geometry.frink

View or download geometry.frink in plain text format


// Formulae for simple plane geometry.

// Find the length of the hypotenuse of a right triangle given sides a and b.
hypotenuseLength[a, b] := sqrt[a^2 + b^2]

// Calculate area of arbitrary triangle given 3 sides of known length.
// Calculate using Heron's formula.
triangleArea[a, b, c] :=
{
   // Semiperimeter s
   s = 1/2 (a + b + c)

   area = sqrt[s (s-a) (s-b) (s-c)]

   return area
}


// Find the distance between a point (x0,y0) and a line specified by
// two points (x1, y1), (x2, y2)
pointToLineDistance[x0, y0, x1, y1, x2, y2] :=
{
   abs[ (x2-x1)(y1-y0) - (x1-x0)(y2-y1)] / sqrt[(x2-x1)^2+(y2-y1)^2 ]
}


// Returns true if a point is inside a polygon.
// Adapted from
// http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
//
// points is an array of two-dimensional arrays of points.
// x,y is the point to test.
//
// Returns:
//   true if [x,y] is inside the polygon, false otherwise.
pointInPolygon[points, x, y] :=
{
   c = false
   i = 0
   npol = length[points]
   j = npol-1
   while i < npol
   {
      xi = points@i@0
      xj = points@j@0
      yi = points@i@1
      yj = points@j@1
      if ((((yi<=y) and (y<yj)) or ((yj<=y) and (y<yi))) and (x < (xj - xi) * (y - yi) / (yj - yi) + xi))
         c = !c
      j = i
      i = i + 1
   }

   return c;
}

// Calculates the area of a two-dimensional polygon.
//
// points is an array of [x,y] coordinates
//
// Returns: [cx, cy]
//
//  See:  http://en.wikipedia.org/wiki/Polygon#Area_and_centroid
//
//  The Frink polygon classes now implement area and centroid calculations
// using this code.
polygonArea[points] :=
{
   size = length[points]
   selfClosed = (points@0 == points@(size-1))

   if selfClosed
      n = size - 1
   else
      n = size

   if size == 0
      return 0

   if size <= 2
      return 0 * points@0@0 * points@0@1

   // Maintain units of measure.
   xs = 0 * points@0@0
   ys = 0 * points@0@1

   ps = xs ys
   xs = xs * ps
   ys = ys * ps
   
   for i = 0 to n-1
   {
      xi = points@i@0
      yi = points@i@1

      ii = i+1
      if ii == n and !selfClosed  // Make last point first point
         ii = 0
      
      xi1 = points@(ii)@0
      yi1 = points@(ii)@1
      pp = (xi yi1 - xi1 yi)
      ps = ps + pp
   }

   A = 1/2 ps
   return abs[A]
}

// Calculates the centroid of a two-dimensional polygon.
//
// points is an array of [x,y] coordinates
//
// Returns: [cx, cy]
//
//  See:  http://en.wikipedia.org/wiki/Polygon#Area_and_centroid
polygonCentroid[points] :=
{
   size = length[points]
   selfClosed = (points@0 == points@(size-1))

   if selfClosed
      n = size - 1
   else
      n = size

   if size == 0
      return undef

   if size == 1
      return points@0

   // Maintain units of measure.
   xs = 0 * points@0@0
   ys = 0 * points@0@1

   ps = xs ys
   xs = xs * ps
   ys = ys * ps
   
   for i = 0 to n-1
   {
      xi = points@i@0
      yi = points@i@1

      ii = i+1
      if ii == n and !selfClosed  // Make last point first point
         ii = 0
      
      xi1 = points@(ii)@0
      yi1 = points@(ii)@1

      pp = (xi yi1 - xi1 yi)
      xs = xs + (xi + xi1) pp
      ys = ys + (yi + yi1) pp
      ps = ps + pp
   }

   A = 1/2 ps
   cx = xs/(6 A)
   cy = ys/(6 A)

   return[cx, cy]
}

// Rotate the point [x,y] around the point [cx,cy] by the specified angle
// and return the coordinates of the new point, [x2, y2]
rotateAroundPoint[x,y, cx, cy, angle] :=
{
   ca = cos[angle]
   sa = sin[-angle]
   dx = x-cx
   dy = y-cy
   return [cx + ca dx - sa dy, cy + sa dx + ca dy]
}

