Commit a01fc64a authored by Sean McBride's avatar Sean McBride Committed by Kitware Robot

Merge topic 'vtkMathReview'

b889b631 Introduce VTK's first use of C++11's std::lround() (in a test)
27579ded Deprecated non-thread-safe colour conversion functions
f98c02bb Code review of vtkMath clamping
b4f84fe1 Added const to various vtkMath API
a7ceef61 Replace some int with vtkTypeBool (if public) or bool (if private)
a440133b General code review of vtkMath: trivial & conservative
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: David Gobbi's avatarDavid Gobbi <david.gobbi@gmail.com>
Merge-request: !3299
parents 6144e9eb b889b631
......@@ -607,7 +607,6 @@ static int TestColorConvert(const Triple &rgb, const Triple &hsv,
cout << " CIE-L*ab: " << lab << endl;
Triple result1;
double *result2;
#define COMPARE(testname, target, dest) \
if ((target) != (dest)) \
......@@ -623,11 +622,6 @@ static int TestColorConvert(const Triple &rgb, const Triple &hsv,
vtkMath::HSVToRGB(hsv(), result1());
COMPARE(HSVToRGB, rgb, result1);
result2 = vtkMath::RGBToHSV(rgb());
COMPARE(RGBToHSV, hsv, result2);
result2 = vtkMath::HSVToRGB(hsv());
COMPARE(HSVToRGB, rgb, result2);
vtkMath::RGBToHSV(rgb[0], rgb[1], rgb[2],
&result1[0], &result1[1], &result1[2]);
COMPARE(RGBToHSV, hsv, result1);
......@@ -641,11 +635,6 @@ static int TestColorConvert(const Triple &rgb, const Triple &hsv,
vtkMath::XYZToRGB(xyz(), result1());
COMPARE(XYZToRGB, rgb, result1);
result2 = vtkMath::RGBToXYZ(rgb());
COMPARE(RGBToXYZ, xyz, result2);
result2 = vtkMath::XYZToRGB(xyz());
COMPARE(XYZToRGB, rgb, result2);
vtkMath::RGBToXYZ(rgb[0], rgb[1], rgb[2],
&result1[0], &result1[1], &result1[2]);
COMPARE(RGBToXYZ, xyz, result1);
......@@ -659,11 +648,6 @@ static int TestColorConvert(const Triple &rgb, const Triple &hsv,
vtkMath::XYZToLab(xyz(), result1());
COMPARE(XYZToLab, lab, result1);
result2 = vtkMath::LabToXYZ(lab());
COMPARE(LabToXYZ, xyz, result2);
result2 = vtkMath::XYZToLab(xyz());
COMPARE(XYZToLab, lab, result2);
vtkMath::LabToXYZ(lab[0], lab[1], lab[2],
&result1[0], &result1[1], &result1[2]);
COMPARE(LabToXYZ, xyz, result1);
......@@ -677,11 +661,6 @@ static int TestColorConvert(const Triple &rgb, const Triple &hsv,
vtkMath::RGBToLab(rgb(), result1());
COMPARE(RGBToLab, lab, result1);
result2 = vtkMath::LabToRGB(lab());
COMPARE(LabToRGB, rgb, result2);
result2 = vtkMath::RGBToLab(rgb());
COMPARE(RGBToLab, lab, result2);
vtkMath::LabToRGB(lab[0], lab[1], lab[2],
&result1[0], &result1[1], &result1[2]);
COMPARE(LabToRGB, rgb, result1);
......
......@@ -23,10 +23,6 @@
#include <vector>
static int TestPi();
#if 0
static int TestDoublePi();
static int TestDoubleTwoPi();
#endif
static int TestDegreesFromRadians();
static int TestRound();
static int TestFloor();
......@@ -98,11 +94,6 @@ int UnitTestMath(int,char *[])
status += TestPi();
#if 0
status += TestDoublePi(); // legacy
status += TestDoubleTwoPi(); // legacy
#endif
status += TestDegreesFromRadians();
status += TestRound();
status += TestFloor();
......@@ -202,56 +193,6 @@ int TestPi()
return status;
}
#if 0
// Validate by comparing to atan/4
int TestDoublePi()
{
int status = 0;
std::cout << "DoublePi..";
if (vtkMath::DoublePi() != std::atan(1.0) * 4.0)
{
std::cout << "Expected " << vtkMath::Pi()
<< " but got " << std::atan(1.0) * 4.0;
++status;
}
if (status)
{
std::cout << "..FAILED" << std::endl;
}
else
{
std::cout << ".PASSED" << std::endl;
}
return status;
}
// Validate by comparing to atan/4 * 2
int TestDoubleTwoPi()
{
int status = 0;
std::cout << "DoubleTwoPi..";
if (vtkMath::DoubleTwoPi() != std::atan(1.0) * 4.0 * 2.0)
{
std::cout << "Expected " << vtkMath::Pi() * 2.0
<< " but got " << std::atan(1.0) * 4.0 * 2.0;
++status;
}
if (status)
{
std::cout << "..FAILED" << std::endl;
}
else
{
std::cout << ".PASSED" << std::endl;
}
return status;
}
#endif
// Validate against RadiansFromDegress
int TestDegreesFromRadians()
{
......
......@@ -67,7 +67,7 @@ vtkMathInternal::vtkMathInternal()
// an initial vtkMinimalStandardRandomSequence is created.
this->Uniform=static_cast<vtkMinimalStandardRandomSequence *>(
this->Gaussian->GetUniformSequence());
this->Uniform->SetSeedOnly(1177); // One authors home address
this->Uniform->SetSeedOnly(1177); // One author's home address
this->MemoizeFactorial.resize(21, 0);
}
......@@ -138,7 +138,7 @@ int vtkMath::CeilLog2(vtkTypeUInt64 x)
int y = (((x & (x - 1)) == 0) ? 0 : 1);
// loop through the table (this unrolls nicely)
for (int i = 0; i < 6; i++)
for (int i = 0; i < 6; ++i)
{
int k = (((x & t[i]) == 0) ? 0 : j);
y += k;
......@@ -165,7 +165,6 @@ double vtkMath::Random()
// is proportional to the seed value! To help solve this, call
// RandomSeed() a few times inside seed. This doesn't ruin the
// repeatability of Random().
//
void vtkMath::RandomSeed(int s)
{
vtkMath::Internal->Uniform->SetSeed(s);
......@@ -226,7 +225,6 @@ vtkTypeInt64 vtkMath::Factorial( int N )
//----------------------------------------------------------------------------
// The number of combinations of n objects from a pool of m objects (m>n).
//
vtkTypeInt64 vtkMath::Binomial( int m, int n )
{
double r = 1;
......@@ -241,7 +239,6 @@ vtkTypeInt64 vtkMath::Binomial( int m, int n )
// Start iterating over "m choose n" objects.
// This function returns an array of n integers, each from 0 to m-1.
// These integers represent the n items chosen from the set [0,m[.
//
int* vtkMath::BeginCombination( int m, int n )
{
if ( m < n )
......@@ -265,7 +262,6 @@ int* vtkMath::BeginCombination( int m, int n )
// If the \a combination is the last item in the sequence on input,
// then \a combination is unaltered and 0 is returned.
// Otherwise, 1 is returned and \a combination is updated.
//
int vtkMath::NextCombination( int m, int n, int* r )
{
int status = 0;
......@@ -300,14 +296,13 @@ template<class T1, class T2, class T3>
inline void vtkMathPerpendiculars(const T1 v1[3], T2 v2[3], T3 v3[3],
double theta)
{
int dv1,dv2,dv3;
double v1sq = v1[0]*v1[0];
double v2sq = v1[1]*v1[1];
double v3sq = v1[2]*v1[2];
double r = sqrt(v1sq + v2sq + v3sq);
// transpose the vector to avoid divide-by-zero error
int dv1, dv2, dv3;
if (v1sq > v2sq && v1sq > v3sq)
{
dv1 = 0; dv2 = 1; dv3 = 2;
......@@ -327,7 +322,7 @@ inline void vtkMathPerpendiculars(const T1 v1[3], T2 v2[3], T3 v3[3],
double tmp = sqrt(a*a+c*c);
if (theta != 0)
if (theta != 0.0)
{
double sintheta = sin(theta);
double costheta = cos(theta);
......@@ -382,17 +377,15 @@ void vtkMath::Perpendiculars(const float v1[3], float v2[3], float v3[3],
// Solve linear equations Ax = b using Crout's method. Input is square matrix A
// and load vector x. Solution x is written over load vector. The dimension of
// the matrix is specified in size. If error is found, method returns a 0.
int vtkMath::SolveLinearSystem(double **A, double *x, int size)
vtkTypeBool vtkMath::SolveLinearSystem(double **A, double *x, int size)
{
// if we solving something simple, just solve it
//
if (size == 2)
{
double det, y[2];
double det = vtkMath::Determinant2x2(A[0][0], A[0][1], A[1][0], A[1][1]);
det = vtkMath::Determinant2x2(A[0][0], A[0][1], A[1][0], A[1][1]);
static const double eps = 256*std::numeric_limits<double>::epsilon();
static const double eps = 256.0 * std::numeric_limits<double>::epsilon();
if (std::fabs(det) < eps)
{
......@@ -400,6 +393,7 @@ int vtkMath::SolveLinearSystem(double **A, double *x, int size)
return 0;
}
double y[2];
y[0] = (A[1][1]*x[0] - A[0][1]*x[1]) / det;
y[1] = (-A[1][0]*x[0] + A[0][0]*x[1]) / det;
......@@ -435,9 +429,12 @@ int vtkMath::SolveLinearSystem(double **A, double *x, int size)
{
return 0;
}
vtkMath::LUSolveLinearSystem(A,index,x,size);
vtkMath::LUSolveLinearSystem(A, index, x, size);
if (size >= 10 ) delete [] index;
if (size >= 10 )
{
delete [] index;
}
return 1;
}
......@@ -445,7 +442,7 @@ int vtkMath::SolveLinearSystem(double **A, double *x, int size)
// Invert input square matrix A into matrix AI. Note that A is modified during
// the inversion. The size variable is the dimension of the matrix. Returns 0
// if inverse not computed.
int vtkMath::InvertMatrix(double **A, double **AI, int size)
vtkTypeBool vtkMath::InvertMatrix(double **A, double **AI, int size)
{
int *index, iScratch[10];
double *column, dScratch[10];
......@@ -463,7 +460,7 @@ int vtkMath::InvertMatrix(double **A, double **AI, int size)
column = new double[size];
}
int retVal = vtkMath::InvertMatrix(A, AI, size, index, column);
vtkTypeBool retVal = vtkMath::InvertMatrix(A, AI, size, index, column);
if ( size > 10 )
{
......@@ -480,7 +477,7 @@ int vtkMath::InvertMatrix(double **A, double **AI, int size)
// square matrix A, integer array of pivot indices index[0->n-1], and size
// of square matrix n. Output factorization LU is in matrix A. If error is
// found, method returns 0.
int vtkMath::LUFactorLinearSystem(double **A, int *index, int size)
vtkTypeBool vtkMath::LUFactorLinearSystem(double **A, int *index, int size)
{
double scratch[10];
double *scale = (size<10 ? scratch : new double[size]);
......@@ -492,9 +489,9 @@ int vtkMath::LUFactorLinearSystem(double **A, int *index, int size)
//
// Loop over rows to get implicit scaling information
//
for ( i = 0; i < size; i++ )
for ( i = 0; i < size; ++i )
{
for ( largest = 0.0, j = 0; j < size; j++ )
for ( largest = 0.0, j = 0; j < size; ++j )
{
if ( (temp2 = fabs(A[i][j])) > largest )
{
......@@ -507,17 +504,17 @@ int vtkMath::LUFactorLinearSystem(double **A, int *index, int size)
vtkGenericWarningMacro(<<"Unable to factor linear system");
return 0;
}
scale[i] = 1.0 / largest;
scale[i] = 1.0 / largest;
}
//
// Loop over all columns using Crout's method
//
for ( j = 0; j < size; j++ )
for ( j = 0; j < size; ++j )
{
for (i = 0; i < j; i++)
for (i = 0; i < j; ++i)
{
sum = A[i][j];
for ( k = 0; k < i; k++ )
for ( k = 0; k < i; ++k )
{
sum -= A[i][k] * A[k][j];
}
......@@ -526,10 +523,10 @@ int vtkMath::LUFactorLinearSystem(double **A, int *index, int size)
//
// Begin search for largest pivot element
//
for ( largest = 0.0, i = j; i < size; i++ )
for ( largest = 0.0, i = j; i < size; ++i )
{
sum = A[i][j];
for ( k = 0; k < j; k++ )
for ( k = 0; k < j; ++k )
{
sum -= A[i][k] * A[k][j];
}
......@@ -546,7 +543,7 @@ int vtkMath::LUFactorLinearSystem(double **A, int *index, int size)
//
if ( j != maxI )
{
for ( k = 0; k < size; k++ )
for ( k = 0; k < size; ++k )
{
temp1 = A[maxI][k];
A[maxI][k] = A[j][k];
......@@ -568,14 +565,17 @@ int vtkMath::LUFactorLinearSystem(double **A, int *index, int size)
if ( j != (size-1) )
{
temp1 = 1.0 / A[j][j];
for ( i = j + 1; i < size; i++ )
for ( i = j + 1; i < size; ++i )
{
A[i][j] *= temp1;
}
}
}
if (size >= 10 ) delete [] scale;
if (size >= 10 )
{
delete [] scale;
}
return 1;
}
......@@ -596,7 +596,7 @@ void vtkMath::LUSolveLinearSystem(double **A, int *index,
// Proceed with forward and backsubstitution for L and U
// matrices. First, forward substitution.
//
for ( ii = -1, i = 0; i < size; i++ )
for ( ii = -1, i = 0; i < size; ++i )
{
idx = index[i];
sum = x[idx];
......@@ -604,7 +604,7 @@ void vtkMath::LUSolveLinearSystem(double **A, int *index,
if ( ii >= 0 )
{
for ( j = ii; j <= (i-1); j++ )
for ( j = ii; j <= (i-1); ++j )
{
sum -= A[i][j]*x[j];
}
......@@ -622,7 +622,7 @@ void vtkMath::LUSolveLinearSystem(double **A, int *index,
for ( i = size-1; i >= 0; i-- )
{
sum = x[i];
for ( j = i + 1; j < size; j++ )
for ( j = i + 1; j < size; ++j )
{
sum -= A[i][j]*x[j];
}
......@@ -648,7 +648,7 @@ void vtkMath::LUSolveLinearSystem(double **A, int *index,
// normalized.
// It assumes a is symmetric and uses only its upper right triangular part.
template<class T>
int vtkJacobiN(T **a, int n, T *w, T **v)
vtkTypeBool vtkJacobiN(T **a, int n, T *w, T **v)
{
int i, j, k, iq, ip, numPos;
T tresh, theta, tau, t, sm, s, h, g, c, tmp;
......@@ -679,7 +679,7 @@ int vtkJacobiN(T **a, int n, T *w, T **v)
}
// begin rotation sequence
for (i=0; i<VTK_MAX_ROTATIONS; i++)
for (i=0; i<VTK_MAX_ROTATIONS; ++i)
{
sm = 0.0;
for (ip=0; ip<n-1; ip++)
......@@ -694,7 +694,7 @@ int vtkJacobiN(T **a, int n, T *w, T **v)
break;
}
if (i < 3) // first 3 sweeps
if (i < 3) // first 3 sweeps
{
tresh = 0.2*sm/(n*n);
}
......@@ -742,21 +742,21 @@ int vtkJacobiN(T **a, int n, T *w, T **v)
a[ip][iq]=0.0;
// ip already shifted left by 1 unit
for (j = 0;j <= ip-1;j++)
for (j = 0;j <= ip-1;++j)
{
VTK_ROTATE(a,j,ip,j,iq);
}
// ip and iq already shifted left by 1 unit
for (j = ip+1;j <= iq-1;j++)
for (j = ip+1;j <= iq-1;++j)
{
VTK_ROTATE(a,ip,j,j,iq);
}
// iq already shifted left by 1 unit
for (j=iq+1; j<n; j++)
for (j=iq+1; j<n; ++j)
{
VTK_ROTATE(a,ip,j,iq,j);
}
for (j=0; j<n; j++)
for (j=0; j<n; ++j)
{
VTK_ROTATE(v,j,ip,j,iq);
}
......@@ -781,11 +781,11 @@ int vtkJacobiN(T **a, int n, T *w, T **v)
}
// sort eigenfunctions these changes do not affect accuracy
for (j=0; j<n-1; j++) // boundary incorrect
for (j=0; j<n-1; ++j) // boundary incorrect
{
k = j;
tmp = w[k];
for (i=j+1; i<n; i++) // boundary incorrect, shifted already
for (i=j+1; i<n; ++i) // boundary incorrect, shifted already
{
if (w[i] >= tmp) // why exchange if same?
{
......@@ -797,7 +797,7 @@ int vtkJacobiN(T **a, int n, T *w, T **v)
{
w[k] = w[j];
w[j] = tmp;
for (i=0; i<n; i++)
for (i=0; i<n; ++i)
{
tmp = v[i][j];
v[i][j] = v[i][k];
......@@ -810,19 +810,18 @@ int vtkJacobiN(T **a, int n, T *w, T **v)
// reek havoc in hyperstreamline/other stuff. We will select the most
// positive eigenvector.
int ceil_half_n = (n >> 1) + (n & 1);
for (j=0; j<n; j++)
for (j=0; j<n; ++j)
{
for (numPos=0, i=0; i<n; i++)
for (numPos=0, i=0; i<n; ++i)
{
if ( v[i][j] >= 0.0 )
{
numPos++;
}
}
// if ( numPos < ceil(double(n)/double(2.0)) )
if ( numPos < ceil_half_n)
{
for(i=0; i<n; i++)
for(i=0; i<n; ++i)
{
v[i][j] *= -1.0;
}
......@@ -841,13 +840,13 @@ int vtkJacobiN(T **a, int n, T *w, T **v)
#undef VTK_MAX_ROTATIONS
//----------------------------------------------------------------------------
int vtkMath::JacobiN(float **a, int n, float *w, float **v)
vtkTypeBool vtkMath::JacobiN(float **a, int n, float *w, float **v)
{
return vtkJacobiN(a,n,w,v);
}
//----------------------------------------------------------------------------
int vtkMath::JacobiN(double **a, int n, double *w, double **v)
vtkTypeBool vtkMath::JacobiN(double **a, int n, double *w, double **v)
{
return vtkJacobiN(a,n,w,v);
}
......@@ -858,13 +857,13 @@ int vtkMath::JacobiN(double **a, int n, double *w, double **v)
// real symmetric matrix. Square 3x3 matrix a; output eigenvalues in w;
// and output eigenvectors in v. Resulting eigenvalues/vectors are sorted
// in decreasing order; eigenvectors are normalized.
int vtkMath::Jacobi(float **a, float *w, float **v)
vtkTypeBool vtkMath::Jacobi(float **a, float *w, float **v)
{
return vtkMath::JacobiN(a, 3, w, v);
}
//----------------------------------------------------------------------------
int vtkMath::Jacobi(double **a, double *w, double **v)
vtkTypeBool vtkMath::Jacobi(double **a, double *w, double **v)
{
return vtkMath::JacobiN(a, 3, w, v);
}
......@@ -883,9 +882,9 @@ double vtkMath::EstimateMatrixCondition(const double *const *A, int size)
double min=VTK_FLOAT_MAX, max=(-VTK_FLOAT_MAX);
// find the maximum value
for (i=0; i < size; i++)
for (i=0; i < size; ++i)
{
for (j=i; j < size; j++)
for (j=i; j < size; ++j)
{
if ( fabs(A[i][j]) > max )
{
......@@ -895,7 +894,7 @@ double vtkMath::EstimateMatrixCondition(const double *const *A, int size)
}
// find the minimum diagonal value
for (i=0; i < size; i++)
for (i=0; i < size; ++i)
{
if ( fabs(A[i][i]) < min )
{
......@@ -924,8 +923,9 @@ double vtkMath::EstimateMatrixCondition(const double *const *A, int size)
// M' dimension is xOrder by 1.
// M' should be pre-allocated. All matrices are row major. The resultant
// matrix M' should be pre-multiplied to X' to get 0', or transposed and
// then post multiplied to X to get 0
int vtkMath::SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int xOrder,
// then post multiplied to X to get 0.
// Returns success/fail.
vtkTypeBool vtkMath::SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int xOrder,
double **mt)
{
// check dimensional consistency
......@@ -935,8 +935,6 @@ int vtkMath::SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int
return 0;
}
int i, j, k;
// set up intermediate variables
// Allocate matrix to hold X times transpose of X
double **XXt = new double *[xOrder]; // size x by x
......@@ -944,24 +942,25 @@ int vtkMath::SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int
double *eigenvals = new double [xOrder];
double **eigenvecs = new double *[xOrder];
int i, j, k;
// Clear the upper triangular region (and btw, allocate the eigenvecs as well)
for (i = 0; i < xOrder; i++)
for (i = 0; i < xOrder; ++i)
{
eigenvecs[i] = new double[xOrder];
XXt[i] = new double[xOrder];
for (j = 0; j < xOrder; j++)
for (j = 0; j < xOrder; ++j)
{
XXt[i][j] = 0.0;
}
}
// Calculate XXt upper half only, due to symmetry
for (k = 0; k < numberOfSamples; k++)
for (k = 0; k < numberOfSamples; ++k)
{
for (i = 0; i < xOrder; i++)
for (i = 0; i < xOrder; ++i)
{
for (j = i; j < xOrder; j++)
for (j = i; j < xOrder; ++j)
{
XXt[i][j] += xt[k][i] * xt[k][j];
}
......@@ -969,9 +968,9 @@ int vtkMath::SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int
}
// now fill in the lower half of the XXt matrix
for (i = 0; i < xOrder; i++)
for (i = 0; i < xOrder; ++i)
{
for (j = 0; j < i; j++)
for (j = 0; j < i; ++j)
{
XXt[i][j] = XXt[j][i];
}
......@@ -982,13 +981,13 @@ int vtkMath::SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int
// Smallest eigenval is at the end of the list (xOrder-1), and solution is
// corresponding eigenvec.
for (i=0; i<xOrder; i++)
for (i=0; i<xOrder; ++i)
{
mt[i][0] = eigenvecs[i][xOrder-1];
}
// Clean up:
for (i=0; i<xOrder; i++)
for (i=0; i<xOrder; ++i)
{
delete [] XXt[i];
delete [] eigenvecs[i];
......@@ -999,7 +998,8 @@ int vtkMath::SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int
return 1;
}
#define VTK_SMALL_NUMBER 1.0e-12
static const double VTK_SMALL_NUMBER = 1.0e-12;
//----------------------------------------------------------------------------
// Solves for the least squares best fit matrix for the equation X'M' = Y'.
......@@ -1014,8 +1014,9 @@ int vtkMath::SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int
// By default, this method checks for the homogeneous condition where Y==0, and
// if so, invokes SolveHomogeneousLeastSquares. For better performance when
// the system is known not to be homogeneous, invoke with checkHomogeneous=0.
int vtkMath::SolveLeastSquares(int numberOfSamples, double **xt, int xOrder,
double **yt, int yOrder, double **mt, int checkHomogeneous)
// Returns success/fail.
vtkTypeBool vtkMath::SolveLeastSquares(int numberOfSamples, double **xt, int xOrder,
double **yt, int yOrder, double **mt, int checkHomogeneous)
{
// check dimensional consistency
if ((numberOfSamples < xOrder) || (numberOfSamples < yOrder))
......@@ -1026,12 +1027,12 @@ int vtkMath::SolveLeastSquares(int numberOfSamples, double **xt, int xOrder,
int i, j, k;
int someHomogeneous = 0;
int allHomogeneous = 1;
bool someHomogeneous = 0;
bool allHomogeneous = 1;
double **hmt = nullptr;
int homogRC = 0;
vtkTypeBool homogRC = 0;
int *homogenFlags = new int[yOrder];
int successFlag;
vtkTypeBool successFlag;
// Ok, first init some flags check and see if all the systems are homogeneous
if (checkHomogeneous)
......@@ -1045,13 +1046,13 @@ int vtkMath::SolveLeastSquares(int numberOfSamples, double **xt, int xOrder,
// Initialize homogeneous flags on a per-right-hand-side basis
for (j=0; j<yOrder; j++)
for (j=0; j<yOrder; ++j)
{
homogenFlags[j] = 1;
}
for (i=0; i<numberOfSamples; i++)
for (i=0; i<numberOfSamples; ++i)
{
for (j=0; j<yOrder; j++)
for (j=0; j<yOrder; ++j)
{
if (fabs(yt[i][j]) > VTK_SMALL_NUMBER)
{
......@@ -1079,7 +1080,7 @@ int vtkMath::SolveLeastSquares(int numberOfSamples, double **xt, int xOrder,
}
else
{
for (j=0; j<yOrder; j++)
for (j=0; j<yOrder; ++j)
{
if (homogenFlags[j])
{
......@@ -1094,7 +1095,7 @@ int vtkMath::SolveLeastSquares(int numberOfSamples, double **xt, int xOrder,
{
// hmt is the homogeneous equation version of mt, the general solution.
hmt = new double *[xOrder];
for (j=0; j<xOrder; j++)
for (j=0; j<xOrder; ++j)
{
// Only allocate 1 here, not yOrder, because here we're going to solve
// just the one homogeneous equation subset of the entire problem
......@@ -1110,37 +1111,37 @@ int vtkMath::SolveLeastSquares(int numberOfSamples, double **xt, int xOrder,
double **XXt = new double *[xOrder]; // size x by x
double **XXtI = new double *[xOrder]; // size x by x
double **XYt = new double *[xOrder]; // size x by y
for (i = 0; i < xOrder; i++)
for (i = 0; i < xOrder; ++i)
{
XXt[i] = new double[xOrder];
XXtI[i] = new double[xOrder];
for (j = 0; j < xOrder; j++)
for (j = 0; j < xOrder; ++j)
{
XXt[i][j] = 0.0;
XXtI[i][j] = 0.0;
}
XYt[i] = new double[yOrder];
for (j = 0; j < yOrder; j++)
for (j = 0; j < yOrder; ++j)
{
XYt[i][j] = 0.0;
}
}
// first find the pseudoinverse matrix
for (k = 0; k < numberOfSamples; k++)
for (k = 0; k < numberOfSamples; ++k)
{
for (i = 0; i < xOrder; i++)
for (i = 0; i < xOrder; ++i)
{
// first calculate the XXt matrix, only do the upper half (symmetrical)
for (j = i; j < xOrder; j++)
for (j = i; j < xOrder; ++j)
{
XXt[i][j] += xt[k][i] * xt[k][j];
}
// now calculate the XYt matrix
for (j = 0; j < yOrder; j++)
for (j = 0; j < yOrder; ++j)
{
XYt[i][j] += xt[k][i] * yt[k][j];
}
......@@ -1148,9 +1149,9 @@ int vtkMath::SolveLeastSquares(int numberOfSamples, double **xt, int xOrder,
}
// now fill in the lower half of the XXt matrix
for (i = 0; i < xOrder; i++)
for (i = 0; i < xOrder; ++i)
{
for (j = 0; j < i; j++)
for (j = 0; j < i; ++j)