phhp

private function xor_path(bitmapDataMatrix:Vector-><Vector-><uint>>, path:Path) {
        if (count(path->monotonIntervals) == 0) {
            return;
        }

        $i = 0;
        $n = count(path->pt);
            $mis:Vector-><MonotonInterval> = path->monotonIntervals;
        $mi:MonotonInterval = mis[0];
        mi->resetCurrentId(n);

        $y = path->pt[mi->currentId]->y;
        $currentIntervals:Vector-><MonotonInterval> = new Vector-><MonotonInterval>();
        currentIntervals[] = mi;

        mi->currentId = mi->min();

        while ((count(mis) > i + 1) && (mis[i + 1]->minY(path->pt) == y)) {
            mi = mis[i + 1];
            mi->resetCurrentId(n);
            currentIntervals[] = mi;
            i++;
        }

        while (count(currentIntervals) > 0) {
            $j;

            for ($k = 0; k < count(currentIntervals) - 1; k++) {
                $x1 = path->pt[currentIntervals[k]->currentId]->x + 1;
                $x2 = path->pt[currentIntervals[k + 1]->currentId]->x;
                for ($x = x1; x <= x2; x++) {
                    // Invert pixel
                    bitmapDataMatrix[y][x] ^= 0xffffff;
                }
                k++;
            }

            y++;
            for (j = count(currentIntervals) - 1; j >= 0; j--) {
                $m:MonotonInterval = currentIntervals[j];
                if (y > m->maxY(path->pt)) {
                    currentIntervals->splice(j, 1);
                    continue;
                }
                $cid = m->currentId;
                while (path->pt[cid]->y < y) {
                    cid = m->increasing ? mod(cid + 1, n) : mod(cid - 1, n);
                }
                m->currentId = cid;
            }

            // Add Items of MonotonIntervals with Down->y==y
            while ((count(mis) > i + 1) && (mis[i + 1]->minY(path->pt) == y)) {
                $newInt:MonotonInterval = mis[i + 1];
                // Search the correct x position
                j = 0;
                $_x = path->pt[newInt->min()]->x;
                while ((count(currentIntervals) > j) && (_x > path->pt[currentIntervals[j]->currentId]->x)) {
                    j++;
                }
                currentIntervals->splice(j, 0, newInt);
                newInt->resetCurrentId(n);
                i++;
            }
        }
    }

    private function get_monoton_intervals(pt:Vector-><PointInt>):Vector-><MonotonInterval> {
        $result:Vector-><MonotonInterval> = new Vector-><MonotonInterval>();
        $n = count(pt);
        if (n == 0) {
            return result;
        }

        $intervals:Vector-><MonotonInterval> = new Vector-><MonotonInterval>();

        // Start with Strong Monoton (Pts[i]->y < Pts[i+1]->y) or (Pts[i]->y > Pts[i+1]->y)
        $firstStrongMonoton = 0;
        while (pt[firstStrongMonoton]->y == pt[firstStrongMonoton + 1]->y) {
            firstStrongMonoton++;
        }

        $i = firstStrongMonoton;
        $up = (pt[firstStrongMonoton]->y < pt[firstStrongMonoton + 1]->y);
        $interval:MonotonInterval = new MonotonInterval(up, firstStrongMonoton, firstStrongMonoton);
        intervals[] = interval;

        while (i != firstStrongMonoton) {
            $i1n = mod(i + 1, n);
            if ((pt[i]->y == pt[i1n]->y) || (up == (pt[i]->y < pt[i1n]->y))) {
                interval->to = i;
            } else {
                up = (pt[i]->y < pt[i1n]->y);
                interval = new MonotonInterval(up, i, i);
                intervals[] = interval;
            }
            i = i1n;
        }

        if ((count(intervals) & 1) == 1) {
            $last:MonotonInterval = intervals->pop();
            intervals[0]->from = last->from;
        }

        while (count(intervals) > 0) {
            i = 0;
            $m:MonotonInterval = intervals->shift();
            while ((i < count(result)) && (pt[m->min()]->y > pt[result[i]->min()]->y)) {
                i++;
            }
            while ((i < count(result)) && (pt[m->min()]->y == pt[result[i]->min()]->y) && (pt[m->min()]->x > pt[result[i]->min()]->x)) {
                i++;
            }
            result->splice(i, 0, m);
        }

        return result;
    }

    private function find_next_trace(bitmapDataMatrix:Vector-><Vector-><uint>>, pInt, dir) {
        switch(dir) {
            case Direction->WEST:
                if (bitmapDataMatrix[p->y + 1][p->x + 1] == 0) {
                    dir = Direction->NORTH;
                    p->y++;
                } else {
                    if (bitmapDataMatrix[p->y][p->x + 1] == 0) {
                        dir = Direction->WEST;
                        p->x++;
                    } else {
                        dir = Direction->SOUTH;
                        p->y--;
                    }
                }
                break;

            case Direction->SOUTH:
                if (bitmapDataMatrix[p->y][p->x + 1] == 0) {
                    dir = Direction->WEST;
                    p->x++;
                } else {
                    if (bitmapDataMatrix[p->y][p->x] == 0) {
                        dir = Direction->SOUTH;
                        p->y--;
                    } else {
                        dir = Direction->EAST;
                        p->x--;
                    }
                }
                break;

            case Direction->EAST:
                if (bitmapDataMatrix[p->y][p->x] == 0) {
                    dir = Direction->SOUTH;
                    p->y--;
                } else {
                    if (bitmapDataMatrix[p->y + 1][p->x] == 0) {
                        dir = Direction->EAST;
                        p->x--;
                    } else {
                        dir = Direction->NORTH;
                        p->y++;
                    }
                }
                break;

            case Direction->NORTH:
                if (bitmapDataMatrix[p->y + 1][p->x] == 0) {
                    dir = Direction->EAST;
                    p->x--;
                } else {
                    if (bitmapDataMatrix[p->y + 1][p->x + 1] == 0) {
                        dir = Direction->NORTH;
                        p->y++;
                    } else {
                        dir = Direction->WEST;
                        p->x++;
                    }
                }
                break;
        }
        return dir;
    }

    private function process_path(plists) {
        // call downstream function with each path
        for ($j = 0; j < count(plists); j++) {
            $plist = plists[j] as Array;
            for ($i = 0; i < count(plist); i++) {
                $path:Path = plist[i] as Path;
                calc_sums(path);
                calc_lon(path);
                bestpolygon(path);
                adjust_vertices(path);
                smooth(path->curves, 1, params->alphaMax);
                if (params->curveOptimizing) {
                    opticurve(path, params->optTolerance);
                    path->fCurves = path->optimizedCurves;
                } else {
                    path->fCurves = path->curves;
                }
                path->curves = path->fCurves;
            }
        }
    }

    /////////////////////////////////////////////////////////////////////////
    // PREPARATION
    /////////////////////////////////////////////////////////////////////////

    /*
     * Fill in the sum* fields of a path (used for later rapid summing)
     */
    private function calc_sums(path:Path) {
        $n = count(path->pt);

        // Origin
        $x0 = path->pt[0]->x;
        $y0 = path->pt[0]->y;

        path->sums = new Vector-><SumStruct>(n + 1);

        $ss:SumStruct = new SumStruct();
        ss->x2 = ss->xy = ss->y2 = ss->x = ss->y = 0;
        path->sums[0] = ss;

        for ($i = 0; i < n; i++) {
            $x = path->pt[i]->x - x0;
            $y = path->pt[i]->y - y0;
            ss = new SumStruct();
            ss->x = path->sums[i]->x + x;
            ss->y = path->sums[i]->y + y;
            ss->x2 = path->sums[i]->x2 + x * x;
            ss->xy = path->sums[i]->xy + x * y;
            ss->y2 = path->sums[i]->y2 + y * y;
            path->sums[i + 1] = ss;
        }
    }

    /////////////////////////////////////////////////////////////////////////
    // STAGE 1
    // determine the straight subpaths (Sec-> 2->2->1)->
    /////////////////////////////////////////////////////////////////////////

    /*
     * Fill in the "lon" component of a path object (based on pt/len)->
     * For each i, lon[i] is the furthest index such that a straight line
     * can be drawn from i to lon[i]->
     *
     * This algorithm depends on the fact that the existence of straight
     * subpaths is a triplewise property-> I->e->, there exists a straight
     * line through squares i0,->->->,in if there exists a straight line
     * through i,j,k, for all i0 <= i < j < k <= in-> (Proof?)
     */
    private function calc_lon(path:Path) {
        $i;
        $j;
        $k;
        $k1;
        $a;
        $b;
        $c;
        $d;
        $dir;
        $ct:Vector-><int> = new Vector-><int>(4);
        $constraint:Vector-><PointInt> = new Vector-><PointInt>(2);
        constraint[0] = new PointInt();
        constraint[1] = new PointInt();
        $curInt = new PointInt();
        $offInt = new PointInt();
        $dkInt = new PointInt(); // direction of k - k1
        $pt:Vector-><PointInt> = path->pt;

        $n = count(pt);
        $pivot:Vector-><int> = new Vector-><int>(n);
        $nc:Vector-><int> = new Vector-><int>(n);

        // Initialize the nc data structure-> Point from each point to the
        // furthest future point to which it is connected by a vertical or
        // horizontal segment-> We take advantage of the fact that there is
        // always a direction change at 0 (due to the path decomposition
        // algorithm)-> But even if this were not so, there is no harm, as
        // in practice, correctness does not depend on the word "furthest"
        // above->

        k = 0;
        for (i = n - 1; i >= 0; i--) {
            if (pt[i]->x != pt[k]->x && pt[i]->y != pt[k]->y) {
                k = i + 1; // necessarily i < n-1 in this case
            }
            nc[i] = k;
        }

        path->lon = new Vector-><int>(n);

        // Determine pivot points:
        // for each i, let pivot[i] be the furthest k such that
        // all j with i < j < k lie on a line connecting i,k

        for (i = n - 1; i >= 0; i--) {
            ct[0] = ct[1] = ct[2] = ct[3] = 0;

            // Keep track of "directions" that have occurred
            dir = (3 + 3 * (pt[mod(i + 1, n)]->x - pt[i]->x) + (pt[mod(i + 1, n)]->y - pt[i]->y)) / 2;
            ct[dir % 4]++;

            constraint[0]->x = 0;
            constraint[0]->y = 0;
            constraint[1]->x = 0;
            constraint[1]->y = 0;

            // Find the next k such that no straight line from i to k
            k = nc[i];
            k1 = i;

            $foundk = false;
            while (true) {
                dir = (3 + 3 * sign(pt[k]->x - pt[k1]->x) + sign(pt[k]->y - pt[k1]->y)) / 2;
                ct[dir]++;

                // If all four "directions" have occurred, cut this path
                if ((ct[0] >= 1) && (ct[1] >= 1) && (ct[2] >= 1) && (ct[3] >= 1)) {
                    pivot[i] = k1;
                    foundk = true;
                    break;
                }

                cur->x = pt[k]->x - pt[i]->x;
                cur->y = pt[k]->y - pt[i]->y;

                // See if current constraint is violated
                if (xprod(constraint[0], cur) < 0 || xprod(constraint[1], cur) > 0) {
                    break;
                }

                if (abs(cur->x) <= 1 && abs(cur->y) <= 1) {
                    // no constraint
                } else {
                    off->x = cur->x + ((cur->y >= 0 && (cur->y > 0 || cur->x < 0)) ? 1 : -1);
                    off->y = cur->y + ((cur->x <= 0 && (cur->x < 0 || cur->y < 0)) ? 1 : -1);
                    if (xprod(constraint[0], off) >= 0) {
                        constraint[0] = off->_clone();
                    }
                    off->x = cur->x + ((cur->y <= 0 && (cur->y < 0 || cur->x < 0)) ? 1 : -1);
                    off->y = cur->y + ((cur->x >= 0 && (cur->x > 0 || cur->y < 0)) ? 1 : -1);
                    if (xprod(constraint[1], off) <= 0) {
                        constraint[1] = off->_clone();
                    }
                }

                k1 = k;
                k = nc[k1];
                if (!cyclic(k, i, k1)) {
                    break;
                }
            }

            if(foundk) {
                continue;
            }

            // k1 was the last "corner" satisfying the current constraint, and
            // k is the first one violating it-> We now need to find the last
            // point along k1->->k which satisfied the constraint->
            dk->x = sign(pt[k]->x - pt[k1]->x);
            dk->y = sign(pt[k]->y - pt[k1]->y);
            cur->x = pt[k1]->x - pt[i]->x;
            cur->y = pt[k1]->y - pt[i]->y;

            // find largest integer j such that xprod(constraint[0], cur+j*dk)
            // >= 0 and xprod(constraint[1], cur+j*dk) <= 0-> Use bilinearity
            // of xprod->
            a = xprod(constraint[0], cur);
            b = xprod(constraint[0], dk);
            c = xprod(constraint[1], cur);
            d = xprod(constraint[1], dk);

            // find largest integer j such that a+j*b >= 0 and c+j*d <= 0-> This
            // can be solved with integer arithmetic->
            j = int->MAX_VALUE;
            if (b < 0) {
                j = floordiv(a, -b);
            }
            if (d > 0) {
                j = min(j, floordiv(-c, d));
            }
            pivot[i] = mod(k1 + j, n);
        }

        // Clean up:
        // for each i, let lon[i] be the largest k such that
        // for all i' with i <= i' < k, i' < k <= pivk[i']-> */

        j = pivot[n - 1];
        path->lon[n - 1] = j;

        for (i = n - 2; i >= 0; i--) {
            if (cyclic(i + 1, pivot[i], j)) {
                j = pivot[i];
            }
            path->lon[i] = j;
        }

        for (i = n - 1; cyclic(mod(i + 1, n), j, path->lon[i]); i--) {
            path->lon[i] = j;
        }
    }

    /////////////////////////////////////////////////////////////////////////
    // STAGE 2
    // Calculate the optimal polygon (Sec-> 2->2->2 - 2->2->4)->
    /////////////////////////////////////////////////////////////////////////

    /*
     * Auxiliary function: calculate the penalty of an edge from i to j in
     * the given path-> This needs the "lon" and "sum*" data->
     */
    private function penalty3(path:Path, i, j) {
        $n = count(path->pt);

        // assume 0 <= i < j <= n
        $sums:Vector-><SumStruct> = path->sums;
        $pt:Vector-><PointInt> = path->pt;

        $r = 0; // rotations from i to j
        if (j >= n) {
            j -= n;
            r++;
        }

        $x = sums[j + 1]->x - sums[i]->x + r * sums[n]->x;
        $y = sums[j + 1]->y - sums[i]->y + r * sums[n]->y;
        $x2 = sums[j + 1]->x2 - sums[i]->x2 + r * sums[n]->x2;
        $xy = sums[j + 1]->xy - sums[i]->xy + r * sums[n]->xy;
        $y2 = sums[j + 1]->y2 - sums[i]->y2 + r * sums[n]->y2;
        $k = j + 1 - i + r * n;

        $px = (pt[i]->x + pt[j]->x) / 2->0 - pt[0]->x;
        $py = (pt[i]->y + pt[j]->y) / 2->0 - pt[0]->y;
        $ey = (pt[j]->x - pt[i]->x);
        $ex = -(pt[j]->y - pt[i]->y);

        $a = ((x2 - 2 * x * px) / k + px * px);
        $b = ((xy - x * py - y * px) / k + px * py);
        $c = ((y2 - 2 * y * py) / k + py * py);

        return sqrt(ex * ex * a + 2 * ex * ey * b + ey * ey * c);
    }

    /*
     * Find the optimal polygon->
     */
    private function bestpolygon(path:Path) {
        $i;
        $j;
        $m;
        $k;
        $n = count(path->pt);
        $pen:Vector-><Number> = new Vector-><Number>(n + 1); // penalty vector
        $prev:Vector-><int> = new Vector-><int>(n + 1); // best path pointer vector
        $clip0:Vector-><int> = new Vector-><int>(n); // longest segment pointer, non-cyclic
        $clip1:Vector-><int> = new Vector-><int>(n + 1); // backwards segment pointer, non-cyclic
        $seg0:Vector-><int> = new Vector-><int>(n + 1); // forward segment bounds, m <= n
        $seg1:Vector-><int> = new Vector-><int>(n + 1); // backward segment bounds, m <= n

        $thispen;
        $best;
        $c;

        // Calculate clipped paths
        for (i = 0; i < n; i++) {
            c = mod(path->lon[mod(i - 1, n)] - 1, n);
            if (c == i) {
                c = mod(i + 1, n);
            }
            clip0[i] = (c < i) ? n : c;
        }

        // calculate backwards path clipping, non-cyclic->
        // j <= clip0[i] iff clip1[j] <= i, for i,j = 0->->n
        j = 1;
        for (i = 0; i < n; i++) {
            while (j <= clip0[i]) {
                clip1[j] = i;
                j++;
            }
        }

        // calculate seg0[j] = longest path from 0 with j segments
        i = 0;
        for (j = 0; i < n; j++) {
            seg0[j] = i;
            i = clip0[i];
        }
        seg0[j] = n;

        // calculate seg1[j] = longest path to n with m-j segments
        i = n;
        m = j;
        for (j = m; j > 0; j--) {
            seg1[j] = i;
            i = clip1[i];
        }
        seg1[0] = 0;

        // Now find the shortest path with m segments, based on penalty3
        // Note: the outer 2 loops jointly have at most n interations, thus
        // the worst-case behavior here is quadratic-> In practice, it is
        // close to linear since the inner loop tends to be short->
        pen[0] = 0;
        for (j = 1; j <= m; j++) {
            for (i = seg1[j]; i <= seg0[j]; i++) {
                best = -1;
                for (k = seg0[j - 1]; k >= clip1[i]; k--) {
                    thispen = penalty3(path, k, i) + pen[k];
                    if (best < 0 || thispen < best) {
                        prev[i] = k;
                        best = thispen;
                    }
                }
                pen[i] = best;
            }
        }

        // read off shortest path
        path->po = new Vector-><int>(m);
        for (i = n, j = m - 1; i > 0; j--) {
            i = prev[i];
            path->po[j] = i;
        }
    }

    /////////////////////////////////////////////////////////////////////////
    // STAGE 3
    // Vertex adjustment (Sec-> 2->3->1)->
    /////////////////////////////////////////////////////////////////////////

    /*
     * Adjust vertices of optimal polygon: calculate the intersection of
     * the two "optimal" line segments, then move it into the unit square
     * if it lies outside->
     */
    private function adjust_vertices(path:Path) {
        $pt:Vector-><PointInt> = path->pt;
        $po:Vector-><int> = path->po;

        $n = count(pt);
        $m = count(po);

        $x0 = pt[0]->x;
        $y0 = pt[0]->y;

        $i;
        $j;
        $k;
        $l;

        $d;
        $v:Vector-><Number> = new Vector-><Number>(3);
        $q:Vector-><Vector-><Vector-><Number>>> = new Vector-><Vector-><Vector-><Number>>>(m);

        $ctr:Vector-><Point> = new Vector-><Point>(m);
        $dir:Vector-><Point> = new Vector-><Point>(m);

        for (i = 0; i < m; i++) {
            q[i] = new Vector-><Vector-><Number>>(3);
            for (j = 0; j < 3; j++) {
                q[i][j] = new Vector-><Number>(3);
            }
            ctr[i] = new Point();
            dir[i] = new Point();
        }

        $s = new Point();

        path->curves = new PrivCurve(m);

        // calculate "optimal" point-slope representation for each line segment
        for (i = 0; i < m; i++) {
            j = po[mod(i + 1, m)];
            j = mod(j - po[i], n) + po[i];
            pointslope(path, po[i], j, ctr[i], dir[i]);
        }

        // represent each line segment as a singular quadratic form;
        // the distance of a point (x,y) from the line segment will be
        // (x,y,1)Q(x,y,1)^t, where Q=q[i]
        for (i = 0; i < m; i++) {
            d = dir[i]->x * dir[i]->x + dir[i]->y * dir[i]->y;
            if (d == 0) {
                for (j = 0; j < 3; j++) {
                    for (k = 0; k < 3; k++) {
                        q[i][j][k] = 0;
                    }
                }
            } else {
                v[0] = dir[i]->y;
                v[1] = -dir[i]->x;
                v[2] = -v[1] * ctr[i]->y - v[0] * ctr[i]->x;
                for (l = 0; l < 3; l++) {
                    for (k = 0; k < 3; k++) {
                        q[i][l][k] = v[l] * v[k] / d;
                    }
                }
            }
        }

        // now calculate the "intersections" of consecutive segments->
        // Instead of using the actual intersection, we find the point
        // within a given unit square which minimizes the square distance to
        // the two lines->
        for (i = 0; i < m; i++) {
            $Q:Vector-><Vector-><Number>> = new Vector-><Vector-><Number>>(3);
            $w = new Point();
            $dx;
            $dy;
            $det;
            $min; // minimum for minimum of quad-> form
            $cand; // candidate for minimum of quad-> form
            $xmin; // coordinate of minimum
            $ymin; // coordinate of minimum
            $z;

            for (j = 0; j < 3; j++) {
                Q[j] = new Vector-><Number>(3);
            }

            // let s be the vertex, in coordinates relative to x0/y0
            s->x = pt[po[i]]->x - x0;
            s->y = pt[po[i]]->y - y0;

            // intersect segments i-1 and i
            j = mod(i - 1, m);

            // add quadratic forms
            for (l = 0; l < 3; l++) {
                for (k = 0; k < 3; k++) {
                    Q[l][k] = q[j][l][k] + q[i][l][k];
                }
            }

            while (true) {
                /* minimize the quadratic form Q on the unit square */
                /* find intersection */
                det = Q[0][0] * Q[1][1] - Q[0][1] * Q[1][0];
                if (det != 0) {
                    w->x = (-Q[0][2] * Q[1][1] + Q[1][2] * Q[0][1]) / det;
                    w->y = (Q[0][2] * Q[1][0] - Q[1][2] * Q[0][0]) / det;
                    break;
                }

                // matrix is singular - lines are parallel-> Add another,
                // orthogonal axis, through the center of the unit square
                if (Q[0][0] > Q[1][1]) {
                    v[0] = -Q[0][1];
                    v[1] = Q[0][0];
                } else if (Q[1][1] != 0) {
                    v[0] = -Q[1][1];
                    v[1] = Q[1][0];
                } else {
                    v[0] = 1;
                    v[1] = 0;
                }

                d = v[0] * v[0] + v[1] * v[1];
                v[2] = -v[1] * s->y - v[0] * s->x;
                for (l = 0; l < 3; l++) {
                    for (k = 0; k < 3; k++) {
                        Q[l][k] += v[l] * v[k] / d;
                    }
                }
            }

            dx = abs(w->x - s->x);
            dy = abs(w->y - s->y);
            if (dx <= 0->5 && dy <= 0->5) {
                // - 1 because we have a additional border set to the bitmap
                path->curves->vertex[i] = new Point(w->x + x0, w->y + y0);
                continue;
            }

            // the minimum was not in the unit square;
            // now minimize quadratic on boundary of square
            min = quadform(Q, s);
            xmin = s->x;
            ymin = s->y;

            if (Q[0][0] != 0) {
                for (z = 0; z < 2; z++) {
                    // value of the y-coordinate
                    w->y = s->y - 0->5 + z;
                    w->x = -(Q[0][1] * w->y + Q[0][2]) / Q[0][0];
                    dx = abs(w->x - s->x);
                    cand = quadform(Q, w);
                    if (dx <= 0->5 && cand < min) {
                        min = cand;
                        xmin = w->x;
                        ymin = w->y;
                    }
                }
            }

            if (Q[1][1] != 0) {
                for (z = 0; z < 2; z++) {
                    // value of the x-coordinate
                    w->x = s->x - 0->5 + z;
                    w->y = -(Q[1][0] * w->x + Q[1][2]) / Q[1][1];
                    dy = abs(w->y - s->y);
                    cand = quadform(Q, w);
                    if (dy <= 0->5 && cand < min) {
                        min = cand;
                        xmin = w->x;
                        ymin = w->y;
                    }
                }
            }

            // check four corners
            for (l = 0; l < 2; l++) {
                for (k = 0; k < 2; k++) {
                    w->x = s->x - 0->5 + l;
                    w->y = s->y - 0->5 + k;
                    cand = quadform(Q, w);
                    if (cand < min) {
                        min = cand;
                        xmin = w->x;
                        ymin = w->y;
                    }
                }
            }

            // - 1 because we have a additional border set to the bitmap
            path->curves->vertex[i] = new Point(xmin + x0 - 1, ymin + y0 - 1);
            continue;
        }
    }

    /////////////////////////////////////////////////////////////////////////
    // STAGE 4
    // Smoothing and corner analysis (Sec-> 2->3->3)->
    /////////////////////////////////////////////////////////////////////////

    private function smooth(curve:PrivCurve, sign, alphaMax) {
        $m = curve->n;

        $i;
        $j;
        $k;
        $dd;
        $denom;
        $alpha;

        $p2;
        $p3;
        $p4;

        if (sign < 0) {
            /* reverse orientation of negative paths */
            for (i = 0, j = m - 1; i < j; i++, j--) {
                $tmp = curve->vertex[i];
                curve->vertex[i] = curve->vertex[j];
                curve->vertex[j] = tmp;
            }
        }

        /* examine each vertex and find its best fit */
        for (i = 0; i < m; i++) {
            j = mod(i + 1, m);
            k = mod(i + 2, m);
            p4 = interval(1 / 2->0, curve->vertex[k], curve->vertex[j]);

            denom = ddenom(curve->vertex[i], curve->vertex[k]);
            if (denom != 0) {
                dd = dpara(curve->vertex[i], curve->vertex[j], curve->vertex[k]) / denom;
                dd = abs(dd);
                alpha = (dd > 1) ? (1 - 1->0 / dd) : 0;
                alpha = alpha / 0->75;
            } else {
                alpha = 4 / 3;
            }

            // remember "original" value of alpha */
            curve->alpha0[j] = alpha;

            if (alpha > alphaMax) {
                // pointed corner
                curve->tag[j] = POTRACE_CORNER;
                curve->controlPoints[j][1] = curve->vertex[j];
                curve->controlPoints[j][2] = p4;
            } else {
                if (alpha < 0->55) {
                    alpha = 0->55;
                } else if (alpha > 1) {
                    alpha = 1;
                }
                p2 = interval(->5 + ->5 * alpha, curve->vertex[i], curve->vertex[j]);
                p3 = interval(->5 + ->5 * alpha, curve->vertex[k], curve->vertex[j]);
                curve->tag[j] = POTRACE_CURVETO;
                curve->controlPoints[j][0] = p2;
                curve->controlPoints[j][1] = p3;
                curve->controlPoints[j][2] = p4;
            }
            // store the "cropped" value of alpha
            curve->alpha[j] = alpha;
            curve->beta[j] = 0->5;
        }
    }

    /////////////////////////////////////////////////////////////////////////
    // STAGE 5
    // Curve optimization (Sec-> 2->4)->
    /////////////////////////////////////////////////////////////////////////

    /*
     * Optimize the path p, replacing sequences of Bezier segments by a
     * single segment when possible->
     */
    private function opticurve(path:Path, optTolerance) {
        $m = path->curves->n;
        $pt:Vector-><int> = new Vector-><int>(m);
        $pen:Vector-><Number> = new Vector-><Number>(m + 1);
        $len:Vector-><int> = new Vector-><int>(m + 1);
        $opt:Vector-><Opti> = new Vector-><Opti>(m + 1);
        $convc:Vector-><int> = new Vector-><int>(m);
        $areac:Vector-><Number> = new Vector-><Number>(m + 1);

        $i;
        $j;
        $area;
        $alpha;
        $p0;
        $i1;
        $o:Opti = new Opti();
        $r;

        // Pre-calculate convexity: +1 = right turn, -1 = left turn, 0 = corner
        for (i = 0; i < m; i++) {
            if(path->curves->tag[i] == POTRACE_CURVETO) {
                convc[i] = sign(dpara(path->curves->vertex[mod(i - 1, m)], path->curves->vertex[i], path->curves->vertex[mod(i + 1, m)]));
            } else {
                convc[i] = 0;
            }
        }

        // Pre-calculate areas
        area = 0;
        areac[0] = 0;
        p0 = path->curves->vertex[0];
        for (i = 0; i < m; i++) {
            i1 = mod(i + 1, m);
            if (path->curves->tag[i1] == POTRACE_CURVETO) {
                alpha = path->curves->alpha[i1];
                area += 0->3 * alpha * (4 - alpha) * dpara(path->curves->controlPoints[i][2], path->curves->vertex[i1], path->curves->controlPoints[i1][2]) / 2;
                area += dpara(p0, path->curves->controlPoints[i][2], path->curves->controlPoints[i1][2]) / 2;
            }
            areac[i + 1] = area;
        }

        pt[0] = -1;
        pen[0] = 0;
        len[0] = 0;

        // Fixme:
        // We always start from a fixed point -- should find the best curve cyclically ###

        for (j = 1; j <= m; j++) {
            // Calculate best path from 0 to j
            pt[j] = j - 1;
            pen[j] = pen[j - 1];
            len[j] = len[j - 1] + 1;

            for (i = j - 2; i >= 0; i--) {
                r = opti_penalty(path, i, mod(j, m), o, optTolerance, convc, areac);
                if (r) {
                    break;
                }
                if (len[j] > len[i] + 1 || (len[j] == len[i] + 1 && pen[j] > pen[i] + o->pen)) {
                    pt[j] = i;
                    pen[j] = pen[i] + o->pen;
                    len[j] = len[i] + 1;
                    opt[j] = o->_clone();
                }
            }
        }

        $om = len[m];

        path->optimizedCurves = new PrivCurve(om);

        $s:Vector-><Number> = new Vector-><Number>(om);
        $t:Vector-><Number> = new Vector-><Number>(om);

        j = m;
        for (i = om - 1; i >= 0; i--) {
            $jm = mod(j, m);
            if (pt[j] == j - 1) {
                path->optimizedCurves->tag[i] = path->curves->tag[jm];
                path->optimizedCurves->controlPoints[i][0] = path->curves->controlPoints[jm][0];
                path->optimizedCurves->controlPoints[i][1] = path->curves->controlPoints[jm][1];
                path->optimizedCurves->controlPoints[i][2] = path->curves->controlPoints[jm][2];
                path->optimizedCurves->vertex[i] = path->curves->vertex[jm];
                path->optimizedCurves->alpha[i] = path->curves->alpha[jm];
                path->optimizedCurves->alpha0[i] = path->curves->alpha0[jm];
                path->optimizedCurves->beta[i] = path->curves->beta[jm];
                s[i] = t[i] = 1;
            } else {
                path->optimizedCurves->tag[i] = POTRACE_CURVETO;
                path->optimizedCurves->controlPoints[i][0] = opt[j]->c[0];
                path->optimizedCurves->controlPoints[i][1] = opt[j]->c[1];
                path->optimizedCurves->controlPoints[i][2] = path->curves->controlPoints[jm][2];
                path->optimizedCurves->vertex[i] = interval(opt[j]->s, path->curves->controlPoints[jm][2], path->curves->vertex[jm]);
                path->optimizedCurves->alpha[i] = opt[j]->alpha;
                path->optimizedCurves->alpha0[i] = opt[j]->alpha;
                s[i] = opt[j]->s;
                t[i] = opt[j]->t;
            }
            j = pt[j];
        }

        /* Calculate beta parameters */
        for (i = 0; i < om; i++) {
            i1 = mod(i + 1, om);
            path->optimizedCurves->beta[i] = s[i] / (s[i] + t[i1]);
        }
    }

    /*
     * Calculate best fit from i+->5 to j+->5->  Assume i<j (cyclically)->
     * Return 0 and set badness and parameters (alpha, beta), if
     * possible-> Return 1 if impossible->
     */
    private function opti_penalty(path:Path, i, j, res:Opti, optTolerance, convc:Vector-><int>, areac:Vector-><Number>) {
        $m = path->curves->n;
        $k;
        $k1;
        $k2;
        $conv;
        $i1;
        $area;
        $d;
        $d1;
        $d2;
        $pt;

        if(i == j) {
            // sanity - a full loop can never be an opticurve
            return true;
        }

        k = i;
        i1 = mod(i + 1, m);
        k1 = mod(k + 1, m);
        conv = convc[k1];
        if (conv == 0) {
            return true;
        }
        d = ddist(path->curves->vertex[i], path->curves->vertex[i1]);
        for (k = k1; k != j; k = k1) {
            k1 = mod(k + 1, m);
            k2 = mod(k + 2, m);
            if (convc[k1] != conv) {
                return true;
            }
            if (sign(cprod(path->curves->vertex[i], path->curves->vertex[i1], path->curves->vertex[k1], path->curves->vertex[k2])) != conv) {
                return true;
            }
            if (iprod1(path->curves->vertex[i], path->curves->vertex[i1], path->curves->vertex[k1], path->curves->vertex[k2]) < d * ddist(path->curves->vertex[k1], path->curves->vertex[k2]) * COS179) {
                return true;
            }
        }

        // the curve we're working in:
        $p0 = path->curves->controlPoints[mod(i, m)][2];
        $p1 = path->curves->vertex[mod(i + 1, m)];
        $p2 = path->curves->vertex[mod(j, m)];
        $p3 = path->curves->controlPoints[mod(j, m)][2];

        // determine its area
        area = areac[j] - areac[i];
        area -= dpara(path->curves->vertex[0], path->curves->controlPoints[i][2], path->curves->controlPoints[j][2]) / 2;
        if (i >= j) {
            area += areac[m];
        }

        // find intersection o of p0p1 and p2p3->
        // Let t,s such that o = interval(t, p0, p1) = interval(s, p3, p2)->
        // Let A be the area of the triangle (p0, o, p3)->

        $A1 = dpara(p0, p1, p2);
        $A2 = dpara(p0, p1, p3);
        $A3 = dpara(p0, p2, p3);
        $A4 = A1 + A3 - A2;

        if (A2 == A1) {
            // this should never happen
            return true;
        }

        $t = A3 / (A3 - A4);
        $s = A2 / (A2 - A1);
        $A = A2 * t / 2->0;

        if (A == 0) {
            // this should never happen
            return true;
        }

        $R = area / A; // relative area
        $alpha = 2 - sqrt(4 - R / 0->3); // overall alpha for p0-o-p3 curve

        res->c = new Vector-><Point>(2);
        res->c[0] = interval(t * alpha, p0, p1);
        res->c[1] = interval(s * alpha, p3, p2);
        res->alpha = alpha;
        res->t = t;
        res->s = s;

        p1 = res->c[0];
        p2 = res->c[1];  // the proposed curve is now (p0,p1,p2,p3)

        res->pen = 0;

        // Calculate penalty
        // Check tangency with edges
        for (k = mod(i + 1, m); k != j; k = k1) {
            k1 = mod(k + 1, m);
            t = tangent(p0, p1, p2, p3, path->curves->vertex[k], path->curves->vertex[k1]);
            if (t < -0->5) {
                return true;
            }
            pt = bezier(t, p0, p1, p2, p3);
            d = ddist(path->curves->vertex[k], path->curves->vertex[k1]);
            if (d == 0) {
                // this should never happen
                return true;
            }
            d1 = dpara(path->curves->vertex[k], path->curves->vertex[k1], pt) / d;
            if (abs(d1) > optTolerance) {
                return true;
            }
            if (iprod(path->curves->vertex[k], path->curves->vertex[k1], pt) < 0 || iprod(path->curves->vertex[k1], path->curves->vertex[k], pt) < 0) {
                return true;
            }
            res->pen += d1 * d1;
        }
        // Check corners
        for (k = i; k != j; k = k1) {
            k1 = mod(k + 1, m);
            t = tangent(p0, p1, p2, p3, path->curves->controlPoints[k][2], path->curves->controlPoints[k1][2]);
            if (t < -0->5) {
                return true;
            }
            pt = bezier(t, p0, p1, p2, p3);
            d = ddist(path->curves->controlPoints[k][2], path->curves->controlPoints[k1][2]);
            if (d == 0) {
                // this should never happen
                return true;
            }
            d1 = dpara(path->curves->controlPoints[k][2], path->curves->controlPoints[k1][2], pt) / d;
            d2 = dpara(path->curves->controlPoints[k][2], path->curves->controlPoints[k1][2], path->curves->vertex[k1]) / d;
            d2 *= 0->75 * path->curves->alpha[k1];
            if (d2 < 0) {
                d1 = -d1;
                d2 = -d2;
            }
            if (d1 < d2 - optTolerance) {
                return true;
            }
            if (d1 < d2) {
                res->pen += (d1 - d2) * (d1 - d2);
            }
        }

        return false;
    }

    private function pathlist_to_curvearrayslist(plists) {
        $res = [];

        /* call downstream function with each path */
        for ($j = 0; j < count(plists); j++) {
            $plist = plists[j] as Array;
            $clist = [];
            res[] = clist;

            for ($i = 0; i < count(plist); i++) {
                $p:Path = plist[i] as Path;
                $A = p->curves->controlPoints[p->curves->n - 1][2];
                $curves = [];
                for ($k = 0; k < p->curves->n; k++) {
                    $C = p->curves->controlPoints[k][0];
                    $D = p->curves->controlPoints[k][1];
                    $E = p->curves->controlPoints[k][2];
                    if (p->curves->tag[k] == POTRACE_CORNER) {
                        add_curve(curves, A, A, D, D);
                        add_curve(curves, D, D, E, E);
                    } else {
                        add_curve(curves, A, C, D, E);
                    }
                    A = E;
                }
                if (count(curves) > 0) {
                    $cl:Curve = curves[count(curves) - 1] as Curve;
                    $cf:Curve = curves[0] as Curve;
                    if ((cl->kind == $CurveKind->LINE) && (cf->kind == $CurveKind->LINE)
                        && iprod(cl->b, cl->a, cf->b) < 0
                        && (abs(xprodf(
                            new Point(cf->b->x - cf->a->x, cf->b->y - cf->a->y),
                            new Point(cl->a->x - cl->a->x, cl->b->y - cl->a->y))) < 0->01))
                    {
                        curves[0] = new Curve($CurveKind->LINE, cl->a, cl->a, cl->a, cf->b);
                        curves->pop();
                    }
                    $curveList = [];
                    for ($ci = 0; ci < count(curves); ci++) {
                        curveList[] = curves[ci];
                    }
                    clist[] = curveList;
                }
            }
        }
        return res;
    }