// Finds the intersection[s] between a line and a circle given two points on
// the line [x1, y1], [x2, y2], the center of the circle [cx, cy],
// and the radius of the circle.
//
//  Returns:
//   an array, each element containing an array with [x,y] coordinates for
//  each intersection point.  There can 
//  be 0, 1, or 2 points of intersection.  For example, this may return
//  something like  [[1,1], [2,1]]
//
// See:
// http://mathworld.wolfram.com/Circle-LineIntersection.html
circleLineIntersections[x1, y1, x2, y2, cx, cy, r] :=
{
   // First shift coordinates so the center of the circle is at 0,0
   x1a = x1 - cx
   y1a = y1 - cy
   x2a = x2 - cx
   y2a = y2 - cy

   dx = x2 - x1
   dy = y2 - y1
   dr2 = dx^2 + dy^2
   D = x1a y2a - x2a y1a
   delta = r^2 dr2 - D^2

   // If delta < 0 there is no intersection.
   if delta < 0
      return new array

   sgndy = dy < 0 ? -1 : 1
   sqrtDelta = sqrt[delta]

   Ddy = D dy
   Ddx = - D dx
   
   // if delta == 0 there is only one point of intersecion
   if delta == 0
      return [[Ddy / dr2 + cx, Ddx / dr2 + cy]]
   
   rightx = sgndy dx sqrtDelta
   righty = abs[dy] sqrtDelta

   return [[(Ddy + rightx) / dr2 + cx, (Ddx + righty) / dr2 + cy],
           [(Ddy - rightx) / dr2 + cx, (Ddx - righty) / dr2 + cy]]
}


/** This performs an arbitrary affine transformation on the point [x,y]
    using the affine coordinates a,b,c,d,e,f.

    where the transformation matrix is defined by:

    x'    a c e     x
    y' =  b d f  *  y
    1     0 0 1     1

    (This is the order of the coordinates in the SVG and HTML5 standards.)

    which gives the equation:

    x' = a x + c y + e
    y' = b x + d y + f

    returns:  [x', y']
*/

affineTransform[x, y, a, b, c, d, e, f] := [a x + c y + e, b x + d y + f]


/** Performs an affineTransform with the coefficients specified as a list. */
affineTransform[x, y, coeffs] :=
{
   [a,b,c,d,e,f] = coeffs
   return affineTransform[x, y, a, b, c, d, e, f]
}


/** This performs the inverse of the arbitrary affine transformation above,
    converting the point [x',y'] back into the original coordinates [x,y]
    using the same affine transform coefficients a,b,c,d,e,f as passed to
    the affineTransform function above.

    returns:  [x, y]
*/

inverseAffineTransform[xp, yp, a, b, c, d, e, f] :=
{
   determinant = a d - b c
   if determinant == 0
      return [undef, undef]

   x =   (c (f - yp) - d (e - xp)) / determinant
   y =  -(a (f - yp) - b (e - xp)) / determinant
   
   return [x, y]
}

/* This performs the inverse of the arbitrary affine transformation above,
    converting the point [x',y'] back into the original coordinates [x,y]
    using the same affine transform coefficients coeffs=[a,b,c,d,e,f] as
    passed to the affineTransform function above.

    returns:  [x, y]
*/

inverseAffineTransform[xp, yp, coeffs] :=
{
   [a,b,c,d,e,f] = coeffs
   return inverseAffineTransform[xp, yp, a, b, c, d, e, f]
}

/** This composes a set of transformation coefficients with another set of
    transformation coefficients.  The inputs should be an array of 6 elements
    containing the coefficients [a,b,c,d,e,f] as described in the documentation
    for the affineTransformation function above.

    Returns a vector [a,b,c,d,e,f] of the
    new coefficients. */

composeTransformation[coeffsA, coeffsB] :=
{
   [a1, b1, c1, d1, e1, f1] = coeffsA
   [a2, b2, c2, d2, e2, f2] = coeffsB
   return composeTransformation[a1, b1, c1, d1, e1, f1, a2, b2, c2, d2, e2, f2]
}

/** This composes a set of transformation coefficients with another set of
    transformation coefficients.  Returns a vector [a,b,c,d,e,f] of the
    new coefficients.  This is just an unrolled, optimized matrix
    multiplication, eliminating always-zero terms. */

composeTransformation[a1, b1, c1, d1, e1, f1, a2, b2, c2, d2, e2, f2] :=
{
   return [a1 a2 + b2 c1,        // a
           a2 b1 + b2 d1,        // b
           a1 c2 + c1 d2,        // c
           b1 c2 + d1 d2,        // d
           a1 e2 + c1 f2 + e1,   // e
           b1 e2 + d1 f2 + f1]   // f
}

