]> Creatis software - clitk.git/blob - tools/clitkPointRigidRegistration.cxx
Update Landmark coordinates with transformation matrix
[clitk.git] / tools / clitkPointRigidRegistration.cxx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
5   - University of LYON              http://www.universite-lyon.fr/
6   - Léon Bérard cancer center       http://www.centreleonberard.fr
7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
8
9   This software is distributed WITHOUT ANY WARRANTY; without even
10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11   PURPOSE.  See the copyright notices for more information.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18
19 /* =================================================
20  * @file   clitkPointRigidRegistrationGenericFilter.txx
21  * @author xxx <xxx@creatis.insa-lyon.fr>
22  * @date   29 June 2029
23  *
24  * @brief PointRigidRegistration an image
25  *
26  ===================================================*/
27
28 // clitk
29 #include "clitkPointRigidRegistration_ggo.h"
30
31
32 //paste from RigidRegistration
33 #include "itkImageFileReader.h"
34 #include "itkImageFileWriter.h"
35 #include "itkImage.h"
36 #include "itkVector.h"
37 #include "itkResampleImageFilter.h"
38 #include "itkLandmarkBasedTransformInitializer.h"
39 #include "itkRigid2DTransform.h"
40 #include "itkVersorRigid3DTransform.h"
41 #include <iostream>
42
43
44 #include "itkMatrix.h"
45 #include "itkArray.h"
46 //#include "itkPoint.h"
47 #include "clitkImageCommon.h"
48 #include "clitkCommon.h"
49 //#define VTK_EXCLUDE_STRSTREAM_HEADERS
50 #include <vtkMatrix4x4.h>
51 #include <vtkSmartPointer.h>
52
53 //for open file for reading
54 #include "clitkIO.h"
55 #include "clitkImageCommon.h"
56 #include "clitkCommon.h"
57
58
59 //--------------------------------------------------------------------
60 int main(int argc, char * argv[])
61 {
62
63   // Init command line
64   GGO(clitkPointRigidRegistration, args_info);
65   CLITK_INIT;
66
67
68
69
70   // Iinit itk
71   // read input file
72   //open stream to reading
73   std::ifstream is;
74   clitk::openFileForReading(is, args_info.input_arg);
75
76
77
78   //reading first line of input file to chck thw dimension of data
79   double x = 0;
80   //clitk::skipComment(is);
81   is >> x;
82
83     typedef   unsigned char  PixelType;
84
85
86
87
88
89     unsigned int Dimension_temp = (unsigned int)x;
90
91
92
93         if (Dimension_temp==2)
94         {
95             const     unsigned int   Dimension = 2;
96             typedef   itk::Image< PixelType, Dimension > ImageType;
97             typedef   float          VectorComponentType;
98             typedef   itk::Vector< VectorComponentType, Dimension >    VectorType;
99             //Typ LandmarkBasedTransormInitializer
100             typedef itk::Rigid2DTransform< double > Rigid2DTransformType;
101             typedef itk::LandmarkBasedTransformInitializer< Rigid2DTransformType,ImageType, ImageType>
102                 LandmarkBasedTransformInitializerType;
103
104             LandmarkBasedTransformInitializerType::Pointer landmarkBasedTransformInitializer = LandmarkBasedTransformInitializerType::New();
105
106             //  Create source and target landmarks.
107             typedef LandmarkBasedTransformInitializerType::LandmarkPointContainer     LandmarkContainerType;
108             typedef LandmarkBasedTransformInitializerType::LandmarkPointType          LandmarkPointType;
109
110             LandmarkContainerType imageLandmarks;
111             LandmarkContainerType trackerLandmarks;
112
113             LandmarkPointType imagePoint;
114             LandmarkPointType trackerPoint;
115
116             is >> x;
117
118
119             while (is && !is.eof()) {
120                 imagePoint[0] = x;
121                 is >> imagePoint[1];
122
123
124                 is >> trackerPoint[0];
125                 is >> trackerPoint[1];
126
127                 imageLandmarks.push_back(imagePoint );
128                 trackerLandmarks.push_back(trackerPoint );
129
130
131               is >> x;
132             }
133
134             is.close();
135
136             landmarkBasedTransformInitializer->SetFixedLandmarks( imageLandmarks);
137             landmarkBasedTransformInitializer->SetMovingLandmarks( trackerLandmarks);
138
139             Rigid2DTransformType::Pointer transform = Rigid2DTransformType::New();
140
141             transform->SetIdentity();
142
143             landmarkBasedTransformInitializer->SetTransform(transform);
144             landmarkBasedTransformInitializer->InitializeTransform();
145
146             Rigid2DTransformType::MatrixType matrix = transform->GetMatrix();
147             Rigid2DTransformType::OffsetType offset = transform->GetOffset();
148
149             // Write result
150               std::ofstream out;
151               clitk::openFileForWriting(out, args_info.output_arg);
152
153               out << matrix[0][0] << ' ' << matrix[0][1] << ' ' << offset[0] << std::endl;
154               out << matrix[1][0] << ' ' << matrix[1][1] << ' ' << offset[1] << std::endl;
155               out << 0.0  << ' ' << 0.0 << ' ' << 1.0;
156               out.close();
157
158
159
160               // Write result
161              if (args_info.error_arg != 0)
162              {
163                 std::ofstream error;
164                 clitk::openFileForWriting(error, args_info.error_arg);
165
166
167               //imagePoints
168               int nImagePoints = imageLandmarks.size();
169               double  imagePoints[3][nImagePoints];
170               double  trackerPoints[3][nImagePoints];
171               double  imagePointstoTrackerCoordinate[3][nImagePoints];
172               double FRE = 0;
173               double PointEuclideanDistance = 0;
174               double SumofPointsEuclideanDistance = 0;
175
176
177               for(int n= 0; n < imageLandmarks.size(); n++)
178               {
179                   for(int j= 0; j < 2; j++)
180                   {
181                       imagePoints[j][n] = imageLandmarks[n][j];
182                       trackerPoints[j][n] = trackerLandmarks[n][j];
183
184
185
186
187 //                      imagePoints[j][n] = trackerLandmarks[n][j];
188 //                      trackerPoints[j][n] = imageLandmarks[n][j];
189                       imagePointstoTrackerCoordinate[j][n] = 0;
190
191                   }
192                   imagePoints[2][n] = 1.0;
193                   trackerPoints[2][n] = 1.0;
194                   imagePointstoTrackerCoordinate[2][n] = 0;
195               }
196
197
198               //transformation matrix
199               double RigidTransformationImagetoTracker[4][4];
200
201
202               for(int i = 0; i < 2; i++)
203                 for(int j = 0; j < 2; j++)
204                     RigidTransformationImagetoTracker[i][j] = matrix[i][j];
205
206
207               for(int i = 0; i < 2; i++)
208                   RigidTransformationImagetoTracker[i][2] = offset[i];
209
210               RigidTransformationImagetoTracker[2][0] = 0.0;
211               RigidTransformationImagetoTracker[2][1] = 0.0;
212               RigidTransformationImagetoTracker[2][2] = 1.0;
213
214
215               //writing results
216               error << "ImagePoints:    CorrespondingTrackerPoints:     ImagePointstoTrackerCoordinateSystem:  PointsDistanceafterRigidRegistration:" << std::endl;
217
218
219
220
221
222               //Calculate FRE
223
224               for(int n= 0; n < imageLandmarks.size(); n++)
225               {
226
227
228                    //Calculate imagePointstoTrackerCoordinate
229                    for(int i = 0; i < 3; i++)
230                    {
231
232                        for(int k = 0; k < 3; k++)
233                        {
234                             imagePointstoTrackerCoordinate[i][n] = imagePointstoTrackerCoordinate[i][n] + RigidTransformationImagetoTracker[i][k]*imagePoints[k][n];
235
236                        }
237
238                    }
239
240
241                    //writing results
242
243                    for(int i = 0; i < 2; i++)
244                             error << imagePoints[i][n] << ' ';
245
246                    error << "    ";
247
248                    for(int i = 0; i < 2; i++)
249                             error << trackerPoints[i][n] << ' ';
250
251                    error << "   ";
252
253                    for(int i = 0; i < 2; i++)
254                             error << imagePointstoTrackerCoordinate[i][n] << ' ';
255
256
257
258                    //calculate PointEuclideanDistance
259                    PointEuclideanDistance = 0;
260                    for(int i=0; i < 2;i++)
261                    {
262                        PointEuclideanDistance = PointEuclideanDistance + std::pow((trackerPoints[i][n] - imagePointstoTrackerCoordinate[i][n]),2);
263                    }
264                    PointEuclideanDistance = std::sqrt(PointEuclideanDistance);
265                    error << "   " <<   PointEuclideanDistance   << "   ";
266
267                    SumofPointsEuclideanDistance = SumofPointsEuclideanDistance +  std::pow(PointEuclideanDistance,2);
268
269
270
271
272
273
274                    error << std::endl;
275
276
277
278
279
280               } //end loop via points
281
282               error << std::endl;
283               SumofPointsEuclideanDistance = std::sqrt(SumofPointsEuclideanDistance/nImagePoints);
284               error << "Fiducial Registration Error:   " <<  SumofPointsEuclideanDistance;
285               error.close();
286
287             } //end if (args_info.error_arg != null)
288         } //end Dimension == 2 case
289         else if (Dimension_temp==3)
290         {
291
292             const     unsigned int   Dimension = 3;
293             typedef   itk::Image< PixelType, Dimension > ImageType;
294             typedef   float          VectorComponentType;
295             typedef   itk::Vector< VectorComponentType, Dimension >    VectorType;
296             //Typ LandmarkBasedTransormInitializer
297             typedef itk::VersorRigid3DTransform< double > Rigid3DTransformType;
298             typedef itk::LandmarkBasedTransformInitializer< Rigid3DTransformType,ImageType, ImageType>
299                 LandmarkBasedTransformInitializerType;
300
301             LandmarkBasedTransformInitializerType::Pointer landmarkBasedTransformInitializer = LandmarkBasedTransformInitializerType::New();
302
303             //  Create source and target landmarks.
304             typedef LandmarkBasedTransformInitializerType::LandmarkPointContainer     LandmarkContainerType;
305             typedef LandmarkBasedTransformInitializerType::LandmarkPointType          LandmarkPointType;
306
307             LandmarkContainerType imageLandmarks;
308             LandmarkContainerType trackerLandmarks;
309             LandmarkContainerType imageLandmarksTotrackerCoordinate;
310
311
312             LandmarkPointType imagePoint;
313             LandmarkPointType trackerPoint;
314
315             is >> x;
316
317             while (is && !is.eof()) {
318                 imagePoint[0] = x;
319                 is >> imagePoint[1];
320                 is >> imagePoint[2];
321
322                 is >> trackerPoint[0];
323                 is >> trackerPoint[1];
324                 is >> trackerPoint[2];
325
326                 imageLandmarks.push_back(imagePoint );
327                 trackerLandmarks.push_back(trackerPoint );
328
329                 is >> x;
330             }
331
332             is.close();
333
334             landmarkBasedTransformInitializer->SetFixedLandmarks( imageLandmarks);
335             landmarkBasedTransformInitializer->SetMovingLandmarks( trackerLandmarks);
336
337             Rigid3DTransformType::Pointer transform = Rigid3DTransformType::New();
338
339             transform->SetIdentity();
340
341             landmarkBasedTransformInitializer->SetTransform(transform);
342             landmarkBasedTransformInitializer->InitializeTransform();
343
344             Rigid3DTransformType::MatrixType matrix = transform->GetMatrix();
345             Rigid3DTransformType::OffsetType offset = transform->GetOffset();
346
347             // Write result
348               std::ofstream out;
349               clitk::openFileForWriting(out, args_info.output_arg);
350
351               out << matrix[0][0] << ' ' << matrix[0][1] << ' ' << matrix[0][2] << ' '<< offset[0] << std::endl;
352               out << matrix[1][0] << ' ' << matrix[1][1] << ' ' << matrix[1][2] << ' '<< offset[1] << std::endl;
353               out << matrix[2][0] << ' ' << matrix[2][1] << ' ' << matrix[2][2] << ' '<< offset[2] << std::endl;
354               out << 0.0  << ' ' << 0.0 << ' ' << 0.0 << ' ' << 1.0;
355               out.close();
356
357
358               if (args_info.error_arg != 0)
359               {
360               // Write result
361                 std::ofstream error;
362                 clitk::openFileForWriting(error, args_info.error_arg);
363
364
365               //imagePoints
366               int nImagePoints = imageLandmarks.size();
367               double  imagePoints[4][nImagePoints];
368               double  trackerPoints[4][nImagePoints];
369               double  imagePointstoTrackerCoordinate[4][nImagePoints];
370               double FRE = 0;
371               double PointEuclideanDistance = 0;
372               double SumofPointsEuclideanDistance = 0;
373
374
375               for(int n= 0; n < imageLandmarks.size(); n++)
376               {
377                   for(int j= 0; j < 3; j++)
378                   {
379                       imagePoints[j][n] = imageLandmarks[n][j];
380                       trackerPoints[j][n] = trackerLandmarks[n][j];
381
382
383
384
385 //                      imagePoints[j][n] = trackerLandmarks[n][j];
386 //                      trackerPoints[j][n] = imageLandmarks[n][j];
387                       imagePointstoTrackerCoordinate[j][n] = 0;
388
389                   }
390                   imagePoints[3][n] = 1.0;
391                   trackerPoints[3][n] = 1.0;
392                   imagePointstoTrackerCoordinate[3][n] = 0;
393               }
394
395
396               //transformation matrix
397               double RigidTransformationImagetoTracker[4][4];
398
399
400               for(int i = 0; i < 3; i++)
401                 for(int j = 0; j < 3; j++)
402                     RigidTransformationImagetoTracker[i][j] = matrix[i][j];
403
404
405               for(int i = 0; i < 3; i++)
406                   RigidTransformationImagetoTracker[i][3] = offset[i];
407
408               RigidTransformationImagetoTracker[3][0] = 0.0;
409               RigidTransformationImagetoTracker[3][1] = 0.0;
410               RigidTransformationImagetoTracker[3][2] = 0.0;
411               RigidTransformationImagetoTracker[3][3] = 1.0;
412
413               //writing results
414               error << "ImagePoints:    CorrespondingTrackerPoints:     ImagePointstoTrackerCoordinateSystem:  PointsDistanceafterRigidRegistration:" << std::endl;
415
416
417
418
419
420               //Calculate FRE
421
422               for(int n= 0; n < imageLandmarks.size(); n++)
423               {
424
425
426                    //Calculate imagePointstoTrackerCoordinate
427                    for(int i = 0; i < 4; i++)
428                    {
429
430                        for(int k = 0; k < 4; k++)
431                        {
432                             imagePointstoTrackerCoordinate[i][n] = imagePointstoTrackerCoordinate[i][n] + RigidTransformationImagetoTracker[i][k]*imagePoints[k][n];
433
434                        }
435
436                    }
437
438
439                    //writing results
440
441                    for(int i = 0; i < 3; i++)
442                             error << imagePoints[i][n] << ' ';
443
444                    error << "    ";
445
446                    for(int i = 0; i < 3; i++)
447                             error << trackerPoints[i][n] << ' ';
448
449                    error << "   ";
450
451                    for(int i = 0; i < 3; i++)
452                             error << imagePointstoTrackerCoordinate[i][n] << ' ';
453
454
455
456                    //calculate PointEuclideanDistance
457                    PointEuclideanDistance = 0;
458                    for(int i=0; i < 3;i++)
459                    {
460                        PointEuclideanDistance = PointEuclideanDistance + std::pow((trackerPoints[i][n] - imagePointstoTrackerCoordinate[i][n]),2);
461                    }
462                    PointEuclideanDistance = std::sqrt(PointEuclideanDistance);
463                    error << "   " <<   PointEuclideanDistance   << "   ";
464
465                    SumofPointsEuclideanDistance = SumofPointsEuclideanDistance +  std::pow(PointEuclideanDistance,2);
466
467
468
469
470
471
472                    error << std::endl;
473
474
475
476
477
478               } //end loop via points
479
480               error << std::endl;
481               SumofPointsEuclideanDistance = std::sqrt(SumofPointsEuclideanDistance/nImagePoints);
482               error << "Fiducial Registration Error:   " <<  SumofPointsEuclideanDistance;
483               error.close();
484
485
486               } //end if (args_info.error_arg != null)
487           } //end Dimension == 3 case
488         else
489         {
490               is.close();
491               return EXIT_FAILURE;
492
493         }
494
495
496
497
498
499
500
501   return EXIT_SUCCESS;
502 }// end main
503 //--------------------------------------------------------------------