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 //--------------------------------------------------------------------