/* Copyright (C) 2001-2024 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.,  39 Mesa Street, Suite 108A, San Francisco,

   CA 94129, USA, for further information.

*/

/* Rendering for Coons and tensor patch shadings */

/*

   A contiguous non-overlapping decomposition

   of a tensor shading into linear color triangles.

 */

#include "memory_.h"

#include "gx.h"

#include "gserrors.h"

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

#include "gscoord.h"

#include "gscicach.h"

#include "gxcspace.h"

#include "gxdcolor.h"

#include "gxgstate.h"

#include "gxshade.h"

#include "gxdevcli.h"

#include "gxshade4.h"

#include "gxarith.h"

#include "gzpath.h"

#include "stdint_.h"

#include "math_.h"

#include "gsicc_cache.h"

#include "gxdevsop.h"

/* The original version of the shading code 'decompose's shadings into

 * smaller and smaller regions until they are smaller than 1 pixel, and then

 * fills them. (Either with a constant colour, or with a linear filled trap).

 *

 * Previous versions of the code (from svn revision 7936 until June 2011

 * (shortly after the switch to git)) changed this limit to be 1 point or 1

 * pixel (whichever is larger) (on the grounds that as resolution increases,

 * we are unlikely to be able to notice increasingly small inaccuracies in

 * the shading). Given how people abuse the resolution at which things are

 * rendered (especially when rendering to images that can subsequently be

 * zoomed), this seems a poor trade off. See bug 691152.

 *

 * The code has now been changed back to operate with the proper 1 pixel

 * decomposition, which will cost us performance in some cases. Should

 * people want to restore the previous operation, they should build with

 * MAX_SHADING_RESOLUTION predefined to 72. In general, this symbol can be

 * set to the maximum dpi that shading should ever be performed at. Leaving

 * it undefined will leave the default (1 pixel limit) in place.

 *

 * A possible future optimisation here may be to use different max shading

 * resolutions for linear and non-linear shadings; linear shadings appear to

 * result in calls to "fill linear traps", and non linear ones appear to

 * result in calls to "fill constant color". As such linear shadings are much

 * more forgiving of a higher decomposition threshold.

 */

#if NOFILL_TEST

static bool dbg_nofill = false;

#endif

#if SKIP_TEST

static int dbg_patch_cnt = 0;

static int dbg_quad_cnt = 0;

static int dbg_triangle_cnt = 0;

static int dbg_wedge_triangle_cnt = 0;

#endif

enum {

    min_linear_grades = 255 /* The minimal number of device color grades,

            required to apply linear color device functions. */

};

/* ================ Utilities ================ */

static int

allocate_color_stack(patch_fill_state_t *pfs, gs_memory_t *memory)

{

    if (pfs->color_stack != NULL)

        return 0;

    pfs->color_stack_step = offset_of(patch_color_t, cc.paint.values[pfs->num_components]);

    pfs->color_stack_step = (pfs->color_stack_step + sizeof(void *) - 1) / sizeof(void *) * sizeof(void *); /* Alignment */

    pfs->color_stack_size = pfs->color_stack_step * SHADING_COLOR_STACK_SIZE;

    pfs->color_stack = gs_alloc_bytes(memory, pfs->color_stack_size, "allocate_color_stack");

    if (pfs->color_stack == NULL)

        return_error(gs_error_VMerror);

    pfs->color_stack_limit = pfs->color_stack + pfs->color_stack_size;

    pfs->color_stack_ptr = pfs->color_stack;

    pfs->memory = memory;

    return 0;

}

static inline byte *

reserve_colors_inline(patch_fill_state_t *pfs, patch_color_t *c[], int n)

{

    int i;

    byte *ptr0 = pfs->color_stack_ptr, *ptr = ptr0;

    for (i = 0; i < n; i++, ptr += pfs->color_stack_step)

        c[i] = (patch_color_t *)ptr;

    if (ptr > pfs->color_stack_limit) {

        c[0] = NULL; /* safety. */

        return NULL;

    }

    pfs->color_stack_ptr = ptr;

    return ptr0;

}

byte *

reserve_colors(patch_fill_state_t *pfs, patch_color_t *c[], int n)

{

    return reserve_colors_inline(pfs, c, n);

}

static inline void

release_colors_inline(patch_fill_state_t *pfs, byte *ptr, int n)

{

#if 0 /* Saving the invariant for records. */

    pfs->color_stack_ptr = pfs->color_stack_step * n;

    assert((byte *)pfs->color_stack_ptr == ptr);

#else

    pfs->color_stack_ptr = ptr;

#endif

}

void

release_colors(patch_fill_state_t *pfs, byte *ptr, int n)

{

    release_colors_inline(pfs, ptr, n);

}

/* Get colors for patch vertices. */

static int

shade_next_colors(shade_coord_stream_t * cs, patch_curve_t * curves,

                  int num_vertices)

