Eigen Vector Calculations for OBB crashing on second pass through "Assertion Failed"

Question

Sorry if the title isn't clear enough, I wasn't really sure how to summarize this problem since it's very specific and odd. Basically, I have a mesh that I've drawn and I am attempting to calculate an Oriented Bounding Box. In order to do this I have used Principle Component Analysis. After computing the Eigenvectors, I find the rotation matrix that would transform the eigen vectors onto the x, y and z axis respectively. I then transform the mesh by this rotation matrix, calculate an axis aligned bounding box and then rotate both the mesh and the box by the inverse of the rotation matrix. When I run the program, it successfuly runs through the bounding box function once, although no box is visible, but on the second run of the display function it crashed when trying to compute the bounding box. I have identified the line that causes the crash and it is part of the eigenvector calculations. However, because I am using a library that I was given to compute the eigenvectors I'm not completely sure what is going on in the function and I cant identify why it's not working. When it crashes it gives the error "assertion failed" followed by the file path of my program and the line.

My function for calculating the covariances is here:

std::array<std::array<double, 3>, 3> covarianceCalc2()
{
   std::array<std::array<double, 3>, 3> sum = {{{0, 0, 0}, {0, 0, 0}, {0, 0, 0,}}};
   std::array<double, 3> tempVec;
   double mean;
   for(int i = 0; i < meshVertices.size(); i++)
   {
     mean = (meshVertices[i].x + meshVertices[i].y + meshVertices[i].z)/3;
     tempVec[0] = meshVertices[i].x - mean;
     tempVec[1] = meshVertices[i].y - mean;
     tempVec[2] = meshVertices[i].z - mean;
     sum = matrixAdd(sum, vectorTranposeMult(tempVec));
   }
 sum = matrixMultNum(sum, 1/(meshVertices.size()));
 return sum;
}

Here is the code that computes the eigenvectors and rotates them:

 Compute_EigenV(covarianceCalc2(), eigenValues, eigenVectors_1, eigenVectors_2, eigenVectors_3);
 std::array<double, 3> x = {1, 0, 0};
 std::array<double, 3> y = {0, 1, 0};
 std::array<double, 3> z = {0, 0, 1};
 std::array<std::array<double, 3>, 3> transformX = findRotation(eigenVectors_1, x);
 std::array<std::array<double, 3>, 3> transformY = findRotation(eigenVectors_2, y);
 std::array<std::array<double, 3>, 3> transformZ = findRotation(eigenVectors_3, z);


 memcpy(transformX, axisRotateX, sizeof(double) * 3 * 3);
 memcpy(transformY, axisRotateY, sizeof(double) * 3 * 3);
 memcpy(transformZ, axisRotateZ, sizeof(double) * 3 * 3);

The first line of this section is the line that causes the crash. The Compute_Eigen function is as follows:

void Compute_EigenV(std::array<std::array<double, 3>, 3> covariance, std::array<double, 3> eigenValues, std::array<double, 3> eigenVectors_1, std::array<double, 3> eigenVectors_2, std::array<double, 3> eigenVectors_3)
{
    printf("Matrix Stuff\n");
    MatrixXd m(3, 3);
    m << covariance[0][0], covariance[0][1], covariance[0][2],
         covariance[1][0], covariance[1][1], covariance[1][2],
         covariance[2][0], covariance[2][1], covariance[2][2];

    // volving SVD
    printf("EigenSolver\n");
    EigenSolver<MatrixXd> solver(m);
    MatrixXd all_eigenVectors = solver.eigenvectors().real();
    MatrixXd all_eigenValues = solver.eigenvalues().real();

    // find the max index
    printf("Find Max Index\n");
    int INDEX[3];
    double max;
    max=all_eigenValues(0,0);
    int index=0;
    for (int i=1;i<3;i++){
        if (max<all_eigenValues(i,0)){
            max=all_eigenValues(i,0);
            index=i;
        }
    }
    INDEX[0]=index;

    // find the min index
    printf("Find Min Index\n");
    double min;
    min=all_eigenValues(0,0);

    index=0;
    for (int i=1;i<3;i++){
        if (min>all_eigenValues(i,0)){
            min=all_eigenValues(i,0);
            index=i;
        }
    }
    INDEX[1]=3-index-INDEX[0];
    INDEX[2]=index;

    // giave eigenvalues and eien vectors to matrix
    printf("Give values and vector to matrix\n");
    eigenValues[0]=all_eigenValues(INDEX[0],0);
    printf("1");
    eigenValues[1]=all_eigenValues(INDEX[1],0);
    printf("1\n");
    eigenValues[2]=all_eigenValues(INDEX[2],0);

    printf("Vector 1\n");
    VectorXd featureVector_1 = all_eigenVectors.col(INDEX[0]);
    eigenVectors_1[0]=featureVector_1(0);
    eigenVectors_1[1]=featureVector_1(1);
    eigenVectors_1[2]=featureVector_1(2);

    printf("Vector 2\n");
    VectorXd featureVector_2 = all_eigenVectors.col(INDEX[1]);
    eigenVectors_2[0]=featureVector_2(0);
    eigenVectors_2[1]=featureVector_2(1);
    eigenVectors_2[2]=featureVector_2(2);

    printf("Vector 3\n");
    VectorXd featureVector_3 = all_eigenVectors.col(INDEX[2]);
    eigenVectors_3[0]=featureVector_3(0);
    eigenVectors_3[1]=featureVector_3(1);
    eigenVectors_3[2]=featureVector_3(2);

}

