Hello Jim.
Thanks for taking the time to think about this.

I implemented a few different offsetting schemes. On well behaved
curves, my original "parallel first and last control vectors and go
through B(0.5)" algorithm was by far the best. Theoretically it is
also somewhat better justified than the others (some of which were
like Apache Harmony's way - offsetting p2 and p3), since we're 
interpolating instead of just using some heuristic. It is also
the only one I've encountered that widens quarter circles optimally,
and it only needs one curve per side too (i.e. if we have to widen
the path of a PathIterator of a circle with radius r, using width w,
Pisces' output will be identical to the output of the PathIterators
of 2 circles with radii r+w and r-w (perhaps not identical identical,
since PathIterators can use doubles, and we only use floats in pisces,
but... that's nitpicking)).

As I've said before, the only 2 problems with it were that it was bad
on high curvatures (but this was fixed by breaking down the curves into
monotonic ones), and it was bad on curves like 
p.moveTo(0,0);p.curveTo(100,100,0,100,100,1). I thought the latter was
fixable with the "d1/(d1+d2)" algorithm I suggested in my previous e-mail.
I implemented it, and it turns out I was wrong. Then I came up with my own
variation of that algorithm (the original one I used was from Don Lancaster's
website) that did sort of work. But then I implemented your suggestion of
splitting the curve at the inflection points as well as breaking it into
monotonic pieces, and the problem was gone. I didn't even have to use the
fancy variation of the "d1/(d1 + d2)" algorithm - just the plain old
interpolation one worked (I should say that the fancy one is still probably
better, but undetectably so, now that we break at inflection points and at
xy direction changes, so the added complexity is probably not worth it).

I've attached my latest Stroker.java, if you want to take a look at these
algorithms (they're in computeOffsetWay1 (for the old interpolation) and
computeOffsetWay3 (for the fancy version). There are 2 more, but these
aren't as good, and one is just shameful). I didn't make a webrev because
I still need to fix how it handles cusps (which should be easy), and I
need to implement a way to avoid subdividing a rotated quarter circle.
I can do this either by rotating curves so that p2-p1 is parallel to the
x axis and then subdivide like I do now (i.e. make monotonic, break at
inflections) or get rid of the monotonic subdivision, and instead subdivide
by checking the angles of the control polygon. I could make it so it
subdivides whenever the angles between p2-p1 and p3-p2 or p3-p2 and p4-p3
are greater than 45 degrees. I think this would actually ensure that
every resulting curve can be rotated to be made monotonic in both x and
y (at least after inflections are eliminated), which means it's at least 
as good as what I have now.

Very well. I've convinced myself. I'll implement the latter.

Thanks,
Denis.

----- "Jim Graham" <james.gra...@oracle.com> wrote:

> Hi Denis,
> 
> Just to let you know that I've seen this and I'm thinking.
> 
> But it will take a bit of "more thinking" when I have time for more. 
> Hopefully in a couple of days.
> 
> For one thing, it sounds like you already understand the Apache
> Harmony 
> approach quite a bit better than I ever did and, in particular, why it
> 
> didn't work well - which is encouraging.  Hopefully your solution will
> 
> sound pretty good when I get around to wrapping my head around it...
> 
>                       ...jim
> 
> On 8/30/2010 3:03 PM, Denis Lila wrote:
> > Hello Jim.
> >
> >> One way in which they may not break enough is that I think that
> >> inflections also need to be broken in order to find a parallel
> curve
> >> (though I suppose a very tiny inflection might still be
> approximated by
> >> a parallel curve easily) and inflections can happen at any angle
> without
> >> going horizontal or vertical.
> >
> >      It wouldn't be hard to add additional breaks at the inflection
> points.
> >
> >> Secondly, although a circle tends to be represented by quadrant
> sections
> >> which are monotonic, a circle rotated by 45 degrees would have
> >> horizontal and vertical sections in the middle of each quadrant and
> you
> >> would split those, but really they can be made parallel regardless
> of
> >> angle so these would be unnecessary splits.
> >
> >      That's true, and it's one reason I didn't like my method very
> much when I sent
> > my previous e-mail. However, what if we rotated curves so that B'(0)
> is
> > parallel to the x-axis before splitting at points where dx/dt or
> dy/dt are 0?
> > This would certainly get rid of this problem for circles, and
> probably for
> > other curves. All it would cost is 1 Math.cos and 1 Math.sin and a
> few
> > multiplications and additions per curve.
> >
> > The biggest problem with my implementation was that some curves that
> were almost
> > lines were drawn horribly. I computed the parallel (p1', p2', p3',
> p4') of
> > a given curve (p1, p2, p3, p4) by making p2'-p1' parallel to p2-p1,
> > making p4'-p3' parallel to p4-p3. This leaves a 2 unknowns, so to
> get 2 more
> > equations I made the computed parallel at t=0.5 equal to the ideal
> parallel
> > at t=0.5. The problem was that for some curves (ones with very high
> curvature)
> > p2'-p1' and p4'-p3' were parallel to p2-p1 and p4-p3 but their
> directions were
> > opposite, so you can imagine what the computed "parallel" looked
> like.
> > However, I am almost certain that the problem was caused by making
> C(0.5) = P(0.5)
> > (where C is the computed curve, and P is the ideal parallel). A much
> better
> > restriction, that I think would eliminate the problem would be
> C(d1/(d1 + d2)) = P(0.5).
> > where d1 = ||P(0.5)-P(0)|| and d2 = ||P(1)-P(0.5)||.
> >
> >> My belief is that lengths and angles of the control polygon help
> >> determine if it is well-behaved and can be made parallel simply by
> >> offsetting.  Some formula involving those values would likely be
> happy
> >> with circle sections regardless of their angle of rotation.  I
> believe
> >> the Apache Harmony version of BasicStroke used those criteria...
> >
> >      I've been reading the Apache Harmony file. The way they do it
> is compute
> > the tangent of an angle using the line width and a predefined
> constant
> > (which doesn't seem to be related to the curve's polygon - it's just
> a decreasing
> > function with the range (-1,1)), and if some angles in the polygon
> are smaller than
> > the computed angle, offset curves are computed. Otherwise the curve
> is subdivided.
> > However, I can't understand how offsets for the control points are
> computed (i.e.
> > why they use the equations they use, and why they work).
> > If we're using the widening of quarter circles as a yard stick, then
> Apache Harmony's
> > BasicStroke does not do very well. If we have a quarter circle from
> (1,0) to (0,1),
> > then the control points c1 and c2 should be (1,k) and (k,1) where k
> = 4*(sqrt(2)-1)/3.
> > A parallel curve w away from this quarter circle should have control
> points
> > (1+w,k+k*w) and (k+k*w,1+w). I traced Apache Harmony's computations
> on a quarter
> > circle, and this is not what it outputs. Furthermore, all quarter
> circles with line
> > width<  9.65 will be subdivided. My method with the modifications
> suggested above
> > doesn't have these problems.
> >
> >      But perhaps it's better to not interpolate through P(0.5) after
> all.
> > In addition to making p4'-p3' and p2'-p1' parallel to p4-p3 and
> p2-p1, we could
> > also make p3'-p2' parallel to p3-p2. These constraints leave just
> one unknown, which
> > needs to be eliminated. I am not sure how to do this. I thought
> about making the line
> > (p2',p3') be a distance of w from (p2,p3) (which is exactly what
> needs to be done for
> > (p1',p2') and (p3',p4')) but this doesn't get us what we want for
> quarter circles. So, finding
> > the control points would boil down to finding intersections of 3
> lines.
> >
> > Do you have any suggestions on how to do the offsetting for the
> control points?
> >
> > Thank you,
> > Denis.
> >
> > ----- "Jim Graham"<james.gra...@oracle.com>  wrote:
> >
> >> Hi Denis,
> >>
> >> At the bottom-most rendering level monotonic curves can be cool to
> deal
> >> with, but I'm dubious that they help with widening.  For one
> things, I
> >> think you need more breaks than they would give you and also they
> might
> >> sometimes break a curve when it doesn't need it.
> >>
> >>            ...jim
> >>
> >> On 8/25/2010 2:36 PM, Denis Lila wrote:
> >>> Hello Jim.
> >>>
> >>>> I think a more dynamic approach that looked at how "regular" the
> >> curve
> >>>> was would do better.  Regular is hard to define, but for
> instance
> >> a
> >>>> bezier section of a circle could have parallel curves computed
> >> very
> >>>> easily without having to flatten or subdivide further.  Curves
> >> with
> >>>> inflections probably require subdividing to get an accurate
> >> parallel
> >>>> curve.
> >>>
> >>> I'm not sure if you read it, but after the email with the webrev
> >> link
> >>> I sent another describing a different method of widening: split
> the
> >>> curve at every t where dx/dt == 0 and dy/dt == 0. This guarantees
> >> that
> >>> there will be no more than 5 curves per side, and since each
> curve
> >> will
> >>> be monotonic in both x and y the curve that interpolates its
> >> parallel
> >>> should do a pretty good job.
> >>>
> >>> This is far better than interpolating at regular t intervals, but
> >> I'm
> >>> trying to find a better way. I don't like this because the split
> >>> depend not only on the curve itself, but also on the axes. The
> axes
> >> are
> >>> arbitrary, so this is not good. For example a curve like this
> >>> |
> >>> \_ would get widened by 1 curve per side (which is good and
> >> optimal), but
> >>> if we rotate this curve by, say, 30 degrees it would be widened by
> 2
> >> curves
> >>> per side.
> >>> It also doesn't handle cusps and locations of high curvature very
> >> well (although
> >>> I think the latter is a numerical issue that could be made better
> by
> >> using
> >>> doubles).
> >>>
> >>> Regards,
> >>> Denis.
> >>>
> >>> ----- "Jim Graham"<james.gra...@oracle.com>   wrote:
> >>>
> >>>> Hi Denis,
> >>>>
> >>>> On 8/23/2010 4:18 PM, Denis Lila wrote:
> >>>>>        To widen cubic curves, I use a cubic spline with a fixed
> >> number
> >>>> of curves for
> >>>>> each curve to be widened. This was meant to be temporary, until
> I
> >>>> could find a
> >>>>> better algorithm for determining the number of curves in the
> >> spline,
> >>>> but I
> >>>>> discovered today that that won't do it.
> >>>>>        For example, the curve p.moveTo(0,0),p.curveTo(84.0,
> 62.0,
> >>>> 32.0, 34.0, 28.0, 5.0)
> >>>>> looks bad all the way up to ~200 curves. Obviously, this is
> >>>> unacceptable.
> >>>>>
> >>>>> It would be great if anyone has any better ideas for how to go
> >> about
> >>>> this.
> >>>>> To me it seems like the problem is that in the webrev I chop up
> >> the
> >>>> curve to be
> >>>>> interpolated at equal intervals of the parameter.
> >>>
> >>>>
> >>>> Perhaps looking at the rate of change of the slope (2nd and/or
> 3rd
> >>>> derivatives) would help to figure out when a given section of
> >> curve
> >>>> can
> >>>> be approximated with a parallel version?
> >>>>
> >>>> I believe that the BasicStroke class of Apache Harmony returns
> >> widened
> >>>>
> >>>> curves, but when I tested it it produced a lot more curves than
> >> Ductus
> >>>>
> >>>> (still, probably a lot less than the lines that would be
> produced
> >> by
> >>>> flattening) and it had some numerical problems.  In the end I
> >> decided
> >>>> to
> >>>> leave our Ductus stuff in there until someone could come up with
> a
> >>>> more
> >>>> reliable Open Source replacement, but hopefully that code is
> close
> >>>> enough to be massaged into working order.
> >>>>
> >>>> You can search the internet for "parallel curves" and find lots
> of
> >>>> literature on the subject.  I briefly looked through the web
> >> sites,
> >>>> but
> >>>> didn't have enough time to remember enough calculus and
> >> trigonometry
> >>>> to
> >>>> garner a solution out of it all with the time that I had.  Maybe
> >>>> you'll
> >>>> have better luck following the algorithms... ;-)
> >>>>
> >>>> As far as flattening at the lowest level when doing scanline
> >>>> conversion,
> >>>> I like the idea of using forward differencing as it can create
> an
> >>>> algorithm that doesn't require all of the intermediate storage
> that
> >> a
> >>>>
> >>>> subdividing flattener requires.  One reference that describes
> the
> >>>> technique is on Dr. Dobbs site, though I imagine there are many
> >>>> others:
> >>>>
> >>>>
> >>
> http://www.drdobbs.com/184403417;jsessionid=O5N5MDJRDMIKHQE1GHOSKH4ATMY32JVN
> >>>>
> >>>> You can also look at the code in
> >>>> src/share/native/sun/java2d/pipe/ProcessPath.c for some examples
> >> of
> >>>> forward differencing in use (and that code also has similar
> >> techniques
> >>>>
> >>>> to what you did to first chop the path into vertical pieces). 
> BUT
> >>>> (*Nota Bene*), I must warn you that the geometry of the path is
> >>>> somewhat
> >>>> perturbed in that code since it tries to combine "path
> >> normalization"
> >>>>
> >>>> and rasterization into a single process.  It's not rendering
> pure
> >>>> geometry, it is rendering tweaked geometry which tries to make
> >> non-AA
> >>>>
> >>>> circles look round and other such aesthetics-targeted
> impurities.
> >>>> While
> >>>> the code does have good examples of how to set up and evaluate
> >> forward
> >>>>
> >>>> differencing equations, don't copy too many of the details or
> you
> >>>> might
> >>>> end up copying some of the tweaks and the results will look
> >> strange
> >>>> under AA.  The Dr. Dobbs article should be your numerical
> >> reference
> >>>> and
> >>>> that reference code a practical, but incompatible, example...
> >>>>
> >>>>                  ...jim
/*
 * Copyright (c) 2007, Oracle and/or its affiliates. All rights reserved.
 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
 *
 * This code is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License version 2 only, as
 * published by the Free Software Foundation.  Oracle designates this
 * particular file as subject to the "Classpath" exception as provided
 * by Oracle in the LICENSE file that accompanied this code.
 *
 * This code is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * version 2 for more details (a copy is included in the LICENSE file that
 * accompanied this code).
 *
 * You should have received a copy of the GNU General Public License version
 * 2 along with this work; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
 *
 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
 * or visit www.oracle.com if you need additional information or have any
 * questions.
 */

package pisces;

import java.util.Arrays;

import sun.awt.geom.PathConsumer2D;

public class Stroker implements PathConsumer2D {

    private static final int MOVE_TO = 0;
    private static final int LINE_TO = 1;
    private static final int CLOSE = 2;
    private static final int CURVE_TO = 3;
    private static final int QUAD_TO = 4;

    /**
     * Constant value for join style.
     */
    public static final int JOIN_MITER = 0;

    /**
     * Constant value for join style.
     */
    public static final int JOIN_ROUND = 1;

    /**
     * Constant value for join style.
     */
    public static final int JOIN_BEVEL = 2;

    /**
     * Constant value for end cap style.
     */
    public static final int CAP_BUTT = 0;

    /**
     * Constant value for end cap style.
     */
    public static final int CAP_ROUND = 1;

    /**
     * Constant value for end cap style.
     */
    public static final int CAP_SQUARE = 2;

    private final PathConsumer2D output;

    private final int capStyle;
    private final int joinStyle;

    private final float m00, m01, m10, m11, det;

    private final float lineWidth2;
    private final float scaledLineWidth2;

    private final float[][] offset = new float[3][2];
    private final PolyStack reverse = new PolyStack();
    private final float[] miter = new float[2];
    private final float miterLimitSq;

    private int prev;
    private boolean started;

    private float sx0, sy0, sx1, sy1, x0, y0, px0, py0;
    private float mx0, my0, omx, omy;

    private final float m00_2_m01_2;
    private final float m10_2_m11_2;
    private final float m00_m10_m01_m11;

    /**
     * Constructs a <code>Stroker</code>.
     *
     * @param output an output <code>LineSink</code>.
     * @param lineWidth the desired line width in pixels
     * @param capStyle the desired end cap style, one of
     * <code>CAP_BUTT</code>, <code>CAP_ROUND</code> or
     * <code>CAP_SQUARE</code>.
     * @param joinStyle the desired line join style, one of
     * <code>JOIN_MITER</code>, <code>JOIN_ROUND</code> or
     * <code>JOIN_BEVEL</code>.
     * @param miterLimit the desired miter limit
     * @param transform a <code>Transform4</code> object indicating
     * the transform that has been previously applied to all incoming
     * coordinates.  This is required in order to produce consistently
     * shaped end caps and joins.
     */
    public Stroker(PathConsumer2D output,
                   float lineWidth,
                   int capStyle,
                   int joinStyle,
                   float miterLimit,
                   float m00, float m01, float m10, float m11) {
        this.output = output;

        this.lineWidth2 = lineWidth / 2;
        this.scaledLineWidth2 = m00 * lineWidth2;
        this.capStyle = capStyle;
        this.joinStyle = joinStyle;

        this.m00_2_m01_2 = m00*m00 + m01*m01;
        this.m10_2_m11_2 = m10*m10 + m11*m11;
        this.m00_m10_m01_m11 = m00*m10 + m01*m11;

        this.m00 = m00;
        this.m01 = m01;
        this.m10 = m10;
        this.m11 = m11;
        this.det = m00*m11 - m01*m10;

        float limit = miterLimit * lineWidth2 * det;
        this.miterLimitSq = limit*limit;

        prev = CLOSE;
        started = false;
    }

    private void computeOffset(float lx, float ly, float[] m) {
        float dx, dy;
        if (m00 > 0 && m00 == m11 && m01 == 0 & m10 == 0) {
            float ilen = (float)Math.hypot(lx, ly);
            if (ilen == 0) {
                dx = dy = 0;
            } else {
                dx = (ly * scaledLineWidth2)/ilen;
                dy = -(lx * scaledLineWidth2)/ilen;
            }
        } else {
            int sdet = (det > 0) ? 1 : -1;
            float a = ly * m00 - lx * m10;
            float b = ly * m01 - lx * m11;
            float dh = (float)Math.hypot(a, b);
            float div = sdet * lineWidth2/dh;

            float ddx = ly * m00_2_m01_2 - lx * m00_m10_m01_m11;
            float ddy = ly * m00_m10_m01_m11 - lx * m10_2_m11_2;
            dx = ddx*div;
            dy = ddy*div;
        }

        m[0] = dx;
        m[1] = dy;
    }

    private boolean isCCW(float x0, float y0,
                          float x1, float y1,
                          float x2, float y2) {
        return (x1 - x0) * (y2 - y1) < (y1 - y0) * (x2 - x1);
    }

    private boolean side(float x,  float y,
                         float x0, float y0,
                         float x1, float y1) {
        return (y0 - y1)*x + (x1 - x0)*y + (x0*y1 - x1*y0) > 0;
    }

    // pisces used to use fixed point arithmetic with 16 decimal digits. I
    // didn't want to change the values of the constants below when I converted
    // it to floating point, so that's why the divisions by 2^16 are there.
    private static final float ROUND_JOIN_THRESHOLD = 1000/65536f;

    private void drawRoundJoin(float x, float y,
                               float omx, float omy, float mx, float my,
                               int side,
                               boolean flip,
                               boolean rev,
                               float threshold) {
        if ((omx == 0 && omy == 0) || (mx == 0 && my == 0)) {
            return;
        }

        float domx = omx - mx;
        float domy = omy - my;
        float len = domx*domx + domy*domy;
        if (len < threshold) {
            return;
        }

        if (rev) {
            omx = -omx;
            omy = -omy;
            mx = -mx;
            my = -my;
        }

        float bx0 = x + omx;
        float by0 = y + omy;
        float bx1 = x + mx;
        float by1 = y + my;

        drawBezRoundJoin(x, y, bx0, by0, bx1, by1, side, flip, rev);
    }

    private void drawBezRoundJoin(float cx, float cy,
    		float ax, float ay,
    		float bx, float by,
    		int oppositeSide,
    		boolean flip,
    		boolean rev) {
    	// the side of the line on which the arc will be drawn. true means left
    	boolean side; 
    	if (oppositeSide == 0) {
    		side = !side(cx, cy, ax, ay, bx, by);
    	} else {
    		side = (oppositeSide != 1);
    	}
        float[] bezJoin = new float[8];
    	bezJoin[0] = cx; bezJoin[1] = cy;
    	bezJoin[2] = ax; bezJoin[3] = ay;
    	bezJoin[4] = bx; bezJoin[5] = by;

    	getUntransformedPoints(bezJoin, m00, m01, m10, m11);

    	float pcx = bezJoin[0], pcy = bezJoin[1];
    	float x0 = bezJoin[2] - pcx, y0 = bezJoin[3] - pcy;
    	float x3 = bezJoin[4] - pcx, y3 = bezJoin[5] - pcy;

    	double a = Math.atan2((double)y0, (double)x0);
    	double t2 = Math.atan2((double)y3, (double)x3);
    	double twoPi = Math.PI * 2;
    	// We need the to be angles between 0 and 2pi.
    	if (a < 0) {a += twoPi;}
    	if (t2 < 0) {t2 = t2 + twoPi;}

    	// ext is the angle that added to a gets us to the angle
    	// of the second end point of the arc being drawn. Its sign 
    	// indicates the direction in which the arc will be drawn.
    	double ext = (t2 > a) ? (t2 - a - twoPi) : (t2 - a);
    	if (!side) { // drawing on the right side
    		ext += twoPi;
    	}

    	// If the arc is more than 90 degrees we approximate it using 2
    	// bezier curves, to get higher accuracy.
    	// Math.max is there in case ceil returns a double slightly smaller than
    	// 1, which would get cast to 0.
    	int numCurves = Math.max((int)Math.ceil((2 * Math.abs(ext)) / Math.PI), 1);
    	ext /= numCurves;
    	// cv is the length of P1-P0 and P2-P3 divided by the radius of the arc
    	// (so, cv assumes the arc has radius 1). P0, P1, P2, P3 are the points that
    	// define the bezier curve we're computing.
    	// It is computed using the constraints that P1-P0 and P3-P2 are parallel
    	// to the arc tangents at a and a+t respectively, and that |P1-P0|=|P3-P2|.
    	double cv = 4.0 / 3.0 * Math.sin(ext/2) / (1.0 + Math.cos(ext/2));

    	// I am not sure why this line and the last are necessary, but without
    	// them certain joins are not drawn properly. In particular, polylines.
    	// eg: g2.drawPolyLine(new int[]{300,300,800}, new int[]{50,500,50}, 3);
    	emitLineTo(ax, ay, rev);
    	double relx = Math.cos(a);
    	double rely = Math.sin(a);
    	for (int i = 0; i < numCurves; i++) {
    		bezJoin[0] = (float) (pcx + relx * lineWidth2);
    		bezJoin[1] = (float) (pcy + rely * lineWidth2);
    		bezJoin[2] = (float) (pcx + (relx - cv*rely) * lineWidth2);
    		bezJoin[3] = (float) (pcy + (rely + cv*relx) * lineWidth2);
    		a += ext;
    		relx = Math.cos(a);
    		rely = Math.sin(a);
    		bezJoin[4] = (float)(pcx + (relx + cv*rely) * lineWidth2);
    		bezJoin[5] = (float)(pcy + (rely - cv*relx) * lineWidth2);
    		bezJoin[6] = (float)(pcx + relx * lineWidth2);
    		bezJoin[7] = (float)(pcy + rely * lineWidth2);

    		getTransformedPoints(bezJoin, m00, m01, m10, m11);
    		if (!rev) {
    			emitCurveTo(bezJoin[2], bezJoin[3],
    					    bezJoin[4], bezJoin[5],
    					    bezJoin[6], bezJoin[7]);
    		} else {
    			reverse.pushCubic(bezJoin[0], bezJoin[1], bezJoin[2], bezJoin[3], bezJoin[4], bezJoin[5]);
    		}
    	}
    	emitLineTo(bx, by, rev);
    }

    private static void getUntransformedPoints(float[] xy, float m00, float m01, float m10, float m11) {
    	float det = (m00 * m11 - m01 * m10);
    	for (int i = 0; i < xy.length; i += 2) {
    		float x = (m11 * xy[i] - m01 * xy[i + 1]) / det;
    		float y = (m00 * xy[i + 1] - m10 * xy[i]) / det;
    		xy[i] = x;
    		xy[i + 1] = y;
    	}
    }

    private static void getTransformedPoints(float[] xy, float m00, float m01, float m10, float m11) {
    	for (int i = 0; i < xy.length; i += 2) {
    	float x = (m00 * xy[i] + m01 * xy[i + 1]);
    	float y = (m10 * xy[i] + m11 * xy[i + 1]);
        xy[i] = x;
        xy[i + 1] = y;
    	}
    }


    // Return the intersection point of the lines (ix0, iy0) -> (ix1, iy1)
    // and (ix0p, iy0p) -> (ix1p, iy1p) in m[0] and m[1]
    private void computeMiter(float x0, float y0, float x1, float y1,
                              float x0p, float y0p, float x1p, float y1p,
                              float[] m) {
        float x10 = x1 - x0;
        float y10 = y1 - y0;
        float x10p = x1p - x0p;
        float y10p = y1p - y0p;

        float den = x10*y10p - x10p*y10;
        if (den == 0) {
            m[0] = x0;
            m[1] = y0;
            return;
        }

        float t = x1p*(y0 - y0p) - x0*y10p + x0p*(y1p - y0);
        m[0] = x0 + (t*x10)/den;
        m[1] = y0 + (t*y10)/den;
    }

    private void drawMiter(float px0, float py0,
                           float x0, float y0,
                           float x1, float y1,
                           float omx, float omy, float mx, float my,
                           boolean rev) {
        if (mx == omx && my == omy) {
            return;
        }
        if (px0 == x0 && py0 == y0) {
            return;
        }
        if (x0 == x1 && y0 == y1) {
            return;
        }

        if (rev) {
            omx = -omx;
            omy = -omy;
            mx = -mx;
            my = -my;
        }

        computeMiter(px0 + omx, py0 + omy, x0 + omx, y0 + omy,
                     x0 + mx, y0 + my, x1 + mx, y1 + my,
                     miter);

        // Compute miter length in untransformed coordinates
        float dx = miter[0] - x0;
        float dy = miter[1] - y0;
        float a = dy*m00 - dx*m10;
        float b = dy*m01 - dx*m11;
        float lenSq = a*a + b*b;

        if (lenSq < miterLimitSq) {
            emitLineTo(miter[0], miter[1], rev);
        }
    }

    public void moveTo(float x0, float y0) {
        if (prev == LINE_TO) {
            finish();
        }

        this.sx0 = this.x0 = x0;
        this.sy0 = this.y0 = y0;
        this.started = false;
        this.prev = MOVE_TO;
    }

    public void lineTo(float x1, float y1) {
    	if (x1 == x0 && y1 == y0) {
            return;
        }
        lineToImpl(x1, y1);
    }

    private void lineToImpl(float x1, float y1) {
        computeOffset(x1 - x0, y1 - y0, offset[0]);
        float mx = offset[0][0];
        float my = offset[0][1];

        drawJoin(px0, py0, x0, y0, x1, y1, omx, omy, mx, my);

        emitLineTo(x0 + mx, y0 + my, false);
        emitLineTo(x1 + mx, y1 + my, false);

        emitLineTo(x0 - mx, y0 - my, true);
        // suppose the next command is a curveTo. We still need to put the 
        // following lineTo on the stack because if that curve is counter
        // clockwise, then a lineTo(x0, y0) will be put on the stack in
        // the function handling that curve (the x0,y0 in that function are
        // our x1, y1). And so, when we're going through the reverse stack,
        // we will emit a line to our x1,y1, and then we will emit a line
        // to our x0+mx,y0+mx, and the line will be a trapezoid instead of
        // a rectangle. Therefore, the following line is necessary. If that
        // curve is clockwise, or if this is the last line, then we will
        // emit a 0-length line, which isn't a problem because it'll
        // probably just be ignored by our consumer.
        emitLineTo(x1 - mx, y1 - my, true);

        this.omx = mx;
        this.omy = my;
        this.px0 = x0;
        this.py0 = y0;
        this.x0 = x1;
        this.y0 = y1;
        this.prev = LINE_TO;
    }

    public void closePath() {
        if (!started) {
            finish();
            return;
        }

        if (x0 != sx0 || y0 != sy0) {
        	lineToImpl(sx0, sy0);
        }

        drawJoin(px0, py0, x0, y0, sx1, sy1, omx, omy, mx0, my0);

        emitLineTo(x0 + mx0, y0 + my0, false);
        
        emitMoveTo(sx0 - mx0, sy0 - my0);
        emitReverse();

        this.started = false;
        this.prev = CLOSE;
        emitClose();
    }

    private void emitReverse() {
    	float[] buf = new float[8];
        while(!reverse.isEmpty()) {
        	int type = reverse.pop(buf);
            switch (type) {
            // TODO: all these can be put in a function that takes the type
            // of curve as an argument, and a float array with the points.
            // then we don't have to have a separate function for each 
            // curve type. But, of course, we'd have to convert everything to
            // arrays. urgh :(
            case LINE_TO:
                emitLineTo(buf[0], buf[1]);
                break;
            case QUAD_TO:
            	emitQuadTo(buf[0], buf[1], buf[2], buf[3]);
            	break;
            case CURVE_TO:
                emitCurveTo(buf[0], buf[1],
                		    buf[2], buf[3],
                		    buf[4], buf[5]);
                break;
            }
        }
    }

    public void pathDone() {
        if (prev == LINE_TO) {
            finish();
        }

        output.pathDone();
        this.prev = MOVE_TO; // TODO: probably remove
    }

    double userSpaceLineLength(double dx, double dy) {
        double a = (dy*m00 - dx*m10)/det;
        double b = (dy*m01 - dx*m11)/det;
        return Math.hypot(a, b);
    }

    private void finish() {
        if (capStyle == CAP_ROUND) {
            drawRoundJoin(x0, y0,
                          omx, omy, -omx, -omy, 1, false, false,
                          ROUND_JOIN_THRESHOLD);
        } else if (capStyle == CAP_SQUARE) {
            float dx = px0 - x0;
            float dy = py0 - y0;
            float len = (float)userSpaceLineLength(dx, dy);
            float s = lineWidth2/len;

            float capx = x0 - dx*s;
            float capy = y0 - dy*s;

            emitLineTo(capx + omx, capy + omy);
            emitLineTo(capx - omx, capy - omy);
        }

        emitReverse();

        if (capStyle == CAP_ROUND) {
            drawRoundJoin(sx0, sy0,
                          -mx0, -my0, mx0, my0, 1, false, false,
                          ROUND_JOIN_THRESHOLD);
        } else if (capStyle == CAP_SQUARE) {
            float dx = sx1 - sx0;
            float dy = sy1 - sy0;
            float len = (float)userSpaceLineLength(dx, dy);
            float s = lineWidth2/len;

            float capx = sx0 - dx*s;
            float capy = sy0 - dy*s;

            emitLineTo(capx - mx0, capy - my0);
            emitLineTo(capx + mx0, capy + my0);
        }

        emitClose();
    }

    private void emitMoveTo(float x0, float y0) {
        output.moveTo(x0, y0);
    }

    private void emitLineTo(float x1, float y1) {
        output.lineTo(x1, y1);
    }
    
    private void emitLineTo(float x1, float y1, boolean rev) {
    	if (rev) { 
    		reverse.pushLine(x1, y1);
    	} else {
    		emitLineTo(x1, y1);
    	}
    }
    
    private void emitQuadTo(float x1, float y1, 
    		                float x2, float y2) {
        output.quadTo(x1, y1, x2, y2);
    }
    
    private void emitCurveTo(float x1, float y1, 
                             float x2, float y2, 
                             float x3, float y3) {
//    	System.out.printf("curveTo: (%f, %f), (%f, %f), (%f, %f)%n", x1, y1, x2, y2, x3, y3);
    	output.curveTo(x1, y1, x2, y2, x3, y3);
    }

    private void emitClose() {
        output.closePath();
    }
   
    private void drawJoin(float px0, float py0,
    		              float x0, float y0,
    		              float x1, float y1,
    		              float omx, float omy,
    		              float mx, float my) {
    	if (!started) {
            emitMoveTo(x0 + mx, y0 + my);
            // TODO: in the flattening version sx1, sy1 was the first (x1,y1)
            // Now that's not true, but this is still correct because it was
            // only used for slopes, and for that purpose it's still good.
            // It should be changed however, because it's confusing.
            this.sx1 = x1;
            this.sy1 = y1;
            this.mx0 = mx;
            this.my0 = my;
            started = true;
        } else {
        	// TODO: px0, py0 have a similar problem to sx1, sy1. Now they should
        	// be set to the second control point of the last curveTo.
            boolean ccw = isCCW(px0, py0, x0, y0, x1, y1);
            if (joinStyle == JOIN_MITER) {
            	drawMiter(px0, py0, x0, y0, x1, y1, omx, omy, mx, my,
            			ccw);
            } else if (joinStyle == JOIN_ROUND) {
            	drawRoundJoin(x0, y0,
            			omx, omy,
            			mx, my, 0, false, ccw,
            			ROUND_JOIN_THRESHOLD);
            }
            emitLineTo(x0, y0, !ccw);
        }
    }

    // compute offset curves using bezier spline through t=0.5 (i.e. 
    // ComputedCurve(0.5) == IdealParallelCurve(0.5))
    private void computeOffsetWay1(float[] pts, float[] right, float[] left) {
    	float x1 = pts[0], y1 = pts[1];
    	float x2 = pts[2], y2 = pts[3];
    	float x3 = pts[4], y3 = pts[5];
    	float x4 = pts[6], y4 = pts[7];
    	BezCurve cur = BezCurve.makeCheapBezCurve(null, pts, 8);
    	
        float x = cur.xat(0.5f);
        float y = cur.yat(0.5f);
        float dx4 = x4 - x3;
        float dy4 = y4 - y3;
        float dx1 = x2 - x1;
        float dy1 = y2 - y1;
        float dxm = cur.dxat(0.5f);
        float dym = cur.dyat(0.5f);

        // this computes the offsets at t=0, 0.5, 1, using the property that
        // for any bezier curve the vectors p2-p1 and p4-p3 are parallel to
        // the (dx/dt, dy/dt) vectors at the endpoints.
        computeOffset(dx1, dy1, offset[0]);
        computeOffset(dxm, dym, offset[1]);
        computeOffset(dx4, dy4, offset[2]);
        float x1p = x1 + offset[0][0]; // start
        float y1p = y1 + offset[0][1]; // point
        float xi  = x + offset[1][0]; // interpolation
        float yi  = y + offset[1][1]; // point
        float x4p = x4 + offset[2][0]; // end
        float y4p = y4 + offset[2][1]; // point

        float det43 = 4f / (3f * (dx1 * dy4 - dy1 * dx4));
        
        float two_pi_m_p1_m_p4x = 2*xi - x1p - x4p;
        float two_pi_m_p1_m_p4y = 2*yi - y1p - y4p;
        float c1 = det43 * (dy4 * two_pi_m_p1_m_p4x - dx4 * two_pi_m_p1_m_p4y);
        float c2 = det43 * (dx1 * two_pi_m_p1_m_p4y - dy1 * two_pi_m_p1_m_p4x);
        
        float x2p, y2p, x3p, y3p;
        x2p = x1p + c1*dx1;
        y2p = y1p + c1*dy1;
        x3p = x4p + c2*dx4;
        y3p = y4p + c2*dy4;

        right[0] = x1p; right[1] = y1p;
        right[2] = x2p; right[3] = y2p;
        right[4] = x3p; right[5] = y3p;
        right[6] = x4p; right[7] = y4p;
        
        x1p = x1 - offset[0][0]; y1p = y1 - offset[0][1];
        xi = xi - 2 * offset[1][0]; yi = yi - 2 * offset[1][1];
        x4p = x4 - offset[2][0]; y4p = y4 - offset[2][1];

        two_pi_m_p1_m_p4x = 2*xi - x1p - x4p;
        two_pi_m_p1_m_p4y = 2*yi - y1p - y4p;
        c1 = det43 * (dy4 * two_pi_m_p1_m_p4x - dx4 * two_pi_m_p1_m_p4y);
        c2 = det43 * (dx1 * two_pi_m_p1_m_p4y - dy1 * two_pi_m_p1_m_p4x);
        
        x2p = x1p + c1*dx1;
        y2p = y1p + c1*dy1;
        x3p = x4p + c2*dx4;
        y3p = y4p + c2*dy4;
        
        left[0] = x1p; left[1] = y1p;
        left[2] = x2p; left[3] = y2p;
        left[4] = x3p; left[5] = y3p;
        left[6] = x4p; left[7] = y4p;
        System.out.println("right: " + Arrays.toString(right));
        System.out.println("left: " + Arrays.toString(left));
    }
    
    private void computeOffsetWay2(float[] pts, float[] right, float[] left) {
    	float x1 = pts[0], y1 = pts[1];
    	float x2 = pts[2], y2 = pts[3];
    	float x3 = pts[4], y3 = pts[5];
    	float x4 = pts[6], y4 = pts[7];
    	BezCurve cur = BezCurve.makeCheapBezCurve(null, pts, 8);

        final float dx4 = x4 - x3;
        final float dy4 = y4 - y3;
        float dx1 = x2 - x1;
        float dy1 = y2 - y1;

        // this computes the offsets at t=0, 0.5, 1, using the property that
        // for any bezier curve the vectors p2-p1 and p4-p3 are parallel to
        // the (dx/dt, dy/dt) vectors at the endpoints.
        computeOffset(dx1, dy1, offset[0]);
        computeOffset(dx4, dy4, offset[2]);
        float x1p = x1 + offset[0][0]; // start
        float y1p = y1 + offset[0][1]; // point
        float x4p = x4 + offset[2][0]; // end
        float y4p = y4 + offset[2][1]; // point

        float x2p, y2p, x3p, y3p;
        computeOffset(cur.dxat(1f/3f), cur.dyat(1f/3f), offset[1]);
        float mx2 = offset[1][0]; float my2 = offset[1][1];
        x2p = x2 + mx2;    y2p = y2 + my2;
        computeOffset(cur.dxat(2f/3f), cur.dyat(2f/3f), offset[1]);
        float mx3 = offset[1][0]; float my3 = offset[1][1];
        x3p = x3 + mx3;    y3p = y3 + my3;

        right[0] = x1p; right[1] = y1p;
        right[2] = x2p; right[3] = y2p;
        right[4] = x3p; right[5] = y3p;
        right[6] = x4p; right[7] = y4p;
        
        x1p = x1 - offset[0][0]; y1p = y1 - offset[0][1];
        x2p = x2 - mx2;    y2p = y2 - my2;
        x3p = x3 - mx3;    y3p = y3 - my3;
        x4p = x4 - offset[2][0]; y4p = y4 - offset[2][1];
        
        left[0] = x1p; left[1] = y1p;
        left[2] = x2p; left[3] = y2p;
        left[4] = x3p; left[5] = y3p;
        left[6] = x4p; left[7] = y4p;
    }
    
    // exactly like way 1, but instead of making ComputedCurve(0.5) == IdealParallel(0.5)
    // we do ComputedCurve(d1/(d1 + d2)) = IdealParallel(0.5).
    private void computeOffsetWay3(float[] pts, float[] right, float[] left) {
    	float x1 = pts[0], y1 = pts[1];
    	float x2 = pts[2], y2 = pts[3];
    	float x3 = pts[4], y3 = pts[5];
    	float x4 = pts[6], y4 = pts[7];
    	BezCurve cur = BezCurve.makeCheapBezCurve(null, pts, 8);
    	
        float x = cur.xat(0.5f);
        float y = cur.yat(0.5f);
        float dx4 = x4 - x3;
        float dy4 = y4 - y3;
        float dx1 = x2 - x1;
        float dy1 = y2 - y1;
        float dxm = cur.dxat(0.5f);
        float dym = cur.dyat(0.5f);

        // this computes the offsets at t=0, 0.5, 1, using the property that
        // for any bezier curve the vectors p2-p1 and p4-p3 are parallel to
        // the (dx/dt, dy/dt) vectors at the endpoints.
        computeOffset(dx1, dy1, offset[0]);
        computeOffset(dxm, dym, offset[1]);
        computeOffset(dx4, dy4, offset[2]);
        float x1p = x1 + offset[0][0]; // start
        float y1p = y1 + offset[0][1]; // point
        float xi  = x + offset[1][0]; // interpolation
        float yi  = y + offset[1][1]; // point
        float x4p = x4 + offset[2][0]; // end
        float y4p = y4 + offset[2][1]; // point

    	double d1 = Math.hypot(x - x1, y - y1);
    	double d2 = Math.hypot(x4 - x, y4 - y);
        
        float[] curarr = right;
        for (int i = 0; i < 2; i++) {
        	double d1p = Math.hypot(xi - x1p, yi - y1p);
        	double d2p = Math.hypot(x4p - xi, y4p - yi);
        	float t = (float) ((d1p * (d1 + d2)) / (2 * d1 * (d1p + d2p)));
        	System.out.println(t);
        	float t1 = 1 - t;
        	float c1 = t1 * t1 * t1;
        	float c2 = 3 * t1 * t1 * t;
        	float c3 = 3 * t1 * t * t;
        	float c4 = t * t * t;
        	float A00 = c2 * (x2 - x1);
        	float A01 = c3 * (x3 - x4);
        	float A10 = c2 * (y2 - y1);
        	float A11 = c3 * (y3 - y4);
        	float det = A00 * A11 - A10 * A01;

        	A00 /= det; A01 /= (-det); A10 /= (-det); A11 /= det;

        	float rhsx = xi - x1p * (c1 + c2) - x4p * (c3 + c4);
        	float rhsy = yi - y1p * (c1 + c2) - y4p * (c3 + c4);

        	float a1 = A11 * rhsx + A01 * rhsy;
        	float a2 = A10 * rhsx + A00 * rhsy;

        	float x2p, y2p, x3p, y3p;
        	x2p = a1 * (x2 - x1) + x1p;
        	y2p = a1 * (y2 - y1) + y1p;
        	x3p = a2 * (x3 - x4) + x4p;
        	y3p = a2 * (y3 - y4) + y4p;

//        	right[0] = x1p; right[1] = y1p;
//        	right[2] = x2p; right[3] = y2p;
//        	right[4] = x3p; right[5] = y3p;
//        	right[6] = x4p; right[7] = y4p;
        	curarr[0] = x1p; curarr[1] = y1p;
        	curarr[2] = x2p; curarr[3] = y2p;
        	curarr[4] = x3p; curarr[5] = y3p;
        	curarr[6] = x4p; curarr[7] = y4p;

        	x1p = x1 - offset[0][0]; y1p = y1 - offset[0][1];
        	xi = xi - 2 * offset[1][0]; yi = yi - 2 * offset[1][1];
        	x4p = x4 - offset[2][0]; y4p = y4 - offset[2][1];
        	curarr = left;
        }
        System.out.println("right: " + Arrays.toString(right));
        System.out.println("left: " + Arrays.toString(left));
    }
    
    private void computeOffsetWay4(float[] pts, float[] right, float[] left) {
    	float x1 = pts[0], y1 = pts[1];
    	float x2 = pts[2], y2 = pts[3];
    	float x3 = pts[4], y3 = pts[5];
    	float x4 = pts[6], y4 = pts[7];
    	BezCurve cur = BezCurve.makeCheapBezCurve(null, pts, 8);

        final float dx4 = x4 - x3;
        final float dy4 = y4 - y3;
        float dx1 = x2 - x1;
        float dy1 = y2 - y1;

        // this computes the offsets at t=0, 0.5, 1, using the property that
        // for any bezier curve the vectors p2-p1 and p4-p3 are parallel to
        // the (dx/dt, dy/dt) vectors at the endpoints.
        computeOffset(dx1, dy1, offset[0]);
        computeOffset(cur.dxat(0.5f), cur.dyat(0.5f), offset[1]);
        computeOffset(dx4, dy4, offset[2]);
        float x1p = x1 + offset[0][0]; // start
        float y1p = y1 + offset[0][1]; // point
        float xi = cur.xat(0.5f) + offset[1][0];
        float yi = cur.yat(0.5f) + offset[1][1];
        float x4p = x4 + offset[2][0]; // end
        float y4p = y4 + offset[2][1]; // point

        float c = (4 * (2 * xi - x1p - x4p)) / (3f * (x2 - x1 - x4 + x3));
        float c2 = (4 * (2 * yi - y1p - y4p)) / (3f * (y2 - y1 - y4 + y3));
        System.out.println("c is : " + c);
        System.out.println("c2 is : " + c2);
        
        float x2p, y2p, x3p, y3p;
        x2p = c * (x2 - x1) + x1p; 
        y2p = c2 * (y2 - y1) + y1p;
        x3p = x4p - c * (x4 - x3);
        y3p = y4p - c2 * (y4 - y3);

        right[0] = x1p; right[1] = y1p;
        right[2] = x2p; right[3] = y2p;
        right[4] = x3p; right[5] = y3p;
        right[6] = x4p; right[7] = y4p;
        
        x1p = x1 - offset[0][0]; y1p = y1 - offset[0][1];
        xi = xi - 2 * offset[1][0]; yi = yi - 2 * offset[1][1];
        x4p = x4 - offset[2][0]; y4p = y4 - offset[2][1];
        
        c = (4 * (2 * xi - x1p - x4p)) / (3f * (x2 - x1 - x4 + x3));
        c2 = (4 * (2 * yi - y1p - y4p)) / (3f * (y2 - y1 - y4 + y3));
        
        x2p = c * (x2 - x1) + x1p; 
        y2p = c2 * (y2 - y1) + y1p;
        x3p = x4p - c * (x4 - x3);
        y3p = y4p - c2 * (y4 - y3);
        
        left[0] = x1p; left[1] = y1p;
        left[2] = x2p; left[3] = y2p;
        left[4] = x3p; left[5] = y3p;
        left[6] = x4p; left[7] = y4p;
    }
    
    BezCurve c = null;
    
    @Override public void curveTo(float x1o, float y1o, 
    		float x2o, float y2o, 
    		float x3o, float y3o) {
    	computeOffset(x1o - x0, y1o - y0, offset[0]);
    	float mx = offset[0][0];
    	float my = offset[0][1];

    	c = BezCurve.makeCheapBezCurve(c, new float[] {x0, y0, x1o, y1o, x2o, y2o, x3o, y3o}, 8);

    	drawJoin(px0, py0, x0, y0, x1o, y1o, omx, omy, mx, my);

    	emitLineTo(x0 + mx, y0 + my);

    	int nCurves = c.breakIntoMonotonicAndInflectionpts();
//    	int nCurves = c.breakIntoMonotonicPieces();
//    	int nCurves = c.breakIntoNoInfs();

    	for (int i = 0; i < nCurves; i++) {
//    		final float[] pts = new float[] {x0, y0, x1o, y1o, x2o, y2o, x3o, y3o};
    		final float[] pts = c.getCurve(i);

    		float[] r = new float[8];
    		float[] l = new float[8];
    		int way = 3;
    		switch(way) {
    		case 1:
    			computeOffsetWay1(pts, r, l);
    			break;
    		case 2:
    			computeOffsetWay2(pts, r, l);
    			break;
    		case 3:
    			computeOffsetWay3(pts, r, l);
    			break;
    		case 4:
    			computeOffsetWay4(pts, r, l);
    			break;
    		}
    		emitCurveTo(r[2], r[3], r[4], r[5], r[6], r[7]);
    		
    		reverse.pushCubic(l[0], l[1], l[2], l[3], l[4], l[5]);
    	}

    	computeOffset(x3o - x2o, y3o - y2o, offset[0]);
    	mx = offset[0][0];
    	my = offset[0][1];
    	emitLineTo(x3o - mx, y3o - my, true);

    	this.omx = mx;
    	this.omy = my;
    	this.px0 = x2o;
    	this.py0 = y2o;
    	this.x0 = x3o;
    	this.y0 = y3o;
    	// not technically true, but it shouldn't matter TODO: look into it
    	this.prev = LINE_TO;
    }
    
    // TODO: must check for cusps (places where |dx,dy| is 0 or very small.
    // For example, the curve (0,0),(1,1),(0,1),(1,0) has |dx,dy| == 0 at t=0.5
    // and so computeOffset will return (0,0).
/*    @Override public void curveTo(float x1o, float y1o, 
                        float x2o, float y2o, 
                        float x3o, float y3o) {
    	computeOffset(x1o - x0, y1o - y0, offset[0]);
    	float mx = offset[0][0];
    	float my = offset[0][1];
    	
//    	BezCurve c = new BezCurve(x0, y0, x1o, y1o, x2o, y2o, x3o, y3o);
    	c = BezCurve.makeCheapBezCurve(c, new float[] {x0, y0, x1o, y1o, x2o, y2o, x3o, y3o}, 8);

        drawJoin(px0, py0, x0, y0, x1o, y1o, omx, omy, mx, my);

        emitLineTo(x0 + mx, y0 + my);

        
//        int nCurves = getNumCurvesCubic(x1o, y1o, x2o, y2o, x3o, y3o);
        int nCurves = c.breakIntoMonotonicPieces();

        // We split the curve to be widened into nCurves curves. For each curve,
        // we compute 2 offset curves. Let p1,p2,p3,p4 be the control points of
        // the middle (non offset) curve, and p1', p2', p3', p4' be the control
        // points of whichever offset curve is being currently computed. (we 
        // can't have "'" in variable names so we denote these with p1p, p2p,
        // p3p, p4p) This loop constructs 2 bezier splines to approximate the
        // ideal offset curves. The splines are composed of nCurves bezier 
        // curves. Each is computed using the constraints that p2-p1 is parallel
        // to p2p-p1p, p3-p4 is parallel to p3p-p4p, and that the bezier defined
        // by p1p,p2p,p3p,p4p must pass through B(0.5)+/-offset (where B(0.5) 
        // is the curve p1,p2,p3,p4 at t=0.5)

//        final float t_inc = 1f / nCurves;
//        float t_end = t_inc, t_middle = t_inc / 2;
//        float x4 = c.xat(0), y4 = c.yat(0);
        float dx4 = c.dxat(0); // this is very ugly. TODO: make a mutable Point/Vector class?
        float dy4 = c.dyat(0);
        for (int i = 0; i < nCurves; i++) {
        	final float[] pts = c.getCurve(i);//new float[] {x0, y0, x1o, y1o, x2o, y2o, x3o, y3o};//c.getCurve(i);
        	final float t_start = c.getLowT(i);//0;//c.getLowT(i);
        	final float t_end = c.getHiT(i);//1;//c.getHiT(i);
        	final float x1 = pts[0], y1 = pts[1];
        	final float x2 = pts[2], y2 = pts[3];
        	final float x3 = pts[4], y3 = pts[5];
            final float x4 = pts[6], y4 = pts[7];

            final float dx1 = dx4;
            final float dy1 = dy4;
            dx4 = c.dxat(t_end);
            dy4 = c.dyat(t_end);

            // this computes the offsets at t=0, 0.5, 1, using the property that
            // for any bezier curve the vectors p2-p1 and p4-p3 are parallel to
            // the (dx/dt, dy/dt) vectors at the endpoints.
            computeOffset(dx1, dy1, offset[0]);
            computeOffset(dx4, dy4, offset[2]);
            float x1p = x1 + offset[0][0]; // start
            float y1p = y1 + offset[0][1]; // point
            float x4p = x4 + offset[2][0]; // end
            float y4p = y4 + offset[2][1]; // point

            float x2p, y2p, x3p, y3p;
            computeOffset(c.dxat((1 - (1f/3)) * t_start + (1f/3)*t_end), c.dyat((1 - (1f/3)) * t_start + (1f/3)*t_end), offset[1]);
            float mx2 = offset[1][0]; float my2 = offset[1][1];
            x2p = x2 + offset[1][0];    y2p = y2 + offset[1][1];
            computeOffset(c.dxat((1 - (2f/3)) * t_start + (2f/3)*t_end), c.dyat((1 - (2f/3)) * t_start + (2f/3)*t_end), offset[1]);
            float mx3 = offset[1][0]; float my3 = offset[1][1];
            x3p = x3 + offset[1][0];    y3p = y3 + offset[1][1];

            emitCurveTo(x2p, y2p, x3p, y3p, x4p, y4p);
            
            x1p = x1 - offset[0][0]; y1p = y1 - offset[0][1];
            x2p = x2 - mx2;    y2p = y2 - my2;
            x3p = x3 - mx3;    y3p = y3 - my3;
            x4p = x4 - offset[2][0]; y4p = y4 - offset[2][1];
            
            reverse.pushCubic(x1p, y1p, x2p, y2p, x3p, y3p);
//            t_middle += t_inc;
//            t_end += t_inc;
        }
        
        computeOffset(x3o - x2o, y3o - y2o, offset[0]);
        mx = offset[0][0];
        my = offset[0][1];
        emitLineTo(x3o - mx, y3o - my, true);
        
        this.omx = mx;
        this.omy = my;
        this.px0 = x2o;
        this.py0 = y2o;
        this.x0 = x3o;
        this.y0 = y3o;
        // not technically true, but it shouldn't matter TODO: look into it
        this.prev = LINE_TO;
    }*/

    @Override
    public long getNativeConsumer() {
        throw new InternalError("Stroker doesn't use a native consumer");
    }

    private int numCurvesQuad(float x0, float y0, float x1, float y1, float x2, float y2) {
    	return 13;
    }
    
    @Override
    public void quadTo(float x1o, float y1o, float x2o, float y2o) {
    	computeOffset(x1o - x0, y1o - y0, offset[0]);
    	float mx = offset[0][0];
    	float my = offset[0][1];

        drawJoin(px0, py0, x0, y0, x1o, y1o, omx, omy, mx, my);

        emitLineTo(x0 + mx, y0 + my);

        int nCurves = numCurvesQuad(x0, y0, x1o, y1o, x2o, y2o);
        BezCurve c = new BezCurve(x0, y0, x1o, y1o, x2o, y2o);
        final float t_inc = 1f / nCurves;
        float t_end = t_inc;
        float x3 = c.xat(0), y3 = c.yat(0);
        float dx3 = c.dxat(0);
        float dy3 = c.dyat(0);
        for (int i = 0; i < nCurves; i++) {
            float x1 = x3, y1 = y3;
            x3 = c.xat(t_end);
            y3 = c.yat(t_end);
            
            float dx1 = dx3;
            float dy1 = dy3;
            dx3 = c.dxat(t_end);
            dy3 = c.dyat(t_end);

            // this computes the offsets at t=0, 0.5, 1, using the property that
            // for any bezier curve the vectors p2-p1 and p4-p3 are parallel to
            // the (dx/dt, dy/dt) vectors at the endpoints.
            computeOffset(dx1, dy1, offset[0]);
            computeOffset(dx3, dy3, offset[1]);
            float x1p = x1 + offset[0][0]; // start
            float y1p = y1 + offset[0][1]; // point
            float x3p = x3 + offset[1][0]; // end
            float y3p = y3 + offset[1][1]; // point
            
            // reciprocal of the determinant L
            float detr = 1f / (dx1 * dy3 - dy1 * dx3);

            float p3pmp1px = x3p - x1p;
            float p3pmp1py = y3p - y1p;
            float c1 = detr * (dy3 * p3pmp1px - dx3 * p3pmp1py);
            
            float x2p, y2p;
            x2p = x1p + c1*dx1;
            y2p = y1p + c1*dy1;

            emitQuadTo(x2p, y2p, x3p, y3p);
            
            x1p = x1 - offset[0][0]; y1p = y1 - offset[0][1];
            x3p = x3 - offset[1][0]; y3p = y3 - offset[1][1];
            
            p3pmp1px = x3p - x1p;
            p3pmp1py = y3p - y1p;
            c1 = detr * (dy3 * p3pmp1px - dx3 * p3pmp1py);
            
            x2p = x1p + c1*dx1;
            y2p = y1p + c1*dy1;
            
            reverse.pushQuad(x1p, y1p, x2p, y2p);
            t_end += t_inc;
        }
        
        computeOffset(x2o - x1o, y2o - y1o, offset[0]);
        mx = offset[0][0];
        my = offset[0][1];
        emitLineTo(x2o - mx, y2o - my, true);
        
        this.omx = mx;
        this.omy = my;
        this.px0 = x2o;
        this.py0 = y2o;
        this.x0 = x2o;
        this.y0 = y2o;
        // not true, but it shouldn't matter TODO: look into it
        this.prev = LINE_TO;
    }
    
    private static int sizeOfPoly(int type) {
    	switch(type) {
    	case LINE_TO: return 4;
    	case QUAD_TO: return 6;
    	case CURVE_TO: return 8;
    	default: throw new InternalError("Unrecognized curve type");
    	}
    }
    
    // a stack of polynomial curves where each curve shares endpoints with
    // adjacent ones.
    private static final class PolyStack {
        float[] curves;
        int end;
        int[] curveTypes;
        int numCurves;
        
        private static final int INIT_SIZE = 50;
        
        PolyStack() {
            curves = new float[8 * INIT_SIZE];
            curveTypes = new int[INIT_SIZE];
            end = 0;
            numCurves = 0;
        }
        
        public boolean isEmpty() {
        	return numCurves == 0;
        }
        
        private void ensureSpace(int n) {
            if (end + n >= curves.length) {
                int newSize = (end + n) * 2;
                curves = Arrays.copyOf(curves, newSize);
            }
            if (numCurves >= curveTypes.length) {
                int newSize = numCurves * 2;
                curveTypes = Arrays.copyOf(curveTypes, newSize);
            }
        }
        
        public void pushCubic(float x0, float y0,
                              float x1, float y1,
                              float x2, float y2) {
            ensureSpace(6);
            curveTypes[numCurves++] = CURVE_TO;
            // assert(x0 == lastX && y0 == lastY)
            
            // we reverse the coordinate order to make popping easier
            curves[end++] = x2;    curves[end++] = y2;
            curves[end++] = x1;    curves[end++] = y1;
            curves[end++] = x0;    curves[end++] = y0;
        }
        
        public void pushQuad(float x0, float y0,
                             float x1, float y1) {
            ensureSpace(4);
            curveTypes[numCurves++] = QUAD_TO;
            // assert(x0 == lastX && y0 == lastY)
            curves[end++] = x1;    curves[end++] = y1;
            curves[end++] = x0;    curves[end++] = y0;
        }
        
        public void pushLine(float x, float y) {
            ensureSpace(2);
            curveTypes[numCurves++] = LINE_TO;
            // assert(x0 == lastX && y0 == lastY)
            curves[end++] = x;    curves[end++] = y;
        }

        public int pop(float[] pts) {
        	int ret = curveTypes[numCurves - 1];
        	int numCoords = sizeOfPoly(ret) - 2;
        	numCurves--;
        	end -= numCoords;
        	System.arraycopy(curves, end, pts, 0, numCoords);
        	return ret;
        }
    }
}

Reply via email to