/usr/include/ITK-4.5/itkTriangleThresholdCalculator.hxx is in libinsighttoolkit4-dev 4.5.0-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 | /*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef __itkTriangleThresholdCalculator_hxx
#define __itkTriangleThresholdCalculator_hxx
#include "itkTriangleThresholdCalculator.h"
#include "itkProgressReporter.h"
#include "vnl/vnl_math.h"
#include <algorithm>
namespace itk
{
/*
* Compute the Triangle's threshold
*/
template<typename THistogram, typename TOutput>
void
TriangleThresholdCalculator<THistogram, TOutput>
::GenerateData(void)
{
const HistogramType * histogram = this->GetInput();
// histogram->Print(std::cout);
if ( histogram->GetTotalFrequency() == 0 )
{
itkExceptionMacro(<< "Histogram is empty");
}
ProgressReporter progress(this, 0, histogram->GetSize(0) );
if( histogram->GetSize(0) == 1 )
{
this->GetOutput()->Set( static_cast<OutputType>(histogram->GetMeasurement(0,0)) );
}
SizeValueType size = histogram->GetSize(0);
// create a histogram
std::vector<double> cumSum(size, 0.0);
std::vector<double> triangle(size, 0.0);
// Triangle method needs the maximum and minimum indexes
// Minimum indexes for this purpose are poorly defined - can't just
// take a index with zero entries.
double Mx = itk::NumericTraits<double>::min();
IndexValueType MxIdx = 0;
for ( SizeValueType j = 0; j < size; j++ )
{
if ( histogram->GetFrequency(j, 0) > Mx )
{
MxIdx = j;
Mx = histogram->GetFrequency(j, 0);
}
}
cumSum[0] = histogram->GetFrequency(0, 0);
for ( SizeValueType j = 1; j < size; j++ )
{
cumSum[j] = histogram->GetFrequency(j, 0) + cumSum[j-1];
}
typename HistogramType::MeasurementVectorType onePC(1), nnPC(1);
onePC.Fill(histogram->Quantile(0, 0.01));
typename HistogramType::IndexType localIndex;
histogram->GetIndex(onePC,localIndex);
const IndexValueType onePCIdx = localIndex[0];
nnPC.Fill(histogram->Quantile(0, 0.99));
histogram->GetIndex(nnPC,localIndex);
const IndexValueType nnPCIdx = localIndex[0];
// figure out which way we are looking - we want to construct our
// line between the max index and the further of 1% and 99%
IndexValueType ThreshIdx = 0;
if (fabs((float)MxIdx - (float)onePCIdx) > fabs((float)MxIdx - (float)nnPCIdx))
{
// line to 1 %
double slope = Mx / ( MxIdx - onePCIdx );
for (IndexValueType k = onePCIdx; k < MxIdx; k++)
{
float line = (slope*(k-onePCIdx));
triangle[k]= line - histogram->GetFrequency(k);
}
ThreshIdx = onePCIdx + std::distance(&(triangle[onePCIdx]), std::max_element(&(triangle[onePCIdx]), &(triangle[MxIdx])));
}
else
{
// line to 99 %
double slope = -Mx / ( nnPCIdx - MxIdx );
for (IndexValueType k = MxIdx; k < nnPCIdx; k++)
{
float line = (slope*(k-MxIdx) + Mx);
triangle[k]= line - histogram->GetFrequency(k);
}
ThreshIdx = MxIdx + std::distance(&(triangle[MxIdx]), std::max_element(&(triangle[MxIdx]), &(triangle[nnPCIdx])));
}
this->GetOutput()->Set( static_cast<OutputType>( histogram->GetMeasurement( ThreshIdx + 1, 0 ) ) );
}
} // end namespace itk
#endif
|