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 
19 /**
20  * This tests the elastic body spline and thin plate spline
21  * transform classes by warping a unit cube into a cube with side length 3.
22  * It performs the test for 2D, 3D, and 4D to ensure that the
23  * class works in N dimensions
24  */
25 #include "itkElasticBodySplineKernelTransform.h"
26 #include "itkElasticBodyReciprocalSplineKernelTransform.h"
27 #include "itkThinPlateSplineKernelTransform.h"
28 #include "itkThinPlateR2LogRSplineKernelTransform.h"
29 #include "itkVolumeSplineKernelTransform.h"
30 #include "itkMath.h"
31 
32 
itkSplineKernelTransformTest(int,char * [])33 int itkSplineKernelTransformTest(int , char* [] )
34 {
35 
36   const double epsilon = 1e-12;
37 
38   // 2-D case
39   int i, j;
40 
41   using EBSTransform2DType = itk::ElasticBodySplineKernelTransform<double, 2>;
42   using EBRSTransform2DType = itk::ElasticBodyReciprocalSplineKernelTransform<double, 2>;
43   using TPSTransform2DType = itk::ThinPlateSplineKernelTransform<double, 2>;
44   using TPR2LRSTransform2DType = itk::ThinPlateR2LogRSplineKernelTransform<double, 2>;
45   using VSTransform2DType = itk::VolumeSplineKernelTransform<double, 2>;
46 
47   using PointType2D = EBSTransform2DType::InputPointType;
48   using Points2DIteratorType = EBSTransform2DType::PointsIterator;
49   using PointSetType2D = EBSTransform2DType::PointSetType;
50 
51   PointType2D sourcePoint2D;
52   PointType2D targetPoint2D;
53   PointType2D mappedPoint2D;
54 
55   EBSTransform2DType::Pointer     ebs2D       = EBSTransform2DType::New();
56   EBRSTransform2DType::Pointer    ebrs2D      = EBRSTransform2DType::New();
57   TPSTransform2DType::Pointer     tps2D       = TPSTransform2DType::New();
58   TPR2LRSTransform2DType::Pointer tpr2lrs2D   = TPR2LRSTransform2DType::New();
59   VSTransform2DType::Pointer      vs2D        = VSTransform2DType::New();
60 
61   // Reserve memory for the number of points
62   PointSetType2D::Pointer sourceLandmarks2D = PointSetType2D::New();
63   PointSetType2D::Pointer targetLandmarks2D = PointSetType2D::New();
64 
65   sourceLandmarks2D->GetPoints()->Reserve( 4 );
66   targetLandmarks2D->GetPoints()->Reserve( 4 );
67 
68   // Create landmark sets
69   Points2DIteratorType source2Dit = sourceLandmarks2D->GetPoints()->Begin();
70   Points2DIteratorType target2Dit = targetLandmarks2D->GetPoints()->Begin();
71 
72   Points2DIteratorType source2Dend = sourceLandmarks2D->GetPoints()->End();
73 
74   for (i = 0; i < 2; i++)
75     {
76     for (j = 0; j < 2; j++)
77       {
78       sourcePoint2D[0] = j;
79       sourcePoint2D[1] = i;
80       source2Dit.Value() = sourcePoint2D;
81       targetPoint2D[0] = 3*j;
82       targetPoint2D[1] = 3*i;
83       target2Dit.Value() = targetPoint2D;
84       source2Dit++;
85       target2Dit++;
86       }
87     }
88 
89 
90   std::cout << "EBS 2D Test:" << std::endl;
91   // Poisson's ration = 0.25, Alpha = 12.0 * ( 1 - \nu ) - 1
92   ebs2D->SetSourceLandmarks( sourceLandmarks2D );
93   ebs2D->SetTargetLandmarks( targetLandmarks2D );
94   ebs2D->SetAlpha( 12.0 * ( 1 -  0.25) - 1.0 );
95   ebs2D->ComputeWMatrix();
96 
97   { // Testing the number of parameters
98   EBSTransform2DType::ParametersType parameters1 = ebs2D->GetParameters();
99   const unsigned int numberOfParameters = parameters1.Size();
100   if( numberOfParameters != 4 * 2 )
101     {
102     std::cerr << "Number of parameters was not updated after" << std::endl;
103     std::cerr << "invoking SetSourceLandmarks and SetTargetLandmarks" << std::endl;
104     std::cerr << "Number of parameters is = " << numberOfParameters << std::endl;
105     std::cerr << "While we were expecting = " << 4 * 2 << std::endl;
106     return EXIT_FAILURE;
107     }
108   }
109 
110   source2Dit = sourceLandmarks2D->GetPoints()->Begin();
111   target2Dit = targetLandmarks2D->GetPoints()->Begin();
112 
113   source2Dend = sourceLandmarks2D->GetPoints()->End();
114   while( source2Dit != source2Dend )
115     {
116     sourcePoint2D = source2Dit.Value();
117     targetPoint2D = target2Dit.Value();
118     mappedPoint2D = ebs2D->TransformPoint(sourcePoint2D);
119     std::cout << sourcePoint2D << " : " << targetPoint2D;
120     std::cout << " warps to: " << mappedPoint2D << std::endl;
121     if( mappedPoint2D.EuclideanDistanceTo( targetPoint2D ) > epsilon )
122       {
123       return EXIT_FAILURE;
124       }
125     source2Dit++;
126     target2Dit++;
127     }
128   std::cout << std::endl;
129 
130 
131   std::cout << "EBRS 2D Test:" << std::endl;
132   ebrs2D->SetSourceLandmarks( sourceLandmarks2D );
133   ebrs2D->SetTargetLandmarks( targetLandmarks2D );
134   ebrs2D->SetAlpha( 12.0 * ( 1 -  0.25) - 1.0 );
135   ebrs2D->ComputeWMatrix();
136 
137   source2Dit = sourceLandmarks2D->GetPoints()->Begin();
138   target2Dit = targetLandmarks2D->GetPoints()->Begin();
139 
140   source2Dend = sourceLandmarks2D->GetPoints()->End();
141   while( source2Dit != source2Dend )
142     {
143     sourcePoint2D = source2Dit.Value();
144     targetPoint2D = target2Dit.Value();
145     mappedPoint2D = ebrs2D->TransformPoint(sourcePoint2D);
146     std::cout << sourcePoint2D << " : " << targetPoint2D;
147     std::cout << " warps to: " << mappedPoint2D << std::endl;
148     if( mappedPoint2D.EuclideanDistanceTo( targetPoint2D ) > epsilon )
149       {
150       return EXIT_FAILURE;
151       }
152     source2Dit++;
153     target2Dit++;
154     }
155   std::cout << std::endl;
156 
157   std::cout << "TPS 2D Test:" << std::endl;
158   tps2D->SetSourceLandmarks( sourceLandmarks2D );
159   tps2D->SetTargetLandmarks( targetLandmarks2D );
160 
161   tps2D->ComputeWMatrix();
162 
163   source2Dit = sourceLandmarks2D->GetPoints()->Begin();
164   target2Dit = targetLandmarks2D->GetPoints()->Begin();
165 
166   source2Dend = sourceLandmarks2D->GetPoints()->End();
167   while( source2Dit != source2Dend )
168     {
169     sourcePoint2D = source2Dit.Value();
170     targetPoint2D = target2Dit.Value();
171     mappedPoint2D = tps2D->TransformPoint(sourcePoint2D);
172     std::cout << sourcePoint2D << " : " << targetPoint2D;
173     std::cout << " warps to: " << mappedPoint2D << std::endl;
174     if( mappedPoint2D.EuclideanDistanceTo( targetPoint2D ) > epsilon )
175       {
176       return EXIT_FAILURE;
177       }
178     source2Dit++;
179     target2Dit++;
180     }
181   if ( tps2D->IsLinear() == true ) //NOTE TPS is never linear!
182     {
183     std::cout << "ERROR:  2D TPS reports as being a linear transform." << std::endl;
184     return EXIT_FAILURE;
185     }
186 
187 
188   //NOTE: The following should set the default values explicitly
189     {
190     constexpr double TestValue = 0.012345;
191     tps2D->SetStiffness(TestValue); //This value should not change the result at all.
192 
193     if ( itk::Math::NotExactlyEquals(tps2D->GetStiffness(), TestValue) )
194       {
195       std::cout << "ERROR:  Explicitly set stiffness value not retained." << std::endl;
196       return EXIT_FAILURE;
197       }
198     }
199     { //Just for code coverage
200     TPSTransform2DType::VectorSetType::ConstPointer tempDisplacements=tps2D->GetDisplacements();
201 
202       {
203       TPSTransform2DType::InputVectorType  testVector;
204       testVector[0] = 0.0;
205       testVector[1] = 1.0;
206       bool exceptionCaught=false;
207       try
208         {
209         tps2D->TransformVector(testVector);
210         }
211       catch(...)
212         {
213         exceptionCaught = true;
214         }
215       if ( exceptionCaught != true )
216         {
217         return EXIT_FAILURE;
218         }
219       }
220       {
221       TPSTransform2DType::InputVnlVectorType  testVector;
222       testVector[0] = 0.0;
223       testVector[1] = 1.0;
224       bool exceptionCaught=false;
225       try
226         {
227         tps2D->TransformVector(testVector);
228         }
229       catch(...)
230         {
231         exceptionCaught = true;
232         }
233       if ( exceptionCaught != true )
234         {
235         return EXIT_FAILURE;
236         }
237       }
238       {
239       TPSTransform2DType::InputCovariantVectorType testVector;
240       testVector[0] = 0.0;
241       testVector[1] = 1.0;
242       bool exceptionCaught=false;
243       try
244         {
245         tps2D->TransformCovariantVector(testVector);
246         }
247       catch(...)
248         {
249         exceptionCaught = true;
250         }
251       if ( exceptionCaught != true )
252         {
253         return EXIT_FAILURE;
254         }
255       }
256       {
257       TPSTransform2DType::JacobianPositionType   testJacobian;
258       TPSTransform2DType::InputPointType testVector;
259       testVector[0] = 0.0;
260       testVector[1] = 1.0;
261       bool exceptionCaught=false;
262       try
263         {
264         tps2D->ComputeJacobianWithRespectToPosition(testVector,testJacobian);
265         }
266       catch(...)
267         {
268         exceptionCaught = true;
269         }
270       if ( exceptionCaught != true )
271         {
272         return EXIT_FAILURE;
273         }
274       }
275     }
276 
277   std::cout << std::endl;
278 
279   std::cout << "TPR2LR 2D Test:" << std::endl;
280   tpr2lrs2D->SetSourceLandmarks( sourceLandmarks2D );
281   tpr2lrs2D->SetTargetLandmarks( targetLandmarks2D );
282 
283   tpr2lrs2D->ComputeWMatrix();
284 
285   source2Dit = sourceLandmarks2D->GetPoints()->Begin();
286   target2Dit = targetLandmarks2D->GetPoints()->Begin();
287 
288   source2Dend = sourceLandmarks2D->GetPoints()->End();
289   while( source2Dit != source2Dend )
290     {
291     sourcePoint2D = source2Dit.Value();
292     targetPoint2D = target2Dit.Value();
293     mappedPoint2D = tpr2lrs2D->TransformPoint(sourcePoint2D);
294     std::cout << sourcePoint2D << " : " << targetPoint2D;
295     std::cout << " warps to: " << mappedPoint2D << std::endl;
296     if( mappedPoint2D.EuclideanDistanceTo( targetPoint2D ) > epsilon )
297       {
298       return EXIT_FAILURE;
299       }
300     source2Dit++;
301     target2Dit++;
302     }
303   std::cout << std::endl;
304 
305   // volume spline transform
306   std::cout << "VS 2D Test:" << std::endl;
307   vs2D->SetSourceLandmarks( sourceLandmarks2D );
308   vs2D->SetTargetLandmarks( targetLandmarks2D );
309 
310   vs2D->ComputeWMatrix();
311 
312   source2Dit = sourceLandmarks2D->GetPoints()->Begin();
313   target2Dit = targetLandmarks2D->GetPoints()->Begin();
314 
315   source2Dend = sourceLandmarks2D->GetPoints()->End();
316   while( source2Dit != source2Dend )
317     {
318     sourcePoint2D = source2Dit.Value();
319     targetPoint2D = target2Dit.Value();
320     mappedPoint2D = vs2D->TransformPoint(sourcePoint2D);
321     std::cout << sourcePoint2D << " : " << targetPoint2D;
322     std::cout << " warps to: " << mappedPoint2D << std::endl;
323     if( mappedPoint2D.EuclideanDistanceTo( targetPoint2D ) > epsilon )
324       {
325       return EXIT_FAILURE;
326       }
327     source2Dit++;
328     target2Dit++;
329     }
330   std::cout << std::endl;
331 
332 
333    // 3-D case
334   int k;
335   using EBSTransform3DType = itk::ElasticBodySplineKernelTransform<double, 3>;
336   using TPSTransform3DType = itk::ThinPlateSplineKernelTransform<double, 3>;
337 
338   using PointType3D = EBSTransform3DType::InputPointType;
339   using Points3DIteratorType = EBSTransform3DType::PointsIterator;
340 
341   PointType3D sourcePoint3D;
342   PointType3D targetPoint3D;
343   PointType3D mappedPoint3D;
344 
345   // Reserve memory for the number of points
346   EBSTransform3DType::Pointer ebs3D = EBSTransform3DType::New();
347   ebs3D->GetModifiableTargetLandmarks()->GetPoints()->Reserve( 8 );
348   ebs3D->GetModifiableSourceLandmarks()->GetPoints()->Reserve( 8 );
349 
350   TPSTransform3DType::Pointer tps3D = TPSTransform3DType::New();
351   tps3D->GetModifiableTargetLandmarks()->GetPoints()->Reserve( 8 );
352   tps3D->GetModifiableSourceLandmarks()->GetPoints()->Reserve( 8 );
353 
354 
355   // Create landmark sets
356   Points3DIteratorType ebs3Ds = ebs3D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
357   Points3DIteratorType ebs3Dt = ebs3D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
358   Points3DIteratorType tps3Ds = tps3D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
359   Points3DIteratorType tps3Dt = tps3D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
360 
361   Points3DIteratorType ebs3DsEnd  = ebs3D->GetModifiableSourceLandmarks()->GetPoints()->End();
362   Points3DIteratorType tps3DsEnd  = tps3D->GetModifiableSourceLandmarks()->GetPoints()->End();
363 
364   for (i = 0; i < 2; i++)
365     {
366     for (j = 0; j < 2; j++)
367       {
368       for (k = 0; k < 2; k++)
369         {
370         sourcePoint3D[0] = k;
371         sourcePoint3D[1] = j;
372         sourcePoint3D[2] = i;
373         ebs3Ds.Value() = sourcePoint3D;
374         tps3Ds.Value() = sourcePoint3D;
375         targetPoint3D[0] = 3*k;
376         targetPoint3D[1] = 3*j;
377         targetPoint3D[2] = 3*i;
378         ebs3Dt.Value() = targetPoint3D;
379         tps3Dt.Value() = targetPoint3D;
380         ebs3Ds++;
381         ebs3Dt++;
382         tps3Ds++;
383         tps3Dt++;
384         }
385       }
386     }
387 
388   std::cout << "EBS 3D Test:" << std::endl;
389   // Poisson's ration = 0.25, Alpha = 12.0 * ( 1 - \nu ) - 1
390   ebs3D->SetAlpha( 12.0 * ( 1 -  0.25) - 1.0 );
391   ebs3D->ComputeWMatrix();
392 
393   ebs3Ds     = ebs3D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
394   ebs3Dt     = ebs3D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
395   ebs3DsEnd  = ebs3D->GetModifiableSourceLandmarks()->GetPoints()->End();
396 
397   while( ebs3Ds != ebs3DsEnd )
398   {
399     sourcePoint3D = ebs3Ds.Value();
400     targetPoint3D = ebs3Dt.Value();
401     mappedPoint3D = ebs3D->TransformPoint(sourcePoint3D);
402     std::cout << sourcePoint3D << " : " << targetPoint3D;
403     std::cout << " warps to: " << mappedPoint3D << std::endl;
404     if( mappedPoint3D.EuclideanDistanceTo( targetPoint3D ) > epsilon )
405     {
406       return EXIT_FAILURE;
407     }
408     ebs3Ds++;
409     ebs3Dt++;
410   }
411   std::cout << std::endl;
412 
413   std::cout << "TPS 3D Test:" << std::endl;
414 
415   tps3D->ComputeWMatrix();
416 
417   tps3Ds = tps3D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
418   tps3Dt = tps3D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
419   tps3DsEnd  = tps3D->GetModifiableSourceLandmarks()->GetPoints()->End();
420 
421   while( tps3Ds != tps3DsEnd )
422   {
423     sourcePoint3D = tps3Ds.Value();
424     targetPoint3D = tps3Dt.Value();
425     mappedPoint3D = tps3D->TransformPoint(sourcePoint3D);
426     std::cout << sourcePoint3D << " : " << targetPoint3D;
427     std::cout << " warps to: " << mappedPoint3D << std::endl;
428     if( mappedPoint3D.EuclideanDistanceTo( targetPoint3D ) > epsilon )
429     {
430       return EXIT_FAILURE;
431     }
432     tps3Ds++;
433     tps3Dt++;
434   }
435   std::cout << std::endl;
436 
437   std::cout << "Get/Set Parameters test " << std::endl;
438   TPSTransform3DType::ParametersType parameters1 = tps3D->GetParameters();
439   tps3D->SetParameters( parameters1 );
440   TPSTransform3DType::ParametersType parameters2 = tps3D->GetParameters();
441   const unsigned int numberOfParameters = parameters1.Size();
442   const double tolerance = 1e-7;
443   for(unsigned int pr = 0; pr < numberOfParameters; pr++)
444     {
445     if( itk::Math::abs( parameters1[pr] - parameters2[pr] ) > tolerance )
446       {
447       std::cout << "Parameters were not correctly recovered " << std::endl;
448       return EXIT_FAILURE;
449       }
450     }
451   std::cout << "Get/Set Parameters Passed" << std::endl << std::endl;
452 
453 
454   // 4-D case
455   int l;
456   using EBSTransform4DType = itk::ElasticBodySplineKernelTransform<double, 4>;
457   using TPSTransform4DType = itk::ThinPlateSplineKernelTransform<double, 4>;
458 
459   using PointType4D = EBSTransform4DType::InputPointType;
460   using Points4DIteratorType = EBSTransform4DType::PointsIterator;
461 
462   PointType4D sourcePoint4D;
463   PointType4D targetPoint4D;
464   PointType4D mappedPoint4D;
465 
466   EBSTransform4DType::Pointer ebs4D = EBSTransform4DType::New();
467   TPSTransform4DType::Pointer tps4D = TPSTransform4DType::New();
468 
469   // Reserve memory for the number of points
470   ebs4D->GetModifiableTargetLandmarks()->GetPoints()->Reserve( 16 );
471   tps4D->GetModifiableTargetLandmarks()->GetPoints()->Reserve( 16 );
472 
473   ebs4D->GetModifiableSourceLandmarks()->GetPoints()->Reserve( 16 );
474   tps4D->GetModifiableSourceLandmarks()->GetPoints()->Reserve( 16 );
475 
476   // Create landmark sets
477   Points4DIteratorType ebs4Ds = ebs4D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
478   Points4DIteratorType ebs4Dt = ebs4D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
479   Points4DIteratorType tps4Ds = tps4D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
480   Points4DIteratorType tps4Dt = tps4D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
481 
482   Points4DIteratorType ebs4DsEnd  = ebs4D->GetModifiableSourceLandmarks()->GetPoints()->End();
483   Points4DIteratorType tps4DsEnd  = tps4D->GetModifiableSourceLandmarks()->GetPoints()->End();
484 
485   for (i = 0; i < 2; i++)
486     {
487     for (j = 0; j < 2; j++)
488       {
489       for (k = 0; k < 2; k++)
490         {
491         for (l = 0; l < 2; l++)
492           {
493           sourcePoint4D[0] = l;
494           sourcePoint4D[1] = k;
495           sourcePoint4D[2] = j;
496           sourcePoint4D[3] = i;
497           ebs4Ds.Value() = sourcePoint4D;
498           tps4Ds.Value() = sourcePoint4D;
499           targetPoint4D[0] = 3*l;
500           targetPoint4D[1] = 3*k;
501           targetPoint4D[2] = 3*j;
502           targetPoint4D[3] = 3*i;
503           ebs4Dt.Value() = targetPoint4D;
504           tps4Dt.Value() = targetPoint4D;
505           ebs4Ds++;
506           ebs4Dt++;
507           tps4Ds++;
508           tps4Dt++;
509           }
510         }
511       }
512     }
513   std::cout << "EBS 4D Test:" << std::endl;
514   // Poisson's ration = 0.25, Alpha = 12.0 * ( 1 - \nu ) - 1
515   ebs4D->SetAlpha( 12.0 * ( 1 -  0.25) - 1.0 );
516   ebs4D->ComputeWMatrix();
517 
518   ebs4Ds = ebs4D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
519   ebs4Dt = ebs4D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
520   ebs4DsEnd  = ebs4D->GetModifiableSourceLandmarks()->GetPoints()->End();
521 
522   while( ebs4Ds != ebs4DsEnd )
523   {
524     sourcePoint4D = ebs4Ds.Value();
525     targetPoint4D = ebs4Dt.Value();
526     mappedPoint4D = ebs4D->TransformPoint(sourcePoint4D);
527     std::cout << sourcePoint4D << " : " << targetPoint4D;
528     std::cout << " warps to: " << mappedPoint4D << std::endl;
529     if( mappedPoint4D.EuclideanDistanceTo( targetPoint4D ) > epsilon )
530     {
531       return EXIT_FAILURE;
532     }
533     ebs4Ds++;
534     ebs4Dt++;
535   }
536   std::cout << std::endl;
537 
538   std::cout << "TPS 4D Test:" << std::endl;
539   tps4D->ComputeWMatrix();
540 
541   tps4Ds = tps4D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
542   tps4Dt = tps4D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
543   tps4DsEnd  = tps4D->GetModifiableSourceLandmarks()->GetPoints()->End();
544   while( tps4Ds != tps4DsEnd )
545   {
546     sourcePoint4D = tps4Ds.Value();
547     targetPoint4D = tps4Dt.Value();
548     mappedPoint4D = tps4D->TransformPoint(sourcePoint4D);
549     std::cout << sourcePoint4D << " : " << targetPoint4D;
550     std::cout << " warps to: " << mappedPoint4D << std::endl;
551     if( mappedPoint4D.EuclideanDistanceTo( targetPoint4D ) > epsilon )
552     {
553       return EXIT_FAILURE;
554     }
555     tps4Ds++;
556     tps4Dt++;
557   }
558   std::cout << std::endl;
559 
560   std::cout << ebs2D << std::endl;
561 
562   std::cout << "TEST DONE" << std::endl;
563 
564   return EXIT_SUCCESS;
565 
566 }
567