00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 #include <stdlib.h>
00286 #include <grass/gis.h>
00287
00288
00289 static int double_comp(const void *, const void *);
00290
00291 #define USE_LOOKUP 1
00292 #define MAX_LOOKUP_TABLE_SIZE 2048
00293 #define NO_DATA (G_set_c_null_value (&tmp, 1), (CELL) tmp)
00294
00295 #undef MIN
00296 #undef MAX
00297 #define MIN(a,b) ((a) < (b) ? (a) : (b))
00298 #define MAX(a,b) ((a) > (b) ? (a) : (b))
00299
00300 #define NO_LEFT_INFINITE_RULE (! q->infiniteLeftSet)
00301 #define NO_RIGHT_INFINITE_RULE (! q->infiniteRightSet)
00302 #define NO_FINITE_RULE (q->nofRules <= 0)
00303 #define NO_EXPLICIT_RULE (NO_FINITE_RULE && \
00304 NO_LEFT_INFINITE_RULE && NO_RIGHT_INFINITE_RULE)
00305
00306
00307
00308 void G_quant_clear(struct Quant *q)
00309 {
00310 q->nofRules = 0;
00311 q->infiniteRightSet = q->infiniteLeftSet = 0;
00312 }
00313
00314
00315
00316 void G_quant_free(struct Quant *q)
00317 {
00318 G_quant_clear(q);
00319
00320 if (q->maxNofRules > 0)
00321 G_free(q->table);
00322 if (q->fp_lookup.active) {
00323 G_free(q->fp_lookup.vals);
00324 G_free(q->fp_lookup.rules);
00325 q->fp_lookup.nalloc = 0;
00326 q->fp_lookup.active = 0;
00327 }
00328 q->maxNofRules = 0;
00329 }
00330
00331
00332
00333 int G__quant_organize_fp_lookup(struct Quant *q)
00334 {
00335 int i;
00336 DCELL val;
00337 CELL tmp;
00338 struct Quant_table *p;
00339
00340 if (q->nofRules * 2 > MAX_LOOKUP_TABLE_SIZE)
00341 return -1;
00342 if (q->nofRules == 0)
00343 return -1;
00344 q->fp_lookup.vals = (DCELL *)
00345 G_calloc(q->nofRules * 2, sizeof(DCELL));
00346
00347 q->fp_lookup.rules = (struct Quant_table **)
00348 G_calloc(q->nofRules * 2, sizeof(struct Quant_table *));
00349
00350
00351 if (!NO_FINITE_RULE) {
00352 i = 0;
00353
00354
00355
00356
00357
00358 for (p = &(q->table[q->nofRules - 1]); p >= q->table; p--) {
00359
00360 if (i == 0 || p->dLow != q->fp_lookup.vals[i - 1])
00361 q->fp_lookup.vals[i++] = p->dLow;
00362 q->fp_lookup.vals[i++] = p->dHigh;
00363 }
00364 q->fp_lookup.nalloc = i;
00365
00366
00367 qsort((char *)q->fp_lookup.vals, q->fp_lookup.nalloc,
00368 sizeof(DCELL), double_comp);
00369
00370
00371 for (i = 0; i < q->fp_lookup.nalloc - 1; i++) {
00372
00373
00374
00375 val = (q->fp_lookup.vals[i] + q->fp_lookup.vals[i + 1]) / 2.;
00376 q->fp_lookup.rules[i] =
00377 G__quant_get_rule_for_d_raster_val(q, val);
00378
00379
00380
00381
00382
00383
00384 }
00385 }
00386
00387 if (!NO_LEFT_INFINITE_RULE) {
00388 q->fp_lookup.inf_dmin = q->infiniteDLeft;
00389 q->fp_lookup.inf_min = q->infiniteCLeft;
00390 }
00391 else {
00392 if (q->fp_lookup.nalloc)
00393 q->fp_lookup.inf_dmin = q->fp_lookup.vals[0];
00394 q->fp_lookup.inf_min = NO_DATA;
00395 }
00396
00397 if (!NO_RIGHT_INFINITE_RULE) {
00398 if (q->fp_lookup.nalloc)
00399 q->fp_lookup.inf_dmax = q->infiniteDRight;
00400 q->fp_lookup.inf_max = q->infiniteCRight;
00401 }
00402 else {
00403 q->fp_lookup.inf_dmax = q->fp_lookup.vals[q->fp_lookup.nalloc - 1];
00404 q->fp_lookup.inf_max = NO_DATA;
00405 }
00406 q->fp_lookup.active = 1;
00407 return 1;
00408 }
00409
00410
00411
00412
00422 int G_quant_init(struct Quant *quant)
00423 {
00424 quant->fp_lookup.active = 0;
00425 quant->maxNofRules = 0;
00426 quant->truncate_only = 0;
00427 quant->round_only = 0;
00428 G_quant_clear(quant);
00429
00430 return 1;
00431 }
00432
00433
00434
00435 int G_quant_is_truncate(const struct Quant *quant)
00436 {
00437 return quant->truncate_only;
00438 }
00439
00440
00441
00442 int G_quant_is_round(const struct Quant *quant)
00443 {
00444 return quant->round_only;
00445 }
00446
00447
00448
00449
00471 int G_quant_truncate(struct Quant *quant)
00472 {
00473 quant->truncate_only = 1;
00474 return 1;
00475 }
00476
00477
00478
00479 int G_quant_round(struct Quant *quant)
00480 {
00481 quant->round_only = 1;
00482 return 1;
00483 }
00484
00485
00486
00487 static void quant_set_limits(struct Quant *q,
00488 DCELL dLow, DCELL dHigh, CELL cLow, CELL cHigh)
00489 {
00490 q->dMin = dLow;
00491 q->dMax = dHigh;
00492 q->cMin = cLow;
00493 q->cMax = cHigh;
00494 }
00495
00496
00497
00498 static void quant_update_limits(struct Quant *q,
00499 DCELL dLow, DCELL dHigh,
00500 CELL cLow, DCELL cHigh)
00501 {
00502 if (NO_EXPLICIT_RULE) {
00503 quant_set_limits(q, dLow, dHigh, cLow, cHigh);
00504 return;
00505 }
00506
00507 q->dMin = MIN(q->dMin, MIN(dLow, dHigh));
00508 q->dMax = MAX(q->dMax, MAX(dLow, dHigh));
00509 q->cMin = MIN(q->cMin, MIN(cLow, cHigh));
00510 q->cMax = MAX(q->cMax, MAX(cLow, cHigh));
00511 }
00512
00513
00514
00515
00534 int G_quant_get_limits(const struct Quant *q,
00535 DCELL * dMin, DCELL * dMax, CELL * cMin, CELL * cMax)
00536 {
00537 if (NO_EXPLICIT_RULE) {
00538 G_set_c_null_value(cMin, 1);
00539 G_set_c_null_value(cMax, 1);
00540 G_set_d_null_value(dMin, 1);
00541 G_set_d_null_value(dMax, 1);
00542 return -1;
00543 }
00544
00545 *dMin = q->dMin;
00546 *dMax = q->dMax;
00547 *cMin = q->cMin;
00548 *cMax = q->cMax;
00549
00550 return 1;
00551 }
00552
00553
00554
00555 int G_quant_nof_rules(const struct Quant *q)
00556 {
00557 return q->nofRules;
00558 }
00559
00560
00561
00562 void G_quant_get_ith_rule(const struct Quant *q,
00563 int i,
00564 DCELL * dLow, DCELL * dHigh,
00565 CELL * cLow, CELL * cHigh)
00566 {
00567 *dLow = q->table[i].dLow;
00568 *dHigh = q->table[i].dHigh;
00569 *cLow = q->table[i].cLow;
00570 *cHigh = q->table[i].cHigh;
00571 }
00572
00573
00574
00575 static void quant_table_increase(struct Quant *q)
00576 {
00577 if (q->nofRules < q->maxNofRules)
00578 return;
00579
00580 if (q->maxNofRules == 0) {
00581 q->maxNofRules = 50;
00582 q->table = (struct Quant_table *)
00583 G_malloc(q->maxNofRules * sizeof(struct Quant_table));
00584 }
00585 else {
00586 q->maxNofRules += 50;
00587 q->table = (struct Quant_table *)
00588 G_realloc((char *)q->table,
00589 q->maxNofRules * sizeof(struct Quant_table));
00590 }
00591 }
00592
00593
00594
00595 void G_quant_set_neg_infinite_rule(struct Quant *q, DCELL dLeft, CELL c)
00596 {
00597 q->infiniteDLeft = dLeft;
00598 q->infiniteCLeft = c;
00599 quant_update_limits(q, dLeft, dLeft, c, c);
00600
00601
00602 if (q->fp_lookup.active) {
00603 q->fp_lookup.inf_dmin = q->infiniteDLeft;
00604 q->fp_lookup.inf_min = q->infiniteCLeft;
00605 }
00606 q->infiniteLeftSet = 1;
00607 }
00608
00609
00610
00611 int G_quant_get_neg_infinite_rule(const struct Quant *q,
00612 DCELL * dLeft, CELL * c)
00613 {
00614 if (q->infiniteLeftSet == 0)
00615 return 0;
00616
00617 *dLeft = q->infiniteDLeft;
00618 *c = q->infiniteCLeft;
00619
00620 return 1;
00621 }
00622
00623
00624
00625 void G_quant_set_pos_infinite_rule(struct Quant *q, DCELL dRight, CELL c)
00626 {
00627 q->infiniteDRight = dRight;
00628 q->infiniteCRight = c;
00629 quant_update_limits(q, dRight, dRight, c, c);
00630
00631
00632 if (q->fp_lookup.active) {
00633 q->fp_lookup.inf_dmax = q->infiniteDRight;
00634 q->fp_lookup.inf_max = q->infiniteCRight;
00635 }
00636 q->infiniteRightSet = 1;
00637 }
00638
00639
00640
00641 int G_quant_get_pos_infinite_rule(const struct Quant *q,
00642 DCELL * dRight, CELL * c)
00643 {
00644 if (q->infiniteRightSet == 0)
00645 return 0;
00646
00647 *dRight = q->infiniteDRight;
00648 *c = q->infiniteCRight;
00649
00650 return 1;
00651 }
00652
00653
00654
00655 void G_quant_add_rule(struct Quant *q,
00656 DCELL dLow, DCELL dHigh, CELL cLow, CELL cHigh)
00657 {
00658 int i;
00659 struct Quant_table *p;
00660
00661 quant_table_increase(q);
00662
00663 i = q->nofRules;
00664
00665 p = &(q->table[i]);
00666 if (dHigh >= dLow) {
00667 p->dLow = dLow;
00668 p->dHigh = dHigh;
00669 p->cLow = cLow;
00670 p->cHigh = cHigh;
00671 }
00672 else {
00673 p->dLow = dHigh;
00674 p->dHigh = dLow;
00675 p->cLow = cHigh;
00676 p->cHigh = cLow;
00677 }
00678
00679
00680 if (q->fp_lookup.active) {
00681 G_free(q->fp_lookup.vals);
00682 G_free(q->fp_lookup.rules);
00683 q->fp_lookup.active = 0;
00684 q->fp_lookup.nalloc = 0;
00685 }
00686
00687 quant_update_limits(q, dLow, dHigh, cLow, cHigh);
00688
00689 q->nofRules++;
00690 }
00691
00692
00693
00694 void G_quant_reverse_rule_order(struct Quant *q)
00695 {
00696 struct Quant_table tmp;
00697 struct Quant_table *pLeft, *pRight;
00698
00699 pLeft = q->table;
00700 pRight = &(q->table[q->nofRules - 1]);
00701
00702 while (pLeft < pRight) {
00703 tmp.dLow = pLeft->dLow;
00704 tmp.dHigh = pLeft->dHigh;
00705 tmp.cLow = pLeft->cLow;
00706 tmp.cHigh = pLeft->cHigh;
00707
00708 pLeft->dLow = pRight->dLow;
00709 pLeft->dHigh = pRight->dHigh;
00710 pLeft->cLow = pRight->cLow;
00711 pLeft->cHigh = pRight->cHigh;
00712
00713 pRight->dLow = tmp.dLow;
00714 pRight->dHigh = tmp.dHigh;
00715 pRight->cLow = tmp.cLow;
00716 pRight->cHigh = tmp.cHigh;
00717
00718 pLeft++;
00719 pRight--;
00720 }
00721 }
00722
00723
00724
00725 static CELL quant_interpolate(DCELL dLow, DCELL dHigh,
00726 CELL cLow, CELL cHigh, DCELL dValue)
00727 {
00728 if (cLow == cHigh)
00729 return cLow;
00730 if (dLow == dHigh)
00731 return cLow;
00732
00733 return (CELL) ((dValue - dLow) / (dHigh - dLow) * (DCELL) (cHigh - cLow) +
00734 (DCELL) cLow);
00735 }
00736
00737
00738 static int less_or_equal(double x, double y)
00739 {
00740 if (x <= y)
00741 return 1;
00742 else
00743 return 0;
00744 }
00745
00746 static int less(double x, double y)
00747 {
00748 if (x < y)
00749 return 1;
00750 else
00751 return 0;
00752 }
00753
00754
00775 CELL G_quant_get_cell_value(struct Quant * q, DCELL dcellVal)
00776 {
00777 CELL tmp;
00778 DCELL dtmp;
00779 int try, min_ind, max_ind;
00780 struct Quant_table *p;
00781 int (*lower) ();
00782
00783 dtmp = dcellVal;
00784
00785
00786 if (G_is_d_null_value(&dtmp))
00787 return NO_DATA;
00788
00789 if (q->truncate_only)
00790 return (CELL) dtmp;
00791
00792 if (q->round_only) {
00793 if (dcellVal > 0)
00794 return (CELL) (dcellVal + .5);
00795 return (CELL) (dcellVal - .5);
00796 }
00797
00798 if (NO_EXPLICIT_RULE)
00799 return NO_DATA;
00800 if (NO_EXPLICIT_RULE)
00801 return NO_DATA;
00802
00803 if (USE_LOOKUP &&
00804 (q->fp_lookup.active || G__quant_organize_fp_lookup(q) > 0)) {
00805
00806
00807 if (dcellVal < q->fp_lookup.vals[0]) {
00808 if (dcellVal <= q->fp_lookup.inf_dmin)
00809 return q->fp_lookup.inf_min;
00810 else
00811 return NO_DATA;
00812 }
00813
00814 if (dcellVal > q->fp_lookup.vals[q->fp_lookup.nalloc - 1]) {
00815 if (dcellVal >= q->fp_lookup.inf_dmax)
00816 return q->fp_lookup.inf_max;
00817 else
00818 return NO_DATA;
00819 }
00820
00821
00822 try = (q->fp_lookup.nalloc - 1) / 2;
00823 min_ind = 0;
00824 max_ind = q->fp_lookup.nalloc - 2;
00825 while (1) {
00826
00827
00828
00829
00830
00831 if (q->fp_lookup.rules[try])
00832 lower = less;
00833 else
00834 lower = less_or_equal;
00835
00836 if (lower(q->fp_lookup.vals[try + 1], dcellVal)) {
00837 min_ind = try + 1;
00838
00839 try = (max_ind + min_ind) / 2;
00840 continue;
00841 }
00842 if (lower(dcellVal, q->fp_lookup.vals[try])) {
00843 max_ind = try - 1;
00844
00845 try = (max_ind + min_ind) / 2;
00846 continue;
00847 }
00848
00849 p = q->fp_lookup.rules[try];
00850 if (p)
00851 return quant_interpolate(p->dLow, p->dHigh, p->cLow, p->cHigh,
00852 dcellVal);
00853
00854 else {
00855 if (dcellVal <= q->fp_lookup.inf_dmin)
00856 return q->fp_lookup.inf_min;
00857 if (dcellVal >= q->fp_lookup.inf_dmax)
00858 return q->fp_lookup.inf_max;
00859 else
00860 return NO_DATA;
00861 }
00862 }
00863 }
00864
00865 if (!NO_FINITE_RULE) {
00866 p = G__quant_get_rule_for_d_raster_val(q, dcellVal);
00867 if (!p)
00868 return NO_DATA;
00869 return quant_interpolate(p->dLow, p->dHigh, p->cLow, p->cHigh,
00870 dcellVal);
00871 }
00872
00873 if ((!NO_LEFT_INFINITE_RULE) && (dcellVal <= q->infiniteDLeft))
00874 return q->infiniteCLeft;
00875
00876 if ((NO_RIGHT_INFINITE_RULE) || (dcellVal < q->infiniteDRight))
00877 return NO_DATA;
00878
00879 return q->infiniteCRight;
00880 }
00881
00882
00883
00884 void G_quant_perform_d(struct Quant *q,
00885 const DCELL * dcell, CELL * cell, int n)
00886 {
00887 int i;
00888
00889 for (i = 0; i < n; i++, dcell++)
00890 if (!G_is_d_null_value(dcell))
00891 *cell++ = G_quant_get_cell_value(q, *dcell);
00892 else
00893 G_set_c_null_value(cell++, 1);
00894 }
00895
00896
00897
00898 void G_quant_perform_f(struct Quant *q,
00899 const FCELL * fcell, CELL * cell, int n)
00900 {
00901 int i;
00902
00903 for (i = 0; i < n; i++, fcell++)
00904 if (!G_is_f_null_value(fcell))
00905 *cell++ = G_quant_get_cell_value(q, (DCELL) * fcell);
00906 else
00907 G_set_c_null_value(cell++, 1);
00908 }
00909
00910
00911
00912 static int double_comp(const void *xx, const void *yy)
00913 {
00914 const DCELL *x = xx;
00915 const DCELL *y = yy;
00916
00917 if (G_is_d_null_value(x))
00918 return 0;
00919 if (*x < *y)
00920 return -1;
00921 else if (*x == *y)
00922 return 0;
00923 else
00924 return 1;
00925 }
00926
00927
00928
00929 struct Quant_table *G__quant_get_rule_for_d_raster_val(const struct Quant *q,
00930 DCELL val)
00931 {
00932 const struct Quant_table *p;
00933
00934 for (p = &(q->table[q->nofRules - 1]); p >= q->table; p--)
00935 if ((val >= p->dLow) && (val <= p->dHigh))
00936 break;
00937 if (p >= q->table)
00938 return (struct Quant_table *)p;
00939 else
00940 return (struct Quant_table *)NULL;
00941 }
00942
00943
00944
00945
00946
00947