src/transformfloat.c

   1 /*
   2  *  transformfloat.c
   3  *
   4  *  Floating point image transformations
   5  *
   6  *  Copyright (C) Georg Martius - June 2011
   7  *   georg dot martius at web dot de
   8  *
   9  *  This file is part of vid.stab video stabilization library
  10  *
  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.
  15  *
  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.
  20  *
  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.
  24  *
  25  */
  26 #include "transformfloat.h"
  27 #include "transform.h"
  28 #include "transformtype_operations.h"
  29
  30
  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)
  36 {
  37   int x_f = myfloor(x);
  38   int x_c = x_f+1;
  39   int y_f = myfloor(y);
  40   int y_c = y_f+1;
  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);
  47   *rv = (uint8_t)s;
  48 }
  49
  50 /** taken from http://en.wikipedia.org/wiki/Bicubic_interpolation for alpha=-0.5
  51     in matrix notation:
  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
  54     | 0, 2, 0, 0 |  |a0|
  55     |-1, 0, 1, 0 |  |a1|
  56     (1,t,t^2,t^3) | 2,-5, 4,-1 |  |a2|
  57     |-1, 3,-3, 1 |  |a3|
  58 */
  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;
  61 }
  62
  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)
  67 {
  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);
  71   } else {
  72     int x_f = myfloor(x);
  73     int y_f = myfloor(y);
  74     float tx = x-x_f;
  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);
  96   }
  97 }
  98
  99
 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)
 104 {
 105   if (x < 0 || x > width - 1 || y < 0 || y > height - 1) {
 106     _FLT(interpolateBiLinBorder)(rv, x, y, img, img_linesize, width, height, def);
 107   } else {
 108     int x_f = myfloor(x);
 109     int x_c = x_f+1;
 110     int y_f = myfloor(y);
 111     int y_c = y_f+1;
 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);
 118     *rv = (uint8_t)s;
 119   }
 120 }
 121
 122
 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)
 127 {
 128   int x_f = myfloor(x);
 129   int x_c = x_f+1;
 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);
 134   *rv = (uint8_t)s;
 135 }
 136
 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)
 141 {
 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);
 145 }
 146
 147
 148 /**
 149  * interpolateN: Bi-linear interpolation function for N channel image.
 150  *
 151  * Parameters:
 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
 155  *            img: source image
 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
 160  * Return value:  None
 161  */
 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,
 166                         uint8_t def)
 167 {
 168   if (x < - 1 || x > width || y < -1 || y > height) {
 169     *rv = def;
 170   } else {
 171     int x_f = myfloor(x);
 172     int x_c = x_f+1;
 173     int y_f = myfloor(y);
 174     int y_c = y_f+1;
 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);
 181     *rv = (uint8_t)s;
 182   }
 183 }
 184
 185
 186 /**
 187  * transformPacked: applies current transformation to frame
 188  * Parameters:
 189  *         td: private data structure of this filter
 190  * Return value:
 191  *         0 for failture, 1 for success
 192  * Preconditions:
 193  *  The frame must be in Packed format
 194  /// TODO Add zoom!
 195  /// Add bytes per pixel usage
 196  */
 197 int _FLT(transformPacked)(VSTransformData* td, VSTransform t)
 198 {
 199   int x = 0, y = 0, z = 0;
 200   uint8_t *D_1, *D_2;
 201   char crop = td->conf.crop;
 202
 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;
 209
 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
 217    */
 218   int channels = td->fiSrc.bytesPerPixel;
 219   /* All channels */
 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);
 234         }
 235       }
 236     }
 237   }else {
 238     /* no rotation, just translation
 239      *(also no interpolation, since no size change (so far)
 240      */
 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);
 248           if (p == -1) {
 249             if (crop == 1)
 250               D_2[(x + y * td->destbuf.linesize[0])*channels+z] = 16;
 251           } else {
 252             D_2[(x + y * td->destbuf.linesize[0])*channels+z] = (uint8_t)p;
 253           }
 254         }
 255       }
 256     }
 257   }
 258   return 1;
 259 }
 260
 261 /**
 262  * transformPlanar: applies current transformation to frame
 263  *
 264  * Parameters:
 265  *         td: private data structure of this filter
 266  * Return value:
 267  *         0 for failture, 1 for success
 268  * Preconditions:
 269  *  The frame must be in Planar format
 270  */
 271 int _FLT(transformPlanar)(VSTransformData* td, VSTransform t)
 272 {
 273   int x = 0, y = 0;
 274   uint8_t *dat_1, *dat_2;
 275   char crop = td->conf.crop;
 276
 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
 280     else {
 281       vsFrameCopy(&td->destbuf, &td->src, &td->fiSrc);
 282       return VS_OK;
 283     }
 284   }
 285   int plane;
 286   for(plane=0; plane< td->fiSrc.planes; plane++){
 287     dat_1  = td->src.data[plane];
 288     dat_2  = td->destbuf.data[plane];
 289
 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;
 297
 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);
 303
 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
 311      */
 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);
 327       }
 328     }
 329   }
 330   return VS_OK;
 331 }
 332
 333 /*
 334  * Local variables:
 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
 339  * End:
 340  *
 341  * vim: expandtab shiftwidth=2:
 342  */