Foreword
C++ complete code for calculating obb bounding box, rendering using osg. Refer to how to generate an OBB (OrientedboundingBox) directional bounding box. The code flow is similar, but my code is smoother and includes the calling part. The matrix class in the original link is not suitable for C++, so some of it is replaced by arrays in this article. In addition, Math::Vec3f in the article is equivalent to osg::Vec3f. I haven’t had time to sort out the code carefully yet. I can optimize it again when I have time later. I’m not very proficient in mathematics. Please point out any errors or omissions.
The first step is to calculate the covariance matrix
/*Calculate covariance matrix */
std::vector<std::vector<float>> computeConvarienceMatrix(const std::vector<Math::Vec3f> & amp; pVertices)
{<!-- -->
int numVertices = pVertices.size();
std::vector<Math::Vec3f> pVectors(numVertices, Math::Vec3f(0.0, 0.0, 0.0));
//Compute the average x,y,z
Math::Vec3f avg(0.0f, 0.0f, 0.0f);
for (int i = 0; i < numVertices; i + + )
{<!-- -->
pVectors[i] = (Math::Vec3f)pVertices[i];
avg + = pVertices[i];
}
avg /= numVertices;
for (int i = 0; i < numVertices; i + + )
pVectors[i] -= (avg);
//Compute the covariance
std::vector<std::vector<float>> covarianceMat(3, {<!-- --> 0.0,0.0,0.0 });
for (int c = 0; c < 3; c + + )
{<!-- -->
float data[3] = {<!-- --> 0 };
for (int r = c; r < 3; r + + )
{<!-- -->
float acc = 0.0f;
for (int i = 0; i < numVertices; i + + )
{<!-- -->
data[0] = pVectors[i].x();
data[1] = pVectors[i].y();
data[2] = pVectors[i].z();
acc + = data[c] * data[r];
}
acc /= numVertices;
covarianceMat[c][r] = acc;
covarianceMat[r][c] = acc; //Symmetric matrix
}
}
return covarianceMat;
}
The second step is to calculate the feature vector
The calculated eigenvector is the axis vector;
/*Find the eigenvector */
void jacobiSolver(std::vector<std::vector<float>> & amp; matrix, std::vector<Math::Vec3f> & amp; eVectors)
{<!-- -->
float eps = (float)1.0E-5;
float p, q, spq;
float cosa = 0, sina = 0; //cos(alpha) and sin(alpha)
floattemp;
float s1 = 0.0f; //sums of squares of diagonal
float s2; //elements
bool flag = true; //determine whether to iterate again
Math::Vec3f eValue(0, 0, 0);
float data[3];
float t[3][3] = {<!-- -->
1.0f,0,0,
0,1.0f,0,
0,0,1.0f };//To store the product of the rotation matrices.
do
{<!-- -->
for (int i = 0; i < 2; i + + )
{<!-- -->
for (int j = i + 1; j < 3; j + + )
{<!-- -->
if (std::abs(matrix[i][j]) < eps)
matrix[i][j] = 0.0f;
else
{<!-- -->
q = std::abs(matrix[i][i] - matrix[j][j]);
if (q > eps)
{<!-- -->
p = (2.0f * matrix[i][j] * q) / (matrix[i][i] - matrix[j][j]);
spq = (float)std::sqrt(p * p + q * q);
cosa = (float)std::sqrt((1.0f + q / spq) / 2.0f);
sina = p / (2.0f * cosa * spq);
}
else
sina = cosa = (std::sqrt(2.0f)) / 2; //INV_SQRT_TWO root of 2
for (int k = 0; k < 3; k + + )
{<!-- -->
temp = t[k][i];
t[k][i] = temp * cosa + t[k][j] * sina;
t[k][j] = temp * sina - t[k][j] * cosa;
}
for (int k = i; k < 3; k + + )
{<!-- -->
if (k > j)
{<!-- -->
temp = matrix[i][k];
matrix[i][k] = cosa * temp + sina * matrix[j][k];
matrix[j][k] = sina * temp - cosa * matrix[j][k];
}
else
{<!-- -->
data[k] = matrix[i][k];
matrix[i][k] = cosa * data[k] + sina * matrix[k][j];
if (k == j)
matrix[j][k] = sina * data[k] - cosa * matrix[j][k];
}
}
data[j] = sina * data[i] - cosa * data[j];
for (int k = 0; k <= j; k + + )
{<!-- -->
if (k <= i)
{<!-- -->
temp = matrix[k][i];
matrix[k][i] = cosa * temp + sina * matrix[k][j];
matrix[k][j] = sina * temp - cosa * matrix[k][j];
}
else
matrix[k][j] = sina * data[k] - cosa * matrix[k][j];
}
}
}
}
s2 = 0.0f;
for (int i = 0; i < 3; i + + )
{<!-- -->
eValue[i] = matrix[i][i];
s2 + = eValue[i] * eValue[i];
}
if (std::abs(s2) < (float)1.e-5 || std::abs(1 - s1 / s2) < eps)
flag = false;
else
s1 = s2;
} while (flag);
eVectors[0].set(t[0][0], t[0][1], t[0][2]);
eVectors[1].set(t[1][0], t[1][1], t[1][2]);
eVectors[2].set(t[2][0], t[2][1], t[2][2]);
Math::Vec3f cross;
cross = eVectors[0].cross(eVectors[1]);
if (cross.dot(eVectors[2]) < 0.0f)
eVectors[2] = -eVectors[2];
}
The third step: Schmidt orthogonalization of eigenvectors
/*Orthogonalization */
void schmidtOrthogonal(Math::Vec3f & amp; v0, Math::Vec3f & amp; v1, Math::Vec3f & amp; v2)
{<!-- -->
v0.normalize();
v1 = v1 - ((v1 * v0) * v0);
v1.normalize();
v2 = v0.cross(v1);
}
Step 4: Calculate the center point and half-length section
With the center point and half length and width, each vertex of the bounding box can be calculated
//Calculate center point and half-field width
void computeCenterNHalflength(std::vector<Math::Vec3f> & vecVertex, osg::Vec3f & amp; leftbottom, osg::Vec3f & amp; rightTop, osg::Vec3f & amp; rightbottom, osg::Vec3f & amp; leftTop)
{<!-- -->
std::vector<std::vector<float>> conMatrix3f = computeConvarienceMatrix(vecVertex);
std::vector<Math::Vec3f> eVectors(3, Math::Vec3f(0, 0, 0)); //Feature vector
jacobiSolver(conMatrix3f, eVectors);
schmidtOrthogonal(eVectors[0], eVectors[1], eVectors[2]);
float infinity = Base::POS_INFINITY32;
Math::Vec3f minExtents(infinity, infinity, infinity);
Math::Vec3f maxExtents(-infinity, -infinity, -infinity);
std::vector<Math::Vec3f> vertices = vecVertex;
int numVertices = vecVertex.size();
//Find the center point position
Math::Vec3f position(0.0f, 0.0f, 0.0f);
for (int i = 0; i < numVertices; i + + )
{<!-- -->
position + = vertices[i];
}
position /= numVertices;
for (int index = 0; index < numVertices; index + + )
{<!-- -->
Math::Vec3f displacement = vertices[index] - position;
//Calculate the dot product of each point on the direction axis as the length
minExtents[0] = std::min(minExtents.x(), displacement.dot(eVectors[0]));
minExtents[1] = std::min(minExtents.y(), displacement.dot(eVectors[1]));
minExtents[2] = std::min(minExtents.z(), displacement.dot(eVectors[2]));
maxExtents.x() = std::max(maxExtents.x(), displacement.dot(eVectors[0]));
maxExtents.y() = std::max(maxExtents.y(), displacement.dot(eVectors[1]));
maxExtents.z() = std::max(maxExtents.z(), displacement.dot(eVectors[2]));
}
Math::Vec3f offset = maxExtents - minExtents;
offset /= (2.0f);
offset + = (minExtents);
//Correction center point
position + = (eVectors[0] * (offset.x()));
position + = (eVectors[1] * (offset.y()));
position + = (eVectors[2] * (offset.z()));
Math::Vec3f halfExtents;
halfExtents.x() = (maxExtents.x() - minExtents.x()) / 2.0f;
halfExtents.y() = (maxExtents.y() - minExtents.y()) / 2.0f;
halfExtents.z() = (maxExtents.z() - minExtents.z()) / 2.0f;
Math::Vec3f mleftbtm = position - halfExtents.x() * eVectors[0] - halfExtents.z() * eVectors[2];
Math::Vec3f mrightbottom = position + halfExtents.x() * eVectors[0] - halfExtents.z() * eVectors[2];
Math::Vec3f mrightTop = position + halfExtents.x() * eVectors[0] + halfExtents.z() * eVectors[2];
Math::Vec3f mleftTop = position - halfExtents.x() * eVectors[0] + halfExtents.z() * eVectors[2];
leftbottom.set(mleftbtm.x(), mleftbtm.y(), mleftbtm.z());
rightTop.set(mrightTop.x(), mrightTop.y(), mrightTop.z());
rightbottom.set(mrightbottom.x(), mrightbottom.y(), mrightbottom.z());
leftTop.set(mleftTop.x(), mleftTop.y(), mleftTop.z());
}
Calling code
void updateRectBox(std::vector<Math::Vec3f> vecVertex)
{<!-- -->
osg::Vec3f leftbottom, rightTop, rightbottom, leftTop;
computeCenterNHalflength(vecVertex, leftbottom, rightTop, rightbottom, leftTop);
//Draw the border
osg::ref_ptr<osg::Geometry> geomLine = new osg::Geometry;
osg::ref_ptr<osg::Vec3Array> temVertex = new osg::Vec3Array();
temVertex->push_back(leftbottom);
temVertex->push_back(rightbottom);
temVertex->push_back(rightTop);
temVertex->push_back(leftTop);
geomLine->setVertexAttribArray(0, temVertex);
geomLine->setVertexAttribBinding(0, osg::Geometry::BIND_PER_VERTEX);
//color array
osg::Vec4Array* vColors = new osg::Vec4Array;
vColors->push_back(COLOR_RECT);
geomLine->setColorArray(vColors);
geomLine->setColorBinding(osg::Geometry::BIND_OVERALL);
//Draw primitive attributes
osg::ref_ptr<osg::DrawArrays> vPri = new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP, 0, temVertex->size());
geomLine->addPrimitiveSet(vPri.get());
osg::Geode* pRectGeode = new osg::Geode;
pRectGeode->getOrCreateStateSet()->setMode(GL_DEPTH_TEST, osg::StateAttribute::OFF);
pRectGeode->addDrawable(geomLine);
}