Subversion Repositories planix.SVN

/* Copyright (C) 1992, 1995, 1996, 1997, 1998, 1999 Aladdin Enterprises.  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.

  For more information about licensing, please refer to

  http://www.ghostscript.com/licensing/. For information on

  commercial licensing, go to http://www.artifex.com/licensing/ or

  contact Artifex Software, Inc., 101 Lucas Valley Road #110,

  San Rafael, CA  94903, U.S.A., +1(415)492-9861.

*/

/* $Id: gscie.c,v 1.16 2004/03/16 02:16:20 dan Exp $ */

/* CIE color rendering cache management */

#include "math_.h"

#include "memory_.h"

#include "gx.h"

#include "gserrors.h"

#include "gsstruct.h"

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

#include "gxcspace.h"		/* for gxcie.c */

#include "gscolor2.h"		/* for gs_set/currentcolorrendering */

#include "gxarith.h"

#include "gxcie.h"

#include "gxdevice.h"		/* for gxcmap.h */

#include "gxcmap.h"

#include "gzstate.h"

#include "gsicc.h"

/*

 * Define whether to optimize the CIE mapping process by combining steps.

 * This should only be disabled (commented out) for debugging.

 */

#define OPTIMIZE_CIE_MAPPING

/* Forward references */

private int cie_joint_caches_init(gx_cie_joint_caches *,

				  const gs_cie_common *,

				  gs_cie_render *);

private void cie_joint_caches_complete(gx_cie_joint_caches *,

				       const gs_cie_common *,

				       const gs_cie_abc *,

				       const gs_cie_render *);

private void cie_cache_restrict(cie_cache_floats *, const gs_range *);

private void cie_mult3(const gs_vector3 *, const gs_matrix3 *,

		       gs_vector3 *);

private void cie_matrix_mult3(const gs_matrix3 *, const gs_matrix3 *,

			      gs_matrix3 *);

private void cie_invert3(const gs_matrix3 *, gs_matrix3 *);

private void cie_matrix_init(gs_matrix3 *);

/* Allocator structure types */

private_st_joint_caches();

extern_st(st_imager_state);

#define RESTRICTED_INDEX(v, n, itemp)\

  ((uint)(itemp = (int)(v)) >= (n) ?\

   (itemp < 0 ? 0 : (n) - 1) : itemp)

/* Define the template for loading a cache. */

/* If we had parameterized types, or a more flexible type system, */

/* this could be done with a single procedure. */

#define CIE_LOAD_CACHE_BODY(pcache, domains, rprocs, dprocs, pcie, cname)\

  BEGIN\

	int j;\

\

	for (j = 0; j < countof(pcache); j++) {\

	    cie_cache_floats *pcf = &(pcache)[j].floats;\

	    int i;\

	    gs_sample_loop_params_t lp;\

\

	    gs_cie_cache_init(&pcf->params, &lp, &(domains)[j], cname);\

	    for (i = 0; i <= lp.N; ++i) {\

		float v = SAMPLE_LOOP_VALUE(i, lp);\

		pcf->values[i] = (*(rprocs)->procs[j])(v, pcie);\

		if_debug5('C', "[C]%s[%d,%d] = %g => %g\n",\

			  cname, j, i, v, pcf->values[i]);\

	    }\

	    pcf->params.is_identity =\

		(rprocs)->procs[j] == (dprocs).procs[j];\

	}\

  END

/* Define cache interpolation threshold values. */

#ifdef CIE_CACHE_INTERPOLATE

#  ifdef CIE_INTERPOLATE_THRESHOLD

#    define CACHE_THRESHOLD CIE_INTERPOLATE_THRESHOLD

#  else

#    define CACHE_THRESHOLD 0	/* always interpolate */

#  endif

#else

#  define CACHE_THRESHOLD 1.0e6	/* never interpolate */

#endif

#ifdef CIE_RENDER_TABLE_INTERPOLATE

#  define RENDER_TABLE_THRESHOLD 0

