-
Notifications
You must be signed in to change notification settings - Fork 0
/
Program.cs
358 lines (289 loc) · 18.8 KB
/
Program.cs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
using System;
using System.Collections.Generic;
using System.Linq;
using System.Threading.Tasks;
using System.Windows.Forms;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
namespace WindowsFormsApp1
{
static class Globals
{
// global int
/*
public static string path_nf;
public static string path_f;
public static string path_tra;
public static Matrix<double> matrix_nf;
public static Matrix<double> matrix_f;
public static Matrix<double> matrix_tra;
*/
}
public class Transform3D
{
#region Fields
public Matrix<double> actualsMatrix;
public Matrix<double> nominalsMatrix;
#endregion
#region Constructors
/// <summary>
/// Creates Transform3D object from two arrays, converts to Math.NET matrices and assigns to actualsMatrix and nominalsMatrix
/// </summary>
/// <param name="actuals"></param>
/// <param name="nominals"></param>
public Transform3D(double[,] actuals, double[,] nominals)
{
//convert array to Math.NET Matrices
this.actualsMatrix = Matrix<double>.Build.DenseOfArray(actuals);
this.nominalsMatrix = Matrix<double>.Build.DenseOfArray(nominals);
}
/// <summary>
/// Creates Transform3D object from two Math.NET matrices and assigns to actualsMatrix and nominalsMatrix
/// </summary>
/// <param name="actuals"></param>
/// <param name="nominals"></param>
public Transform3D(Matrix<double> actuals, Matrix<double> nominals)
{
//convert array to Math.NET Matrices
this.actualsMatrix = actuals;
this.nominalsMatrix = nominals;
}
#endregion
#region Properties
/// <summary>
/// Gets and sets the 4X4 transform matrix (of type double[,])
/// </summary>
public double[,] TransformMatrix { get; set; }
/// <summary>
/// Gets and sets the Transform Vector of the form [XTrans, YTrans, ZTrans, ZRot, YRot, XRot] with rotations in radians
/// </summary>
public double[] Transform6DOFVector { get; set; }
public double[] Transform6DOFVectorSeparate { get; set; }
/// <summary>
/// Gets and sets the Transform Vector in Siemens Convention of the form [XTransSiemens, YTransSiemens, ZTransSiemens, ZRotSiemens, YRotSiemens, XRotSiemens] with rotations in radians
/// </summary>
public double[] TransformSiemens6DOFVector { get; set; }
public double[] TransformSiemens6DOFVectorSeparate { get; set; }
/// <summary>
/// Gets and sets the 4X4 Rotation matrix following the conventional matrix multiplication order -> Rz*Ry*Rx
/// returns as an array of type double[,]
/// </summary>
public double[,] RotationMatrix { get; set; }
/// <summary>
/// Gets and sets the 4X4 Translation matrix (of type double[,])
/// </summary>
public double[,] TranslationMatrix { get; set; }
/// <summary>
/// Gets and sets the RMS of the error between nominal points and best-fit-transformed actual points
/// </summary>
public double ErrorRMS { get; set; }
/// <summary>
/// Gets and sets the determinant of the rotation matrix to determine presence of reflection between point sets (-1 -> reflection present, 1-> no reflection present)
/// </summary>
public double RotMatrixDeterminant { get; set; }
#endregion
#region Public methods
/// <summary>
/// Calculates the 4x4 transformation from actualsMatrix to nominalsMatrix by the Singular Value Decomposition method as well as its component Rotation/Translation
/// matrices, transform fit error, Rotation matrix determinant and Transform vector od the form [XTrans, YTrans, ZTrans, ZRot, YRot, XRot]
/// /// </summary>
/// <param name="actualsMatrix"></param>
/// <param name="nominalsMatrix"></param>
/// <returns></returns>
public bool CalcTransform(Matrix<double> actualsMatrix, Matrix<double> nominalsMatrix)
{
//check to see if actuals and nominals matrices are same dimension and that they have only 3 columns
if (actualsMatrix.ColumnCount != nominalsMatrix.ColumnCount || actualsMatrix.RowCount != nominalsMatrix.RowCount || actualsMatrix.ColumnCount != 3)
{
if (actualsMatrix.ColumnCount != 3)
{
Console.WriteLine("Incoming arrays have more than 3 columns (X,Y,Z)");
}
else
{
Console.WriteLine("Dimensions of actuals and nominals don't match.");
}
//throw new ArgumentException(string.Format("Check that incomiang arrays have matching dimensions and only 3 columns."));
return false;
}
else
{
//find centroids of pointCloudStart and pointCloudFinish
Matrix<double> actualsCentroid = GetCentroid(actualsMatrix).ToRowMatrix();
Matrix<double> nominalsCentroid = GetCentroid(nominalsMatrix).ToRowMatrix();
//make centroid matrix with centroid repeated for every row of actuals/nominalsMatrix
Matrix<double> actualsCentroidNxM = actualsCentroid;
Matrix<double> nominalsCentroidNxM = nominalsCentroid;
int count = 0;
while (count < actualsMatrix.RowCount - 1)
{
actualsCentroidNxM = actualsCentroidNxM.Stack(actualsCentroid);
nominalsCentroidNxM = nominalsCentroidNxM.Stack(nominalsCentroid);
count = count + 1;
}
//center pointCloudStart/pointCloudFinishMatrix on centroids
Matrix<double> actualsMatrixCentered = actualsMatrix - actualsCentroidNxM;
Matrix<double> nominalsMatrixCentered = nominalsMatrix - nominalsCentroidNxM;
//calculate covariance matrix
Matrix<double> covarianceMatrix = actualsMatrixCentered.Transpose() * nominalsMatrixCentered;
//singular value decomposition: http://en.wikipedia.org/wiki/Singular_value_decomposition
var svd = covarianceMatrix.Svd();
Matrix<double> svdU = svd.U;
Matrix<double> svdS = svd.S.ToRowMatrix();
Matrix<double> svdVt = svd.VT;
//calculate rotation matrix of the form -> nominal_i = RotationMatrix * actual_i + TranslationComponentVector
Matrix<double> rotationMatrix = svdVt.Transpose() * svdU.Transpose();
//calclate determinant of rotation matrix
RotMatrixDeterminant = rotationMatrix.Determinant();
//Check rotation matrix for reflection
if (RotMatrixDeterminant < 0)
{
Console.WriteLine("Refelction Detected");
//If refelecction exists, multiply by correction matrix such that reflectionCorrection is similar to an Identity Matrix but with value [N,M] = -1
//See equation 23,24 at http://igl.ethz.ch/projects/ARAP/svd_rot.pdf
Matrix<double> reflectionCorrection = new DiagonalMatrix(3, 3, 1.0);
reflectionCorrection[2, 2] = -1.0;
//Python (not Matlab) code from http://nghiaho.com/?page_id=671 contridicts with equation 23,24 at http://igl.ethz.ch/projects/ARAP/svd_rot.pdf
//svdVt = svdVt.InsertRow(3, svdVt.Row(2) * -1).RemoveRow(2);
//recalculate rotation matrix
rotationMatrix = svdVt.Transpose() * reflectionCorrection * svdU.Transpose();
}
//calculate translation matrix of the form -> nominal_i = RotationMatrix * actual_i + TranslationComponentVector
Matrix<double> TranslationComponentVector = -rotationMatrix * actualsCentroid.Transpose() + nominalsCentroid.Transpose();
//Transform matrix comprised of Rotation and Translation matrices such that -> nominal_i = RotationMatrix * actual_i + TranslationComponentVector
//append translation to rotation and append last row (0,0,0,1) to make to matrix suitable for matrix multiplication (4x4)
Matrix<double> transformMatrix = rotationMatrix.Append(TranslationComponentVector);
transformMatrix = transformMatrix.InsertRow(3, new DenseVector(new[] { 0.0, 0.0, 0.0, 1.0 }));
////append 4th column (0,0,0) and 4th row (0,0,0,1) to make to matrix suitable for matrix multiplication (4x4) multiply rotation and translation matrices
rotationMatrix = rotationMatrix.InsertColumn(3, new DenseVector(new[] { 0.0, 0.0, 0.0 }));
rotationMatrix = rotationMatrix.InsertRow(3, new DenseVector(new[] { 0.0, 0.0, 0.0, 1.0 }));
//Build Translation Matrix suitable for matrix multiplication (of the form transformMatrix = rotationMatrix * translationMatrix) from TranslationComponentVector
Matrix<double> translationMatrix = rotationMatrix.Inverse().Multiply(transformMatrix);
//Error Calculation
//Calcualte actuals' = tranformMatrix * actualsMatrix (made possible by first appending a column matrix of 1's to actuals)
Matrix<double> columnVectorof1s = new DenseMatrix(actualsMatrix.RowCount, 1) + 1;
Matrix<double> actualsMatrixNx4 = actualsMatrix.InsertColumn(3, columnVectorof1s.Column(0));
Matrix<double> actualsMatrixPrime = (transformMatrix * actualsMatrixNx4.Transpose()).Transpose();
Matrix<double> errorMatrix = actualsMatrixPrime.RemoveColumn(3) - nominalsMatrix;
Matrix<double> errorMatrixSquared = errorMatrix.PointwiseMultiply(errorMatrix);
double errorMatrixSum = errorMatrixSquared.ColumnAbsoluteSums().Sum();
//Method outputs
RotationMatrix = rotationMatrix.ToArray();
TranslationMatrix = translationMatrix.ToArray();
ErrorRMS = Math.Sqrt(errorMatrixSum / actualsMatrix.RowCount);
TransformMatrix = transformMatrix.ToArray();
Transform6DOFVector = TransformMatrixTo6DOFVector(transformMatrix);
Transform6DOFVectorSeparate = TransformMatrixTo6DOFVectorSeparate(rotationMatrix, translationMatrix);
TransformSiemens6DOFVector = TransformMatrixToSiemens6DOFVector(transformMatrix);
TransformSiemens6DOFVectorSeparate = TransformMatrixToSiemens6DOFVectorSeparate(rotationMatrix, translationMatrix);
return true;
}
}
#endregion
#region Private methods
/// <summary>
/// Returns the centroid of an NxM matrix
/// </summary>
/// <param name="matrix"></param>
/// <returns></returns>
private Vector<double> GetCentroid(Matrix<double> matrix)
{
return matrix.ColumnSums() / matrix.RowCount;
}
/// <summary>
/// Returns the 6DoF transform vector of the form [XTrans, YTrans, ZTrans, ZRot, YRot, XRot] with rotations in radians
/// </summary>
/// <param name="TransformMatrix"></param>
/// <returns></returns>
public double[] TransformMatrixTo6DOFVector(Matrix<double> TransformMatrix)
{
Vector<double> tempVector = new DenseVector(6);
//Set the X,Y, and Z translation components of the transform vector
tempVector[0] = TransformMatrix[0, 3]; //X translation
tempVector[1] = TransformMatrix[1, 3]; //Y translation
tempVector[2] = TransformMatrix[2, 3]; //Z translation
//Set the Z, Y, and X rotation components of the transform vector - see equation derivation at page 5 of http://www.usna.edu/Users/cs/taylor/courses/si475/class/3dkinematics.pdf
double beta = Math.Atan2(-TransformMatrix[2, 0], Math.Sqrt(Math.Pow(TransformMatrix[0, 0], 2) + Math.Pow(TransformMatrix[1, 0], 2))); //Y rotation - beta
double gamma = Math.Atan2(TransformMatrix[1, 0] / Math.Cos(beta), TransformMatrix[0, 0] / Math.Cos(beta)); //Z rotation - gamma
double alpha = Math.Atan2(TransformMatrix[2, 1] / Math.Cos(beta), TransformMatrix[2, 2] / Math.Cos(beta)); //X rotation - alpha
tempVector[3] = gamma; //Z rotation - kappa
tempVector[4] = beta; //Y rotation - phi
tempVector[5] = alpha; //X rotation - omega
return tempVector.ToArray();
}
private double[] TransformMatrixTo6DOFVectorSeparate(Matrix<double> RotationMatrix, Matrix<double> TranslationMatrix)
{
Vector<double> tempVector = new DenseVector(6);
//Set the X,Y, and Z translation components of the transform vector
tempVector[0] = TranslationMatrix[0, 3]; //X translation
tempVector[1] = TranslationMatrix[1, 3]; //Y translation
tempVector[2] = TranslationMatrix[2, 3]; //Z translation
//Set the Z, Y, and X rotation components of the transform vector - see equation derivation at page 5 of http://www.usna.edu/Users/cs/taylor/courses/si475/class/3dkinematics.pdf
double beta = Math.Atan2(-RotationMatrix[2, 0], Math.Sqrt(Math.Pow(RotationMatrix[0, 0], 2) + Math.Pow(RotationMatrix[1, 0], 2))); //Y rotation - beta
double gamma = Math.Atan2(RotationMatrix[1, 0] / Math.Cos(beta), RotationMatrix[0, 0] / Math.Cos(beta)); //Z rotation - gamma
double alpha = Math.Atan2(RotationMatrix[2, 1] / Math.Cos(beta), RotationMatrix[2, 2] / Math.Cos(beta)); //X rotation - alpha
tempVector[3] = gamma; //Z rotation - kappa
tempVector[4] = beta; //Y rotation - phi
tempVector[5] = alpha; //X rotation - omega
return tempVector.ToArray();
}
/// <summary>
/// Returns the 6DoF transform vector in Siemens Convention of the form [XTransSiemens, YTransSiemens, ZTransSiemens, ZRotSiemens, YRotSiemens, XRotSiemens] with rotations in radians
/// </summary>
/// <param name="RotationMatrix"></param>
/// <param name="TranslationMatrix"></param>
/// <returns></returns>
private double[] TransformMatrixToSiemens6DOFVector(Matrix<double> TransformMatrix)
{
Vector<double> tempVector = new DenseVector(6);
//Set the X,Y, and Z translation components of the transform vector for Siemens convention
//Translation components in Siemens convention is simply the negative of the typical translation components
tempVector[0] = -TransformMatrix[0, 3]; //X translation
tempVector[1] = -TransformMatrix[1, 3]; //Y translation
tempVector[2] = -TransformMatrix[2, 3]; //Z translation
//Set the Z, Y, and X rotation components of the transform vector - see equation derivation at page 5 of http://www.usna.edu/Users/cs/taylor/courses/si475/class/3dkinematics.pdf
//Siemens convention was determined experimentally. To find the typical 3X3 rotation matrix (RotationMatrix in this class), the following matrix multiplication order is used -> Rz*Ry*Rx
//To achieve the same rotation matrix with angles provided by Siemens MEAFrame, the following matrix multiplication must be used -> transpose(Rx)*transpose(Ry)*transpose(Rz)
double beta = Math.Atan2(-TransformMatrix[0, 2], Math.Sqrt(Math.Pow(TransformMatrix[0, 0], 2) + Math.Pow(TransformMatrix[0, 1], 2))); //Y rotation - beta
double gamma = Math.Atan2(TransformMatrix[0, 1] / Math.Cos(beta), TransformMatrix[0, 0] / Math.Cos(beta)); //Z rotation - gamma
double alpha = Math.Atan2(TransformMatrix[1, 2] / Math.Cos(beta), TransformMatrix[2, 2] / Math.Cos(beta)); //X rotation - alpha
tempVector[3] = gamma; //Z rotation - kappa
tempVector[4] = beta; //Y rotation - phi
tempVector[5] = alpha; //X rotation - omega
return tempVector.ToArray();
}
private double[] TransformMatrixToSiemens6DOFVectorSeparate(Matrix<double> RotationMatrix, Matrix<double> TranslationMatrix)
{
Vector<double> tempVector = new DenseVector(6);
//Set the X,Y, and Z translation components of the transform vector for Siemens convention
//Translation components in Siemens convention is simply the negative of the typical translation components
tempVector[0] = -TranslationMatrix[0, 3]; //X translation
tempVector[1] = -TranslationMatrix[1, 3]; //Y translation
tempVector[2] = -TranslationMatrix[2, 3]; //Z translation
//Set the Z, Y, and X rotation components of the transform vector - see equation derivation at page 5 of http://www.usna.edu/Users/cs/taylor/courses/si475/class/3dkinematics.pdf
//Siemens convention was determined experimentally. To find the typical 3X3 rotation matrix (RotationMatrix in this class), the following matrix multiplication order is used -> Rz*Ry*Rx
//To achieve the same rotation matrix with angles provided by Siemens MEAFrame, the following matrix multiplication must be used -> transpose(Rx)*transpose(Ry)*transpose(Rz)
double beta = Math.Atan2(-RotationMatrix[0, 2], Math.Sqrt(Math.Pow(RotationMatrix[0, 0], 2) + Math.Pow(RotationMatrix[0, 1], 2))); //Y rotation - beta
double gamma = Math.Atan2(RotationMatrix[0, 1] / Math.Cos(beta), RotationMatrix[0, 0] / Math.Cos(beta)); //Z rotation - gamma
double alpha = Math.Atan2(RotationMatrix[1, 2] / Math.Cos(beta), RotationMatrix[2, 2] / Math.Cos(beta)); //X rotation - alpha
tempVector[3] = gamma; //Z rotation - kappa
tempVector[4] = beta; //Y rotation - phi
tempVector[5] = alpha; //X rotation - omega
return tempVector.ToArray();
}
#endregion
}
static class Program
{
/// <summary>
/// The main entry point for the application.
/// </summary>
[STAThread]
static void Main()
{
Application.EnableVisualStyles();
Application.SetCompatibleTextRenderingDefault(false);
Application.Run(new Form1());
}
}
}