← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.resamp.bspline/resamp.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.resamp.bspline/resamp.c
 
1
 
 
2
/***********************************************************************
 
3
 *
 
4
 * MODULE:       r.resamp.bspline
 
5
 *
 
6
 * AUTHOR(S):    Markus Metz
 
7
 *
 
8
 * PURPOSE:      Spline Interpolation and cross correlation
 
9
 *
 
10
 * COPYRIGHT:    (C) 2010 GRASS development team
 
11
 *
 
12
 *               This program is free software under the
 
13
 *               GNU General Public License (>=v2).
 
14
 *               Read the file COPYING that comes with GRASS
 
15
 *               for details.
 
16
 *
 
17
 **************************************************************************/
 
18
 
 
19
#include <grass/config.h>
 
20
#include <stdlib.h>
 
21
#include <string.h>
 
22
#include <math.h>
 
23
#include "bspline.h"
 
24
 
 
25
struct Point *P_Read_Raster_Region_masked(SEGMENT *mask_seg,
 
26
                                       struct Cell_head *Original,
 
27
                                       struct bound_box output_box,
 
28
                                       struct bound_box General,
 
29
                                       int *num_points, int dim_vect,
 
30
                                       double mean)
 
31
{
 
32
    int col, row, startcol, endcol, startrow, endrow, nrows, ncols;
 
33
    int pippo, npoints;
 
34
    double X, Y;
 
35
    struct Point *obs;
 
36
    char mask_val;
 
37
 
 
38
    pippo = dim_vect;
 
39
    obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
 
40
 
 
41
    /* Reading points inside output box and inside General box */
 
42
 
 
43
    npoints = 0;
 
44
    nrows = Original->rows;
 
45
    ncols = Original->cols;
 
46
 
 
47
    /* original region = output region
 
48
     * -> General box is somewhere inside output region */
 
49
    if (Original->north > General.N) {
 
50
        startrow = (double)((Original->north - General.N) / Original->ns_res - 1);
 
51
        if (startrow < 0)
 
52
            startrow = 0;
 
53
    }
 
54
    else
 
55
        startrow = 0;
 
56
    if (Original->south < General.S) {
 
57
        endrow = (double)((Original->north - General.S) / Original->ns_res + 1);
 
58
        if (endrow > nrows)
 
59
            endrow = nrows;
 
60
    }
 
61
    else
 
62
        endrow = nrows;
 
63
    if (General.W > Original->west) {
 
64
        startcol = (double)((General.W - Original->west) / Original->ew_res - 1);
 
65
        if (startcol < 0)
 
66
            startcol = 0;
 
67
    }
 
68
    else
 
69
        startcol = 0;
 
70
    if (General.E < Original->east) {
 
71
        endcol = (double)((General.E - Original->west) / Original->ew_res + 1);
 
72
        if (endcol > ncols)
 
73
            endcol = ncols;
 
74
    }
 
75
    else
 
76
        endcol = ncols;
 
77
 
 
78
    for (row = startrow; row < endrow; row++) {
 
79
        for (col = startcol; col < endcol; col++) {
 
80
 
 
81
            Segment_get(mask_seg, &mask_val, row, col);
 
82
            if (!mask_val)
 
83
                continue;
 
84
 
 
85
            X = Rast_col_to_easting((double)(col) + 0.5, Original);
 
86
            Y = Rast_row_to_northing((double)(row) + 0.5, Original);
 
87
 
 
88
            /* Here, mean is just for asking if obs point is in box */
 
89
            if (Vect_point_in_box(X, Y, mean, &General)) { /* General */
 
90
                if (npoints >= pippo) {
 
91
                    pippo += dim_vect;
 
92
                    obs =
 
93
                        (struct Point *)G_realloc((void *)obs,
 
94
                                                  (signed int)pippo *
 
95
                                                  sizeof(struct Point));
 
96
                }
 
97
 
 
98
                /* Storing observation vector */
 
99
                obs[npoints].coordX = X;
 
100
                obs[npoints].coordY = Y;
 
101
                obs[npoints].coordZ = 0;
 
102
                npoints++;
 
103
            }
 
104
        }
 
105
    }
 
106
 
 
107
    *num_points = npoints;
 
108
    return obs;
 
109
}
 