#else

#  define RENDER_TABLE_THRESHOLD 1.0e6

#endif

/*

 * Determine whether a function is a linear transformation of the form

 * f(x) = scale * x + origin.

 */

private bool

cache_is_linear(cie_linear_params_t *params, const cie_cache_floats *pcf)

{

    double origin = pcf->values[0];

    double diff = pcf->values[countof(pcf->values) - 1] - origin;

    double scale = diff / (countof(pcf->values) - 1);

    int i;

    double test = origin + scale;

    for (i = 1; i < countof(pcf->values) - 1; ++i, test += scale)

	if (fabs(pcf->values[i] - test) >= 0.5 / countof(pcf->values))

	    return (params->is_linear = false);

    params->origin = origin - pcf->params.base;

    params->scale =

	diff * pcf->params.factor / (countof(pcf->values) - 1);

    return (params->is_linear = true);

}

private void

cache_set_linear(cie_cache_floats *pcf)

{

	if (pcf->params.is_identity) {

	    if_debug1('c', "[c]is_linear(0x%lx) = true (is_identity)\n",

		      (ulong)pcf);

	    pcf->params.linear.is_linear = true;

	    pcf->params.linear.origin = 0;

	    pcf->params.linear.scale = 1;

	} else if (cache_is_linear(&pcf->params.linear, pcf)) {

	    if (pcf->params.linear.origin == 0 &&

		fabs(pcf->params.linear.scale - 1) < 0.00001)

		pcf->params.is_identity = true;

	    if_debug4('c',

		      "[c]is_linear(0x%lx) = true, origin = %g, scale = %g%s\n",

		      (ulong)pcf, pcf->params.linear.origin,

		      pcf->params.linear.scale,

		      (pcf->params.is_identity ? " (=> is_identity)" : ""));

	}

#ifdef DEBUG

	else

	    if_debug1('c', "[c]linear(0x%lx) = false\n", (ulong)pcf);

#endif

}

private void

cache3_set_linear(gx_cie_vector_cache3_t *pvc)

{

    cache_set_linear(&pvc->caches[0].floats);

    cache_set_linear(&pvc->caches[1].floats);

    cache_set_linear(&pvc->caches[2].floats);

}

#ifdef DEBUG

private void

if_debug_vector3(const char *str, const gs_vector3 *vec)

{

    if_debug4('c', "%s[%g %g %g]\n", str, vec->u, vec->v, vec->w);

}

private void

if_debug_matrix3(const char *str, const gs_matrix3 *mat)

{

    if_debug10('c', "%s [%g %g %g] [%g %g %g] [%g %g %g]\n", str,

	       mat->cu.u, mat->cu.v, mat->cu.w,

	       mat->cv.u, mat->cv.v, mat->cv.w,

	       mat->cw.u, mat->cw.v, mat->cw.w);

}

#else

#  define if_debug_vector3(str, vec) DO_NOTHING

#  define if_debug_matrix3(str, mat) DO_NOTHING

#endif

/* ------ Default values for CIE dictionary elements ------ */

/* Default transformation procedures. */

private float

a_identity(floatp in, const gs_cie_a * pcie)

{

    return in;

}

private float

a_from_cache(floatp in, const gs_cie_a * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches.DecodeA.floats);

}

private float

abc_identity(floatp in, const gs_cie_abc * pcie)

{

    return in;

}

private float

abc_from_cache_0(floatp in, const gs_cie_abc * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches.DecodeABC.caches[0].floats);

}

private float

abc_from_cache_1(floatp in, const gs_cie_abc * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches.DecodeABC.caches[1].floats);

}

private float

abc_from_cache_2(floatp in, const gs_cie_abc * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches.DecodeABC.caches[2].floats);

}

private float

def_identity(floatp in, const gs_cie_def * pcie)

{

    return in;

}

private float

def_from_cache_0(floatp in, const gs_cie_def * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[0].floats);

}

private float

