3 points define a circle

Greetings members of Second Life! My geometry and trigonometry have been out of use for several years, and I was hoping someone on the forums could help shed some insight for me.

I know that with 3 points, you can define a circle. Would anyone here happen to have some code that has a function in which you feed the function the three points (vectors), and it would calculate the centre of the circle and radius?

Thanks in advance!

bisect any two of the sides with perpendicular lines,, where those two lines meet is the center of the circle, distance from there to any point is the radius

Or, only slightly more complicated, a perpendicular to any chord at its midpoint also goes through the center of the circle. With any three points on a circle's circumference, you automatically have three chords defined. Take each one, halve its length, and draw a perpendicular at that point. See where they cross. Translating that geometric solution into an algebraic one is a nice puzzle for an afternoon.

Three points define a PLANE, not a circle.

That's true, but they also automatically lie on the circumference of a unique circle in that plane (and at the corners of a unique triangle as well, for that matter), so the OP's question is a perfectly good one.

I think the Looprez scripts have that formula. (o.o)

In space 3 points will define multiple spheres, if not on a line.
In the plane they will define a circle if they are not on a line.
In the plane you may use the circle equation: (x-a)^2+(y-b)^2=r^2
to find center and radius.
I did it once
Use the circle equation 3 times with the 3 points inserted:
P1=(a1,b1)
P2=(a2,b2)
P3=(a3,b3)
Then solve the 3 equations for x, y and r
Thats all
(if my memory serves me well)

Wins today's award for clarity and brevity!

Now, write a code that does that!

llSay("Wins today's award for clarity and brevity!";

Rest of code is an exercise for the student.

Average the points in space to find the center. The center to any of the three points forms the radius... draw an arc around the center of that radius. Er, on the plane. Which I assumed the points where on.

Plane Geometry flashback, whoa! I knew that class would come n handy in the real world, NOT!

Here you go. Given the points P1, P2, and P3, determine the center point C:

All points are the same distance from C:

   (C-P1)^2 = (C-P2)^2 = (C-P3)^2

   Define vector U = C-P1

   (C-P1)^2 = U^2
   (C-P2)^2 = [U-(P2-P1)]^2 = U^2 + (P2-P1)^2 - 2(P2-P1)*U
   (C-P3)^2 = [U-(P3-P1)]^2 = U^2 + (P3-P1)^2 - 2(P3-P1)*U

   1.) 2(P2-P1)*U = (P2-P1)^2
   2.) 2(P3-P1)*U = (P3-P1)^2

C is in the same plane as P1, P2, and P3 (otherwise there'd be an infinite line satisfying 1. and 2.):

   N = (P2-P1)%(P3-P1)

   3.) N*U = 0

By 1., 2., and 3. we have three simultaneous linear equations for U, expressible in matrix form as:

   [    2(P2-P1)     ]       [ (P2-P1)^2  ]
   [    2(P3-P1)     ] * U = [ (P3-P1)^2  ]
   [ (P2-P1)%(P3-P1) ]       [     0      ]

   C = P1+U

Now over in /54/82/298475/1.html you'll find I just published a library of functions useful for solving such systems of linear equations. With those functions we should be able to solve that matrix equation as follows:

// #include "MHLSL 3D Matrix Library"

// Finds a circle with the given points on its circumference.
// params p1, p2, p3 - The points on the circumference of the circle.
// returns - A list composed of [ centerPoint, planeNormalVector, radius ], or
//           an empty list if the center cannot be determined (e.g. if the
//           points are colinear).
list getCommonCircle(vector p1, vector p2, vector p3)
{
   vector side1 = p2-p1;
   vector side2 = p3-p1;
   vector normal = side1%side2;

   // setup the matrix equation

   list m =
      [ 2*side1,
        2*side2,
        normal ];
   vector v = <side1*side1, side2*side2, 0.0>;

   // solve:

   list mInv = matrixInverse(m);
   if (llGetListLength(mInv) == 0)
   {
      // No inverse; points are colinear
      return [];
   }

   vector u = matrixVectorProduct(mInv, v);
   vector c = p1+u;
   float r = llSqrt(u*u);

   return [ c, llVecNorm(normal), r ];
}

Average won't give you the center of a circle going through those three points, unfortunately. It's just a tad more complicated than that.

The radius computation is the easy part.

and do all the work? nah, I'm not at the computer with all the reference links. just make sure you get all the angles of the rays right, because otherwise youll have problems... just because a circle lies in a plane doesnt mean it lies in a convenient plane like x,y,z, or the inscribed triangle is even aligned to axis either. assuming both of those are true it's dog simple, asumming both are false, it's a major pain in lsl

small shorcut: the center will never lie outside the bounds of a right triangle, and never inside the bounds of a triangle with 2 angles < 45deg

