/* Copyright (C) 2001-2021 Artifex Software, Inc.

   All Rights Reserved.

   This software is provided AS-IS with no warranty, either express or

   implied.

   This software is distributed under license and may not be copied,

   modified or distributed except as expressly authorized under the terms

   of the license contained in the file LICENSE in this distribution.

   Refer to licensing information at http://www.artifex.com or contact

   Artifex Software, Inc.,  1305 Grant Avenue - Suite 200, Novato,

   CA 94945, U.S.A., +1(415)492-9861, for further information.

*/

/* Rendering for non-mesh shadings */

#include "math_.h"

#include "memory_.h"

#include "gx.h"

#include "gserrors.h"

#include "gsmatrix.h"   /* for gscoord.h */

#include "gscoord.h"

#include "gspath.h"

#include "gsptype2.h"

#include "gxcspace.h"

#include "gxdcolor.h"

#include "gxfarith.h"

#include "gxfixed.h"

#include "gxgstate.h"

#include "gxpath.h"

#include "gxshade.h"

#include "gxdevcli.h"

#include "gxshade4.h"

#include "gsicc_cache.h"

/* ---------------- Function-based shading ---------------- */

typedef struct Fb_frame_s { /* A rudiment of old code. */

    gs_rect region;

    gs_client_color cc[4];  /* colors at 4 corners */

    int state;

} Fb_frame_t;

typedef struct Fb_fill_state_s {

    shading_fill_state_common;

    const gs_shading_Fb_t *psh;

    gs_matrix_fixed ptm;  /* parameter space -> device space */

    Fb_frame_t frame;

} Fb_fill_state_t;

/****** NEED GC DESCRIPTOR ******/

static inline void

make_other_poles(patch_curve_t curve[4])

{

    int i, j;

    for (i = 0; i < 4; i++) {

        j = (i + 1) % 4;

        curve[i].control[0].x = (curve[i].vertex.p.x * 2 + curve[j].vertex.p.x) / 3;

        curve[i].control[0].y = (curve[i].vertex.p.y * 2 + curve[j].vertex.p.y) / 3;

        curve[i].control[1].x = (curve[i].vertex.p.x + curve[j].vertex.p.x * 2) / 3;

        curve[i].control[1].y = (curve[i].vertex.p.y + curve[j].vertex.p.y * 2) / 3;

        curve[i].straight = true;

    }

}

/* Transform a point with a fixed-point result. */

static void

gs_point_transform2fixed_clamped(const gs_matrix_fixed * pmat,

                         double x, double y, gs_fixed_point * ppt)

{

    gs_point fpt;

    gs_point_transform(x, y, (const gs_matrix *)pmat, &fpt);

    ppt->x = clamp_coord(fpt.x);

    ppt->y = clamp_coord(fpt.y);

}

static int

Fb_fill_region(Fb_fill_state_t * pfs, const gs_fixed_rect *rect)

{

    patch_fill_state_t pfs1;

    patch_curve_t curve[4];

    Fb_frame_t * fp = &pfs->frame;

    int code;

    memcpy(&pfs1, (shading_fill_state_t *)pfs, sizeof(shading_fill_state_t));

    pfs1.Function = pfs->psh->params.Function;

    code = init_patch_fill_state(&pfs1);

    if (code < 0)

        return code;

    pfs1.maybe_self_intersecting = false;

    pfs1.n_color_args = 2;

    pfs1.rect = *rect;

    gs_point_transform2fixed(&pfs->ptm, fp->region.p.x, fp->region.p.y, &curve[0].vertex.p);

    gs_point_transform2fixed(&pfs->ptm, fp->region.q.x, fp->region.p.y, &curve[1].vertex.p);

    gs_point_transform2fixed(&pfs->ptm, fp->region.q.x, fp->region.q.y, &curve[2].vertex.p);

    gs_point_transform2fixed(&pfs->ptm, fp->region.p.x, fp->region.q.y, &curve[3].vertex.p);

    make_other_poles(curve);

    curve[0].vertex.cc[0] = fp->region.p.x;   curve[0].vertex.cc[1] = fp->region.p.y;

    curve[1].vertex.cc[0] = fp->region.q.x;   curve[1].vertex.cc[1] = fp->region.p.y;

    curve[2].vertex.cc[0] = fp->region.q.x;   curve[2].vertex.cc[1] = fp->region.q.y;

    curve[3].vertex.cc[0] = fp->region.p.x;   curve[3].vertex.cc[1] = fp->region.q.y;

    code = patch_fill(&pfs1, curve, NULL, NULL);

    if (term_patch_fill_state(&pfs1))

        return_error(gs_error_unregistered); /* Must not happen. */

    return code;

}

int

gs_shading_Fb_fill_rectangle(const gs_shading_t * psh0, const gs_rect * rect,

                             const gs_fixed_rect * rect_clip,

                             gx_device * dev, gs_gstate * pgs)