110
 
 
111
int P_Sparse_Raster_Points(SEGMENT *out_seg, struct Cell_head *Elaboration,
 
112
                struct Cell_head *Original, struct bound_box General, struct bound_box Overlap,
 
113
                struct Point *obs, double *param, double pe, double pn,
 
114
                double overlap, int nsplx, int nsply, int num_points,
 
115
                int bilin, double mean)
 
116
{
 
117
    int i, row, col;
 
118
    double X, Y, interpolation, csi, eta, weight, dval;
 
119
    int points_in_box = 0;
 
120
 
 
121
    /* Reading points inside output region and inside general box */
 
122
    /* all points available here are inside the output box,
 
123
     * selected by P_Read_Raster_Region_Nulls(), no check needed */
 
124
 
 
125
    for (i = 0; i < num_points; i++) {
 
126
 
 
127
        X = obs[i].coordX;
 
128
        Y = obs[i].coordY;
 
129
 
 
130
        /* X,Y are cell center cordinates, MUST be inside General box */
 
131
        row = (int) (floor(Rast_northing_to_row(Y, Original)) + 0.1);
 
132
        col = (int) (floor((X - Original->west) / Original->ew_res) + 0.1);
 
133
 
 
134
        if (row < 0 || row >= Original->rows) {
 
135
            G_fatal_error("row index out of range");
 
136
            continue;
 
137
        }
 
138
 
 
139
        if (col < 0 || col >= Original->cols) {
 
140
            G_fatal_error("col index out of range");
 
141
            continue;
 
142
        }
 
143
        points_in_box++;
 
144
 
 
145
        G_debug(3, "P_Sparse_Raster_Points: interpolate point %d...", i);
 
146
        if (bilin)
 
147
            interpolation =
 
148
                dataInterpolateBilin(X, Y, pe, pn, nsplx,
 
149
                                     nsply, Elaboration->west,
 
150
                                     Elaboration->south, param);
 
151
        else
 
152
            interpolation =
 
153
                dataInterpolateBicubic(X, Y, pe, pn, nsplx,
 
154
                                       nsply, Elaboration->west,
 
155
                                       Elaboration->south, param);
 
156
 
 
157
        interpolation += mean;
 
158
 
 
159
        if (Vect_point_in_box(X, Y, interpolation, &Overlap)) { /* (5) */
 
160
            dval = interpolation;
 
161
        }
 
162
        else {
 
163
            Segment_get(out_seg, &dval, row, col);
 
164
            if ((X > Overlap.E) && (X < General.E)) {
 
165
                if ((Y > Overlap.N) && (Y < General.N)) {       /* (3) */
 
166
                    csi = (General.E - X) / overlap;
 
167
                    eta = (General.N - Y) / overlap;
 
168
                    weight = csi * eta;
 
169
                    interpolation *= weight;
 
170
                    dval += interpolation;
 
171
                }
 
172
                else if ((Y < Overlap.S) && (Y > General.S)) {  /* (1) */
 
173
                    csi = (General.E - X) / overlap;
 
174
                    eta = (Y - General.S) / overlap;
 
175
                    weight = csi * eta;
 
176
                    interpolation *= weight;
 
177
                    dval = interpolation;
 
178
                }
 
179
                else if ((Y >= Overlap.S) && (Y <= Overlap.N)) {        /* (1) */
 
180
                    weight = (General.E - X ) / overlap;
 
181
                    interpolation *= weight;
 
182
                    dval = interpolation;
 
183
                }
 
184
            }
 
185
            else if ((X < Overlap.W) && (X > General.W)) {
 
186
                if ((Y > Overlap.N) && (Y < General.N)) {       /* (4) */
 
187
                    csi = (X - General.W) / overlap;
 
188
                    eta = (General.N - Y) / overlap;
 
189
                    weight = eta * csi;
 
190
                    interpolation *= weight;
 
191
                    dval += interpolation;
 
192
                }
 
193
                else if ((Y < Overlap.S) && (Y > General.S)) {  /* (2) */
 
194
                    csi = (X - General.W) / overlap;
 
195
                    eta = (Y - General.S) / overlap;
 
196
                    weight = csi * eta;
 
197
                    interpolation *= weight;
 
198
                    dval += interpolation;
 
199
                }
 
200
                else if ((Y >= Overlap.S) && (Y <= Overlap.N)) {        /* (2) */
 
201
                    weight = (X - General.W) / overlap;
 
202
                    interpolation *= weight;
 
203
                    dval += interpolation;
 
204
                }
 
205
            }
 
206
            else if ((X >= Overlap.W) && (X <= Overlap.E)) {
 
207
                if ((Y > Overlap.N) && (Y < General.N)) {       /* (3) */
 
208
                    weight = (General.N - Y) / overlap;
 
209
                    interpolation *= weight;
 
210
                    dval += interpolation;
 
211
                }
 
212
                else if ((Y < Overlap.S) && (Y > General.S)) {  /* (1) */
 
213
                    weight = (Y - General.S) / overlap;
 
214
                    interpolation *= weight;
 
215
                    dval = interpolation;
 
216
                }
 
217
            }
 
218
        } /* end not in overlap */
 
219
        Segment_put(out_seg, &dval, row, col);
 
220
    }  /* for num_points */
 
221
 
 
222
    return 1;
 
223
}
 
