view release on metacpan or  search on metacpan

lib/Math/PlanePath.pm  view on Meta::CPAN

equilateral triangle.  That sqrt(3) factor can be applied if desired,

    X, Y*sqrt(3)          side length 2

    X/2, Y*sqrt(3)/2      side length 1

Integer Y values have the advantage of fitting pixels on the usual kind of
raster computer screen, and not losing precision in floating point results.

If doing a general-purpose coordinate rotation then be sure to apply the
sqrt(3) scale factor before rotating or the result will be skewed.  60
degree rotations can be made within the integer X,Y coordinates directly as
follows, all giving integer X,Y results.

    ( X-3Y)/2, ( X+Y)/2     rotate +60   (anti-clockwise)
    ( X+3Y)/2, (-X+Y)/2     rotate -60   (clockwise)
    (-X-3Y)/2, ( X-Y)/2     rotate +120
    (-X+3Y)/2, (-X-Y)/2     rotate -120
    -X,-Y                   rotate 180

    (X+3Y)/2, (X-Y)/2       mirror across the X=3*Y twelfth (30deg)

=cut

# GP-DEFINE  sqrt3i = quadgen(-12);
# GP-Test  sqrt3i^2 == -3
# GP-DEFINE  w6 =  1/2 + 1/2*sqrt3i;
# GP-DEFINE  w3 = -1/2 + 1/2*sqrt3i;
# GP-DEFINE  b12 = 3/2 + 1/2*sqrt3i;

# rot +/-60
# GP-Test   ('x/2 + 'y*sqrt3i/2)*w6 \
# GP-Test   == (('x - 3*'y)/2)*1/2  + (( 'x + 'y)/2)*sqrt3i/2

# GP-Test   ('x/2 + 'y*sqrt3i/2)/w6 \
# GP-Test   == (('x + 3*'y)/2)*1/2  + ((-'x + 'y)/2)*sqrt3i/2

# rot +/-120
# GP-Test   ('x/2 + 'y*sqrt3i/2)*w3 \
# GP-Test   == ((-'x - 3*'y)/2)*1/2  + (( 'x - 'y)/2)*sqrt3i/2

# GP-Test   ('x/2 + 'y*sqrt3i/2)/w3 \
# GP-Test   == ((-'x + 3*'y)/2)*1/2  + ((-'x - 'y)/2)*sqrt3i/2

# mirror across 30deg
# GP-Test   conj(('x/2 + 'y*sqrt3i/2)/b12)*b12 \
# GP-Test   == (('x + 3*'y)/2)*1/2  + (('x - 'y)/2)*sqrt3i/2

=pod

The sqrt(3) factor can be worked into a hypotenuse radial distance
calculation as follows if comparing distances from the origin.

    hypot = sqrt(X*X + 3*Y*Y)

See for instance C<TriangularHypot> which is triangular points ordered by
this radial distance.

=head1 FORMULAS

The formulas section in the POD of each class describes some of the
calculations.  This might be of interest even if the code is not.

=head2 Triangular Calculations

For a triangular lattice, the rotation formulas above allow calculations to
be done in the rectangular X,Y coordinates which are the inputs and outputs
of the PlanePath functions.  Another way is to number vertically on a 60
degree angle with coordinates i,j,

          ...
          *   *   *      2
        *   *   *       1
      *   *   *      j=0
    i=0  1   2

These coordinates are sometimes used for hexagonal grids in board games etc.
Using this internally can simplify rotations a little,

    -j, i+j         rotate +60   (anti-clockwise)
    i+j, -i         rotate -60   (clockwise)
    -i-j, i         rotate +120
    j, -i-j         rotate -120
    -i, -j          rotate 180

Conversions between i,j and the rectangular X,Y are

    X = 2*i + j         i = (X-Y)/2
    Y = j               j = Y

A third coordinate k at a +120 degrees angle can be used too,

     k=0  k=1 k=2
        *   *   *
          *   *   *
            *   *   *
             0   1   2