def_from_cache_1(floatp in, const gs_cie_def * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[1].floats);

}

private float

def_from_cache_2(floatp in, const gs_cie_def * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[2].floats);

}

private float

defg_identity(floatp in, const gs_cie_defg * pcie)

{

    return in;

}

private float

defg_from_cache_0(floatp in, const gs_cie_defg * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[0].floats);

}

private float

defg_from_cache_1(floatp in, const gs_cie_defg * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[1].floats);

}

private float

defg_from_cache_2(floatp in, const gs_cie_defg * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[2].floats);

}

private float

defg_from_cache_3(floatp in, const gs_cie_defg * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[3].floats);

}

private float

common_identity(floatp in, const gs_cie_common * pcie)

{

    return in;

}

private float

lmn_from_cache_0(floatp in, const gs_cie_common * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[0].floats);

}

private float

lmn_from_cache_1(floatp in, const gs_cie_common * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[1].floats);

}

private float

lmn_from_cache_2(floatp in, const gs_cie_common * pcie)

{

    return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[2].floats);

}

/* Transformation procedures for accessing an already-loaded cache. */

float

gs_cie_cached_value(floatp in, const cie_cache_floats *pcache)

{

    /*

     * We need to get the same results when we sample an already-loaded

     * cache, so we need to round the index just a tiny bit.

     */

    int index =

	(int)((in - pcache->params.base) * pcache->params.factor + 0.0001);

    CIE_CLAMP_INDEX(index);

    return pcache->values[index];

}

/* Default vectors and matrices. */

const gs_range3 Range3_default = {

    { {0, 1}, {0, 1}, {0, 1} }

};

const gs_range4 Range4_default = {

    { {0, 1}, {0, 1}, {0, 1}, {0, 1} }

};

const gs_cie_defg_proc4 DecodeDEFG_default = {

    {defg_identity, defg_identity, defg_identity, defg_identity}

};

const gs_cie_defg_proc4 DecodeDEFG_from_cache = {

    {defg_from_cache_0, defg_from_cache_1, defg_from_cache_2, defg_from_cache_3}

};

const gs_cie_def_proc3 DecodeDEF_default = {

    {def_identity, def_identity, def_identity}

};

const gs_cie_def_proc3 DecodeDEF_from_cache = {

    {def_from_cache_0, def_from_cache_1, def_from_cache_2}

};

const gs_cie_abc_proc3 DecodeABC_default = {

    {abc_identity, abc_identity, abc_identity}

};

const gs_cie_abc_proc3 DecodeABC_from_cache = {

    {abc_from_cache_0, abc_from_cache_1, abc_from_cache_2}

};

const gs_cie_common_proc3 DecodeLMN_default = {

    {common_identity, common_identity, common_identity}

};

const gs_cie_common_proc3 DecodeLMN_from_cache = {

    {lmn_from_cache_0, lmn_from_cache_1, lmn_from_cache_2}

};

const gs_matrix3 Matrix3_default = {

    {1, 0, 0},

    {0, 1, 0},

    {0, 0, 1},

    1 /*true */

};

const gs_range RangeA_default = {0, 1};

const gs_cie_a_proc DecodeA_default = a_identity;

const gs_cie_a_proc DecodeA_from_cache = a_from_cache;

const gs_vector3 MatrixA_default = {1, 1, 1};

const gs_vector3 BlackPoint_default = {0, 0, 0};

/* Initialize a CIE color. */

/* This only happens on setcolorspace. */

void

gx_init_CIE(gs_client_color * pcc, const gs_color_space * pcs)

{

    gx_init_paint_4(pcc, pcs);

    /* (0...) may not be within the range of allowable values. */

    (*pcs->type->restrict_color)(pcc, pcs);

}

/* Restrict CIE colors. */

inline private void

cie_restrict(float *pv, const gs_range *range)

{

    if (*pv <= range->rmin)

	*pv = range->rmin;

    else if (*pv >= range->rmax)

	*pv = range->rmax;

}

