git.sesse.net Git - mlt/blob - src/modules/videostab/stabilize.c

   1 /*
   2  *  filter_stabilize.c
   3  *
   4  *  Copyright (C) Georg Martius - June 2007
   5  *   georg dot martius at web dot de
   6  *
   7  *  This file is part of transcode, a video stream processing tool
   8  *
   9  *  transcode is free software; you can redistribute it and/or modify
  10  *  it under the terms of the GNU General Public License as published by
  11  *  the Free Software Foundation; either version 2, or (at your option)
  12  *  any later version.
  13  *
  14  *  transcode is distributed in the hope that it will be useful,
  15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  *  GNU General Public License for more details.
  18  *
  19  *  You should have received a copy of the GNU General Public License
  20  *  along with GNU Make; see the file COPYING.  If not, write to
  21  *  the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  22  *
  23  */
  24
  25 /* Typical call:
  26  *  transcode -V -J stabilize=shakiness=5:show=1,preview
  27  *         -i inp.mpeg -y null,null -o dummy
  28  *  all parameters are optional
  29 */
  30
  31 #define MOD_NAME    "filter_stabilize.so"
  32 #define MOD_VERSION "v0.75 (2010-04-07)"
  33 #define MOD_CAP     "extracts relative transformations of \n\
  34     subsequent frames (used for stabilization together with the\n\
  35     transform filter in a second pass)"
  36 #define MOD_AUTHOR  "Georg Martius"
  37
  38 /* Ideas:
  39  - Try OpenCL/Cuda, this should work great
  40  - use smoothing on the frames and then use gradient decent!
  41  - stepsize could be adapted (maybe to check only one field with large
  42    stepsize and use the maximally required for the other fields
  43 */
  44
  45 #define MOD_FEATURES \
  46     TC_MODULE_FEATURE_FILTER|TC_MODULE_FEATURE_VIDEO
  47 #define MOD_FLAGS  \
  48     TC_MODULE_FLAG_RECONFIGURABLE | TC_MODULE_FLAG_DELAY
  49 #define MAX(a,b)         ((a < b) ?  (b) : (a))
  50 #define MIN(a,b)         ((a < b) ?  (a) : (b))
  51 #include "stabilize.h"
  52 #include <stdlib.h>
  53 #include <string.h>
  54 #include <framework/mlt_types.h>
  55
  56
  57 void addTrans(StabData* sd, Transform sl)
  58 {
  59     if (!sd->transs) {
  60        sd->transs = tlist_new(0);
  61     }
  62     tlist_append(sd->transs, &sl,sizeof(Transform) );
  63 }
  64
  65
  66
  67 /** initialise measurement fields on the frame.
  68     The size of the fields and the maxshift is used to
  69     calculate an optimal distribution in the frame.
  70 */
  71 int initFields(StabData* sd)
  72 {
  73     int size     = sd->field_size;
  74     int rows = MAX(3,(sd->height - sd->maxshift*2)/size-1);
  75     int cols = MAX(3,(sd->width  - sd->maxshift*2)/size-1);
  76     // make sure that the remaining rows have the same length
  77     sd->field_num  = rows*cols;
  78     sd->field_rows = rows;
  79     printf("field setup: rows: %i cols: %i Total: %i fields",
  80                 rows, cols, sd->field_num);
  81
  82     if (!(sd->fields = malloc(sizeof(Field) * sd->field_num))) {
  83         printf( "malloc failed!\n");
  84         return 0;
  85     } else {
  86         int i, j;
  87         // the border is the amount by which the field centers
  88         // have to be away from the image boundary
  89         // (stepsize is added in case shift is increased through stepsize)
  90         int border   = size/2 + sd->maxshift + sd->stepsize;
  91         int step_x   = (sd->width  - 2*border)/MAX(cols-1,1);
  92         int step_y   = (sd->height - 2*border) / MAX(rows-1,1);
  93         for (j = 0; j < rows; j++) {
  94             for (i = 0; i < cols; i++) {
  95                 int idx = j*cols+i;
  96                 sd->fields[idx].x = border + i*step_x;
  97                 sd->fields[idx].y = border + j*step_y;
  98                 sd->fields[idx].size = size;
  99             }
 100         }
 101     }
 102     return 1;
 103 }
 104
 105
 106 /**
 107    compares the two given images and returns the average absolute difference
 108    \param d_x shift in x direction
 109    \param d_y shift in y direction
 110 */
 111 double compareImg(unsigned char* I1, unsigned char* I2,
 112                   int width, int height,  int bytesPerPixel, int d_x, int d_y)
 113 {
 114     int i, j;
 115     unsigned char* p1 = NULL;
 116     unsigned char* p2 = NULL;
 117     long int sum = 0;
 118     int effectWidth = width - abs(d_x);
 119     int effectHeight = height - abs(d_y);
 120
 121 /*   DEBUGGING code to export single frames */
 122 /*   char buffer[100]; */
 123 /*   sprintf(buffer, "pic_%02ix%02i_1.ppm", d_x, d_y); */
 124 /*   FILE *pic1 = fopen(buffer, "w"); */
 125 /*   sprintf(buffer, "pic_%02ix%02i_2.ppm", d_x, d_y); */
 126 /*   FILE *pic2 = fopen(buffer, "w"); */
 127 /*   fprintf(pic1, "P6\n%i %i\n255\n", effectWidth, effectHeight); */
 128 /*   fprintf(pic2, "P6\n%i %i\n255\n", effectWidth, effectHeight); */
 129
 130     for (i = 0; i < effectHeight; i++) {
 131         p1 = I1;
 132         p2 = I2;
 133         if (d_y > 0 ){
 134             p1 += (i + d_y) * width * bytesPerPixel;
 135             p2 += i * width * bytesPerPixel;
 136         } else {
 137             p1 += i * width * bytesPerPixel;
 138             p2 += (i - d_y) * width * bytesPerPixel;
 139         }
 140         if (d_x > 0) {
 141             p1 += d_x * bytesPerPixel;
 142         } else {
 143             p2 -= d_x * bytesPerPixel;
 144         }
 145         // TODO: use some mmx or sse stuff here
 146         for (j = 0; j < effectWidth * bytesPerPixel; j++) {
 147             /* debugging code continued */
 148             /* fwrite(p1,1,1,pic1);fwrite(p1,1,1,pic1);fwrite(p1,1,1,pic1);
 149                fwrite(p2,1,1,pic2);fwrite(p2,1,1,pic2);fwrite(p2,1,1,pic2);
 150              */
 151             sum += abs((int)*p1 - (int)*p2);
 152             p1++;
 153             p2++;
 154         }
 155     }
 156     /*  fclose(pic1);
 157         fclose(pic2);
 158      */
 159     return sum/((double) effectWidth * effectHeight * bytesPerPixel);
 160 }
 161
 162 /**
 163    compares a small part of two given images
 164    and returns the average absolute difference.
 165    Field center, size and shift have to be choosen,
 166    so that no clipping is required
 167
 168    \param field Field specifies position(center) and size of subimage
 169    \param d_x shift in x direction
 170    \param d_y shift in y direction
 171 */
 172 double compareSubImg(unsigned char* const I1, unsigned char* const I2,
 173                      const Field* field,
 174                      int width, int height, int bytesPerPixel, int d_x, int d_y)
 175 {
 176     int k, j;
 177     unsigned char* p1 = NULL;
 178     unsigned char* p2 = NULL;
 179     int s2 = field->size / 2;
 180     double sum = 0;
 181
 182     p1=I1 + ((field->x - s2) + (field->y - s2)*width)*bytesPerPixel;
 183     p2=I2 + ((field->x - s2 + d_x) + (field->y - s2 + d_y)*width)*bytesPerPixel;
 184     // TODO: use some mmx or sse stuff here
 185     for (j = 0; j < field->size; j++){
 186
 187         for (k = 0; k < field->size * bytesPerPixel; k++) {
 188 #if HAVE_MMX
 189 #error "wo"
 190 #else
 191             sum += abs((int)*p1 - (int)*p2);
 192                         /*asm {
 193                                 push edx
 194
 195                                 pop edx
 196                                 emms
 197                         }*/
 198 #endif
 199             p1++;
 200             p2++;
 201         }
 202         p1 += (width - field->size) * bytesPerPixel;
 203         p2 += (width - field->size) * bytesPerPixel;
 204     }
 205     return sum/((double) field->size *field->size* bytesPerPixel);
 206 }
 207
 208 /** \see contrastSubImg called with bytesPerPixel=1*/
 209 double contrastSubImgYUV(StabData* sd, const Field* field){
 210     return contrastSubImg(sd->curr,field,sd->width,sd->height,1);
 211 }
 212
 213 /**
 214     \see contrastSubImg three times called with bytesPerPixel=3
 215     for all channels
 216  */
 217 double contrastSubImgRGB(StabData* sd, const Field* field){
 218     unsigned char* const I = sd->curr;
 219     return (  contrastSubImg(I,  field,sd->width,sd->height,3)
 220             + contrastSubImg(I+1,field,sd->width,sd->height,3)
 221             + contrastSubImg(I+2,field,sd->width,sd->height,3))/3;
 222 }
 223
 224 /**
 225    calculates Michelson-contrast in the given small part of the given image
 226
 227    \param I pointer to framebuffer
 228    \param field Field specifies position(center) and size of subimage
 229    \param width width of frame
 230    \param height height of frame
 231    \param bytesPerPixel calc contrast for only for first channel
 232 */
 233 double contrastSubImg(unsigned char* const I, const Field* field,
 234                      int width, int height, int bytesPerPixel)
 235 {
 236     int k, j;
 237     unsigned char* p = NULL;
 238     int s2 = field->size / 2;
 239     unsigned char mini = 255;
 240     unsigned char maxi = 0;
 241
 242     p = I + ((field->x - s2) + (field->y - s2)*width)*bytesPerPixel;
 243     // TODO: use some mmx or sse stuff here
 244     for (j = 0; j < field->size; j++){
 245         for (k = 0; k < field->size * bytesPerPixel; k++) {
 246             mini = (mini < *p) ? mini : *p;
 247             maxi = (maxi > *p) ? maxi : *p;
 248             p += bytesPerPixel;
 249         }
 250         p += (width - field->size) * bytesPerPixel;
 251     }
 252     return (maxi-mini)/(maxi+mini+0.1); // +0.1 to avoid division by 0
 253 }
 254
 255 /** tries to register current frame onto previous frame.
 256     This is the most simple algorithm:
 257     shift images to all possible positions and calc summed error
 258     Shift with minimal error is selected.
 259 */
 260 Transform calcShiftRGBSimple(StabData* sd)
 261 {
 262         printf("calc shoft rgb simple\n");
 263     int x = 0, y = 0;
 264     int i, j;
 265     double minerror = 1e20;
 266     for (i = -sd->maxshift; i <= sd->maxshift; i++) {
 267         for (j = -sd->maxshift; j <= sd->maxshift; j++) {
 268             double error = compareImg(sd->curr, sd->prev,
 269                                       sd->width, sd->height, 3, i, j);
 270             if (error < minerror) {
 271                 minerror = error;
 272                 x = i;
 273                 y = j;
 274            }
 275         }
 276     }
 277     return new_transform(x, y, 0, 0, 0);
 278 }
 279
 280
 281 /** tries to register current frame onto previous frame.
 282     (only the luminance is used)
 283     This is the most simple algorithm:
 284     shift images to all possible positions and calc summed error
 285     Shift with minimal error is selected.
 286 */
 287 Transform calcShiftYUVSimple(StabData* sd)
 288 {
 289         printf("calc shoft yuv\n");
 290     int x = 0, y = 0;
 291     int i, j;
 292     unsigned char *Y_c, *Y_p;// , *Cb, *Cr;
 293 #ifdef STABVERBOSE
 294     FILE *f = NULL;
 295     char buffer[32];
 296     tc_snprintf(buffer, sizeof(buffer), "f%04i.dat", sd->t);
 297     f = fopen(buffer, "w");
 298     fprintf(f, "# splot \"%s\"\n", buffer);
 299 #endif
 300
 301     // we only use the luminance part of the image
 302     Y_c  = sd->curr;
 303     //  Cb_c = sd->curr + sd->width*sd->height;
 304     //Cr_c = sd->curr + 5*sd->width*sd->height/4;
 305     Y_p  = sd->prev;
 306     //Cb_p = sd->prev + sd->width*sd->height;
 307     //Cr_p = sd->prev + 5*sd->width*sd->height/4;
 308
 309     double minerror = 1e20;
 310     for (i = -sd->maxshift; i <= sd->maxshift; i++) {
 311         for (j = -sd->maxshift; j <= sd->maxshift; j++) {
 312             double error = compareImg(Y_c, Y_p,
 313                                       sd->width, sd->height, 1, i, j);
 314 #ifdef STABVERBOSE
 315             fprintf(f, "%i %i %f\n", i, j, error);
 316 #endif
 317             if (error < minerror) {
 318                 minerror = error;
 319                 x = i;
 320                 y = j;
 321             }
 322         }
 323     }
 324 #ifdef STABVERBOSE
 325     fclose(f);
 326     tc_log_msg(MOD_NAME, "Minerror: %f\n", minerror);
 327 #endif
 328     return new_transform(x, y, 0, 0, 0);
 329 }
 330
 331
 332
 333 /* calculates rotation angle for the given transform and
 334  * field with respect to the given center-point
 335  */
 336 double calcAngle(StabData* sd, Field* field, Transform* t,
 337                  int center_x, int center_y)
 338 {
 339     // we better ignore fields that are to close to the rotation center
 340     if (abs(field->x - center_x) + abs(field->y - center_y) < sd->maxshift) {
 341         return 0;
 342     } else {
 343         // double r = sqrt(field->x*field->x + field->y*field->y);
 344         double a1 = atan2(field->y - center_y, field->x - center_x);
 345         double a2 = atan2(field->y - center_y + t->y,
 346                           field->x - center_x + t->x);
 347         double diff = a2 - a1;
 348         return (diff>M_PI) ? diff - 2*M_PI
 349             : ( (diff<-M_PI) ? diff + 2*M_PI : diff);
 350     }
 351 }
 352
 353
 354 /* calculates the optimal transformation for one field in YUV frames
 355  * (only luminance)
 356  */
 357 Transform calcFieldTransYUV(StabData* sd, const Field* field, int fieldnum)
 358 {
 359     Transform t = null_transform();
 360     unsigned char *Y_c = sd->curr, *Y_p = sd->prev;
 361     // we only use the luminance part of the image
 362     int i, j;
 363
 364 /*     // check contrast in sub image */
 365 /*     double contr = contrastSubImg(Y_c, field, sd->width, sd->height, 1); */
 366 /*     if(contr < sd->contrast_threshold) { */
 367 /*         t.extra=-1; */
 368 /*         return t; */
 369 /*     } */
 370 #ifdef STABVERBOSE
 371     // printf("%i %i %f\n", sd->t, fieldnum, contr);
 372     FILE *f = NULL;
 373     char buffer[32];
 374     snprintf(buffer, sizeof(buffer), "f%04i_%02i.dat", sd->t, fieldnum);
 375     f = fopen(buffer, "w");
 376     fprintf(f, "# splot \"%s\"\n", buffer);
 377 #endif
 378
 379     double minerror = 1e10;
 380     double error = 1e10;
 381     for (i = -sd->maxshift; i <= sd->maxshift; i += sd->stepsize) {
 382         for (j = -sd->maxshift; j <= sd->maxshift; j += sd->stepsize) {
 383             error = compareSubImg(Y_c, Y_p, field,
 384                                          sd->width, sd->height, 1, i, j);
 385 #ifdef STABVERBOSE
 386             fprintf(f, "%i %i %f\n", i, j, error);
 387 #endif
 388             if (error < minerror) {
 389                 minerror = error;
 390                 t.x = i;
 391                 t.y = j;
 392             }
 393         }
 394     }
 395
 396     if (sd->stepsize > 1) {    // make fine grain check around the best match
 397         int r = sd->stepsize - 1;
 398         for (i = t.x - r; i <= t.x + r; i += 1) {
 399             for (j = -t.y - r; j <= t.y + r; j += 1) {
 400                 if (i == t.x && j == t.y)
 401                     continue; //no need to check this since already done
 402                 error = compareSubImg(Y_c, Y_p, field,
 403                                       sd->width, sd->height, 1, i, j);
 404 #ifdef STABVERBOSE
 405                 fprintf(f, "%i %i %f\n", i, j, error);
 406 #endif
 407                 if (error < minerror){
 408                     minerror = error;
 409                     t.x = i;
 410                     t.y = j;
 411                 }
 412             }
 413         }
 414     }
 415 #ifdef STABVERBOSE
 416     fclose(f);
 417     printf( "Minerror: %f\n", minerror);
 418 #endif
 419
 420     if (!sd->allowmax && fabs(t.x) == sd->maxshift) {
 421 #ifdef STABVERBOSE
 422         printf( "maximal x shift ");
 423 #endif
 424         t.x = 0;
 425     }
 426     if (!sd->allowmax && fabs(t.y) == sd->maxshift) {
 427 #ifdef STABVERBOSE
 428         printf("maximal y shift ");
 429 #endif
 430         t.y = 0;
 431     }
 432     return t;
 433 }
 434
 435 /* calculates the optimal transformation for one field in RGB
 436  *   slower than the YUV version because it uses all three color channels
 437  */
 438 Transform calcFieldTransRGB(StabData* sd, const Field* field, int fieldnum)
 439 {
 440     Transform t = null_transform();
 441     unsigned char *I_c = sd->curr, *I_p = sd->prev;
 442     int i, j;
 443
 444     double minerror = 1e20;
 445     for (i = -sd->maxshift; i <= sd->maxshift; i += 2) {
 446         for (j=-sd->maxshift; j <= sd->maxshift; j += 2) {
 447             double error = compareSubImg(I_c, I_p, field,
 448                                          sd->width, sd->height, 3, i, j);
 449             if (error < minerror) {
 450                 minerror = error;
 451                 t.x = i;
 452                 t.y = j;
 453             }
 454         }
 455     }
 456     for (i = t.x - 1; i <= t.x + 1; i += 2) {
 457         for (j = -t.y - 1; j <= t.y + 1; j += 2) {
 458             double error = compareSubImg(I_c, I_p, field,
 459                                          sd->width, sd->height, 3, i, j);
 460             if (error < minerror) {
 461                 minerror = error;
 462                 t.x = i;
 463                 t.y = j;
 464             }
 465         }
 466     }
 467     if (!sd->allowmax && fabs(t.x) == sd->maxshift) {
 468         t.x = 0;
 469     }
 470     if (!sd->allowmax && fabs(t.y) == sd->maxshift) {
 471         t.y = 0;
 472     }
 473     return t;
 474 }
 475
 476 /* compares contrast_idx structures respect to the contrast
 477    (for sort function)
 478 */
 479 int cmp_contrast_idx(const void *ci1, const void* ci2)
 480 {
 481     double a = ((contrast_idx*)ci1)->contrast;
 482     double b = ((contrast_idx*)ci2)->contrast;
 483     return a < b ? 1 : ( a > b ? -1 : 0 );
 484 }
 485
 486 /* select only the best 'maxfields' fields
 487    first calc contrasts then select from each part of the
 488    frame a some fields
 489 */
 490 tlist* selectfields(StabData* sd, contrastSubImgFunc contrastfunc)
 491 {
 492     int i,j;
 493     tlist* goodflds = tlist_new(0);
 494     contrast_idx *ci = malloc(sizeof(contrast_idx) * sd->field_num);
 495
 496     // we split all fields into row+1 segments and take from each segment
 497     // the best fields
 498     int numsegms = (sd->field_rows+1);
 499     int segmlen = sd->field_num/(sd->field_rows+1)+1;
 500     // split the frame list into rows+1 segments
 501     contrast_idx *ci_segms = malloc(sizeof(contrast_idx) * sd->field_num);
 502     int remaining   = 0;
 503     // calculate contrast for each field
 504     for (i = 0; i < sd->field_num; i++) {
 505         ci[i].contrast = contrastfunc(sd, &sd->fields[i]);
 506         ci[i].index=i;
 507         if(ci[i].contrast < sd->contrast_threshold) ci[i].contrast = 0;
 508         // else printf("%i %lf\n", ci[i].index, ci[i].contrast);
 509     }
 510
 511     memcpy(ci_segms, ci, sizeof(contrast_idx) * sd->field_num);
 512     // get best fields from each segment
 513     for(i=0; i<numsegms; i++){
 514         int startindex = segmlen*i;
 515         int endindex   = segmlen*(i+1);
 516         endindex       = endindex > sd->field_num ? sd->field_num : endindex;
 517         //printf("Segment: %i: %i-%i\n", i, startindex, endindex);
 518
 519         // sort within segment
 520         qsort(ci_segms+startindex, endindex-startindex,
 521               sizeof(contrast_idx), cmp_contrast_idx);
 522         // take maxfields/numsegms
 523         for(j=0; j<sd->maxfields/numsegms; j++){
 524             if(startindex+j >= endindex) continue;
 525             // printf("%i %lf\n", ci_segms[startindex+j].index,
 526             //                    ci_segms[startindex+j].contrast);
 527             if(ci_segms[startindex+j].contrast > 0){
 528                tlist_append(goodflds, &ci[ci_segms[startindex+j].index],sizeof(contrast_idx));
 529                 // don't consider them in the later selection process
 530                 ci_segms[startindex+j].contrast=0;
 531             }
 532         }
 533     }
 534     // check whether enough fields are selected
 535     // printf("Phase2: %i\n", tc_list_size(goodflds));
 536     remaining = sd->maxfields - tlist_size(goodflds);
 537     if(remaining > 0){
 538         // take the remaining from the leftovers
 539         qsort(ci_segms, sd->field_num,
 540               sizeof(contrast_idx), cmp_contrast_idx);
 541         for(j=0; j < remaining; j++){
 542             if(ci_segms[j].contrast > 0){
 543                 tlist_append(goodflds, &ci_segms[j], sizeof(contrast_idx));
 544             }
 545         }
 546     }
 547     // printf("Ende: %i\n", tc_list_size(goodflds));
 548     free(ci);
 549     free(ci_segms);
 550     return goodflds;
 551 }
 552
 553
 554
 555 /* tries to register current frame onto previous frame.
 556  *   Algorithm:
 557  *   check all fields for vertical and horizontal transformation
 558  *   use minimal difference of all possible positions
 559  *   discards fields with low contrast
 560  *   select maxfields field according to their contrast
 561  *   calculate shift as cleaned mean of all remaining fields
 562  *   calculate rotation angle of each field in respect to center of fields
 563  *   after shift removal
 564  *   calculate rotation angle as cleaned mean of all angles
 565  *   compensate for possibly off-center rotation
 566 */
 567 Transform calcTransFields(StabData* sd, calcFieldTransFunc fieldfunc,
 568                           contrastSubImgFunc contrastfunc)
 569 {
 570     Transform* ts  = malloc(sizeof(Transform) * sd->field_num);
 571     Field** fs     = malloc(sizeof(Field*) * sd->field_num);
 572     double *angles = malloc(sizeof(double) * sd->field_num);
 573     int i, index=0, num_trans;
 574     Transform t;
 575 #ifdef STABVERBOSE
 576     FILE *f = NULL;
 577     char buffer[32];
 578     tc_snprintf(buffer, sizeof(buffer), "k%04i.dat", sd->t);
 579     f = fopen(buffer, "w");
 580     fprintf(f, "# plot \"%s\" w l, \"\" every 2:1:0\n", buffer);
 581 #endif
 582
 583
 584     tlist* goodflds = selectfields(sd, contrastfunc);
 585
 586     // use all "good" fields and calculate optimal match to previous frame
 587     contrast_idx* f;
 588     while((f = (contrast_idx*)tlist_pop(goodflds,0) ) != 0){
 589         int i = f->index;
 590         t =  fieldfunc(sd, &sd->fields[i], i); // e.g. calcFieldTransYUV
 591 #ifdef STABVERBOSE
 592         fprintf(f, "%i %i\n%f %f %i\n \n\n", sd->fields[i].x, sd->fields[i].y,
 593                 sd->fields[i].x + t.x, sd->fields[i].y + t.y, t.extra);
 594 #endif
 595         if (t.extra != -1){ // ignore if extra == -1 (unused at the moment)
 596             ts[index] = t;
 597             fs[index] = sd->fields+i;
 598             index++;
 599         }
 600     }
 601     tlist_fini(goodflds);
 602
 603     t = null_transform();
 604     num_trans = index; // amount of transforms we actually have
 605     if (num_trans < 1) {
 606         printf( "too low contrast! No field remains.\n \
 607                     (no translations are detected in frame %i)", sd->t);
 608         return t;
 609     }
 610
 611     int center_x = 0;
 612     int center_y = 0;
 613     // calc center point of all remaining fields
 614     for (i = 0; i < num_trans; i++) {
 615         center_x += fs[i]->x;
 616         center_y += fs[i]->y;
 617     }
 618     center_x /= num_trans;
 619     center_y /= num_trans;
 620
 621     if (sd->show){ // draw fields and transforms into frame.
 622         // this has to be done one after another to handle possible overlap
 623         if (sd->show > 1) {
 624             for (i = 0; i < num_trans; i++)
 625                 drawFieldScanArea(sd, fs[i], &ts[i]);
 626         }
 627         for (i = 0; i < num_trans; i++)
 628             drawField(sd, fs[i], &ts[i]);
 629         for (i = 0; i < num_trans; i++)
 630             drawFieldTrans(sd, fs[i], &ts[i]);
 631     }
 632     /* median over all transforms
 633        t= median_xy_transform(ts, sd->field_num);*/
 634     // cleaned mean
 635     t = cleanmean_xy_transform(ts, num_trans);
 636
 637     // substract avg
 638     for (i = 0; i < num_trans; i++) {
 639         ts[i] = sub_transforms(&ts[i], &t);
 640     }
 641     // figure out angle
 642     if (sd->field_num < 6) {
 643         // the angle calculation is inaccurate for 5 and less fields
 644         t.alpha = 0;
 645     } else {
 646         for (i = 0; i < num_trans; i++) {
 647             angles[i] = calcAngle(sd, fs[i], &ts[i], center_x, center_y);
 648         }
 649         double min,max;
 650         t.alpha = -cleanmean(angles, num_trans, &min, &max);
 651         if(max-min>sd->maxanglevariation){
 652             t.alpha=0;
 653             printf( "too large variation in angle(%f)\n",
 654                         max-min);
 655         }
 656     }
 657
 658     // compensate for off-center rotation
 659     double p_x = (center_x - sd->width/2);
 660     double p_y = (center_y - sd->height/2);
 661     t.x += (cos(t.alpha)-1)*p_x  - sin(t.alpha)*p_y;
 662     t.y += sin(t.alpha)*p_x  + (cos(t.alpha)-1)*p_y;
 663
 664 #ifdef STABVERBOSE
 665     fclose(f);
 666 #endif
 667     return t;
 668 }
 669
 670 /** draws the field scanning area */
 671 void drawFieldScanArea(StabData* sd, const Field* field, const Transform* t)
 672 {
 673     if (!sd->pixelformat == mlt_image_yuv420p) {
 674                 printf("kein format\n");
 675         return;
 676         }
 677     drawBox(sd->curr, sd->width, sd->height, 1, field->x, field->y,
 678             field->size+2*sd->maxshift, field->size+2*sd->maxshift, 80);
 679 }
 680
 681 /** draws the field */
 682 void drawField(StabData* sd, const Field* field, const Transform* t)
 683 {
 684     if (!sd->pixelformat == mlt_image_yuv420p){
 685                 printf("kein format\n");
 686         return;
 687         }
 688     drawBox(sd->curr, sd->width, sd->height, 1, field->x, field->y,
 689             field->size, field->size, t->extra == -1 ? 100 : 40);
 690 }
 691
 692 /** draws the transform data of this field */
 693 void drawFieldTrans(StabData* sd, const Field* field, const Transform* t)
 694 {
 695     if (!sd->pixelformat == mlt_image_yuv420p){
 696                 printf("kein format\n");
 697         return;
 698         }
 699     drawBox(sd->curr, sd->width, sd->height, 1,
 700             field->x, field->y, 5, 5, 128);     // draw center
 701     drawBox(sd->curr, sd->width, sd->height, 1,
 702             field->x + t->x, field->y + t->y, 8, 8, 250); // draw translation
 703 }
 704
 705 /**
 706  * draws a box at the given position x,y (center) in the given color
 707    (the same for all channels)
 708  */
 709 void drawBox(unsigned char* I, int width, int height, int bytesPerPixel,
 710              int x, int y, int sizex, int sizey, unsigned char color){
 711
 712     unsigned char* p = NULL;
 713     int j,k;
 714     p = I + ((x - sizex/2) + (y - sizey/2)*width)*bytesPerPixel;
 715     for (j = 0; j < sizey; j++){
 716         for (k = 0; k < sizex * bytesPerPixel; k++) {
 717             *p = color;
 718             p++;
 719         }
 720         p += (width - sizex) * bytesPerPixel;
 721     }
 722 }
 723
 724 struct iterdata {
 725     FILE *f;
 726     int  counter;
 727 };
 728 #if 0
 729 static int stabilize_dump_trans(tlist*item, void *userdata)
 730 {
 731     struct iterdata *ID = userdata;
 732
 733     if (item->data) {
 734         Transform* t = (Transform*) item->data;
 735         fprintf(ID->f, "%i %6.4lf %6.4lf %8.5lf %6.4lf %i\n",
 736                 ID->counter, t->x, t->y, t->alpha, t->zoom, t->extra);
 737         ID->counter++;
 738     }
 739     return 0; /* never give up */
 740 }
 741
 742 #endif
 743 /*************************************************************************/
 744
 745 /* Module interface routines and data. */
 746
 747 /*************************************************************************/
 748
 749 /**
 750  * stabilize_init:  Initialize this instance of the module.  See
 751  * tcmodule-data.h for function details.
 752  */
 753 int stabilize_init(StabData* instance)
 754 {
 755
 756     instance = malloc(sizeof(StabData)); // allocation with zero values
 757         memset(instance,sizeof(StabData),0);
 758     if (!instance) {
 759         return -1;
 760     }
 761
 762         //stab->vob=pixelfmt;
 763     /**** Initialise private data structure */
 764
 765     //self->userdata = sd;
 766
 767     return 0;
 768 }
 769
 770 /*
 771  * stabilize_configure:  Configure this instance of the module.  See
 772  * tcmodule-data.h for function details.
 773  */
 774 int stabilize_configure(StabData* instance
 775 /*const char *options,
 776                                int pixelfmt */
 777                                /*TCModuleExtraData *xdata[]*/)
 778 {
 779     StabData *sd = instance;
 780     /*    sd->framesize = sd->vob->im_v_width * MAX_PLANES *
 781           sizeof(char) * 2 * sd->vob->im_v_height * 2;     */
 782     /*TODO sd->framesize = sd->vob->im_v_size;    */
 783     sd->prev = calloc(1,sd->framesize);
 784     if (!sd->prev) {
 785         printf( "malloc failed");
 786         return -1;
 787     }
 788     sd->currcopy = 0;
 789
 790     /*sd->width  = sd->vob->ex_v_width;
 791     sd->height = sd->vob->ex_v_height;
 792 */
 793     sd->hasSeenOneFrame = 0;
 794     sd->transs = 0;
 795
 796     // Options
 797     // done in filter : sd->stepsize   = 6;
 798     sd->allowmax   = 0;
 799     // done in filter :sd->algo = 1;
 800 //    sd->field_num   = 64;
 801     // done in filter : sd->accuracy    = 4;
 802     // done in filter : sd->shakiness   = 4;
 803     sd->field_size  = MIN(sd->width, sd->height)/12;
 804     // done in filter : sd->show        = 0;
 805     // done in filter : sd->contrast_threshold = 0.3;
 806     sd->maxanglevariation = 1;
 807
 808     /*if (options != NULL) {
 809         // for some reason this plugin is called in the old fashion
 810         //  (not with inspect). Anyway we support both ways of getting help.
 811         if(optstr_lookup(options, "help")) {
 812             printf(stabilize_help);
 813             return -2;
 814         }
 815
 816         optstr_get(options, "shakiness",  "%d", &sd->shakiness);
 817         optstr_get(options, "accuracy",   "%d", &sd->accuracy);
 818         optstr_get(options, "stepsize",   "%d", &sd->stepsize);
 819         optstr_get(options, "algo",       "%d", &sd->algo);
 820         optstr_get(options, "mincontrast","%lf",&sd->contrast_threshold);
 821         optstr_get(options, "show",       "%d", &sd->show);
 822     }*/
 823     sd->shakiness = MIN(10,MAX(1,sd->shakiness));
 824     sd->accuracy  = MAX(sd->shakiness,MIN(15,MAX(1,sd->accuracy)));
 825     if (1) {
 826         printf( "Image Stabilization Settings:\n");
 827         printf( "     shakiness = %d\n", sd->shakiness);
 828         printf( "      accuracy = %d\n", sd->accuracy);
 829         printf( "      stepsize = %d\n", sd->stepsize);
 830         printf( "          algo = %d\n", sd->algo);
 831         printf("   mincontrast = %f\n", sd->contrast_threshold);
 832         printf( "          show = %d\n", sd->show);
 833     }
 834
 835     // shift and size: shakiness 1: height/40; 10: height/4
 836     sd->maxshift    = MIN(sd->width, sd->height)*sd->shakiness/40;
 837     sd->field_size   = MIN(sd->width, sd->height)*sd->shakiness/40;
 838
 839     printf( "Fieldsize: %i, Maximal translation: %i pixel\n",
 840                 sd->field_size, sd->maxshift);
 841     if (sd->algo==1) {
 842         // initialize measurement fields. field_num is set here.
 843         if (!initFields(sd)) {
 844             return -1;
 845         }
 846         sd->maxfields = (sd->accuracy) * sd->field_num / 15;
 847         printf( "Number of used measurement fields: %i out of %i\n",
 848                     sd->maxfields, sd->field_num);
 849     }
 850     if (sd->show){
 851         sd->currcopy = malloc(sd->framesize);
 852                 memset ( sd->currcopy, sd->framesize, 0 );
 853         }
 854
 855     /* load unsharp filter to smooth the frames. This allows larger stepsize.*/
 856     char unsharp_param[128];
 857     int masksize = MIN(13,sd->stepsize*1.8); // only works up to 13.
 858     sprintf(unsharp_param,"luma=-1:luma_matrix=%ix%i:pre=1",
 859             masksize, masksize);
 860     return 0;
 861 }
 862
 863
 864 /**
 865  * stabilize_filter_video: performs the analysis of subsequent frames
 866  * See tcmodule-data.h for function details.
 867  */
 868
 869 int stabilize_filter_video(StabData* instance,
 870                                   unsigned char *frame,mlt_image_format pixelformat)
 871 {
 872     StabData *sd = instance;
 873         sd->pixelformat=pixelformat;
 874
 875     if(sd->show)  // save the buffer to restore at the end for prev
 876         memcpy(sd->currcopy, frame, sd->framesize);
 877     if (sd->hasSeenOneFrame) {
 878         sd->curr = frame;
 879         if (pixelformat == mlt_image_rgb24) {
 880             if (sd->algo == 0)
 881                 addTrans(sd, calcShiftRGBSimple(sd));
 882             else if (sd->algo == 1)
 883                 addTrans(sd, calcTransFields(sd, calcFieldTransRGB,
 884                          contrastSubImgRGB));
 885         } else if (pixelformat == mlt_image_yuv420p ) {
 886             if (sd->algo == 0)
 887                 addTrans(sd, calcShiftYUVSimple(sd));
 888             else if (sd->algo == 1)
 889                 addTrans(sd, calcTransFields(sd, calcFieldTransYUV,
 890                                              contrastSubImgYUV));
 891         } else {
 892             printf("unsupported Codec: %i\n",
 893                         pixelformat);
 894             return 0;
 895         }
 896     } else {
 897         sd->hasSeenOneFrame = 1;
 898         addTrans(sd, null_transform());
 899     }
 900
 901     if(!sd->show) { // copy current frame to prev for next frame comparison
 902         memcpy(sd->prev, frame, sd->framesize);
 903     } else { // use the copy because we changed the original frame
 904         memcpy(sd->prev, sd->currcopy, sd->framesize);
 905     }
 906     sd->t++;
 907     return 0;
 908 }
 909
 910 /**
 911  * stabilize_stop:  Reset this instance of the module.  See tcmodule-data.h
 912  * for function details.
 913  */
 914
 915 int stabilize_stop(StabData* instance)
 916 {
 917     StabData *sd = instance;
 918
 919     // print transs
 920         /*struct iterdata ID;
 921         ID.counter = 0;
 922         ID.f       = sd->f;
 923         // write parameters as comments to file
 924         fprintf(sd->f, "#      accuracy = %d\n", sd->accuracy);
 925         fprintf(sd->f, "#     shakiness = %d\n", sd->shakiness);
 926         fprintf(sd->f, "#      stepsize = %d\n", sd->stepsize);
 927         fprintf(sd->f, "#          algo = %d\n", sd->algo);
 928         fprintf(sd->f, "#   mincontrast = %f\n", sd->contrast_threshold);
 929         // write header line
 930         fprintf(sd->f, "# Transforms\n#C FrameNr x y alpha zoom extra\n");
 931         // and all transforms
 932         tc_list_foreach(sd->transs, stabilize_dump_trans, &ID);
 933 */
 934         tlist* transf=sd->transs;
 935         int num=0;
 936         while (transf ){
 937                 Transform* t=transf->data;
 938                 if (t)
 939                 printf("%d %f %f %f %f %d\n",num++,t->x,t->y,t->alpha,t->zoom,t->extra);
 940                 transf=transf->next;
 941         }
 942     //tc_list_del(sd->transs, 1 );
 943     if (sd->prev) {
 944         free(sd->prev);
 945         sd->prev = NULL;
 946     }
 947     return 0;
 948 }
 949
 950
 951