free hint for alternate method. compress your triangle by dropping two axis seperately and comparing.

Seriously, the code snippet I posted above works. Here is a full script utilizing it. Drop it in a cylinder, rez 3 prims called "CircumferencePoint" around it, and touch it when you want it to automatically move.

// "MHLSL 3D Matrix Library" by Hewee "MH" Zetkin
// License:
//    The author hereby grants unlimited permission to copy, modify, and
//    redistribute this script in any form, with the exception that any
//    source code listing must retain this header giving credit to the author,
//    similar licensing, and a basic summary of the changes made along with the
//    contributors who made them.
// Description:
//    Provides basic 3D (3x3, 3x1) matrix and matrix-vector operations in LSL,
//    similar to the functionality found in BLAS levels 2 and 3 (with level 1
//    provided by the LSL language itself) plus matrix transposition and
//    inversion.
// Release Notes:
//    v0.1.1 (Dec. 18, 2008)
//       - List of row-vectors matrix encoding.
//       - Basic matrix functions, matrix product, matrix-vector product
//       - Matrix inverse
//       - Matrix inverse test harness
//       - Published to SL Scripting Tips forums in hopes of testing for
//         version 1.0.
//    v0.1.2 (Dec. 18, 2008)
//       - Fixed matrix-vector product return type.

// "Close enough" floating point value for zero comparison and such.
float EPSILON = 0.001;

// The 3x3 identity matrix in row-major order, with each row represented as a
// vector.
list MATRIX_IDENTITY =
   [ <1.0, 0.0, 0.0>,
     <0.0, 1.0, 0.0>,
     <0.0, 0.0, 1.0> ];

// Computes the sum of two 3x3 matrix arguments.
// See http://en.wikipedia.org/wiki/Matrix_product
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// param n - A 3x3 matrix encoded like 'm';
// returns - A matrix encoded like 'm' that is the matrix sum of 'm' and
//           'n'.
list matrixSum(list m, list n)
{
   vector mr1 = llList2Vector(m, 0);
   vector mr2 = llList2Vector(m, 1);
   vector mr3 = llList2Vector(m, 2);

   vector nr1 = llList2Vector(n, 0);
   vector nr2 = llList2Vector(n, 1);
   vector nr3 = llList2Vector(n, 2);

   return [ mr1+nr1, mr2+nr2, mr3+nr3 ];
}

// Computes the difference of two 3x3 matrix arguments.
// See http://en.wikipedia.org/wiki/Matrix_product
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// param n - A 3x3 matrix encoded like 'm';
// returns - A matrix encoded like 'm' that is the matrix difference of 'm' and
//           'n'.
list matrixDiff(list m, list n)
{
   vector mr1 = llList2Vector(m, 0);
   vector mr2 = llList2Vector(m, 1);
   vector mr3 = llList2Vector(m, 2);

   vector nr1 = llList2Vector(n, 0);
   vector nr2 = llList2Vector(n, 1);
   vector nr3 = llList2Vector(n, 2);

   return [ mr1-nr1, mr2-nr2, mr3-nr3 ];
}

// Computes the product of a 3x3 matrix and a scalar.
// See http://en.wikipedia.org/wiki/Matrix_(mathematics)#Basic_operations
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// param s - A scalar.
// returns - A matrix encoded like 'm' that is the product of 'm' and 's'.
list matrixScalarProduct(list m, float s)
{
   vector r1 = llList2Vector(m, 0);
   vector r2 = llList2Vector(m, 1);
   vector r3 = llList2Vector(m, 2);

   return [ s*r1, s*r2, s*r3 ];
}

// Transposes a 3x3 matrix.
// See http://en.wikipedia.org/wiki/Transpose
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// returns - A matrix encoded like 'm' that is the transpose of 'm'.
list matrixTranspose(list m)
{
   vector r1 = llList2Vector(m, 0);
   vector r2 = llList2Vector(m, 1);
   vector r3 = llList2Vector(m, 2);

   return [ <r1.x, r2.x, r3.x>,
            <r1.y, r2.y, r3.y>,
            <r1.z, r2.z, r3.z> ];
}