void

gx_restrict_CIEDEFG(gs_client_color * pcc, const gs_color_space * pcs)

{

    const gs_cie_defg *pcie = pcs->params.defg;

    cie_restrict(&pcc->paint.values[0], &pcie->RangeDEFG.ranges[0]);

    cie_restrict(&pcc->paint.values[1], &pcie->RangeDEFG.ranges[1]);

    cie_restrict(&pcc->paint.values[2], &pcie->RangeDEFG.ranges[2]);

    cie_restrict(&pcc->paint.values[3], &pcie->RangeDEFG.ranges[3]);

}

void

gx_restrict_CIEDEF(gs_client_color * pcc, const gs_color_space * pcs)

{

    const gs_cie_def *pcie = pcs->params.def;

    cie_restrict(&pcc->paint.values[0], &pcie->RangeDEF.ranges[0]);

    cie_restrict(&pcc->paint.values[1], &pcie->RangeDEF.ranges[1]);

    cie_restrict(&pcc->paint.values[2], &pcie->RangeDEF.ranges[2]);

}

void

gx_restrict_CIEABC(gs_client_color * pcc, const gs_color_space * pcs)

{

    const gs_cie_abc *pcie = pcs->params.abc;

    cie_restrict(&pcc->paint.values[0], &pcie->RangeABC.ranges[0]);

    cie_restrict(&pcc->paint.values[1], &pcie->RangeABC.ranges[1]);

    cie_restrict(&pcc->paint.values[2], &pcie->RangeABC.ranges[2]);

}

void

gx_restrict_CIEA(gs_client_color * pcc, const gs_color_space * pcs)

{

    const gs_cie_a *pcie = pcs->params.a;

    cie_restrict(&pcc->paint.values[0], &pcie->RangeA);

}

/* ================ Table setup ================ */

/* ------ Install a CIE color space ------ */

private void cie_cache_mult(gx_cie_vector_cache *, const gs_vector3 *,

			    const cie_cache_floats *, floatp);

private bool cie_cache_mult3(gx_cie_vector_cache3_t *,

			     const gs_matrix3 *, floatp);

private int

gx_install_cie_abc(gs_cie_abc *pcie, gs_state * pgs)

{

    if_debug_matrix3("[c]CIE MatrixABC =", &pcie->MatrixABC);

    cie_matrix_init(&pcie->MatrixABC);

    CIE_LOAD_CACHE_BODY(pcie->caches.DecodeABC.caches, pcie->RangeABC.ranges,

			&pcie->DecodeABC, DecodeABC_default, pcie,

			"DecodeABC");

    gx_cie_load_common_cache(&pcie->common, pgs);

    gs_cie_abc_complete(pcie);

    return gs_cie_cs_complete(pgs, true);

}

int

gx_install_CIEDEFG(const gs_color_space * pcs, gs_state * pgs)

{

    gs_cie_defg *pcie = pcs->params.defg;

    CIE_LOAD_CACHE_BODY(pcie->caches_defg.DecodeDEFG, pcie->RangeDEFG.ranges,

			&pcie->DecodeDEFG, DecodeDEFG_default, pcie,

			"DecodeDEFG");

    return gx_install_cie_abc((gs_cie_abc *)pcie, pgs);

}

int

gx_install_CIEDEF(const gs_color_space * pcs, gs_state * pgs)

{

    gs_cie_def *pcie = pcs->params.def;

    CIE_LOAD_CACHE_BODY(pcie->caches_def.DecodeDEF, pcie->RangeDEF.ranges,

			&pcie->DecodeDEF, DecodeDEF_default, pcie,

			"DecodeDEF");

    return gx_install_cie_abc((gs_cie_abc *)pcie, pgs);

}

int

gx_install_CIEABC(const gs_color_space * pcs, gs_state * pgs)

{

    return gx_install_cie_abc(pcs->params.abc, pgs);

}

int

gx_install_CIEA(const gs_color_space * pcs, gs_state * pgs)

