git.sesse.net Git - ffmpeg/blob - libavcodec/imgresample.c

   1 /*
   2  * High quality image resampling with polyphase filters
   3  * Copyright (c) 2001 Gerard Lantau.
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 2 of the License, or
   8  * (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18  */
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <string.h>
  22 #include <math.h>
  23 #include "dsputil.h"
  24 #include "avcodec.h"
  25
  26 #ifdef USE_FASTMEMCPY
  27 #include "fastmemcpy.h"
  28 #endif
  29
  30
  31 #define NB_COMPONENTS 3
  32
  33 #define PHASE_BITS 4
  34 #define NB_PHASES  (1 << PHASE_BITS)
  35 #define NB_TAPS    4
  36 #define FCENTER    1  /* index of the center of the filter */
  37
  38 #define POS_FRAC_BITS 16
  39 #define POS_FRAC      (1 << POS_FRAC_BITS)
  40 /* 6 bits precision is needed for MMX */
  41 #define FILTER_BITS   8
  42
  43 #define LINE_BUF_HEIGHT (NB_TAPS * 4)
  44
  45 struct ImgReSampleContext {
  46     int iwidth, iheight, owidth, oheight;
  47     int h_incr, v_incr;
  48     INT16 h_filters[NB_PHASES][NB_TAPS] __align8; /* horizontal filters */
  49     INT16 v_filters[NB_PHASES][NB_TAPS] __align8; /* vertical filters */
  50     UINT8 *line_buf;
  51 };
  52
  53 static inline int get_phase(int pos)
  54 {
  55     return ((pos) >> (POS_FRAC_BITS - PHASE_BITS)) & ((1 << PHASE_BITS) - 1);
  56 }
  57
  58 /* This function must be optimized */
  59 static void h_resample_fast(UINT8 *dst, int dst_width, UINT8 *src, int src_width,
  60                             int src_start, int src_incr, INT16 *filters)
  61 {
  62     int src_pos, phase, sum, i;
  63     UINT8 *s;
  64     INT16 *filter;
  65
  66     src_pos = src_start;
  67     for(i=0;i<dst_width;i++) {
  68 #ifdef TEST
  69         /* test */
  70         if ((src_pos >> POS_FRAC_BITS) < 0 ||
  71             (src_pos >> POS_FRAC_BITS) > (src_width - NB_TAPS))
  72             abort();
  73 #endif
  74         s = src + (src_pos >> POS_FRAC_BITS);
  75         phase = get_phase(src_pos);
  76         filter = filters + phase * NB_TAPS;
  77 #if NB_TAPS == 4
  78         sum = s[0] * filter[0] +
  79             s[1] * filter[1] +
  80             s[2] * filter[2] +
  81             s[3] * filter[3];
  82 #else
  83         {
  84             int j;
  85             sum = 0;
  86             for(j=0;j<NB_TAPS;j++)
  87                 sum += s[j] * filter[j];
  88         }
  89 #endif
  90         sum = sum >> FILTER_BITS;
  91         if (sum < 0)
  92             sum = 0;
  93         else if (sum > 255)
  94             sum = 255;
  95         dst[0] = sum;
  96         src_pos += src_incr;
  97         dst++;
  98     }
  99 }
 100
 101 /* This function must be optimized */
 102 static void v_resample(UINT8 *dst, int dst_width, UINT8 *src, int wrap,
 103                        INT16 *filter)
 104 {
 105     int sum, i;
 106     UINT8 *s;
 107
 108     s = src;
 109     for(i=0;i<dst_width;i++) {
 110 #if NB_TAPS == 4
 111         sum = s[0 * wrap] * filter[0] +
 112             s[1 * wrap] * filter[1] +
 113             s[2 * wrap] * filter[2] +
 114             s[3 * wrap] * filter[3];
 115 #else
 116         {
 117             int j;
 118             UINT8 *s1 = s;
 119
 120             sum = 0;
 121             for(j=0;j<NB_TAPS;j++) {
 122                 sum += s1[0] * filter[j];
 123                 s1 += wrap;
 124             }
 125         }
 126 #endif
 127         sum = sum >> FILTER_BITS;
 128         if (sum < 0)
 129             sum = 0;
 130         else if (sum > 255)
 131             sum = 255;
 132         dst[0] = sum;
 133         dst++;
 134         s++;
 135     }
 136 }
 137
 138 #ifdef HAVE_MMX
 139
 140 #include "i386/mmx.h"
 141
 142 #define FILTER4(reg) \
 143 {\
 144         s = src + (src_pos >> POS_FRAC_BITS);\
 145         phase = get_phase(src_pos);\
 146         filter = filters + phase * NB_TAPS;\
 147         movq_m2r(*s, reg);\
 148         punpcklbw_r2r(mm7, reg);\
 149         movq_m2r(*filter, mm6);\
 150         pmaddwd_r2r(reg, mm6);\
 151         movq_r2r(mm6, reg);\
 152         psrlq_i2r(32, reg);\
 153         paddd_r2r(mm6, reg);\
 154         psrad_i2r(FILTER_BITS, reg);\
 155         src_pos += src_incr;\
 156 }
 157
 158 #define DUMP(reg) movq_r2m(reg, tmp); printf(#reg "=%016Lx\n", tmp.uq);
 159
 160 /* XXX: do four pixels at a time */
 161 static void h_resample_fast4_mmx(UINT8 *dst, int dst_width, UINT8 *src, int src_width,
 162                                  int src_start, int src_incr, INT16 *filters)
 163 {
 164     int src_pos, phase;
 165     UINT8 *s;
 166     INT16 *filter;
 167     mmx_t tmp;
 168
 169     src_pos = src_start;
 170     pxor_r2r(mm7, mm7);
 171
 172     while (dst_width >= 4) {
 173
 174         FILTER4(mm0);
 175         FILTER4(mm1);
 176         FILTER4(mm2);
 177         FILTER4(mm3);
 178
 179         packuswb_r2r(mm7, mm0);
 180         packuswb_r2r(mm7, mm1);
 181         packuswb_r2r(mm7, mm3);
 182         packuswb_r2r(mm7, mm2);
 183         movq_r2m(mm0, tmp);
 184         dst[0] = tmp.ub[0];
 185         movq_r2m(mm1, tmp);
 186         dst[1] = tmp.ub[0];
 187         movq_r2m(mm2, tmp);
 188         dst[2] = tmp.ub[0];
 189         movq_r2m(mm3, tmp);
 190         dst[3] = tmp.ub[0];
 191         dst += 4;
 192         dst_width -= 4;
 193     }
 194     while (dst_width > 0) {
 195         FILTER4(mm0);
 196         packuswb_r2r(mm7, mm0);
 197         movq_r2m(mm0, tmp);
 198         dst[0] = tmp.ub[0];
 199         dst++;
 200         dst_width--;
 201     }
 202     emms();
 203 }
 204
 205 static void v_resample4_mmx(UINT8 *dst, int dst_width, UINT8 *src, int wrap,
 206                             INT16 *filter)
 207 {
 208     int sum, i, v;
 209     UINT8 *s;
 210     mmx_t tmp;
 211     mmx_t coefs[4];
 212
 213     for(i=0;i<4;i++) {
 214         v = filter[i];
 215         coefs[i].uw[0] = v;
 216         coefs[i].uw[1] = v;
 217         coefs[i].uw[2] = v;
 218         coefs[i].uw[3] = v;
 219     }
 220
 221     pxor_r2r(mm7, mm7);
 222     s = src;
 223     while (dst_width >= 4) {
 224         movq_m2r(s[0 * wrap], mm0);
 225         punpcklbw_r2r(mm7, mm0);
 226         movq_m2r(s[1 * wrap], mm1);
 227         punpcklbw_r2r(mm7, mm1);
 228         movq_m2r(s[2 * wrap], mm2);
 229         punpcklbw_r2r(mm7, mm2);
 230         movq_m2r(s[3 * wrap], mm3);
 231         punpcklbw_r2r(mm7, mm3);
 232
 233         pmullw_m2r(coefs[0], mm0);
 234         pmullw_m2r(coefs[1], mm1);
 235         pmullw_m2r(coefs[2], mm2);
 236         pmullw_m2r(coefs[3], mm3);
 237
 238         paddw_r2r(mm1, mm0);
 239         paddw_r2r(mm3, mm2);
 240         paddw_r2r(mm2, mm0);
 241         psraw_i2r(FILTER_BITS, mm0);
 242
 243         packuswb_r2r(mm7, mm0);
 244         movq_r2m(mm0, tmp);
 245
 246         *(UINT32 *)dst = tmp.ud[0];
 247         dst += 4;
 248         s += 4;
 249         dst_width -= 4;
 250     }
 251     while (dst_width > 0) {
 252         sum = s[0 * wrap] * filter[0] +
 253             s[1 * wrap] * filter[1] +
 254             s[2 * wrap] * filter[2] +
 255             s[3 * wrap] * filter[3];
 256         sum = sum >> FILTER_BITS;
 257         if (sum < 0)
 258             sum = 0;
 259         else if (sum > 255)
 260             sum = 255;
 261         dst[0] = sum;
 262         dst++;
 263         s++;
 264         dst_width--;
 265     }
 266     emms();
 267 }
 268 #endif
 269
 270 /* slow version to handle limit cases. Does not need optimisation */
 271 static void h_resample_slow(UINT8 *dst, int dst_width, UINT8 *src, int src_width,
 272                             int src_start, int src_incr, INT16 *filters)
 273 {
 274     int src_pos, phase, sum, j, v, i;
 275     UINT8 *s, *src_end;
 276     INT16 *filter;
 277
 278     src_end = src + src_width;
 279     src_pos = src_start;
 280     for(i=0;i<dst_width;i++) {
 281         s = src + (src_pos >> POS_FRAC_BITS);
 282         phase = get_phase(src_pos);
 283         filter = filters + phase * NB_TAPS;
 284         sum = 0;
 285         for(j=0;j<NB_TAPS;j++) {
 286             if (s < src)
 287                 v = src[0];
 288             else if (s >= src_end)
 289                 v = src_end[-1];
 290             else
 291                 v = s[0];
 292             sum += v * filter[j];
 293             s++;
 294         }
 295         sum = sum >> FILTER_BITS;
 296         if (sum < 0)
 297             sum = 0;
 298         else if (sum > 255)
 299             sum = 255;
 300         dst[0] = sum;
 301         src_pos += src_incr;
 302         dst++;
 303     }
 304 }
 305
 306 static void h_resample(UINT8 *dst, int dst_width, UINT8 *src, int src_width,
 307                        int src_start, int src_incr, INT16 *filters)
 308 {
 309     int n, src_end;
 310
 311     if (src_start < 0) {
 312         n = (0 - src_start + src_incr - 1) / src_incr;
 313         h_resample_slow(dst, n, src, src_width, src_start, src_incr, filters);
 314         dst += n;
 315         dst_width -= n;
 316         src_start += n * src_incr;
 317     }
 318     src_end = src_start + dst_width * src_incr;
 319     if (src_end > ((src_width - NB_TAPS) << POS_FRAC_BITS)) {
 320         n = (((src_width - NB_TAPS + 1) << POS_FRAC_BITS) - 1 - src_start) /
 321             src_incr;
 322     } else {
 323         n = dst_width;
 324     }
 325 #ifdef HAVE_MMX
 326     if ((mm_flags & MM_MMX) && NB_TAPS == 4)
 327         h_resample_fast4_mmx(dst, n,
 328                              src, src_width, src_start, src_incr, filters);
 329     else
 330 #endif
 331         h_resample_fast(dst, n,
 332                         src, src_width, src_start, src_incr, filters);
 333     if (n < dst_width) {
 334         dst += n;
 335         dst_width -= n;
 336         src_start += n * src_incr;
 337         h_resample_slow(dst, dst_width,
 338                         src, src_width, src_start, src_incr, filters);
 339     }
 340 }
 341
 342 static void component_resample(ImgReSampleContext *s,
 343                                UINT8 *output, int owrap, int owidth, int oheight,
 344                                UINT8 *input, int iwrap, int iwidth, int iheight)
 345 {
 346     int src_y, src_y1, last_src_y, ring_y, phase_y, y1, y;
 347     UINT8 *new_line, *src_line;
 348
 349     last_src_y = - FCENTER - 1;
 350     /* position of the bottom of the filter in the source image */
 351     src_y = (last_src_y + NB_TAPS) * POS_FRAC;
 352     ring_y = NB_TAPS; /* position in ring buffer */
 353     for(y=0;y<oheight;y++) {
 354         /* apply horizontal filter on new lines from input if needed */
 355         src_y1 = src_y >> POS_FRAC_BITS;
 356         while (last_src_y < src_y1) {
 357             if (++ring_y >= LINE_BUF_HEIGHT + NB_TAPS)
 358                 ring_y = NB_TAPS;
 359             last_src_y++;
 360             /* handle limit conditions : replicate line (slighly
 361                inefficient because we filter multiple times */
 362             y1 = last_src_y;
 363             if (y1 < 0) {
 364                 y1 = 0;
 365             } else if (y1 >= iheight) {
 366                 y1 = iheight - 1;
 367             }
 368             src_line = input + y1 * iwrap;
 369             new_line = s->line_buf + ring_y * owidth;
 370             /* apply filter and handle limit cases correctly */
 371             h_resample(new_line, owidth,
 372                        src_line, iwidth, - FCENTER * POS_FRAC, s->h_incr,
 373                        &s->h_filters[0][0]);
 374             /* handle ring buffer wraping */
 375             if (ring_y >= LINE_BUF_HEIGHT) {
 376                 memcpy(s->line_buf + (ring_y - LINE_BUF_HEIGHT) * owidth,
 377                        new_line, owidth);
 378             }
 379         }
 380         /* apply vertical filter */
 381         phase_y = get_phase(src_y);
 382 #ifdef HAVE_MMX
 383         /* desactivated MMX because loss of precision */
 384         if ((mm_flags & MM_MMX) && NB_TAPS == 4 && 0)
 385             v_resample4_mmx(output, owidth,
 386                             s->line_buf + (ring_y - NB_TAPS + 1) * owidth, owidth,
 387                             &s->v_filters[phase_y][0]);
 388         else
 389 #endif
 390             v_resample(output, owidth,
 391                        s->line_buf + (ring_y - NB_TAPS + 1) * owidth, owidth,
 392                        &s->v_filters[phase_y][0]);
 393
 394         src_y += s->v_incr;
 395         output += owrap;
 396     }
 397 }
 398
 399 /* XXX: the following filter is quite naive, but it seems to suffice
 400    for 4 taps */
 401 static void build_filter(INT16 *filter, float factor)
 402 {
 403     int ph, i, v;
 404     float x, y, tab[NB_TAPS], norm, mult;
 405
 406     /* if upsampling, only need to interpolate, no filter */
 407     if (factor > 1.0)
 408         factor = 1.0;
 409
 410     for(ph=0;ph<NB_PHASES;ph++) {
 411         norm = 0;
 412         for(i=0;i<NB_TAPS;i++) {
 413
 414             x = M_PI * ((float)(i - FCENTER) - (float)ph / NB_PHASES) * factor;
 415             if (x == 0)
 416                 y = 1.0;
 417             else
 418                 y = sin(x) / x;
 419             tab[i] = y;
 420             norm += y;
 421         }
 422
 423         /* normalize so that an uniform color remains the same */
 424         mult = (float)(1 << FILTER_BITS) / norm;
 425         for(i=0;i<NB_TAPS;i++) {
 426             v = (int)(tab[i] * mult);
 427             filter[ph * NB_TAPS + i] = v;
 428         }
 429     }
 430 }
 431
 432 ImgReSampleContext *img_resample_init(int owidth, int oheight,
 433                                       int iwidth, int iheight)
 434 {
 435     ImgReSampleContext *s;
 436
 437     s = av_mallocz(sizeof(ImgReSampleContext));
 438     if (!s)
 439         return NULL;
 440     s->line_buf = av_mallocz(owidth * (LINE_BUF_HEIGHT + NB_TAPS));
 441     if (!s->line_buf)
 442         goto fail;
 443
 444     s->owidth = owidth;
 445     s->oheight = oheight;
 446     s->iwidth = iwidth;
 447     s->iheight = iheight;
 448
 449     s->h_incr = (iwidth * POS_FRAC) / owidth;
 450     s->v_incr = (iheight * POS_FRAC) / oheight;
 451
 452     build_filter(&s->h_filters[0][0], (float)owidth / (float)iwidth);
 453     build_filter(&s->v_filters[0][0], (float)oheight / (float)iheight);
 454
 455     return s;
 456  fail:
 457     free(s);
 458     return NULL;
 459 }
 460
 461 void img_resample(ImgReSampleContext *s,
 462                   AVPicture *output, AVPicture *input)
 463 {
 464     int i, shift;
 465
 466     for(i=0;i<3;i++) {
 467         shift = (i == 0) ? 0 : 1;
 468         component_resample(s, output->data[i], output->linesize[i],
 469                            s->owidth >> shift, s->oheight >> shift,
 470                            input->data[i], input->linesize[i],
 471                            s->iwidth >> shift, s->iheight >> shift);
 472     }
 473 }
 474
 475 void img_resample_close(ImgReSampleContext *s)
 476 {
 477     free(s->line_buf);
 478     free(s);
 479 }
 480
 481 #ifdef TEST
 482
 483 void *av_mallocz(int size)
 484 {
 485     void *ptr;
 486     ptr = malloc(size);
 487     memset(ptr, 0, size);
 488     return ptr;
 489 }
 490
 491 /* input */
 492 #define XSIZE 256
 493 #define YSIZE 256
 494 UINT8 img[XSIZE * YSIZE];
 495
 496 /* output */
 497 #define XSIZE1 512
 498 #define YSIZE1 512
 499 UINT8 img1[XSIZE1 * YSIZE1];
 500 UINT8 img2[XSIZE1 * YSIZE1];
 501
 502 void save_pgm(const char *filename, UINT8 *img, int xsize, int ysize)
 503 {
 504     FILE *f;
 505     f=fopen(filename,"w");
 506     fprintf(f,"P5\n%d %d\n%d\n", xsize, ysize, 255);
 507     fwrite(img,1, xsize * ysize,f);
 508     fclose(f);
 509 }
 510
 511 static void dump_filter(INT16 *filter)
 512 {
 513     int i, ph;
 514
 515     for(ph=0;ph<NB_PHASES;ph++) {
 516         printf("%2d: ", ph);
 517         for(i=0;i<NB_TAPS;i++) {
 518             printf(" %5.2f", filter[ph * NB_TAPS + i] / 256.0);
 519         }
 520         printf("\n");
 521     }
 522 }
 523
 524 #ifdef HAVE_MMX
 525 int mm_flags;
 526 #endif
 527
 528 int main(int argc, char **argv)
 529 {
 530     int x, y, v, i, xsize, ysize;
 531     ImgReSampleContext *s;
 532     float fact, factors[] = { 1/2.0, 3.0/4.0, 1.0, 4.0/3.0, 16.0/9.0, 2.0 };
 533     char buf[256];
 534
 535     /* build test image */
 536     for(y=0;y<YSIZE;y++) {
 537         for(x=0;x<XSIZE;x++) {
 538             if (x < XSIZE/2 && y < YSIZE/2) {
 539                 if (x < XSIZE/4 && y < YSIZE/4) {
 540                     if ((x % 10) <= 6 &&
 541                         (y % 10) <= 6)
 542                         v = 0xff;
 543                     else
 544                         v = 0x00;
 545                 } else if (x < XSIZE/4) {
 546                     if (x & 1)
 547                         v = 0xff;
 548                     else
 549                         v = 0;
 550                 } else if (y < XSIZE/4) {
 551                     if (y & 1)
 552                         v = 0xff;
 553                     else
 554                         v = 0;
 555                 } else {
 556                     if (y < YSIZE*3/8) {
 557                         if ((y+x) & 1)
 558                             v = 0xff;
 559                         else
 560                             v = 0;
 561                     } else {
 562                         if (((x+3) % 4) <= 1 &&
 563                             ((y+3) % 4) <= 1)
 564                             v = 0xff;
 565                         else
 566                             v = 0x00;
 567                     }
 568                 }
 569             } else if (x < XSIZE/2) {
 570                 v = ((x - (XSIZE/2)) * 255) / (XSIZE/2);
 571             } else if (y < XSIZE/2) {
 572                 v = ((y - (XSIZE/2)) * 255) / (XSIZE/2);
 573             } else {
 574                 v = ((x + y - XSIZE) * 255) / XSIZE;
 575             }
 576             img[y * XSIZE + x] = v;
 577         }
 578     }
 579     save_pgm("/tmp/in.pgm", img, XSIZE, YSIZE);
 580     for(i=0;i<sizeof(factors)/sizeof(float);i++) {
 581         fact = factors[i];
 582         xsize = (int)(XSIZE * fact);
 583         ysize = (int)(YSIZE * fact);
 584         s = img_resample_init(xsize, ysize, XSIZE, YSIZE);
 585         printf("Factor=%0.2f\n", fact);
 586         dump_filter(&s->h_filters[0][0]);
 587         component_resample(s, img1, xsize, xsize, ysize,
 588                            img, XSIZE, XSIZE, YSIZE);
 589         img_resample_close(s);
 590
 591         sprintf(buf, "/tmp/out%d.pgm", i);
 592         save_pgm(buf, img1, xsize, ysize);
 593     }
 594
 595     /* mmx test */
 596 #ifdef HAVE_MMX
 597     printf("MMX test\n");
 598     fact = 0.72;
 599     xsize = (int)(XSIZE * fact);
 600     ysize = (int)(YSIZE * fact);
 601     mm_flags = MM_MMX;
 602     s = img_resample_init(xsize, ysize, XSIZE, YSIZE);
 603     component_resample(s, img1, xsize, xsize, ysize,
 604                        img, XSIZE, XSIZE, YSIZE);
 605
 606     mm_flags = 0;
 607     s = img_resample_init(xsize, ysize, XSIZE, YSIZE);
 608     component_resample(s, img2, xsize, xsize, ysize,
 609                        img, XSIZE, XSIZE, YSIZE);
 610     if (memcmp(img1, img2, xsize * ysize) != 0) {
 611         fprintf(stderr, "mmx error\n");
 612         exit(1);
 613     }
 614     printf("MMX OK\n");
 615 #endif
 616     return 0;
 617 }
 618
 619 #endif