X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkPointRigidRegistration.cxx;h=065b7d904847ba3081f054505ab1e68e9b092b5a;hb=15ebed734d7ef4727dff24b21d00d929b9162f90;hp=e66edd25a70abf6a3d7b842c594207f7dd5b3c72;hpb=54e318c028d7dc8441a263ee71b6757e03d57ecf;p=clitk.git diff --git a/tools/clitkPointRigidRegistration.cxx b/tools/clitkPointRigidRegistration.cxx index e66edd2..065b7d9 100644 --- a/tools/clitkPointRigidRegistration.cxx +++ b/tools/clitkPointRigidRegistration.cxx @@ -27,7 +27,7 @@ // clitk #include "clitkPointRigidRegistration_ggo.h" -#include "clitkPointRigidRegistrationGenericFilter.h" + //paste from RigidRegistration #include "itkImageFileReader.h" @@ -40,10 +40,10 @@ #include "itkVersorRigid3DTransform.h" #include -//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 //--------------------------------------------------------------------