{

    gs_cie_a *pcie = pcs->params.a;

    gs_sample_loop_params_t lp;

    int i;

    gs_cie_cache_init(&pcie->caches.DecodeA.floats.params, &lp,

		      &pcie->RangeA, "DecodeA");

    for (i = 0; i <= lp.N; ++i) {

	float in = SAMPLE_LOOP_VALUE(i, lp);

	pcie->caches.DecodeA.floats.values[i] = (*pcie->DecodeA)(in, pcie);

	if_debug3('C', "[C]DecodeA[%d] = %g => %g\n",

		  i, in, pcie->caches.DecodeA.floats.values[i]);

    }

    gx_cie_load_common_cache(&pcie->common, pgs);

    gs_cie_a_complete(pcie);

    return gs_cie_cs_complete(pgs, true);

}

/* Load the common caches when installing the color space. */

/* This routine is exported for the benefit of gsicc.c */

void

gx_cie_load_common_cache(gs_cie_common * pcie, gs_state * pgs)

{

    if_debug_matrix3("[c]CIE MatrixLMN =", &pcie->MatrixLMN);

    cie_matrix_init(&pcie->MatrixLMN);

    CIE_LOAD_CACHE_BODY(pcie->caches.DecodeLMN, pcie->RangeLMN.ranges,

			&pcie->DecodeLMN, DecodeLMN_default, pcie,

			"DecodeLMN");

}

/* Complete loading the common caches. */

/* This routine is exported for the benefit of gsicc.c */

void

gx_cie_common_complete(gs_cie_common *pcie)

{

    int i;

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

	cache_set_linear(&pcie->caches.DecodeLMN[i].floats);

}

/*

 * Restrict the DecodeDEF[G] cache according to RangeHIJ[K], and scale to

 * the dimensions of Table.

 */

private void

gs_cie_defx_scale(float *values, const gs_range *range, int dim)

{

    double scale = (dim - 1.0) / (range->rmax - range->rmin);

    int i;

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

	float value = values[i];

	values[i] =

	    (value <= range->rmin ? 0 :

	     value >= range->rmax ? dim - 1 :

	     (value - range->rmin) * scale);

    }

}

/* Complete loading a CIEBasedDEFG color space. */

/* This routine is NOT idempotent. */

void

gs_cie_defg_complete(gs_cie_defg * pcie)

{

    int j;

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

	gs_cie_defx_scale(pcie->caches_defg.DecodeDEFG[j].floats.values,

			  &pcie->RangeHIJK.ranges[j], pcie->Table.dims[j]);

    gs_cie_abc_complete((gs_cie_abc *)pcie);

}

/* Complete loading a CIEBasedDEF color space. */

/* This routine is NOT idempotent. */

void

gs_cie_def_complete(gs_cie_def * pcie)

{

    int j;

    for (j = 0; j < 3; ++j)

	gs_cie_defx_scale(pcie->caches_def.DecodeDEF[j].floats.values,

			  &pcie->RangeHIJ.ranges[j], pcie->Table.dims[j]);

    gs_cie_abc_complete((gs_cie_abc *)pcie);

}

/* Complete loading a CIEBasedABC color space. */

/* This routine is idempotent. */

void

gs_cie_abc_complete(gs_cie_abc * pcie)

{

    cache3_set_linear(&pcie->caches.DecodeABC);

    pcie->caches.skipABC =

	cie_cache_mult3(&pcie->caches.DecodeABC, &pcie->MatrixABC,

			CACHE_THRESHOLD);

    gx_cie_common_complete((gs_cie_common *)pcie);

}

/* Complete loading a CIEBasedA color space. */

/* This routine is idempotent. */

void

gs_cie_a_complete(gs_cie_a * pcie)

{

    cie_cache_mult(&pcie->caches.DecodeA, &pcie->MatrixA,

		   &pcie->caches.DecodeA.floats,

		   CACHE_THRESHOLD);

    cache_set_linear(&pcie->caches.DecodeA.floats);

    gx_cie_common_complete((gs_cie_common *)pcie);

}

