Hello !
So I solved my problems these ways, but I am not sure whether these are the
correct ways (I seems to work thought).
1) This bug was due to the outputlowpass which was never allocated when
direction != 0 (direction was 1 at some time during the execution because
it was a 3D file) and m_SubsampleImageFactor = 1. Thus I moved the lines
previously in the if (direction == 0) clause
outputLowPass->SetRegions(
outputImageRegion);
outputLowPass->Allocate();
outputLowPass->FillBuffer(0);
before the lines
if (m_SubsampleImageFactor > 1)
{
input = m_InternalImages[direction][idx / (1 << (direction
+ 1))];
if (direction != 0)
{
outputLowPass = m_InternalImages[direction - 1][idx /
(1 << direction)];
outputHighPass = m_InternalImages[direction - 1][idx /
(1 << direction) + 1];
}
}
2) I changed
for (unsigned int i = 0; i < m_InternalImages[direction].size(); ++i)
to
for (unsigned int i = 0; i < m_InternalImages[InputImageDimension - 2 -
direction].size(); ++i)
3) I changed
static_cast<OutputImageType*>(m_InternalImages[direction - 2][idx / (1 <<
(direction - 1))])
to
static_cast<OutputImageType*>(m_InternalImages[direction][idx / (1 <<
(direction + 1))])
Attached is the itk version of the otbWaveletFilterBank.txx with my
modifications (it works with both 2D and 3D images, in single and
multiresolution).
I have another problem with the threads, it does not seem to work if I set
more than 1 thread (on 3D output image, I have a black plane in the midle
of the z dimension but the other parts of the image are good).
Any ideas on how to solve this problem?
Thanks in advance !
Regards,
Flo
--
--
Check the OTB FAQ at
http://www.orfeo-toolbox.org/FAQ.html
You received this message because you are subscribed to the Google
Groups "otb-users" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/otb-users?hl=en
---
You received this message because you are subscribed to the Google Groups
"otb-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
Copyright (c) Institut Mines-Telecom. All rights reserved.
See IMTCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef itkWaveletFilterBank_hxx
#define itkWaveletFilterBank_hxx
#include "itkWaveletFilterBank.h"
#include "itkSubsampleImageFilter.h"
#include "itkSubsampledImageRegionConstIterator.h"
#include <itkPeriodicBoundaryCondition.h>
// FIXME
#define __myDebug__ 0
namespace itk {
/**
* Template Specialization for the Wavelet::FORWARD case
*/
template <class TInputImage, class TOutputImage, class TWaveletOperator>
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::WaveletFilterBank()
{
this->SetNumberOfRequiredInputs(1);
this->SetNumberOfRequiredInputs(1);
unsigned int numOfOutputs = 1 << InputImageDimension;
this->SetNumberOfRequiredOutputs(numOfOutputs);
for (unsigned int i = 0; i < numOfOutputs; ++i)
{
this->SetNthOutput(i, OutputImageType::New());
}
m_UpSampleFilterFactor = 0;
m_SubsampleImageFactor = 1;
this->SetNumberOfThreads(1);//
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::GenerateOutputInformation()
{
Superclass::GenerateOutputInformation();
if (GetSubsampleImageFactor() == 1) return;
#if __myDebug__
itkDebugMacro(<< " down sampling output regions by a factor of " << GetSubsampleImageFactor());
itkDebugMacro(<< "initial region " << this->GetInput()->GetLargestPossibleRegion().GetSize()[0]
<< "," << this->GetInput()->GetLargestPossibleRegion().GetSize()[1]);
#endif
OutputImageRegionType newRegion;
this->CallCopyInputRegionToOutputRegion(newRegion, this->GetInput()->GetLargestPossibleRegion());
for (unsigned int i = 0; i < this->GetNumberOfOutputs(); ++i)
{
this->GetOutput(i)->SetRegions(newRegion);
}
#if __myDebug__
itkDebugMacro(<< "new region output " << newRegion.GetSize()[0] << "," << newRegion.GetSize()[1]);
#endif
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::GenerateInputRequestedRegion()
throw (itk::InvalidRequestedRegionError)
{
Superclass::GenerateInputRequestedRegion();
// Filter length calculation
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(0);
lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
lowPassOperator.CreateDirectional();
unsigned int radius = lowPassOperator.GetRadius()[0];
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(0);
highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
highPassOperator.CreateDirectional();
if (radius < highPassOperator.GetRadius()[0]) radius = highPassOperator.GetRadius()[0];
// Get the requested region and pad it
InputImagePointerType input = const_cast<InputImageType*>(this->GetInput());
InputImageRegionType inputRegion = input->GetRequestedRegion();
inputRegion.PadByRadius(radius);
if (inputRegion.Crop(input->GetLargestPossibleRegion()))
{
input->SetRequestedRegion(inputRegion);
}
else
{
input->SetRequestedRegion(inputRegion);
itk::InvalidRequestedRegionError err(__FILE__, __LINE__);
err.SetLocation(ITK_LOCATION);
err.SetDescription("Requested region is (at least partially) outside the largest possible region.");
err.SetDataObject(input);
throw err;
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::BeforeThreadedGenerateData()
{
if (m_SubsampleImageFactor > 1)
{
// Check the dimension
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
if ((m_SubsampleImageFactor
* (this->GetInput()->GetRequestedRegion().GetSize()[i] / m_SubsampleImageFactor))
!= this->GetInput()->GetRequestedRegion().GetSize()[i])
{
itk::InvalidRequestedRegionError err(__FILE__, __LINE__);
err.SetLocation(ITK_LOCATION);
err.SetDescription("Requested region dimension cannot be used in multiresolution analysis (crop it).");
err.SetDataObject(const_cast<InputImageType*>(this->GetInput()));
throw err;
}
}
if (InputImageDimension > 1)
{
// Internal images will be used only if m_SubsampledInputImages != 1
m_InternalImages.resize(InputImageDimension - 1);
for (unsigned int i = 0; i < m_InternalImages.size(); ++i)
{
// the size is linked to the SubsampleImageFactor that is assume to be 2!!!
m_InternalImages[InputImageDimension - 2 - i].resize(1 << (i + 1));
}
OutputImageRegionType intermediateRegion;
this->Superclass::CallCopyInputRegionToOutputRegion(intermediateRegion,
this->GetInput()->GetLargestPossibleRegion());
AllocateInternalData(intermediateRegion);
}
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::AllocateInternalData
(const OutputImageRegionType& outputRegion)
{
OutputImageRegionType smallerRegion;
OutputImageRegionType largerRegion = outputRegion;
for (unsigned int direction = 0; direction < InputImageDimension - 1; direction++)
{
this->CallCopyInputRegionToOutputRegion(InputImageDimension - 1 - direction,
smallerRegion, largerRegion);
// for (unsigned int i = 0; i < m_InternalImages[direction].size(); ++i)
for (unsigned int i = 0; i < m_InternalImages[InputImageDimension - 2 - direction].size(); ++i)
{
m_InternalImages[InputImageDimension - 2 - direction][i] = OutputImageType::New();
m_InternalImages[InputImageDimension - 2 - direction][i]->SetRegions(smallerRegion);
m_InternalImages[InputImageDimension - 2 - direction][i]->Allocate();
m_InternalImages[InputImageDimension - 2 - direction][i]->FillBuffer(0);
}
largerRegion = smallerRegion;
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::AfterThreadedGenerateData()
{
if (m_SubsampleImageFactor > 1 && InputImageDimension > 1)
{
m_InternalImages.clear();
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::CallCopyOutputRegionToInputRegion
(InputImageRegionType& destRegion, const OutputImageRegionType& srcRegion)
{
Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
if (GetSubsampleImageFactor() > 1)
{
OutputIndexType srcIndex = srcRegion.GetIndex();
OutputSizeType srcSize = srcRegion.GetSize();
InputIndexType destIndex;
InputSizeType destSize;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
destIndex[i] = srcIndex[i] * GetSubsampleImageFactor();
destSize[i] = srcSize[i] * GetSubsampleImageFactor();
}
destRegion.SetIndex(destIndex);
destRegion.SetSize(destSize);
#if 0
// Contrairement a Wavelet::INVERSE, ici ca ne sera a rien apparemment...
// Region Padding
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(0);
lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
lowPassOperator.CreateDirectional();
unsigned long radius[InputImageDimension];
radius[0] = lowPassOperator.GetRadius()[0];
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(0);
highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
highPassOperator.CreateDirectional();
if (radius[0] < highPassOperator.GetRadius()[0]) radius[0] = highPassOperator.GetRadius()[0];
for (unsigned int i = 1; i < InputImageDimension; ++i)
radius[i] = 0;
InputImageRegionType paddedRegion = destRegion;
paddedRegion.PadByRadius(radius);
if (paddedRegion.Crop(this->GetInput()->GetLargestPossibleRegion()))
{
destRegion = paddedRegion;
}
#endif
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::CallCopyOutputRegionToInputRegion
(unsigned int direction,
InputImageRegionType& destRegion, const OutputImageRegionType& srcRegion)
{
Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
if (GetSubsampleImageFactor() > 1)
{
OutputIndexType srcIndex = srcRegion.GetIndex();
OutputSizeType srcSize = srcRegion.GetSize();
InputIndexType destIndex;
InputSizeType destSize;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
if (i == direction)
{
destIndex[i] = srcIndex[i] * GetSubsampleImageFactor();
destSize[i] = srcSize[i] * GetSubsampleImageFactor();
}
else
{
destIndex[i] = srcIndex[i];
destSize[i] = srcSize[i];
}
}
destRegion.SetIndex(destIndex);
destRegion.SetSize(destSize);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::CallCopyInputRegionToOutputRegion
(OutputImageRegionType& destRegion, const InputImageRegionType& srcRegion)
{
Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
if (GetSubsampleImageFactor() > 1)
{
typename InputImageRegionType::IndexType srcIndex = srcRegion.GetIndex();
typename InputImageRegionType::SizeType srcSize = srcRegion.GetSize();
typename OutputImageRegionType::IndexType destIndex;
typename OutputImageRegionType::SizeType destSize;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
// TODO: This seems not right in odd index cases
destIndex[i] = srcIndex[i] / GetSubsampleImageFactor();
destSize[i] = srcSize[i] / GetSubsampleImageFactor();
}
destRegion.SetIndex(destIndex);
destRegion.SetSize(destSize);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::CallCopyInputRegionToOutputRegion
(unsigned int direction,
OutputImageRegionType& destRegion, const InputImageRegionType& srcRegion)
{
Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
if (GetSubsampleImageFactor() > 1)
{
typename InputImageRegionType::IndexType srcIndex = srcRegion.GetIndex();
typename InputImageRegionType::SizeType srcSize = srcRegion.GetSize();
typename OutputImageRegionType::IndexType destIndex;
typename OutputImageRegionType::SizeType destSize;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
if (i == direction)
{
// TODO: This seems not right in odd index cases
destIndex[i] = srcIndex[i] / GetSubsampleImageFactor();
destSize[i] = srcSize[i] / GetSubsampleImageFactor();
}
else
{
destIndex[i] = srcIndex[i];
destSize[i] = srcSize[i];
}
}
destRegion.SetIndex(destIndex);
destRegion.SetSize(destSize);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::ThreadedGenerateData
(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
{
unsigned int dir = InputImageDimension - 1;
if ((1 << dir) >= static_cast<int>(this->GetNumberOfOutputs()))
{
std::ostringstream msg;
msg << "Output number 1<<" << dir << " = " << (1 << dir) << " not allocated\n";
msg << "Number of excpected outputs " << this->GetNumberOfOutputs() << "\n";
throw itk::ExceptionObject(__FILE__, __LINE__, msg.str().c_str(), ITK_LOCATION);
}
itk::ProgressReporter reporter(this, threadId,
outputRegionForThread.GetNumberOfPixels() * this->GetNumberOfOutputs() * 2);
const InputImageType * input = this->GetInput();
InputImageRegionType inputRegionForThread;
this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
OutputImagePointerType outputLowPass = this->GetOutput(0);
OutputImagePointerType outputHighPass = this->GetOutput(1 << dir);
// On multidimensional case, if m_SubsampleImageFactor != 1, we need internal images of different size
if (dir != 0 && m_SubsampleImageFactor > 1)
{
outputLowPass = m_InternalImages[dir - 1][0];
outputHighPass = m_InternalImages[dir - 1][1];
}
// typedef for the iterations over the input image
typedef itk::ConstNeighborhoodIterator<InputImageType> NeighborhoodIteratorType;
typedef itk::NeighborhoodInnerProduct<InputImageType> InnerProductType;
typedef itk::ImageRegionIterator<OutputImageType> IteratorType;
// Prepare the subsampling image factor, if any.
typedef SubsampledImageRegionConstIterator<InputImageType> SubsampleIteratorType;
typename SubsampleIteratorType::IndexType subsampling;
subsampling.Fill(1);
subsampling[dir] = GetSubsampleImageFactor();
// Inner product
InnerProductType innerProduct;
// High pass part calculation
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(dir);
highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
highPassOperator.CreateDirectional();
SubsampleIteratorType subItHighPass(input, inputRegionForThread);
subItHighPass.SetSubsampleFactor(subsampling);
subItHighPass.GoToBegin();
NeighborhoodIteratorType itHighPass(highPassOperator.GetRadius(), input, inputRegionForThread);
itk::PeriodicBoundaryCondition<InputImageType> boundaryCondition;
itHighPass.OverrideBoundaryCondition(&boundaryCondition);
IteratorType outHighPass(outputHighPass, subItHighPass.GenerateOutputInformation());
outHighPass.GoToBegin();
while (!subItHighPass.IsAtEnd() && !outHighPass.IsAtEnd())
{
itHighPass.SetLocation(subItHighPass.GetIndex());
outHighPass.Set(innerProduct(itHighPass, highPassOperator));
++subItHighPass;
++outHighPass;
reporter.CompletedPixel();
}
// Low pass part calculation
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(dir);
lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
lowPassOperator.CreateDirectional();
SubsampleIteratorType subItLowPass(input, inputRegionForThread);
subItLowPass.SetSubsampleFactor(subsampling);
subItLowPass.GoToBegin();
NeighborhoodIteratorType itLowPass(lowPassOperator.GetRadius(), input, inputRegionForThread);
itLowPass.OverrideBoundaryCondition(&boundaryCondition);
IteratorType outLowPass(outputLowPass, subItLowPass.GenerateOutputInformation());
outLowPass.GoToBegin();
while (!subItLowPass.IsAtEnd() && !outLowPass.IsAtEnd())
{
itLowPass.SetLocation(subItLowPass.GetIndex());
outLowPass.Set(innerProduct(itLowPass, lowPassOperator));
++subItLowPass;
++outLowPass;
reporter.CompletedPixel();
}
if (dir > 0)
{
// Note that outputImageRegion correspond to the actual region of (local) input !
OutputImageRegionType outputImageRegion;
this->CallCopyInputRegionToOutputRegion(dir, outputImageRegion, inputRegionForThread);
ThreadedGenerateDataAtDimensionN(0, dir - 1, reporter, outputImageRegion, threadId);
ThreadedGenerateDataAtDimensionN(1 << dir, dir - 1, reporter, outputImageRegion, threadId);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
::ThreadedGenerateDataAtDimensionN
(unsigned int idx, unsigned int direction, itk::ProgressReporter& reporter,
const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
{
// Note that outputRegionForThread correspond to the actual region of input !
OutputImagePointerType input = this->GetOutput(idx);
OutputImagePointerType outputHighPass = this->GetOutput(idx + (1 << direction));
OutputImagePointerType outputLowPass = OutputImageType::New();
OutputImageRegionType outputImageRegion;
this->CallCopyInputRegionToOutputRegion(direction, outputImageRegion, outputRegionForThread);
// The output image has to be allocated
// May be not valid on multithreaded process ???
outputLowPass->SetRegions(outputImageRegion);
outputLowPass->Allocate(); // FIXME Check this line...
outputLowPass->FillBuffer(0);
if (m_SubsampleImageFactor > 1)
{
input = m_InternalImages[direction][idx / (1 << (direction + 1))];
if (direction != 0)
{
outputLowPass = m_InternalImages[direction - 1][idx / (1 << direction)];
outputHighPass = m_InternalImages[direction - 1][idx / (1 << direction) + 1];
}
}
/*if (direction == 0)
{
// The output image has to be allocated
// May be not valid on multithreaded process ???
outputLowPass->SetRegions(outputImageRegion);
outputLowPass->Allocate(); // FIXME Check this line...
outputLowPass->FillBuffer(0);
}*/
/** Inner product iterators */
typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
typedef itk::NeighborhoodInnerProduct<OutputImageType> InnerProductType;
typedef itk::ImageRegionIterator<OutputImageType> IteratorType;
// Prepare the subsampling image factor, if any.
typedef SubsampledImageRegionConstIterator<InputImageType> SubsampleIteratorType;
typename SubsampleIteratorType::IndexType subsampling;
subsampling.Fill(1);
subsampling[direction] = GetSubsampleImageFactor();
// Inner products
InnerProductType innerProduct;
// High pass part calculation
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(direction);
highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
highPassOperator.CreateDirectional();
SubsampleIteratorType subItHighPass(input, outputRegionForThread);
subItHighPass.SetSubsampleFactor(subsampling);
subItHighPass.GoToBegin();
NeighborhoodIteratorType itHighPass(highPassOperator.GetRadius(), input, outputRegionForThread);
itk::PeriodicBoundaryCondition<OutputImageType> boundaryCondition;
itHighPass.OverrideBoundaryCondition(&boundaryCondition);
IteratorType outHighPass(outputHighPass, subItHighPass.GenerateOutputInformation());
outHighPass.GoToBegin();
while (!subItHighPass.IsAtEnd() && !outHighPass.IsAtEnd())
{
itHighPass.SetLocation(subItHighPass.GetIndex());
outHighPass.Set(innerProduct(itHighPass, highPassOperator));
++subItHighPass;
++outHighPass;
reporter.CompletedPixel();
}
// Low pass part calculation
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(direction);
lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
lowPassOperator.CreateDirectional();
SubsampleIteratorType subItLowPass(input, outputRegionForThread);
subItLowPass.SetSubsampleFactor(subsampling);
subItLowPass.GoToBegin();
NeighborhoodIteratorType itLowPass(lowPassOperator.GetRadius(), input, outputRegionForThread);
itLowPass.OverrideBoundaryCondition(&boundaryCondition);
IteratorType outLowPass(outputLowPass, subItLowPass.GenerateOutputInformation());
outLowPass.GoToBegin();
while (!subItLowPass.IsAtEnd() && !outLowPass.IsAtEnd())
{
itLowPass.SetLocation(subItLowPass.GetIndex());
outLowPass.Set(innerProduct(itLowPass, lowPassOperator));
++subItLowPass;
++outLowPass;
reporter.CompletedPixel();
}
// Swap input and lowPassOutput
itk::ImageRegionConstIterator<OutputImageType> inIt(outputLowPass, outputImageRegion);
IteratorType outIt((direction != 0 && m_SubsampleImageFactor > 1) ?
//static_cast<OutputImageType*>(m_InternalImages[direction - 2][idx / (1 << (direction - 1))])
static_cast<OutputImageType*>(m_InternalImages[direction][idx / (1 << (direction + 1))]) //?
: this->GetOutput(idx),
outputImageRegion);
for (inIt.GoToBegin(), outIt.GoToBegin(); !inIt.IsAtEnd(); ++inIt, ++outIt)
{
outIt.Set(inIt.Get());
}
if (direction != 0)
{
ThreadedGenerateDataAtDimensionN(idx, direction - 1, reporter, outputImageRegion, threadId);
ThreadedGenerateDataAtDimensionN(idx + (1 << direction), direction - 1, reporter, outputImageRegion, threadId);
}
}
/**
* Template Specialization for the Wavelet::INVERSE case
*/
template <class TInputImage, class TOutputImage, class TWaveletOperator>
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::WaveletFilterBank()
{
this->SetNumberOfRequiredInputs(1 << InputImageDimension);
m_UpSampleFilterFactor = 0;
m_SubsampleImageFactor = 1;
// TODO: For now, we force the number threads to 1 because there is a bug with multithreading in INVERSE transform
// Resulting in discontinuities in the reconstructed images
this->SetNumberOfThreads(1);
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::GenerateOutputInformation()
{
Superclass::GenerateOutputInformation();
for (unsigned int i = 1; i < this->GetNumberOfInputs(); ++i)
{
for (unsigned int dim = 0; dim < InputImageDimension; dim++)
{
if (this->GetInput(0)->GetLargestPossibleRegion().GetSize()[dim]
!= this->GetInput(i)->GetLargestPossibleRegion().GetSize()[dim])
{
throw itk::ExceptionObject(__FILE__, __LINE__,
"Input images are not of the same dimension", ITK_LOCATION);
}
}
}
#if __myDebug__
itkDebugMacro(<< " up sampling output regions by a factor of " << GetSubsampleImageFactor());
itkDebugMacro(<< "initial region "
<< this->GetInput(0)->GetLargestPossibleRegion().GetSize()[0]
<< "," << this->GetInput(0)->GetLargestPossibleRegion().GetSize()[1]);
#endif
OutputImageRegionType newRegion;
this->CallCopyInputRegionToOutputRegion(newRegion, this->GetInput(0)->GetLargestPossibleRegion());
this->GetOutput()->SetRegions(newRegion);
#if __myDebug__
itkDebugMacro(<< "new region output " << newRegion.GetSize()[0] << "," << newRegion.GetSize()[1]);
#endif
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::GenerateInputRequestedRegion()
throw (itk::InvalidRequestedRegionError)
{
Superclass::GenerateInputRequestedRegion();
// Filter length calculation
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(0);
lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
lowPassOperator.CreateDirectional();
unsigned int radius = lowPassOperator.GetRadius()[0];
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(0);
highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
highPassOperator.CreateDirectional();
if (radius < highPassOperator.GetRadius()[0]) radius = highPassOperator.GetRadius()[0];
// Get the requested regionand pad it
for (unsigned int idx = 0; idx < this->GetNumberOfInputs(); ++idx)
{
InputImagePointerType input = const_cast<InputImageType*>(this->GetInput(idx));
InputImageRegionType inputRegion = input->GetRequestedRegion();
inputRegion.PadByRadius(radius);
if (inputRegion.Crop(input->GetLargestPossibleRegion()))
{
input->SetRequestedRegion(inputRegion);
}
else
{
input->SetRequestedRegion(inputRegion);
itk::InvalidRequestedRegionError err(__FILE__, __LINE__);
err.SetLocation(ITK_LOCATION);
err.SetDescription("Requested region is (at least partially) outside the largest possible region.");
err.SetDataObject(input);
throw err;
}
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::BeforeThreadedGenerateData()
{
if (InputImageDimension > 1)
{
// Internal images will be used only if m_SubsampleImageFactor != 1
m_InternalImages.resize(InputImageDimension - 1);
for (unsigned int i = 0; i < m_InternalImages.size(); ++i)
{
// the size is linked to the SubsampleImageFactor that is assume to be 2!!!
//m_InternalImages[i].resize(1 << (i + 1));
m_InternalImages[InputImageDimension - 2 - i].resize(1 << (i + 1));
}
OutputImageRegionType intermediateRegion;
Superclass::CallCopyInputRegionToOutputRegion(intermediateRegion,
this->GetInput(0)->GetLargestPossibleRegion());
AllocateInternalData(intermediateRegion);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::AllocateInternalData
(const OutputImageRegionType& outputRegion)
{
OutputImageRegionType largerRegion;
OutputImageRegionType smallerRegion = outputRegion;
for (unsigned int direction = 0; direction < InputImageDimension - 1; direction++)
{
this->CallCopyInputRegionToOutputRegion(direction,
largerRegion, smallerRegion);
for (unsigned int i = 0; i < m_InternalImages[direction].size(); ++i)
{
m_InternalImages[direction][i] = OutputImageType::New();
m_InternalImages[direction][i]->SetRegions(largerRegion);
m_InternalImages[direction][i]->Allocate();
m_InternalImages[direction][i]->FillBuffer(0);
}
smallerRegion = largerRegion;
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::AfterThreadedGenerateData()
{
if (m_SubsampleImageFactor > 1 && InputImageDimension > 1)
{
m_InternalImages.clear();
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::CallCopyOutputRegionToInputRegion
(InputImageRegionType& destRegion, const OutputImageRegionType& srcRegion)
{
Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
if (GetSubsampleImageFactor() > 1)
{
OutputIndexType srcIndex = srcRegion.GetIndex();
OutputSizeType srcSize = srcRegion.GetSize();
InputIndexType destIndex;
InputSizeType destSize;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
// TODO: This seems not right in odd index cases
destIndex[i] = srcIndex[i] / GetSubsampleImageFactor();
destSize[i] = srcSize[i] / GetSubsampleImageFactor();
}
destRegion.SetIndex(destIndex);
destRegion.SetSize(destSize);
#if 1
// Region Padding
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(0);
lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
lowPassOperator.CreateDirectional();
typename InputImageRegionType::SizeType radius;
radius[0] = lowPassOperator.GetRadius()[0];
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(0);
highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
highPassOperator.CreateDirectional();
if (radius[0] < highPassOperator.GetRadius()[0]) radius[0] = highPassOperator.GetRadius()[0];
for (unsigned int i = 1; i < InputImageDimension; ++i)
radius[i] = 0;
// for ( unsigned int i = 0; i < InputImageDimension; ++i )
// {
// radius[i] = lowPassOperator.GetRadius()[i];
// if ( radius[i] < highPassOperator.GetRadius()[i] )
// radius[i] = highPassOperator.GetRadius()[i];
// }
InputImageRegionType paddedRegion = destRegion;
paddedRegion.PadByRadius(radius);
if (paddedRegion.Crop(this->GetInput(0)->GetLargestPossibleRegion()))
{
destRegion = paddedRegion;
}
#endif
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::CallCopyInputRegionToOutputRegion
(OutputImageRegionType& destRegion, const InputImageRegionType& srcRegion)
{
Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
if (GetSubsampleImageFactor() > 1)
{
OutputIndexType srcIndex = srcRegion.GetIndex();
OutputSizeType srcSize = srcRegion.GetSize();
InputIndexType destIndex;
InputSizeType destSize;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
destIndex[i] = srcIndex[i] * GetSubsampleImageFactor();
destSize[i] = srcSize[i] * GetSubsampleImageFactor();
}
destRegion.SetIndex(destIndex);
destRegion.SetSize(destSize);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::CallCopyOutputRegionToInputRegion
(unsigned int direction,
InputImageRegionType& destRegion, const OutputImageRegionType& srcRegion)
{
Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
if (GetSubsampleImageFactor() > 1)
{
OutputIndexType srcIndex = srcRegion.GetIndex();
OutputSizeType srcSize = srcRegion.GetSize();
InputIndexType destIndex;
InputSizeType destSize;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
if (i == direction)
{
// TODO: This seems not right in odd index cases
destIndex[i] = srcIndex[i] / GetSubsampleImageFactor();
destSize[i] = srcSize[i] / GetSubsampleImageFactor();
}
else
{
destIndex[i] = srcIndex[i];
destSize[i] = srcSize[i];
}
}
destRegion.SetIndex(destIndex);
destRegion.SetSize(destSize);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::CallCopyInputRegionToOutputRegion
(unsigned int direction,
OutputImageRegionType& destRegion, const InputImageRegionType& srcRegion)
{
Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
if (GetSubsampleImageFactor() > 1)
{
typename InputImageRegionType::IndexType srcIndex = srcRegion.GetIndex();
typename InputImageRegionType::SizeType srcSize = srcRegion.GetSize();
typename OutputImageRegionType::IndexType destIndex;
typename OutputImageRegionType::SizeType destSize;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
if (i == direction)
{
destIndex[i] = srcIndex[i] * GetSubsampleImageFactor();
destSize[i] = srcSize[i] * GetSubsampleImageFactor();
}
else
{
destIndex[i] = srcIndex[i];
destSize[i] = srcSize[i];
}
}
destRegion.SetIndex(destIndex);
destRegion.SetSize(destSize);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::ThreadedGenerateData
(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
{
itk::ProgressReporter reporter(this, threadId,
outputRegionForThread.GetNumberOfPixels() * this->GetNumberOfInputs());
InputImageRegionType inputRegionForThread;
this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
unsigned int dir = 0;
// Low pass part calculation
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(dir);
lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
lowPassOperator.CreateDirectional();
// High pass part calculation
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(dir);
highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
highPassOperator.CreateDirectional();
// typedef for the iterations over the input image
typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
typedef itk::NeighborhoodInnerProduct<OutputImageType> InnerProductType;
// Faces iterations
typename NeighborhoodIteratorType::RadiusType radiusMax;
for (unsigned int idx = 0; idx < OutputImageDimension; ++idx)
{
radiusMax[idx] = lowPassOperator.GetRadius(idx);
if (radiusMax[idx] < highPassOperator.GetRadius(idx)) radiusMax[idx] = highPassOperator.GetRadius(idx);
}
// The multiresolution case requires a SubsampleImageFilter step
if (m_SubsampleImageFactor > 1)
{
for (unsigned int i = 0; i < this->GetNumberOfInputs(); i += 2)
{
int test = this->GetNumberOfInputs();
InputImagePointerType imgLowPass = const_cast<InputImageType*>(this->GetInput(i));
InputImagePointerType imgHighPass = const_cast<InputImageType*>(this->GetInput(i + 1));
OutputImagePointerType outputImage = this->GetOutput();
if (dir != InputImageDimension - 1)
{
outputImage = m_InternalImages[0][i / 2];
}
typedef itk::SubsampleImageFilter<InputImageType, OutputImageType, Wavelet::INVERSE> FilterType;
typename FilterType::InputImageIndexType delta;
delta.Fill(1);
delta[dir] = this->GetSubsampleImageFactor();
InputImagePointerType cropedLowPass = InputImageType::New();
cropedLowPass->SetRegions(inputRegionForThread);
cropedLowPass->Allocate();
cropedLowPass->FillBuffer(0.);
itk::ImageRegionIterator<InputImageType> cropedLowPassIt(cropedLowPass, inputRegionForThread);
itk::ImageRegionIterator<InputImageType> imgLowPassIt(imgLowPass, inputRegionForThread);
for (cropedLowPassIt.GoToBegin(), imgLowPassIt.GoToBegin();
!cropedLowPassIt.IsAtEnd() && !imgLowPassIt.IsAtEnd();
++cropedLowPassIt, ++imgLowPassIt)
{
cropedLowPassIt.Set(imgLowPassIt.Get());
}
typename FilterType::Pointer overSampledLowPass = FilterType::New();
overSampledLowPass->SetInput(cropedLowPass);
overSampledLowPass->SetSubsampleFactor(delta);
overSampledLowPass->Update();
InputImagePointerType cropedHighPass = InputImageType::New();
cropedHighPass->SetRegions(inputRegionForThread);
cropedHighPass->Allocate();
cropedHighPass->FillBuffer(0.);
itk::ImageRegionIterator<InputImageType> cropedHighPassIt(cropedHighPass, inputRegionForThread);
itk::ImageRegionIterator<InputImageType> imgHighPassIt(imgHighPass, inputRegionForThread);
for (cropedHighPassIt.GoToBegin(), imgHighPassIt.GoToBegin();
!cropedHighPassIt.IsAtEnd() && !imgHighPassIt.IsAtEnd();
++cropedHighPassIt, ++imgHighPassIt)
{
cropedHighPassIt.Set(imgHighPassIt.Get());
}
typename FilterType::Pointer overSampledHighPass = FilterType::New();
overSampledHighPass->SetInput(cropedHighPass);
overSampledHighPass->SetSubsampleFactor(delta);
overSampledHighPass->Update();
InnerProductType innerProduct;
itk::ImageRegionIterator<OutputImageType> out
(outputImage, overSampledLowPass->GetOutput()->GetRequestedRegion());
NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(),
overSampledLowPass->GetOutput(),
overSampledLowPass->GetOutput()->GetRequestedRegion());
itk::PeriodicBoundaryCondition<OutputImageType> boundaryCondition;
lowIter.OverrideBoundaryCondition(&boundaryCondition);
NeighborhoodIteratorType highIter(highPassOperator.GetRadius(),
overSampledHighPass->GetOutput(),
overSampledHighPass->GetOutput()->GetRequestedRegion());
highIter.OverrideBoundaryCondition(&boundaryCondition);
lowIter.GoToBegin();
highIter.GoToBegin();
out.GoToBegin();
while (!out.IsAtEnd())
{
out.Set(innerProduct(lowIter, lowPassOperator)
+ innerProduct(highIter, highPassOperator));
++lowIter;
++highIter;
++out;
reporter.CompletedPixel();
}
} // end for each overSampledLowPass/overSampledhighPass pair of entry
}
else // multiscale case
{
for (unsigned int i = 0; i < this->GetNumberOfInputs(); i += 2)
{
InputImagePointerType imgLowPass = const_cast<InputImageType*>(this->GetInput(i));
InputImagePointerType imgHighPass = const_cast<InputImageType*>(this->GetInput(i + 1));
OutputImagePointerType outputImage = this->GetOutput();
if (dir != InputImageDimension - 1)
{
outputImage = m_InternalImages[0][i / 2];
}
InnerProductType innerProduct;
itk::ImageRegionIterator<OutputImageType> out(outputImage, imgLowPass->GetRequestedRegion());
NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(), imgLowPass, imgLowPass->GetRequestedRegion());
itk::PeriodicBoundaryCondition<OutputImageType> boundaryCondition;
lowIter.OverrideBoundaryCondition(&boundaryCondition);
NeighborhoodIteratorType highIter(highPassOperator.GetRadius(), imgHighPass, imgLowPass->GetRequestedRegion());
highIter.OverrideBoundaryCondition(&boundaryCondition);
lowIter.GoToBegin();
highIter.GoToBegin();
out.GoToBegin();
while (!out.IsAtEnd())
{
out.Set((innerProduct(lowIter, lowPassOperator)
+ innerProduct(highIter, highPassOperator)) / 2.);
++lowIter;
++highIter;
++out;
reporter.CompletedPixel();
}
} // end for each imgLowPass/imghighPass pair of entry
} // end multiscale case
if (dir != InputImageDimension - 1)
{
// Note that outputImageRegion correspond to the actual region of (local) input !
OutputImageRegionType outputImageRegion;
this->CallCopyInputRegionToOutputRegion(dir, outputImageRegion, inputRegionForThread);
ThreadedGenerateDataAtDimensionN(dir + 1, reporter, outputImageRegion, threadId);
}
}
template <class TInputImage, class TOutputImage, class TWaveletOperator>
void
WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
::ThreadedGenerateDataAtDimensionN
(unsigned int direction,
itk::ProgressReporter& reporter,
const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
{
OutputImageRegionType outputImageRegion;
this->CallCopyInputRegionToOutputRegion(direction, outputImageRegion, outputRegionForThread);
// Low pass part calculation
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(direction);
lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
lowPassOperator.CreateDirectional();
// High pass part calculation
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(direction);
highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
highPassOperator.CreateDirectional();
// typedef for the iterations over the input image
typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
typedef itk::NeighborhoodInnerProduct<OutputImageType> InnerProductType;
// Faces iterations
typename NeighborhoodIteratorType::RadiusType radiusMax;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
radiusMax[i] = lowPassOperator.GetRadius(i);
if (radiusMax[i] < highPassOperator.GetRadius(i)) radiusMax[i] = highPassOperator.GetRadius(i);
}
// The multiresolution case requires a SubsampleImageFilter step
if (m_SubsampleImageFactor > 1)
{
for (unsigned int i = 0; i < m_InternalImages[direction - 1].size(); i += 2)
{
InputImagePointerType imgLowPass = m_InternalImages[direction - 1][i];
InputImagePointerType imgHighPass = m_InternalImages[direction - 1][i + 1];
OutputImagePointerType outputImage = this->GetOutput();
if (direction < InputImageDimension - 1)
{
outputImage = m_InternalImages[direction][i / 2];
}
typedef itk::SubsampleImageFilter<OutputImageType, OutputImageType, Wavelet::INVERSE> FilterType;
typename FilterType::InputImageIndexType delta;
delta.Fill(1);
delta[direction] = this->GetSubsampleImageFactor();
InputImagePointerType cropedLowPass = InputImageType::New();
cropedLowPass->SetRegions(outputRegionForThread);
cropedLowPass->Allocate();
cropedLowPass->FillBuffer(0.);
itk::ImageRegionIterator<InputImageType> cropedLowPassIt(cropedLowPass, outputRegionForThread);
itk::ImageRegionIterator<InputImageType> imgLowPassIt(imgLowPass, outputRegionForThread);
for (cropedLowPassIt.GoToBegin(), imgLowPassIt.GoToBegin();
!cropedLowPassIt.IsAtEnd() && !imgLowPassIt.IsAtEnd();
++cropedLowPassIt, ++imgLowPassIt)
{
cropedLowPassIt.Set(imgLowPassIt.Get());
}
typename FilterType::Pointer overSampledLowPass = FilterType::New();
overSampledLowPass->SetInput(cropedLowPass);
overSampledLowPass->SetSubsampleFactor(delta);
overSampledLowPass->Update();
InputImagePointerType cropedHighPass = InputImageType::New();
cropedHighPass->SetRegions(outputRegionForThread);
cropedHighPass->Allocate();
cropedHighPass->FillBuffer(0.);
itk::ImageRegionIterator<InputImageType> cropedHighPassIt(cropedHighPass, outputRegionForThread);
itk::ImageRegionIterator<InputImageType> imgHighPassIt(imgHighPass, outputRegionForThread);
for (cropedHighPassIt.GoToBegin(), imgHighPassIt.GoToBegin();
!cropedHighPassIt.IsAtEnd() && !imgHighPassIt.IsAtEnd();
++cropedHighPassIt, ++imgHighPassIt)
{
cropedHighPassIt.Set(imgHighPassIt.Get());
}
typename FilterType::Pointer overSampledHighPass = FilterType::New();
overSampledHighPass->SetInput(cropedHighPass);
overSampledHighPass->SetSubsampleFactor(delta);
overSampledHighPass->Update();
InnerProductType innerProduct;
itk::ImageRegionIterator<OutputImageType> out(outputImage,
overSampledLowPass->GetOutput()->GetRequestedRegion());
// TODO: This might be the cause of the multithreading bug : we use a neighborhood iterator on cropped data
// Are we sure that we have cropped enough data to access the neighborhood ?
NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(),
overSampledLowPass->GetOutput(),
overSampledLowPass->GetOutput()->GetRequestedRegion());
itk::PeriodicBoundaryCondition<OutputImageType> boundaryCondition;
lowIter.OverrideBoundaryCondition(&boundaryCondition);
NeighborhoodIteratorType highIter(highPassOperator.GetRadius(),
overSampledHighPass->GetOutput(),
overSampledHighPass->GetOutput()->GetRequestedRegion());
highIter.OverrideBoundaryCondition(&boundaryCondition);
lowIter.GoToBegin();
highIter.GoToBegin();
out.GoToBegin();
while (!out.IsAtEnd())
{
out.Set(innerProduct(lowIter, lowPassOperator)
+ innerProduct(highIter, highPassOperator));
++lowIter;
++highIter;
++out;
reporter.CompletedPixel();
}
} // end for each overSampledLowPass/overSampledhighPass pair of entry
}
else // multiscale case
{
for (unsigned int i = 0; i < m_InternalImages[direction - 1].size(); i += 2)
{
InputImagePointerType imgLowPass = m_InternalImages[direction - 1][i];
InputImagePointerType imgHighPass = m_InternalImages[direction - 1][i + 1];
OutputImagePointerType outputImage = this->GetOutput();
if (direction < InputImageDimension - 1)
{
outputImage = m_InternalImages[direction][i / 2];
}
InnerProductType innerProduct;
itk::ImageRegionIterator<OutputImageType> out(outputImage, imgLowPass->GetRequestedRegion());
NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(), imgLowPass, imgLowPass->GetRequestedRegion());
itk::PeriodicBoundaryCondition<OutputImageType> boundaryCondition;
lowIter.OverrideBoundaryCondition(&boundaryCondition);
NeighborhoodIteratorType highIter(highPassOperator.GetRadius(), imgHighPass, imgLowPass->GetRequestedRegion());
highIter.OverrideBoundaryCondition(&boundaryCondition);
lowIter.GoToBegin();
highIter.GoToBegin();
out.GoToBegin();
while (!out.IsAtEnd())
{
out.Set((innerProduct(lowIter, lowPassOperator)
+ innerProduct(highIter, highPassOperator)) / 2.);
++lowIter;
++highIter;
++out;
reporter.CompletedPixel();
}
} // end for each imgLowPass/imghighPass pair of entry
}
if (direction < InputImageDimension - 1)
{
ThreadedGenerateDataAtDimensionN(direction + 1, reporter, outputImageRegion, threadId);
}
}
} // end of namespace
#endif
