path: root/rspl/rev.h

#ifndef RSPL_REV_H
#define RSPL_REV_H

/* 
 * Argyll Color Correction System
 * Multi-dimensional regularized spline data structure
 *
 * Reverse interpolation support code.
 *
 * Author: Graeme W. Gill
 * Date:   30/1/00
 *
 * Copyright 1999 - 2008 Graeme W. Gill
 * All rights reserved.
 *
 * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
 * see the License.txt file for licencing details.
 *
 * Latest simplex/linear equation version.
 */

#undef STATS		/* Collect and print reverse cach stats */

/* Data structures used by reverse lookup code. */
/* Note that the reverse lookup code only supports a more limited */
/* dimension range than the general rspl code. */

/*
 * Note on simplex parameter space. 
 *
 * Simplex interpolation is normaly done in Baricentric space, where there
 * is one more baricentric coordinate than dimensions, and the sum of all
 * the baricentric coordinates must be 1. 
 *
 * To simplify things, we work in a "Simplex parameter" space, in which
 * there are only dimension parameters, and each directly corresponds
 * with a cartesian coordinate, but the parameter order corresponds with
 * the baricentric order.
 *
 * For example, given cartesian sub-coordinates D0, D1, D2
 * into a (3D) forward interpolation cube, these should be sorted
 * from smallest to largest, thereby choosing a particular
 * simplex within a cube, and allowing a correspondence to
 * the parameter coordinates, ie:
 *
 * D2 D0 D1		Smallest -> Largest cartesian sort
 * P0 P1 P2 	Corresponding Parameter coordinates
 * 
 * B0 = P0		Conversion to Baricentric weighting/coordinates
 * B1 = P1 - P0
 * B2 = P2 - P1
 * B3 = 1  - P2
 *
 * The 4 (tetrahedron) vertex values directly correspond to Baricentric
 * weighting/coordinates, giving the usual interpolation equation of:
 *
 *   VV0 * B0
 * + VV1 * B1
 * + VV2 * B2
 * + VV3 * B3
 *
 * Reversing the Parameter -> Baricentric equations gives the
 * following interpolation equation using Parameter coordinates:
 *
 *   VV0 - VV1 * P0
 * + VV1 - VV2 * P1
 * + VV2 - VV3 * P2
 * + VV3
 *
 * Note that withing the simplex, 0 <= P0 && P0 <= P1 && P1 <= P2 && P2 <= 1
 *
 * It is this which is used in rev.c for solving the reverse
 * interpolation problem.
 */

/* - - - - - - - - - - - - - - - - - - - - - */
/* Group size information for nn LCh weighted quick accept/reject testing */
typedef struct {
	double bcent[MXRO];	/* Group center location in output space */ 
	double brad;		/* Output value bounding shere radius */
	double bradsq;		/* Output value bounding shere radius squared */

	double maxDlc;		/* Maximum members weighted delta LC (+extra dims) squared */
	double maxDh;		/* Maximum members delta h squared */
	double maxDh_;		/* Maximum members delta h (not squared) */
	double sratio;		/* Minimum members C ratio to Gc squared */
	double bratio;		/* Maximum members C ratio to Gc squared */
	double Wsratio;		/* Minimum members C ratio to Gc squared - pre-weighted */
	double Wbratio;		/* Maximum members C ratio to Gc squared - pre-weighted */
	double Gc;			/* Group center Chrominance squared */
	double Gc_;			/* Group center Chrominance (not squared) */
} nn_grp;

/* - - - - - - - - - - - - - - - - - - - - - */
/* A structure to hold per simplex coordinate and vertex mapping for ssxinfo. */
/* This is relative to the construction cube. A face sub-simplex */
/* that is common between cubes, will have a different psxinfo */
/* depending on which cube created it. */
typedef struct {
	int face;			/* Flag, nz if simplex lies on cube surface */
	int icomb[MXDI];    /* icomb[] specifies the equation to convert simplex space */
						/* coordinates back into cartesian space. */
	                    /* Index by Absolute[di] -> Simplex Parameter[sdi], */
						/*                          -1 == value 0, -2 == value 1 */
						/* Note that many Abs can map to one Param to form a sum. */
	int offs[MXDI+1];	/* Offsets to simplex verticies within cube == bit per dim */
	int goffs[MXDI+1];	/* Offsets to simplex verticies within grid */
	int foffs[MXDI+1];	/* Fwd grid floating offsets to simplex verticies from cube base */
	int pmino[MXDI], pmaxo[MXDI]; /* Cube verticy offsets to setup simplex pmin[] and */
						/* pmax[] bounding box pointers. */
} psxinfo;