{

    const gs_shading_Fb_t * const psh = (const gs_shading_Fb_t *)psh0;

    gs_matrix save_ctm;

    int xi, yi, code;

    float x[2], y[2];

    Fb_fill_state_t state;

    code = shade_init_fill_state((shading_fill_state_t *) & state, psh0, dev, pgs);

    if (code < 0)

        return code;

    state.psh = psh;

    /****** HACK FOR FIXED-POINT MATRIX MULTIPLY ******/

    gs_currentmatrix((gs_gstate *) pgs, &save_ctm);

    gs_concat((gs_gstate *) pgs, &psh->params.Matrix);

    state.ptm = pgs->ctm;

    gs_setmatrix((gs_gstate *) pgs, &save_ctm);

    /* Compute the parameter X and Y ranges. */

    {

        gs_rect pbox;

        code = gs_bbox_transform_inverse(rect, &psh->params.Matrix, &pbox);

        if (code < 0)

            return code;

        x[0] = max(pbox.p.x, psh->params.Domain[0]);

        x[1] = min(pbox.q.x, psh->params.Domain[1]);

        y[0] = max(pbox.p.y, psh->params.Domain[2]);

        y[1] = min(pbox.q.y, psh->params.Domain[3]);

    }

    if (x[0] > x[1] || y[0] > y[1]) {

        /* The region is outside the shading area. */

        if (state.icclink != NULL) gsicc_release_link(state.icclink);

        return 0;

    }

    for (xi = 0; xi < 2; ++xi)

        for (yi = 0; yi < 2; ++yi) {

            float v[2];

            v[0] = x[xi], v[1] = y[yi];

            gs_function_evaluate(psh->params.Function, v,

                                 state.frame.cc[yi * 2 + xi].paint.values);

        }

    state.frame.region.p.x = x[0];

    state.frame.region.p.y = y[0];

    state.frame.region.q.x = x[1];

    state.frame.region.q.y = y[1];

    code = Fb_fill_region(&state, rect_clip);

    if (state.icclink != NULL) gsicc_release_link(state.icclink);

    return code;

}

/* ---------------- Axial shading ---------------- */

typedef struct A_fill_state_s {

    const gs_shading_A_t *psh;

    gs_point delta;

    double length;

    double t0, t1;

    double v0, v1, u0, u1;

} A_fill_state_t;

/****** NEED GC DESCRIPTOR ******/

/* Note t0 and t1 vary over [0..1], not the Domain. */

typedef struct

{

    patch_curve_t curve[4];

    gs_point corners[4];

} corners_and_curves;

/* Ghostscript cannot possibly render any patch whose bounds aren't

 * representable in fixed's. In fact, this is a larger limit than

 * we need. We notionally have an area defined by coordinates

 * that can be represented in fixed point with at least 1 bit to

 * spare.

 *

 * Any patch that lies completely outside this region can be clipped

 * away. Any patch that isn't representable by fixed points can be

 * subdivided into 4.

 *

 * This avoids us subdividing patches huge numbers of times because

 * one side is just outside the region we will accept.

 */

#define MIN_CLIP_LIMIT ((int)(fixed2int(min_fixed)/2))

#define MAX_CLIP_LIMIT ((int)(fixed2int(max_fixed)/2))

static int not_clipped_away(const gs_point *p)

{

    if (p[0].x < MIN_CLIP_LIMIT &&

        p[1].x < MIN_CLIP_LIMIT &&

        p[2].x < MIN_CLIP_LIMIT &&

        p[3].x < MIN_CLIP_LIMIT)

        return 0; /* Clipped away! */

    if (p[0].x > MAX_CLIP_LIMIT &&

        p[1].x > MAX_CLIP_LIMIT &&

        p[2].x > MAX_CLIP_LIMIT &&

        p[3].x > MAX_CLIP_LIMIT)

        return 0; /* Clipped away! */

    if (p[0].y < MIN_CLIP_LIMIT &&

        p[1].y < MIN_CLIP_LIMIT &&

        p[2].y < MIN_CLIP_LIMIT &&

        p[3].y < MIN_CLIP_LIMIT)

        return 0; /* Clipped away! */

    if (p[0].y > MAX_CLIP_LIMIT &&

        p[1].y > MAX_CLIP_LIMIT &&

        p[2].y > MAX_CLIP_LIMIT &&

        p[3].y > MAX_CLIP_LIMIT)

        return 0; /* Clipped away! */

    return 1;

}

#define f_fits_in_fixed(f) f_fits_in_bits(f, fixed_int_bits)

static int

A_fill_region_floats(patch_fill_state_t *pfs1, corners_and_curves *cc, int depth)

