August 2014

Sun Mon Tue Wed Thu Fri Sat
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            








« An interesting challenge: generating variable density fill patterns for 3D printing | Main | Generating hyperbolic tessellations inside AutoCAD using .NET »

January 23, 2012

Generating hyperbolic geometry on a Poincaré disk in AutoCAD using .NET

As mentioned in the last post, today’s post looks at how to draw hyperbolic geometry using the Poincaré disk model inside AutoCAD.

This is an interesting exercise, but probably won’t ultimately help us with the hyperbolic tessellations we’re aiming to create: it’s interesting as it will end up with us having a mechanism for mapping hyperbolic geometry onto an arbitrary circle inside AutoCAD, but as we’re probably going to end up creating straight-lined segments for the polygons we use to create our pattern, it’s mostly an exercise (and will probably be largely irrelevant to our end result). The tough part actually comes in the next post, when we move on to the tessellation algorithm.

To get started, I managed to find a straightforward C++ implementation of a set of functions to create hyperbolic geometry on a Poincaré disk. This library was developed by by David Coeurjolly, Research Director for the team focused on Multiresolution, Discrete and Combinatorial Models at LIRIS (Laboratoire d'InfoRmatique en Image et Systèmes d'information), part of the CNRS (Centre National de la Recherche Scientifique) in France. David’s implementation makes use of the Cairo graphics library, so there was a modest amount of work needed to generate geometry inside AutoCAD, instead. Indeed this was the main focus of my efforts: I fully admit I haven’t spent much time looking at the underlying mathematics (in particular, please don’t ask me to explain how the radius is calculated in the ComputeCircleParameters() method… ;-).

Here’s the resultant C# file containing the library class I’ve named HyperbolicEngine (for want of a better term):

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.Geometry;

using System;

 

namespace HyperbolicGeometry

{

  public class HyperbolicEngine

  {

    // We'll hold some AutoCAD objects to provide quick access

    // to the drawing database

 

    Transaction _tr;

    Database _db;

    BlockTableRecord _btr;

 

    // We'll also hold information about the target circle

 

    double _radius;

    Point3d _center;

    Matrix3d _disp;

 

    // Whether we'll draw curved or straight lines

 

    bool _straight;

 

    // Target layers for points, supporting lines and our

    // actual geometry

 

    ObjectId _ptLay, _lnLay, _hgLay;

 

    // Constructor

 

    public HyperbolicEngine(

      Transaction tr,

      Database db,

      Circle cir,

      bool straight,

      ObjectId? ptLay = null,

      ObjectId? lnLay = null,

      ObjectId? hgLay = null

    )

    {

      _tr = tr;

      _db = db;

 

      if (cir == null)

      {

        _radius = 400;

        _center = Point3d.Origin;

      }

      else

      {

        _center = cir.Center;

        _radius = cir.Radius;

      }

      _disp = Matrix3d.Displacement(_center - Point3d.Origin);

 

      _btr =

        (BlockTableRecord)tr.GetObject(

          db.CurrentSpaceId, OpenMode.ForWrite

        );

 

      _straight = straight;

 

      _ptLay = ptLay.HasValue ? ptLay.Value : db.Clayer;

      _lnLay = lnLay.HasValue ? lnLay.Value : db.Clayer;

      _hgLay = hgLay.HasValue ? hgLay.Value : db.Clayer;

    }

 

    // Functions to map a point from the unit disk to WCS

 

    private Point3d WcsPoint(double x, double y)

    {

      Point3d pt = new Point3d(cX(x), cY(y), 0);

      return pt.TransformBy(_disp);

    }

 

    private double cX(double x)

    {

      return x * _radius;

    }

 

    private double cY(double y)

    {

      return y * _radius;

    }

 

    /*

     * Compute the circle containing points a and b, and

     * orthogonal to the unit circle (a.k.a shortest path between

     * a and b in the hyperbolic sense).

     *

     * Return true if the output is a circle. False if the

     * shortest path is a straight line (diameter of the unit

     * circle).

     */

 

    private static bool ComputeCircleParameters(

      Point2d a, Point2d b,

      out double cx, out double cy, out double radius

    )

