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