› ›

inscribing ellipses in arbitrary polygons?

terekitaterekita
in Programming Questions

Hello, I'm looking to inscribe ellipses in arbitrary polygons, along the lines of http://faculty.mae.carleton.ca/John_Hayes/Papers/InscribingEllipse.pdf , but not constrained to quadrilaterals. Kind of a long shot, but has anyone happened to code up this particular concept? cheers...

Answers

  • HorseofCaligulusHorseofCaligulus
    edited November 2015

    I'll give it a try, just to clarify though, you want the largest area of an ellipse that can fit in a given polygon?

  • terekitaterekita

    thanks! Yes, the largest area of an ellipse that can fit in a given polygon is exactly what I'm looking for. I would be super excited to see what you might come up with—this problem takes significantly more computational geometry than I have. At the same time, it feels tantalizingly solvable...

  • terekitaterekita

    also, this would be for convex shapes only...

  • HorseofCaligulusHorseofCaligulus

    haha, that makes it much easier, I had made an entire pathfinding program to find my way through screwy polygons, haven't had much time recently but I should be able to finish it this weekend... maybe, but anyway basically i run a function several times in different areas around the shape, and the function expands an ellipse, tilts moves and shifts it until the amount of precision needed is past anything useful. then it returns the coordinates of the ellipse, I'll post the code as soon as it's done

  • terekitaterekita

    oh man, that sounds awesome! sorry I wasn't more precise in the initial description!

  • quarkquark

    i run a function several times in different areas around the shape, and the function expands an ellipse,

    The article shows how to calculate the centre of the ellipse, see fig 2 and the first paragraph of part 4

  • quarkquark

    It has taken me about 4-5 hours research and programming but EUREKA the sketch below will calculate the largest inscribed ellipse to fit inside a convex quadrilateral.

    Here is the output

    image

    Here is the code

     /**
 * This sketch was inspired buy a question on the Processing forum
 * https://forum.processing.org/two/discussion/13399/inscribing-ellipses-in-arbitrary-polygons#latest
 * 
 * REFERENCES
 * http://faculty.mae.carleton.ca/John_Hayes/Papers/InscribingEllipse.pdf
 *    How to calculate the largest ellipse to fit inside a convex quadrilateral
 *    
 * http://chrisjones.id.au/Ellipses/ellipse.html
 *    How to calculate the 'projective collineation' matrix mention in the
 *    first reference and from that the coefficients of the general equation of an ellipse.
 *    ax2 + 2bxy + cy2 + 2dx + 2fy + g = 0
 *    
 * http://mathworld.wolfram.com/Ellipse.html
 *    How to calculate the centre, semi-axis lengths and angle of rotation of the ellipse
 *    from the coefficients a,b,c,d,f,g  
 * 
 * @author Peter Lager
 *
 */

Quadrilateral q;
Ellipse e;

public void settings() {
  size(400, 400);
}

public void setup() {
  e = new Ellipse();
  q = new Quadrilateral(149, 43, 53, 271, 340, 333, 301, 199);
  calcLargestInscribedEllipse(q, e);
  background(220, 250, 220);
  q.render();
  e.render();
}

