71 (AAT && (pResultMatrix->
columns != pMatrix->
rows || pResultMatrix->
rows != pMatrix->
rows)) ||
75 printf(
"error: matrix dimensions do not match for LinearAlgebra::SelfProduct\n");
83 pMatrix = pTransposedMatrix;
86 const int columns = pMatrix->
columns;
87 const int rows = pMatrix->
rows;
90 float *result = pResultMatrix->
data;
94 for (
int i = 0, input_offset1 = 0; i < rows; i++, input_offset1 += columns)
96 for (
int j = i, input_offset2 = input_offset1; j < rows; j++, input_offset2 += columns)
98 const float *data1 = data + input_offset1;
99 const float *data2 = data + input_offset2;
102 for (
int k = 0; k < columns; k++)
103 sum += data1[k] * data2[k];
105 result[i * rows + j] = result[j * rows + i] = sum;
112 float *temp_data = temp.
data;
114 for (
int i = 0, input_offset1 = 0; i < rows; i++, input_offset1 += columns)
116 for (
int j = i, input_offset2 = input_offset1; j < rows; j++, input_offset2 += columns)
118 const float *data1 = data + input_offset1;
119 const float *data2 = data + input_offset2;
122 for (
int k = 0; k < columns; k++)
123 sum += data1[k] * data2[k];
125 temp_data[i * rows + j] = temp_data[j * rows + i] = sum;
129 memcpy(pResultMatrix->
data, temp_data, pResultMatrix->
rows * pResultMatrix->
columns *
sizeof(
float));
140 printf(
"error: matrices A, B, and pResultMatrix do not satisfy requirements for LinearAlgebra::Multiply\n");
146 printf(
"error: matrices A, B, and pResultMatrix do not satisfy requirements for LinearAlgebra::Multiply\n");
150 const float *data1 = A->
data;
151 const float *data2 = B->
data;
152 float *result = pResultMatrix->
data;
154 const int nSize = pResultMatrix->
columns * pResultMatrix->
rows;
158 const int rowsA = A->
rows;
159 const int rowsB = B->
rows;
160 const int columns = A->
columns;
162 if (data1 != result && data2 != result)
164 for (
int i = 0, input_offset1 = 0, output_offset = 0; i < rowsA; i++, input_offset1 += columns)
166 for (
int j = 0, input_offset2 = 0; j < rowsB; j++, input_offset2 += columns, output_offset++)
168 const float *data1_ = data1 + input_offset1;
169 const float *data2_ = data2 + input_offset2;
172 for (
int k = 0; k < columns; k++)
173 sum += data1_[k] * data2_[k];
175 result[output_offset] = sum;
182 float *temp_data = temp.
data;
186 for (l = 0; l < nSize; l++)
189 for (
int i = 0, input_offset1 = 0, output_offset = 0; i < rowsA; i++, input_offset1 += columns)
191 for (
int j = 0, input_offset2 = input_offset1; j < rowsB; j++, input_offset2 += columns, output_offset++)
193 const float *data1_ = data1 + input_offset1;
194 const float *data2_ = data2 + input_offset2;
197 for (
int k = 0; k < columns; k++)
198 sum += data1_[k] * data2_[k];
200 temp_data[output_offset] = sum;
204 for (l = 0; l < nSize; l++)
205 result[l] = temp_data[l];
210 const int i_max = A->
rows;
215 if (data1 != result && data2 != result)
217 for (
int l = 0; l < nSize; l++)
220 for (
int i = 0, input1_offset = 0; i < i_max; i++)
222 const int result_offset_start = i * j_max;
224 for (
int k = 0, input2_offset = 0; k < k_max; k++, input1_offset++)
225 for (
int j = 0, result_offset = result_offset_start; j < j_max; j++, result_offset++, input2_offset++)
226 result[result_offset] += data1[input1_offset] * data2[input2_offset];
232 float *temp_data = temp.
data;
236 for (l = 0; l < nSize; l++)
239 for (
int i = 0, input1_offset = 0; i < i_max; i++)
241 const int result_offset_start = i * j_max;
243 for (
int k = 0, input2_offset = 0; k < k_max; k++, input1_offset++)
244 for (
int j = 0, result_offset = result_offset_start; j < j_max; j++, result_offset++, input2_offset++)
245 temp_data[result_offset] += data1[input1_offset] * data2[input2_offset];
248 for (l = 0; l < nSize; l++)
249 result[l] = temp_data[l];
258 printf(
"error: input and output matrix do not match LinearAlgebra::Transpose\n");
262 const int rows = pMatrix->
rows;
263 const int columns = pMatrix->
columns;
265 float *output = pResultMatrix->
data;
269 for (
int i = 0, input_offset = 0; i < rows; i++)
270 for (
int j = 0, output_offset = i; j < columns; j++, input_offset++, output_offset += rows)
271 output[output_offset] = input[input_offset];
277 float *temp_data = temp.
data;
281 for (i = 0, input_offset = 0; i < rows; i++)
282 for (
int j = 0, output_offset = i; j < columns; j++, input_offset++, output_offset += rows)
283 temp_data[output_offset] = input[input_offset];
285 const int nSize = rows * columns;
287 for (i = 0; i < nSize; i++)
288 output[i] = temp_data[i];
296 printf(
"error: input matrix and vector and output vector do not match for CFloatMatrix::Multiply\n");
300 const int rows = pMatrix->
rows;
301 const int columns = pMatrix->
columns;
302 const float *data_m = pMatrix->
data;
303 const float *data_v = pVector->
data;
304 float *data_r = pResultVector->
data;
306 if (data_r != data_v)
308 for (
int i = 0,
offset = 0; i < rows; i++)
312 for (
int j = 0; j < columns; j++,
offset++)
313 sum += data_m[
offset] * data_v[j];
320 float *temp =
new float[rows];
323 for (i = 0, offset = 0; i < rows; i++)
327 for (
int j = 0; j < columns; j++, offset++)
328 sum += data_m[offset] * data_v[j];
333 for (i = 0; i < rows; i++)
342 MulMatVec(pMatrix, pVector1, pResultVector);
351 printf(
"error: matrix dimensions do not match in LinearAlgebra::AddMatMat\n");
355 const int nSize = pMatrix1->
rows * pMatrix1->
columns;
356 const float *data1 = pMatrix1->
data;
357 const float *data2 = pMatrix2->
data;
358 float *result = pResultMatrix->
data;
360 for (
int i = 0; i < nSize; i++)
361 result[i] = data1[i] + data2[i];
369 printf(
"error: matrix dimensions do not match in LinearAlgebra::SubtractMatMat\n");
373 const int nSize = pMatrix1->
rows * pMatrix1->
columns;
374 const float *data1 = pMatrix1->
data;
375 const float *data2 = pMatrix2->
data;
376 float *result = pResultMatrix->
data;
378 for (
int i = 0; i < nSize; i++)
379 result[i] = data1[i] - data2[i];
386 printf(
"error: vector dimensions do not match in LinearAlgebra::AddVecVec\n");
390 const int nDimension = pVector1->
dimension;
391 const float *data1 = pVector1->
data;
392 const float *data2 = pVector2->
data;
393 float *result = pResultVector->
data;
395 for (
int i = 0; i < nDimension; i++)
396 result[i] = data1[i] + data2[i];
403 printf(
"error: vector dimensions do not match in LinearAlgebra::SubtractVecVec\n");
407 const int nDimension = pVector1->
dimension;
408 const float *data1 = pVector1->
data;
409 const float *data2 = pVector2->
data;
410 float *result = pResultVector->
data;
412 for (
int i = 0; i < nDimension; i++)
413 result[i] = data1[i] - data2[i];
420 printf(
"error: matrix dimensions do not match in LinearAlgebra::AddToMat\n");
424 const int nSize = pMatrix->
rows * pMatrix->
columns;
426 const float *data_to_add = pMatrixToAdd->
data;
428 for (
int i = 0; i < nSize; i++)
429 data[i] += data_to_add[i];
436 printf(
"error: matrix dimensions do not match in LinearAlgebra::SubtractFromMat\n");
440 const int nSize = pMatrix->
rows * pMatrix->
columns;
442 const float *data_to_subtract = pMatrixToAdd->
data;
444 for (
int i = 0; i < nSize; i++)
445 data[i] -= data_to_subtract[i];
452 printf(
"error: vector dimensions do not match in LinearAlgebra::AddToVec\n");
456 const int nDimension = pVector->
dimension;
458 const float *data_to_add = pVectorToAdd->
data;
460 for (
int i = 0; i < nDimension; i++)
461 data[i] += data_to_add[i];
468 printf(
"error: vector dimensions do not match in LinearAlgebra::SubtractFromVec\n");
472 const int nDimension = pVector->
dimension;
474 const float *data_to_subtract = pVectorToSubtract->
data;
476 for (
int i = 0; i < nDimension; i++)
477 data[i] -= data_to_subtract[i];
486 const int columns = pMatrix->
columns;
487 const int rows = pMatrix->
rows;
488 const int max = rows * columns;
490 float *output = pResultMatrix->
data;
492 for (
int i = 0; i < columns; i++)
497 for (j = 0; j < max; j += columns)
502 for (j = 0; j < max; j += columns)
503 output[j] = data[j] - mean;
515 const int columns = pMatrix->
columns;
516 const int rows = pMatrix->
rows;
518 float *output = pResultMatrix->
data;
520 for (
int i = 0; i < rows; i++)
525 for (j = 0; j < columns; j++)
530 for (j = 0; j < columns; j++)
531 output[j] = data[j] - mean;
542 printf(
"error: A, x do not match LinearAlgebra::SolveLinearLeastSquaresHomogeneousSVD\n");
548 printf(
"error: equation system is underdetermined in LinearAlgebra::SolveLinearLeastSquaresHomogeneousSVD\n");
552 const int m = A->
rows;
557 SVD(A, &W, 0, &V,
false,
false,
false);
567 printf(
"error: A, b, x do not match LinearAlgebra::SolveLinearLeastSquaresSVD\n");
573 printf(
"error: equation system is underdetermined in LinearAlgebra::SolveLinearLeastSquaresSVD\n");
586 printf(
"error: A, b, x do not match LinearAlgebra::SolveLinearLeastSquaresSimple\n");
592 printf(
"error: equation system is underdetermined in LinearAlgebra::SolveLinearLeastSquaresSimple\n");
610 printf(
"error: input and output matrix do not match LinearAlgebra::CalculatePseudoInverseSVD\n");
620 const int m = pInputMatrix->
rows;
621 const int n = pInputMatrix->
columns;
626 SVD(A, &W, &UT, &V,
false,
true,
false);
631 const int min =
MY_MIN(m, n);
634 float fThreshold = 0.0f;
636 for (i = 0; i < min; i++)
637 fThreshold += WT(i, i);
639 fThreshold *= 2 * FLT_EPSILON;
642 for (i = 0; i < min; i++)
643 if (WT(i, i) < fThreshold)
646 WT(i, i) = 1.0f / WT(i, i);
658 printf(
"error: input and output matrix do not match LinearAlgebra::CalculatePseudoInverseSimple\n");
666 const int m = pInputMatrix->
rows;
667 const int n = pInputMatrix->
columns;
674 if (!
Invert(&ATA, &ATA_inverted))
677 MulMatMat(&ATA_inverted, &AT, pResultMatrix);
684 if (nTargetDimension > pData->
columns)
686 printf(
"error: target dimension is greater than number of columns in training data matrix in LinearAlgebra::PCA\n");
690 const int samples = pData->
rows;
691 const int dimension = pData->
columns;
693 if (pTransformationMatrix->
columns != dimension || pTransformationMatrix->
rows != nTargetDimension ||
694 pTransformedData->
columns != samples || pTransformedData->
rows != nTargetDimension)
696 printf(
"error: input to LinearAlgebra::PCA does not match\n");
712 SVD(&covarianceMatrix, &eigenValues, &eigenVectors, 0,
true,
true);
716 for (
int i = 0,
offset = 0; i < nTargetDimension; i++)
718 for (
int j = 0; j < dimension; j++)
720 pTransformationMatrix->
data[
offset] = eigenVectors.
data[i * dimension + j];
725 MulMatMat(pTransformationMatrix, pData, pTransformedData,
true);
730 const int dimension = pData->
columns;
732 if (pTransformationMatrix->
columns != dimension || pTransformationMatrix->
rows != dimension || pEigenValues->
columns != 1 || pEigenValues->
rows != dimension)
734 printf(
"error: input to LinearAlgebra::PCA does not match\n");
750 SVD(&covarianceMatrix, &eigenValues, pTransformationMatrix, 0,
true,
true);
754 for (
int i = 0; i < dimension; i++)
755 pEigenValues->
data[i] = eigenValues.
data[i * dimension + i];
760 if (pInputMatrix->
data == pResultMatrix->
data)
762 printf(
"error: pInputMatrix and pResultMatrix may not point to the same data in LinearAlgebra::Invert\n");
768 printf(
"error: input matrix must be square matrix for LinearAlgebra::Invert\n");
774 printf(
"error: input and output matrix do not match LinearAlgebra::Invert\n");
778 const int n = pInputMatrix->
columns;
786 for (i = 0; i <
n; i++)
787 tempMatrix[i][i] = 1;
789 int *pPivotRows =
new int[
n];
790 for (i = 0; i <
n; i++)
794 for (i = 0; i <
n; i++)
796 int j, nPivotColumn = 0;
798 float *helper1 = copiedMatrix[i];
799 float *helper2 = tempMatrix[i];
803 for (j = 0; j <
n; j++)
804 if (fabsf(helper1[j]) > max)
806 max = fabsf(helper1[j]);
810 pPivotRows[nPivotColumn] = i;
812 const float fPivotElement = copiedMatrix[i][nPivotColumn];
814 if (fabsf(fPivotElement) < 0.00001f)
817 delete [] pPivotRows;
821 const float fFactor = 1.0f / fPivotElement;
823 for (j = 0; j <
n; j++)
825 helper1[j] *= fFactor;
826 helper2[j] *= fFactor;
829 for (j = 0; j <
n; j++)
833 const float v = copiedMatrix[j][nPivotColumn];
836 helper1 = copiedMatrix[j];
837 helper2 = copiedMatrix[i];
838 for (k = 0; k <
n; k++)
839 helper1[k] -= v * helper2[k];
840 helper1[nPivotColumn] = 0;
842 helper1 = tempMatrix[j];
843 helper2 = tempMatrix[i];
844 for (k = 0; k <
n; k++)
845 helper1[k] -= v * helper2[k];
851 for (i = 0; i <
n; i++)
853 float *helper1 = (*pResultMatrix)[i];
854 const float *helper2 = tempMatrix[pPivotRows[i]];
856 for (
int j = 0; j <
n; j++)
857 helper1[j] = helper2[j];
860 delete [] pPivotRows;
867 memset(pMatrix->
data, 0, pMatrix->
columns * pMatrix->
rows *
sizeof(
float));
872 memset(pVector->
data, 0, pVector->
dimension *
sizeof(
float));
885 (AAT && (pResultMatrix->
columns != pMatrix->
rows || pResultMatrix->
rows != pMatrix->
rows)) ||
889 printf(
"error: matrix dimensions do not match for LinearAlgebra::SelfProduct\n");
897 pMatrix = pTransposedMatrix;
900 const int columns = pMatrix->
columns;
901 const int rows = pMatrix->
rows;
904 double *result = pResultMatrix->
data;
908 for (
int i = 0, input_offset1 = 0; i < rows; i++, input_offset1 += columns)
910 for (
int j = i, input_offset2 = input_offset1; j < rows; j++, input_offset2 += columns)
912 const double *data1 = data + input_offset1;
913 const double *data2 = data + input_offset2;
916 for (
int k = 0; k < columns; k++)
917 sum += data1[k] * data2[k];
919 result[i * rows + j] = result[j * rows + i] = sum;
926 double *temp_data = temp.
data;
928 for (
int i = 0, input_offset1 = 0; i < rows; i++, input_offset1 += columns)
930 for (
int j = i, input_offset2 = input_offset1; j < rows; j++, input_offset2 += columns)
932 const double *data1 = data + input_offset1;
933 const double *data2 = data + input_offset2;
936 for (
int k = 0; k < columns; k++)
937 sum += data1[k] * data2[k];
939 temp_data[i * rows + j] = temp_data[j * rows + i] = sum;
943 memcpy(pResultMatrix->
data, temp_data, pResultMatrix->
rows * pResultMatrix->
columns *
sizeof(
double));
954 printf(
"error: matrices A, B, and pResultMatrix do not satisfy requirements for LinearAlgebra::Multiply\n");
960 printf(
"error: matrices A, B, and pResultMatrix do not satisfy requirements for LinearAlgebra::Multiply\n");
964 const double *data1 = A->
data;
965 const double *data2 = B->
data;
966 double *result = pResultMatrix->
data;
968 const int nSize = pResultMatrix->
columns * pResultMatrix->
rows;
972 const int rowsA = A->
rows;
973 const int rowsB = B->
rows;
974 const int columns = A->
columns;
976 if (data1 != result && data2 != result)
978 for (
int i = 0, input_offset1 = 0, output_offset = 0; i < rowsA; i++, input_offset1 += columns)
980 for (
int j = 0, input_offset2 = input_offset1; j < rowsB; j++, input_offset2 += columns, output_offset++)
982 const double *data1_ = data1 + input_offset1;
983 const double *data2_ = data2 + input_offset2;
986 for (
int k = 0; k < columns; k++)
987 sum += data1_[k] * data2_[k];
989 result[output_offset] = sum;
996 double *temp_data = temp.
data;
1000 for (l = 0; l < nSize; l++)
1003 for (
int i = 0, input_offset1 = 0, output_offset = 0; i < rowsA; i++, input_offset1 += columns)
1005 for (
int j = 0, input_offset2 = input_offset1; j < rowsB; j++, input_offset2 += columns, output_offset++)
1007 const double *data1_ = data1 + input_offset1;
1008 const double *data2_ = data2 + input_offset2;
1011 for (
int k = 0; k < columns; k++)
1012 sum += data1_[k] * data2_[k];
1014 temp_data[output_offset] = sum;
1018 for (l = 0; l < nSize; l++)
1019 result[l] = temp_data[l];
1024 const int i_max = A->
rows;
1029 if (data1 != result && data2 != result)
1031 for (
int l = 0; l < nSize; l++)
1034 for (
int i = 0, input1_offset = 0; i < i_max; i++)
1036 const int result_offset_start = i * j_max;
1038 for (
int k = 0, input2_offset = 0; k < k_max; k++, input1_offset++)
1039 for (
int j = 0, result_offset = result_offset_start; j < j_max; j++, result_offset++, input2_offset++)
1040 result[result_offset] += data1[input1_offset] * data2[input2_offset];
1046 double *temp_data = temp.
data;
1050 for (l = 0; l < nSize; l++)
1053 for (
int i = 0, input1_offset = 0; i < i_max; i++)
1055 const int result_offset_start = i * j_max;
1057 for (
int k = 0, input2_offset = 0; k < k_max; k++, input1_offset++)
1058 for (
int j = 0, result_offset = result_offset_start; j < j_max; j++, result_offset++, input2_offset++)
1059 temp_data[result_offset] += data1[input1_offset] * data2[input2_offset];
1062 for (l = 0; l < nSize; l++)
1063 result[l] = temp_data[l];
1072 printf(
"error: input and output matrix do not match LinearAlgebra::Transpose\n");
1076 const int rows = pMatrix->
rows;
1077 const int columns = pMatrix->
columns;
1079 double *output = pResultMatrix->
data;
1081 if (input != output)
1083 for (
int i = 0, input_offset = 0; i < rows; i++)
1084 for (
int j = 0, output_offset = i; j < columns; j++, input_offset++, output_offset += rows)
1085 output[output_offset] = input[input_offset];
1091 double *temp_data = temp.
data;
1093 int i, input_offset;
1095 for (i = 0, input_offset = 0; i < rows; i++)
1096 for (
int j = 0, output_offset = i; j < columns; j++, input_offset++, output_offset += rows)
1097 temp_data[output_offset] = input[input_offset];
1099 const int nSize = rows * columns;
1101 for (i = 0; i < nSize; i++)
1102 output[i] = temp_data[i];
1110 printf(
"error: input matrix and vector and output vector do not match for CDoubleMatrix::Multiply\n");
1114 const int rows = pMatrix->
rows;
1115 const int columns = pMatrix->
columns;
1116 const double *data_m = pMatrix->
data;
1117 const double *data_v = pVector->
data;
1118 double *data_r = pResultVector->
data;
1120 if (data_r != data_v)
1122 for (
int i = 0,
offset = 0; i < rows; i++)
1126 for (
int j = 0; j < columns; j++,
offset++)
1127 sum += data_m[
offset] * data_v[j];
1134 double *temp =
new double[rows];
1137 for (i = 0, offset = 0; i < rows; i++)
1141 for (
int j = 0; j < columns; j++, offset++)
1142 sum += data_m[offset] * data_v[j];
1147 for (i = 0; i < rows; i++)
1148 data_r[i] = temp[i];
1156 MulMatVec(pMatrix, pVector1, pResultVector);
1165 printf(
"error: matrix dimensions do not match in LinearAlgebra::AddMatMat\n");
1169 const int nSize = pMatrix1->
rows * pMatrix1->
columns;
1170 const double *data1 = pMatrix1->
data;
1171 const double *data2 = pMatrix2->
data;
1172 double *result = pResultMatrix->
data;
1174 for (
int i = 0; i < nSize; i++)
1175 result[i] = data1[i] + data2[i];
1183 printf(
"error: matrix dimensions do not match in LinearAlgebra::SubtractMatMat\n");
1187 const int nSize = pMatrix1->
rows * pMatrix1->
columns;
1188 const double *data1 = pMatrix1->
data;
1189 const double *data2 = pMatrix2->
data;
1190 double *result = pResultMatrix->
data;
1192 for (
int i = 0; i < nSize; i++)
1193 result[i] = data1[i] - data2[i];
1200 printf(
"error: vector dimensions do not match in LinearAlgebra::AddVecVec\n");
1204 const int nDimension = pVector1->
dimension;
1205 const double *data1 = pVector1->
data;
1206 const double *data2 = pVector2->
data;
1207 double *result = pResultVector->
data;
1209 for (
int i = 0; i < nDimension; i++)
1210 result[i] = data1[i] + data2[i];
1217 printf(
"error: vector dimensions do not match in LinearAlgebra::SubtractVecVec\n");
1221 const int nDimension = pVector1->
dimension;
1222 const double *data1 = pVector1->
data;
1223 const double *data2 = pVector2->
data;
1224 double *result = pResultVector->
data;
1226 for (
int i = 0; i < nDimension; i++)
1227 result[i] = data1[i] - data2[i];
1234 printf(
"error: matrix dimensions do not match in LinearAlgebra::AddToMat\n");
1238 const int nSize = pMatrix->
rows * pMatrix->
columns;
1240 const double *data_to_add = pMatrixToAdd->
data;
1242 for (
int i = 0; i < nSize; i++)
1243 data[i] += data_to_add[i];
1250 printf(
"error: matrix dimensions do not match in LinearAlgebra::SubtractFromMat\n");
1254 const int nSize = pMatrix->
rows * pMatrix->
columns;
1256 const double *data_to_subtract = pMatrixToAdd->
data;
1258 for (
int i = 0; i < nSize; i++)
1259 data[i] -= data_to_subtract[i];
1266 printf(
"error: vector dimensions do not match in LinearAlgebra::AddToVec\n");
1270 const int nDimension = pVector->
dimension;
1272 const double *data_to_add = pVectorToAdd->
data;
1274 for (
int i = 0; i < nDimension; i++)
1275 data[i] += data_to_add[i];
1282 printf(
"error: vector dimensions do not match in LinearAlgebra::SubtractFromVec\n");
1286 const int nDimension = pVector->
dimension;
1288 const double *data_to_subtract = pVectorToSubtract->
data;
1290 for (
int i = 0; i < nDimension; i++)
1291 data[i] -= data_to_subtract[i];
1300 const int columns = pMatrix->
columns;
1301 const int rows = pMatrix->
rows;
1302 const int max = rows * columns;
1303 const double *
data = pMatrix->
data;
1304 double *output = pResultMatrix->
data;
1306 for (
int i = 0; i < columns; i++)
1311 for (j = 0; j < max; j += columns)
1316 for (j = 0; j < max; j += columns)
1317 output[j] = data[j] - mean;
1329 const int columns = pMatrix->
columns;
1330 const int rows = pMatrix->
rows;
1331 const double *
data = pMatrix->
data;
1332 double *output = pResultMatrix->
data;
1334 for (
int i = 0; i < rows; i++)
1339 for (j = 0; j < columns; j++)
1344 for (j = 0; j < columns; j++)
1345 output[j] = data[j] - mean;
1354 if (pInputMatrix->
data == pResultMatrix->
data)
1356 printf(
"error: pInputMatrix and pResultMatrix may not point to the same data in LinearAlgebra::Invert\n");
1362 printf(
"error: input matrix must be square matrix for LinearAlgebra::Invert\n");
1368 printf(
"error: input and output matrix do not match LinearAlgebra::Invert\n");
1372 const int n = pInputMatrix->
columns;
1380 for (i = 0; i <
n; i++)
1381 tempMatrix[i][i] = 1;
1383 int *pPivotRows =
new int[
n];
1384 for (i = 0; i <
n; i++)
1388 for (i = 0; i <
n; i++)
1390 int j, nPivotColumn = 0;
1392 double *helper1 = copiedMatrix[i];
1393 double *helper2 = tempMatrix[i];
1397 for (j = 0; j <
n; j++)
1398 if (fabs(helper1[j]) > max)
1400 max = fabs(helper1[j]);
1404 pPivotRows[nPivotColumn] = i;
1406 const double dPivotElement = copiedMatrix[i][nPivotColumn];
1408 if (fabs(dPivotElement) < 0.00001)
1411 delete [] pPivotRows;
1415 const double dFactor = 1.0 / dPivotElement;
1417 for (j = 0; j <
n; j++)
1419 helper1[j] *= dFactor;
1420 helper2[j] *= dFactor;
1423 for (j = 0; j <
n; j++)
1427 const double v = copiedMatrix[j][nPivotColumn];
1430 helper1 = copiedMatrix[j];
1431 helper2 = copiedMatrix[i];
1432 for (k = 0; k <
n; k++)
1433 helper1[k] -= v * helper2[k];
1434 helper1[nPivotColumn] = 0;
1436 helper1 = tempMatrix[j];
1437 helper2 = tempMatrix[i];
1438 for (k = 0; k <
n; k++)
1439 helper1[k] -= v * helper2[k];
1445 for (i = 0; i <
n; i++)
1447 double *helper1 = (*pResultMatrix)[i];
1448 const double *helper2 = tempMatrix[pPivotRows[i]];
1450 for (
int j = 0; j <
n; j++)
1451 helper1[j] = helper2[j];
1454 delete [] pPivotRows;
1463 printf(
"error: A, x do not match LinearAlgebra::SolveLinearLeastSquaresHomogeneousSVD\n");
1469 printf(
"error: equation system is underdetermined in LinearAlgebra::SolveLinearLeastSquaresHomogeneousSVD\n");
1473 const int m = A->
rows;
1478 SVD(A, &W, 0, &V,
false,
false,
false);
1488 printf(
"error: A, b, x do not match LinearAlgebra::SolveLinearLeastSquaresSVD\n");
1494 printf(
"error: equation system is underdetermined in LinearAlgebra::SolveLinearLeastSquaresSVD\n");
1507 printf(
"error: A, b, x do not match LinearAlgebra::SolveLinearLeastSquaresSimple\n");
1513 printf(
"error: equation system is underdetermined in LinearAlgebra::SolveLinearLeastSquaresSimple\n");
1531 printf(
"error: input and output matrix do not match LinearAlgebra::CalculatePseudoInverseSVD\n");
1541 const int m = pInputMatrix->
rows;
1542 const int n = pInputMatrix->
columns;
1547 SVD(A, &W, &UT, &V,
false,
true,
false);
1552 const int min =
MY_MIN(m, n);
1555 double dThreshold = 0.0;
1557 for(i = 0; i < min; i++)
1558 dThreshold += WT(i, i);
1560 dThreshold *= 2 * DBL_EPSILON;
1563 for (i = 0; i < min; i++)
1564 if (WT(i, i) < dThreshold)
1567 WT(i, i) = 1.0 / WT(i, i);
1579 printf(
"error: input and output matrix do not match LinearAlgebra::CalculatePseudoInverseSimple\n");
1587 const int m = pInputMatrix->
rows;
1588 const int n = pInputMatrix->
columns;
1595 if (!
Invert(&ATA, &ATA_inverted))
1598 MulMatMat(&ATA_inverted, &AT, pResultMatrix);
1605 memset(pVector->
data, 0, pVector->
dimension *
sizeof(
double));
1610 memset(pMatrix->
data, 0, pMatrix->
columns * pMatrix->
rows *
sizeof(
double));
1620 printf(
"error: not enough input point pairs for LinearAlgebra::DetermineAffineTransformation (must be at least 3)\n");
1629 for (
int i = 0,
offset = 0; i < nPoints; i++,
offset += 12)
1631 data[
offset] = pSourcePoints[i].
x;
1632 data[offset + 1] = pSourcePoints[i].
y;
1633 data[offset + 2] = 1;
1634 data[offset + 3] = 0;
1635 data[offset + 4] = 0;
1636 data[offset + 5] = 0;
1638 data[offset + 6] = 0;
1639 data[offset + 7] = 0;
1640 data[offset + 8] = 0;
1641 data[offset + 9] = pSourcePoints[i].
x;
1642 data[offset + 10] = pSourcePoints[i].
y;
1643 data[offset + 11] = 1;
1645 const int index = 2 * i;
1647 b.
data[index + 1] = pTargetPoints[i].
y;
1652 for (
int i = 0; i < 6; i++)
1674 printf(
"error: not enough input point pairs for LinearAlgebra::DetermineHomography (must be at least 4)\n");
1685 for (
int i = 0,
offset = 0; i < nPoints; i++,
offset += 16)
1687 data[
offset] = pSourcePoints[i].
x;
1688 data[offset + 1] = pSourcePoints[i].
y;
1689 data[offset + 2] = 1;
1690 data[offset + 3] = 0;
1691 data[offset + 4] = 0;
1692 data[offset + 5] = 0;
1693 data[offset + 6] = -pSourcePoints[i].
x * pTargetPoints[i].
x;
1694 data[offset + 7] = -pSourcePoints[i].
y * pTargetPoints[i].
x;
1696 data[offset + 8] = 0;
1697 data[offset + 9] = 0;
1698 data[offset + 10] = 0;
1699 data[offset + 11] = pSourcePoints[i].
x;
1700 data[offset + 12] = pSourcePoints[i].
y;
1701 data[offset + 13] = 1;
1702 data[offset + 14] = -pSourcePoints[i].
x * pTargetPoints[i].
y;
1703 data[offset + 15] = -pSourcePoints[i].
y * pTargetPoints[i].
y;
1705 const int index = 2 * i;
1707 b.
data[index + 1] = pTargetPoints[i].
y;
1712 for (
int i = 0; i < 8; i++)
void SubtractFromMat(CFloatMatrix *pMatrix, const CFloatMatrix *pMatrixToSubtract)
bool DetermineHomography(const Vec2d *pSourcePoints, const Vec2d *pTargetPoints, int nPoints, Mat3d &A, bool bUseSVD=false)
Determines a homography based on a set of 2d-2d point correspondences.
void AddToMat(CFloatMatrix *pMatrix, const CFloatMatrix *pMatrixToAdd)
void SolveLinearLeastSquaresHomogeneousSVD(const CFloatMatrix *A, CFloatVector *x)
bool Invert(const CByteImage *pInputImage, CByteImage *pOutputImage)
Calculates the inverted image of a CByteImage and writes the result to a CByteImage.
bool CalculatePseudoInverseSimple(const CFloatMatrix *pInputMatrix, CFloatMatrix *pResultMatrix)
void SolveLinearLeastSquaresSVD(const CFloatMatrix *A, const CFloatVector *b, CFloatVector *x)
Data structure for the representation of a vector of values of the data type double.
void PCA(const CFloatMatrix *pData, CFloatMatrix *pTransformationMatrix, CFloatMatrix *pTransformedData, int nTargetDimension)
void SubtractFromVec(CFloatVector *pVector, const CFloatVector *pVectorToSubtract)
void AddVecVec(const CFloatVector *pVector1, const CFloatVector *pVector2, CFloatVector *pResultVector)
void SVD(const CFloatMatrix *A, CFloatMatrix *W, CFloatMatrix *U=0, CFloatMatrix *V=0, bool bAllowModifyA=false, bool bReturnUTransposed=false, bool bReturnVTransposed=false)
bool DetermineAffineTransformation(const Vec2d *pSourcePoints, const Vec2d *pTargetPoints, int nPoints, Mat3d &A, bool bUseSVD=false)
Determines an affine transformation based on a set of 2d-2d point correspondences.
bool CopyMatrix(const CFloatMatrix *pInputImage, CFloatMatrix *pOutputImage)
Copies one CFloatMatrix to another.
Data structure for the representation of a matrix of values of the data type float.
void SubtractMeanFromRows(const CFloatMatrix *pMatrix, CFloatMatrix *pResultMatrix)
GLsizei GLsizei GLenum GLenum const GLvoid * data
void AddMatMat(const CFloatMatrix *pMatrix1, const CFloatMatrix *pMatrix2, CFloatMatrix *pResultMatrix)
void SubtractMeanFromColumns(const CFloatMatrix *pMatrix, CFloatMatrix *pResultMatrix)
void MulMatVec(const CFloatMatrix *pMatrix, const CFloatVector *pVector, CFloatVector *pResultVector)
bool Invert(const CFloatMatrix *pInputMatrix, CFloatMatrix *pResultMatrix)
void MulMatMat(const CFloatMatrix *A, const CFloatMatrix *B, CFloatMatrix *pResultMatrix, bool bTransposeB=false)
void Zero(CFloatMatrix *pMatrix)
void CalculatePseudoInverseSVD(const CFloatMatrix *pInputMatrix, CFloatMatrix *pOutputMatrix)
Data structure for the representation of a vector of values of the data type float.
GLenum GLenum GLenum input
bool SolveLinearLeastSquaresSimple(const CFloatMatrix *A, const CFloatVector *b, CFloatVector *x)
Data structure for the representation of a 2D vector.
void SelfProduct(const CFloatMatrix *pMatrix, CFloatMatrix *pResultMatrix, bool AAT=false)
Data structure for the representation of a matrix of values of the data type double.
void Transpose(const CFloatMatrix *pMatrix, CFloatMatrix *pResultMatrix)
void SubtractMatMat(const CFloatMatrix *pMatrix1, const CFloatMatrix *pMatrix2, CFloatMatrix *pResultMatrix)
void AddToVec(CFloatVector *pVector, const CFloatVector *pVectorToAdd)
void SubtractVecVec(const CFloatVector *pVector1, const CFloatVector *pVector2, CFloatVector *pResultVector)
void Zero(CByteImage *pImage, const MyRegion *pROI=0)
Sets all values in a CByteImage to zero.
void SetMat(Mat3d &matrix, float r1, float r2, float r3, float r4, float r5, float r6, float r7, float r8, float r9)
Data structure for the representation of a 3x3 matrix.