    {

      double ax = a.X, ay = a.Y, bx = b.X, by = b.Y;

 

      cx = 0; cy = 0; radius = 0;

 

      if (Math.Abs(ax * by - ay * bx) < Tolerance.Global.EqualPoint)

        return false;

 

      cx =

        (-1 / 2.0 *

          (ay * bx * bx + ay * by * by -

            (ax * ax + ay * ay + 1) * by + ay

          )

          / (ax * by - ay * bx)

        );

 

      cy =

        (1 / 2.0 *

          (ax * bx * bx + ax * by * by -

            (ax * ax + ay * ay + 1) * bx + ax

          )

          / (ax * by - ay * bx)

        );

 

      radius =

        -1 / 2.0 * Math.Sqrt(

          ax * ax * ax * ax * bx * bx +

          ax * ax * ax * ax * by * by -

          2 * ax * ax * ax * bx * bx * bx -

          2 * ax * ax * ax * bx * by * by +

          2 * ax * ax * ay * ay * bx * bx +

          2 * ax * ax * ay * ay * by * by -

          2 * ax * ax * ay * bx * bx * by -

          2 * ax * ax * ay * by * by * by +

          ax * ax * bx * bx * bx * bx +

          2 * ax * ax * bx * bx * by * by +

          ax * ax * by * by * by * by -

          2 * ax * ay * ay * bx * bx * bx -

          2 * ax * ay * ay * bx * by * by +

          ay * ay * ay * ay * bx * bx +

          ay * ay * ay * ay * by * by -

          2 * ay * ay * ay * bx * bx * by -

          2 * ay * ay * ay * by * by * by +

          ay * ay * bx * bx * bx * bx +

          2 * ay * ay * bx * bx * by * by +

          ay * ay * by * by * by * by -

          2 * ax * ax * ax * bx -

          2 * ax * ax * ay * by +

          4 * ax * ax * bx * bx -

          2 * ax * ay * ay * bx +

          8 * ax * ay * bx * by -

          2 * ax * bx * bx * bx -

          2 * ax * bx * by * by -

          2 * ay * ay * ay * by +

          4 * ay * ay * by * by -

          2 * ay * bx * bx * by -

          2 * ay * by * by * by +

          ax * ax - 2 * ax * bx +

          ay * ay - 2 * ay * by +

          bx * bx + by * by

        )

        / (ax * by - ay * bx);

 

      if (radius < 0.0)

        radius =

          1 / 2.0 * Math.Sqrt(

            ax * ax * ax * ax * bx * bx +

            ax * ax * ax * ax * by * by -

            2 * ax * ax * ax * bx * bx * bx -

            2 * ax * ax * ax * bx * by * by +

            2 * ax * ax * ay * ay * bx * bx +

            2 * ax * ax * ay * ay * by * by -

            2 * ax * ax * ay * bx * bx * by -

            2 * ax * ax * ay * by * by * by +

            ax * ax * bx * bx * bx * bx +

            2 * ax * ax * bx * bx * by * by +

            ax * ax * by * by * by * by -

            2 * ax * ay * ay * bx * bx * bx -

            2 * ax * ay * ay * bx * by * by +

            ay * ay * ay * ay * bx * bx +

            ay * ay * ay * ay * by * by -

            2 * ay * ay * ay * bx * bx * by -

            2 * ay * ay * ay * by * by * by +

            ay * ay * bx * bx * bx * bx +

            2 * ay * ay * bx * bx * by * by +

            ay * ay * by * by * by * by -

            2 * ax * ax * ax * bx -

            2 * ax * ax * ay * by +

            4 * ax * ax * bx * bx -

            2 * ax * ay * ay * bx +

            8 * ax * ay * bx * by -

            2 * ax * bx * bx * bx -

            2 * ax * bx * by * by -

            2 * ay * ay * ay * by +

            4 * ay * ay * by * by -

            2 * ay * bx * bx * by -

            2 * ay * by * by * by +

            ax * ax - 2 * ax * bx +

            ay * ay - 2 * ay * by +

            bx * bx + by * by

          )

          / (ax * by - ay * bx);

 

      return true;

    }

 

    // Draw the unit circle

 

    public void DrawUnitCircle()

    {

      if (_tr != null && _btr != null)

      {

        Circle c = new Circle(_center, Vector3d.ZAxis, _radius);

 

        c.LayerId = _hgLay;

 

        _btr.AppendEntity(c);

        _tr.AddNewlyCreatedDBObject(c, true);

      }

    }

 

    // Draw a point in the Poincare disc

 

    public void DrawPoint(Point2d a)

