Skip to content

Commit

Permalink
ENH: Update ImageToHistogram to use modern dynamic threading.
Browse files Browse the repository at this point in the history 
Change-Id: I979069f0ae1fdde9175a1f86853378db4c86ca5b
  • Loading branch information
@blowekamp @hjmjohnson
blowekamp authored and hjmjohnson committed Nov 6, 2018
1 parent d12b531 commit 95640f7
Show file tree
Hide file tree
Showing 4 changed files with 113 additions and 125 deletions.
14 changes: 8 additions & 6 deletions Modules/Numerics/Statistics/include/itkImageToHistogramFilter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,23 +122,25 @@ class ITK_TEMPLATE_EXPORT ImageToHistogramFilter:public ImageTransformer<TImage>

void GenerateData() override;
void BeforeThreadedGenerateData() override;
void ThreadedGenerateData(const RegionType & inputRegionForThread, ThreadIdType threadId) override;
void AfterThreadedGenerateData() override;
static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION ThreaderMinMaxCallback(void *arg);

/** Method that construct the outputs */
using DataObjectPointerArraySizeType = ProcessObject::DataObjectPointerArraySizeType;
using Superclass::MakeOutput;
DataObject::Pointer MakeOutput(DataObjectPointerArraySizeType) override;

virtual void ThreadedComputeMinimumAndMaximum( const RegionType & inputRegionForThread, ThreadIdType threadId );
virtual void ThreadedComputeHistogram(const RegionType &);
virtual void ThreadedComputeMinimumAndMaximum( const RegionType & inputRegionForThread );

std::vector< HistogramPointer > m_Histograms;
std::vector< HistogramMeasurementVectorType > m_Minimums;
std::vector< HistogramMeasurementVectorType > m_Maximums;
std::mutex m_Mutex;

HistogramMeasurementVectorType m_Minimum;
HistogramMeasurementVectorType m_Maximum;

private:
void ApplyMarginalScale( HistogramMeasurementVectorType & min, HistogramMeasurementVectorType & max, HistogramSizeType & size );


};
} // end of namespace Statistics
} // end of namespace itk
Expand Down
178 changes: 68 additions & 110 deletions Modules/Numerics/Statistics/include/itkImageToHistogramFilter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,19 +103,12 @@ ImageToHistogramFilter< TImage >
this->AllocateOutputs();
this->BeforeThreadedGenerateData();

m_Histograms[0] = this->GetOutput(); // just use the main one
m_Histograms[0]->SetClipBinsAtEnds(true);
for (unsigned i=1; i<m_Histograms.size(); i++)
{
m_Histograms[i] = HistogramType::New();
m_Histograms[i]->SetClipBinsAtEnds(true);
}
HistogramType *outputHistogram = this->GetOutput();
outputHistogram->SetClipBinsAtEnds(true);

// the parameter needed to initialize the histogram
unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
const unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
HistogramSizeType size(nbOfComponents);
HistogramMeasurementVectorType min(nbOfComponents);
HistogramMeasurementVectorType max(nbOfComponents);
if( this->GetHistogramSizeInput() )
{
// user provided value
Expand All @@ -128,145 +121,83 @@ ImageToHistogramFilter< TImage >
}

this->UpdateProgress(0.01f);
//HistogramType * hist = m_Histograms[threadId];