/** Makes a set of transformation coefficients that correspond to a translation
    by dx, dy.

    This corresponds to the matrix:

    x'    a=1  c=0  e=dx     x
    y' =  b=0  d=1  f=dy  *  y
    1       0    0    1      1

    which gives the equation:

    x' = x + dx
    y' = y + dy
*/

makeTranslate[dx, dy] :=
{
   return [1,   // a
           0,   // b
           0,   // c
           1,   // d
           dx,  // e
           dy]  // f
}


/** Create a scaling transformation that scales the coordinate system
    by sx, sy around 0,0.  This corresponds to the matrix:

    x'    a=sx  c=0   e=0     x
    y' =  b=0   d=sy  f=0  *  y
    1       0     0     1     1

    which gives the equation:

    x' = x sx
    y' = y sy
*/

makeScale[sx, sy] :=
{
   return [sx,  // a
           0,   // b
           0,   // c
           sy,  // d
           0,   // e
           0]   // f
}


/** Create a rotating transformation that rotates the coordinate system
    around 0,0 by the angle theta.  This corresponds to the matrix:

    x'    a=cos[theta]  c=-sin[theta]   e=0     x
    y' =  b=sin[theta]  d= cos[theta]   f=0  *  y
    1       0              0              1     1

    which gives the equations:

    x' = x cos[theta] - y sin[theta]
    y' = x sin[theta] + y cos[theta]
*/

makeRotate[theta] :=
{
   c = cos[theta]
   s = sin[theta]
   return [c,   // a
           s,   // b
           -s,  // c
           c,   // d
           0,   // e
           0]   // f
}

/** Create a rotating transformation that rotates the coordinate system
    around cx, cy by the angle theta.  This is equivalent to:

    translate[cx,cy]
    rotate[theta]
    translate[-cx, -cy]

    This corresponds to the matrix:

x'    a=cos[theta]  c=-sin[theta]  e=cx (1 - cos[theta]) + cy sin[theta]     x
y' =  b=sin[theta]  d= cos[theta]  f=cy (1 - cos[theta]) - cx sin[theta] *   y 
1       0              0             1                                       1

    which gives the equations:

    x' = (x - cx) cos[theta] - (y - cy) sin[theta] + cx
    y' = (x - cx) sin[theta] + (y - cy) cos[theta] + cy
*/

makeRotateAround[cx, cy, theta] :=
{
   transform = makeTranslate[cx, cy]
   transform = composeTransformation[transform, makeRotate[theta]]
   transform = composeTransformation[transform, makeTranslate[-cx, -cy]]
   return transform
}

/** Create a scaling transformation that scales the coordinate system
    around cx, cy by the scales sx, sy.  This is equivalent to:

    translate[cx, cy]
    scale[sx, sy]
    translate[-cx, -cy]

    This corresponds to the matrix:

    x'    a=sx   c=0    e=cx (1-sx)        x
    y' =  b=0    d=sy   f=cy (1-sy)    *   y 
    1       0      0      1                1

    which gives the equations:

    x' = sx (x-cx) + cx
    y' = sy (y-cy) + cy

   or

   x' = sx x + cx (1 - sx)
   y' = sy y + cy (1 - sy)
*/

makeScaleAround[cx, cy, sx, sy] :=
{
   transform = makeTranslate[cx, cy]
   transform = composeTransformation[transform, makeScale[sx, sy]]
   transform = composeTransformation[transform, makeTranslate[-cx, -cy]]
   return transform
}

/** Make an arbitrary transform with the specified coordinates. */
makeArbitrary[a, b, c, d, e, f] :=
{
   return [a, b, c, d, e, f]
}

/** Creates the inverse transform of a general transform.  For this to be
    valid, the determinant (a d - b c) must be nonzero. */

makeInverse[a, b, c, d, e, f] :=
{
   det = 1 / (a d - b c)
   return [ d det,
           -b det,
           -c det,
            a det,
           (c f - d e) det,
           (b e - a f) det]
}

/** Creates the inverse transform of a general transform.  For this to be
    valid, the determinant (a d - b c) must be nonzero. */

makeInverse[array] :=
{
   [a,b,c,d,e,f] = array
   makeInverse[a,b,c,d,e,f]
}

"Geometry.frink included OK"


View or download geometry.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 17646 days, 9 hours, 39 minutes ago.