/*

 * Set the ranges where interpolation is required in a vector cache.

 * This procedure is idempotent.

 */

typedef struct cie_cache_range_temp_s {

    cie_cached_value prev;

    int imin, imax;

} cie_cache_range_temp_t;

private void

check_interpolation_required(cie_cache_range_temp_t *pccr,

			     cie_cached_value cur, int i, floatp threshold)

{

    cie_cached_value prev = pccr->prev;

    if (any_abs(cur - prev) > threshold * min(any_abs(prev), any_abs(cur))) {

	if (i - 1 < pccr->imin)

	    pccr->imin = i - 1;

	if (i > pccr->imax)

	    pccr->imax = i;

    }

    pccr->prev = cur;

}

private void

cie_cache_set_interpolation(gx_cie_vector_cache *pcache, floatp threshold)

{

    cie_cached_value base = pcache->vecs.params.base;

    cie_cached_value factor = pcache->vecs.params.factor;

    cie_cache_range_temp_t temp[3];

    int i, j;

    for (j = 0; j < 3; ++j)

	temp[j].imin = gx_cie_cache_size, temp[j].imax = -1;

    temp[0].prev = pcache->vecs.values[0].u;

    temp[1].prev = pcache->vecs.values[0].v;

    temp[2].prev = pcache->vecs.values[0].w;

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

	check_interpolation_required(&temp[0], pcache->vecs.values[i].u, i,

				     threshold);

	check_interpolation_required(&temp[1], pcache->vecs.values[i].v, i,

				     threshold);

	check_interpolation_required(&temp[2], pcache->vecs.values[i].w, i,

				     threshold);

    }

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

	pcache->vecs.params.interpolation_ranges[j].rmin =

	    base + (cie_cached_value)((double)temp[j].imin / factor);

	pcache->vecs.params.interpolation_ranges[j].rmax =

	    base + (cie_cached_value)((double)temp[j].imax / factor);

	if_debug3('c', "[c]interpolation_ranges[%d] = %g, %g\n", j,

		  cie_cached2float(pcache->vecs.params.interpolation_ranges[j].rmin),

		  cie_cached2float(pcache->vecs.params.interpolation_ranges[j].rmax));

    }

}

/*

 * Convert a scalar cache to a vector cache by multiplying the scalar

 * values by a vector.  Also set the range where interpolation is needed.

 * This procedure is idempotent.

 */

private void

cie_cache_mult(gx_cie_vector_cache * pcache, const gs_vector3 * pvec,

	       const cie_cache_floats * pcf, floatp threshold)

{

    float u = pvec->u, v = pvec->v, w = pvec->w;

    int i;

    pcache->vecs.params.base = float2cie_cached(pcf->params.base);

    pcache->vecs.params.factor = float2cie_cached(pcf->params.factor);

    pcache->vecs.params.limit =

	float2cie_cached((gx_cie_cache_size - 1) / pcf->params.factor +

			 pcf->params.base);

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

	float f = pcf->values[i];

	pcache->vecs.values[i].u = float2cie_cached(f * u);

	pcache->vecs.values[i].v = float2cie_cached(f * v);

	pcache->vecs.values[i].w = float2cie_cached(f * w);

    }

    cie_cache_set_interpolation(pcache, threshold);

}

/*

 * Set the interpolation ranges in a 3-vector cache, based on the ranges in

 * the individual vector caches.  This procedure is idempotent.

 */

private void

cie_cache3_set_interpolation(gx_cie_vector_cache3_t * pvc)

