4 * Floating point image transformations
6 * Copyright (C) Georg Martius - June 2011
7 * georg dot martius at web dot de
9 * This file is part of vid.stab video stabilization library
11 * vid.stab is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License,
13 * as published by the Free Software Foundation; either version 2, or
14 * (at your option) any later version.
16 * vid.stab is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
21 * You should have received a copy of the GNU General Public License
22 * along with GNU Make; see the file COPYING. If not, write to
23 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
26 #include "transformfloat.h"
27 #include "transform.h"
28 #include "transformtype_operations.h"
31 /** interpolateBiLinBorder: bi-linear interpolation function that also works at the border.
32 This is used by many other interpolation methods at and outsize the border, see interpolate */
33 void _FLT(interpolateBiLinBorder
)(uint8_t *rv
, float x
, float y
,
34 const uint8_t *img
, int img_linesize
,
35 int width
, int height
, uint8_t def
)
41 short v1
= PIXEL(img
, img_linesize
, x_c
, y_c
, width
, height
, def
);
42 short v2
= PIXEL(img
, img_linesize
, x_c
, y_f
, width
, height
, def
);
43 short v3
= PIXEL(img
, img_linesize
, x_f
, y_c
, width
, height
, def
);
44 short v4
= PIXEL(img
, img_linesize
, x_f
, y_f
, width
, height
, def
);
45 float s
= (v1
*(x
- x_f
)+v3
*(x_c
- x
))*(y
- y_f
) +
46 (v2
*(x
- x_f
) + v4
*(x_c
- x
))*(y_c
- y
);
50 /** taken from http://en.wikipedia.org/wiki/Bicubic_interpolation for alpha=-0.5
52 a0-a3 are the neigthboring points where the target point is between a1 and a2
53 t is the point of interpolation (position between a1 and a2) value between 0 and 1
56 (1,t,t^2,t^3) | 2,-5, 4,-1 | |a2|
59 static short _FLT(bicub_kernel
)(float t
, short a0
, short a1
, short a2
, short a3
){
60 return (2*a1
+ t
*((-a0
+a2
) + t
*((2*a0
-5*a1
+4*a2
-a3
) + t
*(-a0
+3*a1
-3*a2
+a3
) )) ) / 2;
63 /** interpolateBiCub: bi-cubic interpolation function using 4x4 pixel, see interpolate */
64 void _FLT(interpolateBiCub
)(uint8_t *rv
, float x
, float y
,
65 const uint8_t *img
, int img_linesize
,
66 int width
, int height
, uint8_t def
)
68 // do a simple linear interpolation at the border
69 if (x
< 1 || x
> width
- 2 || y
< 1 || y
> height
- 2) {
70 _FLT(interpolateBiLinBorder
)(rv
, x
, y
, img
, img_linesize
, width
, height
, def
);
75 short v1
= _FLT(bicub_kernel
)(tx
,
76 PIX(img
, img_linesize
, x_f
-1, y_f
-1),
77 PIX(img
, img_linesize
, x_f
, y_f
-1),
78 PIX(img
, img_linesize
, x_f
+1, y_f
-1),
79 PIX(img
, img_linesize
, x_f
+2, y_f
-1));
80 short v2
= _FLT(bicub_kernel
)(tx
,
81 PIX(img
, img_linesize
, x_f
-1, y_f
),
82 PIX(img
, img_linesize
, x_f
, y_f
),
83 PIX(img
, img_linesize
, x_f
+1, y_f
),
84 PIX(img
, img_linesize
, x_f
+2, y_f
));
85 short v3
= _FLT(bicub_kernel
)(tx
,
86 PIX(img
, img_linesize
, x_f
-1, y_f
+1),
87 PIX(img
, img_linesize
, x_f
, y_f
+1),
88 PIX(img
, img_linesize
, x_f
+1, y_f
+1),
89 PIX(img
, img_linesize
, x_f
+2, y_f
+1));
90 short v4
= _FLT(bicub_kernel
)(tx
,
91 PIX(img
, img_linesize
, x_f
-1, y_f
+2),
92 PIX(img
, img_linesize
, x_f
, y_f
+2),
93 PIX(img
, img_linesize
, x_f
+1, y_f
+2),
94 PIX(img
, img_linesize
, x_f
+2, y_f
+2));
95 *rv
= (uint8_t)_FLT(bicub_kernel
)(y
-y_f
, v1
, v2
, v3
, v4
);
100 /** interpolateBiLin: bi-linear interpolation function, see interpolate */
101 void _FLT(interpolateBiLin
)(uint8_t *rv
, float x
, float y
,
102 const uint8_t *img
, int img_linesize
,
103 int width
, int height
, uint8_t def
)
105 if (x
< 0 || x
> width
- 1 || y
< 0 || y
> height
- 1) {
106 _FLT(interpolateBiLinBorder
)(rv
, x
, y
, img
, img_linesize
, width
, height
, def
);
108 int x_f
= myfloor(x
);
110 int y_f
= myfloor(y
);
112 short v1
= PIX(img
, img_linesize
, x_c
, y_c
);
113 short v2
= PIX(img
, img_linesize
, x_c
, y_f
);
114 short v3
= PIX(img
, img_linesize
, x_f
, y_c
);
115 short v4
= PIX(img
, img_linesize
, x_f
, y_f
);
116 float s
= (v1
*(x
- x_f
)+v3
*(x_c
- x
))*(y
- y_f
) +
117 (v2
*(x
- x_f
) + v4
*(x_c
- x
))*(y_c
- y
);
123 /** interpolateLin: linear (only x) interpolation function, see interpolate */
124 void _FLT(interpolateLin
)(uint8_t *rv
, float x
, float y
,
125 const uint8_t *img
, int img_linesize
,
126 int width
, int height
, uint8_t def
)
128 int x_f
= myfloor(x
);
130 int y_n
= myround(y
);
131 float v1
= PIXEL(img
, img_linesize
, x_c
, y_n
, width
, height
, def
);
132 float v2
= PIXEL(img
, img_linesize
, x_f
, y_n
, width
, height
, def
);
133 float s
= v1
*(x
- x_f
) + v2
*(x_c
- x
);
137 /** interpolateZero: nearest neighbor interpolation function, see interpolate */
138 void _FLT(interpolateZero
)(uint8_t *rv
, float x
, float y
,
139 const uint8_t *img
, int img_linesize
,
140 int width
, int height
, uint8_t def
)
142 int x_n
= myround(x
);
143 int y_n
= myround(y
);
144 *rv
= (uint8_t) PIXEL(img
, img_linesize
, x_n
, y_n
, width
, height
, def
);
149 * interpolateN: Bi-linear interpolation function for N channel image.
152 * rv: destination pixel (call by reference)
153 * x,y: the source coordinates in the image img. Note this
154 * are real-value coordinates, that's why we interpolate
156 * width,height: dimension of image
157 * N: number of channels
158 * channel: channel number (0..N-1)
159 * def: default value if coordinates are out of range
162 void _FLT(interpolateN
)(uint8_t *rv
, float x
, float y
,
163 const uint8_t *img
, int img_linesize
,
164 int width
, int height
,
165 uint8_t N
, uint8_t channel
,
168 if (x
< - 1 || x
> width
|| y
< -1 || y
> height
) {
171 int x_f
= myfloor(x
);
173 int y_f
= myfloor(y
);
175 short v1
= PIXELN(img
, img_linesize
, x_c
, y_c
, width
, height
, N
, channel
, def
);
176 short v2
= PIXELN(img
, img_linesize
, x_c
, y_f
, width
, height
, N
, channel
, def
);
177 short v3
= PIXELN(img
, img_linesize
, x_f
, y_c
, width
, height
, N
, channel
, def
);
178 short v4
= PIXELN(img
, img_linesize
, x_f
, y_f
, width
, height
, N
, channel
, def
);
179 float s
= (v1
*(x
- x_f
)+v3
*(x_c
- x
))*(y
- y_f
) +
180 (v2
*(x
- x_f
) + v4
*(x_c
- x
))*(y_c
- y
);
187 * transformPacked: applies current transformation to frame
189 * td: private data structure of this filter
191 * 0 for failture, 1 for success
193 * The frame must be in Packed format
195 /// Add bytes per pixel usage
197 int _FLT(transformPacked
)(VSTransformData
* td
, VSTransform t
)
199 int x
= 0, y
= 0, z
= 0;
201 char crop
= td
->conf
.crop
;
203 D_1
= td
->src
.data
[0];
204 D_2
= td
->destbuf
.data
[0];
205 float c_s_x
= td
->fiSrc
.width
/2.0;
206 float c_s_y
= td
->fiSrc
.height
/2.0;
207 float c_d_x
= td
->fiDest
.width
/2.0;
208 float c_d_y
= td
->fiDest
.height
/2.0;
210 /* for each pixel in the destination image we calc the source
211 * coordinate and make an interpolation:
212 * p_d = c_d + M(p_s - c_s) + t
213 * where p are the points, c the center coordinate,
214 * _s source and _d destination,
215 * t the translation, and M the rotation matrix
216 * p_s = M^{-1}(p_d - c_d - t) + c_s
218 int channels
= td
->fiSrc
.bytesPerPixel
;
220 if (fabs(t
.alpha
) > 0.1*M_PI
/180.0) { // 0.1 deg
221 for (x
= 0; x
< td
->fiDest
.width
; x
++) {
222 for (y
= 0; y
< td
->fiDest
.height
; y
++) {
223 float x_d1
= (x
- c_d_x
);
224 float y_d1
= (y
- c_d_y
);
225 float x_s
= cos(-t
.alpha
) * x_d1
226 + sin(-t
.alpha
) * y_d1
+ c_s_x
-t
.x
;
227 float y_s
= -sin(-t
.alpha
) * x_d1
228 + cos(-t
.alpha
) * y_d1
+ c_s_y
-t
.y
;
229 for (z
= 0; z
< channels
; z
++) { // iterate over colors
230 uint8_t *dest
= &D_2
[x
+ y
* td
->destbuf
.linesize
[0]+z
];
231 _FLT(interpolateN
)(dest
, x_s
, y_s
, D_1
, td
->src
.linesize
[0],
232 td
->fiSrc
.width
, td
->fiSrc
.height
,
233 channels
, z
, crop
? 16 : *dest
);
238 /* no rotation, just translation
239 *(also no interpolation, since no size change (so far)
241 int round_tx
= myround(t
.x
);
242 int round_ty
= myround(t
.y
);
243 for (x
= 0; x
< td
->fiDest
.width
; x
++) {
244 for (y
= 0; y
< td
->fiDest
.height
; y
++) {
245 for (z
= 0; z
< channels
; z
++) { // iterate over colors
246 short p
= PIXELN(D_1
, td
->src
.linesize
[0], x
- round_tx
, y
- round_ty
,
247 td
->fiSrc
.width
, td
->fiSrc
.height
, channels
, z
, -1);
250 D_2
[(x
+ y
* td
->destbuf
.linesize
[0])*channels
+z
] = 16;
252 D_2
[(x
+ y
* td
->destbuf
.linesize
[0])*channels
+z
] = (uint8_t)p
;
262 * transformPlanar: applies current transformation to frame
265 * td: private data structure of this filter
267 * 0 for failture, 1 for success
269 * The frame must be in Planar format
271 int _FLT(transformPlanar
)(VSTransformData
* td
, VSTransform t
)
274 uint8_t *dat_1
, *dat_2
;
275 char crop
= td
->conf
.crop
;
277 if (t
.alpha
==0 && t
.x
==0 && t
.y
==0 && t
.zoom
== 0){
278 if(vsFramesEqual(&td
->src
,&td
->destbuf
))
279 return VS_OK
; // noop
281 vsFrameCopy(&td
->destbuf
, &td
->src
, &td
->fiSrc
);
286 for(plane
=0; plane
< td
->fiSrc
.planes
; plane
++){
287 dat_1
= td
->src
.data
[plane
];
288 dat_2
= td
->destbuf
.data
[plane
];
290 int wsub
= vsGetPlaneWidthSubS(&td
->fiSrc
,plane
);
291 int hsub
= vsGetPlaneHeightSubS(&td
->fiSrc
,plane
);
292 float c_s_x
= (td
->fiSrc
.width
>> wsub
)/2.0;
293 float c_s_y
= (td
->fiSrc
.height
>> hsub
)/2.0;
294 float c_d_x
= (td
->fiDest
.width
>> wsub
)/2.0;
295 float c_d_y
= (td
->fiDest
.height
>> hsub
)/2.0;
296 uint8_t black
= plane
==0 ? 0 : 0x80;
298 float z
= 1.0-t
.zoom
/100;
299 float zcos_a
= z
*cos(-t
.alpha
); // scaled cos
300 float zsin_a
= z
*sin(-t
.alpha
); // scaled sin
301 float tx
= t
.x
/ (float)(1 << wsub
);
302 float ty
= t
.y
/ (float)(1 << hsub
);
304 /* for each pixel in the destination image we calc the source
305 * coordinate and make an interpolation:
306 * p_d = c_d + M(p_s - c_s) + t
307 * where p are the points, c the center coordinate,
308 * _s source and _d destination,
309 * t the translation, and M the rotation and scaling matrix
310 * p_s = M^{-1}(p_d - c_d - t) + c_s
312 int w
= CHROMA_SIZE(td
->fiDest
.width
,wsub
);
313 int h
= CHROMA_SIZE(td
->fiDest
.height
,hsub
);
314 int sw
= CHROMA_SIZE(td
->fiSrc
.width
,wsub
);
315 int sh
= CHROMA_SIZE(td
->fiSrc
.height
,hsub
);
316 for (x
= 0; x
< w
; x
++) {
317 for (y
= 0; y
< h
; y
++) {
318 float x_d1
= (x
- c_d_x
);
319 float y_d1
= (y
- c_d_y
);
320 float x_s
= zcos_a
* x_d1
321 + zsin_a
* y_d1
+ c_s_x
-tx
;
322 float y_s
= -zsin_a
* x_d1
323 + zcos_a
* y_d1
+ c_s_y
-ty
;
324 uint8_t *dest
= &dat_2
[x
+ y
* td
->destbuf
.linesize
[plane
]];
325 td
->_FLT(interpolate
)(dest
, x_s
, y_s
, dat_1
, td
->src
.linesize
[plane
],
326 sw
, sh
, crop
? black
: *dest
);
335 * c-file-style: "stroustrup"
336 * c-file-offsets: ((case-label . *) (statement-case-intro . *))
337 * indent-tabs-mode: nil
338 * c-basic-offset: 2 t
341 * vim: expandtab shiftwidth=2: