git.sesse.net Git - ffmpeg/blob - libavfilter/vf_deshake.c

   1 /*
   2  * Copyright (C) 2010 Georg Martius <georg.martius@web.de>
   3  * Copyright (C) 2010 Daniel G. Taylor <dan@programmer-art.org>
   4  *
   5  * This file is part of FFmpeg.
   6  *
   7  * FFmpeg is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * FFmpeg is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with FFmpeg; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21
  22 /**
  23  * @file
  24  * fast deshake / depan video filter
  25  *
  26  * SAD block-matching motion compensation to fix small changes in
  27  * horizontal and/or vertical shift. This filter helps remove camera shake
  28  * from hand-holding a camera, bumping a tripod, moving on a vehicle, etc.
  29  *
  30  * Algorithm:
  31  *   - For each frame with one previous reference frame
  32  *       - For each block in the frame
  33  *           - If contrast > threshold then find likely motion vector
  34  *       - For all found motion vectors
  35  *           - Find most common, store as global motion vector
  36  *       - Find most likely rotation angle
  37  *       - Transform image along global motion
  38  *
  39  * TODO:
  40  *   - Fill frame edges based on previous/next reference frames
  41  *   - Fill frame edges by stretching image near the edges?
  42  *       - Can this be done quickly and look decent?
  43  *
  44  * Dark Shikari links to http://wiki.videolan.org/SoC_x264_2010#GPU_Motion_Estimation_2
  45  * for an algorithm similar to what could be used here to get the gmv
  46  * It requires only a couple diamond searches + fast downscaling
  47  *
  48  * Special thanks to Jason Kotenko for his help with the algorithm and my
  49  * inability to see simple errors in C code.
  50  */
  51
  52 #include "avfilter.h"
  53 #include "libavutil/common.h"
  54 #include "libavutil/mem.h"
  55 #include "libavutil/pixdesc.h"
  56 #include "libavcodec/dsputil.h"
  57
  58 #include "transform.h"
  59
  60 #define CHROMA_WIDTH(link)  -((-link->w) >> av_pix_fmt_descriptors[link->format].log2_chroma_w)
  61 #define CHROMA_HEIGHT(link) -((-link->h) >> av_pix_fmt_descriptors[link->format].log2_chroma_h)
  62
  63 enum SearchMethod {
  64     EXHAUSTIVE,        ///< Search all possible positions
  65     SMART_EXHAUSTIVE,  ///< Search most possible positions (faster)
  66     SEARCH_COUNT
  67 };
  68
  69 typedef struct {
  70     int x;             ///< Horizontal shift
  71     int y;             ///< Vertical shift
  72 } IntMotionVector;
  73
  74 typedef struct {
  75     double x;             ///< Horizontal shift
  76     double y;             ///< Vertical shift
  77 } MotionVector;
  78
  79 typedef struct {
  80     MotionVector vector;  ///< Motion vector
  81     double angle;         ///< Angle of rotation
  82     double zoom;          ///< Zoom percentage
  83 } Transform;
  84
  85 typedef struct {
  86     AVClass av_class;
  87     AVFilterBufferRef *ref;    ///< Previous frame
  88     int rx;                    ///< Maximum horizontal shift
  89     int ry;                    ///< Maximum vertical shift
  90     enum FillMethod edge;      ///< Edge fill method
  91     int blocksize;             ///< Size of blocks to compare
  92     int contrast;              ///< Contrast threshold
  93     enum SearchMethod search;  ///< Motion search method
  94     AVCodecContext *avctx;
  95     DSPContext c;              ///< Context providing optimized SAD methods
  96     Transform last;            ///< Transform from last frame
  97     int refcount;              ///< Number of reference frames (defines averaging window)
  98     FILE *fp;
  99     Transform avg;
 100     int cw;                    ///< Crop motion search to this box
 101     int ch;
 102     int cx;
 103     int cy;
 104 } DeshakeContext;
 105
 106 static int cmp(const double *a, const double *b)
 107 {
 108     return *a < *b ? -1 : ( *a > *b ? 1 : 0 );
 109 }
 110
 111 /**
 112  * Cleaned mean (cuts off 20% of values to remove outliers and then averages)
 113  */
 114 static double clean_mean(double *values, int count)
 115 {
 116     double mean = 0;
 117     int cut = count / 5;
 118     int x;
 119
 120     qsort(values, count, sizeof(double), (void*)cmp);
 121
 122     for (x = cut; x < count - cut; x++) {
 123         mean += values[x];
 124     }
 125
 126     return mean / (count - cut * 2);
 127 }
 128
 129 /**
 130  * Find the most likely shift in motion between two frames for a given
 131  * macroblock. Test each block against several shifts given by the rx
 132  * and ry attributes. Searches using a simple matrix of those shifts and
 133  * chooses the most likely shift by the smallest difference in blocks.
 134  */
 135 static void find_block_motion(DeshakeContext *deshake, uint8_t *src1,
 136                               uint8_t *src2, int cx, int cy, int stride,
 137                               IntMotionVector *mv)
 138 {
 139     int x, y;
 140     int diff;
 141     int smallest = INT_MAX;
 142     int tmp, tmp2;
 143
 144     #define CMP(i, j) deshake->c.sad[0](deshake, src1 + cy * stride + cx, \
 145                                         src2 + (j) * stride + (i), stride, \
 146                                         deshake->blocksize)
 147
 148     if (deshake->search == EXHAUSTIVE) {
 149         // Compare every possible position - this is sloooow!
 150         for (y = -deshake->ry; y <= deshake->ry; y++) {
 151             for (x = -deshake->rx; x <= deshake->rx; x++) {
 152                 diff = CMP(cx - x, cy - y);
 153                 if (diff < smallest) {
 154                     smallest = diff;
 155                     mv->x = x;
 156                     mv->y = y;
 157                 }
 158             }
 159         }
 160     } else if (deshake->search == SMART_EXHAUSTIVE) {
 161         // Compare every other possible position and find the best match
 162         for (y = -deshake->ry + 1; y < deshake->ry - 2; y += 2) {
 163             for (x = -deshake->rx + 1; x < deshake->rx - 2; x += 2) {
 164                 diff = CMP(cx - x, cy - y);
 165                 if (diff < smallest) {
 166                     smallest = diff;
 167                     mv->x = x;
 168                     mv->y = y;
 169                 }
 170             }
 171         }
 172
 173         // Hone in on the specific best match around the match we found above
 174         tmp = mv->x;
 175         tmp2 = mv->y;
 176
 177         for (y = tmp2 - 1; y <= tmp2 + 1; y++) {
 178             for (x = tmp - 1; x <= tmp + 1; x++) {
 179                 if (x == tmp && y == tmp2)
 180                     continue;
 181
 182                 diff = CMP(cx - x, cy - y);
 183                 if (diff < smallest) {
 184                     smallest = diff;
 185                     mv->x = x;
 186                     mv->y = y;
 187                 }
 188             }
 189         }
 190     }
 191
 192     if (smallest > 512) {
 193         mv->x = -1;
 194         mv->y = -1;
 195     }
 196     emms_c();
 197     //av_log(NULL, AV_LOG_ERROR, "%d\n", smallest);
 198     //av_log(NULL, AV_LOG_ERROR, "Final: (%d, %d) = %d x %d\n", cx, cy, mv->x, mv->y);
 199 }
 200
 201 /**
 202  * Find the contrast of a given block. When searching for global motion we
 203  * really only care about the high contrast blocks, so using this method we
 204  * can actually skip blocks we don't care much about.
 205  */
 206 static int block_contrast(uint8_t *src, int x, int y, int stride, int blocksize)
 207 {
 208     int highest = 0;
 209     int lowest = 0;
 210     int i, j, pos;
 211
 212     for (i = 0; i <= blocksize * 2; i++) {
 213         // We use a width of 16 here to match the libavcodec sad functions
 214         for (j = 0; i <= 15; i++) {
 215             pos = (y - i) * stride + (x - j);
 216             if (src[pos] < lowest)
 217                 lowest = src[pos];
 218             else if (src[pos] > highest) {
 219                 highest = src[pos];
 220             }
 221         }
 222     }
 223
 224     return highest - lowest;
 225 }
 226
 227 /**
 228  * Find the rotation for a given block.
 229  */
 230 static double block_angle(int x, int y, int cx, int cy, IntMotionVector *shift)
 231 {
 232     double a1, a2, diff;
 233
 234     a1 = atan2(y - cy, x - cx);
 235     a2 = atan2(y - cy + shift->y, x - cx + shift->x);
 236
 237     diff = a2 - a1;
 238
 239     return (diff > M_PI)  ? diff - 2 * M_PI :
 240            (diff < -M_PI) ? diff + 2 * M_PI :
 241            diff;
 242 }
 243
 244 /**
 245  * Find the estimated global motion for a scene given the most likely shift
 246  * for each block in the frame. The global motion is estimated to be the
 247  * same as the motion from most blocks in the frame, so if most blocks
 248  * move one pixel to the right and two pixels down, this would yield a
 249  * motion vector (1, -2).
 250  */
 251 static void find_motion(DeshakeContext *deshake, uint8_t *src1, uint8_t *src2,
 252                         int width, int height, int stride, Transform *t)
 253 {
 254     int x, y;
 255     IntMotionVector mv = {0, 0};
 256     int counts[128][128];
 257     int count_max_value = 0;
 258     int contrast;
 259
 260     int pos;
 261     double *angles = av_malloc(sizeof(*angles) * width * height / (16 * deshake->blocksize));
 262     int center_x = 0, center_y = 0;
 263     double p_x, p_y;
 264
 265     // Reset counts to zero
 266     for (x = 0; x < deshake->rx * 2 + 1; x++) {
 267         for (y = 0; y < deshake->ry * 2 + 1; y++) {
 268             counts[x][y] = 0;
 269         }
 270     }
 271
 272     pos = 0;
 273     // Find motion for every block and store the motion vector in the counts
 274     for (y = deshake->ry; y < height - deshake->ry - (deshake->blocksize * 2); y += deshake->blocksize * 2) {
 275         // We use a width of 16 here to match the libavcodec sad functions
 276         for (x = deshake->rx; x < width - deshake->rx - 16; x += 16) {
 277             // If the contrast is too low, just skip this block as it probably
 278             // won't be very useful to us.
 279             contrast = block_contrast(src2, x, y, stride, deshake->blocksize);
 280             if (contrast > deshake->contrast) {
 281                 //av_log(NULL, AV_LOG_ERROR, "%d\n", contrast);
 282                 find_block_motion(deshake, src1, src2, x, y, stride, &mv);
 283                 if (mv.x != -1 && mv.y != -1) {
 284                     counts[mv.x + deshake->rx][mv.y + deshake->ry] += 1;
 285                     if (x > deshake->rx && y > deshake->ry)
 286                         angles[pos++] = block_angle(x, y, 0, 0, &mv);
 287
 288                     center_x += mv.x;
 289                     center_y += mv.y;
 290                 }
 291             }
 292         }
 293     }
 294
 295     pos = FFMAX(1, pos);
 296
 297     center_x /= pos;
 298     center_y /= pos;
 299
 300     t->angle = clean_mean(angles, pos);
 301     if (t->angle < 0.001)
 302         t->angle = 0;
 303
 304     // Find the most common motion vector in the frame and use it as the gmv
 305     for (y = deshake->ry * 2; y >= 0; y--) {
 306         for (x = 0; x < deshake->rx * 2 + 1; x++) {
 307             //av_log(NULL, AV_LOG_ERROR, "%5d ", counts[x][y]);
 308             if (counts[x][y] > count_max_value) {
 309                 t->vector.x = x - deshake->rx;
 310                 t->vector.y = y - deshake->ry;
 311                 count_max_value = counts[x][y];
 312             }
 313         }
 314         //av_log(NULL, AV_LOG_ERROR, "\n");
 315     }
 316
 317     p_x = (center_x - width / 2);
 318     p_y = (center_y - height / 2);
 319     t->vector.x += (cos(t->angle)-1)*p_x  - sin(t->angle)*p_y;
 320     t->vector.y += sin(t->angle)*p_x  + (cos(t->angle)-1)*p_y;
 321
 322     // Clamp max shift & rotation?
 323     t->vector.x = av_clipf(t->vector.x, -deshake->rx * 2, deshake->rx * 2);
 324     t->vector.y = av_clipf(t->vector.y, -deshake->ry * 2, deshake->ry * 2);
 325     t->angle = av_clipf(t->angle, -0.1, 0.1);
 326
 327     //av_log(NULL, AV_LOG_ERROR, "%d x %d\n", avg->x, avg->y);
 328     av_free(angles);
 329 }
 330
 331 static av_cold int init(AVFilterContext *ctx, const char *args, void *opaque)
 332 {
 333     DeshakeContext *deshake = ctx->priv;
 334     char filename[256] = {0};
 335
 336     deshake->rx = 16;
 337     deshake->ry = 16;
 338     deshake->edge = FILL_MIRROR;
 339     deshake->blocksize = 8;
 340     deshake->contrast = 125;
 341     deshake->search = EXHAUSTIVE;
 342     deshake->refcount = 20;
 343
 344     deshake->cw = -1;
 345     deshake->ch = -1;
 346     deshake->cx = -1;
 347     deshake->cy = -1;
 348
 349     if (args) {
 350         sscanf(args, "%d:%d:%d:%d:%d:%d:%d:%d:%d:%d:%255s",
 351                &deshake->cx, &deshake->cy, &deshake->cw, &deshake->ch,
 352                &deshake->rx, &deshake->ry, (int *)&deshake->edge,
 353                &deshake->blocksize, &deshake->contrast, (int *)&deshake->search, filename);
 354
 355         deshake->blocksize /= 2;
 356
 357         deshake->rx = av_clip(deshake->rx, 0, 64);
 358         deshake->ry = av_clip(deshake->ry, 0, 64);
 359         deshake->edge = av_clip(deshake->edge, FILL_BLANK, FILL_COUNT - 1);
 360         deshake->blocksize = av_clip(deshake->blocksize, 4, 128);
 361         deshake->contrast = av_clip(deshake->contrast, 1, 255);
 362         deshake->search = av_clip(deshake->search, EXHAUSTIVE, SEARCH_COUNT - 1);
 363
 364     }
 365     if (*filename)
 366         deshake->fp = fopen(filename, "w");
 367     if (deshake->fp)
 368         fwrite("Ori x, Avg x, Fin x, Ori y, Avg y, Fin y, Ori angle, Avg angle, Fin angle, Ori zoom, Avg zoom, Fin zoom\n", sizeof(char), 104, deshake->fp);
 369
 370     // Quadword align left edge of box for MMX code, adjust width if necessary
 371     // to keep right margin
 372     if (deshake->cx > 0) {
 373         deshake->cw += deshake->cx - (deshake->cx & ~15);
 374         deshake->cx &= ~15;
 375     }
 376
 377     av_log(ctx, AV_LOG_INFO, "cx: %d, cy: %d, cw: %d, ch: %d, rx: %d, ry: %d, edge: %d blocksize: %d contrast: %d search: %d\n",
 378            deshake->cx, deshake->cy, deshake->cw, deshake->ch,
 379            deshake->rx, deshake->ry, deshake->edge, deshake->blocksize * 2, deshake->contrast, deshake->search);
 380
 381     return 0;
 382 }
 383
 384 static int query_formats(AVFilterContext *ctx)
 385 {
 386     enum PixelFormat pix_fmts[] = {
 387         PIX_FMT_YUV420P,  PIX_FMT_YUV422P,  PIX_FMT_YUV444P,  PIX_FMT_YUV410P,
 388         PIX_FMT_YUV411P,  PIX_FMT_YUV440P,  PIX_FMT_YUVJ420P, PIX_FMT_YUVJ422P,
 389         PIX_FMT_YUVJ444P, PIX_FMT_YUVJ440P, PIX_FMT_NONE
 390     };
 391
 392     avfilter_set_common_pixel_formats(ctx, avfilter_make_format_list(pix_fmts));
 393
 394     return 0;
 395 }
 396
 397 static int config_props(AVFilterLink *link)
 398 {
 399     DeshakeContext *deshake = link->dst->priv;
 400
 401     deshake->ref = NULL;
 402     deshake->last.vector.x = 0;
 403     deshake->last.vector.y = 0;
 404     deshake->last.angle = 0;
 405     deshake->last.zoom = 0;
 406
 407     deshake->avctx = avcodec_alloc_context3(NULL);
 408     dsputil_init(&deshake->c, deshake->avctx);
 409
 410     return 0;
 411 }
 412
 413 static av_cold void uninit(AVFilterContext *ctx)
 414 {
 415     DeshakeContext *deshake = ctx->priv;
 416
 417     avfilter_unref_buffer(deshake->ref);
 418     if (deshake->fp)
 419         fclose(deshake->fp);
 420 }
 421
 422 static void end_frame(AVFilterLink *link)
 423 {
 424     DeshakeContext *deshake = link->dst->priv;
 425     AVFilterBufferRef *in  = link->cur_buf;
 426     AVFilterBufferRef *out = link->dst->outputs[0]->out_buf;
 427     Transform t = {{0},0}, orig = {{0},0};
 428     float matrix[9];
 429     float alpha = 2.0 / deshake->refcount;
 430     char tmp[256];
 431
 432     if (deshake->cx < 0 || deshake->cy < 0 || deshake->cw < 0 || deshake->ch < 0) {
 433         // Find the most likely global motion for the current frame
 434         find_motion(deshake, (deshake->ref == NULL) ? in->data[0] : deshake->ref->data[0], in->data[0], link->w, link->h, in->linesize[0], &t);
 435     } else {
 436         uint8_t *src1 = (deshake->ref == NULL) ? in->data[0] : deshake->ref->data[0];
 437         uint8_t *src2 = in->data[0];
 438
 439         deshake->cx = FFMIN(deshake->cx, link->w);
 440         deshake->cy = FFMIN(deshake->cy, link->h);
 441
 442         if ((unsigned)deshake->cx + (unsigned)deshake->cw > link->w) deshake->cw = link->w - deshake->cx;
 443         if ((unsigned)deshake->cy + (unsigned)deshake->ch > link->h) deshake->ch = link->h - deshake->cy;
 444
 445         // Quadword align right margin
 446         deshake->cw &= ~15;
 447
 448         src1 += deshake->cy * in->linesize[0] + deshake->cx;
 449         src2 += deshake->cy * in->linesize[0] + deshake->cx;
 450
 451         find_motion(deshake, src1, src2, deshake->cw, deshake->ch, in->linesize[0], &t);
 452     }
 453
 454
 455     // Copy transform so we can output it later to compare to the smoothed value
 456     orig.vector.x = t.vector.x;
 457     orig.vector.y = t.vector.y;
 458     orig.angle = t.angle;
 459     orig.zoom = t.zoom;
 460
 461     // Generate a one-sided moving exponential average
 462     deshake->avg.vector.x = alpha * t.vector.x + (1.0 - alpha) * deshake->avg.vector.x;
 463     deshake->avg.vector.y = alpha * t.vector.y + (1.0 - alpha) * deshake->avg.vector.y;
 464     deshake->avg.angle = alpha * t.angle + (1.0 - alpha) * deshake->avg.angle;
 465     deshake->avg.zoom = alpha * t.zoom + (1.0 - alpha) * deshake->avg.zoom;
 466
 467     // Remove the average from the current motion to detect the motion that
 468     // is not on purpose, just as jitter from bumping the camera
 469     t.vector.x -= deshake->avg.vector.x;
 470     t.vector.y -= deshake->avg.vector.y;
 471     t.angle -= deshake->avg.angle;
 472     t.zoom -= deshake->avg.zoom;
 473
 474     // Invert the motion to undo it
 475     t.vector.x *= -1;
 476     t.vector.y *= -1;
 477     t.angle *= -1;
 478
 479     // Write statistics to file
 480     if (deshake->fp) {
 481         snprintf(tmp, 256, "%f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f\n", orig.vector.x, deshake->avg.vector.x, t.vector.x, orig.vector.y, deshake->avg.vector.y, t.vector.y, orig.angle, deshake->avg.angle, t.angle, orig.zoom, deshake->avg.zoom, t.zoom);
 482         fwrite(tmp, sizeof(char), strlen(tmp), deshake->fp);
 483     }
 484
 485     // Turn relative current frame motion into absolute by adding it to the
 486     // last absolute motion
 487     t.vector.x += deshake->last.vector.x;
 488     t.vector.y += deshake->last.vector.y;
 489     t.angle += deshake->last.angle;
 490     t.zoom += deshake->last.zoom;
 491
 492     // Shrink motion by 10% to keep things centered in the camera frame
 493     t.vector.x *= 0.9;
 494     t.vector.y *= 0.9;
 495     t.angle *= 0.9;
 496
 497     // Store the last absolute motion information
 498     deshake->last.vector.x = t.vector.x;
 499     deshake->last.vector.y = t.vector.y;
 500     deshake->last.angle = t.angle;
 501     deshake->last.zoom = t.zoom;
 502
 503     // Generate a luma transformation matrix
 504     avfilter_get_matrix(t.vector.x, t.vector.y, t.angle, 1.0 + t.zoom / 100.0, matrix);
 505
 506     // Transform the luma plane
 507     avfilter_transform(in->data[0], out->data[0], in->linesize[0], out->linesize[0], link->w, link->h, matrix, INTERPOLATE_BILINEAR, deshake->edge);
 508
 509     // Generate a chroma transformation matrix
 510     avfilter_get_matrix(t.vector.x / (link->w / CHROMA_WIDTH(link)), t.vector.y / (link->h / CHROMA_HEIGHT(link)), t.angle, 1.0 + t.zoom / 100.0, matrix);
 511
 512     // Transform the chroma planes
 513     avfilter_transform(in->data[1], out->data[1], in->linesize[1], out->linesize[1], CHROMA_WIDTH(link), CHROMA_HEIGHT(link), matrix, INTERPOLATE_BILINEAR, deshake->edge);
 514     avfilter_transform(in->data[2], out->data[2], in->linesize[2], out->linesize[2], CHROMA_WIDTH(link), CHROMA_HEIGHT(link), matrix, INTERPOLATE_BILINEAR, deshake->edge);
 515
 516     // Store the current frame as the reference frame for calculating the
 517     // motion of the next frame
 518     if (deshake->ref != NULL)
 519         avfilter_unref_buffer(deshake->ref);
 520
 521     // Cleanup the old reference frame
 522     deshake->ref = in;
 523
 524     // Draw the transformed frame information
 525     avfilter_draw_slice(link->dst->outputs[0], 0, link->h, 1);
 526     avfilter_end_frame(link->dst->outputs[0]);
 527     avfilter_unref_buffer(out);
 528 }
 529
 530 static void draw_slice(AVFilterLink *link, int y, int h, int slice_dir)
 531 {
 532 }
 533
 534 AVFilter avfilter_vf_deshake = {
 535     .name      = "deshake",
 536     .description = NULL_IF_CONFIG_SMALL("Stabilize shaky video."),
 537
 538     .priv_size = sizeof(DeshakeContext),
 539
 540     .init = init,
 541     .uninit = uninit,
 542     .query_formats = query_formats,
 543
 544     .inputs    = (const AVFilterPad[]) {{ .name       = "default",
 545                                     .type             = AVMEDIA_TYPE_VIDEO,
 546                                     .draw_slice       = draw_slice,
 547                                     .end_frame        = end_frame,
 548                                     .config_props     = config_props,
 549                                     .min_perms        = AV_PERM_READ, },
 550                                   { .name = NULL}},
 551
 552     .outputs   = (const AVFilterPad[]) {{ .name       = "default",
 553                                     .type             = AVMEDIA_TYPE_VIDEO, },
 554                                   { .name = NULL}},
 555 };