1 /*
2 * Name
3 * lanczos.c -- use lanczos interpolation for given row, col
4 *
5 * Description
6 * lanczos returns the value in the buffer that is the result of
7 * lanczos interpolation for the given row, column indices.
8 * If the given row or column is outside the bounds of the input map,
9 * the corresponding point in the output map is set to NULL.
10 *
11 * If single surrounding points in the interpolation matrix are
12 * NULL they where filled with their neighbor
13 */
14
15 #include <grass/gis.h>
16 #include <grass/raster.h>
17 #include <math.h>
18 #include "global.h"
19
p_lanczos(struct cache * ibuffer,void * obufptr,int cell_type,double * row_idx,double * col_idx,struct Cell_head * cellhd)20 void p_lanczos(struct cache *ibuffer, /* input buffer */
21 void *obufptr, /* ptr in output buffer */
22 int cell_type, /* raster map type of obufptr */
23 double *row_idx, /* row index (decimal) */
24 double *col_idx, /* column index (decimal) */
25 struct Cell_head *cellhd /* information of output map */
26 )
27 {
28 int row; /* row indices for interp */
29 int col; /* column indices for interp */
30 int i, j, k;
31 DCELL t, u; /* intermediate slope */
32 DCELL result; /* result of interpolation */
33 DCELL cell[25];
34
35 /* cut indices to integer */
36 row = (int)floor(*row_idx);
37 col = (int)floor(*col_idx);
38
39 /* check for out of bounds of map - if out of bounds set NULL value */
40 if (row - 2 < 0 || row + 2 >= cellhd->rows ||
41 col - 2 < 0 || col + 2 >= cellhd->cols) {
42 Rast_set_null_value(obufptr, 1, cell_type);
43 return;
44 }
45
46 k = 0;
47 for (i = 0; i < 5; i++) {
48 for (j = 0; j < 5; j++) {
49 const DCELL *cellp = CPTR(ibuffer, row - 2 + i, col - 2 + j);
50 if (Rast_is_d_null_value(cellp)) {
51 Rast_set_null_value(obufptr, 1, cell_type);
52 return;
53 }
54 cell[k++] = *cellp;
55 }
56 }
57
58 /* do the interpolation */
59 t = *col_idx - 0.5 - col;
60 u = *row_idx - 0.5 - row;
61
62 result = Rast_interp_lanczos(t, u, cell);
63
64 Rast_set_d_value(obufptr, result, cell_type);
65 }
66
p_lanczos_f(struct cache * ibuffer,void * obufptr,int cell_type,double * row_idx,double * col_idx,struct Cell_head * cellhd)67 void p_lanczos_f(struct cache *ibuffer, /* input buffer */
68 void *obufptr, /* ptr in output buffer */
69 int cell_type, /* raster map type of obufptr */
70 double *row_idx, /* row index (decimal) */
71 double *col_idx, /* column index (decimal) */
72 struct Cell_head *cellhd /* information of output map */
73 )
74 {
75 int row; /* row indices for interp */
76 int col; /* column indices for interp */
77 DCELL *cellp, cell;
78
79 /* cut indices to integer */
80 row = (int)floor(*row_idx);
81 col = (int)floor(*col_idx);
82
83 /* check for out of bounds - if out of bounds set NULL value */
84 if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
85 Rast_set_null_value(obufptr, 1, cell_type);
86 return;
87 }
88
89 cellp = CPTR(ibuffer, row, col);
90 /* if nearest is null, all the other interps will be null */
91 if (Rast_is_d_null_value(cellp)) {
92 Rast_set_null_value(obufptr, 1, cell_type);
93 return;
94 }
95 cell = *cellp;
96
97 p_lanczos(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
98 /* fallback to bicubic if lanczos is null */
99 if (Rast_is_d_null_value(obufptr)) {
100 p_cubic(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
101 /* fallback to bilinear if cubic is null */
102 if (Rast_is_d_null_value(obufptr)) {
103 p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
104 /* fallback to nearest if bilinear is null */
105 if (Rast_is_d_null_value(obufptr))
106 Rast_set_d_value(obufptr, cell, cell_type);
107 }
108 }
109 }
110