{

    int j, k;

    /* Iterate over output components. */

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

	/* Iterate over sub-caches. */

	cie_interpolation_range_t *p = 

		&pvc->caches[0].vecs.params.interpolation_ranges[j];

        cie_cached_value rmin = p->rmin, rmax = p->rmax;

	for (k = 1; k < 3; ++k) {

	    p = &pvc->caches[k].vecs.params.interpolation_ranges[j];

	    rmin = min(rmin, p->rmin), rmax = max(rmax, p->rmax);

	}

	pvc->interpolation_ranges[j].rmin = rmin;

	pvc->interpolation_ranges[j].rmax = rmax;

	if_debug3('c', "[c]Merged interpolation_ranges[%d] = %g, %g\n",

		  j, rmin, rmax);

    }

}

/*

 * Convert 3 scalar caches to vector caches by multiplying by a matrix.

 * Return true iff the resulting cache is an identity transformation.

 * This procedure is idempotent.

 */

private bool

cie_cache_mult3(gx_cie_vector_cache3_t * pvc, const gs_matrix3 * pmat,

		floatp threshold)

{

    cie_cache_mult(&pvc->caches[0], &pmat->cu, &pvc->caches[0].floats, threshold);

    cie_cache_mult(&pvc->caches[1], &pmat->cv, &pvc->caches[1].floats, threshold);

    cie_cache_mult(&pvc->caches[2], &pmat->cw, &pvc->caches[2].floats, threshold);

    cie_cache3_set_interpolation(pvc);

    return pmat->is_identity & pvc->caches[0].floats.params.is_identity &

	pvc->caches[1].floats.params.is_identity &

	pvc->caches[2].floats.params.is_identity;

}

/* ------ Install a rendering dictionary ------ */

/* setcolorrendering */

int

gs_setcolorrendering(gs_state * pgs, gs_cie_render * pcrd)

{

    int code = gs_cie_render_complete(pcrd);

    const gs_cie_render *pcrd_old = pgs->cie_render;

    bool joint_ok;

    if (code < 0)

	return code;

    if (pcrd_old != 0 && pcrd->id == pcrd_old->id)

	return 0;		/* detect needless reselecting */

    joint_ok =

	pcrd_old != 0 &&

#define CRD_SAME(elt) !memcmp(&pcrd->elt, &pcrd_old->elt, sizeof(pcrd->elt))

	CRD_SAME(points.WhitePoint) && CRD_SAME(points.BlackPoint) &&

	CRD_SAME(MatrixPQR) && CRD_SAME(RangePQR) &&

	CRD_SAME(TransformPQR);

#undef CRD_SAME

    rc_assign(pgs->cie_render, pcrd, "gs_setcolorrendering");

    /* Initialize the joint caches if needed. */

    if (!joint_ok)

	code = gs_cie_cs_complete(pgs, true);

    gx_unset_dev_color(pgs);

    return code;

}

/* currentcolorrendering */

const gs_cie_render *

gs_currentcolorrendering(const gs_state * pgs)

{

    return pgs->cie_render;

}

/* Unshare (allocating if necessary) the joint caches. */

gx_cie_joint_caches *

gx_currentciecaches(gs_state * pgs)

{

    gx_cie_joint_caches *pjc = pgs->cie_joint_caches;

    rc_unshare_struct(pgs->cie_joint_caches, gx_cie_joint_caches,

		      &st_joint_caches, pgs->memory,

		      return 0, "gx_currentciecaches");

    if (pgs->cie_joint_caches != pjc) {

	pjc = pgs->cie_joint_caches;

	pjc->cspace_id = pjc->render_id = gs_no_id;

	pjc->id_status = pjc->status = CIE_JC_STATUS_BUILT;

    }

    return pjc;

}

/* Compute the parameters for loading a cache, setting base and factor. */

/* This procedure is idempotent. */

void

gs_cie_cache_init(cie_cache_params * pcache, gs_sample_loop_params_t * pslp,

		  const gs_range * domain, client_name_t cname)