{

    corners_and_curves sub[4];

    int code;

    if (depth == 32)

        return gs_error_limitcheck;

    if (depth > 0 &&

        f_fits_in_fixed(cc->corners[0].x) &&

        f_fits_in_fixed(cc->corners[0].y) &&

        f_fits_in_fixed(cc->corners[1].x) &&

        f_fits_in_fixed(cc->corners[1].y) &&

        f_fits_in_fixed(cc->corners[2].x) &&

        f_fits_in_fixed(cc->corners[2].y) &&

        f_fits_in_fixed(cc->corners[3].x) &&

        f_fits_in_fixed(cc->corners[3].y))

    {

        cc->curve[0].vertex.p.x = float2fixed(cc->corners[0].x);

        cc->curve[0].vertex.p.y = float2fixed(cc->corners[0].y);

        cc->curve[1].vertex.p.x = float2fixed(cc->corners[1].x);

        cc->curve[1].vertex.p.y = float2fixed(cc->corners[1].y);

        cc->curve[2].vertex.p.x = float2fixed(cc->corners[2].x);

        cc->curve[2].vertex.p.y = float2fixed(cc->corners[2].y);

        cc->curve[3].vertex.p.x = float2fixed(cc->corners[3].x);

        cc->curve[3].vertex.p.y = float2fixed(cc->corners[3].y);

        cc->curve[0].vertex.cc[1] = cc->curve[1].vertex.cc[1] =

                                    cc->curve[2].vertex.cc[1] =

                                    cc->curve[3].vertex.cc[1] = 0;

        make_other_poles(cc->curve);

        return patch_fill(pfs1, cc->curve, NULL, NULL);

    }

    /* We have patches with corners:

     *  0  1

     *  3  2

     * We subdivide these into 4 smaller patches:

     *

     *  0   10   1     Where 0123 are corners

     *   [0]  [1]      [0][1][2][3] are patches.

     *  3   23   2

     *  0   10   1

     *   [3]  [2]

     *  3   23   2

     */

    sub[0].corners[0].x = cc->corners[0].x;

    sub[0].corners[0].y = cc->corners[0].y;

    sub[1].corners[1].x = cc->corners[1].x;

    sub[1].corners[1].y = cc->corners[1].y;

    sub[2].corners[2].x = cc->corners[2].x;

    sub[2].corners[2].y = cc->corners[2].y;

    sub[3].corners[3].x = cc->corners[3].x;

    sub[3].corners[3].y = cc->corners[3].y;

    sub[1].corners[0].x = sub[0].corners[1].x = (cc->corners[0].x + cc->corners[1].x)/2;

    sub[1].corners[0].y = sub[0].corners[1].y = (cc->corners[0].y + cc->corners[1].y)/2;

    sub[3].corners[2].x = sub[2].corners[3].x = (cc->corners[2].x + cc->corners[3].x)/2;

    sub[3].corners[2].y = sub[2].corners[3].y = (cc->corners[2].y + cc->corners[3].y)/2;

    sub[3].corners[0].x = sub[0].corners[3].x = (cc->corners[0].x + cc->corners[3].x)/2;

    sub[3].corners[0].y = sub[0].corners[3].y = (cc->corners[0].y + cc->corners[3].y)/2;

    sub[2].corners[1].x = sub[1].corners[2].x = (cc->corners[1].x + cc->corners[2].x)/2;

    sub[2].corners[1].y = sub[1].corners[2].y = (cc->corners[1].y + cc->corners[2].y)/2;

    sub[0].corners[2].x = sub[1].corners[3].x =

                          sub[2].corners[0].x =

                          sub[3].corners[1].x = (sub[0].corners[3].x + sub[1].corners[2].x)/2;

    sub[0].corners[2].y = sub[1].corners[3].y =

                          sub[2].corners[0].y =

                          sub[3].corners[1].y = (sub[0].corners[3].y + sub[1].corners[2].y)/2;

    sub[0].curve[0].vertex.cc[0] = sub[0].curve[3].vertex.cc[0] =

                                   sub[3].curve[0].vertex.cc[0] =

                                   sub[3].curve[3].vertex.cc[0] = cc->curve[0].vertex.cc[0];

    sub[1].curve[1].vertex.cc[0] = sub[1].curve[2].vertex.cc[0] =

                                   sub[2].curve[1].vertex.cc[0] =

                                   sub[2].curve[2].vertex.cc[0] = cc->curve[1].vertex.cc[0];

    sub[0].curve[1].vertex.cc[0] = sub[0].curve[2].vertex.cc[0] =

                                   sub[1].curve[0].vertex.cc[0] =

                                   sub[1].curve[3].vertex.cc[0] =

                                   sub[2].curve[0].vertex.cc[0] =

                                   sub[2].curve[3].vertex.cc[0] =

                                   sub[3].curve[1].vertex.cc[0] =

                                   sub[3].curve[2].vertex.cc[0] = (cc->curve[0].vertex.cc[0] + cc->curve[1].vertex.cc[0])/2;

    depth++;

    if (not_clipped_away(sub[0].corners)) {

        code = A_fill_region_floats(pfs1, &sub[0], depth);

        if (code < 0)

            return code;

    }

    if (not_clipped_away(sub[1].corners)) {

        code = A_fill_region_floats(pfs1, &sub[1], depth);

        if (code < 0)

            return code;

    }

    if (not_clipped_away(sub[2].corners)) {

        code = A_fill_region_floats(pfs1, &sub[2], depth);

        if (code < 0)

            return code;

    }

    if (not_clipped_away(sub[3].corners)) {

        code = A_fill_region_floats(pfs1, &sub[3], depth);

        if (code < 0)

            return code;

    }

    return 0;

}

static int

A_fill_region(A_fill_state_t * pfs, patch_fill_state_t *pfs1)

