- //=========================================================================================================================
- //update the output for the outputRegionForThread
- template<class InputImageType, class OutputImageType>
- void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
- {
-
- //Get pointer to the input
- typename InputImageType::ConstPointer inputPtr = this->GetInput();
-
- //Get pointer to the output
- typename OutputImageType::Pointer outputPtr = this->GetOutput();
- typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
-
- //Iterator over input
- typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
-
- //define them over the outputRegionForThread
- InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
-
- //Initialize
- typename InputImageType::IndexType index;
- itk::ContinuousIndex<double,ImageDimension> contIndex;
- typename InputImageType::PointType ipoint;
- typename OutputImageType::PointType opoint;
- typedef typename OutputImageType::PixelType DisplacementType;
- DisplacementType displacement;
- inputIt.GoToBegin();
-
- //define some temp variables
- signed long baseIndex[ImageDimension];
- double distance[ImageDimension];
- unsigned int dim, counter, upper;
- double totalOverlap,overlap;
- typename OutputImageType::IndexType neighIndex;
-
- //Find the number of neighbors
- unsigned int neighbors = 1 << ImageDimension;
-
- //==================================================================================================
- //Loop over the region and add the intensities to the output and the weight to the weights
- //==================================================================================================
- while( !inputIt.IsAtEnd() )
- {
- // get the input image index
- index = inputIt.GetIndex();
- inputPtr->TransformIndexToPhysicalPoint( index,ipoint );
-
- // get the required displacement
- displacement = inputIt.Get();
-
- // compute the required output image point
- for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j];
-
- // Update the output and the weights
- if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) )
- {
- for(dim = 0; dim < ImageDimension; dim++)
- {
- // The following block is equivalent to the following line without
- // having to call floor. (Only for positive inputs, we already now that is in the image)
- // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
-
- baseIndex[dim] = (long) contIndex[dim];
- distance[dim] = contIndex[dim] - double( baseIndex[dim] );
- }
-
- //Add contribution for each neighbor
- totalOverlap = itk::NumericTraits<double>::Zero;
- for( counter = 0; counter < neighbors ; counter++ )
- {
- overlap = 1.0; // fraction overlap
- upper = counter; // each bit indicates upper/lower neighbour
-
- // get neighbor index and overlap fraction
- for( dim = 0; dim < 3; dim++ )
- {
- if ( upper & 1 )
- {
- neighIndex[dim] = baseIndex[dim] + 1;
- overlap *= distance[dim];
- }
- else
- {
- neighIndex[dim] = baseIndex[dim];
- overlap *= 1.0 - distance[dim];
- }
- upper >>= 1;
- }
-
-
-
- //Set neighbor value only if overlap is not zero
- if( (overlap>0.0)) // &&
- // (static_cast<unsigned int>(neighIndex[0])<size[0]) &&
- // (static_cast<unsigned int>(neighIndex[1])<size[1]) &&
- // (static_cast<unsigned int>(neighIndex[2])<size[2]) &&
- // (neighIndex[0]>=0) &&
- // (neighIndex[1]>=0) &&
- // (neighIndex[2]>=0) )
- {
- //what to store? the original displacement vector?
- if (! m_ThreadSafe)
- {
- //Set the pixel and weight at neighIndex
- outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
- m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
- }
-
- else
- {
- //Entering critilal section: shared memory
- m_MutexImage->GetPixel(neighIndex).Lock();
-
- //Set the pixel and weight at neighIndex
- outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
- m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
-
- //Unlock
- m_MutexImage->GetPixel(neighIndex).Unlock();
-
- }
- //Add to total overlap
- totalOverlap += overlap;
- }
-
- if( totalOverlap == 1.0 )
- {
- // finished
- break;
- }
- }
- }
-
- ++inputIt;
+ /** Typedef to describe the output image region type. */
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+protected:
+ HelperClass1();
+ ~HelperClass1() {};
+
+ //the actual processing
+ void BeforeThreadedGenerateData();
+ void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
+
+ //member data
+ typename WeightsImageType::Pointer m_Weights;
+ typename MutexImageType::Pointer m_MutexImage;
+ bool m_ThreadSafe;
+
+};
+
+
+
+//=========================================================================================================================
+//Member functions of the helper class 1
+//=========================================================================================================================
+
+
+//=========================================================================================================================
+//Empty constructor
+template<class InputImageType, class OutputImageType >
+HelperClass1<InputImageType, OutputImageType>::HelperClass1()
+{
+ m_ThreadSafe=false;
+}
+
+//=========================================================================================================================
+//Before threaded data
+template<class InputImageType, class OutputImageType >
+void HelperClass1<InputImageType, OutputImageType>::BeforeThreadedGenerateData()
+{
+ //Since we will add, put to zero!
+ this->GetOutput()->FillBuffer(itk::NumericTraits<double>::Zero);
+ this->GetWeights()->FillBuffer(itk::NumericTraits<double>::Zero);
+}
+
+//=========================================================================================================================
+//update the output for the outputRegionForThread
+template<class InputImageType, class OutputImageType>
+void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
+{
+
+ //Get pointer to the input
+ typename InputImageType::ConstPointer inputPtr = this->GetInput();
+
+ //Get pointer to the output
+ typename OutputImageType::Pointer outputPtr = this->GetOutput();
+ typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
+
+ //Iterator over input
+ typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
+
+ //define them over the outputRegionForThread
+ InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
+
+ //Initialize
+ typename InputImageType::IndexType index;
+ itk::ContinuousIndex<double,ImageDimension> contIndex;
+ typename InputImageType::PointType ipoint;
+ typename OutputImageType::PointType opoint;
+ typedef typename OutputImageType::PixelType DisplacementType;
+ DisplacementType displacement;
+ inputIt.GoToBegin();
+
+ //define some temp variables
+ signed long baseIndex[ImageDimension];
+ double distance[ImageDimension];
+ unsigned int dim, counter, upper;
+ double totalOverlap,overlap;
+ typename OutputImageType::IndexType neighIndex;
+
+ //Find the number of neighbors
+ unsigned int neighbors = 1 << ImageDimension;
+
+ //==================================================================================================
+ //Loop over the region and add the intensities to the output and the weight to the weights
+ //==================================================================================================
+ while( !inputIt.IsAtEnd() ) {
+ // get the input image index
+ index = inputIt.GetIndex();
+ inputPtr->TransformIndexToPhysicalPoint( index,ipoint );
+
+ // get the required displacement
+ displacement = inputIt.Get();
+
+ // compute the required output image point
+ for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j];
+
+ // Update the output and the weights
+ if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) ) {
+ for(dim = 0; dim < ImageDimension; dim++) {
+ // The following block is equivalent to the following line without
+ // having to call floor. (Only for positive inputs, we already now that is in the image)
+ // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
+
+ baseIndex[dim] = (long) contIndex[dim];
+ distance[dim] = contIndex[dim] - double( baseIndex[dim] );