/* Sub simplexes of a cube information structure */
typedef struct {
	int sdi;			/* Sub-simplex dimensionality */
	int nospx;			/* Number of sub-simplexs per cube */
	psxinfo *spxi;		/* Per sub-simplex info array, NULL if not initialised */
} ssxinfo;

/* - - - - - - - - - - - - - - - - - - - - - */
/* NOTE :- This should really be re-arranged to be per-sub-simplex caching, */
/* rather than fxcell caching. Rather than stash the simplex info in the fxcells, */
/* create a separate cache or some other way of sharing the common simplexes. */
/* The code that ignores common face simplexes in fxcells could then be removed (?). */

/* Simplex definition. Each top level fwd interpolation cell, */
/* is decomposed into sub-simplexes. Sub-simplexes are of equal or */
/* lower dimensionality (ie. faces, edges, verticies) to the cube. */
struct _simplex {
	int refcount;				/* reference count */
	struct _rspl *s;			/* Pointer to parent rspl */
	int ix;						/* Construction Fwd cell index */
	int si;						/* Diagnostic - simplex number within level */
	int sdi;					/* Sub-simplex dimensionality */
	int efdi;					/* Effective fdi. This will be = fdi for a non clip */
								/* plane simplex, and fdi+1 for a clip plane simplex. */
								/* The DOF (Degress of freedom) within this simplex = sdi - efdi */

	psxinfo *psxi;				/* Generic per simplex info (construction cube relative) */
	int vix[MXRI+1];			/* fwd cell vertex indexes of this simplex [sdi+1] */
								/* This is a universal identification of this simplex. */
	struct _simplex *hlink;		/* Link to other cells with this hash */
	unsigned int touch;			/* Last touch count. */
	short flags;				/* Various flags */

#define SPLX_CLIPSX  0x01		/* This is a clip plane simplex */

#define SPLX_FLAG_1  0x04		/* v, linmin/max  initialised */
#define SPLX_FLAG_2  0x08		/* lu/svd initialised */
#define SPLX_FLAG_2F 0x10		/* lu/svd init. failed */
#define SPLX_FLAG_4  0x20		/* locus found */
#define SPLX_FLAG_5  0x40		/* auxiliary lu/svd initialised */
#define SPLX_FLAG_5F 0x80		/* auxiliary lu/svd init. failed */

#define SPLX_FLAGS  (SPLX_FLAG_1 | SPLX_FLAG_2 | SPLX_FLAG_2F \
                   | SPLX_FLAG_4 | SPLX_FLAG_5 | SPLX_FLAG_5F)

	double v[MXRI+1][MXRO+1]; 	/* Simplex Vertex values */
								/* v[0..sdi][0..fdi-1] are the output interpolation values */
								/* v[0..sdi][fdi] are the ink limit interpolation values */

								/* Baricentric vv[x][y] = (v[y][x] - v[y+1][x]) */
								/* and         vv[x][sdi] = v[sdi][x]           */

								/* Note that #num indicates appropriate flag number */
								/* and *num indicates a validator */

	double p0[MXRI]; 			/* Simplex base position = construction cube p[0] */
	double pmin[MXRI];			/* Simplex vertex input space min and */
	double pmax[MXRI];			/* max values [di] */

	double min[MXRO+1], max[MXRO+1]; /* Simplex vertex output space [fdi+1] and */
								/* ink limit bounding values at minmax[fdi] */

	/* If sdi == efdi, this holds the LU decomposition, */
	/* else this holds the SVD and solution locus info */

	char *aloc2;		/* Memory allocation for #2 & #4 */

	/* double **d_u;	   LU decomp of vv, U[0..efdi-1][0..sdi-1]		#2 */
	/* int *d_w;		   LU decomp of vv, W[0..sdi-1]					#2 */

	double **d_u;		/* SVD decomp of vv, U[0..efdi-1][0..sdi-1]		#2 */
	double *d_w;		/* SVD decomp of vv, W[0..sdi-1]				#2 */
	double **d_v;		/* SVD decomp of vv, V[0..sdi-1][0..sdi-1]		#2 */

						/* Degrees of freedom = dof = sdi - efdi */
	double **lo_l;		/* Locus coefficients,  [0..sdi-1][0..dof-1]	#2 */

	double *lo_xb;		/* RHS used to compute lo_bd [0..efdi-1]		*4 */
	double *lo_bd;		/* Locus base solution, [0..sdi-1]				#4 */

	unsigned auxbm;		/* aux bitmap mask for ax_lu and ax_svd 		*5 */
	int      aaux;		/* naux count for allocation					*5 */
	int      naux;		/* naux for calculation (may be < aaux ?)		*5 */

	/* if (sdi-efdi = dof) == naux this holds LU of lo_l */
	/* else this holds the SVD of lo_l */

	char *aloc5;		/* Memory allocation for #5 */

	/* double **ax_u;	   LU decomp of lo_l							#5 */
	/* int *ax_w;		   Pivot record for ax_lu decomp				#5 */

	double **ax_u;		/* SVD decomp of lo_l, U[0..naux-1][0..dof-1]	#5 */
	double *ax_w;		/* SVD decomp of lo_l, W[0..dof-1]				#5 */
	double **ax_v;		/* SVD decomp of lo_l, V[0..dof-1][0..dof-1]	#5 */

}; typedef struct _simplex simplex;

