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 MIN(a,b) (a < b ? a : b)
00294 #define MAX(a,b) (a > b ? a : b)
00295 #define NO_DATA (G_set_c_null_value (&tmp, 1), (CELL) tmp)
00296
00297 #define NO_LEFT_INFINITE_RULE (! q->infiniteLeftSet)
00298 #define NO_RIGHT_INFINITE_RULE (! q->infiniteRightSet)
00299 #define NO_FINITE_RULE (q->nofRules <= 0)
00300 #define NO_EXPLICIT_RULE (NO_FINITE_RULE && \
00301 NO_LEFT_INFINITE_RULE && NO_RIGHT_INFINITE_RULE)
00302
00303
00304
00305 void G_quant_clear (struct Quant *q)
00306
00307 {
00308 q->nofRules = 0;
00309 q->infiniteRightSet = q->infiniteLeftSet = 0;
00310 }
00311
00312
00313
00314 void G_quant_free (struct Quant *q)
00315 {
00316 G_quant_clear (q);
00317
00318 if (q->maxNofRules > 0) G_free (q->table);
00319 if (q->fp_lookup.active)
00320 {
00321 G_free (q->fp_lookup.vals);
00322 G_free (q->fp_lookup.rules);
00323 q->fp_lookup.nalloc = 0;
00324 q->fp_lookup.active = 0;
00325 }
00326 q->maxNofRules = 0;
00327 }
00328
00329
00330
00331 int G__quant_organize_fp_lookup (struct Quant *q)
00332 {
00333 int i;
00334 DCELL val;
00335 CELL tmp;
00336 struct Quant_table *p;
00337
00338 if(q->nofRules * 2 > MAX_LOOKUP_TABLE_SIZE) return -1;
00339 if(q->nofRules == 0) return -1;
00340 q->fp_lookup.vals = (DCELL *)
00341 G_calloc(q->nofRules * 2, sizeof(DCELL));
00342
00343 q->fp_lookup.rules = (struct Quant_table **)
00344 G_calloc(q->nofRules * 2, sizeof(struct Quant_table *));
00345
00346
00347 if (! NO_FINITE_RULE)
00348 {
00349 i = 0;
00350
00351
00352
00353
00354
00355 for (p = &(q->table [q->nofRules - 1]); p >= q->table; p--)
00356 {
00357
00358 if(i==0 || p->dLow != q->fp_lookup.vals[i-1])
00359 q->fp_lookup.vals[i++] = p->dLow;
00360 q->fp_lookup.vals[i++] = p->dHigh;
00361 }
00362 q->fp_lookup.nalloc = i;
00363
00364
00365 qsort((char *) q->fp_lookup.vals, q->fp_lookup.nalloc,
00366 sizeof(DCELL), double_comp);
00367
00368
00369 for(i=0;i<q->fp_lookup.nalloc-1;i++)
00370 {
00371
00372
00373
00374 val = (q->fp_lookup.vals[i] + q->fp_lookup.vals[i+1])/2.;
00375 q->fp_lookup.rules[i] = G__quant_get_rule_for_d_raster_val (q, val);
00376
00377
00378
00379
00380
00381
00382 }
00383 }
00384
00385 if (! NO_LEFT_INFINITE_RULE)
00386 {
00387 q->fp_lookup.inf_dmin = q->infiniteDLeft;
00388 q->fp_lookup.inf_min = q->infiniteCLeft;
00389 }
00390 else
00391 {
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 {
00399 if(q->fp_lookup.nalloc)
00400 q->fp_lookup.inf_dmax = q->infiniteDRight;
00401 q->fp_lookup.inf_max = q->infiniteCRight;
00402 }
00403 else
00404 {
00405 q->fp_lookup.inf_dmax = q->fp_lookup.vals[q->fp_lookup.nalloc-1];
00406 q->fp_lookup.inf_max = NO_DATA;
00407 }
00408 q->fp_lookup.active = 1;
00409 return 1;
00410 }
00411
00412
00413
00414
00424 int G_quant_init (struct Quant *quant)
00425 {
00426 quant->fp_lookup.active = 0;
00427 quant->maxNofRules = 0;
00428 quant->truncate_only = 0;
00429 quant->round_only = 0;
00430 G_quant_clear (quant);
00431
00432 return 1;
00433 }
00434
00435
00436
00437 int G_quant_is_truncate (struct Quant *quant)
00438 {
00439 return quant->truncate_only;
00440 }
00441
00442
00443
00444 int G_quant_is_round (struct Quant *quant)
00445 {
00446 return quant->round_only;
00447 }
00448
00449
00450
00451
00473 int G_quant_truncate (struct Quant *quant)
00474 {
00475 quant->truncate_only = 1;
00476 return 1;
00477 }
00478
00479
00480
00481 int G_quant_round (struct Quant *quant)
00482 {
00483 quant->round_only = 1;
00484 return 1;
00485 }
00486
00487
00488
00489 static void quant_set_limits (
00490 struct Quant *q,
00491 DCELL dLow,DCELL dHigh,
00492 CELL cLow,CELL cHigh)
00493 {
00494 q->dMin = dLow;
00495 q->dMax = dHigh;
00496 q->cMin = cLow;
00497 q->cMax = cHigh;
00498 }
00499
00500
00501
00502 static void quant_update_limits (
00503 struct Quant *q,
00504 DCELL dLow,DCELL dHigh,
00505 CELL cLow,DCELL cHigh)
00506
00507 {
00508 if (NO_EXPLICIT_RULE) {
00509 quant_set_limits (q, dLow, dHigh, cLow, cHigh);
00510 return;
00511 }
00512
00513 q->dMin = MIN (q->dMin, MIN (dLow, dHigh));
00514 q->dMax = MAX (q->dMax, MAX (dLow, dHigh));
00515 q->cMin = MIN (q->cMin, MIN (cLow, cHigh));
00516 q->cMax = MAX (q->cMax, MAX (cLow, cHigh));
00517 }
00518
00519
00520
00521
00540 int G_quant_get_limits (
00541 struct Quant *q,
00542 DCELL *dMin,DCELL *dMax,
00543 CELL *cMin,CELL *cMax)
00544
00545 {
00546 if (NO_EXPLICIT_RULE)
00547 {
00548 G_set_c_null_value(cMin, 1);
00549 G_set_c_null_value(cMax, 1);
00550 G_set_d_null_value(dMin, 1);
00551 G_set_d_null_value(dMax, 1);
00552 return -1;
00553 }
00554
00555 *dMin = q->dMin; *dMax = q->dMax;
00556 *cMin = q->cMin; *cMax = q->cMax;
00557
00558 return 1;
00559 }
00560
00561
00562
00563 int G_quant_nof_rules (struct Quant *q)
00564
00565 {
00566 return q->nofRules;
00567 }
00568
00569
00570
00571 void G_quant_get_ith_rule (
00572 struct Quant *q,
00573 int i,
00574 DCELL *dLow,DCELL *dHigh,
00575 CELL *cLow,CELL *cHigh)
00576
00577 {
00578 *dLow = q->table[i].dLow;
00579 *dHigh = q->table[i].dHigh;
00580 *cLow = q->table[i].cLow;
00581 *cHigh = q->table[i].cHigh;
00582 }
00583
00584
00585
00586 static void quant_table_increase (struct Quant *q)
00587
00588 {
00589 if (q->nofRules < q->maxNofRules) return;
00590
00591 if (q->maxNofRules == 0) {
00592 q->maxNofRules = 50;
00593 q->table = (struct Quant_table *)
00594 G_malloc (q->maxNofRules * sizeof (struct Quant_table));
00595 } else {
00596 q->maxNofRules += 50;
00597 q->table = (struct Quant_table *)
00598 G_realloc ((char *) q->table,
00599 q->maxNofRules * sizeof (struct Quant_table));
00600 }
00601 }
00602
00603
00604
00605 void G_quant_set_neg_infinite_rule (
00606
00607 struct Quant *q,
00608 DCELL dLeft,
00609 CELL c)
00610
00611 {
00612 q->infiniteDLeft = dLeft;
00613 q->infiniteCLeft = c;
00614 quant_update_limits (q, dLeft, dLeft, c, c);
00615
00616
00617 if (q->fp_lookup.active)
00618 {
00619 q->fp_lookup.inf_dmin = q->infiniteDLeft;
00620 q->fp_lookup.inf_min = q->infiniteCLeft;
00621 }
00622 q->infiniteLeftSet = 1;
00623 }
00624
00625
00626
00627 int G_quant_get_neg_infinite_rule (
00628
00629 struct Quant *q,
00630 DCELL *dLeft,
00631 CELL *c)
00632
00633 {
00634 if (q->infiniteLeftSet == 0) return 0;
00635
00636 *dLeft = q->infiniteDLeft;
00637 *c = q->infiniteCLeft;
00638
00639 return 1;
00640 }
00641
00642
00643
00644 void G_quant_set_pos_infinite_rule (
00645
00646 struct Quant *q,
00647 DCELL dRight,
00648 CELL c)
00649
00650 {
00651 q->infiniteDRight = dRight;
00652 q->infiniteCRight = c;
00653 quant_update_limits (q, dRight, dRight, c, c);
00654
00655
00656 if (q->fp_lookup.active)
00657 {
00658 q->fp_lookup.inf_dmax = q->infiniteDRight;
00659 q->fp_lookup.inf_max = q->infiniteCRight;
00660 }
00661 q->infiniteRightSet = 1;
00662 }
00663
00664
00665
00666 int G_quant_get_pos_infinite_rule (
00667
00668 struct Quant *q,
00669 DCELL *dRight,
00670 CELL *c)
00671
00672 {
00673 if (q->infiniteRightSet == 0) return 0;
00674
00675 *dRight = q->infiniteDRight;
00676 *c = q->infiniteCRight;
00677
00678 return 1;
00679 }
00680
00681
00682
00683 void G_quant_add_rule (
00684
00685 struct Quant *q,
00686 DCELL dLow,DCELL dHigh,
00687 CELL cLow,CELL cHigh)
00688
00689 {
00690 int i;
00691 struct Quant_table *p;
00692
00693 quant_table_increase (q);
00694
00695 i = q->nofRules;
00696
00697 p = &(q->table[i]);
00698 if (dHigh >= dLow) {
00699 p->dLow = dLow; p->dHigh = dHigh; p->cLow = cLow; p->cHigh = cHigh;
00700 } else {
00701 p->dLow = dHigh; p->dHigh = dLow; p->cLow = cHigh; p->cHigh = cLow;
00702 }
00703
00704
00705 if (q->fp_lookup.active)
00706 {
00707 G_free (q->fp_lookup.vals);
00708 G_free (q->fp_lookup.rules);
00709 q->fp_lookup.active = 0;
00710 q->fp_lookup.nalloc = 0;
00711 }
00712
00713 quant_update_limits (q, dLow, dHigh, cLow, cHigh);
00714
00715 q->nofRules++;
00716 }
00717
00718
00719
00720 void G_quant_reverse_rule_order (struct Quant *q)
00721
00722 {
00723 struct Quant_table tmp;
00724 struct Quant_table *pLeft, *pRight;
00725
00726 pLeft = q->table;
00727 pRight = &(q->table [q->nofRules - 1]);
00728
00729 while (pLeft < pRight) {
00730 tmp.dLow = pLeft->dLow; tmp.dHigh = pLeft->dHigh;
00731 tmp.cLow = pLeft->cLow; tmp.cHigh = pLeft->cHigh;
00732
00733 pLeft->dLow = pRight->dLow; pLeft->dHigh = pRight->dHigh;
00734 pLeft->cLow = pRight->cLow; pLeft->cHigh = pRight->cHigh;
00735
00736 pRight->dLow = tmp.dLow; pRight->dHigh = tmp.dHigh;
00737 pRight->cLow = tmp.cLow; pRight->cHigh = tmp.cHigh;
00738
00739 pLeft++; pRight--;
00740 }
00741 }
00742
00743
00744
00745 static CELL quant_interpolate (
00746 DCELL dLow,DCELL dHigh,
00747 CELL cLow,CELL cHigh,
00748 DCELL dValue)
00749
00750 {
00751 if (cLow == cHigh) return cLow;
00752 if (dLow == dHigh) return cLow;
00753
00754 return (CELL) ((dValue - dLow) / (dHigh - dLow) * (DCELL) (cHigh - cLow) +
00755 (DCELL) cLow);
00756 }
00757
00758
00759 static int less_or_equal(double x,double y)
00760 {
00761 if(x<=y) return 1;
00762 else return 0;
00763 }
00764
00765 static int less( double x,double y)
00766 {
00767 if(x<y) return 1;
00768 else return 0;
00769 }
00770
00771
00792 CELL G_quant_get_cell_value ( struct Quant *q, DCELL dcellVal)
00793 {
00794 CELL tmp;
00795 DCELL dtmp;
00796 int try, min_ind, max_ind;
00797 struct Quant_table *p;
00798 int (*lower)();
00799
00800 dtmp = dcellVal;
00801
00802
00803 if(G_is_d_null_value(&dtmp))
00804 return NO_DATA;
00805
00806 if(q->truncate_only)
00807 return (CELL) dtmp;
00808
00809 if(q->round_only)
00810 {
00811 if(dcellVal > 0) return (CELL) (dcellVal + .5);
00812 return (CELL) (dcellVal - .5);
00813 }
00814
00815 if (NO_EXPLICIT_RULE) return NO_DATA;
00816 if (NO_EXPLICIT_RULE) return NO_DATA;
00817
00818 if(USE_LOOKUP && (q->fp_lookup.active || G__quant_organize_fp_lookup(q) > 0))
00819 {
00820
00821
00822 if(dcellVal < q->fp_lookup.vals[0])
00823 {
00824 if(dcellVal <= q->fp_lookup.inf_dmin)
00825 return q->fp_lookup.inf_min;
00826 else return NO_DATA;
00827 }
00828
00829 if(dcellVal > q->fp_lookup.vals[q->fp_lookup.nalloc-1])
00830 {
00831 if(dcellVal >= q->fp_lookup.inf_dmax)
00832 return q->fp_lookup.inf_max;
00833 else return NO_DATA;
00834 }
00835
00836
00837 try = (q->fp_lookup.nalloc-1)/2;
00838 min_ind = 0;
00839 max_ind = q->fp_lookup.nalloc-2;
00840 while(1)
00841 {
00842
00843
00844
00845
00846
00847 if(q->fp_lookup.rules[try])
00848 lower = less;
00849 else
00850 lower = less_or_equal;
00851
00852 if (lower(q->fp_lookup.vals[try+1], dcellVal))
00853 {
00854 min_ind = try+1;
00855
00856 try = (max_ind + min_ind)/2;
00857 continue;
00858 }
00859 if (lower(dcellVal, q->fp_lookup.vals[try]))
00860 {
00861 max_ind = try-1;
00862
00863 try = (max_ind + min_ind)/2;
00864 continue;
00865 }
00866
00867 p = q->fp_lookup.rules[try];
00868 if(p)
00869 return quant_interpolate (p->dLow, p->dHigh, p->cLow, p->cHigh,
00870 dcellVal);
00871
00872 else
00873 {
00874 if(dcellVal <= q->fp_lookup.inf_dmin)
00875 return q->fp_lookup.inf_min;
00876 if(dcellVal >= q->fp_lookup.inf_dmax)
00877 return q->fp_lookup.inf_max;
00878 else return NO_DATA;
00879 }
00880 }
00881 }
00882
00883 if (! NO_FINITE_RULE)
00884 {
00885 p = G__quant_get_rule_for_d_raster_val (q, dcellVal);
00886 if (!p)
00887 return NO_DATA;
00888 return quant_interpolate (p->dLow, p->dHigh, p->cLow, p->cHigh,
00889 dcellVal);
00890 }
00891
00892 if ((! NO_LEFT_INFINITE_RULE) && (dcellVal <= q->infiniteDLeft))
00893 return q->infiniteCLeft;
00894
00895 if ((NO_RIGHT_INFINITE_RULE) || (dcellVal < q->infiniteDRight))
00896 return NO_DATA;
00897
00898 return q->infiniteCRight;
00899 }
00900
00901
00902
00903 void G_quant_perform_d (
00904 struct Quant *q,
00905 DCELL *dcell,
00906 CELL *cell,
00907 int n)
00908
00909 {
00910 int i;
00911
00912 for (i = 0; i < n; i++, dcell++)
00913 if (! G_is_d_null_value (dcell))
00914 *cell++ = G_quant_get_cell_value (q, *dcell);
00915 else
00916 G_set_c_null_value (cell++, 1);
00917 }
00918
00919
00920
00921 void G_quant_perform_f (
00922 struct Quant *q,
00923 FCELL *fcell,
00924 CELL *cell,
00925 int n)
00926
00927 {
00928 int i;
00929
00930 for (i = 0; i < n; i++, fcell++)
00931 if (! G_is_f_null_value (fcell))
00932 *cell++ = G_quant_get_cell_value (q, (DCELL) *fcell);
00933 else
00934 G_set_c_null_value (cell++, 1);
00935 }
00936
00937
00938
00939 static int double_comp(const void *xx,const void *yy)
00940 {
00941 const DCELL *x = xx;
00942 const DCELL *y = yy;
00943 if(G_is_d_null_value(x)) return 0;
00944 if(*x < *y) return -1;
00945 else if(*x==*y) return 0;
00946 else return 1;
00947 }
00948
00949
00950
00951 struct Quant_table *G__quant_get_rule_for_d_raster_val (
00952 struct Quant *q,
00953 DCELL val)
00954 {
00955 struct Quant_table *p;
00956 for (p = &(q->table[q->nofRules - 1]); p >= q->table; p--)
00957 if ((val >= p->dLow) && (val <= p->dHigh))
00958 break;
00959 if(p >= q->table) return p;
00960 else return (struct Quant_table *) NULL;
00961 }
00962
00963
00964
00965