1 /*
2 * Name
3 * cubic.c -- use cubic convolution interpolation for given row, col
4 *
5 * Description
6 * cubic returns the value in the buffer that is the result of cubic
7 * convolution 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 <math.h>
16 #include "global.h"
17
p_cubic(struct cache * ibuffer,void * obufptr,int cell_type,double * row_idx,double * col_idx,struct Cell_head * cellhd)18 void p_cubic(struct cache *ibuffer, /* input buffer */
19 void *obufptr, /* ptr in output buffer */
20 int cell_type, /* raster map type of obufptr */
21 double *row_idx, /* row index (decimal) */
22 double *col_idx, /* column index (decimal) */
23 struct Cell_head *cellhd /* information of output map */
24 )
25 {
26 int row; /* row indices for interp */
27 int col; /* column indices for interp */
28 int i, j;
29 double t, u; /* intermediate slope */
30 DCELL result; /* result of interpolation */
31 DCELL val[4]; /* buffer for temporary values */
32 DCELL cell[4][4];
33
34 /* cut indices to integer */
35 row = (int)floor(*row_idx - 0.5);
36 col = (int)floor(*col_idx - 0.5);
37
38 /* check for out of bounds of map - if out of bounds set NULL value */
39 if (row - 1 < 0 || row + 2 >= cellhd->rows ||
40 col - 1 < 0 || col + 2 >= cellhd->cols) {
41 Rast_set_null_value(obufptr, 1, cell_type);
42 return;
43 }
44
45 for (i = 0; i < 4; i++)
46 for (j = 0; j < 4; j++) {
47 const DCELL *cellp = CPTR(ibuffer, row - 1 + i, col - 1 + j);
48 if (Rast_is_d_null_value(cellp)) {
49 Rast_set_null_value(obufptr, 1, cell_type);
50 return;
51 }
52 cell[i][j] = *cellp;
53 }
54
55 /* do the interpolation */
56 t = *col_idx - 0.5 - col;
57 u = *row_idx - 0.5 - row;
58
59 for (i = 0; i < 4; i++) {
60 val[i] = Rast_interp_cubic(t, cell[i][0], cell[i][1], cell[i][2], cell[i][3]);
61 }
62
63 result = Rast_interp_cubic(u, val[0], val[1], val[2], val[3]);
64
65 Rast_set_d_value(obufptr, result, cell_type);
66 }
67