/* A candidate search (fwd) fxcell (cell cache entry structure) */
struct _fxcell {
	struct _rspl *s;		/* Pointer to parent rspl */

	/* Cache information */
	int ix;					/* Corresponding fwd cell index */
	struct _fxcell *hlink;	/* Link to other cells with this hash */
	struct _fxcell *mrudown;/* Links to next most recently used fxcell */
	struct _fxcell *mruup;
	int refcount;			/* Reference count */
	int flags;				/* Non-zero if the fxcell has been initialised */
#define CELL_FLAG_1  0x01	/* Basic initialisation, including nn_grp */
#define CELL_FLAG_2  0x02	/* Simplex information initialised */

	/* Use information */
	double sort;			/* Sort key */
	double limmin, limmax;	/* limit() min/max for this fxcell */

	/* Quick nn distance information */
	nn_grp g;	

//	double bcent[MXRO];		/* Output value bounding shere center */
//	double brad;			/* Output value bounding shere radius */
//	double bradsq;			/* Output value bounding shere radius squared */
//	double wbrad;			/* Output value weighted bounding shere radius */

	double p[POW2MXRI][MXRI]; /* Vertex input positions for this cube. */
							/* Copied to x->pmin/pmax[] & ink limit */

	double v[POW2MXRI][MXRO+1]; /* Vertex data for this cube. Copied to x->v[] */
							/* v[][fdi] is the ink limit values, if relevant */

	simplex **sx[MXRI+1];	/* Lists of simplexes that make up this fxcell. */
							/* Each list indexed by the non-limited simplex */
							/* dimensionality (similar to sspxi[]) */
							/* Each list will be NULL if it hasn't been created yet */
	int sxno[MXRI+1];		/* Corresponding count of each list */

}; typedef struct _fxcell fxcell;

/* surface bxcell sl status */
typedef enum {
	bx_uninit  = 0,		/* sl is not initialised */
	bx_filled  = 1,		/* sl has been filled with initial fwd cell vertexes */
	bx_rethinnd = 2,	/* sl vertexes need to be re-thinned */
	bx_thinned = 3,		/* sl vertexes have been thinned */
	bx_conv    = 4		/* sl vertexes have been converted to fwd cell indexes */
} bxstat;

/* Structure to hold bwd cell information for surface list, and also */
/* for seed fill bwd cell propogation. (Cells on surface will have two */
/* of these) */
struct _bxcell{
	int ix;						/* nnrev[] index of this bwd cell */
	int gc[MXRO];				/* coordinate of this bwd cell */
	nn_grp g;					/* Group nn calculation info */
	struct _bxcell *ss;			/* bwd surface cell to start search from */
	double sdist;				/* Est. wtd distance from this cell to ss */
	int tix;					/* target rev[] index being filled (visited check) */

