Creating a circle in 3D space using VEX

[nsabesan]

Hey everyone, 

I am in the process of creating a digital asset and its basic building block is a circle. First I have the user plot three points in 3D space. And from these 3 points which form a triangle, I am creating a circle that circumscribes this triangle. 

The algorithm I have so far from gathered resources are below;

1.) calculate the normal vector to the triangle formed by the three user-defined points.

2.) use a reference plane into which this triangle will rotate onto, say the XY plane

3.) rotate the three points into this reference plane after calculating and building a rotation matrix

4.) find the center of the circumscribing circle from the three points

5.) calculate the radius of the circle which is easy

6.) build the circle on the reference plane

7.) build a reverse rotation matrix to rotate the points back to the original plane

I'm not sure what I am doing wrong but the code I have so far works only if the points lie on a plane in 2D space. The minute the points have values in the 3'rd dimension, it gets wonky. I'm also not sure if the rotation matrix in Houdini does the rotation about the origin or about the pivot of the three points. I can share the code to take a look at it.

If anyone has any solutions or suggestions, kindly let me know! Much appreciated!

[deniz]

code ? 

Edited July 22, 2016 by deniz

[nsabesan]

sry here is the code

vector CalcCircleCenter(vector pt1, pt2, pt3, center) {
    float yDelta_a = pt2.y - pt1.y;
    float xDelta_a = pt2.x - pt1.x;
    float yDelta_b = pt3.y - pt2.y;
    float xDelta_b = pt3.x - pt2.x;
    
    if(abs(xDelta_a) <= 0.000000001 && abs(yDelta_b) <= 0.000000001) {
        center.x = 0.5 * (pt2.x + pt3.x);
        center.y = 0.5 * (pt1.y + pt2.y);
        center.z = pt1.z;
        return 1;
    }
    
    //the IsPerpendicular() assures that xDelta(s) are not zero
    float aSlope = yDelta_a / xDelta_a;
    float bSlope = yDelta_b / xDelta_b;
    
    //checking whether the given points are colinear
    if(abs(aSlope - bSlope) <= 0.000000001) {
        return -1;
    }
    
    //calculate the center
    center.x = (aSlope * bSlope * (pt1.y - pt3.y) + bSlope * (pt1.x + pt2.x) - aSlope * (pt2.x + pt3.x)) / (2 * (bSlope - aSlope));
    center.y = -1 * (center.x - (pt1.x + pt2.x) / 2) / aSlope + (pt1.y + pt2.y) / 2;

    return center;
}

//check if the given point is perpendicular to x or y axis
int IsPerpendicular(vector pt1, pt2, pt3) {
    float yDelta_a= pt2.y - pt1.y;
    float xDelta_a= pt2.x - pt1.x;
    float yDelta_b= pt3.y - pt2.y;
    float xDelta_b= pt3.x - pt2.x;

    //check whether the line of the two pts are vertical
    if (abs(xDelta_a) <= 0.000000001 && abs(yDelta_b) <= 0.000000001){
        return false;
    }

    if (abs(yDelta_a) <= 0.000000001){
        return true;
    }
    else if (abs(yDelta_b) <= 0.000000001){
        return true;
    }
    else if (abs(xDelta_a) <= 0.000000001){
        return true;
    }
    else if (abs(xDelta_b) <= 0.000000001){
        return true;
    }
    else
        return false ;
}

//find circle center based on a check
void FindCircleCenter(vector pt1, pt2, pt3, center) {
    if(!IsPerpendicular(pt1, pt2, pt3)) {
        CalcCircleCenter(pt1, pt2, pt3, center);
    }
    else if(!IsPerpendicular(pt1, pt3, pt2)) {
        CalcCircleCenter(pt1, pt3, pt2, center);
    }
    else if(!IsPerpendicular(pt2, pt1, pt3)) {
        CalcCircleCenter(pt2, pt1, pt3, center);
    }
    else if(!IsPerpendicular(pt2, pt3, pt1)) {
        CalcCircleCenter(pt2, pt3, pt1, center);
    }
    else if(!IsPerpendicular(pt3, pt2, pt1)) {
        CalcCircleCenter(pt3, pt2, pt1, center);
    }
    else if(!IsPerpendicular(pt3, pt1, pt2)) {
        CalcCircleCenter(pt3, pt1, pt2, center);
    }
    else {
        removepoint(0,0);
        removepoint(0,1);
        removepoint(0,2);
        return;
    }

    removepoint(0,0);
    removepoint(0,1);
    removepoint(0,2);
}