{

    int i, code = 0;

    for (i = 0; i < num_vertices && code >= 0; ++i) {

        curves[i].vertex.cc[1] = 0; /* safety. (patch_fill may assume 2 arguments) */

        code = shade_next_color(cs, curves[i].vertex.cc);

    }

    return code;

}

/* Get a Bezier or tensor patch element. */

static int

shade_next_curve(shade_coord_stream_t * cs, patch_curve_t * curve)

{

    int code = shade_next_coords(cs, &curve->vertex.p, 1);

    if (code >= 0)

        code = shade_next_coords(cs, curve->control,

                                 countof(curve->control));

    return code;

}

/*

 * Parse the next patch out of the input stream.  Return 1 if done,

 * 0 if patch, <0 on error.

 */

static int

shade_next_patch(shade_coord_stream_t * cs, int BitsPerFlag,

        patch_curve_t curve[4], gs_fixed_point interior[4] /* 0 for Coons patch */)

{

    int flag = shade_next_flag(cs, BitsPerFlag);

    int num_colors, code;

    if (flag < 0) {

        if (!cs->is_eod(cs))

            return_error(gs_error_rangecheck);

        return 1;               /* no more data */

    }

    if (cs->first_patch && (flag & 3) != 0) {

        return_error(gs_error_rangecheck);

    }

    cs->first_patch = 0;

    switch (flag & 3) {

        default:

            return_error(gs_error_rangecheck);  /* not possible */

        case 0:

            if ((code = shade_next_curve(cs, &curve[0])) < 0 ||

                (code = shade_next_coords(cs, &curve[1].vertex.p, 1)) < 0

                )

                return code;

            num_colors = 4;

            goto vx;

        case 1:

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

            goto v3;

        case 2:

            curve[0] = curve[2], curve[1].vertex = curve[3].vertex;

            goto v3;

        case 3:

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

v3:         num_colors = 2;

vx:         if ((code = shade_next_coords(cs, curve[1].control, 2)) < 0 ||

                (code = shade_next_curve(cs, &curve[2])) < 0 ||

                (code = shade_next_curve(cs, &curve[3])) < 0 ||

                (interior != 0 &&

                 (code = shade_next_coords(cs, interior, 4)) < 0) ||

                (code = shade_next_colors(cs, &curve[4 - num_colors],

                                          num_colors)) < 0

                )

                return code;

            cs->align(cs, 8); /* See shade_next_vertex. */

    }

    return 0;

}

static inline bool

is_linear_color_applicable(const patch_fill_state_t *pfs)

{

    if (!USE_LINEAR_COLOR_PROCS)

        return false;

    if (!colors_are_separable_and_linear(&pfs->dev->color_info))

        return false;

    if (gx_get_cmap_procs(pfs->pgs, pfs->dev)->is_halftoned(pfs->pgs, pfs->dev))

        return false;

    return true;

}

static int

alloc_patch_fill_memory(patch_fill_state_t *pfs, gs_memory_t *memory, const gs_color_space *pcs)

{

    int code;

    pfs->memory = memory;

#   if LAZY_WEDGES

        code = wedge_vertex_list_elem_buffer_alloc(pfs);

        if (code < 0)

            return code;

#   endif

    pfs->max_small_coord = 1 << ((sizeof(int64_t) * 8 - 1/*sign*/) / 3);

    code = allocate_color_stack(pfs, memory);

    if (code < 0)

        return code;

    if (pfs->unlinear || pcs == NULL)

        pfs->pcic = NULL;

    else {

        pfs->pcic = gs_color_index_cache_create(memory, pcs, pfs->dev, pfs->pgs, true, pfs->trans_device);

        if (pfs->pcic == NULL)

            return_error(gs_error_VMerror);

    }

    return 0;

}

int

init_patch_fill_state(patch_fill_state_t *pfs)

{

    /* Warning : pfs->Function must be set in advance. */

    const gs_color_space *pcs = pfs->direct_space;

    gs_client_color fcc0, fcc1;

    int i;

    for (i = 0; i < pfs->num_components; i++) {

        fcc0.paint.values[i] = -1000000;

        fcc1.paint.values[i] = 1000000;

    }

    pcs->type->restrict_color(&fcc0, pcs);

    pcs->type->restrict_color(&fcc1, pcs);

    for (i = 0; i < pfs->num_components; i++)

        pfs->color_domain.paint.values[i] = max(fcc1.paint.values[i] - fcc0.paint.values[i], 1);

    pfs->vectorization = false; /* A stub for a while. Will use with pclwrite. */

    pfs->maybe_self_intersecting = true;

    pfs->monotonic_color = (pfs->Function == NULL);

    pfs->function_arg_shift = 0;

    pfs->linear_color = false;

    pfs->inside = false;

    pfs->n_color_args = 1;

#ifdef MAX_SHADING_RESOLUTION

    pfs->decomposition_limit = float2fixed(min(pfs->dev->HWResolution[0],

                                               pfs->dev->HWResolution[1]) / MAX_SHADING_RESOLUTION);

    pfs->decomposition_limit = max(pfs->decomposition_limit, fixed_1);

#else

    pfs->decomposition_limit = fixed_1;

#endif

    pfs->fixed_flat = float2fixed(pfs->pgs->flatness);

    /* Restrict the pfs->smoothness with 1/min_linear_grades, because cs_is_linear

       can't provide a better precision due to the color

       representation with integers.

     */

    pfs->smoothness = max(pfs->pgs->smoothness, 1.0 / min_linear_grades);

    pfs->color_stack_size = 0;

    pfs->color_stack_step = 0;

    pfs->color_stack_ptr = NULL;

    pfs->color_stack = NULL;

    pfs->color_stack_limit = NULL;

    pfs->unlinear = !is_linear_color_applicable(pfs);

    pfs->pcic = NULL;

    return alloc_patch_fill_memory(pfs, pfs->pgs->memory, pcs);

}

