ANIMA  4.0
animaJacobianMatrixImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionConstIteratorWithIndex.h>
5 #include <itkImageRegionIterator.h>
6 
7 #include <vnl_qr.h>
8 
9 namespace anima
10 {
11 
12 template <typename TPixelType, typename TOutputPixelType, unsigned int Dimension>
13 void
14 JacobianMatrixImageFilter <TPixelType, TOutputPixelType, Dimension>
15 ::BeforeThreadedGenerateData()
16 {
17  this->Superclass::BeforeThreadedGenerateData();
18 
19  m_FieldInterpolator = InterpolatorType::New();
20  m_FieldInterpolator->SetInputImage(this->GetInput());
21 
22  if (m_ComputeDeterminant)
23  {
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());
30 
31  m_DeterminantImage->Allocate();
32  m_DeterminantImage->FillBuffer(0.0);
33  }
34 }
35 
36 template <typename TPixelType, typename TOutputPixelType, unsigned int Dimension>
37 void
39 ::SetNeighborhood(unsigned int val)
40 {
41  if (val == 0)
42  {
43  m_Neighborhood = 1;
44  m_OnlySixConnectivity = true;
45  }
46  else
47  {
48  m_Neighborhood = val;
49  m_OnlySixConnectivity = false;
50  }
51 }
52 
53 template <typename TPixelType, typename TOutputPixelType, unsigned int Dimension>
54 void
56 ::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread)
57 {
58  typedef itk::ImageRegionConstIteratorWithIndex <InputImageType> InputIteratorType;
59  typedef itk::ImageRegionIterator <OutputImageType> OutputIteratorType;
60  typedef itk::ImageRegionIterator <DeterminantImageType> DetIteratorType;
61 
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);
67 
68  IndexType currentIndex;
69  IndexType internalIndex, internalIndexOpposite;
70  vnl_matrix <double> dataMatrix, dataMatrixToSolve;
71  vnl_vector <double> dataVector, dataVectorToSolve, outputVector;
72  RegionType inputRegion;
73  RegionType largestRegion = this->GetInput()->GetLargestPossibleRegion();
74  InputPixelType pixelBeforeValue, pixelAfterValue;
75  OutputPixelType outputValue;
76  PointType pointBefore, pointAfter, unitVector;
77  std::vector <IndexType> indexesUsedVector;
78  vnl_matrix <double> jacMatrix(Dimension,Dimension);
79 
80  double meanSpacing = 0;
81  for (unsigned int i = 0;i < Dimension;++i)
82  meanSpacing += this->GetInput()->GetSpacing()[i];
83  meanSpacing /= Dimension;
84 
85  while (!inputItr.IsAtEnd())
86  {
87  currentIndex = inputItr.GetIndex();
88 
89  // Build the explored region
90  for (unsigned int i = 0;i < Dimension;++i)
91  {
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));
97 
98  inputRegion.SetIndex(i,baseIndex);
99  inputRegion.SetSize(i,upperIndex - baseIndex + 1);
100  }
101 
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);
106 
107  InputIteratorType internalItr(this->GetInput(),inputRegion);
108  unsigned int pos = 0;
109  indexesUsedVector.clear();
110  while (!internalItr.IsAtEnd())
111  {
112  internalIndex = internalItr.GetIndex();
113  if (internalIndex == currentIndex)
114  {
115  ++internalItr;
116  continue;
117  }
118 
119  if (m_OnlySixConnectivity)
120  {
121  bool faceConnex = this->CheckFaceConnectivity(internalIndex,currentIndex);
122  if (!faceConnex)
123  {
124  ++internalItr;
125  continue;
126  }
127  }
128 
129  for (unsigned int i = 0;i < Dimension;++i)
130  internalIndexOpposite[i] = 2 * currentIndex[i] - internalIndex[i];
131 
132  bool skipVoxel = false;
133  for (unsigned int i = 0;i < indexesUsedVector.size();++i)
134  {
135  if ((internalIndexOpposite == indexesUsedVector[i])||(internalIndex == indexesUsedVector[i]))
136  {
137  skipVoxel = true;
138  break;
139  }
140  }
141 
142  if (skipVoxel)
143  {
144  ++internalItr;
145  continue;
146  }
147 
148  indexesUsedVector.push_back(internalIndex);
149  indexesUsedVector.push_back(internalIndexOpposite);
150 
151  pixelBeforeValue = internalItr.Get();
152  pixelAfterValue = m_FieldInterpolator->EvaluateAtIndex(internalIndexOpposite);
153 
154  this->GetInput()->TransformIndexToPhysicalPoint(internalIndex,pointBefore);
155  this->GetInput()->TransformIndexToPhysicalPoint(internalIndexOpposite,pointAfter);
156 
157  double pointDist = 0;
158  for (unsigned int i = 0;i < Dimension;++i)
159  {
160  unitVector[i] = pointAfter[i] - pointBefore[i];
161  pointDist += unitVector[i] * unitVector[i];
162  }
163 
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;
168 
169  for (unsigned int i = 0;i < Dimension;++i)
170  {
171  double diffValue = (pixelAfterValue[i] - pixelBeforeValue[i]) / pointDist;
172  dataVector[pos * Dimension + i] = diffValue * dataWeight;
173 
174  for (unsigned int j = 0;j < Dimension;++j)
175  dataMatrix(pos * Dimension + i,i * Dimension + j) = unitVector[j] * dataWeight;
176  }
177 
178  ++pos;
179  ++internalItr;
180  }
181 
182  dataMatrixToSolve.set_size(pos * Dimension,Dimension * Dimension);
183  dataVectorToSolve.set_size(pos * Dimension);
184  for (unsigned int i = 0;i < pos * Dimension;++i)
185  {
186  dataVectorToSolve[i] = dataVector[i];
187  for (unsigned int j = 0;j < Dimension * Dimension;++j)
188  dataMatrixToSolve(i,j) = dataMatrix(i,j);
189  }
190 
191  outputVector = vnl_qr <double> (dataMatrixToSolve).solve(dataVectorToSolve);
192 
193  for (unsigned int i = 0;i < outputValue.GetVectorDimension();++i)
194  outputValue[i] = outputVector[i];
195 
196  if (!m_NoIdentity)
197  {
198  for (unsigned int i = 0;i < Dimension;++i)
199  outputValue[i * (Dimension + 1)] += 1.0;
200  }
201 
202  outputItr.Set(outputValue);
203 
204  if (m_ComputeDeterminant)
205  {
206  for (unsigned int i = 0;i < Dimension;++i)
207  {
208  for (unsigned int j = 0;j < Dimension;++j)
209  jacMatrix(i,j) = outputValue[i * Dimension + j];
210  }
211 
212  detImageIt.Set(vnl_determinant (jacMatrix));
213  }
214 
215  ++inputItr;
216  ++outputItr;
217  if (m_ComputeDeterminant)
218  ++detImageIt;
219  }
220 }
221 
222 template <typename TPixelType, typename TOutputPixelType, unsigned int Dimension>
223 bool
225 ::CheckFaceConnectivity(const IndexType &internalIndex, const IndexType &currentIndex)
226 {
227  unsigned int numDiff = 0;
228  for (unsigned int i = 0;i < Dimension;++i)
229  {
230  if (internalIndex[i] != currentIndex[i])
231  ++numDiff;
232 
233  if (numDiff > 1)
234  return false;
235  }
236 
237  return true;
238 }
239 
240 } //end of namespace anima
Compute the Jacobian matrix in real coordinates of a displacement field.
Superclass::OutputImageRegionType OutputImageRegionType