224
 
 
225
/* align elaboration box to source region
 
226
 * grow each side */
 
227
int align_elaboration_box(struct Cell_head *elaboration, struct Cell_head *original, int type)
 
228
{
 
229
    int row, col;
 
230
 
 
231
    switch (type) {
 
232
    case GENERAL_ROW:           /* General case N-S direction */
 
233
        /* northern edge */
 
234
        row = (int)((original->north - elaboration->north) / original->ns_res);
 
235
 
 
236
        if (row < 0)
 
237
            row = 0;
 
238
 
 
239
        elaboration->north = original->north - original->ns_res * row;
 
240
        
 
241
        /* southern edge */
 
242
        row = (int)((original->north - elaboration->south) / original->ns_res) + 1;
 
243
 
 
244
        if (row > original->rows + 1)
 
245
            row = original->rows + 1;
 
246
 
 
247
        elaboration->south = original->north - original->ns_res * row;
 
248
 
 
249
        return 1;
 
250
    
 
251
    case GENERAL_COLUMN:        /* General case E-W direction */
 
252
 
 
253
        /* eastern edge */
 
254
        col = (int)((elaboration->east - original->west) / original->ew_res) + 1;
 
255
 
 
256
        if (col > original->cols + 1)
 
257
            col = original->cols + 1;
 
258
 
 
259
        elaboration->east = original->west + original->ew_res * col;
 
260
        
 
261
        /* western edge */
 
262
        col = (int)((elaboration->west - original->west) / original->ew_res);
 
263
 
 
264
        if (col < 0)
 
265
            col = 0;
 
266
 
 
267
        elaboration->west = original->west + original->ew_res * col;
 
268
 
 
269
        return 1;
 
270
    }
 
271
    return 0;
 
272
}
 
273
 
 
274
 
 
275
/* align interpolation boxes to destination region
 
276
 * return 1 on success, 0 on failure */
 
277
 
 
278
int align_interp_boxes(struct bound_box *general, struct bound_box *overlap,
 
279
                       struct Cell_head *original, struct bound_box last_general,
 
280
                       struct bound_box last_overlap, int type)
 