bool

term_patch_fill_state(patch_fill_state_t *pfs)

{

    bool b = (pfs->color_stack_ptr != pfs->color_stack);

#   if LAZY_WEDGES

        wedge_vertex_list_elem_buffer_free(pfs);

#   endif

    if (pfs->color_stack)

        gs_free_object(pfs->memory, pfs->color_stack, "term_patch_fill_state");

    if (pfs->pcic != NULL)

        gs_color_index_cache_destroy(pfs->pcic);

    return b;

}

/* Resolve a patch color using the Function if necessary. */

static inline void

patch_resolve_color_inline(patch_color_t * ppcr, const patch_fill_state_t *pfs)

{

    if (pfs->Function) {

        const gs_color_space *pcs = pfs->direct_space;

        gs_function_evaluate(pfs->Function, ppcr->t, ppcr->cc.paint.values);

        pcs->type->restrict_color(&ppcr->cc, pcs);

    }

}

void

patch_resolve_color(patch_color_t * ppcr, const patch_fill_state_t *pfs)

{

    patch_resolve_color_inline(ppcr, pfs);

}

/*

 * Calculate the interpolated color at a given point.

 * Note that we must do this twice for bilinear interpolation.

 * We use the name ppcr rather than ppc because an Apple compiler defines

 * ppc when compiling on PowerPC systems (!).

 */

static void

patch_interpolate_color(patch_color_t * ppcr, const patch_color_t * ppc0,

       const patch_color_t * ppc1, const patch_fill_state_t *pfs, double t)

{

    /* The old code gives -IND on Intel. */

    if (pfs->Function) {

        ppcr->t[0] = ppc0->t[0] * (1 - t) + t * ppc1->t[0];

        ppcr->t[1] = ppc0->t[1] * (1 - t) + t * ppc1->t[1];

        patch_resolve_color_inline(ppcr, pfs);

    } else {

        int ci;

        for (ci = pfs->num_components - 1; ci >= 0; --ci)

            ppcr->cc.paint.values[ci] =

                ppc0->cc.paint.values[ci] * (1 - t) + t * ppc1->cc.paint.values[ci];

    }

}

/* ================ Specific shadings ================ */

/*

 * The curves are stored in a clockwise or counter-clockwise order that maps

 * to the patch definition in a non-intuitive way.  The documentation

 * (pp. 277-281 of the PostScript Language Reference Manual, Third Edition)

 * says that the curves are in the order D1, C2, D2, C1.

 */

/* The starting points of the curves: */

#define D1START 0

#define C2START 1

#define D2START 3

#define C1START 0

/* The control points of the curves (x means reversed order): */

#define D1CTRL 0

#define C2CTRL 1

#define D2XCTRL 2

#define C1XCTRL 3

/* The end points of the curves: */

#define D1END 1

#define C2END 2

#define D2END 2

#define C1END 3

/* ---------------- Common code ---------------- */

/* Evaluate a curve at a given point. */

static void

curve_eval(gs_fixed_point * pt, const gs_fixed_point * p0,

           const gs_fixed_point * p1, const gs_fixed_point * p2,

           const gs_fixed_point * p3, double t)

{

    fixed a, b, c, d;

    fixed t01, t12;

    d = p0->x;

    curve_points_to_coefficients(d, p1->x, p2->x, p3->x,

                                 a, b, c, t01, t12);

    pt->x = (fixed) (((a * t + b) * t + c) * t + d);

    d = p0->y;

    curve_points_to_coefficients(d, p1->y, p2->y, p3->y,

                                 a, b, c, t01, t12);

    pt->y = (fixed) (((a * t + b) * t + c) * t + d);

    if_debug3('2', "[2]t=%g => (%g,%g)\n", t, fixed2float(pt->x),

              fixed2float(pt->y));

}

/* ---------------- Coons patch shading ---------------- */

