Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

WIP: PERF: Use NeighborhoodRange for metric image computation #98

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

WIP: PERF: Use NeighborhoodRange for metric image computation #98

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 52 additions & 10 deletions include/itkBlockMatchingNormalizedCrossCorrelationNeighborhoodIteratorMetricImageFilter.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,14 @@
#include "itkConstantBoundaryCondition.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkImageKernelOperator.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkNeighborhoodInnerProduct.h"
#include "itkConstantBoundaryImageNeighborhoodPixelAccessPolicy.h"
#include "itkShapedImageNeighborhoodRange.h"
#include "itkImageNeighborhoodOffsets.h"
#include <numeric>

namespace itk
{
Expand Down Expand Up @@ -88,7 +92,8 @@ NormalizedCrossCorrelationNeighborhoodIteratorMetricImageFilter< TFixedImage, TM
typedef ImageRegionIterator< MetricImageType > MetricIteratorType;
MetricImageRegionType metricRegion;
typename MovingImageType::IndexType fitIndex;
typedef ImageRegionConstIterator< MetricImageType > MetricConstIteratorType;
//typedef ImageRegionConstIterator< MetricImageType > MetricConstIteratorType;
typedef ImageRegionConstIteratorWithIndex< MetricImageType > MetricConstIteratorType;
typedef ConstantBoundaryCondition< MetricImageType > BoundaryConditionType;
BoundaryConditionType boundaryCondition;
boundaryCondition.SetConstant( NumericTraits< MetricImagePixelType >::Zero );
Expand Down Expand Up @@ -116,9 +121,10 @@ NormalizedCrossCorrelationNeighborhoodIteratorMetricImageFilter< TFixedImage, TM
typename FaceCalculatorType::FaceListType faceList = faceCalculator(
movingPtr, this->m_MovingImageRegion, radius );
typename FaceCalculatorType::FaceListType::iterator fit;
const MetricImagePixelType negativeOne = -1 * NumericTraits< MetricImagePixelType >::One;
const MetricImagePixelType positiveOne = NumericTraits< MetricImagePixelType >::One;
MetricImagePixelType normXcorr;
constexpr MetricImagePixelType negativeOne = -1 * NumericTraits< MetricImagePixelType >::One;
constexpr MetricImagePixelType positiveOne = NumericTraits< MetricImagePixelType >::One;
const std::vector<Offset<ImageDimension>> offsets =
Experimental::GenerateRectangularImageNeighborhoodOffsets(radius);
for( fit = faceList.begin(); fit != faceList.end(); ++fit )
{
(*fit).Crop( outputRegion );
Expand All @@ -136,24 +142,60 @@ NormalizedCrossCorrelationNeighborhoodIteratorMetricImageFilter< TFixedImage, TM
MetricIteratorType metricIt( metricPtr, metricRegion );

MetricConstIteratorType denomIt( denom, *fit );
NeighborhoodIteratorType movingNeighborIt( radius, movingMinusMean, *fit );
movingNeighborIt.OverrideBoundaryCondition( &boundaryCondition );


// for every pixel in the output
for( metricIt.GoToBegin(), denomIt.GoToBegin(), movingNeighborIt.GoToBegin();
//NeighborhoodIteratorType movingNeighborIt( radius, movingMinusMean, *fit );
//movingNeighborIt.OverrideBoundaryCondition( &boundaryCondition );
//for( metricIt.GoToBegin(), denomIt.GoToBegin(), movingNeighborIt.GoToBegin();
//!metricIt.IsAtEnd();
//++metricIt, ++denomIt, ++movingNeighborIt )
//{
//if( !(denomIt.Get() == NumericTraits< MetricImagePixelType >::Zero) )
//{
//const MetricImagePixelType normXcorr = innerProduct( movingNeighborIt, fixedKernelOperator ) / denomIt.Get();

//// Why does this happen? Bug? Funky floating point behavior?
//if( normXcorr < negativeOne )
//{
//metricIt.Set( negativeOne );
//}
//else if ( normXcorr > positiveOne )
//{
//metricIt.Set( positiveOne );
//}
//else
//{
//metricIt.Set( normXcorr );
//}
//}
//}
for( metricIt.GoToBegin(), denomIt.GoToBegin();
!metricIt.IsAtEnd();
++metricIt, ++denomIt, ++movingNeighborIt )
++metricIt, ++denomIt)
{
if( !(denomIt.Get() == NumericTraits< MetricImagePixelType >::Zero) )
{
normXcorr = innerProduct( movingNeighborIt, fixedKernelOperator ) / denomIt.Get();
using RangeType = Experimental::ShapedImageNeighborhoodRange<const MetricImageType,
Experimental::ConstantBoundaryImageNeighborhoodPixelAccessPolicy<const MetricImageType> >;
const RangeType movingRange{ *movingMinusMean, denomIt.GetIndex(), offsets };
const RangeType kernelRange{ *fixedMinusMean, denomIt.GetIndex(), offsets };
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it really intended that the kernel is set on a new location (denomIt.GetIndex()) with every iteration?


const MetricImagePixelType normXcorr = std::inner_product( movingRange.begin(), movingRange.end(), kernelRange.begin(), 0.0 ) / denomIt.Get();
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have to admit I could not get optimal performance when using std::inner_product at https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Core/ImageFunction/include/itkGaussianDerivativeImageFunction.hxx#L219, so instead I just manually wrote the inner product calculation in a few lines of code. In this case, it might look as follows:

auto normXcorr = NumericTraits<MetricImagePixelType>::Zero;
auto movingRangeIterator = movingRange.cbegin();

for (const auto& kernelValue : kernelRange)
{
  normXcorr += kernelValue * (*movingRangeIterator);
  ++movingRangeIterator;
}
normXcorr /= denomIt.Get();


// Why does this happen? Bug? Funky floating point behavior?
if( normXcorr < negativeOne )
{
metricIt.Set( negativeOne );
}
else if ( normXcorr > positiveOne )
{
metricIt.Set( positiveOne );
}
else
{
metricIt.Set( normXcorr );
}
}
}
}
Expand Down