2 # ---------------------------------------------------------------------
4 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
6 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
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.
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
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 # ------------------------------------------------------------------------
27 #include "bbPackRecalageTransform3D3PointsBox.h"
28 #include "bbPackRecalagePackage.h"
29 namespace bbPackRecalage
32 BBTK_ADD_BLACK_BOX_TO_PACKAGE(PackRecalage,Transform3D3PointsBox)
33 BBTK_BLACK_BOX_IMPLEMENTATION(Transform3D3PointsBox,bbtk::AtomicBlackBox);
34 void Transform3D3PointsBox::Process()
36 if(!bbGetInputInX1().empty() && bbGetInputInX1().size() >= 3 && !bbGetInputInX2().empty() && bbGetInputInX2().size() >= 3)
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)
50 for(i = 0; i < 3; i++)
52 if(bbGetInputLabels1()[i].compare("1") == 0)
54 O1[0] = bbGetInputInX1()[i];
55 O1[1] = bbGetInputInY1()[i];
56 O1[2] = bbGetInputInZ1()[i];
58 else if(bbGetInputLabels1()[i].compare("2") == 0)
60 A1[0] = bbGetInputInX1()[i];
61 A1[1] = bbGetInputInY1()[i];
62 A1[2] = bbGetInputInZ1()[i];
64 else if(bbGetInputLabels1()[i].compare("3") == 0)
66 B1[0] = bbGetInputInX1()[i];
67 B1[1] = bbGetInputInY1()[i];
68 B1[2] = bbGetInputInZ1()[i];
71 for(i = 0; i < 3; i++)
73 if(bbGetInputLabels2()[i].compare("1") == 0)
75 O2[0] = bbGetInputInX2()[i];
76 O2[1] = bbGetInputInY2()[i];
77 O2[2] = bbGetInputInZ2()[i];
79 else if(bbGetInputLabels2()[i].compare("2") == 0)
81 A2[0] = bbGetInputInX2()[i];
82 A2[1] = bbGetInputInY2()[i];
83 A2[2] = bbGetInputInZ2()[i];
85 else if(bbGetInputLabels2()[i].compare("3") == 0)
87 B2[0] = bbGetInputInX2()[i];
88 B2[1] = bbGetInputInY2()[i];
89 B2[2] = bbGetInputInZ2()[i];
95 O1[0] = bbGetInputInX1()[0];
96 O1[1] = bbGetInputInY1()[0];
97 O1[2] = bbGetInputInZ1()[0];
99 A1[0] = bbGetInputInX1()[1];
100 A1[1] = bbGetInputInY1()[1];
101 A1[2] = bbGetInputInZ1()[1];
103 B1[0] = bbGetInputInX1()[2];
104 B1[1] = bbGetInputInY1()[2];
105 B1[2] = bbGetInputInZ1()[2];
107 O2[0] = bbGetInputInX2()[0];
108 O2[1] = bbGetInputInY2()[0];
109 O2[2] = bbGetInputInZ2()[0];
111 A2[0] = bbGetInputInX2()[1];
112 A2[1] = bbGetInputInY2()[1];
113 A2[2] = bbGetInputInZ2()[1];
115 B2[0] = bbGetInputInX2()[2];
116 B2[1] = bbGetInputInY2()[2];
117 B2[2] = bbGetInputInZ2()[2];
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];
126 transformer->SetSecondTranslation(centerPoint);
128 centerPoint[0] = O1[0];
129 centerPoint[1] = O1[1];
130 centerPoint[2] = O1[2];
132 transformer->SetFirstTranslation(centerPoint);
134 //Create the first vector
136 vector1 = planes->makeVector(O1, A1);
138 //create the second vector
140 vector2 = planes->makeVector(O2, A2);
143 double* axis = planes->getCrossProduct(vector2, vector1);
146 axis = planes->getNormal(axis);
148 //Normalize the 2 vectors for calculating the dot product
149 double* vector1N = planes->getNormal(vector1);
150 double* vector2N = planes->getNormal(vector2);
153 double angle = planes->getDotProduct(vector2N, vector1N);
155 //Convert the dot product to radians
158 //Convert from Radians to Degrees (necesary for the transformation)
159 angle = vtkMath::DegreesFromRadians(angle);
161 //Set the cross product and dot product of the first rotation
162 transformer->SetRotationAxis(axis);
163 transformer->SetAngle(angle);
165 //Sets the second rotation axis (defined by the first vector) and the second angle as input
166 transformer->SetSecondRotationAxis(vector1);
168 ///////////////////////////////
169 //Second Rotation calculations start now (adjust the angle between the planes)
170 /////////////////////////////////
172 //Obtain the resultant transformation without moving the image first
173 vtkMatrix4x4* matrix = transformer->GetFirstResult()->GetMatrix();
175 double newMatrix[3][3];
178 for(i = 0; i < 3; i ++)
180 for(j = 0; j < 3; j ++)
182 newMatrix[i][j] = matrix->GetElement(i,j);
186 //Calculates the new position of the points after the first transformation
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[1]*newMatrix[1][1] + A2[1]*newMatrix[1][2];
190 A2N[2] = A2[2]*newMatrix[2][0] + A2[2]*newMatrix[2][1] + A2[2]*newMatrix[2][2];
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[1]*newMatrix[1][1] + B2[1]*newMatrix[1][2];
195 B2N[2] = B2[2]*newMatrix[2][0] + B2[2]*newMatrix[2][1] + B2[2]*newMatrix[2][2];
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[1]*newMatrix[1][1] + O2[1]*newMatrix[1][2];
200 O2N[2] = O2[2]*newMatrix[2][0] + O2[2]*newMatrix[2][1] + O2[2]*newMatrix[2][2];
202 ///////////////////////////////////////////////////////////////////////////////
203 //Creation of the two planes for the second transformation
204 ///////////////////////////////////////////////////////////////////////////////
207 vector1A = planes->makeVector(O1, A1);
210 vector1B = planes->makeVector(O1, B1);
213 vector2AN = planes->makeVector(O2N, A2N);
216 vector2BN = planes->makeVector(O2N, B2N);
219 double* crossPlaneA = planes->getCrossProduct(vector1A, vector1B);
222 double* crossPlaneB = planes->getCrossProduct(vector2AN, vector2BN);
224 //Normalize the vectors
225 crossPlaneA = planes->getNormal(crossPlaneA);
226 crossPlaneB = planes->getNormal(crossPlaneB);
228 //Angle between the planes
229 double anglePlanes = planes->getDotProduct(crossPlaneB, crossPlaneA);
231 //Make necessary changes to the angle
232 anglePlanes = acos(anglePlanes);
234 anglePlanes = vtkMath::DegreesFromRadians(anglePlanes);
236 transformer->SetSecondAngle(-anglePlanes);
240 vector2B = planes->makeVector(O2, B2);
241 double mag1= planes->getMagnitud(vector1B);
242 double mag2= planes->getMagnitud(vector2B);
243 transformer->SetScale(mag1/mag2);
248 vector2B = planes->makeVector(O2, B2);
249 double mag1= planes->getMagnitud(vector1);
250 double mag2= planes->getMagnitud(vector2);
251 transformer->SetScale(mag2/mag1);
254 // The calculation of the transformations are made
257 // We get the results of transformer and set it as result of this box
258 bbSetOutputOut( transformer->GetResult() );
262 bbSetOutputOut( NULL );
265 void Transform3D3PointsBox::bbUserSetDefaultValues()
267 std::vector<int> empty;
268 std::vector<std::string> emptyString;
269 bbSetInputInX1(empty);
270 bbSetInputInY1(empty);
271 bbSetInputInZ1(empty);
272 bbSetInputInX2(empty);
273 bbSetInputInY2(empty);
274 bbSetInputInZ2(empty);
275 bbSetInputLabels1(emptyString);
276 bbSetInputLabels2(emptyString);
277 bbSetOutputOut(NULL);
279 void Transform3D3PointsBox::bbUserInitializeProcessing()
281 //We initialize the transformer
282 transformer=new Transformer3D();
284 //We initialize the plane operator
285 planes=new PlanesOperations();
287 void Transform3D3PointsBox::bbUserFinalizeProcessing()
289 //We delete the transformer
292 //We delete the plane operator
296 // EO namespace bbPackRecalage