{

    const gs_shading_A_t * const psh = pfs->psh;

    double x0 = psh->params.Coords[0] + pfs->delta.x * pfs->v0;

    double y0 = psh->params.Coords[1] + pfs->delta.y * pfs->v0;

    double x1 = psh->params.Coords[0] + pfs->delta.x * pfs->v1;

    double y1 = psh->params.Coords[1] + pfs->delta.y * pfs->v1;

    double h0 = pfs->u0, h1 = pfs->u1;

    corners_and_curves cc;

    int code;

    double dx0 = pfs->delta.x * h0;

    double dy0 = pfs->delta.y * h0;

    double dx1 = pfs->delta.x * h1;

    double dy1 = pfs->delta.y * h1;

    cc.curve[0].vertex.cc[0] = pfs->t0; /* The element cc[1] is set to a dummy value against */

    cc.curve[1].vertex.cc[0] = pfs->t1; /* interrupts while an idle priocessing in gxshade.6.c .  */

    cc.curve[2].vertex.cc[0] = pfs->t1;

    cc.curve[3].vertex.cc[0] = pfs->t0;

    cc.curve[0].vertex.cc[1] = 0; /* The element cc[1] is set to a dummy value against */

    cc.curve[1].vertex.cc[1] = 0; /* interrupts while an idle priocessing in gxshade.6.c .  */

    cc.curve[2].vertex.cc[1] = 0;

    cc.curve[3].vertex.cc[1] = 0;

    cc.corners[0].x = x0 + dy0;

    cc.corners[0].y = y0 - dx0;

    cc.corners[1].x = x1 + dy0;

    cc.corners[1].y = y1 - dx0;

    cc.corners[2].x = x1 + dy1;

    cc.corners[2].y = y1 - dx1;

    cc.corners[3].x = x0 + dy1;

    cc.corners[3].y = y0 - dx1;

    code = gs_point_transform2fixed(&pfs1->pgs->ctm, cc.corners[0].x, cc.corners[0].y, &cc.curve[0].vertex.p);

    if (code < 0)

        goto fail;

    code = gs_point_transform2fixed(&pfs1->pgs->ctm, cc.corners[1].x, cc.corners[1].y, &cc.curve[1].vertex.p);

    if (code < 0)

        goto fail;

    code = gs_point_transform2fixed(&pfs1->pgs->ctm, cc.corners[2].x, cc.corners[2].y, &cc.curve[2].vertex.p);

    if (code < 0)

        goto fail;

    code = gs_point_transform2fixed(&pfs1->pgs->ctm, cc.corners[3].x, cc.corners[3].y, &cc.curve[3].vertex.p);

    if (code < 0)

        goto fail;

    make_other_poles(cc.curve);

    return patch_fill(pfs1, cc.curve, NULL, NULL);

fail:

    if (code != gs_error_limitcheck)

        return code;

    code = gs_point_transform(cc.corners[0].x, cc.corners[0].y, (const gs_matrix *)&pfs1->pgs->ctm, &cc.corners[0]);

    if (code < 0)

        return code;

    code = gs_point_transform(cc.corners[1].x, cc.corners[1].y, (const gs_matrix *)&pfs1->pgs->ctm, &cc.corners[1]);

    if (code < 0)

        return code;

    code = gs_point_transform(cc.corners[2].x, cc.corners[2].y, (const gs_matrix *)&pfs1->pgs->ctm, &cc.corners[2]);

    if (code < 0)

        return code;

    code = gs_point_transform(cc.corners[3].x, cc.corners[3].y, (const gs_matrix *)&pfs1->pgs->ctm, &cc.corners[3]);

    if (code < 0)

        return code;

    return A_fill_region_floats(pfs1, &cc, 0);

}

static inline int

gs_shading_A_fill_rectangle_aux(const gs_shading_t * psh0, const gs_rect * rect,

                            const gs_fixed_rect *clip_rect,

                            gx_device * dev, gs_gstate * pgs)

{

    const gs_shading_A_t *const psh = (const gs_shading_A_t *)psh0;

    gs_function_t * const pfn = psh->params.Function;

    gs_matrix cmat;

    gs_rect t_rect;

    A_fill_state_t state;

    patch_fill_state_t pfs1;

    float d0 = psh->params.Domain[0], d1 = psh->params.Domain[1];

    float dd = d1 - d0;

    double t0, t1;

    gs_point dist;

    int code;

    state.psh = psh;

    code = shade_init_fill_state((shading_fill_state_t *)&pfs1, psh0, dev, pgs);

    if (code < 0)

        return code;

    pfs1.Function = pfn;

    pfs1.rect = *clip_rect;

    code = init_patch_fill_state(&pfs1);

    if (code < 0)

        goto fail;

    pfs1.maybe_self_intersecting = false;

    pfs1.function_arg_shift = 1;

    /*

     * Compute the parameter range.  We construct a matrix in which

     * (0,0) corresponds to t = 0 and (0,1) corresponds to t = 1,

     * and use it to inverse-map the rectangle to be filled.

     */

    cmat.tx = psh->params.Coords[0];

    cmat.ty = psh->params.Coords[1];

    state.delta.x = psh->params.Coords[2] - psh->params.Coords[0];

    state.delta.y = psh->params.Coords[3] - psh->params.Coords[1];

    cmat.yx = state.delta.x;

    cmat.yy = state.delta.y;

    cmat.xx = cmat.yy;

    cmat.xy = -cmat.yx;

    code = gs_bbox_transform_inverse(rect, &cmat, &t_rect);

    if (code < 0) {

        code = 0; /* Swallow this silently */

        goto fail;

    }

    t0 = min(max(t_rect.p.y, 0), 1);

    t1 = max(min(t_rect.q.y, 1), 0);

    state.v0 = t0;

    state.v1 = t1;

    state.u0 = t_rect.p.x;

    state.u1 = t_rect.q.x;

    state.t0 = t0 * dd + d0;

    state.t1 = t1 * dd + d0;

    code = gs_distance_transform(state.delta.x, state.delta.y, &ctm_only(pgs),

                          &dist);

    if (code < 0)

        goto fail;

    state.length = hypot(dist.x, dist.y); /* device space line length */

    code = A_fill_region(&state, &pfs1);

    if (psh->params.Extend[0] && t0 > t_rect.p.y) {

        if (code < 0)

            goto fail;

        /* Use the general algorithm, because we need the trapping. */

        state.v0 = t_rect.p.y;

        state.v1 = t0;

        state.t0 = state.t1 = t0 * dd + d0;

        code = A_fill_region(&state, &pfs1);

    }

    if (psh->params.Extend[1] && t1 < t_rect.q.y) {

        if (code < 0)

            goto fail;

        /* Use the general algorithm, because we need the trapping. */

        state.v0 = t1;

        state.v1 = t_rect.q.y;

        state.t0 = state.t1 = t1 * dd + d0;

        code = A_fill_region(&state, &pfs1);

    }

fail:

    gsicc_release_link(pfs1.icclink);

    if (term_patch_fill_state(&pfs1))

        return_error(gs_error_unregistered); /* Must not happen. */

    return code;

}

