]> Creatis software - clitk.git/blobdiff - tools/clitkPointRigidRegistration.cxx
Add comment to precise the functionality of the inputs
[clitk.git] / tools / clitkPointRigidRegistration.cxx
index e66edd25a70abf6a3d7b842c594207f7dd5b3c72..065b7d904847ba3081f054505ab1e68e9b092b5a 100644 (file)
@@ -27,7 +27,7 @@
 
 // clitk
 #include "clitkPointRigidRegistration_ggo.h"
-#include "clitkPointRigidRegistrationGenericFilter.h"
+
 
 //paste from RigidRegistration
 #include "itkImageFileReader.h"
 #include "itkVersorRigid3DTransform.h"
 #include <iostream>
 
-//paste from /home/dspinczyk/dev/clitk_superbuild_Agata/Source/clitk/common/clitkTransformUtilities.h
+
 #include "itkMatrix.h"
 #include "itkArray.h"
-#include "itkPoint.h"
+//#include "itkPoint.h"
 #include "clitkImageCommon.h"
 #include "clitkCommon.h"
 //#define VTK_EXCLUDE_STRSTREAM_HEADERS
@@ -75,20 +75,17 @@ int main(int argc, char * argv[])
 
 
 
-  //pojedyncza zaczytana wartosc
+  //reading first line of input file to chck thw dimension of data
   double x = 0;
   //clitk::skipComment(is);
   is >> x;
 
-  //typ piksela - bo wykorzystujemy strukture itk::Image
-  typedef   unsigned char  PixelType;
+    typedef   unsigned char  PixelType;
 
 
 
 
 
-    //Ddefinicja typu obrazu: zawiera rozmiar danych, typ piksela
-    //rozmiarowosc danych
     unsigned int Dimension_temp = (unsigned int)x;
 
 
@@ -110,25 +107,22 @@ int main(int argc, char * argv[])
             typedef LandmarkBasedTransformInitializerType::LandmarkPointContainer     LandmarkContainerType;
             typedef LandmarkBasedTransformInitializerType::LandmarkPointType          LandmarkPointType;
 
-            //bufory na zaczyt wspolrzedncyh
             LandmarkContainerType imageLandmarks;
             LandmarkContainerType trackerLandmarks;
 
