2 * High quality image resampling with polyphase filters
3 * Copyright (c) 2001 Fabrice Bellard.
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2 of the License, or (at your option) any later version.
10 * This library 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 GNU
13 * Lesser General Public License for more details.
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with this library; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 #include "fastmemcpy.h"
27 #define NB_COMPONENTS 3
30 #define NB_PHASES (1 << PHASE_BITS)
32 #define FCENTER 1 /* index of the center of the filter */
33 //#define TEST 1 /* Test it */
35 #define POS_FRAC_BITS 16
36 #define POS_FRAC (1 << POS_FRAC_BITS)
37 /* 6 bits precision is needed for MMX */
40 #define LINE_BUF_HEIGHT (NB_TAPS * 4)
42 struct ImgReSampleContext {
43 int iwidth, iheight, owidth, oheight, topBand, bottomBand, leftBand, rightBand;
45 INT16 h_filters[NB_PHASES][NB_TAPS] __align8; /* horizontal filters */
46 INT16 v_filters[NB_PHASES][NB_TAPS] __align8; /* vertical filters */
50 static inline int get_phase(int pos)
52 return ((pos) >> (POS_FRAC_BITS - PHASE_BITS)) & ((1 << PHASE_BITS) - 1);
55 /* This function must be optimized */
56 static void h_resample_fast(UINT8 *dst, int dst_width, UINT8 *src, int src_width,
57 int src_start, int src_incr, INT16 *filters)
59 int src_pos, phase, sum, i;
64 for(i=0;i<dst_width;i++) {
67 if ((src_pos >> POS_FRAC_BITS) < 0 ||
68 (src_pos >> POS_FRAC_BITS) > (src_width - NB_TAPS))
71 s = src + (src_pos >> POS_FRAC_BITS);
72 phase = get_phase(src_pos);
73 filter = filters + phase * NB_TAPS;
75 sum = s[0] * filter[0] +
83 for(j=0;j<NB_TAPS;j++)
84 sum += s[j] * filter[j];
87 sum = sum >> FILTER_BITS;
98 /* This function must be optimized */
99 static void v_resample(UINT8 *dst, int dst_width, UINT8 *src, int wrap,
106 for(i=0;i<dst_width;i++) {
108 sum = s[0 * wrap] * filter[0] +
109 s[1 * wrap] * filter[1] +
110 s[2 * wrap] * filter[2] +
111 s[3 * wrap] * filter[3];
118 for(j=0;j<NB_TAPS;j++) {
119 sum += s1[0] * filter[j];
124 sum = sum >> FILTER_BITS;
137 #include "i386/mmx.h"
139 #define FILTER4(reg) \
141 s = src + (src_pos >> POS_FRAC_BITS);\
142 phase = get_phase(src_pos);\
143 filter = filters + phase * NB_TAPS;\
145 punpcklbw_r2r(mm7, reg);\
146 movq_m2r(*filter, mm6);\
147 pmaddwd_r2r(reg, mm6);\
150 paddd_r2r(mm6, reg);\
151 psrad_i2r(FILTER_BITS, reg);\
152 src_pos += src_incr;\
155 #define DUMP(reg) movq_r2m(reg, tmp); printf(#reg "=%016Lx\n", tmp.uq);
157 /* XXX: do four pixels at a time */
158 static void h_resample_fast4_mmx(UINT8 *dst, int dst_width, UINT8 *src, int src_width,
159 int src_start, int src_incr, INT16 *filters)
169 while (dst_width >= 4) {
176 packuswb_r2r(mm7, mm0);
177 packuswb_r2r(mm7, mm1);
178 packuswb_r2r(mm7, mm3);
179 packuswb_r2r(mm7, mm2);
191 while (dst_width > 0) {
193 packuswb_r2r(mm7, mm0);
202 static void v_resample4_mmx(UINT8 *dst, int dst_width, UINT8 *src, int wrap,
220 while (dst_width >= 4) {
221 movq_m2r(s[0 * wrap], mm0);
222 punpcklbw_r2r(mm7, mm0);
223 movq_m2r(s[1 * wrap], mm1);
224 punpcklbw_r2r(mm7, mm1);
225 movq_m2r(s[2 * wrap], mm2);
226 punpcklbw_r2r(mm7, mm2);
227 movq_m2r(s[3 * wrap], mm3);
228 punpcklbw_r2r(mm7, mm3);
230 pmullw_m2r(coefs[0], mm0);
231 pmullw_m2r(coefs[1], mm1);
232 pmullw_m2r(coefs[2], mm2);
233 pmullw_m2r(coefs[3], mm3);
238 psraw_i2r(FILTER_BITS, mm0);
240 packuswb_r2r(mm7, mm0);
243 *(UINT32 *)dst = tmp.ud[0];
248 while (dst_width > 0) {
249 sum = s[0 * wrap] * filter[0] +
250 s[1 * wrap] * filter[1] +
251 s[2 * wrap] * filter[2] +
252 s[3 * wrap] * filter[3];
253 sum = sum >> FILTER_BITS;
267 /* slow version to handle limit cases. Does not need optimisation */
268 static void h_resample_slow(UINT8 *dst, int dst_width, UINT8 *src, int src_width,
269 int src_start, int src_incr, INT16 *filters)
271 int src_pos, phase, sum, j, v, i;
275 src_end = src + src_width;
277 for(i=0;i<dst_width;i++) {
278 s = src + (src_pos >> POS_FRAC_BITS);
279 phase = get_phase(src_pos);
280 filter = filters + phase * NB_TAPS;
282 for(j=0;j<NB_TAPS;j++) {
285 else if (s >= src_end)
289 sum += v * filter[j];
292 sum = sum >> FILTER_BITS;
303 static void h_resample(UINT8 *dst, int dst_width, UINT8 *src, int src_width,
304 int src_start, int src_incr, INT16 *filters)
309 n = (0 - src_start + src_incr - 1) / src_incr;
310 h_resample_slow(dst, n, src, src_width, src_start, src_incr, filters);
313 src_start += n * src_incr;
315 src_end = src_start + dst_width * src_incr;
316 if (src_end > ((src_width - NB_TAPS) << POS_FRAC_BITS)) {
317 n = (((src_width - NB_TAPS + 1) << POS_FRAC_BITS) - 1 - src_start) /
323 if ((mm_flags & MM_MMX) && NB_TAPS == 4)
324 h_resample_fast4_mmx(dst, n,
325 src, src_width, src_start, src_incr, filters);
328 h_resample_fast(dst, n,
329 src, src_width, src_start, src_incr, filters);
333 src_start += n * src_incr;
334 h_resample_slow(dst, dst_width,
335 src, src_width, src_start, src_incr, filters);
339 static void component_resample(ImgReSampleContext *s,
340 UINT8 *output, int owrap, int owidth, int oheight,
341 UINT8 *input, int iwrap, int iwidth, int iheight)
343 int src_y, src_y1, last_src_y, ring_y, phase_y, y1, y;
344 UINT8 *new_line, *src_line;
346 last_src_y = - FCENTER - 1;
347 /* position of the bottom of the filter in the source image */
348 src_y = (last_src_y + NB_TAPS) * POS_FRAC;
349 ring_y = NB_TAPS; /* position in ring buffer */
350 for(y=0;y<oheight;y++) {
351 /* apply horizontal filter on new lines from input if needed */
352 src_y1 = src_y >> POS_FRAC_BITS;
353 while (last_src_y < src_y1) {
354 if (++ring_y >= LINE_BUF_HEIGHT + NB_TAPS)
357 /* handle limit conditions : replicate line (slightly
358 inefficient because we filter multiple times) */
362 } else if (y1 >= iheight) {
365 src_line = input + y1 * iwrap;
366 new_line = s->line_buf + ring_y * owidth;
367 /* apply filter and handle limit cases correctly */
368 h_resample(new_line, owidth,
369 src_line, iwidth, - FCENTER * POS_FRAC, s->h_incr,
370 &s->h_filters[0][0]);
371 /* handle ring buffer wraping */
372 if (ring_y >= LINE_BUF_HEIGHT) {
373 memcpy(s->line_buf + (ring_y - LINE_BUF_HEIGHT) * owidth,
377 /* apply vertical filter */
378 phase_y = get_phase(src_y);
380 /* desactivated MMX because loss of precision */
381 if ((mm_flags & MM_MMX) && NB_TAPS == 4 && 0)
382 v_resample4_mmx(output, owidth,
383 s->line_buf + (ring_y - NB_TAPS + 1) * owidth, owidth,
384 &s->v_filters[phase_y][0]);
387 v_resample(output, owidth,
388 s->line_buf + (ring_y - NB_TAPS + 1) * owidth, owidth,
389 &s->v_filters[phase_y][0]);
396 /* XXX: the following filter is quite naive, but it seems to suffice
398 static void build_filter(INT16 *filter, float factor)
401 float x, y, tab[NB_TAPS], norm, mult;
403 /* if upsampling, only need to interpolate, no filter */
407 for(ph=0;ph<NB_PHASES;ph++) {
409 for(i=0;i<NB_TAPS;i++) {
411 x = M_PI * ((float)(i - FCENTER) - (float)ph / NB_PHASES) * factor;
420 /* normalize so that an uniform color remains the same */
421 mult = (float)(1 << FILTER_BITS) / norm;
422 for(i=0;i<NB_TAPS;i++) {
423 v = (int)(tab[i] * mult);
424 filter[ph * NB_TAPS + i] = v;
429 ImgReSampleContext *img_resample_init(int owidth, int oheight,
430 int iwidth, int iheight)
432 return img_resample_full_init(owidth, oheight, iwidth, iheight, 0, 0, 0, 0);
435 ImgReSampleContext *img_resample_full_init(int owidth, int oheight,
436 int iwidth, int iheight,
437 int topBand, int bottomBand,
438 int leftBand, int rightBand)
440 ImgReSampleContext *s;
442 s = av_mallocz(sizeof(ImgReSampleContext));
445 s->line_buf = av_mallocz(owidth * (LINE_BUF_HEIGHT + NB_TAPS));
450 s->oheight = oheight;
452 s->iheight = iheight;
453 s->topBand = topBand;
454 s->bottomBand = bottomBand;
455 s->leftBand = leftBand;
456 s->rightBand = rightBand;
458 s->h_incr = ((iwidth - leftBand - rightBand) * POS_FRAC) / owidth;
459 s->v_incr = ((iheight - topBand - bottomBand) * POS_FRAC) / oheight;
461 build_filter(&s->h_filters[0][0], (float) owidth / (float) (iwidth - leftBand - rightBand));
462 build_filter(&s->v_filters[0][0], (float) oheight / (float) (iheight - topBand - bottomBand));
470 void img_resample(ImgReSampleContext *s,
471 AVPicture *output, AVPicture *input)
476 shift = (i == 0) ? 0 : 1;
477 component_resample(s, output->data[i], output->linesize[i],
478 s->owidth >> shift, s->oheight >> shift,
479 input->data[i] + (input->linesize[i] * (s->topBand >> shift)) + (s->leftBand >> shift),
480 input->linesize[i], ((s->iwidth - s->leftBand - s->rightBand) >> shift),
481 (s->iheight - s->topBand - s->bottomBand) >> shift);
485 void img_resample_close(ImgReSampleContext *s)
487 av_free(s->line_buf);
493 void *av_mallocz(int size)
497 memset(ptr, 0, size);
501 void av_free(void *ptr)
503 /* XXX: this test should not be needed on most libcs */
511 UINT8 img[XSIZE * YSIZE];
516 UINT8 img1[XSIZE1 * YSIZE1];
517 UINT8 img2[XSIZE1 * YSIZE1];
519 void save_pgm(const char *filename, UINT8 *img, int xsize, int ysize)
522 f=fopen(filename,"w");
523 fprintf(f,"P5\n%d %d\n%d\n", xsize, ysize, 255);
524 fwrite(img,1, xsize * ysize,f);
528 static void dump_filter(INT16 *filter)
532 for(ph=0;ph<NB_PHASES;ph++) {
534 for(i=0;i<NB_TAPS;i++) {
535 printf(" %5.2f", filter[ph * NB_TAPS + i] / 256.0);
545 int main(int argc, char **argv)
547 int x, y, v, i, xsize, ysize;
548 ImgReSampleContext *s;
549 float fact, factors[] = { 1/2.0, 3.0/4.0, 1.0, 4.0/3.0, 16.0/9.0, 2.0 };
552 /* build test image */
553 for(y=0;y<YSIZE;y++) {
554 for(x=0;x<XSIZE;x++) {
555 if (x < XSIZE/2 && y < YSIZE/2) {
556 if (x < XSIZE/4 && y < YSIZE/4) {
562 } else if (x < XSIZE/4) {
567 } else if (y < XSIZE/4) {
579 if (((x+3) % 4) <= 1 &&
586 } else if (x < XSIZE/2) {
587 v = ((x - (XSIZE/2)) * 255) / (XSIZE/2);
588 } else if (y < XSIZE/2) {
589 v = ((y - (XSIZE/2)) * 255) / (XSIZE/2);
591 v = ((x + y - XSIZE) * 255) / XSIZE;
593 img[(YSIZE - y) * XSIZE + (XSIZE - x)] = v;
596 save_pgm("/tmp/in.pgm", img, XSIZE, YSIZE);
597 for(i=0;i<sizeof(factors)/sizeof(float);i++) {
599 xsize = (int)(XSIZE * fact);
600 ysize = (int)((YSIZE - 100) * fact);
601 s = img_resample_full_init(xsize, ysize, XSIZE, YSIZE, 50 ,50, 0, 0);
602 printf("Factor=%0.2f\n", fact);
603 dump_filter(&s->h_filters[0][0]);
604 component_resample(s, img1, xsize, xsize, ysize,
605 img + 50 * XSIZE, XSIZE, XSIZE, YSIZE - 100);
606 img_resample_close(s);
608 sprintf(buf, "/tmp/out%d.pgm", i);
609 save_pgm(buf, img1, xsize, ysize);
614 printf("MMX test\n");
616 xsize = (int)(XSIZE * fact);
617 ysize = (int)(YSIZE * fact);
619 s = img_resample_init(xsize, ysize, XSIZE, YSIZE);
620 component_resample(s, img1, xsize, xsize, ysize,
621 img, XSIZE, XSIZE, YSIZE);
624 s = img_resample_init(xsize, ysize, XSIZE, YSIZE);
625 component_resample(s, img2, xsize, xsize, ysize,
626 img, XSIZE, XSIZE, YSIZE);
627 if (memcmp(img1, img2, xsize * ysize) != 0) {
628 fprintf(stderr, "mmx error\n");