// 
public void calcLargestInscribedEllipse(Quadrilateral q, Ellipse e) {
  PVector r = q.v[0];
  PVector s = q.v[1];
  PVector t = q.v[2];
  PVector u = q.v[3];
  // Create the 'projective collineation' matrix
  Matrix3x3F mat =  new Matrix3x3F(
    s.x * t.x * u.y - r.x * t.x * u.y - s.x * t.y * u.x + r.x * t.y * u.x - r.x * s.y * u.x + r.y * s.x * u.x + r.x * s.y * t.x - r.y * s.x * t.x, 
    r.x * t.x * u.y - r.x * s.x * u.y - s.x * t.y * u.x + s.y * t.x * u.x - r.y * t.x * u.x + r.y * s.x * u.x + r.x * s.x * t.y - r.x * s.y * t.x, 
    s.x * t.x * u.y - r.x * s.x * u.y - r.x * t.y * u.x - s.y * t.x * u.x + r.y * t.x * u.x + r.x * s.y * u.x + r.x * s.x * t.y - r.y * s.x * t.x, 
    s.y * t.x * u.y - r.y * t.x * u.y - r.x * s.y * u.y + r.y * s.x * u.y - s.y * t.y * u.x + r.y * t.y * u.x + r.x * s.y * t.y - r.y * s.x * t.y, 
    -s.x * t.y * u.y + r.x * t.y * u.y + s.y * t.x * u.y - r.x * s.y * u.y - r.y * t.y * u.x + r.y * s.y * u.x + r.y * s.x * t.y - r.y * s.y * t.x, 
    s.x * t.y * u.y - r.x * t.y * u.y + r.y * t.x * u.y - r.y * s.x * u.y - s.y * t.y * u.x + r.y * s.y * u.x + r.x * s.y * t.y - r.y * s.y * t.x, 
    s.x * u.y - r.x * u.y - s.y * u.x + r.y * u.x - s.x * t.y + r.x * t.y + s.y * t.x - r.y * t.x, 
    t.x * u.y - s.x * u.y - t.y * u.x + s.y * u.x + r.x * t.y - r.y * t.x - r.x * s.y + r.y * s.x, 
    t.x * u.y - r.x * u.y - t.y * u.x + r.y * u.x + s.x * t.y - s.y * t.x + r.x * s.y - r.y * s.x
    );

  // Get the inverse of the 'projective collineation' matrix
  Matrix3x3F matI = mat.inverse();
  // Calculate the coefficients of the ellipse equation
  // ax2 + 2bxy + cy2 + 2dx + 2fy + g = 0
  float a = matI.a11 * matI.a11 + matI.a21 * matI.a21 - matI.a31 * matI.a31;
  float b = matI.a11 * matI.a12 + matI.a21 * matI.a22 - matI.a31 * matI.a32;
  float c = matI.a12 * matI.a12 + matI.a22 * matI.a22 - matI.a32 * matI.a32;
  float d = matI.a11 * matI.a13 + matI.a21 * matI.a23 - matI.a31 * matI.a33;
  float f = matI.a12 * matI.a13 + matI.a22 * matI.a23 - matI.a32 * matI.a33;
  float g = matI.a13 * matI.a13 + matI.a23 * matI.a23 - matI.a33 * matI.a33;
  // Calculate the centre of the ellipse
  float b2ac = b*b - a*c;
  e.cx = (c*d - b*f) / b2ac;
  e.cy = (a*f - b*d) / b2ac;
  // Calculate the semi-axis lengths of the ellipse
  float k = 2*(a*f*f + c*d*d + g*b*b - 2*b*d*f - a*c*g) / b2ac;
  float ac24b2 = sqrt( (a - c)*(a - c) + 4*b*b );
  e.sa = sqrt(abs(k / abs( ac24b2 - a - c)));
  e.sb = sqrt(abs(k / abs(-ac24b2 - a - c)));
  // Calculate the rotation angle of the ellipse
  if (b == 0) {
    e.ang = (a < c) ? 0 : HALF_PI;
  } else {
    float acot = atan(2*b / (a-c));
    e.ang = (a < c) ? 0.5f * acot : HALF_PI + 0.5f * acot;
  }
}

/**
 * Simple class to hold the coordinates of a convex quadrilateral.
 */
public class Quadrilateral {
  public PVector[] v = new PVector[4];

  /**
   * The coordinates of the quads vertices expressed in anti-clockwise order
   */
  public Quadrilateral(float x0, float y0, float x1, float y1, float x2, float y2, float x3, float y3) {
    v[0] = new PVector(x0, y0, 1);
    v[1] = new PVector(x1, y1, 1);
    v[2] = new PVector(x2, y2, 1);
    v[3] = new PVector(x3, y3, 1);
  }  

  public void render() {
    stroke(128);
    strokeWeight(2);
    for (int i = 0; i < 4; i++) {
      int i2 = (i + 1) % 4;
      line(v[i].x, v[i].y, v[i2].x, v[i2].y);
    }
  }
}

/**
 * Simple class to store the details of an ellipse
 * @author peter
 *
 */
public class Ellipse {
  public float cx, cy;
  public float sa, sb;
  public float ang;