int

gs_shading_A_fill_rectangle(const gs_shading_t * psh0, const gs_rect * rect,

                            const gs_fixed_rect * rect_clip,

                            gx_device * dev, gs_gstate * pgs)

{

    return gs_shading_A_fill_rectangle_aux(psh0, rect, rect_clip, dev, pgs);

}

/* ---------------- Radial shading ---------------- */

/* Some notes on what I have struggled to understand about the following

 * function. This function renders the 'tube' given by interpolating one

 * circle to another.

 *

 * The first circle is at (x0, y0) with radius r0, and has 'color' t0.

 * The other circle is at (x1, y1) with radius r1, and has 'color' t1.

 *

 * We perform this rendering by approximating each quadrant of the 'tube'

 * by a tensor patch. The tensor patch is formed by taking a curve along

 * 1/4 of the circumference of the first circle, a straight line to the

 * equivalent point on the circumference of the second circle, a curve

 * back along the circumference of the second circle, and then a straight

 * line back to where we started.

 *

 * There is additional logic in this function that forms the directions of

 * the curves differently for different quadrants. This is done to ensure

 * that we always paint 'around' the tube from the back towards the front,

 * so we don't get unexpected regions showing though. This is explained more

 * below.

 *

 * The original code here examined the position change between the two

 * circles dx and dy. Based upon this vector it would pick which quadrant/

 * tensor patch to draw first. It would draw the quadrants/tensor patches

 * in anticlockwise order. Presumably this was intended to be done so that

 * the 'top' quadrant would be drawn last.

 *

 * Unfortunately this did not always work; see bug 692513. If the quadrants

 * were rendered in the order 0,1,2,3, the rendering of 1 was leaving traces

 * on top of 0, which was unexpected.

 *

 * I have therefore altered the code slightly; rather than picking a start

 * quadrant and moving anticlockwise, we now draw the 'undermost' quadrant,

 * then the two adjacent quadrants, then the topmost quadrant.

 *

 * For the purposes of explanation, we shall label the octants as below:

 *

 *     \2|1/       and Quadrants as:       |

 *     3\|/0                            Q1 | Q0

 *    ---+---                          ----+----

 *     4/|\7                            Q2 | Q3

 *     /5|6\                               |

 *

 * We find (dx,dy), the difference between the centres of the circles.

 * We look to see which octant this falls in. Firstly, this tells us which

 * quadrant of the circle we need to draw first (Octant n, starts with

 * Quadrant floor(n/2)). Secondly, it tells us which direction to form the

 * tensor patch in; we always want to draw from the side 'closest' to

 * dx/dy to the side further away. This ensures that we don't overwrite

 * pixels in the incorrect order as the patch decomposes.

 */

static int

R_tensor_annulus(patch_fill_state_t *pfs,

    double x0, double y0, double r0, double t0,

    double x1, double y1, double r1, double t1)