This is redundant in that it doesn't number anything i,j alone can't
already, but it has the advantage of turning rotations into just sign
changes and swaps,

    -k, i, j        rotate +60
    j, k, -i        rotate -60
    -j, -k, i       rotate +120
    k, -i, -j       rotate -120
    -i, -j, -k      rotate 180

The conversions between i,j,k and the rectangular X,Y are like the i,j above
but with k worked in too.

    X = 2i + j - k        i = (X-Y)/2        i = (X+Y)/2
    Y = j + k             j = Y         or   j = 0
                          k = 0              k = Y

=head2 N to dX,dY -- Fractional

C<n_to_dxdy()> is the change from N to N+1, and is designed both for integer
N and fractional N.  For fractional N it can be convenient to calculate a
dX,dY at floor(N) and at floor(N)+1 and then combine the two in proportion
to frac(N).

                     int+2
                      |
                      |
                      N+1    \
                     /|       |
                    / |       |
                   /  |       | frac
                  /   |       |
                 /    |       |
                /     |      /
       int-----N------int+1
    this_dX  dX,dY     next_dX
    this_dY            next_dY

       |-------|------|
         frac   1-frac


    int = int(N)
    frac = N - int    0 <= frac < 1

    this_dX,this_dY  at int
    next_dX,next_dY  at int+1

    at fractional N
      dX = this_dX * (1-frac) + next_dX * frac
      dY = this_dY * (1-frac) + next_dY * frac

This is combination of this_dX,this_dY and next_dX,next_dY in proportion to
the distances from positions N to int+1 and from int+1 to N+1.

The formulas can be rearranged to

    dX = this_dX + frac*(next_dX - this_dX)
    dY = this_dY + frac*(next_dY - this_dY)

which is like dX,dY at the integer position plus fractional part of a turn
or change to the next dX,dY.

=head2 N to dX,dY -- Self-Similar

For most of the self-similar paths such as C<HilbertCurve>, the change dX,dY
is determined by following the state table transitions down through either
all digits of N, or to the last non-9 digit, ie. drop any low digits equal
to radix-1.

Generally paths which are the edges of some tiling use all digits, and those
which are the centres of a tiling stop at the lowest non-9.  This can be
seen for example in the C<DekkingCurve> using all digits, whereas its
C<DekkingCentres> variant stops at the lowest non-24.

Perhaps this all-digits vs low-non-9 would even characterize path style as
edges or centres of a tiling, when a path is specified in some way that a
tiling is not quite obvious.

=head1 SUBCLASSING

The mandatory methods for a PlanePath subclass are

    n_to_xy()
    xy_to_n()
    xy_to_n_list()     if multiple N's map to an X,Y
    rect_to_n_range()

It sometimes happens that one of C<n_to_xy()> or C<xy_to_n()> is easier than
the other but both should be implemented.

C<n_to_xy()> should do something sensible on fractional N.  The suggestion
is to make it an X,Y proportionally between integer N positions.  It could
be along a straight line or an arc as best suits the path.  A straight line
can be done simply by two calculations of the surrounding integer points,
until it's clear how to work the fraction into the code directly.

C<xy_to_n_list()> has a base implementation calling plain C<xy_to_n()> to
give a single N at X,Y.  If a path has multiple Ns at an X,Y
(eg. C<DragonCurve>) then it must implement C<xy_to_n_list()> to return all
those Ns, and must also implement a plain C<xy_to_n()> returning the first
of them.

C<rect_to_n_range()> can initially be any convenient over-estimate.  It
should give N big enough that from there onwards all points are sure to be
beyond the given X,Y rectangle.

The following descriptive methods have base implementations

    n_start()           1
    class_x_negative()  \ 1, so whole plane
    class_y_negative()  /
    x_negative()        calls class_x_negative()
    y_negative()        calls class_x_negative()
    x_negative_at_n()   undef \ as for no negatives
    y_negative_at_n()   undef /