if( this->GetAutoMinimumMaximumInput() && this->GetAutoMinimumMaximum() )
{
// we have to compute the minimum and maximum values
this->ClassicMultiThread(this->ThreaderMinMaxCallback); //calls ThreadedComputeMinimumAndMaximum
this->GetMultiThreader()-> template ParallelizeImageRegion<ImageType::ImageDimension>(
this->GetInput()->GetRequestedRegion(),
[this](const RegionType & inputRegionForThread)
{this->ThreadedComputeMinimumAndMaximum(inputRegionForThread);},
this);

this->UpdateProgress(0.3f);

min = m_Minimums[0];
max = m_Maximums[0];
for( unsigned int t=1; t<m_Minimums.size(); t++ )
{
for( unsigned int i=0; i<nbOfComponents; i++ )
{
min[i] = std::min( min[i], m_Minimums[t][i] );
max[i] = std::max( max[i], m_Maximums[t][i] );
}
}
this->ApplyMarginalScale( min, max, size );
this->ApplyMarginalScale( m_Minimum, m_Maximum, size );

}
else
{
if( this->GetHistogramBinMinimumInput() )
{
min = this->GetHistogramBinMinimum();
m_Minimum = this->GetHistogramBinMinimum();
}
else
{
min.Fill( NumericTraits<ValueType>::NonpositiveMin() - 0.5 );
m_Minimum.Fill( NumericTraits<ValueType>::NonpositiveMin() - 0.5 );
}
if( this->GetHistogramBinMaximumInput() )
{
max = this->GetHistogramBinMaximum();
m_Maximum = this->GetHistogramBinMaximum();
}
else
{
max.Fill( NumericTraits<ValueType>::max() + 0.5 );
m_Maximum.Fill( NumericTraits<ValueType>::max() + 0.5 );
// this->ApplyMarginalScale( min, max, size );
}
}

// finally, initialize the histograms
for (unsigned i=0; i<m_Histograms.size(); i++)
{
m_Histograms[i]->SetMeasurementVectorSize(nbOfComponents);
m_Histograms[i]->Initialize(size, min, max);
}
outputHistogram->SetMeasurementVectorSize(nbOfComponents);
outputHistogram->Initialize(size, m_Minimum, m_Maximum);

this->ClassicMultiThread(this->ThreaderCallback); //parallelizes ThreadedGenerateData
this->GetMultiThreader()-> template ParallelizeImageRegion<ImageType::ImageDimension>(
this->GetInput()->GetRequestedRegion(),
[this](const RegionType & inputRegionForThread)
{this->ThreadedComputeHistogram(inputRegionForThread);},
this);
this->UpdateProgress(0.8f);

this->AfterThreadedGenerateData();
this->UpdateProgress(1.0f);
}


template< typename TImage >
ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
ImageToHistogramFilter< TImage >
::ThreaderMinMaxCallback(void *arg)
{
using ThreadInfo = MultiThreaderBase::WorkUnitInfo;
ThreadInfo * threadInfo = static_cast<ThreadInfo *>(arg);
ThreadIdType threadId = threadInfo->WorkUnitID;
ThreadIdType threadCount = threadInfo->NumberOfWorkUnits;
using FilterStruct = typename ImageTransformer< TImage >::ThreadStruct;
FilterStruct* str = (FilterStruct *)(threadInfo->UserData);
Self* filter = static_cast<Self*>(str->Filter.GetPointer());

// execute the actual method with appropriate output region
// first find out how many pieces extent can be split into.
typename TImage::RegionType splitRegion;
ThreadIdType total = filter->SplitRequestedRegion(threadId, threadCount, splitRegion);

if ( threadId < total )
{
filter->ThreadedComputeMinimumAndMaximum(splitRegion, threadId);
}
// else don't use this thread. Threads were not split conveniently.
return ITK_THREAD_RETURN_DEFAULT_VALUE;
}

template< typename TImage >
void
ImageToHistogramFilter< TImage >
::BeforeThreadedGenerateData()
{
// find the actual number of threads
long nbOfThreads = this->GetNumberOfWorkUnits();
this->GetMultiThreader()->SetNumberOfWorkUnits( nbOfThreads );
// number of work units could be clamped, so get actual number
nbOfThreads = this->GetMultiThreader()->GetNumberOfWorkUnits();
// number of work units can be constrained by the region size,
// so call the SplitRequestedRegion
// to get the real number of threads which will be used
RegionType splitRegion; // dummy region - just to call the following method
nbOfThreads = this->SplitRequestedRegion( 0, nbOfThreads, splitRegion );

// and allocate one histogram per thread
m_Histograms.resize(nbOfThreads);
m_Minimums.resize(nbOfThreads);
m_Maximums.resize(nbOfThreads);
const unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
m_Minimum = HistogramMeasurementVectorType( nbOfComponents );
m_Maximum = HistogramMeasurementVectorType( nbOfComponents );

m_Minimum.Fill( NumericTraits<ValueType>::max() );
m_Maximum.Fill( NumericTraits<ValueType>::NonpositiveMin() );
}