{

    double dx = x1 - x0, dy = y1 - y0;

    double d = hypot(dx, dy);

    gs_point p0, p1, pc0, pc1;

    int k, j, code, dirn;

    bool inside = 0;

    /* pc0 and pc1 are the centres of the respective circles. */

    pc0.x = x0, pc0.y = y0;

    pc1.x = x1, pc1.y = y1;

    /* Set p0 up so it's a unit vector giving the direction of 90 degrees

     * to the right of the major axis as we move from p0c to p1c. */

    if (r0 + d <= r1 || r1 + d <= r0) {

        /* One circle is inside another one.

           Use any subdivision,

           but don't depend on dx, dy, which may be too small. */

        p0.x = 0, p0.y = -1, dirn = 0;

        /* Align stripes along radii for faster triangulation : */

        inside = 1;

        pfs->function_arg_shift = 1;

    } else {

        /* Must generate canonic quadrangle arcs,

           because we approximate them with curves. */

        if(dx >= 0) {

            if (dy >= 0)

                p0.x = 1, p0.y = 0, dirn = (dx >= dy ? 1 : 0);

            else

                p0.x = 0, p0.y = -1, dirn = (dx >= -dy ? 0 : 1);

        } else {

            if (dy >= 0)

                p0.x = 0, p0.y = 1, dirn = (-dx >= dy ? 1 : 0);

            else

                p0.x = -1, p0.y = 0, dirn = (-dx >= -dy ? 0 : 1);

        }

        pfs->function_arg_shift = 0;

    }

    /* fixme: wish: cut invisible parts off.

       Note : when r0 != r1 the invisible part is not a half circle. */

    for (k = 0; k < 4; k++) {

        gs_point p[12];

        patch_curve_t curve[4];

        /* Set p1 to be 90 degrees anticlockwise from p0 */

        p1.x = -p0.y; p1.y = p0.x;

        if (dirn == 0) { /* Clockwise */

            make_quadrant_arc(p + 0, &pc0, &p1, &p0, r0);

            make_quadrant_arc(p + 6, &pc1, &p0, &p1, r1);

        } else { /* Anticlockwise */

            make_quadrant_arc(p + 0, &pc0, &p0, &p1, r0);

            make_quadrant_arc(p + 6, &pc1, &p1, &p0, r1);

        }

        p[4].x = (p[3].x * 2 + p[6].x) / 3;

        p[4].y = (p[3].y * 2 + p[6].y) / 3;

        p[5].x = (p[3].x + p[6].x * 2) / 3;

        p[5].y = (p[3].y + p[6].y * 2) / 3;

        p[10].x = (p[9].x * 2 + p[0].x) / 3;

        p[10].y = (p[9].y * 2 + p[0].y) / 3;

        p[11].x = (p[9].x + p[0].x * 2) / 3;

        p[11].y = (p[9].y + p[0].y * 2) / 3;

        for (j = 0; j < 4; j++) {

            int jj = (j + inside) % 4;

            if (gs_point_transform2fixed(&pfs->pgs->ctm,         p[j*3 + 0].x, p[j*3 + 0].y, &curve[jj].vertex.p) < 0)

                gs_point_transform2fixed_clamped(&pfs->pgs->ctm, p[j*3 + 0].x, p[j*3 + 0].y, &curve[jj].vertex.p);

            if (gs_point_transform2fixed(&pfs->pgs->ctm,         p[j*3 + 1].x, p[j*3 + 1].y, &curve[jj].control[0]) < 0)

                gs_point_transform2fixed_clamped(&pfs->pgs->ctm, p[j*3 + 1].x, p[j*3 + 1].y, &curve[jj].control[0]);

            if (gs_point_transform2fixed(&pfs->pgs->ctm,         p[j*3 + 2].x, p[j*3 + 2].y, &curve[jj].control[1]) < 0)

                gs_point_transform2fixed_clamped(&pfs->pgs->ctm, p[j*3 + 2].x, p[j*3 + 2].y, &curve[jj].control[1]);

            curve[j].straight = (((j + inside) & 1) != 0);

        }

        curve[(0 + inside) % 4].vertex.cc[0] = t0;

        curve[(1 + inside) % 4].vertex.cc[0] = t0;

        curve[(2 + inside) % 4].vertex.cc[0] = t1;

        curve[(3 + inside) % 4].vertex.cc[0] = t1;

        curve[0].vertex.cc[1] = curve[1].vertex.cc[1] = 0; /* Initialize against FPE. */

        curve[2].vertex.cc[1] = curve[3].vertex.cc[1] = 0; /* Initialize against FPE. */

        code = patch_fill(pfs, curve, NULL, NULL);

        if (code < 0)

            return code;

        /* Move p0 to be ready for the next position */

        if (k == 0) {

            /* p0 moves clockwise */

            p1 = p0;

            p0.x = p1.y; p0.y = -p1.x;

            dirn = 0;

        } else if (k == 1) {

            /* p0 flips sides */

            p0.x = -p0.x; p0.y = -p0.y;

            dirn = 1;

        } else if (k == 2) {

            /* p0 moves anti-clockwise */

            p1 = p0;

            p0.x = -p1.y; p0.y = p1.x;

            dirn = 0;

        }

    }

    return 0;

}

/* Find the control points for two points on the arc of a circle

 * the points must be within the same quadrant.

 */

static int find_arc_control_points(gs_point *from, gs_point *to, gs_point *from_control, gs_point *to_control, gs_point *centre)