    {

      if (_tr != null && _btr != null)

      {

        DBPoint p = new DBPoint(WcsPoint(a.X, a.Y));

 

        p.LayerId = _ptLay;

 

        _btr.AppendEntity(p);

        _tr.AddNewlyCreatedDBObject(p, true);

      }

    }

 

    // Draw a hyperbolic line through a and b

 

    public void DrawSupportingLine(

      Point2d a, Point2d b, bool withPoints

    )

    {

      if (withPoints)

      {

        DrawPoint(a);

        DrawPoint(b);

      }

 

      double ax = a.X, ay = a.Y, bx = b.X, by = b.Y;

      double cx, cy, r;

 

      bool result =

        ComputeCircleParameters(a, b, out cx, out cy, out r);

 

      Entity ent;

 

      if (!result)

      {

        // We project a and b points onto the unit circle

 

        double theta = Math.Atan2(ay, ax);

        double theta2 = Math.Atan2(by, bx);

 

        ent =

          new Line(

            WcsPoint(Math.Cos(theta), Math.Sin(theta)),

            WcsPoint(Math.Cos(theta2), Math.Sin(theta2))

          );

      }

      else

      {

        ent =

          new Circle(WcsPoint(cx, cy), Vector3d.ZAxis, r * _radius);

      }

 

      if (ent != null)

      {

        ent.LayerId = _lnLay;

 

        if (_tr != null && _btr != null)

        {

          _btr.AppendEntity(ent);

          _tr.AddNewlyCreatedDBObject(ent, true);

        }

        else

        {

          ent.Dispose();

        }

      }

    }

 

    // Draw a hyperbolic edge between a and b

 

    public void DrawEdge(

      Point2d a, Point2d b, bool withPoints,

      bool straightSegments = false

    )

    {

      if (withPoints)

      {

        DrawPoint(a);

        DrawPoint(b);

      }

 

      double ax = a.X, ay = a.Y, bx = b.X, by = b.Y;

      double cx = 0.0, cy = 0.0, r = 0.0;

      bool result =

        _straight ? false :

          ComputeCircleParameters(a, b, out cx, out cy, out r);

 

      // Near-aligned points -> straight segment

 

      Entity ent = null;

 

      if (!result)

      {

        ent = new Line(WcsPoint(ax, ay), WcsPoint(bx, by));

      }

      else

      {

        double theta = -Math.Atan2(-ay + cy, ax - cx);

        double theta2 = -Math.Atan2(-by + cy, bx - cx);

 

        // We recenter the angles to [0, 2Pi]

 

        while (theta < 0)

        {

          theta += 2 * Math.PI;

        }

        while (theta2 < 0)

        {

          theta2 += 2 * Math.PI;

        }

 

        Point3d cen = WcsPoint(cx, cy);

 

        double delta = Math.Abs(theta - theta2);

 

        if (

          (theta > theta2 && delta > Math.PI) ||

          (theta < theta2 && delta < Math.PI)

        )

        {

          ent = new Arc(cen, r * _radius, theta, theta2);

        }

        else

        {

          ent = new Arc(cen, r * _radius, theta2, theta);

        }

      }

      if (ent != null)

      {

        ent.LayerId = _hgLay;

 

        if (_tr != null && _btr != null)

        {

          _btr.AppendEntity(ent);

          _tr.AddNewlyCreatedDBObject(ent, true);

        }

        else

        {

          ent.Dispose();

        }

      }

    }

 

    // Draw a hyperbolic triangle with vertices (a,b,c)

 

    public void DrawTriangle(

      Point2d a, Point2d b, Point2d c,

      bool withPoints = false, bool withSupportLines = false

    )

    {

      DrawPolygon(

        new Point2d[] { a, b, c }, withPoints, withSupportLines

      );

    }

 

    // Draw a hyperbolic polygon with the provided vertices

 

    public void DrawPolygon(

      Point2d[] pts, bool withPoints = false,

      bool withSupportLines = false

    )

    {

      if (withPoints)

      {

        foreach (Point2d pt in pts)

        {

          DrawPoint(pt);

        }

      }

 

      int numPts = pts.Length;

 

      if (withSupportLines && !_straight)

      {

        for (int i = 0; i < numPts; i++)

        {

          DrawSupportingLine(

            pts[i], pts[(i + 1) % numPts], withPoints

          );

        }

      }

 

      for (int i = 0; i < numPts; i++)

      {

        DrawEdge(

          pts[i], pts[(i + 1) % numPts], withPoints

        );

      }

    }

  }

}