	bxstat status;				/* State of sl list */					
	int *sl;					/* fwd vertex seed list for surface bxcells */
	int *dl;					/* deleted fwd vertex list for this bxcell */

	int *scell;					/* If non-NULL, this is a (non-surface) */
								/* seeding super bxcell, and scell contains */
								/* the list of bxcells covered */ 

	struct _bxcell *slist;		/* Linked list of all surface bxcells */
	struct _bxcell *hlink;		/* Linked list of surface bxcells with same ix hash */
	struct _bxcell *xlist;		/* Linked list of surface exploration search region */
	double emin;				/* estimated minimum wtd distance of this cell in current search */
	struct _bxcell *tlist;		/* Linked list of solution surface cells for current search. */

	struct _bxcell *flist;		/* Linked list for nnrev[] fill seeds */

	double cc;					/* Distance of group from gamut center */
	double dw;					/* Delta width from ocent of furthest point */
	struct _bxcell *wlist;		/* Linked list for shadow bxcells during thinning */

	int    debug;				/* debug flag - for VRML tagging */
}; typedef struct _bxcell bxcell;

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


/* Enough space is needed to cache all the fxcells/simplexes */
/* for a full aux. locus, or the query will be processed in */
/* several "chunks" and will be slower. */
/* This sets the basic memory usage of the rev code. */

#define REV_ACC_GRES_MUL 2.0		/* 2.0 Reverse accelleration grid resolution */
									/* multiplier over fwd grid resolution */
#define REV_ACC_GRES_LIMIT 43		/* Reverse accelln. grid resolution limit before env. mult. */
#define REV_MAX_MEM_RATIO 0.3		/* 0.3 Proportion of first 1G of Ram to use */
#define REV_MAX_MEM_RATIO2 0.4		/* 0.4 Proportion of rest of Ram to use */
									/* rev as a fraction of the System RAM. */
#define HASH_FILL_RATIO 3			/* 3 Ratio of entries to hash size */

/* The structure where fxcells and simplexes are allocated and cached. */

/* Holds the fxcell and simplex match cache specific information */
typedef struct {
	struct _rspl *s;			/* Pointer to parent rspl */
	int nacells;				/* Number of allocated cells */
	int nunlocked;				/* Number of unlocked cells that could be freed */
	int cell_hash_size;			/* Current size of cell hash list */
	fxcell **hashtop;			/* Top of hash list [cell_hash_size] */
	fxcell *mrutop, *mrubot;	/* Top and bottom pointers of mru list */
								/* that tracks allocated fxcells */

	int spx_hash_size;			/* Current size of shared face simplex hash list */
	simplex **spxhashtop;		/* Shared face simplex hash index list */
	int nspx;					/* Number of simplexes in hash list */
} revcache;

/* common search information */
/* Type of (internal) reverse search */
enum ops {
	exact  = 0,		/* Search for all input values that exactly map to given output value */
	clipv  = 1,		/* Search for input values that map to outermost solution along a vector */
	clipn  = 2,		/* Search for input values that map to closest solution */
	auxil  = 3,		/* Search for input values that map to given output, and closest to auxiliary target */
	locus  = 4 };	/* Return range of auxiliary values that contains solution */

/* + possible exact with clip */

/* Structure to hold clip line state information */
typedef struct {
	struct _rspl *s;	/* Pointer to parent rspl */
	double st[MXRO];	/* start of line - reverse grid base value */
	double de[MXRO];	/* direction of line */
	int    di[MXRO];	/* incerement in line direction */
	int    ci[MXRO];	/* current rev grid coordinate */
	double t;			/* Parameter 0.0 - 1.0, line finished if t > 1.0 */
} line;


/* Structure to hold aux value of an intersection of a */
/* solution locus with a sub-simplex. Used when asegs flag is set */
typedef struct {
	double xval;		/* Auxiliary value */
	int nv;				/* Number of verticies valid */
	int vix[MXRI+1];	/* Verticy indexes of sub-simplex involved */
} axisec;