// Computes the product of two 3x3 matrix arguments.
// See http://en.wikipedia.org/wiki/Matrix_product
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// param n - A 3x3 matrix encoded like 'm';
// returns - A matrix encoded like 'm' that is the matrix product of 'm' and
//           'n'.
list matrixProduct(list m, list n)
{
   vector mr1 = llList2Vector(m, 0);
   vector mr2 = llList2Vector(m, 1);
   vector mr3 = llList2Vector(m, 2);

   vector nr1 = llList2Vector(n, 0);
   vector nr2 = llList2Vector(n, 1);
   vector nr3 = llList2Vector(n, 2);

   // columns
   vector nc1 = <nr1.x, nr2.x, nr3.x>;
   vector nc2 = <nr1.y, nr2.y, nr3.y>;
   vector nc3 = <nr1.z, nr2.z, nr3.z>;

   // product
   float p11 = mr1*nc1; float p12 = mr1*nc2; float p13 = mr1*nc3;
   float p21 = mr2*nc1; float p22 = mr2*nc2; float p23 = mr2*nc3;
   float p31 = mr3*nc1; float p32 = mr3*nc2; float p33 = mr3*nc3;

   return [ <p11, p12, p13>,
            <p21, p22, p23>,
            <p31, p32, p33> ];
}

// Computes the product of a 3x3 matrix and a 3x1 matrix (vector).
// See http://en.wikipedia.org/wiki/Matrix_product
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// param v - A 3x1 matrix encoded as a vector.
// returns - A vector that is the product of 'm' and 'v'.
vector matrixVectorProduct(list m, vector v)
{
   vector mr1 = llList2Vector(m, 0);
   vector mr2 = llList2Vector(m, 1);
   vector mr3 = llList2Vector(m, 2);

   return <mr1*v, mr2*v, mr3*v>;
}

// Computes the determinant of the 3x3 matrix argument.
// See http://en.wikipedia.org/wiki/Determinant
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// returns - The determinant of 'm'.
float matrixDeterminant(list m)
{
   vector r1 = llList2Vector(m, 0);
   vector r2 = llList2Vector(m, 1);
   vector r3 = llList2Vector(m, 2);

   return r1*(r2%r3);
}

// Computes the cofactor matrix of the 3x3 matrix argument.
// See http://en.wikipedia.org/wiki/Cofactor_(linear_algebra)#Matrix_of_cofactors
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// returns - A matrix encoded like 'm' that is composed of the cofactors of
//           'm'.
list matrixCofactors(list m)
{
   vector r1 = llList2Vector(m, 0);
   vector r2 = llList2Vector(m, 1);
   vector r3 = llList2Vector(m, 2);

   float m11 = r1.x; float m12 = r1.y; float m13 = r1.z;
   float m21 = r2.x; float m22 = r2.y; float m23 = r2.z;
   float m31 = r3.x; float m32 = r3.y; float m33 = r3.z;

   // cofactors
   float c11 =  m22*m33-m32*m23;
   float c12 = -m21*m33+m31*m23;
   float c13 =  m21*m32-m31*m22;
   float c21 = -m12*m33+m32*m13;
   float c22 =  m11*m33-m31*m13;
   float c23 = -m11*m32+m31*m12;
   float c31 =  m12*m23-m22*m13;
   float c32 = -m11*m23+m21*m13;
   float c33 =  m11*m22-m21*m12;

   return [ <c11, c12, c13>,
            <c21, c22, c23>,
            <c31, c32, c33> ];
}

// Computes the adjugate matrix of the 3x3 matrix argument.
// See http://en.wikipedia.org/wiki/Cofactor_(linear_algebra)#Adjugate
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// returns - A matrix encoded like 'm' that is the adjugate matrix to 'm'.
list matrixAdjugate(list m)
{
   return matrixTranspose(matrixCofactors(m));
}

// Inverts a 3x3 matrix.
// See http://en.wikipedia.org/wiki/Matrix_inverse
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// returns - A matrix encoded like 'm' that is the inverse of 'm'.  If 'm' is
//           not invertable, returns an empty list instead.
list matrixInverse(list m)
{
   float det = matrixDeterminant(m);
   if (llFabs(det) < EPSILON)
   {
      return [];
   }

   list adj = matrixAdjugate(m);
   vector ar1 = llList2Vector(adj, 0);
   vector ar2 = llList2Vector(adj, 1);
   vector ar3 = llList2Vector(adj, 2);

   return [ ar1/det, ar2/det, ar3/det ];
}

// Converts a matrix to a human-readable string for debugging output.
// param m - A 3x3 matrix in row-major order, each row represented as a vector
//           in the list (so 'm11' is m[0].x and 'm23' is m[1].z).  In other
//           words:
//              [ <m11, m12, m13>,
//              [ <m21, m22, m23>,
//              [ <m31, m32, m33> ]
// returns - A string formatted with rows separated by newlines and row
//           elements separated by ", ".
string matrixToString(list m)
{
   vector r1 = llList2Vector(m, 0);
   vector r2 = llList2Vector(m, 1);
   vector r3 = llList2Vector(m, 2);

   float m11 = r1.x; float m12 = r1.y; float m13 = r1.z;
   float m21 = r2.x; float m22 = r2.y; float m23 = r2.z;
   float m31 = r3.x; float m32 = r3.y; float m33 = r3.z;

   return (string)m11+", "+(string)m12+", "+(string)m13+"\n"+
          (string)m21+", "+(string)m22+", "+(string)m23+"\n"+
          (string)m31+", "+(string)m32+", "+(string)m33;
}


