1 /*!
2  * \file lib/raster/quant.c
3  *
4  * \brief Raster Library - Quantization rules.
5  *
6  * The quantization table is stored as a linear array. Rules are added
7  * starting from index 0. Redundant rules are not eliminated. Rules
8  * are tested from the highest index downto 0. There are two
9  * "infinite" rules. Support is provided to reverse the order of the
10  * rules.
11  *
12  * (C) 1999-2009 by the GRASS Development Team
13  *
14  * This program is free software under the GNU General Public License
15  * (>=v2). Read the file COPYING that comes with GRASS for details.
16  *
17  * \author USACERL and many others
18  */
19 
20 #include <stdlib.h>
21 #include <grass/gis.h>
22 #include <grass/raster.h>
23 
24 static int double_comp(const void *, const void *);
25 
26 #define USE_LOOKUP 1
27 #define MAX_LOOKUP_TABLE_SIZE 2048
28 #define NO_DATA (Rast_set_c_null_value (&tmp, 1), (CELL) tmp)
29 
30 #undef MIN
31 #undef MAX
32 #define MIN(a,b) ((a) < (b) ? (a) : (b))
33 #define MAX(a,b) ((a) > (b) ? (a) : (b))
34 
35 #define NO_LEFT_INFINITE_RULE (! q->infiniteLeftSet)
36 #define NO_RIGHT_INFINITE_RULE (! q->infiniteRightSet)
37 #define NO_FINITE_RULE (q->nofRules <= 0)
38 #define NO_EXPLICIT_RULE (NO_FINITE_RULE && \
39 			  NO_LEFT_INFINITE_RULE && NO_RIGHT_INFINITE_RULE)
40 
41 /*!
42    \brief Resets the number of defined rules and number of infinite rules to 0
43 
44    \param q pointer to Quant structure to be reset
45  */
Rast_quant_clear(struct Quant * q)46 void Rast_quant_clear(struct Quant *q)
47 {
48     q->nofRules = 0;
49     q->infiniteRightSet = q->infiniteLeftSet = 0;
50 }
51 
52 /*!
53    \brief Resets and frees allocated memory
54 
55    Resets the number of defined rules to 0 and free's space allocated
56    for rules. Calls Rast_quant_clear().
57 
58    \param q pointer to Quant structure to be reset
59  */
Rast_quant_free(struct Quant * q)60 void Rast_quant_free(struct Quant *q)
61 {
62     Rast_quant_clear(q);
63 
64     if (q->maxNofRules > 0)
65 	G_free(q->table);
66     if (q->fp_lookup.active) {
67 	G_free(q->fp_lookup.vals);
68 	G_free(q->fp_lookup.rules);
69 	q->fp_lookup.nalloc = 0;
70 	q->fp_lookup.active = 0;
71     }
72     q->maxNofRules = 0;
73 }
74 
75 /*!
76  * \brief Organized fp_lookup table.
77  *
78  *  Organizes fp_lookup table for faster (logarithmic) lookup time
79  *  G_quant_organize_fp_lookup() creates a list of min and max for
80  *  each quant rule, sorts this list, and stores the pointer to quant
81  *  rule that should be used inbetween any 2 numbers in this list.
82  *  Also it stores extreme points for 2 infinite rules, if exist.
83  *  After the call to G_quant_organize_fp_lookup()
84  *  instead of linearly searching through list of rules to find
85  *  a rule to apply, quant lookup will perform a binary search
86  *  to find an interval containing floating point value, and then use
87  *  the rule associated with this interval.
88  *  when the value doesn't fall within any interval, check for the
89  *  infinite rules.
90  *
91  * \param q pointer to Quant structure which holds quant rules info
92  *
93  * \return 1 on success
94  */
Rast__quant_organize_fp_lookup(struct Quant * q)95 int Rast__quant_organize_fp_lookup(struct Quant *q)
96 {
97     int i;
98     DCELL val;
99     CELL tmp;
100     struct Quant_table *p;
101 
102     if (q->nofRules * 2 > MAX_LOOKUP_TABLE_SIZE)
103 	return -1;
104     if (q->nofRules == 0)
105 	return -1;
106     q->fp_lookup.vals = (DCELL *)
107 	G_calloc(q->nofRules * 2, sizeof(DCELL));
108     /* 2 endpoints for each rule */
109     q->fp_lookup.rules = (struct Quant_table **)
110 	G_calloc(q->nofRules * 2, sizeof(struct Quant_table *));
111 
112     /* first we organize finite rules into a table */
113     if (!NO_FINITE_RULE) {
114 	i = 0;
115 	/* get the list of DCELL values from set of all dLows and dHighs
116 	   of all rules */
117 	/* NOTE: if dLow==DHigh in a rule, the value appears twice in a list
118 	   but if dLow==DHigh of the previous, rule the value appears only once */
119 
120 	for (p = &(q->table[q->nofRules - 1]); p >= q->table; p--) {
121 	    /* check if the min is the same as previous maximum */
122 	    if (i == 0 || p->dLow != q->fp_lookup.vals[i - 1])
123 		q->fp_lookup.vals[i++] = p->dLow;
124 	    q->fp_lookup.vals[i++] = p->dHigh;
125 	}
126 	q->fp_lookup.nalloc = i;
127 
128 	/* now sort the values */
129 	qsort((char *)q->fp_lookup.vals, q->fp_lookup.nalloc,
130 	      sizeof(DCELL), double_comp);
131 
132 	/* now find the rule to apply inbetween each 2 values in a list */
133 	for (i = 0; i < q->fp_lookup.nalloc - 1; i++) {
134 	    /*debug
135 	       fprintf (stderr, "%lf %lf ", q->fp_lookup.vals[i], q->fp_lookup.vals[i+1]);
136 	     */
137 	    val = (q->fp_lookup.vals[i] + q->fp_lookup.vals[i + 1]) / 2.;
138 	    q->fp_lookup.rules[i] =
139 		Rast__quant_get_rule_for_d_raster_val(q, val);
140 	    /* debug
141 	       if(q->fp_lookup.rules[i])
142 	       fprintf (stderr, "%lf %lf %d %d\n", q->fp_lookup.rules[i]->dLow, q->fp_lookup.rules[i]->dHigh, q->fp_lookup.rules[i]->cLow, q->fp_lookup.rules[i]->cHigh);
143 	       else fprintf (stderr, "null\n");
144 	     */
145 
146 	}
147     }				/* organizing finite rules */
148 
149     if (!NO_LEFT_INFINITE_RULE) {
150 	q->fp_lookup.inf_dmin = q->infiniteDLeft;
151 	q->fp_lookup.inf_min = q->infiniteCLeft;
152     }
153     else {
154 	if (q->fp_lookup.nalloc)
155 	    q->fp_lookup.inf_dmin = q->fp_lookup.vals[0];
156 	q->fp_lookup.inf_min = NO_DATA;
157     }
158 
159     if (!NO_RIGHT_INFINITE_RULE) {
160 	if (q->fp_lookup.nalloc)
161 	    q->fp_lookup.inf_dmax = q->infiniteDRight;
162 	q->fp_lookup.inf_max = q->infiniteCRight;
163     }
164     else {
165 	q->fp_lookup.inf_dmax = q->fp_lookup.vals[q->fp_lookup.nalloc - 1];
166 	q->fp_lookup.inf_max = NO_DATA;
167     }
168     q->fp_lookup.active = 1;
169     return 1;
170 }
171 
172 /*!
173  * \brief Initialize the structure
174  *
175  * Initializes the <i>q</i> struct.
176  *
177  * \param quant pointer to Quant structure to be initialized
178  */
Rast_quant_init(struct Quant * quant)179 void Rast_quant_init(struct Quant *quant)
180 {
181     quant->fp_lookup.active = 0;
182     quant->maxNofRules = 0;
183     quant->truncate_only = 0;
184     quant->round_only = 0;
185     Rast_quant_clear(quant);
186 }
187 
188 /*!
189    \brief Returns whether or not quant rules are set to truncate map
190 
191    \param quant pointer to Quant structure which holds quant rules info
192 
193    \return 1 if truncate is enable
194    \return 0 if not truncated
195  */
Rast_quant_is_truncate(const struct Quant * quant)196 int Rast_quant_is_truncate(const struct Quant *quant)
197 {
198     return quant->truncate_only;
199 }
200 
201 /*!
202    \brief  Returns whether or not quant rules are set to round map
203    \param quant pointer to Quant structure which holds quant rules info
204 
205    \return 1 is round
206    \return 0 not round
207  */
Rast_quant_is_round(const struct Quant * quant)208 int Rast_quant_is_round(const struct Quant *quant)
209 {
210     return quant->round_only;
211 }
212 
213 /*!
214  * \brief Sets the quant rules to perform simple truncation on floats.
215  *
216  * Sets the quant for <i>q</i> rules to perform simple truncation on
217  * floats.
218  *
219  * \param quant pointer to Quant structure which holds quant rules info
220  */
Rast_quant_truncate(struct Quant * quant)221 void Rast_quant_truncate(struct Quant *quant)
222 {
223     quant->truncate_only = 1;
224 }
225 
226 /*!
227  * \brief Sets the quant rules to perform simple rounding on floats.
228  *
229  * Sets the quant for <i>q</i> rules to perform simple rounding on
230  * floats.
231  *
232  * \param quant pointer to Quant structure which holds quant rules info
233  */
Rast_quant_round(struct Quant * quant)234 void Rast_quant_round(struct Quant *quant)
235 {
236     quant->round_only = 1;
237 }
238 
quant_set_limits(struct Quant * q,DCELL dLow,DCELL dHigh,CELL cLow,CELL cHigh)239 static void quant_set_limits(struct Quant *q,
240 			     DCELL dLow, DCELL dHigh, CELL cLow, CELL cHigh)
241 {
242     q->dMin = dLow;
243     q->dMax = dHigh;
244     q->cMin = cLow;
245     q->cMax = cHigh;
246 }
247 
quant_update_limits(struct Quant * q,DCELL dLow,DCELL dHigh,CELL cLow,DCELL cHigh)248 static void quant_update_limits(struct Quant *q,
249 				DCELL dLow, DCELL dHigh,
250 				CELL cLow, DCELL cHigh)
251 {
252     if (NO_EXPLICIT_RULE) {
253 	quant_set_limits(q, dLow, dHigh, cLow, cHigh);
254 	return;
255     }
256 
257     q->dMin = MIN(q->dMin, MIN(dLow, dHigh));
258     q->dMax = MAX(q->dMax, MAX(dLow, dHigh));
259     q->cMin = MIN(q->cMin, MIN(cLow, cHigh));
260     q->cMax = MAX(q->cMax, MAX(cLow, cHigh));
261 }
262 
263 /*!
264  * \brief Returns the minimum and maximum cell and dcell values of all
265  *  the ranges defined.
266  *
267  * Extracts the minimum and maximum floating-point and integer values
268  * from all the rules (except the "infinite" rules) in <i>q</i> into
269  * <i>dmin</i>, <i>dmax</i>, <i>cmin</i>, and <i>cmax</i>.
270  *
271  * \param quant pointer to Quant structure which holds quant rules info
272  * \param[out] dmin minimum fp value
273  * \param[out] dmax maximum fp value
274  * \param[out] cmin minimum value
275  * \param[out] cmax maximum value
276  *
277  * \return -1 if q->truncate or q->round are true or after
278  * Rast_quant_init (), or any call to Rast_quant_clear () or Rast_quant_free()
279  * no explicit rules have been added. In this case the returned
280  * minimum and maximum CELL and DCELL values are null.
281  * \return 1 if there are any explicit rules
282  * \return 0 if there are no explicit rules (this includes cases when
283  * q is set to truncate or round map), and sets <i>dmin</i>,
284  * <i>dmax</i>, <i>cmin</i>, and <i>cmax</i> to NULL.
285  */
Rast_quant_get_limits(const struct Quant * q,DCELL * dMin,DCELL * dMax,CELL * cMin,CELL * cMax)286 int Rast_quant_get_limits(const struct Quant *q,
287 			  DCELL * dMin, DCELL * dMax, CELL * cMin,
288 			  CELL * cMax)
289 {
290     if (NO_EXPLICIT_RULE) {
291 	Rast_set_c_null_value(cMin, 1);
292 	Rast_set_c_null_value(cMax, 1);
293 	Rast_set_d_null_value(dMin, 1);
294 	Rast_set_d_null_value(dMax, 1);
295 	return -1;
296     }
297 
298     *dMin = q->dMin;
299     *dMax = q->dMax;
300     *cMin = q->cMin;
301     *cMax = q->cMax;
302 
303     return 1;
304 }
305 
306 /*!
307    \brief Returns the number of quantization rules defined.
308 
309    This number does not include the 2 infinite intervals.
310 
311    \param q pointer to Quant structure which holds quant rules info
312 
313    \return number of quantization rules
314  */
Rast_quant_nof_rules(const struct Quant * q)315 int Rast_quant_nof_rules(const struct Quant *q)
316 {
317     return q->nofRules;
318 }
319 
320 /*!
321    \brief Returns the i'th quantization rule.
322 
323    For 0 <= i < Rast_quant_nof_rules(). A larger value for i means that
324    the rule has been added later.
325 
326    \param q pointer to Quant structure which holds quant rules info
327    \param i index
328    \param[out] dLow minimum fp value
329    \param[out] dHigh maximum fp value
330    \param[out] cLow minimum value
331    \param[out] cHigh maximum value
332  */
Rast_quant_get_ith_rule(const struct Quant * q,int i,DCELL * dLow,DCELL * dHigh,CELL * cLow,CELL * cHigh)333 void Rast_quant_get_ith_rule(const struct Quant *q,
334 			     int i,
335 			     DCELL * dLow, DCELL * dHigh,
336 			     CELL * cLow, CELL * cHigh)
337 {
338     *dLow = q->table[i].dLow;
339     *dHigh = q->table[i].dHigh;
340     *cLow = q->table[i].cLow;
341     *cHigh = q->table[i].cHigh;
342 }
343 
quant_table_increase(struct Quant * q)344 static void quant_table_increase(struct Quant *q)
345 {
346     if (q->nofRules < q->maxNofRules)
347 	return;
348 
349     if (q->maxNofRules == 0) {
350 	q->maxNofRules = 50;
351 	q->table = (struct Quant_table *)
352 	    G_malloc(q->maxNofRules * sizeof(struct Quant_table));
353     }
354     else {
355 	q->maxNofRules += 50;
356 	q->table = (struct Quant_table *)
357 	    G_realloc((char *)q->table,
358 		      q->maxNofRules * sizeof(struct Quant_table));
359     }
360 }
361 
362 /*!
363    \brief Defines a rule for values "dLeft" and smaller.
364 
365    Values in this range are mapped to "c" if none of the "finite"
366    quantization rules applies.
367 
368    \param q pointer to Quant structure which holds quant rules info
369 
370    \param dLeft fp value
371    \param c value
372  */
Rast_quant_set_neg_infinite_rule(struct Quant * q,DCELL dLeft,CELL c)373 void Rast_quant_set_neg_infinite_rule(struct Quant *q, DCELL dLeft, CELL c)
374 {
375     q->infiniteDLeft = dLeft;
376     q->infiniteCLeft = c;
377     quant_update_limits(q, dLeft, dLeft, c, c);
378 
379     /* update lookup table */
380     if (q->fp_lookup.active) {
381 	q->fp_lookup.inf_dmin = q->infiniteDLeft;
382 	q->fp_lookup.inf_min = q->infiniteCLeft;
383     }
384     q->infiniteLeftSet = 1;
385 }
386 
387 /*!
388    \brief Returns in "dLeft" and "c" the rule values.
389 
390    For the negative infinite interval (see Rast_quant_set_neg_infinite_rule()).
391 
392    \param q pointer to Quant structure which holds quant rules info
393    \param[out] dLeft fp value
394    \param[out] c value
395 
396    \return 0 if this rule is not defined
397    \return 1 otherwise
398  */
Rast_quant_get_neg_infinite_rule(const struct Quant * q,DCELL * dLeft,CELL * c)399 int Rast_quant_get_neg_infinite_rule(const struct Quant *q,
400 				     DCELL * dLeft, CELL * c)
401 {
402     if (q->infiniteLeftSet == 0)
403 	return 0;
404 
405     *dLeft = q->infiniteDLeft;
406     *c = q->infiniteCLeft;
407 
408     return 1;
409 }
410 
411 /*!
412    \brief Defines a rule for values "dRight" and larger.
413 
414    Values in this range are mapped to "c" if none of the "finite"
415    quantization rules or the negative infinite rule applies.
416 
417    \param q pointer to Quant structure which holds quant rules info
418    \param dRight fp value
419    \param c value
420  */
Rast_quant_set_pos_infinite_rule(struct Quant * q,DCELL dRight,CELL c)421 void Rast_quant_set_pos_infinite_rule(struct Quant *q, DCELL dRight, CELL c)
422 {
423     q->infiniteDRight = dRight;
424     q->infiniteCRight = c;
425     quant_update_limits(q, dRight, dRight, c, c);
426 
427     /* update lookup table */
428     if (q->fp_lookup.active) {
429 	q->fp_lookup.inf_dmax = q->infiniteDRight;
430 	q->fp_lookup.inf_max = q->infiniteCRight;
431     }
432     q->infiniteRightSet = 1;
433 }
434 
435 /*!
436    \brief Returns in "dRight" and "c" the rule values.
437 
438    For the positive infinite interval (see Rast_quant_set_pos_infinite_rule()).
439 
440    \param q pointer to Quant structure which holds quant rules info
441    \param[out] dRight fp value
442    \param[out] c value
443 
444    \return 0 if this rule is not defined
445    \return 1 otherwise
446  */
Rast_quant_get_pos_infinite_rule(const struct Quant * q,DCELL * dRight,CELL * c)447 int Rast_quant_get_pos_infinite_rule(const struct Quant *q,
448 				     DCELL * dRight, CELL * c)
449 {
450     if (q->infiniteRightSet == 0)
451 	return 0;
452 
453     *dRight = q->infiniteDRight;
454     *c = q->infiniteCRight;
455 
456     return 1;
457 }
458 
459 /*!
460    \brief Adds a new rule to the set of quantization rules.
461 
462    If dLow < dHigh the rule will be stored with the low and high values
463    interchanged.
464 
465    Note: currently no cleanup of rules is performed, i.e. redundant
466    rules are not removed. This can't be changed because Categories
467    structure HEAVILY depends of quant rules stored in exactly the same
468    order they are entered. So if the cleanup or rearrangement is done in
469    the future make a flag for add_rule whether or not to do it, then
470    quant will not set this flag.
471 
472    \param q pointer to Quant structure which holds quant rules info
473    \param dLow minimum fp value
474    \param dHigh maximum fp value
475    \param cLow minimum value
476    \param cHigh maximum value
477  */
Rast_quant_add_rule(struct Quant * q,DCELL dLow,DCELL dHigh,CELL cLow,CELL cHigh)478 void Rast_quant_add_rule(struct Quant *q,
479 			 DCELL dLow, DCELL dHigh, CELL cLow, CELL cHigh)
480 {
481     int i;
482     struct Quant_table *p;
483 
484     quant_table_increase(q);
485 
486     i = q->nofRules;
487 
488     p = &(q->table[i]);
489     if (dHigh >= dLow) {
490 	p->dLow = dLow;
491 	p->dHigh = dHigh;
492 	p->cLow = cLow;
493 	p->cHigh = cHigh;
494     }
495     else {
496 	p->dLow = dHigh;
497 	p->dHigh = dLow;
498 	p->cLow = cHigh;
499 	p->cHigh = cLow;
500     }
501 
502     /* destroy lookup table, it has to be rebuilt */
503     if (q->fp_lookup.active) {
504 	G_free(q->fp_lookup.vals);
505 	G_free(q->fp_lookup.rules);
506 	q->fp_lookup.active = 0;
507 	q->fp_lookup.nalloc = 0;
508     }
509 
510     quant_update_limits(q, dLow, dHigh, cLow, cHigh);
511 
512     q->nofRules++;
513 }
514 
515 /*!
516    \brief Rreverses the order in which the qunatization rules are stored.
517 
518    See also Rast_quant_get_ith_rule() and Rast_quant_perform_d()).
519 
520    \param q pointer to Quant rules which holds quant rules info
521  */
Rast_quant_reverse_rule_order(struct Quant * q)522 void Rast_quant_reverse_rule_order(struct Quant *q)
523 {
524     struct Quant_table tmp;
525     struct Quant_table *pLeft, *pRight;
526 
527     pLeft = q->table;
528     pRight = &(q->table[q->nofRules - 1]);
529 
530     while (pLeft < pRight) {
531 	tmp.dLow = pLeft->dLow;
532 	tmp.dHigh = pLeft->dHigh;
533 	tmp.cLow = pLeft->cLow;
534 	tmp.cHigh = pLeft->cHigh;
535 
536 	pLeft->dLow = pRight->dLow;
537 	pLeft->dHigh = pRight->dHigh;
538 	pLeft->cLow = pRight->cLow;
539 	pLeft->cHigh = pRight->cHigh;
540 
541 	pRight->dLow = tmp.dLow;
542 	pRight->dHigh = tmp.dHigh;
543 	pRight->cLow = tmp.cLow;
544 	pRight->cHigh = tmp.cHigh;
545 
546 	pLeft++;
547 	pRight--;
548     }
549 }
550 
quant_interpolate(DCELL dLow,DCELL dHigh,CELL cLow,CELL cHigh,DCELL dValue)551 static CELL quant_interpolate(DCELL dLow, DCELL dHigh,
552 			      CELL cLow, CELL cHigh, DCELL dValue)
553 {
554     if (cLow == cHigh)
555 	return cLow;
556     if (dLow == dHigh)
557 	return cLow;
558 
559     return (CELL) ((dValue - dLow) / (dHigh - dLow) * (DCELL) (cHigh - cLow) +
560 		   (DCELL) cLow);
561 }
562 
less_or_equal(double x,double y)563 static int less_or_equal(double x, double y)
564 {
565     if (x <= y)
566 	return 1;
567     else
568 	return 0;
569 }
570 
less(double x,double y)571 static int less(double x, double y)
572 {
573     if (x < y)
574 	return 1;
575     else
576 	return 0;
577 }
578 
579 /*!
580  * \brief
581  *
582  *
583  * Returns a CELL category for the floating-point <i>value</i> based
584  * on the quantization rules in <i>q</i>. The first rule found that
585  * applies is used. The rules are searched in the reverse order they
586  * are added to <i>q</i>. If no rule is found, the <i>value</i>
587  * is first tested against the negative infinite rule, and finally
588  * against the positive infinite rule. If none of these rules apply,
589  * the NULL-value is returned.
590  *
591  * <b>Note:</b> See G_quant_organize_fp_lookup() for details on how
592  * the values are looked up from fp_lookup table when it is
593  * active. Right now fp_lookup is automatically organized during the
594  * first call to Rast_quant_get_cell_value().
595  *
596  * \param q pointer to Quant structure which holds quant rules info
597  * \param dcellValue fp cell value
598  *
599  * \return cell value (integer)
600  */
Rast_quant_get_cell_value(struct Quant * q,DCELL dcellVal)601 CELL Rast_quant_get_cell_value(struct Quant * q, DCELL dcellVal)
602 {
603     CELL tmp;
604     DCELL dtmp;
605     int try, min_ind, max_ind;
606     struct Quant_table *p;
607     int (*lower) ();
608 
609     dtmp = dcellVal;
610     /* I know the functions which call me already check for null values,
611        but I am a public function, and can be called from outside */
612     if (Rast_is_d_null_value(&dtmp))
613 	return NO_DATA;
614 
615     if (q->truncate_only)
616 	return (CELL) dtmp;
617 
618     if (q->round_only) {
619 	if (dcellVal > 0)
620 	    return (CELL) (dcellVal + .5);
621 	return (CELL) (dcellVal - .5);
622     }
623 
624     if (NO_EXPLICIT_RULE)
625 	return NO_DATA;
626     if (NO_EXPLICIT_RULE)
627 	return NO_DATA;
628 
629     if (USE_LOOKUP &&
630 	(q->fp_lookup.active || Rast__quant_organize_fp_lookup(q) > 0)) {
631 	/* first check if values fall within range */
632 	/* if value is below the range */
633 	if (dcellVal < q->fp_lookup.vals[0]) {
634 	    if (dcellVal <= q->fp_lookup.inf_dmin)
635 		return q->fp_lookup.inf_min;
636 	    else
637 		return NO_DATA;
638 	}
639 	/* if value is below above range */
640 	if (dcellVal > q->fp_lookup.vals[q->fp_lookup.nalloc - 1]) {
641 	    if (dcellVal >= q->fp_lookup.inf_dmax)
642 		return q->fp_lookup.inf_max;
643 	    else
644 		return NO_DATA;
645 	}
646 	/* make binary search to find which interval our value belongs to
647 	   and apply the rule for this interval */
648 	try = (q->fp_lookup.nalloc - 1) / 2;
649 	min_ind = 0;
650 	max_ind = q->fp_lookup.nalloc - 2;
651 	while (1) {
652 	    /* DEBUG
653 	       fprintf (stderr, "%d %d %d\n", min_ind, max_ind, try);
654 	     */
655 	    /* when the ruke for the interval is NULL, we exclude the end points.
656 	       when it exists, we include the end-points */
657 	    if (q->fp_lookup.rules[try])
658 		lower = less;
659 	    else
660 		lower = less_or_equal;
661 
662 	    if (lower(q->fp_lookup.vals[try + 1], dcellVal)) {	/* recurse to the second half */
663 		min_ind = try + 1;
664 		/* must be still < nalloc-1, since number is within the range */
665 		try = (max_ind + min_ind) / 2;
666 		continue;
667 	    }
668 	    if (lower(dcellVal, q->fp_lookup.vals[try])) {	/* recurse to the second half */
669 		max_ind = try - 1;
670 		/* must be still >= 0, since number is within the range */
671 		try = (max_ind + min_ind) / 2;
672 		continue;
673 	    }
674 	    /* the value fits into the interval! */
675 	    p = q->fp_lookup.rules[try];
676 	    if (p)
677 		return quant_interpolate(p->dLow, p->dHigh, p->cLow, p->cHigh,
678 					 dcellVal);
679 	    /* otherwise when finite rule for this interval doesn't exist */
680 	    else {		/* first check if maybe infinite rule applies */
681 		if (dcellVal <= q->fp_lookup.inf_dmin)
682 		    return q->fp_lookup.inf_min;
683 		if (dcellVal >= q->fp_lookup.inf_dmax)
684 		    return q->fp_lookup.inf_max;
685 		else
686 		    return NO_DATA;
687 	    }
688 	}			/* while */
689     }				/* looking up in fp_lookup */
690 
691     if (!NO_FINITE_RULE) {
692 	p = Rast__quant_get_rule_for_d_raster_val(q, dcellVal);
693 	if (!p)
694 	    return NO_DATA;
695 	return quant_interpolate(p->dLow, p->dHigh, p->cLow, p->cHigh,
696 				 dcellVal);
697     }
698 
699     if ((!NO_LEFT_INFINITE_RULE) && (dcellVal <= q->infiniteDLeft))
700 	return q->infiniteCLeft;
701 
702     if ((NO_RIGHT_INFINITE_RULE) || (dcellVal < q->infiniteDRight))
703 	return NO_DATA;
704 
705     return q->infiniteCRight;
706 }
707 
708 /*!
709    \brief Returns in "cell" the quantized CELL values.
710 
711    Returns in "cell" the quantized CELL values corresponding to the
712    DCELL values stored in "dcell". the number of elements quantized
713    is n. quantization is performed by repeated application of
714    Rast_quant_get_cell_value().
715 
716    \param q pointer to Quant structure which holds quant rules info
717    \param dcell pointer to fp cell values array
718    \param[out] cell pointer cell values array
719    \param n number of cells
720  */
Rast_quant_perform_d(struct Quant * q,const DCELL * dcell,CELL * cell,int n)721 void Rast_quant_perform_d(struct Quant *q,
722 			  const DCELL * dcell, CELL * cell, int n)
723 {
724     int i;
725 
726     for (i = 0; i < n; i++, dcell++)
727 	if (!Rast_is_d_null_value(dcell))
728 	    *cell++ = Rast_quant_get_cell_value(q, *dcell);
729 	else
730 	    Rast_set_c_null_value(cell++, 1);
731 }
732 
733 /*!
734    \brief Same as Rast_quant_perform_d(), except the type.
735 
736    \param q pointer to Quant structure which holds quant rules info
737    \param fcell pointer to fp cell values array
738    \param[out] cell pointer cell values array
739    \param n number of cells
740  */
Rast_quant_perform_f(struct Quant * q,const FCELL * fcell,CELL * cell,int n)741 void Rast_quant_perform_f(struct Quant *q,
742 			  const FCELL * fcell, CELL * cell, int n)
743 {
744     int i;
745 
746     for (i = 0; i < n; i++, fcell++)
747 	if (!Rast_is_f_null_value(fcell))
748 	    *cell++ = Rast_quant_get_cell_value(q, (DCELL) * fcell);
749 	else
750 	    Rast_set_c_null_value(cell++, 1);
751 }
752 
double_comp(const void * xx,const void * yy)753 static int double_comp(const void *xx, const void *yy)
754 {
755     const DCELL *x = xx;
756     const DCELL *y = yy;
757 
758     if (Rast_is_d_null_value(x))
759 	return 0;
760     if (*x < *y)
761 	return -1;
762     else if (*x == *y)
763 	return 0;
764     else
765 	return 1;
766 }
767 
768 /*!
769    \brief Returns quant rule which will be applied.
770 
771    Returns quant rule which will be applied when looking up the integer
772    quant value for val (used when organizing fp_lookup).
773 
774    \param q pointer to Quant structure which holds quant rules info
775    \param val fp cell value
776 
777    \return pointer to the Quant_table (color rule)
778    \return NULL otherwise
779  */
Rast__quant_get_rule_for_d_raster_val(const struct Quant * q,DCELL val)780 struct Quant_table *Rast__quant_get_rule_for_d_raster_val(const struct Quant
781 							  *q, DCELL val)
782 {
783     const struct Quant_table *p;
784 
785     for (p = &(q->table[q->nofRules - 1]); p >= q->table; p--)
786 	if ((val >= p->dLow) && (val <= p->dHigh))
787 	    break;
788     if (p >= q->table)
789 	return (struct Quant_table *)p;
790     else
791 	return (struct Quant_table *)NULL;
792 }
793