  public void render() {
    pushMatrix();
    translate(cx, cy);
    rotate(ang);
    stroke(200, 20, 200);
    strokeWeight(1);
    fill(255, 230, 255);
    ellipse(0, 0, 2*sa, 2*sb);
    popMatrix();
  }
}


/**
 * A minimal implementation of a 3x3 matrix for this sketch
 * 
 * This is a modified version of the Matrix3 class found at 
 * http://www.java2s.com/Code/Java/2D-Graphics-GUI/A3x3matriximplementation.htm
 * and released under the GPL license, available at
 * http://www.gnu.org/copyleft/gpl.html.
 */
public class Matrix3x3F {
  public float a11, a12, a13;
  public float a21, a22, a23;
  public float a31, a32, a33;

  public Matrix3x3F() {
    a11=0; 
    a12=0; 
    a13=0;
    a21=0; 
    a22=0; 
    a23=0;
    a31=0; 
    a32=0; 
    a33=0;
  }

  /**
   * Construct a 3x3 matrix using specified fields
   */
  public Matrix3x3F(float a11, float a12, float a13, float a21, float a22, 
    float a23, float a31, float a32, float a33) {  
    this.a11 = a11;
    this.a12 = a12;
    this.a13 = a13;
    this.a21 = a21;
    this.a22 = a22;
    this.a23 = a23;
    this.a31 = a31;
    this.a32 = a32;
    this.a33 = a33;
  }

  /**
   * Compute the inverse of the matrix A, place the result in C
   */
  private Matrix3x3F inverse( final Matrix3x3F A, final Matrix3x3F C ) {
    float d = (A.a31*A.a12*A.a23-A.a31*A.a13*A.a22-A.a21*A.a12*A.a33+A.a21*A.a13*A.a32+A.a11*A.a22*A.a33-A.a11*A.a23*A.a32);
    float t11 =  (A.a22*A.a33-A.a23*A.a32)/d;
    float t12 = -(A.a12*A.a33-A.a13*A.a32)/d;
    float t13 =  (A.a12*A.a23-A.a13*A.a22)/d;
    float t21 = -(-A.a31*A.a23+A.a21*A.a33)/d;
    float t22 =  (-A.a31*A.a13+A.a11*A.a33)/d;
    float t23 = -(-A.a21*A.a13+A.a11*A.a23)/d;
    float t31 =  (-A.a31*A.a22+A.a21*A.a32)/d;
    float t32 = -(-A.a31*A.a12+A.a11*A.a32)/d;
    float t33 =  (-A.a21*A.a12+A.a11*A.a22)/d;

    C.a11 = t11; 
    C.a12 = t12; 
    C.a13 = t13;
    C.a21 = t21; 
    C.a22 = t22; 
    C.a23 = t23;
    C.a31 = t31; 
    C.a32 = t32; 
    C.a33 = t33;
    return C;
  }

  public final Matrix3x3F inverse() {
    return inverse(this, new Matrix3x3F());
  }
}
  • terekitaterekita

    wow, this is awesome! If we ever meet in the real world, I'll buy you lunch, or drinks, your choice. Hell, I'll buy you both! This code will absolutely be extremely useful to me. That said, and feeling sheepish for pointing this out, what I am ultimately looking for is inscribing the ellipse within arbitrary convex polygons, not only quadrilaterals. My original post is probably confusing in that though I mention arbitrary polygons, I linked to a paper specifically about quads. Again, thanks so much for taking all that time to code this up!

  • quarkquark
    edited November 2015

    The paper stated it requires a numerical approach or polygons with >=5 edges which would be much slower and I am not sure where to start, more research probably :)

    Don't worry about the time I spent on this I enjoyed the challenge and I was really pleased with the outcome.

  • terekitaterekita

    again, I really appreciate it, and this will still be extremely helpful :)

  • quarkquark

    OK I did some research of polygons with >= 5 sides. Most info was about inscribing polygons inside ellipses and I found nothing I could use.

    Perhaps someone else can help :)

  • terekitaterekita

    ok, thanks so much for looking into it, much appreciated!

Sign In or Register to comment.