-            //bufory na pojedyncze punkty
-            LandmarkPointType imagePoint; // dane z CT
-            LandmarkPointType trackerPoint; //dane z tracking system
+            LandmarkPointType imagePoint;
+            LandmarkPointType trackerPoint;
 
             is >> x;
 
 
             while (is && !is.eof()) {
-                trackerPoint[0] = x;
-                is >> trackerPoint[1];
-                //is >> trackerPoint[2];
-
-                is >> imagePoint[0];
+                imagePoint[0] = x;
                 is >> imagePoint[1];
-                //is >> imagePoint[2];
+
+
+                is >> trackerPoint[0];
+                is >> trackerPoint[1];
 
                 imageLandmarks.push_back(imagePoint );
                 trackerLandmarks.push_back(trackerPoint );
@@ -146,7 +140,6 @@ int main(int argc, char * argv[])
 
             transform->SetIdentity();
 
-            // Versor rotacji i macierz translacji
             landmarkBasedTransformInitializer->SetTransform(transform);
             landmarkBasedTransformInitializer->InitializeTransform();
 
@@ -163,11 +156,139 @@ int main(int argc, char * argv[])
               out.close();
 
 
-        }
+
+              // Write result
+             if (args_info.error_arg != 0)
+             {
+                std::ofstream error;
+                clitk::openFileForWriting(error, args_info.error_arg);
+
+
+              //imagePoints
+              int nImagePoints = imageLandmarks.size();
+              double  imagePoints[3][nImagePoints];
+              double  trackerPoints[3][nImagePoints];
+              double  imagePointstoTrackerCoordinate[3][nImagePoints];
+              double FRE = 0;
+              double PointEuclideanDistance = 0;
+              double SumofPointsEuclideanDistance = 0;
+
+
+              for(int n= 0; n < imageLandmarks.size(); n++)
+              {
+                  for(int j= 0; j < 2; j++)
+                  {
+                      imagePoints[j][n] = imageLandmarks[n][j];
+                      trackerPoints[j][n] = trackerLandmarks[n][j];
+
+
+
+
+//                      imagePoints[j][n] = trackerLandmarks[n][j];
+//                      trackerPoints[j][n] = imageLandmarks[n][j];
+                      imagePointstoTrackerCoordinate[j][n] = 0;
+
+                  }
+                  imagePoints[2][n] = 1.0;
+                  trackerPoints[2][n] = 1.0;
+                  imagePointstoTrackerCoordinate[2][n] = 0;
+              }
+
+
+              //transformation matrix
+              double RigidTransformationImagetoTracker[4][4];
+
+
+              for(int i = 0; i < 2; i++)
+                for(int j = 0; j < 2; j++)
+                    RigidTransformationImagetoTracker[i][j] = matrix[i][j];
+
+
+              for(int i = 0; i < 2; i++)
+                  RigidTransformationImagetoTracker[i][2] = offset[i];
+
+              RigidTransformationImagetoTracker[2][0] = 0.0;
+              RigidTransformationImagetoTracker[2][1] = 0.0;
+              RigidTransformationImagetoTracker[2][2] = 1.0;
+
+
+              //writing results
+              error << "ImagePoints:    CorrespondingTrackerPoints:     ImagePointstoTrackerCoordinateSystem:  PointsDistanceafterRigidRegistration:" << std::endl;
+
+
+
+
+
+              //Calculate FRE
+
+              for(int n= 0; n < imageLandmarks.size(); n++)
+              {
+
+
+                   //Calculate imagePointstoTrackerCoordinate
+                   for(int i = 0; i < 3; i++)
+                   {
+
+                       for(int k = 0; k < 3; k++)
+                       {
+                            imagePointstoTrackerCoordinate[i][n] = imagePointstoTrackerCoordinate[i][n] + RigidTransformationImagetoTracker[i][k]*imagePoints[k][n];
+
+                       }
+
+                   }
+
+
+                   //writing results
+
+                   for(int i = 0; i < 2; i++)
+                            error << imagePoints[i][n] << ' ';
+
+                   error << "    ";
+
+                   for(int i = 0; i < 2; i++)
+                            error << trackerPoints[i][n] << ' ';
+
+                   error << "   ";
+
+                   for(int i = 0; i < 2; i++)
+                            error << imagePointstoTrackerCoordinate[i][n] << ' ';
+
+
+
+                   //calculate PointEuclideanDistance
+                   PointEuclideanDistance = 0;
+                   for(int i=0; i < 2;i++)
+                   {
+                       PointEuclideanDistance = PointEuclideanDistance + std::pow((trackerPoints[i][n] - imagePointstoTrackerCoordinate[i][n]),2);
+                   }
+                   PointEuclideanDistance = std::sqrt(PointEuclideanDistance);
+                   error << "   " <<   PointEuclideanDistance   << "   ";
+
+                   SumofPointsEuclideanDistance = SumofPointsEuclideanDistance +  std::pow(PointEuclideanDistance,2);
+
+
+
+
+
+
+                   error << std::endl;
+
+
+
+
+
+              } //end loop via points
+
+              error << std::endl;
+              SumofPointsEuclideanDistance = std::sqrt(SumofPointsEuclideanDistance/nImagePoints);
+              error << "Fiducial Registration Error:   " <<  SumofPointsEuclideanDistance;
+              error.close();
+
+            } //end if (args_info.error_arg != null)
+        } //end Dimension == 2 case
         else if (Dimension_temp==3)
         {
 
-            std::cout << "to ja";
             const     unsigned int   Dimension = 3;
             typedef   itk::Image< PixelType, Dimension > ImageType;
             typedef   float          VectorComponentType;
@@ -183,25 +304,25 @@ int main(int argc, char * argv[])
             typedef LandmarkBasedTransformInitializerType::LandmarkPointContainer     LandmarkContainerType;
             typedef LandmarkBasedTransformInitializerType::LandmarkPointType          LandmarkPointType;
 
-            //bufory na zaczyt wspolrzedncyh
             LandmarkContainerType imageLandmarks;
             LandmarkContainerType trackerLandmarks;
+            LandmarkContainerType imageLandmarksTotrackerCoordinate;
 
-            //bufory na pojedyncze punkty
-            LandmarkPointType imagePoint; // dane z CT
-            LandmarkPointType trackerPoint; //dane z tracking system
+
+            LandmarkPointType imagePoint;
+            LandmarkPointType trackerPoint;
 
             is >> x;
 
             while (is && !is.eof()) {
-                trackerPoint[0] = x;
-                is >> trackerPoint[1];
-                is >> trackerPoint[2];
-
-                is >> imagePoint[0];
+                imagePoint[0] = x;
                 is >> imagePoint[1];
                 is >> imagePoint[2];
 
+                is >> trackerPoint[0];
+                is >> trackerPoint[1];
+                is >> trackerPoint[2];
+
                 imageLandmarks.push_back(imagePoint );
                 trackerLandmarks.push_back(trackerPoint );
 
@@ -217,7 +338,6 @@ int main(int argc, char * argv[])
 
             transform->SetIdentity();
 
-            // Versor rotacji i macierz translacji
             landmarkBasedTransformInitializer->SetTransform(transform);
             landmarkBasedTransformInitializer->InitializeTransform();
 
@@ -235,8 +355,136 @@ int main(int argc, char * argv[])
               out.close();
 
 
+              if (args_info.error_arg != 0)
+              {
+              // Write result
+                std::ofstream error;
+                clitk::openFileForWriting(error, args_info.error_arg);
 
-        }
+
+              //imagePoints
+              int nImagePoints = imageLandmarks.size();
+              double  imagePoints[4][nImagePoints];
+              double  trackerPoints[4][nImagePoints];
+              double  imagePointstoTrackerCoordinate[4][nImagePoints];
+              double FRE = 0;
+              double PointEuclideanDistance = 0;
+              double SumofPointsEuclideanDistance = 0;
+
+
+              for(int n= 0; n < imageLandmarks.size(); n++)
+              {
+                  for(int j= 0; j < 3; j++)
+                  {
+                      imagePoints[j][n] = imageLandmarks[n][j];
+                      trackerPoints[j][n] = trackerLandmarks[n][j];
+
+
+
+
+//                      imagePoints[j][n] = trackerLandmarks[n][j];
+//                      trackerPoints[j][n] = imageLandmarks[n][j];
+                      imagePointstoTrackerCoordinate[j][n] = 0;
+
+                  }
+                  imagePoints[3][n] = 1.0;
+                  trackerPoints[3][n] = 1.0;
+                  imagePointstoTrackerCoordinate[3][n] = 0;
+              }
+
+
+              //transformation matrix
+              double RigidTransformationImagetoTracker[4][4];
+
+
+              for(int i = 0; i < 3; i++)
+                for(int j = 0; j < 3; j++)
+                    RigidTransformationImagetoTracker[i][j] = matrix[i][j];
+
+
+              for(int i = 0; i < 3; i++)
+                  RigidTransformationImagetoTracker[i][3] = offset[i];
+
+              RigidTransformationImagetoTracker[3][0] = 0.0;
+              RigidTransformationImagetoTracker[3][1] = 0.0;
+              RigidTransformationImagetoTracker[3][2] = 0.0;
+              RigidTransformationImagetoTracker[3][3] = 1.0;
+
+              //writing results
+              error << "ImagePoints:    CorrespondingTrackerPoints:     ImagePointstoTrackerCoordinateSystem:  PointsDistanceafterRigidRegistration:" << std::endl;
+
+
+
+
+
+              //Calculate FRE
+
+              for(int n= 0; n < imageLandmarks.size(); n++)
+              {
+
+
+                   //Calculate imagePointstoTrackerCoordinate
+                   for(int i = 0; i < 4; i++)
+                   {
+
+                       for(int k = 0; k < 4; k++)
+                       {
+                            imagePointstoTrackerCoordinate[i][n] = imagePointstoTrackerCoordinate[i][n] + RigidTransformationImagetoTracker[i][k]*imagePoints[k][n];
+
+                       }
+
+                   }
+
+
+                   //writing results
+
+                   for(int i = 0; i < 3; i++)
+                            error << imagePoints[i][n] << ' ';
+
+                   error << "    ";
+
+                   for(int i = 0; i < 3; i++)
+                            error << trackerPoints[i][n] << ' ';
+
+                   error << "   ";
+
+                   for(int i = 0; i < 3; i++)
+                            error << imagePointstoTrackerCoordinate[i][n] << ' ';
+
+
+
+                   //calculate PointEuclideanDistance
+                   PointEuclideanDistance = 0;
+                   for(int i=0; i < 3;i++)
+                   {
+                       PointEuclideanDistance = PointEuclideanDistance + std::pow((trackerPoints[i][n] - imagePointstoTrackerCoordinate[i][n]),2);
+                   }
+                   PointEuclideanDistance = std::sqrt(PointEuclideanDistance);
+                   error << "   " <<   PointEuclideanDistance   << "   ";
+
+                   SumofPointsEuclideanDistance = SumofPointsEuclideanDistance +  std::pow(PointEuclideanDistance,2);
+
+
+
+
+
+
+                   error << std::endl;
+
+
+
+
+
+              } //end loop via points
+
+              error << std::endl;
+              SumofPointsEuclideanDistance = std::sqrt(SumofPointsEuclideanDistance/nImagePoints);
+              error << "Fiducial Registration Error:   " <<  SumofPointsEuclideanDistance;
+              error.close();
+
+
+              } //end if (args_info.error_arg != null)
+          } //end Dimension == 3 case
         else
         {
               is.close();
@@ -248,6 +496,8 @@ int main(int argc, char * argv[])
 
 
 
+
+
   return EXIT_SUCCESS;
 }// end main
 //--------------------------------------------------------------------