{

    /*

      We need to map the values in the range [domain->rmin..domain->rmax].

      However, if rmin < 0 < rmax and the function is non-linear, this can

      lead to anomalies at zero, which is the default value for CIE colors.

      The "correct" way to approach this is to run the mapping functions on

      demand, but we don't want to deal with the complexities of the

      callbacks this would involve (especially in the middle of rendering

      images); instead, we adjust the range so that zero maps precisely to a

      cache slot.  Define:

      A = domain->rmin;

      B = domain->rmax;

      N = gx_cie_cache_size - 1;

      R = B - A;

      h(v) = N * (v - A) / R;		// the index of v in the cache

      X = h(0).

      If X is not an integer, we can decrease A and/increase B to make it

      one.  Let A' and B' be the adjusted values of A and B respectively,

      and let K be the integer derived from X (either floor(X) or ceil(X)).

      Define

      f(K) = (K * B' + (N - K) * A') / N).

      We want f(K) = 0.  This occurs precisely when, for any real number

      C != 0,

      A' = -K * C;

      B' = (N - K) * C.

      In order to ensure A' <= A and B' >= B, we require

      C >= -A / K;

      C >= B / (N - K).

      Since A' and B' must be exactly representable as floats, we round C

      upward to ensure that it has no more than M mantissa bits, where

      M = ARCH_FLOAT_MANTISSA_BITS - ceil(log2(N)).

    */

    float A = domain->rmin, B = domain->rmax;

    double R = B - A, delta;

#define NN (gx_cie_cache_size - 1) /* 'N' is a member name, see end of proc */

#define N NN

#define CEIL_LOG2_N CIE_LOG2_CACHE_SIZE

    /* Adjust the range if necessary. */

    if (A < 0 && B >= 0) {

	const double X = -N * A / R;	/* know X > 0 */

	/* Choose K to minimize range expansion. */

	const int K = (int)(A + B < 0 ? floor(X) : ceil(X)); /* know 0 < K < N */

	const double Ca = -A / K, Cb = B / (N - K); /* know Ca, Cb > 0 */

	double C = max(Ca, Cb);	/* know C > 0 */

	const int M = ARCH_FLOAT_MANTISSA_BITS - CEIL_LOG2_N;

	int cexp;

	const double cfrac = frexp(C, &cexp);

	if_debug4('c', "[c]adjusting cache_init(%8g, %8g), X = %8g, K = %d:\n",

		  A, B, X, K);

	/* Round C to no more than M significant bits.  See above. */

	C = ldexp(ceil(ldexp(cfrac, M)), cexp - M);

	/* Finally, compute A' and B'. */

	A = -K * C;

	B = (N - K) * C;

	if_debug2('c', "[c]  => %8g, %8g\n", A, B);

	R = B - A;

    }

    delta = R / N;

#ifdef CIE_CACHE_INTERPOLATE

    pcache->base = A;		/* no rounding */

#else

    pcache->base = A - delta / 2;	/* so lookup will round */

#endif

    /*

     * If size of the domain is zero, then use 1.0 as the scaling

     * factor.  This prevents divide by zero errors in later calculations.

     * This should only occurs with zero matrices.  It does occur with

     * Genoa test file 050-01.ps.

     */

    pcache->factor = (any_abs(delta) < 1e-30 ? 1.0 : N / R);

    if_debug4('c', "[c]cache %s 0x%lx base=%g, factor=%g\n",

	      (const char *)cname, (ulong) pcache,

	      pcache->base, pcache->factor);

    pslp->A = A;

    pslp->B = B;

#undef N

    pslp->N = NN;

#undef NN

}

/* ------ Complete a rendering structure ------ */

/*

 * Compute the derived values in a CRD that don't involve the cached

 * procedure values.  This procedure is idempotent.

 */

private void cie_transform_range3(const gs_range3 *, const gs_matrix3 *,

				  gs_range3 *);

int

gs_cie_render_init(gs_cie_render * pcrd)

{

    gs_matrix3 PQR_inverse;

    if (pcrd->status >= CIE_RENDER_STATUS_INITED)

	return 0;		/* init already done */

    if_debug_matrix3("[c]CRD MatrixLMN =", &pcrd->MatrixLMN);