summaryrefslogtreecommitdiff
path: root/rspl/rev.h
blob: b947cbe59a8dca24a1a3a0d09e5dbb08ee2e0541 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
#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 coordinates D0, D1, D2
 * 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:
 *
 * D1 D0 D2		Smallest -> Largest cartesian sort
 * P2 P1 P0 	Corresponding Parameter coordinates (note reverse order!)
 * 
 * B0 = P0		Conversion to Baricentric coordinates
 * B1 = P1 - P0
 * B2 = P2 - P1
 * B3 = 1  - P2
 *
 * The vertex values directly correspond to Baricentric 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
 *
 * It is this which is used in rev.c for solving the reverse
 * interpolation problem.
 */

/* A structure to hold per simplex coordinate and vertex mapping */
/* 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];	/* 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. */
						/* icomb[] specifies the equation to convert simplex space */
						/* coordinates back into cartesian space. */
	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 cell caching. Rather than stash the simplex info in the cells, */
/* create a separate cache or some other way of sharing the common simplexes. */
/* The code that ignores common face simplexes in cells 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) withing 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 cell (cell cache entry structure) */
struct _cell {
	struct _rspl *s;		/* Pointer to parent rspl */

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

	/* Use information */
	double sort;			/* Sort key */

	double limmin, limmax;	/* limit() min/max for this cell */
    double bcent[MXRO];		/* Output value bounding shere center */
    double brad;			/* Output value bounding shere radius */
    double bradsq;			/* Output value bounding shere radius squared */

	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 cell. */
							/* 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 _cell cell;

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


/* Enough space is needed to cache all the cells/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 cells are allocated and cached. */

/* Holds the cell 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 */
	cell **hashtop;				/* Top of hash list [cell_hash_size] */
	cell *mrutop, *mrubot;		/* Top and bottom pointers of mru list */
								/* that tracks allocated cells */

	int spx_hash_size;			/* Current size of simplex hash list */
	simplex **spxhashtop;		/* 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, cell *c);	/* Function to triage & set cube sort index */
	int (*check)(struct _schbase *b, cell *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) */
	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 cell sort list */
	cell **lclist;		/* Sorted list of pointers to candidate cells */

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

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

	double crad;			/* nn current radius distance */
	double bw;				/* nn current window distance */
	int wex[MXRO * 2];		/* nn current window edge indexes */
	double wed[MXRO * 2];	/* nn current window edge distances */

}; 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 */

	struct _rev_struct *next;	/* Linked list of 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 ares;			/* Reverse grid allocated resolution, = res + 1, */
						/* allows for extra row used during construction */
	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 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. */
						/* rev.no pointers to lists of fwd grid indexes. */
						/* First int is allocation length */
						/* Second int is reference count */
						/* Then follows cube indexes */
						/* Last int is -1 */
	int **nnrev;		/* Nearest reverse lookup starting list. */
						/* rev.no pointers to lists of fwd grid indexes. */
						/* [0] is allocation length */
						/* [1] is the next free entry index */
						/* [2] is reference count */
						/* Then follows cube indexes */
						/* The last entry is marked with -1 */

	/* Third section */
	revcache *cache;	/* Where cells 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 */

}; 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 */