/* Calculate the device-space coordinate corresponding to (u,v). */

static void

Cp_transform(gs_fixed_point * pt, const patch_curve_t curve[4],

             const gs_fixed_point ignore_interior[4], double u, double v)

{

    double co_u = 1.0 - u, co_v = 1.0 - v;

    gs_fixed_point c1u, d1v, c2u, d2v;

    curve_eval(&c1u, &curve[C1START].vertex.p,

               &curve[C1XCTRL].control[1], &curve[C1XCTRL].control[0],

               &curve[C1END].vertex.p, u);

    curve_eval(&d1v, &curve[D1START].vertex.p,

               &curve[D1CTRL].control[0], &curve[D1CTRL].control[1],

               &curve[D1END].vertex.p, v);

    curve_eval(&c2u, &curve[C2START].vertex.p,

               &curve[C2CTRL].control[0], &curve[C2CTRL].control[1],

               &curve[C2END].vertex.p, u);

    curve_eval(&d2v, &curve[D2START].vertex.p,

               &curve[D2XCTRL].control[1], &curve[D2XCTRL].control[0],

               &curve[D2END].vertex.p, v);

#define COMPUTE_COORD(xy)\

    pt->xy = (fixed)\

        ((co_v * c1u.xy + v * c2u.xy) + (co_u * d1v.xy + u * d2v.xy) -\

         (co_v * (co_u * curve[C1START].vertex.p.xy +\

                  u * curve[C1END].vertex.p.xy) +\

          v * (co_u * curve[C2START].vertex.p.xy +\

               u * curve[C2END].vertex.p.xy)))

    COMPUTE_COORD(x);

    COMPUTE_COORD(y);

#undef COMPUTE_COORD

    if_debug4('2', "[2](u=%g,v=%g) => (%g,%g)\n",

              u, v, fixed2float(pt->x), fixed2float(pt->y));

}

int

gs_shading_Cp_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_Cp_t * const psh = (const gs_shading_Cp_t *)psh0;

    patch_fill_state_t state;

    shade_coord_stream_t cs;

    patch_curve_t curve[4];

    int code;

    code = mesh_init_fill_state((mesh_fill_state_t *) &state,

                         (const gs_shading_mesh_t *)psh0, rect_clip, dev, pgs);

    if (code < 0) {

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

        return code;

    }

    state.Function = psh->params.Function;

    code = init_patch_fill_state(&state);

    if(code < 0) {

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

        return code;

    }

    curve[0].straight = curve[1].straight = curve[2].straight = curve[3].straight = false;

    shade_next_init(&cs, (const gs_shading_mesh_params_t *)&psh->params, pgs);

    while ((code = shade_next_patch(&cs, psh->params.BitsPerFlag,

                                    curve, NULL)) == 0 &&

           (code = patch_fill(&state, curve, NULL, Cp_transform)) >= 0

        ) {

        DO_NOTHING;

    }

    if (term_patch_fill_state(&state))

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

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

    return min(code, 0);

}

/* ---------------- Tensor product patch shading ---------------- */

/* Calculate the device-space coordinate corresponding to (u,v). */

static void

Tpp_transform(gs_fixed_point * pt, const patch_curve_t curve[4],

              const gs_fixed_point interior[4], double u, double v)

{

    double Bu[4], Bv[4];

    gs_fixed_point pts[4][4];

    int i, j;

    double x = 0, y = 0;

    /* Compute the Bernstein polynomials of u and v. */

    {

        double u2 = u * u, co_u = 1.0 - u, co_u2 = co_u * co_u;

        double v2 = v * v, co_v = 1.0 - v, co_v2 = co_v * co_v;

        Bu[0] = co_u * co_u2, Bu[1] = 3 * u * co_u2,

            Bu[2] = 3 * u2 * co_u, Bu[3] = u * u2;

        Bv[0] = co_v * co_v2, Bv[1] = 3 * v * co_v2,

            Bv[2] = 3 * v2 * co_v, Bv[3] = v * v2;

    }

    /* Arrange the points into an indexable order. */

    pts[0][0] = curve[0].vertex.p;

    pts[0][1] = curve[0].control[0];

    pts[0][2] = curve[0].control[1];

    pts[0][3] = curve[1].vertex.p;

    pts[1][3] = curve[1].control[0];

    pts[2][3] = curve[1].control[1];

    pts[3][3] = curve[2].vertex.p;

    pts[3][2] = curve[2].control[0];

    pts[3][1] = curve[2].control[1];

    pts[3][0] = curve[3].vertex.p;

    pts[2][0] = curve[3].control[0];

    pts[1][0] = curve[3].control[1];

    pts[1][1] = interior[0];

    pts[2][1] = interior[1];

    pts[2][2] = interior[2];

    pts[1][2] = interior[3];

    /* Now compute the actual point. */

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

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

            double coeff = Bu[i] * Bv[j];

            x += pts[i][j].x * coeff, y += pts[i][j].y * coeff;

        }

    pt->x = (fixed)x, pt->y = (fixed)y;

}