{

    double from_tan_alpha, to_tan_alpha, from_alpha, to_alpha;

    double half_inscribed_angle, intersect_x, intersect_y, intersect_dist;

    double radius = sqrt(((from->x - centre->x) * (from->x - centre->x)) + ((from->y - centre->y) * (from->y - centre->y)));

    double tangent_intersect_dist;

    double F;

    int quadrant;

    /* Quadrant 0 is upper right, numbered anti-clockwise.

     * If the direction of the from->to is atni-clockwise, add 4

     */

    if (from->x > to->x) {

        if (from->y > to->y) {

            if (to->y >= centre->y)

                quadrant = 1 + 4;

            else

                quadrant = 3;

        } else {

            if (to->x >= centre->x)

                quadrant = 0 + 4;

            else

                quadrant = 2;

        }

    } else {

        if (from->y > to->y) {

            if (from->x >= centre->x)

                quadrant = 0;

            else

                quadrant = 2 + 4;

        } else {

            if (from->x >= centre->x)

                quadrant = 3 + 4;

            else

                quadrant = 1;

        }

    }

    switch(quadrant) {

        /* quadrant 0, arc goes clockwise */

        case 0:

            if (from->x == centre->x) {

                from_alpha = M_PI / 2;

            } else {

                from_tan_alpha = (from->y - centre->y) / (from->x - centre->x);

                from_alpha = atan(from_tan_alpha);

            }

            to_tan_alpha = (to->y - centre->y) / (to->x - centre->x);

            to_alpha = atan(to_tan_alpha);

            half_inscribed_angle = (from_alpha - to_alpha) / 2;

            intersect_dist = radius / cos(half_inscribed_angle);

            tangent_intersect_dist = tan(half_inscribed_angle) * radius;

            intersect_x = centre->x + cos(to_alpha + half_inscribed_angle) * intersect_dist;

            intersect_y = centre->y + sin(to_alpha + half_inscribed_angle) * intersect_dist;

            break;

        /* quadrant 1, arc goes clockwise */

        case 1:

            from_tan_alpha = (from->y - centre->y) / (centre->x - from->x);

            from_alpha = atan(from_tan_alpha);

            if (to->x == centre->x) {

                to_alpha = M_PI / 2;

            } else {

                to_tan_alpha = (to->y - centre->y) / (centre->x - to->x);

                to_alpha = atan(to_tan_alpha);

            }

            half_inscribed_angle = (to_alpha - from_alpha) / 2;

            intersect_dist = radius / cos(half_inscribed_angle);

            tangent_intersect_dist = tan(half_inscribed_angle) * radius;

            intersect_x = centre->x - cos(from_alpha + half_inscribed_angle) * intersect_dist;

            intersect_y = centre->y + sin(from_alpha + half_inscribed_angle) * intersect_dist;

            break;

        /* quadrant 2, arc goes clockwise */

        case 2:

            if (from->x == centre->x) {

                from_alpha = M_PI / 2;

            } else {

                from_tan_alpha = (centre->y - from->y) / (centre->x - from->x);

                from_alpha = atan(from_tan_alpha);

            }

            to_tan_alpha = (centre->y - to->y) / (centre->x - to->x);

            to_alpha = atan(to_tan_alpha);

            half_inscribed_angle = (to_alpha - from_alpha) / 2;

            intersect_dist = radius / cos(half_inscribed_angle);

            tangent_intersect_dist = tan(half_inscribed_angle) * radius;

            intersect_x = centre->x - cos(from_alpha + half_inscribed_angle) * intersect_dist;

            intersect_y = centre->y - sin(from_alpha + half_inscribed_angle) * intersect_dist;

            break;

        /* quadrant 3, arc goes clockwise */

        case 3:

            from_tan_alpha = (centre->y - from->y) / (from->x - centre->x);

            from_alpha = atan(from_tan_alpha);

            if (to->x == centre->x) {

                to_alpha = M_PI / 2;

            } else {

                to_tan_alpha = (centre->y - to->y) / (to->x - centre->x);

                to_alpha = atan(to_tan_alpha);

            }

            half_inscribed_angle = (to_alpha - from_alpha) / 2;

            intersect_dist = radius / cos(half_inscribed_angle);

            tangent_intersect_dist = tan(half_inscribed_angle) * radius;

            intersect_x = centre->x + cos(from_alpha + half_inscribed_angle) * intersect_dist;

            intersect_y = centre->y - sin(from_alpha + half_inscribed_angle) * intersect_dist;

            break;

        /* quadrant 0, arc goes anti-clockwise */

        case 4:

            from_tan_alpha = (from->y - centre->y) / (from->x - centre->x);

            from_alpha = atan(from_tan_alpha);

            if (to->y == centre->y)

                to_alpha = M_PI / 2;

            else {

                to_tan_alpha = (to->y - centre->y) / (to->x - centre->x);

                to_alpha = atan(to_tan_alpha);

            }

            half_inscribed_angle = (to_alpha - from_alpha) / 2;

            intersect_dist = radius / cos(half_inscribed_angle);

            tangent_intersect_dist = tan(half_inscribed_angle) * radius;

            intersect_x = centre->x + cos(from_alpha + half_inscribed_angle) * intersect_dist;

            intersect_y = centre->y + sin(from_alpha + half_inscribed_angle) * intersect_dist;

            break;

        /* quadrant 1, arc goes anti-clockwise */

        case 5:

            from_tan_alpha = (centre->x - from->x) / (from->y - centre->y);

            from_alpha = atan(from_tan_alpha);

            if (to->y == centre->y) {

                to_alpha = M_PI / 2;

            }

            else {

                to_tan_alpha = (centre->x - to->x) / (to->y - centre->y);

                to_alpha = atan(to_tan_alpha);

            }

            half_inscribed_angle = (to_alpha - from_alpha) / 2;

            intersect_dist = radius / cos(half_inscribed_angle);

            tangent_intersect_dist = tan(half_inscribed_angle) * radius;

            intersect_x = centre->x - sin(from_alpha + half_inscribed_angle) * intersect_dist;

            intersect_y = centre->y + cos(from_alpha + half_inscribed_angle) * intersect_dist;

            break;

        /* quadrant 2, arc goes anti-clockwise */

        case 6:

            from_tan_alpha = (from->y - centre->y) / (centre->x - from->x);

            from_alpha = atan(from_tan_alpha);

            if (to->x == centre->x) {

                to_alpha = M_PI / 2;

            } else {

                to_tan_alpha = (centre->y - to->y) / (centre->x - to->x);

                to_alpha = atan(to_tan_alpha);

            }

            half_inscribed_angle = (to_alpha - from_alpha) / 2;

            intersect_dist = radius / cos(half_inscribed_angle);

            tangent_intersect_dist = tan(half_inscribed_angle) * radius;

            intersect_x = centre->x - cos(from_alpha + half_inscribed_angle) * intersect_dist;

            intersect_y = centre->y - sin(from_alpha + half_inscribed_angle) * intersect_dist;

            break;

        /* quadrant 3, arc goes anti-clockwise */

        case 7:

            if (from->x == centre->x) {

                from_alpha = M_PI / 2;

            } else {

                from_tan_alpha = (centre->y - from->y) / (from->x - centre->x);

                from_alpha = atan(from_tan_alpha);

            }

            to_tan_alpha = (centre->y - to->y) / (to->x - centre->x);

            to_alpha = atan(to_tan_alpha);

            half_inscribed_angle = (from_alpha - to_alpha) / 2;

            intersect_dist = radius / cos(half_inscribed_angle);

            tangent_intersect_dist = tan(half_inscribed_angle) * radius;

            intersect_x = centre->x + cos(to_alpha + half_inscribed_angle) * intersect_dist;

            intersect_y = centre->y - sin(to_alpha + half_inscribed_angle) * intersect_dist;

            break;

    }

    F = (4.0 / 3.0) / (1 + sqrt(1 + ((tangent_intersect_dist / radius) * (tangent_intersect_dist / radius))));

    from_control->x = from->x - ((from->x - intersect_x) * F);

    from_control->y = from->y - ((from->y - intersect_y) * F);

    to_control->x = to->x - ((to->x - intersect_x) * F);

    to_control->y = to->y - ((to->y - intersect_y) * F);

    return 0;

}