/* -------------------------------------------- */
/* Information needed/cached for reverse lookup */
struct _schbase {
	struct _rspl *s;	/* Pointer to parent rspl */

	int flags;			/* Hint flags */
	enum ops op;		/* Type of reverse search operation */
	int ixc;			/* Cube index of corner that holds maximum input values */

	int snsdi, ensdi;	/* Start and end extra sub-simplex dimensionality */
	int (*setsort)(struct _schbase *b, fxcell *c);	/* Function to triage & set cube sort index */
	int (*check)(struct _schbase *b, fxcell *c);		/* Function to recheck cube worth keeping */
	int (*compute)(struct _schbase *b, simplex *x);	/* Function to compute a simplex solution */

	double v[MXRO+1];	/* Target output value, + ink limit */
	double av[MXRI];	/* Target auxiliary values */
	int auxm[MXRI];		/* aux mask flags */
	unsigned auxbm;		/* aux bitmap mask */
	int naux;			/* Number of auxiliary target input values */
	int auxi[MXRI];		/* aux list of auxiliary target input values */
	double idist;		/* best input distance auxiliary target found (smaller is better) */
	int iabove;			/* Number of auxiliaries at or above zero */

	int canvecclip;		/* Non-zero if vector clip direction usable */
	double cdir[MXRO];	/* Clip vector direction and length wrt. v[] */
	double ncdir[MXRO];	/* Normalised clip vector */
	double **cla;		/* Clip vector LHS implicit equation matrix [fdi][fdi+1] (inc. ink tgt.) */
	double clb[MXRO+1];	/* Clip vector RHS implicit equation vector [fdi+1] (inc. ink tgt.) */
	double cdist;		/* Best clip locus distance found (aim is min +ve) :- weighted for nn */
	int    iclip;		/* NZ if result is above (disabled) ink limit */

	int mxsoln;			/* Maximum number of solutions that we want */
	int nsoln;			/* Current number of solutions found */
	co *cpp;			/* Store solutions here */

	int   lxi;			/* Locus search axiliary index */
	double min, max;	/* current extreme locus values for locus search */
	int    asegs;		/* flag - find all search segments */
	int    axisln;		/* Number of elements used in axisl[] */
	int    axislz;		/* Space allocated to axisl[] */
	axisec *axisl;		/* Auxiliary intersections */

	int lclistz;		/* Allocated space to fxcell sort list */
	fxcell **lclist;	/* Sorted list of pointers to candidate fxcells */

	int pauxcell;		/* Indexe of previous call solution fxcell, -1 if not relevant */
	int plmaxcell;		/* Indexe of previous call solution fxcell, -1 if not relevant */
	int plmincell;		/* Indexe of previous call solution fxcell, -1 if not relevant */

	int lsxfilt;		/* Allocated space of simplex filter list */
	char *sxfilt;		/* Flag for simplexes that should be in an fxcell */

	int rix;			/* Diagnostic - rev[] or nnrev[] index for this point */

}; typedef struct _schbase schbase;

/* ----------------------------------------- */

#ifdef STATS
struct _stats {
	int 	searchcalls;/* Number of top level searches */
	int 	csearched;	/* Cells searched */
	int 	ssearched;	/* Simplexes searched */
	int		sinited;	/* Simplexes initialised to base level */
	int		sinited2a;	/* Simplexes initialised to 2nd level with LU */
	int		sinited2b;	/* Simplexes initialised to 2nd level with SVD */
	int		sinited4i;	/* Simplexes invalidated at 4th level */
	int		sinited4;	/* Simplexes initialised to 4th level */
	int		sinited5i;	/* Simplexes invalidated at 5th level */
	int		sinited5a;	/* Simplexes initialised to 5th level with LU */
	int		sinited5b;	/* Simplexes initialised to 5th level with SVD */
	int		chits;		/* Cells hit in cache */
	int		cmiss;		/* Cells misses in cache */
}; typedef struct _stats stats;
#endif /* STATS */

/* ----------------------------------------- */
/* Reverse info stored in main rspl function */
struct _rev_struct {