int

gs_shading_Tpp_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_Tpp_t * const psh = (const gs_shading_Tpp_t *)psh0;

    patch_fill_state_t state;

    shade_coord_stream_t cs;

    patch_curve_t curve[4];

    gs_fixed_point interior[4];

    int code;

    code = mesh_init_fill_state((mesh_fill_state_t *) & state,

                         (const gs_shading_mesh_t *)psh0, rect_clip, dev, pgs);

    if (code < 0) {

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

        return code;

    }

    state.Function = psh->params.Function;

    code = init_patch_fill_state(&state);

    if(code < 0)

        return code;

    curve[0].straight = curve[1].straight = curve[2].straight = curve[3].straight = false;

    shade_next_init(&cs, (const gs_shading_mesh_params_t *)&psh->params, pgs);

    while ((code = shade_next_patch(&cs, psh->params.BitsPerFlag,

                                    curve, interior)) == 0) {

        /*

         * The order of points appears to be consistent with that for Coons

         * patches, which is different from that documented in Red Book 3.

         */

        gs_fixed_point swapped_interior[4];

        swapped_interior[0] = interior[0];

        swapped_interior[1] = interior[3];

        swapped_interior[2] = interior[2];

        swapped_interior[3] = interior[1];

        code = patch_fill(&state, curve, swapped_interior, Tpp_transform);

        if (code < 0)

            break;

    }

    if (term_patch_fill_state(&state))

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

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

    return min(code, 0);

}

/*

    This algorithm performs a decomposition of the shading area

    into a set of constant color trapezoids, some of which

    may use the transpozed coordinate system.

    The target device assumes semi-open intrvals by X to be painted

    (See PLRM3, 7.5. Scan conversion details), i.e.

    it doesn't paint pixels which falls exactly to the right side.

    Note that with raster devices the algorithm doesn't paint pixels,

    whigh are partially covered by the shading area,

    but which's centers are outside the area.

    Pixels inside a monotonic part of the shading area are painted

    at once, but some exceptions may happen :

        - While flattening boundaries of a subpatch,

        to keep the plane coverage contiguity we insert wedges

        between neighbor subpatches, which use a different

        flattening factor. With non-monotonic curves

        those wedges may overlap or be self-overlapping, and a pixel

        is painted so many times as many wedges cover it. Fortunately

        the area of most wedges is zero or extremily small.

        - Since quazi-horizontal wedges may have a non-constant color,

        they can't decompose into constant color trapezoids with

        keeping the coverage contiguity. To represent them we

        apply the XY plane transposition. But with the transposition

        a semiopen interval can met a non-transposed one,

        so that some lines are not covered. Therefore we emulate

        closed intervals with expanding the transposed trapesoids in

        fixed_epsilon, and pixels at that boundary may be painted twice.

        - A boundary of a monotonic area can't compute in XY

        preciselly due to high order polynomial equations.

        Therefore the subdivision near the monotonity boundary

        may paint some pixels twice within same monotonic part.

    Non-monotonic areas slow down due to a tinny subdivision required.

    The target device may be either raster or vector.

    Vector devices should preciselly pass trapezoids to the output.

    Note that ends of sides of a trapesoid are not necessary

    the trapezoid's vertices. Converting this thing into

    an exact quadrangle may cause an arithmetic error,

    and the rounding must be done so that the coverage

    contiguity is not lost.

    When a device passes a trapezoid to it's output,

    a regular rounding would keep the coverage contiguity,

    except for the transposed trapesoids.

    If a transposed trapezoid is being transposed back,

    it doesn't become a canonic trapezoid, and a further

    decomposition is neccessary. But rounding errors here

    would break the coverage contiguity at boundaries

    of the tansposed part of the area.

    Devices, which have no transposed trapezoids and represent

    trapezoids only with 8 coordinates of vertices of the quadrangle

    (pclwrite is an example) may apply the backward transposition,

    and a clipping instead the further decomposition.

    Note that many clip regions may appear for all wedges.

    Note that in some cases the adjustment of the right side to be

    withdrown before the backward transposition.

 */

 /* We believe that a multiplication of 32-bit integers with a

    64-bit result is performed by modern platforms performs

    in hardware level. Therefore we widely use it here,

    but we minimize the usage of a multiplication of longer integers.

    Unfortunately we do need a multiplication of long integers

    in intersection_of_small_bars, because solving the linear system

    requires tripple multiples of 'fixed'. Therefore we retain

    of it's usage in the algorithm of the main branch.

    Configuration macro QUADRANGLES prevents it.

  */

typedef struct {

    gs_fixed_point pole[4][4]; /* [v][u] */

    patch_color_t *c[2][2];     /* [v][u] */

} tensor_patch;

