1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #ifndef itkMIRegistrationFunction_hxx
19 #define itkMIRegistrationFunction_hxx
20
21 #include "itkMIRegistrationFunction.h"
22 #include "itkImageRandomIteratorWithIndex.h"
23 #include "itkMacro.h"
24 #include "itkMath.h"
25 #include "itkNeighborhoodIterator.h"
26
27 #include "vnl/vnl_matrix.h"
28
29 namespace itk
30 {
31
32 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
33 MIRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
MIRegistrationFunction()34 ::MIRegistrationFunction() :
35 m_TimeStep( 1.0 ),
36 m_FixedImageGradientCalculator( GradientCalculatorType::New() ),
37 m_DenominatorThreshold( 1e-9 ),
38 m_IntensityDifferenceThreshold( 0.001 ),
39 m_MetricTotal( 0.0 ),
40 m_NumberOfSamples( 1 ),
41 m_NumberOfBins( 4 ),
42 m_Minnorm( 1.0 ),
43 m_DoInverse( false )
44 {
45 m_FixedImageSpacing.Fill( 1 );
46 m_FixedImageOrigin.Fill( 0 );
47
48 RadiusType r;
49
50 for( unsigned int j = 0; j < ImageDimension; j++ )
51 {
52 r[j] = 2;
53 m_NumberOfSamples *= ( r[j] * 2 + 1 );
54 }
55 this->SetRadius(r);
56
57 this->SetMovingImage(nullptr);
58 this->SetFixedImage(nullptr);
59
60 if( m_DoInverse )
61 {
62 m_MovingImageGradientCalculator = GradientCalculatorType::New();
63 }
64 typename DefaultInterpolatorType::Pointer interp =
65 DefaultInterpolatorType::New();
66
67 m_MovingImageInterpolator = static_cast< InterpolatorType * >(
68 interp.GetPointer() );
69 }
70
71 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
72 void
73 MIRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
InitializeIteration()74 ::InitializeIteration()
75 {
76 if( !this->m_MovingImage || !this->m_FixedImage || !m_MovingImageInterpolator )
77 {
78 itkExceptionMacro(<< "MovingImage, FixedImage and/or Interpolator not set");
79 }
80
81 // Set up gradient calculator
82 m_FixedImageGradientCalculator->SetInputImage(this->m_FixedImage);
83
84 if( m_DoInverse )
85 {
86 // Set up gradient calculator
87 m_MovingImageGradientCalculator->SetInputImage(this->m_MovingImage);
88 }
89
90 // Set up moving image interpolator
91 m_MovingImageInterpolator->SetInputImage(this->m_MovingImage);
92
93 m_MetricTotal = 0.0;
94 }
95
96 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
97 typename MIRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
98 ::PixelType
99 MIRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
ComputeUpdate(const NeighborhoodType & it,void * itkNotUsed (globalData),const FloatOffsetType & itkNotUsed (offset))100 ::ComputeUpdate( const NeighborhoodType & it, void *itkNotUsed(globalData),
101 const FloatOffsetType & itkNotUsed(offset) )
102 {
103 // We compute the derivative of MI w.r.t. the infinitesimal
104 // displacement, following Viola and Wells.
105
106 // 1) collect samples from M (Moving) and F (Fixed)
107 // 2) compute minimum and maximum values of M and F
108 // 3) discretized M and F into N bins
109 // 4) estimate joint probability P(M,F) and P(F)
110 // 5) derivatives is given as :
111 //
112 // $$ \nabla MI = \frac{1}{N} \sum_i \sum_j (F_i-F_j)
113 // ( W(F_i,F_j) \frac{1}{\sigma_v} -
114 // W((F_i,M_i),(F_j,M_j)) \frac{1}{\sigma_vv} ) \nabla F
115 //
116 // NOTE : must estimate sigma for each pdf
117
118 using sampleContainerType = std::vector< double >;
119 using gradContainerType = std::vector< CovariantVectorType >;
120 using gradMagContainerType = std::vector< double >;
121 using inImageIndexContainerType = std::vector< unsigned int >;
122
123 PixelType update;
124 PixelType derivative;
125
126 const IndexType oindex = it.GetIndex();
127
128 unsigned int indct;
129
130 for( indct = 0; indct < ImageDimension; indct++ )
131 {
132 update[indct] = 0.0;
133 derivative[indct] = 0.0;
134 }
135
136 float thresh2 = 1.0 / 255.; // FIX ME : FOR PET LUNG ONLY !!
137 float thresh1 = 1.0 / 255.;
138 if( this->m_MovingImage->GetPixel(oindex) <= thresh1
139 && this->m_FixedImage->GetPixel(oindex) <= thresh2 )
140 {
141 return update;
142 }
143
144 typename FixedImageType::SizeType hradius = this->GetRadius();
145
146 auto * img = const_cast< FixedImageType * >( this->m_FixedImage.GetPointer() );
147 typename FixedImageType::SizeType imagesize = img->GetLargestPossibleRegion().GetSize();
148
149 bool inimage;
150
151 // Now collect the samples
152 sampleContainerType fixedSamplesA;
153 sampleContainerType movingSamplesA;
154 sampleContainerType fixedSamplesB;
155 sampleContainerType movingSamplesB;
156 inImageIndexContainerType inImageIndicesA;
157 gradContainerType fixedGradientsA;
158 gradMagContainerType fixedGradMagsA;
159 inImageIndexContainerType inImageIndicesB;
160 gradContainerType fixedGradientsB;
161 gradMagContainerType fixedGradMagsB;
162
163 unsigned int samplestep = 2; //m_Radius[0];
164
165 double minf = 1.e9, minm = 1.e9, maxf = 0.0, maxm = 0.0;
166 double movingMean = 0.0;
167 double fixedMean = 0.0;
168 double fixedValue = 0, movingValue = 0;
169
170 unsigned int sampct = 0;
171
172 ConstNeighborhoodIterator< DisplacementFieldType >
173 asamIt( hradius,
174 this->GetDisplacementField(),
175 this->GetDisplacementField()->GetRequestedRegion() );
176 asamIt.SetLocation(oindex);
177 unsigned int hoodlen = asamIt.Size();
178
179 // First get the density-related sample
180 for( indct = 0; indct < hoodlen; indct = indct + samplestep )
181 {
182 IndexType index = asamIt.GetIndex(indct);
183 inimage = true;
184 for( unsigned int dd = 0; dd < ImageDimension; dd++ )
185 {
186 if( index[dd] < 0 || index[dd] >
187 static_cast< typename IndexType::IndexValueType >( imagesize[dd] - 1 ) ) { inimage = false; }
188 }
189 if( inimage )
190 {
191 fixedValue = 0.;
192 movingValue = 0.0;
193 CovariantVectorType fixedGradient;
194
195 // Get fixed image related information
196 fixedValue = (double)this->m_FixedImage->GetPixel(index);
197
198 fixedGradient = m_FixedImageGradientCalculator->EvaluateAtIndex(index);
199
200 // Get moving image related information
201 using DeformationPixelType = typename DisplacementFieldType::PixelType;
202 const DeformationPixelType itvec = this->GetDisplacementField()->GetPixel(index);
203
204 PointType mappedPoint;
205 this->GetFixedImage()->TransformIndexToPhysicalPoint(index, mappedPoint);
206 for( unsigned int j = 0; j < ImageDimension; j++ )
207 {
208 mappedPoint[j] += itvec[j];
209 }
210 if( m_MovingImageInterpolator->IsInsideBuffer(mappedPoint) )
211 {
212 movingValue = m_MovingImageInterpolator->Evaluate(mappedPoint);
213 }
214 else
215 {
216 movingValue = 0.0;
217 }
218
219 if( fixedValue > maxf )
220 {
221 maxf = fixedValue;
222 }
223 else if ( fixedValue < minf )
224 {
225 minf = fixedValue;
226 }
227
228 if( movingValue > maxm )
229 {
230 maxm = movingValue;
231 }
232 else if( movingValue < minm )
233 {
234 minm = movingValue;
235 }
236
237 fixedMean += fixedValue;
238 movingMean += movingValue;
239
240 fixedSamplesA.insert(fixedSamplesA.begin(), (double)fixedValue);
241 fixedGradientsA.insert(fixedGradientsA.begin(), fixedGradient);
242 movingSamplesA.insert(movingSamplesA.begin(), (double)movingValue);
243
244 // fixedSamplesB.insert(fixedSamplesB.begin(),(double)fixedValue);
245 // fixedGradientsB.insert(fixedGradientsB.begin(),fixedGradient);
246 // movingSamplesB.insert(movingSamplesB.begin(),(double)movingValue);
247
248 sampct++;
249 }
250 }
251
252 // Begin random a samples
253 bool getrandasamples = true;
254 if( getrandasamples )
255 {
256 typename FixedImageType::RegionType region = img->GetLargestPossibleRegion();
257
258 ImageRandomIteratorWithIndex< FixedImageType > randasamit(img, region);
259 unsigned int numberOfSamples = 20;
260 randasamit.SetNumberOfSamples(numberOfSamples);
261
262 indct = 0;
263
264 randasamit.GoToBegin();
265 while( !randasamit.IsAtEnd() && indct < numberOfSamples )
266 {
267 IndexType index = randasamit.GetIndex();
268 inimage = true;
269
270 float d = 0.0;
271 for( unsigned int dd = 0; dd < ImageDimension; dd++ )
272 {
273 if( index[dd] < 0 || index[dd] >
274 static_cast< typename IndexType::IndexValueType >( imagesize[dd] - 1 ) )
275 {
276 inimage = false;
277 }
278 d += ( index[dd] - oindex[dd] ) * ( index[dd] - oindex[dd] );
279 }
280
281 if( inimage )
282 {
283 fixedValue = 0.;
284 movingValue = 0.0;
285 CovariantVectorType fixedGradient;
286 double fgm = 0;
287 // Get fixed image related information
288 fixedValue = (double)this->m_FixedImage->GetPixel(index);
289 fixedGradient = m_FixedImageGradientCalculator->EvaluateAtIndex(index);
290
291 for( unsigned int j = 0; j < ImageDimension; j++ )
292 {
293 fgm += fixedGradient[j] * fixedGradient[j];
294 }
295 // Get moving image related information
296 using DeformationPixelType = typename DisplacementFieldType::PixelType;
297 const DeformationPixelType itvec = this->GetDisplacementField()->GetPixel(index);
298 PointType mappedPoint;
299 this->GetFixedImage()->TransformIndexToPhysicalPoint(index, mappedPoint);
300 for( unsigned int j = 0; j < ImageDimension; j++ )
301 {
302 mappedPoint[j] += itvec[j];
303 }
304 if( m_MovingImageInterpolator->IsInsideBuffer(mappedPoint) )
305 {
306 movingValue = m_MovingImageInterpolator->Evaluate(mappedPoint);
307 }
308 else
309 {
310 movingValue = 0.0;
311 }
312
313 // if ( (fixedValue > 0 || movingValue > 0 || fgm > 0) ||
314 // !filtersamples)
315
316 if( fixedValue > 0 || movingValue > 0 || fgm > 0 )
317 {
318 fixedMean += fixedValue;
319 movingMean += movingValue;
320
321 fixedSamplesA.insert(fixedSamplesA.begin(), (double)fixedValue);
322 fixedGradientsA.insert(fixedGradientsA.begin(), fixedGradient);
323 movingSamplesA.insert(movingSamplesA.begin(), (double)movingValue);
324 sampct++;
325 indct++;
326 }
327 }
328 ++randasamit;
329 }
330 }
331 // End random a samples
332
333 const DisplacementFieldType *const field = this->GetDisplacementField();
334
335 for( unsigned int j = 0; j < ImageDimension; j++ )
336 {
337 hradius[j] = 0;
338 }
339 ConstNeighborhoodIterator< DisplacementFieldType >
340 hoodIt( hradius, field, field->GetRequestedRegion() );
341 hoodIt.SetLocation(oindex);
342
343 // Then get the entropy ( and MI derivative ) related sample
344 for( indct = 0; indct < hoodIt.Size(); indct = indct + 1 )
345 {
346 const IndexType index = hoodIt.GetIndex(indct);
347 inimage = true;
348 float d = 0.0;
349 for( unsigned int dd = 0; dd < ImageDimension; dd++ )
350 {
351 if( index[dd] < 0 || index[dd] >
352 static_cast< typename IndexType::IndexValueType >( imagesize[dd] - 1 ) )
353 {
354 inimage = false;
355 }
356 d += ( index[dd] - oindex[dd] ) * ( index[dd] - oindex[dd] );
357 }
358 if( inimage && std::sqrt(d) <= 1.0 )
359 {
360 fixedValue = 0.;
361 movingValue = 0.0;
362 CovariantVectorType fixedGradient;
363
364 // Get fixed image related information
365 fixedValue = (double)this->m_FixedImage->GetPixel(index);
366 fixedGradient = m_FixedImageGradientCalculator->EvaluateAtIndex(index);
367
368 // Get moving image related information
369 const typename DisplacementFieldType::PixelType hooditvec =
370 this->m_DisplacementField->GetPixel(index);
371 PointType mappedPoint;
372 this->GetFixedImage()->TransformIndexToPhysicalPoint(index, mappedPoint);
373 for( unsigned int j = 0; j < ImageDimension; j++ )
374 {
375 mappedPoint[j] += hooditvec[j];
376 }
377 if ( m_MovingImageInterpolator->IsInsideBuffer(mappedPoint) )
378 {
379 movingValue = m_MovingImageInterpolator->Evaluate(mappedPoint);
380 }
381 else
382 {
383 movingValue = 0.0;
384 }
385
386 fixedSamplesB.insert(fixedSamplesB.begin(), (double)fixedValue);
387 fixedGradientsB.insert(fixedGradientsB.begin(), fixedGradient);
388 movingSamplesB.insert(movingSamplesB.begin(), (double)movingValue);
389 }
390 }
391
392 double fsigma = 0.0;
393 double msigma = 0.0;
394 double jointsigma = 0.0;
395
396 const auto numsamplesB = (double)fixedSamplesB.size();
397 const auto numsamplesA = (double)fixedSamplesA.size();
398 double nsamp = numsamplesB;
399 // if (maxf == minf && maxm == minm) return update;
400 // else std::cout << " b samps " << fixedSamplesB.size()
401 // << " a samps " << fixedSamplesA.size() <<
402 // oindex << hoodIt.Size() << it.Size() << std::endl;
403
404 fixedMean /= (double)sampct;
405 movingMean /= (double)sampct;
406
407 bool mattes = false;
408
409 for( indct = 0; indct < (unsigned int)numsamplesA; indct++ )
410 {
411 // Get fixed image related information
412 fixedValue = fixedSamplesA[indct];
413 movingValue = movingSamplesA[indct];
414
415 fsigma += ( fixedValue - fixedMean ) * ( fixedValue - fixedMean );
416 msigma += ( movingValue - movingMean ) * ( movingValue - movingMean );
417 jointsigma += fsigma + msigma;
418
419 if( mattes )
420 {
421 fixedSamplesA[indct] = fixedSamplesA[indct] - minf;
422 movingSamplesA[indct] = movingSamplesA[indct] - minm;
423 if( indct < numsamplesB )
424 {
425 fixedSamplesB[indct] = fixedSamplesB[indct] - minf;
426 movingSamplesB[indct] = movingSamplesB[indct] - minm;
427 }
428 }
429 }
430
431 fsigma = std::sqrt(fsigma / numsamplesA);
432 float sigmaw = 0.8;
433 double m_FixedImageStandardDeviation = fsigma * sigmaw;
434 msigma = std::sqrt(msigma / numsamplesA);
435 double m_MovingImageStandardDeviation = msigma * sigmaw;
436 jointsigma = std::sqrt(jointsigma / numsamplesA);
437
438 if( fsigma < 1.e-7 || msigma < 1.e-7 )
439 {
440 return update;
441 }
442
443 double m_MinProbability = 0.0001;
444 double dLogSumFixed = 0., dLogSumMoving = 0., dLogSumJoint = 0.0;
445 unsigned int bsamples;
446 unsigned int asamples;
447
448 // the B samples estimate the entropy
449 for ( bsamples = 0; bsamples < (unsigned int)numsamplesB; bsamples++ )
450 {
451 double dDenominatorMoving = m_MinProbability;
452 double dDenominatorJoint = m_MinProbability;
453 double dDenominatorFixed = m_MinProbability;
454 double dSumFixed = m_MinProbability;
455
456 // Estimate the density
457 for( asamples = 0; asamples < (unsigned int)numsamplesA; asamples++ )
458 {
459 double valueFixed = ( fixedSamplesB[bsamples] - fixedSamplesA[asamples] )
460 / m_FixedImageStandardDeviation;
461 valueFixed = std::exp(-0.5 * valueFixed * valueFixed);
462
463 double valueMoving = ( movingSamplesB[bsamples] - movingSamplesA[asamples] )
464 / m_MovingImageStandardDeviation;
465 valueMoving = std::exp(-0.5 * valueMoving * valueMoving);
466
467 dDenominatorMoving += valueMoving;
468 dDenominatorFixed += valueFixed;
469 dSumFixed += valueFixed;
470
471 // Everything above here can be pre-computed only once and stored,
472 // assuming const v.f. in small n-hood
473 dDenominatorJoint += valueMoving * valueFixed;
474 } // end of sample A loop
475
476 dLogSumFixed -= std::log(dSumFixed);
477 dLogSumMoving -= std::log(dDenominatorMoving);
478 dLogSumJoint -= std::log(dDenominatorJoint);
479
480 // Estimate the density
481 for( asamples = 0; asamples < (unsigned int)numsamplesA; asamples++ )
482 {
483 double valueFixed = ( fixedSamplesB[bsamples] - fixedSamplesA[asamples] )
484 / m_FixedImageStandardDeviation;
485 valueFixed = std::exp(-0.5 * valueFixed * valueFixed);
486
487 double valueMoving = ( movingSamplesB[bsamples] - movingSamplesA[asamples] )
488 / m_MovingImageStandardDeviation;
489 valueMoving = std::exp(-0.5 * valueMoving * valueMoving);
490 const double weightFixed = valueFixed / dDenominatorFixed;
491 // dDenominatorJoint and weightJoint are what need to be computed each time
492 const double weightJoint = valueMoving * valueFixed / dDenominatorJoint;
493
494 // Begin where we may switch fixed and moving
495 double weight = ( weightFixed - weightJoint );
496 weight *= ( fixedSamplesB[bsamples] - fixedSamplesA[asamples] );
497 // End where we may switch fixed and moving
498
499 // This can also be stored away
500 for ( unsigned int i = 0; i < ImageDimension; i++ )
501 {
502 derivative[i] += ( fixedGradientsB[bsamples][i] - fixedGradientsA[asamples][i] ) * weight;
503 }
504 } // end of sample A loop
505 } // end of sample B loop
506
507 const double threshold = -0.1 *nsamp *std::log(m_MinProbability);
508 if( dLogSumMoving > threshold || dLogSumFixed > threshold
509 || dLogSumJoint > threshold )
510 {
511 // At least half the samples in B did not occur within
512 // the Parzen window width of samples in A
513 return update;
514 }
515
516 double value = 0.0;
517 value = dLogSumFixed + dLogSumMoving - dLogSumJoint;
518 value /= nsamp;
519 value += std::log(nsamp);
520
521 m_MetricTotal += value;
522 this->m_Energy += value;
523
524 derivative /= nsamp;
525 derivative /= itk::Math::sqr(m_FixedImageStandardDeviation);
526
527 double updatenorm = 0.0;
528 for( unsigned int tt = 0; tt < ImageDimension; tt++ )
529 {
530 updatenorm += derivative[tt] * derivative[tt];
531 }
532 updatenorm = std::sqrt(updatenorm);
533
534 if( updatenorm > 1.e-20 && this->GetNormalizeGradient() )
535 {
536 derivative = derivative / updatenorm;
537 }
538
539 return derivative * this->GetGradientStep();
540 }
541
542 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
543 void
544 MIRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
PrintSelf(std::ostream & os,Indent indent) const545 ::PrintSelf(std::ostream & os, Indent indent) const
546 {
547 Superclass::PrintSelf(os, indent);
548
549 os << indent << "TimeStep: " << m_TimeStep << std::endl;
550 os << indent << "FixedImageSpacing: "
551 << static_cast< typename itk::NumericTraits< SpacingType >::PrintType >( m_FixedImageSpacing )
552 << std::endl;
553 os << indent << "FixedImageOrigin: "
554 << static_cast< typename itk::NumericTraits< PointType >::PrintType >( m_FixedImageOrigin )
555 << std::endl;
556
557 itkPrintSelfObjectMacro( FixedImageGradientCalculator );
558 itkPrintSelfObjectMacro( MovingImageGradientCalculator );
559 itkPrintSelfObjectMacro( MovingImageInterpolator );
560
561 os << indent << "DenominatorThreshold: " << m_DenominatorThreshold
562 << std::endl;
563 os << indent << "IntensityDifferenceThreshold: "
564 << m_IntensityDifferenceThreshold << std::endl;
565 os << indent << "MetricTotal: " << m_MetricTotal << std::endl;
566 os << indent << "NumberOfSamples: " << m_NumberOfSamples << std::endl;
567 os << indent << "NumberOfBins: " << m_NumberOfBins << std::endl;
568 os << indent << "Minnorm: " << m_Minnorm << std::endl;
569
570 if( m_DoInverse )
571 {
572 os << indent << "DoInverse: On" << std::endl;
573 }
574 else
575 {
576 os << indent << "DoInverse: Off" << std::endl;
577 }
578 }
579 } // end namespace itk
580
581 #endif
582