//// Main Script

string POINT_NAME = "CircumferencePoint";
float SENSOR_RADIUS = 10.0;
float THICKNESS = 0.1;
string NOT_FOUND_TEXT = "Insufficient points found";


// Finds a circle with the given points on its circumference.
// params p1, p2, p3 - The points on the circumference of the circle.
// returns - A list composed of [ centerPoint, planeNormalVector, radius ], or
//           an empty list if the center cannot be determined (e.g. if the
//           points are colinear).
list getCommonCircle(vector p1, vector p2, vector p3)
{
   vector side1 = p2-p1;
   vector side2 = p3-p1;
   vector normal = side1%side2;

   // setup the matrix equation

   list m =
      [ 2*side1,
        2*side2,
        normal ];
   vector v = <side1*side1, side2*side2, 0.0>;

   // solve:

   list mInv = matrixInverse(m);
   if (llGetListLength(mInv) == 0)
   {
      // No inverse; points are colinear
      return [];
   }

   vector u = matrixVectorProduct(mInv, v);
   vector c = p1+u;
   float r = llSqrt(u*u);

   return [ c, llVecNorm(normal), r ];
}


default
{
    state_entry()
    {}
    
    touch_start(integer nDetected)
    {
        llSensor(POINT_NAME, NULL_KEY, PASSIVE | ACTIVE, SENSOR_RADIUS, PI);
    }
    
    no_sensor()
    {
        llSetText(NOT_FOUND_TEXT, <1.0, 1.0, 1.0>, 1.0);
    }
    
    sensor(integer nDetected)
    {
        if (nDetected < 3)
        {
            llSetText(NOT_FOUND_TEXT, <1.0, 1.0, 1.0>, 1.0);
            return;
        }
        
        vector p1 = llDetectedPos(0);
        vector p2 = llDetectedPos(1);
        vector p3 = llDetectedPos(2);
        
        list circleProps = getCommonCircle(p1, p2, p3);
        if (llGetListLength(circleProps) < 3)
        {
            llSetText(NOT_FOUND_TEXT, <1.0, 1.0, 1.0>, 1.0);
            return;
        }
        
        vector center = llList2Vector(circleProps, 0);
        vector normal = llList2Vector(circleProps, 1);
        float radius = llList2Float(circleProps, 2);
        
        if (radius < 0.01 || radius > 10.0)
        {
            llSetText(NOT_FOUND_TEXT, <1.0, 1.0, 1.0>, 1.0);
            return;
        }
        
        rotation rot = llRotBetween(<0.0, 0.0, 1.0>, normal);
        
        llSetText("", ZERO_VECTOR, 0.0);
        
        llSetPrimitiveParams(
            [ PRIM_POSITION, center,
              PRIM_ROTATION, rot,
              PRIM_SIZE, <2.0*radius, 2.0*radius, THICKNESS> ]);
    }
}

Hewee, thank you for the code! I lack the higher level of mathematics than is required (more specifically, matrix math). I have searched the web and found examples, but all the examples I seen applicable to putting it into code involved matrix math.

I found this youtube link which shows an example of the math, but after a point I was just lost on how the answers found. (http://www.youtube.com/watch?v=2uS6OFpKQzQ).

I also seen numerous examples of the geometry involved (perpendicular bisectors of the chords, etc), but Hewee's soloution is just so elegant. I used to code Hewee provided and it works just like a charm!

Thanks to all for your answers! It has shed light on where I need to grow! Hewee, your contribution of the MHLSL 3D Matrix Library is simply awesome...I plan to it to learn a new thing or two!

I just love this scripting forum. I have been lurking here on and off and the things I learn are just great.

LOL, yeah. I was pondering it earlier... and had a WTF was I thinking moment.

...this is an awful long thread considering it was answered by the first response (Void's).

/esc

It was not answered until Hewee's code.
The OP asked for a code, remember?

props to heewee for the code, though I do think the whole thing could be simplified by planar compression... something to try later when I have my good computer back and set up

Just out of curiosity, what would happen if you found the average of the three points? Would that equal the center?


As far as I know, it would work, but I've never tried anything like that with vectors before.



*EDIT* Never mind, it appears this has been proven incorrect ^^ *EDIT*

Errrr guess none of you skipped class and spent days down at the bus stop smoking.......

average will find the center of the triangle... which can be useful in other ways.... but the only case that it'll find the center of the circle is an equilateral triangle.

For a quick visual proof, imagine 3 points very close together on the perimiter of a large circle. The center would be just inside the permiter, nowhere near the middle of the circle.