PackRecalage/src/bbPackRecalageTransform3D3PointsBox.cxx

   1 /*
   2 # ---------------------------------------------------------------------
   3 #
   4 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
   5 #                        pour la Santé)
   6 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
   7 #
   8 #  This software is governed by the CeCILL-B license under French law and
   9 #  abiding by the rules of distribution of free software. You can  use,
  10 #  modify and/ or redistribute the software under the terms of the CeCILL-B
  11 #  license as circulated by CEA, CNRS and INRIA at the following URL
  12 #  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  13 #  or in the file LICENSE.txt.
  14 #
  15 #  As a counterpart to the access to the source code and  rights to copy,
  16 #  modify and redistribute granted by the license, users are provided only
  17 #  with a limited warranty  and the software's author,  the holder of the
  18 #  economic rights,  and the successive licensors  have only  limited
  19 #  liability.
  20 #
  21 #  The fact that you are presently reading this means that you have had
  22 #  knowledge of the CeCILL-B license and that you accept its terms.
  23 # ------------------------------------------------------------------------
  24 */
  25
  26
  27 #include "bbPackRecalageTransform3D3PointsBox.h"
  28 #include "bbPackRecalagePackage.h"
  29 namespace bbPackRecalage
  30 {
  31
  32 BBTK_ADD_BLACK_BOX_TO_PACKAGE(PackRecalage,Transform3D3PointsBox)
  33 BBTK_BLACK_BOX_IMPLEMENTATION(Transform3D3PointsBox,bbtk::AtomicBlackBox);
  34 void Transform3D3PointsBox::Process()
  35 {
  36         if(!bbGetInputInX1().empty() && bbGetInputInX1().size() >= 3 && !bbGetInputInX2().empty() && bbGetInputInX2().size() >= 3)
  37         {
  38                 double A1[3];
  39                 double B1[3];
  40                 double O1[3];
  41
  42                 double A2[3];
  43                 double B2[3];
  44                 double O2[3];
  45
  46                 //In case there is an order in the points
  47                 if(!bbGetInputLabels1().empty() && !bbGetInputLabels2().empty() && bbGetInputLabels1()[0].compare("1") == 0 && bbGetInputLabels1()[1].compare("2") == 0 && bbGetInputLabels1()[2].compare("3") == 0)
  48                 {
  49                         int i;
  50                         for(i = 0; i < 3; i++)
  51                         {
  52                                 if(bbGetInputLabels1()[i].compare("1") == 0)
  53                                 {
  54                                         O1[0] = bbGetInputInX1()[i];
  55                                         O1[1] = bbGetInputInY1()[i];
  56                                         O1[2] = bbGetInputInZ1()[i];
  57                                 }
  58                                 else if(bbGetInputLabels1()[i].compare("2") == 0)
  59                                 {
  60                                         A1[0] = bbGetInputInX1()[i];
  61                                         A1[1] = bbGetInputInY1()[i];
  62                                         A1[2] = bbGetInputInZ1()[i];
  63                                 }
  64                                 else if(bbGetInputLabels1()[i].compare("3") == 0)
  65                                 {
  66                                         B1[0] = bbGetInputInX1()[i];
  67                                         B1[1] = bbGetInputInY1()[i];
  68                                         B1[2] = bbGetInputInZ1()[i];
  69                                 }
  70                         }
  71                         for(i = 0; i < 3; i++)
  72                         {
  73                                 if(bbGetInputLabels2()[i].compare("1") == 0)
  74                                 {
  75                                         O2[0] = bbGetInputInX2()[i];
  76                                         O2[1] = bbGetInputInY2()[i];
  77                                         O2[2] = bbGetInputInZ2()[i];
  78                                 }
  79                                 else if(bbGetInputLabels2()[i].compare("2") == 0)
  80                                 {
  81                                         A2[0] = bbGetInputInX2()[i];
  82                                         A2[1] = bbGetInputInY2()[i];
  83                                         A2[2] = bbGetInputInZ2()[i];
  84                                 }
  85                                 else if(bbGetInputLabels2()[i].compare("3") == 0)
  86                                 {
  87                                         B2[0] = bbGetInputInX2()[i];
  88                                         B2[1] = bbGetInputInY2()[i];
  89                                         B2[2] = bbGetInputInZ2()[i];
  90                                 }
  91                         }
  92                 }
  93                 else
  94                 {
  95                         O1[0] = bbGetInputInX1()[0];
  96                         O1[1] = bbGetInputInY1()[0];
  97                         O1[2] = bbGetInputInZ1()[0];
  98
  99                         A1[0] = bbGetInputInX1()[1];
 100                         A1[1] = bbGetInputInY1()[1];
 101                         A1[2] = bbGetInputInZ1()[1];
 102
 103                         B1[0] = bbGetInputInX1()[2];
 104                         B1[1] = bbGetInputInY1()[2];
 105                         B1[2] = bbGetInputInZ1()[2];
 106
 107                         O2[0] = bbGetInputInX2()[0];
 108                         O2[1] = bbGetInputInY2()[0];
 109                         O2[2] = bbGetInputInZ2()[0];
 110
 111                         A2[0] = bbGetInputInX2()[1];
 112                         A2[1] = bbGetInputInY2()[1];
 113                         A2[2] = bbGetInputInZ2()[1];
 114
 115                         B2[0] = bbGetInputInX2()[2];
 116                         B2[1] = bbGetInputInY2()[2];
 117                         B2[2] = bbGetInputInZ2()[2];
 118                 }
 119
 120                 //Create the new center point
 121                 double centerPoint[3];  // JPR
 122                 centerPoint[0] = O2[0];
 123                 centerPoint[1] = O2[1];
 124                 centerPoint[2] = O2[2];
 125
 126                 transformer->SetSecondTranslation(centerPoint);
 127
 128                 centerPoint[0] = O1[0];
 129                 centerPoint[1] = O1[1];
 130                 centerPoint[2] = O1[2];
 131
 132                 transformer->SetFirstTranslation(centerPoint);
 133
 134                 //Create the first vector
 135                 double* vector1;
 136                 vector1 = planes->makeVector(O1, A1);
 137
 138                 //create the second vector
 139                 double* vector2;
 140                 vector2 = planes->makeVector(O2, A2);
 141
 142                 //The rotation axis
 143                 double* axis = planes->getCrossProduct(vector2, vector1);
 144
 145                 //Normalize the axis
 146                 axis = planes->getNormal(axis);
 147
 148                 //Normalize the 2 vectors for calculating the dot product
 149                 double* vector1N = planes->getNormal(vector1);
 150                 double* vector2N = planes->getNormal(vector2);
 151
 152                 //Dot product
 153                 double angle = planes->getDotProduct(vector2N, vector1N);
 154
 155                 //Convert the dot product to radians
 156                 angle = acos(angle);
 157
 158                 //Convert from Radians to Degrees (necesary for the transformation)
 159                 angle = vtkMath::DegreesFromRadians(angle);
 160
 161                 //Set the cross product and dot product of the first rotation
 162                 transformer->SetRotationAxis(axis);
 163                 transformer->SetAngle(angle);
 164
 165                 //Sets the second rotation axis (defined by the first vector) and the second angle as input
 166                 transformer->SetSecondRotationAxis(vector1);
 167
 168                 ///////////////////////////////
 169                 //Second Rotation calculations start now (adjust the angle between the planes)
 170                 /////////////////////////////////
 171
 172                 //Obtain the resultant transformation without moving the image first
 173                 vtkMatrix4x4* matrix = transformer->GetFirstResult()->GetMatrix();
 174
 175                 double newMatrix[3][3];
 176
 177                 int i,j;
 178                 for(i = 0; i < 3; i ++)
 179                 {
 180                         for(j = 0; j < 3; j ++)
 181                         {
 182                                 newMatrix[i][j] = matrix->GetElement(i,j);
 183                         }
 184                 }
 185
 186                 //Calculates the new position of the points after the first transformation
 187                 double A2N[3];
 188                 A2N[0] = A2[0]*newMatrix[0][0] + A2[0]*newMatrix[0][1] + A2[0]*newMatrix[0][2];
 189                 A2N[1] = A2[1]*newMatrix[1][0] + A2[0]*newMatrix[1][1] + A2[0]*newMatrix[1][2];
 190                 A2N[2] = A2[2]*newMatrix[2][0] + A2[0]*newMatrix[2][1] + A2[0]*newMatrix[2][2];
 191
 192                 double B2N[3];
 193                 B2N[0] = B2[0]*newMatrix[0][0] + B2[0]*newMatrix[0][1] + B2[0]*newMatrix[0][2];
 194                 B2N[1] = B2[1]*newMatrix[1][0] + B2[0]*newMatrix[1][1] + B2[0]*newMatrix[1][2];
 195                 B2N[2] = B2[2]*newMatrix[2][0] + B2[0]*newMatrix[2][1] + B2[0]*newMatrix[2][2];
 196
 197                 double O2N[3];
 198                 O2N[0] = O2[0]*newMatrix[0][0] + O2[0]*newMatrix[0][1] + O2[0]*newMatrix[0][2];
 199                 O2N[1] = O2[1]*newMatrix[1][0] + O2[0]*newMatrix[1][1] + O2[0]*newMatrix[1][2];
 200                 O2N[2] = O2[2]*newMatrix[2][0] + O2[0]*newMatrix[2][1] + O2[0]*newMatrix[2][2];
 201
 202                 ///////////////////////////////////////////////////////////////////////////////
 203                 //Creation of the two planes for the second transformation
 204                 ///////////////////////////////////////////////////////////////////////////////
 205
 206                 double * vector1A;
 207                 vector1A = planes->makeVector(O1, A1);
 208
 209                 double * vector1B;
 210                 vector1B = planes->makeVector(O1, B1);
 211
 212                 double * vector2AN;
 213                 vector2AN = planes->makeVector(O2N, A2N);
 214
 215                 double * vector2BN;
 216                 vector2BN = planes->makeVector(O2N, B2N);
 217
 218                 //Normal Plane 1
 219                 double* crossPlaneA = planes->getCrossProduct(vector1A, vector1B);
 220
 221                 //Normal Plane 2
 222                 double* crossPlaneB = planes->getCrossProduct(vector2AN, vector2BN);
 223
 224                 //Normalize the vectors
 225                 crossPlaneA = planes->getNormal(crossPlaneA);
 226                 crossPlaneB = planes->getNormal(crossPlaneB);
 227
 228                 //Angle between the planes
 229                 double anglePlanes = planes->getDotProduct(crossPlaneB, crossPlaneA);
 230
 231                 //Make necessary changes to the angle
 232                 anglePlanes = acos(anglePlanes);
 233
 234                 anglePlanes = vtkMath::DegreesFromRadians(anglePlanes);
 235
 236                 transformer->SetSecondAngle(-anglePlanes);
 237
 238                 // The calculation of the transformations are made
 239                 transformer->Run();
 240
 241                 // We get the results of transformer and set it as result of this box
 242                 bbSetOutputOut( transformer->GetResult() );
 243         }
 244         else
 245         {
 246                 bbSetOutputOut( NULL );
 247         }
 248 }
 249 void Transform3D3PointsBox::bbUserSetDefaultValues()
 250 {
 251         std::vector<int> empty;
 252         std::vector<std::string> emptyString;
 253         bbSetInputInX1(empty);
 254         bbSetInputInY1(empty);
 255         bbSetInputInZ1(empty);
 256         bbSetInputInX2(empty);
 257         bbSetInputInY2(empty);
 258         bbSetInputInZ2(empty);
 259         bbSetInputLabels1(emptyString);
 260         bbSetInputLabels2(emptyString);
 261         bbSetOutputOut(NULL);
 262 }
 263 void Transform3D3PointsBox::bbUserInitializeProcessing()
 264 {
 265         //We initialize the transformer
 266         transformer=new Transformer3D();
 267
 268         //We initialize the plane operator
 269         planes=new PlanesOperations();
 270 }
 271 void Transform3D3PointsBox::bbUserFinalizeProcessing()
 272 {
 273         //We delete the transformer
 274         delete transformer;
 275
 276         //We delete the plane operator
 277         delete planes;
 278 }
 279 }
 280 // EO namespace bbPackRecalage
 281