The specific section at which the program crashes is:

    eigenValues[0]=all_eigenValues(INDEX[0],0);
    eigenValues[1]=all_eigenValues(INDEX[1],0);
    eigenValues[2]=all_eigenValues(INDEX[2],0);

Using go to definition on all_eigenValues shows this section of code:

/** \returns a reference to the coefficient at given the given row and column.
  *
  * \sa operator[](Index)
  */

EIGEN_STRONG_INLINE Scalar&
operator()(Index row, Index col)
{
  eigen_assert(row >= 0 && row < rows()
      && col >= 0 && col < cols());
  return derived().coeffRef(row, col);
}

Unfortunately I'm not experienced enough with C to really figure out what is causing the error here. Because it might be relevant, here is my function that computes the rotation matrix:

std::array<std::array<double, 3>, 3> findRotation(std::array<double, 3> vec1, std::array<double, 3> vec2)
{
std::array<std::array<double, 3>, 3> tempMatrix, tempMatrix2;
std::array<double, 3> crossProd = crossProduct(vec1, vec2);
double dotProd = dotProduct(vec1, vec2);
double sinAngle = modVector(crossProd)/(modVector(vec1) * modVector(vec2));
std::array<std::array<double, 3>, 3> xProdMatrix = {{{0, crossProd[2] * -1, crossProd[1]}, {crossProd[2], 0, crossProd[0] * -1}, {crossProd[1] * -1, crossProd[0], 0}}}; 
tempMatrix = matrixAdd(iMatrix, xProdMatrix);
tempMatrix2 = matrixMult(xProdMatrix, xProdMatrix);
tempMatrix2 = matrixMultNum(tempMatrix2, (1 - dotProd)/pow(sinAngle, 2));
tempMatrix = matrixAdd(tempMatrix, tempMatrix2);
return tempMatrix;
}

If you need any more code segments I can provide them, I've hit a wall with this issue and haven't made any progress in hours so any help would be appreciated.

I'm assuming that since the eigen calculation is from a library that the error is mine and not with the function and if I had to guess where my mistake was it would be in the rotation matrix calculation but I've not been able to spot the error myself if that is where it is.

Edit: Ok I've fixed the issue with my covariance calculations I hope but unfortunately the same crash still occurs. I've also changed all my arrays to std::arrays as per the advice of another commenter. I've been at this for 7 hours now.

Answer 1

I believe I found your problem.

You need to check again the theory! As I may recall, you have the covariance defined in the theory as:

C=1/M \Sum ( (p-pmean)*(p-pmean)^t )

Well, you may notice that C is a 3x3 matrix, NOT a value. Therefore, when you call Compute_EigenV and it tries to make the SVD of a 3x3 matrix, it crashes.

EDIT Looks that you solved and explained that part a bit , very good. The assertion is part of eigen library and as you can see:

  eigen_assert(row >= 0 && row < rows()   &&  col >= 0 && col < cols());

Its just making sure that the input is a 3x3 matrix (in your case). Somehow this is not happening. I am pretty sure the code of Compute_EigenV works without modifications because well, my co-worker wrote it and I have tried it, therefore it must be something wrong with how you declared the variables you are imputing to the function.

---

Again, that code looks familiar ;).