typedef enum {

    interpatch_padding = 1, /* A Padding between patches for poorly designed documents. */

    inpatch_wedge = 2  /* Wedges while a patch decomposition. */

} wedge_type_t;

int

wedge_vertex_list_elem_buffer_alloc(patch_fill_state_t *pfs)

{

    const int max_level = LAZY_WEDGES_MAX_LEVEL;

    gs_memory_t *memory = pfs->memory;

    /* We have 'max_level' levels, each of which divides 1 or 3 sides.

       LAZY_WEDGES stores all 2^level divisions until

       the other area of same bnoundary is processed.

       Thus the upper estimation of the buffer size is :

       max_level * (1 << max_level) * 3.

       Likely this estimation to be decreased to

       max_level * (1 << max_level) * 2.

       because 1 side of a triangle is always outside the division path.

       For now we set the smaller estimation for obtaining an experimental data

       from the wild. */

    pfs->wedge_vertex_list_elem_count_max = max_level * (1 << max_level) * 2;

    pfs->wedge_vertex_list_elem_buffer = (wedge_vertex_list_elem_t *)gs_alloc_bytes(memory,

            sizeof(wedge_vertex_list_elem_t) * pfs->wedge_vertex_list_elem_count_max,

            "alloc_wedge_vertex_list_elem_buffer");

    if (pfs->wedge_vertex_list_elem_buffer == NULL)

        return_error(gs_error_VMerror);

    pfs->free_wedge_vertex = NULL;

    pfs->wedge_vertex_list_elem_count = 0;

    return 0;

}

void

wedge_vertex_list_elem_buffer_free(patch_fill_state_t *pfs)

{

    gs_memory_t *memory = pfs->memory;

    gs_free_object(memory, pfs->wedge_vertex_list_elem_buffer,

                "wedge_vertex_list_elem_buffer_free");

    pfs->wedge_vertex_list_elem_buffer = NULL;

    pfs->free_wedge_vertex = NULL;

}

static inline wedge_vertex_list_elem_t *

wedge_vertex_list_elem_reserve(patch_fill_state_t *pfs)

{

    wedge_vertex_list_elem_t *e = pfs->free_wedge_vertex;

    if (e != NULL) {

        pfs->free_wedge_vertex = e->next;

        return e;

    }

    if (pfs->wedge_vertex_list_elem_count < pfs->wedge_vertex_list_elem_count_max)

        return pfs->wedge_vertex_list_elem_buffer + pfs->wedge_vertex_list_elem_count++;

    return NULL;

}

static inline void

wedge_vertex_list_elem_release(patch_fill_state_t *pfs, wedge_vertex_list_elem_t *e)

{

    e->next = pfs->free_wedge_vertex;

    pfs->free_wedge_vertex = e;

}

static inline void

release_wedge_vertex_list_interval(patch_fill_state_t *pfs,

    wedge_vertex_list_elem_t *beg, wedge_vertex_list_elem_t *end)

{

    wedge_vertex_list_elem_t *e = beg->next, *ee;

    beg->next = end;

    end->prev = beg;

    for (; e != end; e = ee) {

        ee = e->next;

        wedge_vertex_list_elem_release(pfs, e);

    }

}

static inline int

release_wedge_vertex_list(patch_fill_state_t *pfs, wedge_vertex_list_t *ll, int n)

{

    int i;

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

        wedge_vertex_list_t *l = ll + i;

        if (l->beg != NULL) {

            if (l->end == NULL)

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

            release_wedge_vertex_list_interval(pfs, l->beg, l->end);

            wedge_vertex_list_elem_release(pfs, l->beg);

            wedge_vertex_list_elem_release(pfs, l->end);

            l->beg = l->end = NULL;

        } else if (l->end != NULL)

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

    }

    return 0;

}

static inline wedge_vertex_list_elem_t *

wedge_vertex_list_find(wedge_vertex_list_elem_t *beg, const wedge_vertex_list_elem_t *end,

            int level)

{

    wedge_vertex_list_elem_t *e = beg;

    if (beg == NULL || end == NULL)

        return NULL; /* Must not happen. */

    for (; e != end; e = e->next)

        if (e->level == level)

            return e;

    return NULL;

}

static inline void

init_wedge_vertex_list(wedge_vertex_list_t *l, int n)

{

    memset(l, 0, sizeof(*l) * n);

}

/* For a given set of poles in the tensor patch (for instance

 * [0][0], [0][1], [0][2], [0][3] or [0][2], [1][2], [2][2], [3][2])

 * return the number of subdivisions required to flatten the bezier

 * given by those poles to the required flatness.

 */

static inline int

curve_samples(patch_fill_state_t *pfs,

                const gs_fixed_point *pole, int pole_step, fixed fixed_flat)