Here’s the command class implementation:

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.Colors;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.EditorInput;

using Autodesk.AutoCAD.Geometry;

using Autodesk.AutoCAD.Runtime;

using System;

 

namespace HyperbolicGeometry

{

  public class Commands

  {

    [CommandMethod("HGT")]

    public static void CreateHyperbolicTriangles()

    {

      Document doc =

        Application.DocumentManager.MdiActiveDocument;

      Database db = doc.Database;

 

      Transaction tr =

        doc.TransactionManager.StartTransaction();     

      using (tr)

      {

        LayerTable lt =

          (LayerTable)tr.GetObject(

            db.LayerTableId, OpenMode.ForWrite

          );

 

        // Points and supporting lines will be on layers that

        // are grey and initially frozen

 

        ObjectId ptLay = CreateLayer(tr, lt, "HG_PTS", 9, false);

        ObjectId lnLay = CreateLayer(tr, lt, "HG_LNS", 9, false);

        ObjectId htLay = CreateLayer(tr, lt, "HG_CNT", 6, true);

 

        // Create our hyperbolic engine

 

        HyperbolicEngine he =

          new HyperbolicEngine(

            tr, db, null, false, ptLay, lnLay, htLay

          );

 

        he.DrawUnitCircle();

 

        he.DrawTriangle(

          new Point2d(0.8, -0.3),

          new Point2d(0.2, -0.6),

          new Point2d(-0.61, -0.6),

          true, true

        );

 

        tr.Commit();

      }

    }

 

    [CommandMethod("HGS")]

    public static void CreateHyperbolicSpiro()

    {

      Document doc =

        Application.DocumentManager.MdiActiveDocument;

      Database db = doc.Database;

 

      Transaction tr =

        doc.TransactionManager.StartTransaction();

      using (tr)

      {

        LayerTable lt =

          (LayerTable)tr.GetObject(

            db.LayerTableId, OpenMode.ForWrite

          );

 

        // Points and supporting lines will be on layers that

        // are grey and initially frozen

 

        ObjectId ptLay = CreateLayer(tr, lt, "HG_PTS", 9, false);

        ObjectId lnLay = CreateLayer(tr, lt, "HG_LNS", 9, false);

        ObjectId htLay = CreateLayer(tr, lt, "HG_CNT", 6, true);

 

        // Create our hyperbolic engine

 

        HyperbolicEngine he =

          new HyperbolicEngine(

            tr, db, null, false, ptLay, lnLay, htLay

          );

 

        he.DrawUnitCircle();

 

        int max = 13;

        for (int k = 0; k < (max - 1); k++)

        {

          he.DrawTriangle(

            new Point2d(1, 0),

            new Point2d(

              Math.Cos(k * 134 / (double)max),

              Math.Sin(k * 134 / (double)max)

            ),

            new Point2d(

              Math.Cos((k + 1) * 134 / (double)max),

              Math.Sin((k + 1) * 134 / (double)max)

            ),

            true, true

          );

        }

 

        tr.Commit();

      }

    }

 

    private static ObjectId CreateLayer(

      Transaction tr, LayerTable lt,

      string lname, short col, bool on

    )

    {

      if (lt.Has(lname))

        return lt[lname];

 

      LayerTableRecord ltr = new LayerTableRecord();

      ltr.Color = Color.FromColorIndex(ColorMethod.ByAci, col);

      ltr.Name = lname;

      ltr.IsFrozen = !on;

 

      ObjectId layId = lt.Add(ltr);

      tr.AddNewlyCreatedDBObject(ltr, true);

 

      return layId;

    }

  }

}

The code implements two commands, both of which have been lifted fairly directly from David’s test cases.

  • HGT – creates a simple hyperbolic triangle
  • HGS – creates a more complex set of hyperbolic triangles

Here are the results of running the HGT command:

A simple hyperbolic triangle

And if you turn on the layer showing the supporting lines (HG_LNS), you’ll see the circles used to generated these arcs:

And with its supporting lines displayed

The HGS command creates more interesting geometry:

A more complex  set of triangles

Which, of course, also has more interesting supporting lines along with it:

Some of their supporting lines

You can certainly see similarities between geometry on a Poincaré disk and the hypocycloids generated by Spiro.

Next time we’ll look at hyperbolic tessellations and how to generate them inside AutoCAD.

blog comments powered by Disqus

Feed/Share

10 Random Posts