//get the position of the three points circumscribing the triangle
vector A = point(@OpInput1,"P",0);
vector B = point(@OpInput1,"P",1);
vector C = point(@OpInput1,"P",2);

//initializing center of circumscribing circle
vector center = {0,0,0};

//calculating normal vector to the triangle formed by the 3 points
vector Delta_a = A - B;
vector Delta_b = C - B;

vector circle_normal = cross(normalize(Delta_b), normalize(Delta_a));
vector xy_normal = {0,0,1};

//rotation axis around which the circle will rotate into XY plane
vector rot_axis = cross(normalize(xy_normal),normalize(circle_normal));
float rot_angle = acos(dot(normalize(xy_normal),normalize(circle_normal))); //in radians

//building quaternion rotation matrix
float angle = cos(rot_angle * 0.5);
rot_axis *= sin(rot_angle * 0.5);

matrix3 rot_quat = ident();
rotate(rot_quat, angle, normalize(rot_axis));
quaternion(rot_quat);

//rotate points A, B and C
A *= rot_quat;
B *= rot_quat;
C *= rot_quat;
center *= rot_quat;

FindCircleCenter(A, B, C, center);

v@offset = center;

int temp_pnt = addpoint(0,center);

//calculating radius of circle
float radius = distance(A,center);
float deg_to_rad = 3.14159 / 180.0;

//making a primitive for the circle
int prim_circle = addprim(0,"poly");

//building reverse quaternion rotation
rot_angle = -(acos(dot(normalize(xy_normal),normalize(circle_normal)))); //in radians
angle = cos(rot_angle * 0.5);

vector rot_back_axis = cross(normalize(xy_normal),normalize(circle_normal));
rot_back_axis *= sin(rot_angle * 0.5);

matrix3 rot_back_quat = ident();
rotate(rot_back_quat, angle, normalize(rot_back_axis));
quaternion(rot_back_quat);

//creating the circle
for(int i = 0; i<360; i++) {
    float deg_in_rad = i * deg_to_rad;
    center.x = cos(deg_in_rad) * radius;
    center.y = sin(deg_in_rad) * radius;
    center.z = 0;
    
    //rotate the point back to original plane
    center *= rot_back_quat;
    
    temp_pnt = addpoint(0,center);
    addvertex(0, prim_circle, temp_pnt);
}

 

[f1480187]

// Detail wrangle.

#define PI 3.1415926535897932384

vector a = point(0, "P", 0);
vector b = point(0, "P", 1);
vector c = point(0, "P", 2);

vector ac = c - a;
vector ab = b - a;
vector ab_x_ac = cross(ab, ac);
vector a_center = (cross(ab_x_ac, ab) * length2(ac) + cross(ac, ab_x_ac) * length2(ab)) / (2.0 * length2(ab_x_ac));

float radius = length(a_center);
vector normal = normalize(ab_x_ac);
vector center = a + a_center;

// Test.
vector x, y, z=normal;
makebasis(x, y, z);
matrix3 r = set(x, y, z);

int count = 50;
for (int i = 0; i < count; i++)
{
    float angle = 2 * PI / count * i;
    vector pos = set(sin(angle), cos(angle), 0) * radius;
    addpoint(0, pos * r + center);
}
addpoint(0, center);

Here is working stuff found there. Add degeneracy testing if you need it.

circumcircle.hipnc

[nsabesan]

@f1480187 

you sir are a JEDI!!!! Thanks so much!!! 

Now please reveal yourself!

Kind regards,
Naveen
SideFX Intern