/* Create a 'patch_curve' element whch is a straight line between two points */

static int patch_lineto(gs_matrix_fixed *ctm, gs_point *from, gs_point *to, patch_curve_t *p, float t)

{

    double x_1third, x_2third, y_1third, y_2third;

    x_1third = (to->x - from->x) / 3;

    x_2third = x_1third * 2;

    y_1third = (to->y - from->y) / 3;

    y_2third = y_1third * 2;

    gs_point_transform2fixed(ctm, from->x, from->y, &p->vertex.p);

    gs_point_transform2fixed(ctm, from->x + x_1third, from->y + y_1third, &p->control[0]);

    gs_point_transform2fixed(ctm, from->x + x_2third, from->y + y_2third, &p->control[1]);

    p->vertex.cc[0] = t;

    p->vertex.cc[1] = t;

    p->straight = 1;

    return 0;

}

static int patch_curveto(gs_matrix_fixed *ctm, gs_point *centre, gs_point *from, gs_point *to, patch_curve_t *p, float t)

{

    gs_point from_control, to_control;

    find_arc_control_points(from, to, &from_control, &to_control, centre);

    gs_point_transform2fixed(ctm, from->x, from->y, &p->vertex.p);

    gs_point_transform2fixed(ctm, from_control.x, from_control.y, &p->control[0]);

    gs_point_transform2fixed(ctm, to_control.x, to_control.y, &p->control[1]);

    p->vertex.cc[0] = t;

    p->vertex.cc[1] = t;

    p->straight = 0;

    return 0;

}

static int draw_quarter_annulus(patch_fill_state_t *pfs, gs_point *centre, double radius, gs_point *corner, float t)

{

    gs_point p0, p1, initial;

    patch_curve_t p[4];

    int code;

    if (corner->x > centre->x) {

        initial.x = centre->x + radius;

    }

    else {

        initial.x = centre->x - radius;

    }

    initial.y = centre->y;

    p1.x = initial.x;

    p1.y = corner->y;

    patch_lineto(&pfs->pgs->ctm, &initial, &p1, &p[0], t);

    p0.x = centre->x;

    p0.y = p1.y;

    patch_lineto(&pfs->pgs->ctm, &p1, &p0, &p[1], t);

    p1.x = centre->x;

    if (centre->y > corner->y) {

        p1.y = centre->y - radius;

    } else {

        p1.y = centre->y + radius;

    }

    patch_lineto(&pfs->pgs->ctm, &p0, &p1, &p[2], t);

    patch_curveto(&pfs->pgs->ctm, centre, &p1, &initial, &p[3], t);

    code = patch_fill(pfs, (const patch_curve_t *)&p, NULL, NULL);

    if (code < 0)

        return code;

    if (corner->x > centre->x)

        initial.x = corner->x - (corner->x - (centre->x + radius));

    else

        initial.x = centre->x - radius;

    initial.y = corner->y;

    patch_lineto(&pfs->pgs->ctm, corner, &initial, &p[0], t);

    p0.x = initial.x;

    p0.y = centre->y;

    patch_lineto(&pfs->pgs->ctm, &initial, &p0, &p[1], t);

    p1.y = p0.y;

    p1.x = corner->x;

    patch_lineto(&pfs->pgs->ctm, &p0, &p1, &p[2], t);

    patch_lineto(&pfs->pgs->ctm, &p1, corner, &p[3], t);

    return (patch_fill(pfs, (const patch_curve_t *)&p, NULL, NULL));

}

static int R_tensor_annulus_extend_tangent(patch_fill_state_t *pfs,

    double x0, double y0, double r0, double t0,

    double x1, double y1, double r1, double t1, double r2)

{

    patch_curve_t curve[4];

    gs_point p0, p1;

    int code = 0, q = 0;

    /* special case axis aligned circles. Its quicker to handle these specially as it

     * avoid lots of trigonometry in the general case code, and avoids us

     * having to watch out for infinity as the result of tan() operations.

     */

    if (x0 == x1 || y0 == y1) {

        if (x0 == x1 && y0 > y1) {

            /* tangent at top of circles */