4 #include <itkImageRegionConstIteratorWithIndex.h> 5 #include <itkImageRegionIteratorWithIndex.h> 8 #include <itkTranslationTransform.h> 9 #include <itkMatrixOffsetTransformBase.h> 10 #include <itkCompositeTransform.h> 17 template <
typename TImageType,
typename TInterpolatorPrecisionType>
24 typename InterpolatorType::Pointer tmpInterpolator = InterpolatorType::New();
25 this->SetInterpolator(tmpInterpolator.GetPointer());
29 template <
typename TImageType,
typename TInterpolatorPrecisionType>
34 Superclass::GenerateInputRequestedRegion();
35 if (!this->GetInput())
39 inputPtr->SetRequestedRegionToLargestPossibleRegion();
42 template <
typename TImageType,
typename TInterpolatorPrecisionType>
47 Superclass::GenerateOutputInformation();
49 this->GetOutput()->SetSpacing(m_OutputSpacing);
50 this->GetOutput()->SetOrigin(m_OutputOrigin);
51 this->GetOutput()->SetDirection(m_OutputDirection);
52 this->GetOutput()->SetRegions(m_OutputLargestPossibleRegion);
53 this->GetOutput()->SetNumberOfComponentsPerPixel(this->GetOutputVectorLength());
56 template <
typename TImageType,
typename TInterpolatorPrecisionType>
61 if (this->GetNumberOfIndexedInputs() > 0)
62 return this->GetInput(0)->GetNumberOfComponentsPerPixel();
67 template <
typename TImageType,
typename TInterpolatorPrecisionType>
72 Superclass::BeforeThreadedGenerateData();
74 if (m_Transform.IsNull())
75 itkExceptionMacro(
"No valid transformation...");
77 if (m_Interpolator.IsNull())
78 this->InitializeInterpolator();
80 m_Interpolator->SetInputImage(this->GetInput(0));
82 if (!m_Transform->IsLinear())
84 m_StartIndDef = m_OutputLargestPossibleRegion.GetIndex();
85 m_EndIndDef = m_StartIndDef + m_OutputLargestPossibleRegion.GetSize();
89 template <
typename TImageType,
typename TInterpolatorPrecisionType>
94 unsigned int threadId = this->GetSafeThreadId();
96 if (m_Transform->IsLinear())
97 this->LinearThreadedGenerateData(outputRegionForThread,threadId);
99 this->NonLinearThreadedGenerateData(outputRegionForThread,threadId);
101 this->SafeReleaseThreadId(threadId);
104 template <
typename TImageType,
typename TInterpolatorPrecisionType>
109 typedef itk::ImageRegionIteratorWithIndex <InputImageType> IteratorType;
111 IteratorType outputItr(this->GetOutput(),outputRegionForThread);
115 unsigned int vectorSize = this->GetOutputVectorLength();
119 vnl_matrix <double> orientationMatrix = this->ComputeLinearJacobianMatrix();
120 vnl_matrix <double> parametersRotationMatrix;
121 this->ComputeRotationParametersFromReorientationMatrix(orientationMatrix,parametersRotationMatrix);
123 bool lastDimensionUseless = (this->GetInput(0)->GetLargestPossibleRegion().GetSize()[ImageDimension - 1] <= 1);
125 while (!outputItr.IsAtEnd())
127 tmpInd = outputItr.GetIndex();
129 this->GetOutput()->TransformIndexToPhysicalPoint(tmpInd,tmpPoint);
131 tmpPoint = m_Transform->TransformPoint(tmpPoint);
133 this->GetInput(0)->TransformPhysicalPointToContinuousIndex(tmpPoint,index);
135 if (lastDimensionUseless)
136 index[ImageDimension - 1] = 0;
138 if (m_Interpolator->IsInsideBuffer(index))
139 tmpRes = m_Interpolator->EvaluateAtContinuousIndex(index);
141 this->InitializeZeroPixel(tmpRes);
145 this->ReorientInterpolatedModel(tmpRes,parametersRotationMatrix,resRotated,threadId);
146 outputItr.Set(resRotated);
149 outputItr.Set(tmpRes);
156 template <
typename TImageType,
typename TInterpolatorPrecisionType>
161 typedef itk::ImageRegionIteratorWithIndex <InputImageType> IteratorType;
163 IteratorType outputItr(this->GetOutput(),outputRegionForThread);
167 unsigned int vectorSize = this->GetOutputVectorLength();
171 vnl_matrix <double> orientationMatrix(ImageDimension,ImageDimension);
172 vnl_matrix <double> parametersRotationMatrix;
174 while (!outputItr.IsAtEnd())
176 tmpInd = outputItr.GetIndex();
177 this->GetOutput()->TransformIndexToPhysicalPoint(tmpInd,tmpPoint);
179 tmpPoint = m_Transform->TransformPoint(tmpPoint);
181 this->GetInput(0)->TransformPhysicalPointToContinuousIndex(tmpPoint,index);
183 if (m_Interpolator->IsInsideBuffer(index))
184 tmpRes = m_Interpolator->EvaluateAtContinuousIndex(index);
186 this->InitializeZeroPixel(tmpRes);
190 this->ComputeLocalJacobianMatrix(tmpInd,orientationMatrix);
191 this->ComputeRotationParametersFromReorientationMatrix(orientationMatrix,parametersRotationMatrix);
192 this->ReorientInterpolatedModel(tmpRes,parametersRotationMatrix,resRotated,threadId);
193 outputItr.Set(resRotated);
196 outputItr.Set(tmpRes);
202 template <
typename TImageType,
typename TInterpolatorPrecisionType>
207 vnl_matrix <double> reorientationMatrix(ImageDimension,ImageDimension);
208 reorientationMatrix.set_identity();
210 typedef itk::TranslationTransform <TInterpolatorPrecisionType,ImageDimension> TranslationType;
211 if (dynamic_cast <TranslationType *> (m_Transform.GetPointer()))
212 return reorientationMatrix;
214 typedef itk::MatrixOffsetTransformBase <TInterpolatorPrecisionType,ImageDimension,ImageDimension> AffineTransformType;
216 AffineTransformType *tmpTrsf = dynamic_cast <AffineTransformType *> (m_Transform.GetPointer());
220 reorientationMatrix = tmpTrsf->GetMatrix().GetVnlMatrix().as_matrix();
221 return reorientationMatrix;
224 typedef itk::CompositeTransform <TInterpolatorPrecisionType,ImageDimension> CompositeTransformType;
225 CompositeTransformType *compositeTrsf = dynamic_cast <CompositeTransformType *> (m_Transform.GetPointer());
227 for (
unsigned int i = 0;i < compositeTrsf->GetNumberOfTransforms();++i)
229 tmpTrsf = dynamic_cast <AffineTransformType *> (compositeTrsf->GetNthTransform(i).GetPointer());
231 reorientationMatrix *= tmpTrsf->GetMatrix().GetVnlMatrix().as_matrix();
234 return reorientationMatrix;
237 template <
typename TImageType,
typename TInterpolatorPrecisionType>
242 vnl_matrix <double> jacMatrix(ImageDimension,ImageDimension);
244 vnl_matrix <double> deltaMatrix(ImageDimension,ImageDimension,0);
248 vnl_matrix <double> resDiff(ImageDimension,ImageDimension,0);
253 for (
unsigned int i = 0;i < ImageDimension;++i)
258 if (posBef[i] < m_StartIndDef[i])
259 posBef[i] = m_StartIndDef[i];
261 outputPtr->TransformIndexToPhysicalPoint(posBef,tmpPosBef);
266 if (posAfter[i] >= m_EndIndDef[i])
267 posAfter[i] = m_EndIndDef[i] - 1;
269 outputPtr->TransformIndexToPhysicalPoint(posAfter,tmpPosAfter);
271 if (posAfter[i] == posBef[i])
273 deltaMatrix(i,i) = 1;
277 double normDelta = 0.0;
278 for (
unsigned int j = 0;j < ImageDimension;++j)
280 deltaMatrix(j,i) = tmpPosAfter[j] - tmpPosBef[j];
281 normDelta += deltaMatrix(j,i) * deltaMatrix(j,i);
284 normDelta = std::sqrt(normDelta);
285 for (
unsigned int j = 0;j < ImageDimension;++j)
286 deltaMatrix(j,i) /= normDelta;
288 tmpPos = m_Transform->TransformPoint(tmpPosAfter);
289 for (
unsigned int j = 0;j < ImageDimension;++j)
290 resDiff(j,i) = tmpPos[j];
292 tmpPos = m_Transform->TransformPoint(tmpPosBef);
293 for (
unsigned int j = 0;j < ImageDimension;++j)
294 resDiff(j,i) = (resDiff(j,i) - tmpPos[j]) / normDelta;
297 vnl_matrix_inverse <double> invDelta(deltaMatrix);
298 deltaMatrix = invDelta.inverse();
300 for (
unsigned int i = 0;i < ImageDimension;++i)
302 for (
unsigned int j = 0;j < ImageDimension;++j)
305 for (
unsigned int k = 0;k < ImageDimension;++k)
306 jacMatrix(i,j) += resDiff(i,k) * deltaMatrix(k,j);
310 if (m_StartIndDef[2] == (m_EndIndDef[2]-1))
313 reorientationMatrix = jacMatrix;
316 template <
typename TImageType,
typename TInterpolatorPrecisionType>
320 vnl_matrix <double> &modelOrientationMatrix)
322 if (m_FiniteStrainReorientation)
324 vnl_matrix <double> tmpMat(ImageDimension,ImageDimension);
328 modelOrientationMatrix = reorientationMatrix;
virtual unsigned int GetOutputVectorLength()
virtual void GenerateOutputInformation() ITK_OVERRIDE
virtual void GenerateInputRequestedRegion() ITK_OVERRIDE
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE
InterpolatorType::PointType PointType
void ComputeLocalJacobianMatrix(InputIndexType &index, vnl_matrix< double > &reorientationMatrix)
virtual void InitializeInterpolator()
Initializes the default interpolator, might change in derived classes.
virtual void BeforeThreadedGenerateData() ITK_OVERRIDE
void LinearThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, unsigned int threadId)
vnl_matrix< double > ComputeLinearJacobianMatrix()
itk::InterpolateImageFunction< InputImageType, TInterpolatorPrecisionType > InterpolatorType
void ExtractRotationFromJacobianMatrix(vnl_matrix< RealType > &jacobianMatrix, vnl_matrix< RealType > &rotationMatrix, vnl_matrix< RealType > &tmpMat)
TOutputImage::Pointer OutputImagePointer
InputImageType::PixelType InputPixelType
InterpolatorType::ContinuousIndexType ContinuousIndexType
InputImageType::Pointer InputImagePointer
void NonLinearThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, unsigned int threadId)
TOutputImage::RegionType OutputImageRegionType
InputImageType::IndexType InputIndexType
virtual void ComputeRotationParametersFromReorientationMatrix(vnl_matrix< double > &reorientationMatrix, vnl_matrix< double > &modelOrientationMatrix)