281
{
 
282
    int row, col;
 
283
 
 
284
    switch (type) {
 
285
    case GENERAL_ROW:           /* General case N-S direction */
 
286
 
 
287
        /* general box */
 
288
        /* grow north */
 
289
        general->N = last_overlap.S;
 
290
 
 
291
        /* shrink south */
 
292
        row = (int)((original->north - general->S) / original->ns_res);
 
293
 
 
294
        if (row > original->rows + 1)
 
295
            row = original->rows + 1;
 
296
 
 
297
        general->S = original->north - original->ns_res * row;
 
298
        
 
299
        /* overlap box */
 
300
        /* grow north */
 
301
        overlap->N = last_general.S;
 
302
 
 
303
        /* shrink south */
 
304
        row = (int)((original->north - overlap->S) / original->ns_res);
 
305
 
 
306
        if (row > original->rows + 1)
 
307
            row = original->rows + 1;
 
308
 
 
309
        overlap->S = original->north - original->ns_res * row;
 
310
 
 
311
        return 1;
 
312
 
 
313
    case GENERAL_COLUMN:        /* General case E-W direction */
 
314
 
 
315
        /* general box */
 
316
        /* grow west */
 
317
        general->W = last_overlap.E;
 
318
 
 
319
        /* shrink east */
 
320
        col = (int)((general->E - original->west) / original->ew_res);
 
321
 
 
322
        if (col > original->cols + 1)
 
323
            col = original->cols + 1;
 
324
 
 
325
        general->E = original->west + original->ew_res * col;
 
326
        
 
327
        /* overlap box */
 
328
        /* grow west */
 
329
        overlap->W = last_general.E;
 
330
 
 
331
        /* shrink east */
 
332
        col = (int)((overlap->E - original->west) / original->ew_res);
 
333
 
 
334
        if (col > original->cols + 1)
 
335
            col = original->cols + 1;
 
336
 
 
337
        overlap->E = original->west + original->ew_res * col;
 
338
 
 
339
        return 1;
 
340
 
 
341
    case FIRST_ROW:             /* Just started with first row */
 
342
        general->N = original->north;
 
343
        overlap->N = original->north;
 
344
 
 
345
        /* shrink south */
 
346
        row = (int)((original->north - general->S) / original->ns_res);
 
347
 
 
348
        if (row > original->rows)
 
349
            row = original->rows;
 
350
 
 
351
        general->S = original->north - original->ns_res * row;
 
352
 
 
353
        row = (int)((original->north - overlap->S) / original->ns_res);
 
354
 
 
355
        if (row > original->rows + 1)
 
356
            row = original->rows + 1;
 
357
 
 
358
        overlap->S = original->north - original->ns_res * row;
 
359
 
 
360
        return 1;
 
361
 
 
362
    case LAST_ROW:              /* Reached last row */
 
363
        general->S = original->south;
 
364
        overlap->S = original->south;
 
365
 
 
366
        return 1;
 
367
 
 
368
    case FIRST_COLUMN:          /* Just started with first column */
 
369
        general->W = original->west;
 
370
        overlap->W = original->west;
 
371
 
 
372
        /* shrink east */
 
373
        col = (int)((general->E - original->west) / original->ew_res);
 
374
 
 
375
        if (col > original->cols + 1)
 
376
            col = original->cols + 1;
 
377
 
 
378
        general->E = original->west + original->ew_res * col;
 
379
 
 
380
        col = (int)((overlap->E - original->west) / original->ew_res);
 
381
 
 
382
        if (col > original->cols + 1)
 
383
            col = original->cols + 1;
 
384
 
 
385
        overlap->E = original->west + original->ew_res * col;
 
386
 
 
387
        return 1;
 
388
 
 
389
    case LAST_COLUMN:           /* Reached last column */
 
390
        general->E = original->east;
 
391
        overlap->E = original->east;
 
392
 
 
393
        return 1;
 
394
    }
 
395
 
 
396
    return 0;
 
397
}