template< typename TImage >
void
ImageToHistogramFilter< TImage >
::AfterThreadedGenerateData()
{
// group the results in the output histogram
HistogramType * hist = m_Histograms[0];
typename HistogramType::IndexType index;
for( unsigned int i=1; i<m_Histograms.size(); i++ )
{
using HistogramIterator = typename HistogramType::ConstIterator;

HistogramIterator hit = m_Histograms[i]->Begin();
HistogramIterator end = m_Histograms[i]->End();
while ( hit != end )
{
hist->GetIndex( hit.GetMeasurementVector(), index);
hist->IncreaseFrequencyOfIndex( index, hit.GetFrequency() );
++hit;
}
}

// and drop the temporary histograms
m_Histograms.clear();
m_Minimums.clear();
m_Maximums.clear();
}


template< typename TImage >
void
ImageToHistogramFilter< TImage >
::ThreadedComputeMinimumAndMaximum(const RegionType & inputRegionForThread, ThreadIdType threadId )
::ThreadedComputeMinimumAndMaximum(const RegionType & inputRegionForThread)
{
unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
const unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
HistogramMeasurementVectorType min( nbOfComponents );
HistogramMeasurementVectorType max( nbOfComponents );

Expand All @@ -287,16 +218,29 @@ ImageToHistogramFilter< TImage >
}
++inputIt;
}
m_Minimums[threadId] = min;
m_Maximums[threadId] = max;

std::lock_guard<std::mutex> mutexHolder( m_Mutex );
for( unsigned int i=0; i<nbOfComponents; i++ )
{
m_Minimum[i] = std::min( m_Minimum[i], min[i] );
m_Maximum[i] = std::max( m_Maximum[i], max[i] );
}
}

template< typename TImage >
void
ImageToHistogramFilter< TImage >
::ThreadedGenerateData(const RegionType & inputRegionForThread, ThreadIdType threadId)
::ThreadedComputeHistogram(const RegionType & inputRegionForThread )
{
unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
const unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
const HistogramType *outputHistogram = this->GetOutput();

HistogramPointer histogram = HistogramType::New();
histogram->SetClipBinsAtEnds(outputHistogram->GetClipBinsAtEnds());
histogram->SetMeasurementVectorSize(nbOfComponents);
histogram->Initialize(outputHistogram->GetSize(), m_Minimum, m_Maximum);


ImageRegionConstIterator< TImage > inputIt( this->GetInput(), inputRegionForThread );
inputIt.GoToBegin();
HistogramMeasurementVectorType m( nbOfComponents );
Expand All @@ -306,18 +250,35 @@ ImageToHistogramFilter< TImage >
{
const PixelType & p = inputIt.Get();
NumericTraits<PixelType>::AssignToArray( p, m );
m_Histograms[threadId]->GetIndex( m, index );
m_Histograms[threadId]->IncreaseFrequencyOfIndex( index, 1 );
histogram->GetIndex( m, index );
histogram->IncreaseFrequencyOfIndex( index, 1 );
++inputIt;
}


std::lock_guard<std::mutex> mutexHolder( m_Mutex );

// reduce the results in the output histogram
HistogramType *writableOutputHistogram = this->GetOutput();

using HistogramIterator = typename HistogramType::ConstIterator;

HistogramIterator hit = histogram->Begin();
HistogramIterator end = histogram->End();
while ( hit != end )
{
writableOutputHistogram->GetIndex( hit.GetMeasurementVector(), index);
writableOutputHistogram->IncreaseFrequencyOfIndex( index, hit.GetFrequency() );
++hit;
}
}

template< typename TImage >
void
ImageToHistogramFilter< TImage >
::ApplyMarginalScale( HistogramMeasurementVectorType & min, HistogramMeasurementVectorType & max, HistogramSizeType & size )
{
unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
const unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
bool clipHistograms = true;
for ( unsigned int i = 0; i < nbOfComponents; i++ )
{
Expand Down Expand Up @@ -376,10 +337,7 @@ ImageToHistogramFilter< TImage >
}
if( clipHistograms == false )
{
for( unsigned int i=0; i<m_Histograms.size(); i++ )
{
m_Histograms[i]->SetClipBinsAtEnds(false);
}
this->GetOutput()->SetClipBinsAtEnds(false);
}
}