{

    curve_segment s;

    int k;

    s.p1.x = pole[pole_step].x;

    s.p1.y = pole[pole_step].y;

    s.p2.x = pole[pole_step * 2].x;

    s.p2.y = pole[pole_step * 2].y;

    s.pt.x = pole[pole_step * 3].x;

    s.pt.y = pole[pole_step * 3].y;

    k = gx_curve_log2_samples(pole[0].x, pole[0].y, &s, fixed_flat);

    {

#       if LAZY_WEDGES || QUADRANGLES

            int k1;

            fixed L = any_abs(pole[pole_step * 1].x - pole[pole_step * 0].x) + any_abs(pole[pole_step * 1].y - pole[pole_step * 0].y) +

                      any_abs(pole[pole_step * 2].x - pole[pole_step * 1].x) + any_abs(pole[pole_step * 2].y - pole[pole_step * 1].y) +

                      any_abs(pole[pole_step * 3].x - pole[pole_step * 2].x) + any_abs(pole[pole_step * 3].y - pole[pole_step * 2].y);

#       endif

#       if LAZY_WEDGES

            /* Restrict lengths for a reasonable memory consumption : */

            k1 = ilog2(L / fixed_1 / (1 << (LAZY_WEDGES_MAX_LEVEL - 1)));

            k = max(k, k1);

#       endif

#       if QUADRANGLES

            /* Restrict lengths for intersection_of_small_bars : */

            k = max(k, ilog2(L) - ilog2(pfs->max_small_coord));

#       endif

    }

    return 1 << k;

}

static inline bool

intersection_of_small_bars(const gs_fixed_point q[4], int i0, int i1, int i2, int i3, fixed *ry, fixed *ey)

{

    /* This function is only used with QUADRANGLES. */

    return gx_intersect_small_bars(q[i0].x, q[i0].y, q[i1].x, q[i1].y, q[i2].x, q[i2].y, q[i3].x, q[i3].y, ry, ey);

}

static inline void

adjust_swapped_boundary(fixed *b, bool swap_axes)

{

    if (swap_axes) {

        /*  Sinse the rasterizer algorithm assumes semi-open interval

            when computing pixel coverage, we should expand

            the right side of the area. Otherwise a dropout can happen :

            if the left neighbour is painted with !swap_axes,

            the left side of this area appears to be the left side

            of the neighbour area, and the side is not included

            into both areas.

         */

        *b += fixed_epsilon;

    }

}

static inline void

make_trapezoid(const gs_fixed_point q[4],

        int vi0, int vi1, int vi2, int vi3, fixed ybot, fixed ytop,

        bool swap_axes, bool orient, gs_fixed_edge *le, gs_fixed_edge *re)

{

    if (!orient) {

        le->start = q[vi0];

        le->end = q[vi1];

        re->start = q[vi2];

        re->end = q[vi3];

    } else {

        le->start = q[vi2];

        le->end = q[vi3];

        re->start = q[vi0];

        re->end = q[vi1];

    }

    adjust_swapped_boundary(&re->start.x, swap_axes);

    adjust_swapped_boundary(&re->end.x, swap_axes);

}

static inline int

gx_shade_trapezoid(patch_fill_state_t *pfs, const gs_fixed_point q[4],

        int vi0, int vi1, int vi2, int vi3, fixed ybot0, fixed ytop0,

        bool swap_axes, const gx_device_color *pdevc, bool orient)

