forked from tbriand/inverse_compositional
-
Notifications
You must be signed in to change notification settings - Fork 0
/
bicubic_interpolation.cpp
executable file
·182 lines (159 loc) · 4.92 KB
/
bicubic_interpolation.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
// This program is free software: you can use, modify and/or redistribute it
// under the terms of the simplified BSD License. You should have received a
// copy of this license along this program. If not, see
// <http://www.opensource.org/licenses/bsd-license.html>.
//
// Copyright (C) 2017, Thibaud Briand <[email protected]>
// Copyright (C) 2015, Javier Sánchez Pérez <[email protected]>
// Copyright (C) 2014, Nelson Monzón López <[email protected]>
// All rights reserved.
#include "bicubic_interpolation.h"
#include "transformation.h"
#include <cmath>
#include <stdio.h>
#include <stdlib.h>
#include "smapa.h"
SMART_PARAMETER(NANIFOUTSIDE,1) //set boundary pixels values to 0 if 0 and NAN if 1
SMART_PARAMETER(EDGEPADDING,5) //control the width of discarded values
/**
*
* Neumann boundary condition test
*
**/
int
neumann_bc (int x, int nx)
{
if (x<0)
x = 0;
else if (x >= nx)
x = nx - 1;
return x;
}
/**
*
* Bicubic interpolation in one dimension
*
**/
double
cubic_interpolation(
double v[4], //interpolation points
double x //point to be interpolated
)
{
return v[1] + 0.5 * x * (v[2] - v[0]
+ x * (2.0 * v[0] - 5.0 * v[1] + 4.0 * v[2] - v[3]
+ x * (3.0 * (v[1] - v[2]) + v[3] - v[0])));
}
/**
*
* Bicubic interpolation in two dimension
*
**/
double
bicubic_interpolation(
double p[4][4], //array containing the interpolation points
double x, //x position to be interpolated
double y //y position to be interpolated
)
{
double v[4];
v[0] = cubic_interpolation (p[0], y);
v[1] = cubic_interpolation (p[1], y);
v[2] = cubic_interpolation (p[2], y);
v[3] = cubic_interpolation (p[3], y);
return cubic_interpolation (v, x);
}
/**
*
* Compute the bicubic interpolation of a point in an image.
* Detects if the point goes outside the image domain
*
**/
double
bicubic_interpolation(
double *input,//image to be interpolated
double uu, //x component of the vector field
double vv, //y component of the vector field
int nx, //width of the image
int ny, //height of the image
int nz, //number of channels of the image
int k //actual channel
)
{
int sx = (uu < 0) ? -1 : 1;
int sy = (vv < 0) ? -1 : 1;
int x, y, mx, my, dx, dy, ddx, ddy;
x = neumann_bc ((int) uu, nx);
y = neumann_bc ((int) vv, ny);
mx = neumann_bc ((int) uu - sx, nx);
my = neumann_bc ((int) vv - sx, ny);
dx = neumann_bc ((int) uu + sx, nx);
dy = neumann_bc ((int) vv + sy, ny);
ddx = neumann_bc ((int) uu + 2 * sx, nx);
ddy = neumann_bc ((int) vv + 2 * sy, ny);
//obtain the interpolation points of the image
double p11 = input[(mx + nx * my) * nz + k];
double p12 = input[(x + nx * my) * nz + k];
double p13 = input[(dx + nx * my) * nz + k];
double p14 = input[(ddx + nx * my) * nz + k];
double p21 = input[(mx + nx * y) * nz + k];
double p22 = input[(x + nx * y) * nz + k];
double p23 = input[(dx + nx * y) * nz + k];
double p24 = input[(ddx + nx * y) * nz + k];
double p31 = input[(mx + nx * dy) * nz + k];
double p32 = input[(x + nx * dy) * nz + k];
double p33 = input[(dx + nx * dy) * nz + k];
double p34 = input[(ddx + nx * dy) * nz + k];
double p41 = input[(mx + nx * ddy) * nz + k];
double p42 = input[(x + nx * ddy) * nz + k];
double p43 = input[(dx + nx * ddy) * nz + k];
double p44 = input[(ddx + nx * ddy) * nz + k];
//create array
double pol[4][4] = {
{p11, p21, p31, p41}, {p12, p22, p32, p42},
{p13, p23, p33, p43}, {p14, p24, p34, p44}
};
//return interpolation
return bicubic_interpolation (pol, (double) uu - x, (double) vv - y);
}
/**
*
* Compute the bicubic interpolation of an image from a parametric transform
*
**/
void bicubic_interpolation(
double *input, //image to be warped
double *output, //warped output image with bicubic interpolation
double *params, //parameters of the transformation
int nparams, //number of parameters of the transform
int nx, //width of the image
int ny, //height of the image
int nz //number of channels of the image
)
{
double out_value = 0;
if ( NANIFOUTSIDE() )
out_value = NAN;
int delta = EDGEPADDING();
for (int i=0; i<ny; i++)
for (int j=0; j<nx; j++) {
int p=i*nx+j;
double x, y;
//transform coordinates using the parametric model
project(j, i, params, x, y, nparams);
bool out = false;
if ( x < delta || x > nx-1-delta
|| y < delta || y > ny-1-delta )
out = true;
if ( out ) {
for(int k=0; k<nz; k++)
output[p*nz+k]=out_value;
}
else {
//obtain the bicubic interpolation at position (uu, vv)
for(int k=0; k<nz; k++)
output[p*nz+k]=bicubic_interpolation(
input, x, y, nx, ny, nz, k);
}
}
}