	/* First section, basic information that doesn't change */
	/* Has been initialised if inited != 0 */

	int inited;			/* Non-zero if first section has been initialised */
						/* All other sections depend on this. */
	int fastsetup;		/* Flag - NZ if fast setup at cost of slow throughput */

	int lchweighted;	/* Non-zero if nearest search is LCh weighted */
	double lchw[MXRO];	/* LCh nearest weighting */
	double lchw_sq[MXRO];	/* LCh nearest weighting squared */
	double lchw_chsq;		/* lchw_sq[1] - lchw_sq[2] */

	struct _rev_struct *next;	/* Linked list of global instances sharing memory */
	size_t max_sz;		/* Maximum size permitted */
	size_t sz;			/* Total memory current allocated by rev */

#ifdef NEVER
	int thissz, lastsz;	/* Debug reporting */
#endif

	/* Reverse grid lookup information */
	int res;			/* Reverse grid resolution == ncells on a side */
	int no;				/* Total number of points in reverse grid = rev.ares ^ fdi */
	int coi[MXRO];		/* Coordinate increments for each dimension */
	int hoi[1 << MXRO];	/* Coordinate increments for progress through cube */
	datao gl,gh,gw;		/* Reverse grid low, high, grid bwd cell width */

	/* Second section, accelleration information that changes with ink limit. */
	int rev_valid;		/* nz if this information in rev[] and nnrev[] is valid */
	int **rev;			/* Exact reverse lookup starting list. */
	int **nnrev;		/* Nearest reverse lookup starting list. */
						/* These lists are of fwd grid indexes. */
						/* [0] is allocation length */
						/* [1] is the next free entry index (length + 3, not counting -1) */
						/* [2] is index into share lists, -1 if not shared. */
						/* Then follows cube indexes */
						/* Last entry is marked with -1 */

	double ocent[MXRO];	/* rev cell gamut "center" point for thinning and shadow testing. */

	bxcell *surflist;	/* Linked list of rev[] bwd cells that contain gamut surface fwd cells. */
						/* Used to speed up fill_nncell() when rev.fastsetup is set, else NULL */
	int surf_hash_size;	/* Current size of bxcell hash list */
	bxcell **surfhash;	/* bxcell hash index list */

	int **sharelist;	/* Array of pointers to shared (fwd grid list sharer) records. */ 
						/* Each record is same format as rev[]/nnrev[], except */
						/* [2] is used to detect scanning the same list. */
	int sharellen;		/* Size of sharelist */ 
	int sharelaloc;		/* Allocation size of sharelist */ 

	/* Third section */
	revcache *cache;	/* Where fxcells and simplexes are allocated and cached */
	/* Sub-dimension simplex information */
	ssxinfo sspxi[MXRI+1];/* One per sub dimenstionality at offset sdi */

	/* Fourth section */
	/* Has been initialise if sb != NULL */
	schbase *sb;		/* Structure holding calculated per-search call information */

	unsigned int stouch; /* Simplex touch count to avoid searching shared simplexs twice */
#ifdef STATS
	stats st[5];	/* Set of stats info indexed by enum ops */
#endif	/* STATS */

	int primsecwarn;	/* Not primary or secondary warning has been issued */

#ifdef CHECK_NNLU
int cknn_no;			/* Number checked */
double cknn_we;			/* Worst DE */
int cknn_noerrs;		/* Number not as good as closet vertex */
int cknn_nobsb;			/* Number not as good as second closest vertex */
int cknn_nonis;			/* Number not in surface */
#endif /* CHECK_NNLU */

}; typedef struct _rev_struct rev_struct;


/* ------------------------------------ */
/* Utility functions used by other parts of rspl implementation */

/* Initialise a static sub-simplex verticy information table */
void rspl_init_ssimplex_info(struct _rspl *s,	/* RSPL object */
ssxinfo *xip,				/* Pointer to sub-simplex info structure to init. */
int sdi);					/* Sub-simplex dimensionality (range 0 - di) */

/* Free the given sub-simplex verticy information */
void rspl_free_ssimplex_info(struct _rspl *s,
ssxinfo *xip);		/* Pointer to sub-simplex info structure */

#endif /* RSPL_REV_H */