{

    gs_fixed_edge le, re;

    int code;

    fixed ybot = max(ybot0, swap_axes ? pfs->rect.p.x : pfs->rect.p.y);

    fixed ytop = min(ytop0, swap_axes ? pfs->rect.q.x : pfs->rect.q.y);

    fixed xleft  = (swap_axes ? pfs->rect.p.y : pfs->rect.p.x);

    fixed xright = (swap_axes ? pfs->rect.q.y : pfs->rect.q.x);

    if (ybot >= ytop)

        return 0;

#   if NOFILL_TEST

        if (dbg_nofill)

            return 0;

#   endif

    make_trapezoid(q, vi0, vi1, vi2, vi3, ybot, ytop, swap_axes, orient, &le, &re);

    if (!pfs->inside) {

        bool clip = false;

        /* We are asked to clip a trapezoid to a rectangle. If the rectangle

         * is entirely contained within the rectangle, then no clipping is

         * actually required. If the left edge is entirely to the right of

         * the rectangle, or the right edge is entirely to the left, we

         * clip away to nothing. If the left edge is entirely to the left of

         * the rectangle, then we can simplify it to a vertical edge along

         * the edge of the rectangle. Likewise with the right edge if it's

         * entirely to the right of the rectangle.*/

        if (le.start.x > xright) {

            if (le.end.x > xright) {

                return 0;

            }

            clip = true;

        } else if (le.end.x > xright) {

            clip = true;

        }

        if (le.start.x < xleft) {

            if (le.end.x < xleft) {

                le.start.x = xleft;

                le.end.x   = xleft;

                le.start.y = ybot;

                le.end.y   = ytop;

            } else {

                clip = true;

            }

        } else if (le.end.x < xleft) {

            clip = true;

        }

        if (re.start.x < xleft) {

            if (re.end.x < xleft) {

                return 0;

            }

            clip = true;

        } else if (re.end.x < xleft) {

            clip = true;

        }

        if (re.start.x > xright) {

            if (re.end.x > xright) {

                re.start.x = xright;

                re.end.x   = xright;

                re.start.y = ybot;

                re.end.y   = ytop;

            } else {

                clip = true;

            }

        } else if (re.end.x > xright) {

            clip = true;

        }

        if (clip)

        {

            /* Some form of clipping seems to be required. A certain amount

             * of horridness goes on here to ensure that we round 'outwards'

             * in all cases. */

            gs_fixed_edge lenew, renew;

            fixed ybl, ybr, ytl, ytr, ymid;

            /* Reduce the clipping region horizontally if possible. */

            if (re.start.x > re.end.x) {

                if (re.start.x < xright)

                    xright = re.start.x;

            } else if (re.end.x < xright)

                xright = re.end.x;

            if (le.start.x > le.end.x) {

                if (le.end.x > xleft)

                    xleft = le.end.x;

            } else if (le.start.x > xleft)

                xleft = le.start.x;

            ybot = max(ybot, min(le.start.y, re.start.y));

            ytop = min(ytop, max(le.end.y, re.end.y));

#if 0

            /* RJW: I've disabled this code a) because it doesn't make any

             * difference in the cluster tests, and b) because I think it's wrong.

             * Taking the first case as an example; just because the le.start.x

             * is > xright, does not mean that we can simply truncate the edge at

             * xright, as this may throw away part of the trap between ybot and

             * the new le.start.y. */

            /* Reduce the edges to the left/right of the clipping region. */

            /* Only in the 4 cases which can bring ytop/ybot closer */

            if (le.start.x > xright) {

                le.start.y += (fixed)((int64_t)(le.end.y-le.start.y)*

                                      (int64_t)(le.start.x-xright)/

                                      (int64_t)(le.start.x-le.end.x));

                le.start.x = xright;

            }

            if (re.start.x < xleft) {

                re.start.y += (fixed)((int64_t)(re.end.y-re.start.y)*

                                      (int64_t)(xleft-re.start.x)/

                                      (int64_t)(re.end.x-re.start.x));

                re.start.x = xleft;

            }

            if (le.end.x > xright) {

                le.end.y -= (fixed)((int64_t)(le.end.y-le.start.y)*

                                    (int64_t)(le.end.x-xright)/

                                    (int64_t)(le.end.x-le.start.x));

                le.end.x = xright;

            }

            if (re.end.x < xleft) {

                re.end.y -= (fixed)((int64_t)(re.end.y-re.start.y)*

                                    (int64_t)(xleft-re.end.x)/

                                    (int64_t)(re.start.x-re.end.x));

                re.end.x = xleft;

            }

#endif

            if (ybot >= ytop)

                return 0;

            /* Follow the edges in, so that le.start.y == ybot etc. */

            if (le.start.y < ybot) {

                int round = ((le.end.x < le.start.x) ?

                             (le.end.y-le.start.y-1) : 0);

                le.start.x += (fixed)(((int64_t)(ybot-le.start.y)*

                                       (int64_t)(le.end.x-le.start.x)-round)/

                                      (int64_t)(le.end.y-le.start.y));

                le.start.y = ybot;

            }

            if (le.end.y > ytop) {

                int round = ((le.end.x > le.start.x) ?

                             (le.end.y-le.start.y-1) : 0);

                le.end.x += (fixed)(((int64_t)(le.end.y-ytop)*

                                     (int64_t)(le.start.x-le.end.x)-round)/

                                    (int64_t)(le.end.y-le.start.y));

                le.end.y = ytop;

            }

            if ((le.start.x < xleft) && (le.end.x < xleft)) {

                le.start.x = xleft;

                le.end.x   = xleft;

                le.start.y = ybot;

                le.end.y   = ytop;

            }

            if (re.start.y < ybot) {

                int round = ((re.end.x > re.start.x) ?

                             (re.end.y-re.start.y-1) : 0);

                re.start.x += (fixed)(((int64_t)(ybot-re.start.y)*

                                       (int64_t)(re.end.x-re.start.x)+round)/

                                      (int64_t)(re.end.y-re.start.y));

                re.start.y = ybot;

            }

            if (re.end.y > ytop) {

                int round = ((re.end.x < re.start.x) ?

                             (re.end.y-re.start.y-1) : 0);

                re.end.x += (fixed)(((int64_t)(re.end.y-ytop)*

                                     (int64_t)(re.start.x-re.end.x)+round)/

                                    (int64_t)(re.end.y-re.start.y));

                re.end.y = ytop;

            }