22 #ifndef __StreamingLargeFeatherMosaicFilter_hxx
23 #define __StreamingLargeFeatherMosaicFilter_hxx
33 template <
class TInputImage,
class TOutputImage,
class TDistanceImage,
class TInternalValueType>
39 itkDebugMacro(<<
"Actually executing thread " << threadId <<
" in region " << outputRegionForThread);
42 itk::ProgressReporter progress(
this, threadId, outputRegionForThread.GetNumberOfPixels());
48 const unsigned int nbOfUsedInputImages = Superclass::GetNumberOfUsedInputImages();
51 const unsigned int nBands = Superclass::GetNumberOfBands();
54 IteratorType outputIt(mosaicImage, outputRegionForThread);
57 typename std::vector<InputImageType*> currentImage;
58 typename std::vector<InterpolatorPointerType> interp;
59 Superclass::PrepareImageAccessors(currentImage, interp);
62 typename std::vector<DistanceImageType*> currentDistanceImage;
63 typename std::vector<DistanceImageInterpolatorPointer> distanceInterpolator;
64 Superclass::PrepareDistanceImageAccessors(currentDistanceImage, distanceInterpolator);
71 interpolatedMathPixel.SetSize(nBands);
72 tempOutputPixel.SetSize(nBands);
76 bool isDataInCurrentOutputPixel;
83 for (outputIt.GoToBegin(); !outputIt.IsAtEnd(); ++outputIt)
86 mosaicImage->TransformIndexToPhysicalPoint(outputIt.GetIndex(), geoPoint);
89 isDataInCurrentOutputPixel =
false;
93 tempOutputPixel.Fill(0.0);
96 for (i = 0; i < nbOfUsedInputImages; i++)
102 if (interp[i]->IsInsideBuffer(geoPoint))
105 interpolatedPixel = interp[i]->Evaluate(geoPoint);
108 if (Superclass::IsPixelNotEmpty(interpolatedPixel))
111 if (distanceInterpolator[i]->IsInsideBuffer(geoPoint))
113 distanceImagePixel = distanceInterpolator[i]->Evaluate(geoPoint);
114 distanceImagePixel -= Superclass::GetDistanceOffset();
116 if (distanceImagePixel > 0)
119 isDataInCurrentOutputPixel =
true;
122 sumDistances += distanceImagePixel;
129 const unsigned int inputImageIndex = Superclass::GetUsedInputImageIndice(i);
130 for (band = 0; band < nBands; band++)
133 interpolatedMathPixel[band] =
static_cast<InternalValueType>(interpolatedPixel[band]);
136 if (this->GetShiftScaleInputImages())
138 this->ShiftScaleValue(interpolatedMathPixel[band], inputImageIndex, band);
142 interpolatedMathPixel[band] *= distanceImagePixel;
145 tempOutputPixel[band] += interpolatedMathPixel[band];
153 itkWarningMacro(<<
"Unable to evaluate distance at point " << geoPoint);
162 if (isDataInCurrentOutputPixel)
166 for (band = 0; band < nBands; band++)
170 pixelValue = tempOutputPixel[band] / sumDistances;
171 Superclass::NormalizePixelValue(pixelValue);
178 outputIt.Set(outputPixel);
181 progress.CompletedPixel();