quant.c 21 KB

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