Expand Down
4 changes: 2 additions & 2 deletions Modules/Numerics/Statistics/include/itkMaskedImageToHistogramFilter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,8 +84,8 @@ class ITK_TEMPLATE_EXPORT MaskedImageToHistogramFilter:public ImageToHistogramFi
MaskedImageToHistogramFilter();
~MaskedImageToHistogramFilter() override = default;

void ThreadedGenerateData( const RegionType & inputRegionForThread, ThreadIdType threadId ) override;
void ThreadedComputeMinimumAndMaximum( const RegionType & inputRegionForThread, ThreadIdType threadId ) override;
void ThreadedComputeHistogram( const RegionType & inputRegionForThread ) override;
void ThreadedComputeMinimumAndMaximum( const RegionType & inputRegionForThread) override;
};
} // end of namespace Statistics
} // end of namespace itk
Expand Down
42 changes: 35 additions & 7 deletions Modules/Numerics/Statistics/include/itkMaskedImageToHistogramFilter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ MaskedImageToHistogramFilter< TImage, TMaskImage >
template< typename TImage, typename TMaskImage >
void
MaskedImageToHistogramFilter< TImage, TMaskImage >
::ThreadedComputeMinimumAndMaximum( const RegionType & inputRegionForThread, ThreadIdType threadId )
::ThreadedComputeMinimumAndMaximum( const RegionType & inputRegionForThread )
{
unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
HistogramMeasurementVectorType min( nbOfComponents );
Expand Down Expand Up @@ -68,16 +68,27 @@ MaskedImageToHistogramFilter< TImage, TMaskImage >
++inputIt;
++maskIt;
}
this->m_Minimums[threadId] = min;
this->m_Maximums[threadId] = max;
std::lock_guard<std::mutex> mutexHolder( this->m_Mutex );
for( unsigned int i=0; i<nbOfComponents; i++ )
{
this->m_Minimum[i] = std::min( this->m_Minimum[i], min[i] );
this->m_Maximum[i] = std::max( this->m_Maximum[i], max[i] );
}
}

template< typename TImage, typename TMaskImage >
void
MaskedImageToHistogramFilter< TImage, TMaskImage >
::ThreadedGenerateData(const RegionType & inputRegionForThread, ThreadIdType threadId)
::ThreadedComputeHistogram(const RegionType & inputRegionForThread)
{
unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
const unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
const HistogramType *outputHistogram = this->GetOutput();

HistogramPointer histogram = HistogramType::New();
histogram->SetClipBinsAtEnds(outputHistogram->GetClipBinsAtEnds());
histogram->SetMeasurementVectorSize(nbOfComponents);
histogram->Initialize(outputHistogram->GetSize(), this->m_Minimum, this->m_Maximum);

ImageRegionConstIterator< TImage > inputIt( this->GetInput(), inputRegionForThread );
ImageRegionConstIterator< TMaskImage > maskIt( this->GetMaskImage(), inputRegionForThread );
inputIt.GoToBegin();
Expand All @@ -92,12 +103,29 @@ MaskedImageToHistogramFilter< TImage, TMaskImage >
{
const PixelType & p = inputIt.Get();
NumericTraits<PixelType>::AssignToArray( p, m );
this->m_Histograms[threadId]->GetIndex( m, index );
this->m_Histograms[threadId]->IncreaseFrequencyOfIndex( index, 1 );
histogram->GetIndex( m, index );
histogram->IncreaseFrequencyOfIndex( index, 1 );
}
++inputIt;
++maskIt;
}


std::lock_guard<std::mutex> mutexHolder( this->m_Mutex );

// reduce the results in the output histogram
HistogramType *writableOutputHistogram = this->GetOutput();

using HistogramIterator = typename HistogramType::ConstIterator;

HistogramIterator hit = histogram->Begin();
HistogramIterator end = histogram->End();
while ( hit != end )
{
writableOutputHistogram->GetIndex( hit.GetMeasurementVector(), index);
writableOutputHistogram->IncreaseFrequencyOfIndex( index, hit.GetFrequency() );
++hit;
}
}

} // end of namespace Statistics
Expand Down

0 comments on commit 95640f7

Please sign in to comment.