4 #include <itkImageRegionConstIteratorWithIndex.h> 5 #include <itkImageRegionIterator.h> 12 template <
typename TPixelType,
typename TOutputPixelType,
unsigned int Dimension>
14 JacobianMatrixImageFilter <TPixelType, TOutputPixelType, Dimension>
15 ::BeforeThreadedGenerateData()
17 this->Superclass::BeforeThreadedGenerateData();
19 m_FieldInterpolator = InterpolatorType::New();
20 m_FieldInterpolator->SetInputImage(this->GetInput());
22 if (m_ComputeDeterminant)
24 m_DeterminantImage = DeterminantImageType::New();
25 m_DeterminantImage->Initialize();
26 m_DeterminantImage->SetOrigin(this->GetOutput()->GetOrigin());
27 m_DeterminantImage->SetSpacing(this->GetOutput()->GetSpacing());
28 m_DeterminantImage->SetDirection(this->GetOutput()->GetDirection());
29 m_DeterminantImage->SetRegions(this->GetOutput()->GetLargestPossibleRegion());
31 m_DeterminantImage->Allocate();
32 m_DeterminantImage->FillBuffer(0.0);
36 template <
typename TPixelType,
typename TOutputPixelType,
unsigned int Dimension>
39 ::SetNeighborhood(
unsigned int val)
44 m_OnlySixConnectivity =
true;
49 m_OnlySixConnectivity =
false;
53 template <
typename TPixelType,
typename TOutputPixelType,
unsigned int Dimension>
58 typedef itk::ImageRegionConstIteratorWithIndex <InputImageType> InputIteratorType;
59 typedef itk::ImageRegionIterator <OutputImageType> OutputIteratorType;
60 typedef itk::ImageRegionIterator <DeterminantImageType> DetIteratorType;
62 InputIteratorType inputItr(this->GetInput(),outputRegionForThread);
63 OutputIteratorType outputItr(this->GetOutput(),outputRegionForThread);
64 DetIteratorType detImageIt;
65 if (m_ComputeDeterminant)
66 detImageIt = DetIteratorType (m_DeterminantImage,outputRegionForThread);
69 IndexType internalIndex, internalIndexOpposite;
70 vnl_matrix <double> dataMatrix, dataMatrixToSolve;
71 vnl_vector <double> dataVector, dataVectorToSolve, outputVector;
73 RegionType largestRegion = this->GetInput()->GetLargestPossibleRegion();
76 PointType pointBefore, pointAfter, unitVector;
77 std::vector <IndexType> indexesUsedVector;
78 vnl_matrix <double> jacMatrix(Dimension,Dimension);
80 double meanSpacing = 0;
81 for (
unsigned int i = 0;i < Dimension;++i)
82 meanSpacing += this->GetInput()->GetSpacing()[i];
83 meanSpacing /= Dimension;
85 while (!inputItr.IsAtEnd())
87 currentIndex = inputItr.GetIndex();
90 for (
unsigned int i = 0;i < Dimension;++i)
92 int largestIndex = largestRegion.GetIndex()[i];
93 int testedIndex = currentIndex[i] - m_Neighborhood;
94 unsigned int baseIndex = std::max(largestIndex,testedIndex);
95 unsigned int maxValue = largestIndex + largestRegion.GetSize()[i] - 1;
96 unsigned int upperIndex = std::min(maxValue,(
unsigned int)(currentIndex[i] + m_Neighborhood));
98 inputRegion.SetIndex(i,baseIndex);
99 inputRegion.SetSize(i,upperIndex - baseIndex + 1);
102 unsigned int maxRowSize = (inputRegion.GetNumberOfPixels() - 1) * Dimension;
103 dataMatrix.set_size(maxRowSize,Dimension * Dimension);
104 dataMatrix.fill(0.0);
105 dataVector.set_size(maxRowSize);
107 InputIteratorType internalItr(this->GetInput(),inputRegion);
108 unsigned int pos = 0;
109 indexesUsedVector.clear();
110 while (!internalItr.IsAtEnd())
112 internalIndex = internalItr.GetIndex();
113 if (internalIndex == currentIndex)
119 if (m_OnlySixConnectivity)
121 bool faceConnex = this->CheckFaceConnectivity(internalIndex,currentIndex);
129 for (
unsigned int i = 0;i < Dimension;++i)
130 internalIndexOpposite[i] = 2 * currentIndex[i] - internalIndex[i];
132 bool skipVoxel =
false;
133 for (
unsigned int i = 0;i < indexesUsedVector.size();++i)
135 if ((internalIndexOpposite == indexesUsedVector[i])||(internalIndex == indexesUsedVector[i]))
148 indexesUsedVector.push_back(internalIndex);
149 indexesUsedVector.push_back(internalIndexOpposite);
151 pixelBeforeValue = internalItr.Get();
152 pixelAfterValue = m_FieldInterpolator->EvaluateAtIndex(internalIndexOpposite);
154 this->GetInput()->TransformIndexToPhysicalPoint(internalIndex,pointBefore);
155 this->GetInput()->TransformIndexToPhysicalPoint(internalIndexOpposite,pointAfter);
157 double pointDist = 0;
158 for (
unsigned int i = 0;i < Dimension;++i)
160 unitVector[i] = pointAfter[i] - pointBefore[i];
161 pointDist += unitVector[i] * unitVector[i];
164 double dataWeight = std::sqrt(std::exp(- pointDist / (meanSpacing * meanSpacing)));
165 pointDist = std::sqrt(pointDist);
166 for (
unsigned int i = 0;i < Dimension;++i)
167 unitVector[i] /= pointDist;
169 for (
unsigned int i = 0;i < Dimension;++i)
171 double diffValue = (pixelAfterValue[i] - pixelBeforeValue[i]) / pointDist;
172 dataVector[pos * Dimension + i] = diffValue * dataWeight;
174 for (
unsigned int j = 0;j < Dimension;++j)
175 dataMatrix(pos * Dimension + i,i * Dimension + j) = unitVector[j] * dataWeight;
182 dataMatrixToSolve.set_size(pos * Dimension,Dimension * Dimension);
183 dataVectorToSolve.set_size(pos * Dimension);
184 for (
unsigned int i = 0;i < pos * Dimension;++i)
186 dataVectorToSolve[i] = dataVector[i];
187 for (
unsigned int j = 0;j < Dimension * Dimension;++j)
188 dataMatrixToSolve(i,j) = dataMatrix(i,j);
191 outputVector = vnl_qr <double> (dataMatrixToSolve).solve(dataVectorToSolve);
193 for (
unsigned int i = 0;i < outputValue.GetVectorDimension();++i)
194 outputValue[i] = outputVector[i];
198 for (
unsigned int i = 0;i < Dimension;++i)
199 outputValue[i * (Dimension + 1)] += 1.0;
202 outputItr.Set(outputValue);
204 if (m_ComputeDeterminant)
206 for (
unsigned int i = 0;i < Dimension;++i)
208 for (
unsigned int j = 0;j < Dimension;++j)
209 jacMatrix(i,j) = outputValue[i * Dimension + j];
212 detImageIt.Set(vnl_determinant (jacMatrix));
217 if (m_ComputeDeterminant)
222 template <
typename TPixelType,
typename TOutputPixelType,
unsigned int Dimension>
227 unsigned int numDiff = 0;
228 for (
unsigned int i = 0;i < Dimension;++i)
230 if (internalIndex[i] != currentIndex[i])
InputImageType::RegionType RegionType
InputImageType::IndexType IndexType
Compute the Jacobian matrix in real coordinates of a displacement field.
InputImageType::PixelType InputPixelType
InputImageType::PointType PointType
Superclass::OutputImageRegionType OutputImageRegionType
OutputImageType::PixelType OutputPixelType