2e67fecb3d093133d43bf39b1f123e3ca8add4ec
[ffmpeg.git] / libavcodec / imgresample.c
1 /*
2  * High quality image resampling with polyphase filters 
3  * Copyright (c) 2001 Fabrice Bellard.
4  *
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.
9  *
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.
14  *
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
18  */
19 #include "avcodec.h"
20 #include "dsputil.h"
21
22 #ifdef USE_FASTMEMCPY
23 #include "fastmemcpy.h"
24 #endif
25
26 #define NB_COMPONENTS 3
27
28 #define PHASE_BITS 4
29 #define NB_PHASES  (1 << PHASE_BITS)
30 #define NB_TAPS    4
31 #define FCENTER    1  /* index of the center of the filter */
32 //#define TEST    1  /* Test it */
33
34 #define POS_FRAC_BITS 16
35 #define POS_FRAC      (1 << POS_FRAC_BITS)
36 /* 6 bits precision is needed for MMX */
37 #define FILTER_BITS   8
38
39 #define LINE_BUF_HEIGHT (NB_TAPS * 4)
40
41 struct ImgReSampleContext {
42     int iwidth, iheight, owidth, oheight, topBand, bottomBand, leftBand, rightBand;
43     int h_incr, v_incr;
44     int16_t h_filters[NB_PHASES][NB_TAPS] __align8; /* horizontal filters */
45     int16_t v_filters[NB_PHASES][NB_TAPS] __align8; /* vertical filters */
46     uint8_t *line_buf;
47 };
48
49 static inline int get_phase(int pos)
50 {
51     return ((pos) >> (POS_FRAC_BITS - PHASE_BITS)) & ((1 << PHASE_BITS) - 1);
52 }
53
54 /* This function must be optimized */
55 static void h_resample_fast(uint8_t *dst, int dst_width, uint8_t *src, int src_width,
56                             int src_start, int src_incr, int16_t *filters)
57 {
58     int src_pos, phase, sum, i;
59     uint8_t *s;
60     int16_t *filter;
61
62     src_pos = src_start;
63     for(i=0;i<dst_width;i++) {
64 #ifdef TEST
65         /* test */
66         if ((src_pos >> POS_FRAC_BITS) < 0 ||
67             (src_pos >> POS_FRAC_BITS) > (src_width - NB_TAPS))
68             av_abort();
69 #endif
70         s = src + (src_pos >> POS_FRAC_BITS);
71         phase = get_phase(src_pos);
72         filter = filters + phase * NB_TAPS;
73 #if NB_TAPS == 4
74         sum = s[0] * filter[0] +
75             s[1] * filter[1] +
76             s[2] * filter[2] +
77             s[3] * filter[3];
78 #else
79         {
80             int j;
81             sum = 0;
82             for(j=0;j<NB_TAPS;j++)
83                 sum += s[j] * filter[j];
84         }
85 #endif
86         sum = sum >> FILTER_BITS;
87         if (sum < 0)
88             sum = 0;
89         else if (sum > 255)
90             sum = 255;
91         dst[0] = sum;
92         src_pos += src_incr;
93         dst++;
94     }
95 }
96
97 /* This function must be optimized */
98 static void v_resample(uint8_t *dst, int dst_width, uint8_t *src, int wrap, 
99                        int16_t *filter)
100 {
101     int sum, i;
102     uint8_t *s;
103
104     s = src;
105     for(i=0;i<dst_width;i++) {
106 #if NB_TAPS == 4
107         sum = s[0 * wrap] * filter[0] +
108             s[1 * wrap] * filter[1] +
109             s[2 * wrap] * filter[2] +
110             s[3 * wrap] * filter[3];
111 #else
112         {
113             int j;
114             uint8_t *s1 = s;
115
116             sum = 0;
117             for(j=0;j<NB_TAPS;j++) {
118                 sum += s1[0] * filter[j];
119                 s1 += wrap;
120             }
121         }
122 #endif
123         sum = sum >> FILTER_BITS;
124         if (sum < 0)
125             sum = 0;
126         else if (sum > 255)
127             sum = 255;
128         dst[0] = sum;
129         dst++;
130         s++;
131     }
132 }
133
134 #ifdef HAVE_MMX
135
136 #include "i386/mmx.h"
137
138 #define FILTER4(reg) \
139 {\
140         s = src + (src_pos >> POS_FRAC_BITS);\
141         phase = get_phase(src_pos);\
142         filter = filters + phase * NB_TAPS;\
143         movq_m2r(*s, reg);\
144         punpcklbw_r2r(mm7, reg);\
145         movq_m2r(*filter, mm6);\
146         pmaddwd_r2r(reg, mm6);\
147         movq_r2r(mm6, reg);\
148         psrlq_i2r(32, reg);\
149         paddd_r2r(mm6, reg);\
150         psrad_i2r(FILTER_BITS, reg);\
151         src_pos += src_incr;\
152 }
153
154 #define DUMP(reg) movq_r2m(reg, tmp); printf(#reg "=%016Lx\n", tmp.uq);
155
156 /* XXX: do four pixels at a time */
157 static void h_resample_fast4_mmx(uint8_t *dst, int dst_width, uint8_t *src, int src_width,
158                                  int src_start, int src_incr, int16_t *filters)
159 {
160     int src_pos, phase;
161     uint8_t *s;
162     int16_t *filter;
163     mmx_t tmp;
164     
165     src_pos = src_start;
166     pxor_r2r(mm7, mm7);
167
168     while (dst_width >= 4) {
169
170         FILTER4(mm0);
171         FILTER4(mm1);
172         FILTER4(mm2);
173         FILTER4(mm3);
174
175         packuswb_r2r(mm7, mm0);
176         packuswb_r2r(mm7, mm1);
177         packuswb_r2r(mm7, mm3);
178         packuswb_r2r(mm7, mm2);
179         movq_r2m(mm0, tmp);
180         dst[0] = tmp.ub[0];
181         movq_r2m(mm1, tmp);
182         dst[1] = tmp.ub[0];
183         movq_r2m(mm2, tmp);
184         dst[2] = tmp.ub[0];
185         movq_r2m(mm3, tmp);
186         dst[3] = tmp.ub[0];
187         dst += 4;
188         dst_width -= 4;
189     }
190     while (dst_width > 0) {
191         FILTER4(mm0);
192         packuswb_r2r(mm7, mm0);
193         movq_r2m(mm0, tmp);
194         dst[0] = tmp.ub[0];
195         dst++;
196         dst_width--;
197     }
198     emms();
199 }
200
201 static void v_resample4_mmx(uint8_t *dst, int dst_width, uint8_t *src, int wrap, 
202                             int16_t *filter)
203 {
204     int sum, i, v;
205     uint8_t *s;
206     mmx_t tmp;
207     mmx_t coefs[4];
208     
209     for(i=0;i<4;i++) {
210         v = filter[i];
211         coefs[i].uw[0] = v;
212         coefs[i].uw[1] = v;
213         coefs[i].uw[2] = v;
214         coefs[i].uw[3] = v;
215     }
216     
217     pxor_r2r(mm7, mm7);
218     s = src;
219     while (dst_width >= 4) {
220         movq_m2r(s[0 * wrap], mm0);
221         punpcklbw_r2r(mm7, mm0);
222         movq_m2r(s[1 * wrap], mm1);
223         punpcklbw_r2r(mm7, mm1);
224         movq_m2r(s[2 * wrap], mm2);
225         punpcklbw_r2r(mm7, mm2);
226         movq_m2r(s[3 * wrap], mm3);
227         punpcklbw_r2r(mm7, mm3);
228
229         pmullw_m2r(coefs[0], mm0);
230         pmullw_m2r(coefs[1], mm1);
231         pmullw_m2r(coefs[2], mm2);
232         pmullw_m2r(coefs[3], mm3);
233
234         paddw_r2r(mm1, mm0);
235         paddw_r2r(mm3, mm2);
236         paddw_r2r(mm2, mm0);
237         psraw_i2r(FILTER_BITS, mm0);
238         
239         packuswb_r2r(mm7, mm0);
240         movq_r2m(mm0, tmp);
241
242         *(uint32_t *)dst = tmp.ud[0];
243         dst += 4;
244         s += 4;
245         dst_width -= 4;
246     }
247     while (dst_width > 0) {
248         sum = s[0 * wrap] * filter[0] +
249             s[1 * wrap] * filter[1] +
250             s[2 * wrap] * filter[2] +
251             s[3 * wrap] * filter[3];
252         sum = sum >> FILTER_BITS;
253         if (sum < 0)
254             sum = 0;
255         else if (sum > 255)
256             sum = 255;
257         dst[0] = sum;
258         dst++;
259         s++;
260         dst_width--;
261     }
262     emms();
263 }
264 #endif
265
266 #ifdef HAVE_ALTIVEC
267 typedef union {
268     vector unsigned char v;
269     unsigned char c[16];
270 } vec_uc_t;
271
272 typedef union {
273     vector signed short v;
274     signed short s[8];
275 } vec_ss_t;
276
277 void v_resample16_altivec(uint8_t *dst, int dst_width, uint8_t *src, int wrap,
278                             int16_t *filter)
279 {
280     int sum, i;
281     uint8_t *s;
282     vector unsigned char *tv, tmp, dstv, zero;
283     vec_ss_t srchv[4], srclv[4], fv[4];
284     vector signed short zeros, sumhv, sumlv;    
285     s = src;
286
287     for(i=0;i<4;i++)
288     {
289         /*
290            The vec_madds later on does an implicit >>15 on the result.
291            Since FILTER_BITS is 8, and we have 15 bits of magnitude in
292            a signed short, we have just enough bits to pre-shift our
293            filter constants <<7 to compensate for vec_madds.
294         */
295         fv[i].s[0] = filter[i] << (15-FILTER_BITS);
296         fv[i].v = vec_splat(fv[i].v, 0);
297     }
298     
299     zero = vec_splat_u8(0);
300     zeros = vec_splat_s16(0);
301
302
303     /*
304        When we're resampling, we'd ideally like both our input buffers,
305        and output buffers to be 16-byte aligned, so we can do both aligned
306        reads and writes. Sadly we can't always have this at the moment, so
307        we opt for aligned writes, as unaligned writes have a huge overhead.
308        To do this, do enough scalar resamples to get dst 16-byte aligned.
309     */
310     i = (-(int)dst) & 0xf;
311     while(i>0) {
312         sum = s[0 * wrap] * filter[0] +
313         s[1 * wrap] * filter[1] +
314         s[2 * wrap] * filter[2] +
315         s[3 * wrap] * filter[3];
316         sum = sum >> FILTER_BITS;
317         if (sum<0) sum = 0; else if (sum>255) sum=255;
318         dst[0] = sum;
319         dst++;
320         s++;
321         dst_width--;
322         i--;
323     }
324     
325     /* Do our altivec resampling on 16 pixels at once. */
326     while(dst_width>=16) {
327         /*
328            Read 16 (potentially unaligned) bytes from each of
329            4 lines into 4 vectors, and split them into shorts.
330            Interleave the multipy/accumulate for the resample
331            filter with the loads to hide the 3 cycle latency
332            the vec_madds have.
333         */
334         tv = (vector unsigned char *) &s[0 * wrap];
335         tmp = vec_perm(tv[0], tv[1], vec_lvsl(0, &s[i * wrap]));
336         srchv[0].v = (vector signed short) vec_mergeh(zero, tmp);
337         srclv[0].v = (vector signed short) vec_mergel(zero, tmp);
338         sumhv = vec_madds(srchv[0].v, fv[0].v, zeros);
339         sumlv = vec_madds(srclv[0].v, fv[0].v, zeros);
340
341         tv = (vector unsigned char *) &s[1 * wrap];
342         tmp = vec_perm(tv[0], tv[1], vec_lvsl(0, &s[1 * wrap]));
343         srchv[1].v = (vector signed short) vec_mergeh(zero, tmp);
344         srclv[1].v = (vector signed short) vec_mergel(zero, tmp);
345         sumhv = vec_madds(srchv[1].v, fv[1].v, sumhv);
346         sumlv = vec_madds(srclv[1].v, fv[1].v, sumlv);
347
348         tv = (vector unsigned char *) &s[2 * wrap];
349         tmp = vec_perm(tv[0], tv[1], vec_lvsl(0, &s[2 * wrap]));
350         srchv[2].v = (vector signed short) vec_mergeh(zero, tmp);
351         srclv[2].v = (vector signed short) vec_mergel(zero, tmp);
352         sumhv = vec_madds(srchv[2].v, fv[2].v, sumhv);
353         sumlv = vec_madds(srclv[2].v, fv[2].v, sumlv);
354
355         tv = (vector unsigned char *) &s[3 * wrap];
356         tmp = vec_perm(tv[0], tv[1], vec_lvsl(0, &s[3 * wrap]));
357         srchv[3].v = (vector signed short) vec_mergeh(zero, tmp);
358         srclv[3].v = (vector signed short) vec_mergel(zero, tmp);
359         sumhv = vec_madds(srchv[3].v, fv[3].v, sumhv);
360         sumlv = vec_madds(srclv[3].v, fv[3].v, sumlv);
361     
362         /*
363            Pack the results into our destination vector,
364            and do an aligned write of that back to memory.
365         */
366         dstv = vec_packsu(sumhv, sumlv) ;
367         vec_st(dstv, 0, (vector unsigned char *) dst);
368         
369         dst+=16;
370         s+=16;
371         dst_width-=16;
372     }
373
374     /*
375        If there are any leftover pixels, resample them
376        with the slow scalar method.
377     */
378     while(dst_width>0) {
379         sum = s[0 * wrap] * filter[0] +
380         s[1 * wrap] * filter[1] +
381         s[2 * wrap] * filter[2] +
382         s[3 * wrap] * filter[3];
383         sum = sum >> FILTER_BITS;
384         if (sum<0) sum = 0; else if (sum>255) sum=255;
385         dst[0] = sum;
386         dst++;
387         s++;
388         dst_width--;
389     }
390 }
391 #endif
392
393 /* slow version to handle limit cases. Does not need optimisation */
394 static void h_resample_slow(uint8_t *dst, int dst_width, uint8_t *src, int src_width,
395                             int src_start, int src_incr, int16_t *filters)
396 {
397     int src_pos, phase, sum, j, v, i;
398     uint8_t *s, *src_end;
399     int16_t *filter;
400
401     src_end = src + src_width;
402     src_pos = src_start;
403     for(i=0;i<dst_width;i++) {
404         s = src + (src_pos >> POS_FRAC_BITS);
405         phase = get_phase(src_pos);
406         filter = filters + phase * NB_TAPS;
407         sum = 0;
408         for(j=0;j<NB_TAPS;j++) {
409             if (s < src)
410                 v = src[0];
411             else if (s >= src_end)
412                 v = src_end[-1];
413             else
414                 v = s[0];
415             sum += v * filter[j];
416             s++;
417         }
418         sum = sum >> FILTER_BITS;
419         if (sum < 0)
420             sum = 0;
421         else if (sum > 255)
422             sum = 255;
423         dst[0] = sum;
424         src_pos += src_incr;
425         dst++;
426     }
427 }
428
429 static void h_resample(uint8_t *dst, int dst_width, uint8_t *src, int src_width,
430                        int src_start, int src_incr, int16_t *filters)
431 {
432     int n, src_end;
433
434     if (src_start < 0) {
435         n = (0 - src_start + src_incr - 1) / src_incr;
436         h_resample_slow(dst, n, src, src_width, src_start, src_incr, filters);
437         dst += n;
438         dst_width -= n;
439         src_start += n * src_incr;
440     }
441     src_end = src_start + dst_width * src_incr;
442     if (src_end > ((src_width - NB_TAPS) << POS_FRAC_BITS)) {
443         n = (((src_width - NB_TAPS + 1) << POS_FRAC_BITS) - 1 - src_start) / 
444             src_incr;
445     } else {
446         n = dst_width;
447     }
448 #ifdef HAVE_MMX
449     if ((mm_flags & MM_MMX) && NB_TAPS == 4)
450         h_resample_fast4_mmx(dst, n, 
451                              src, src_width, src_start, src_incr, filters);
452     else
453 #endif
454         h_resample_fast(dst, n, 
455                         src, src_width, src_start, src_incr, filters);
456     if (n < dst_width) {
457         dst += n;
458         dst_width -= n;
459         src_start += n * src_incr;
460         h_resample_slow(dst, dst_width, 
461                         src, src_width, src_start, src_incr, filters);
462     }
463 }
464
465 static void component_resample(ImgReSampleContext *s, 
466                                uint8_t *output, int owrap, int owidth, int oheight,
467                                uint8_t *input, int iwrap, int iwidth, int iheight)
468 {
469     int src_y, src_y1, last_src_y, ring_y, phase_y, y1, y;
470     uint8_t *new_line, *src_line;
471
472     last_src_y = - FCENTER - 1;
473     /* position of the bottom of the filter in the source image */
474     src_y = (last_src_y + NB_TAPS) * POS_FRAC; 
475     ring_y = NB_TAPS; /* position in ring buffer */
476     for(y=0;y<oheight;y++) {
477         /* apply horizontal filter on new lines from input if needed */
478         src_y1 = src_y >> POS_FRAC_BITS;
479         while (last_src_y < src_y1) {
480             if (++ring_y >= LINE_BUF_HEIGHT + NB_TAPS)
481                 ring_y = NB_TAPS;
482             last_src_y++;
483             /* handle limit conditions : replicate line (slightly
484                inefficient because we filter multiple times) */
485             y1 = last_src_y;
486             if (y1 < 0) {
487                 y1 = 0;
488             } else if (y1 >= iheight) {
489                 y1 = iheight - 1;
490             }
491             src_line = input + y1 * iwrap;
492             new_line = s->line_buf + ring_y * owidth;
493             /* apply filter and handle limit cases correctly */
494             h_resample(new_line, owidth, 
495                        src_line, iwidth, - FCENTER * POS_FRAC, s->h_incr, 
496                        &s->h_filters[0][0]);
497             /* handle ring buffer wraping */
498             if (ring_y >= LINE_BUF_HEIGHT) {
499                 memcpy(s->line_buf + (ring_y - LINE_BUF_HEIGHT) * owidth,
500                        new_line, owidth);
501             }
502         }
503         /* apply vertical filter */
504         phase_y = get_phase(src_y);
505 #ifdef HAVE_MMX
506         /* desactivated MMX because loss of precision */
507         if ((mm_flags & MM_MMX) && NB_TAPS == 4 && 0)
508             v_resample4_mmx(output, owidth, 
509                             s->line_buf + (ring_y - NB_TAPS + 1) * owidth, owidth, 
510                             &s->v_filters[phase_y][0]);
511         else
512 #endif
513 #ifdef HAVE_ALTIVEC
514             if ((mm_flags & MM_ALTIVEC) && NB_TAPS == 4 && FILTER_BITS <= 6)
515                 v_resample16_altivec(output, owidth,
516                                 s->line_buf + (ring_y - NB_TAPS + 1) * owidth, owidth,
517                                 &s->v_filters[phase_y][0]);
518         else
519 #endif
520             v_resample(output, owidth, 
521                        s->line_buf + (ring_y - NB_TAPS + 1) * owidth, owidth, 
522                        &s->v_filters[phase_y][0]);
523             
524         src_y += s->v_incr;
525         output += owrap;
526     }
527 }
528
529 /* XXX: the following filter is quite naive, but it seems to suffice
530    for 4 taps */
531 static void build_filter(int16_t *filter, float factor)
532 {
533     int ph, i, v;
534     float x, y, tab[NB_TAPS], norm, mult;
535
536     /* if upsampling, only need to interpolate, no filter */
537     if (factor > 1.0)
538         factor = 1.0;
539
540     for(ph=0;ph<NB_PHASES;ph++) {
541         norm = 0;
542         for(i=0;i<NB_TAPS;i++) {
543             
544             x = M_PI * ((float)(i - FCENTER) - (float)ph / NB_PHASES) * factor;
545             if (x == 0)
546                 y = 1.0;
547             else
548                 y = sin(x) / x;
549             tab[i] = y;
550             norm += y;
551         }
552
553         /* normalize so that an uniform color remains the same */
554         mult = (float)(1 << FILTER_BITS) / norm;
555         for(i=0;i<NB_TAPS;i++) {
556             v = (int)(tab[i] * mult);
557             filter[ph * NB_TAPS + i] = v;
558         }
559     }
560 }
561
562 ImgReSampleContext *img_resample_init(int owidth, int oheight,
563                                       int iwidth, int iheight)
564 {
565         return img_resample_full_init(owidth, oheight, iwidth, iheight, 0, 0, 0, 0);
566 }
567
568 ImgReSampleContext *img_resample_full_init(int owidth, int oheight,
569                                       int iwidth, int iheight,
570                                       int topBand, int bottomBand,
571                                       int leftBand, int rightBand)
572 {
573     ImgReSampleContext *s;
574
575     s = av_mallocz(sizeof(ImgReSampleContext));
576     if (!s)
577         return NULL;
578     s->line_buf = av_mallocz(owidth * (LINE_BUF_HEIGHT + NB_TAPS));
579     if (!s->line_buf) 
580         goto fail;
581     
582     s->owidth = owidth;
583     s->oheight = oheight;
584     s->iwidth = iwidth;
585     s->iheight = iheight;
586     s->topBand = topBand;
587     s->bottomBand = bottomBand;
588     s->leftBand = leftBand;
589     s->rightBand = rightBand;
590     
591     s->h_incr = ((iwidth - leftBand - rightBand) * POS_FRAC) / owidth;
592     s->v_incr = ((iheight - topBand - bottomBand) * POS_FRAC) / oheight;
593     
594     build_filter(&s->h_filters[0][0], (float) owidth  / (float) (iwidth - leftBand - rightBand));
595     build_filter(&s->v_filters[0][0], (float) oheight / (float) (iheight - topBand - bottomBand));
596
597     return s;
598  fail:
599     av_free(s);
600     return NULL;
601 }
602
603 void img_resample(ImgReSampleContext *s, 
604                   AVPicture *output, AVPicture *input)
605 {
606     int i, shift;
607
608     for(i=0;i<3;i++) {
609         shift = (i == 0) ? 0 : 1;
610         component_resample(s, output->data[i], output->linesize[i], 
611                            s->owidth >> shift, s->oheight >> shift,
612                            input->data[i] + (input->linesize[i] * (s->topBand >> shift)) + (s->leftBand >> shift),
613                            input->linesize[i], ((s->iwidth - s->leftBand - s->rightBand) >> shift),
614                            (s->iheight - s->topBand - s->bottomBand) >> shift);
615     }
616 }
617
618 void img_resample_close(ImgReSampleContext *s)
619 {
620     av_free(s->line_buf);
621     av_free(s);
622 }
623
624 #ifdef TEST
625
626 void *av_mallocz(int size)
627 {
628     void *ptr;
629     ptr = malloc(size);
630     memset(ptr, 0, size);
631     return ptr;
632 }
633
634 void av_free(void *ptr)
635 {
636     /* XXX: this test should not be needed on most libcs */
637     if (ptr)
638         free(ptr);
639 }
640
641 /* input */
642 #define XSIZE 256
643 #define YSIZE 256
644 uint8_t img[XSIZE * YSIZE];
645
646 /* output */
647 #define XSIZE1 512
648 #define YSIZE1 512
649 uint8_t img1[XSIZE1 * YSIZE1];
650 uint8_t img2[XSIZE1 * YSIZE1];
651
652 void save_pgm(const char *filename, uint8_t *img, int xsize, int ysize)
653 {
654     FILE *f;
655     f=fopen(filename,"w");
656     fprintf(f,"P5\n%d %d\n%d\n", xsize, ysize, 255);
657     fwrite(img,1, xsize * ysize,f);
658     fclose(f);
659 }
660
661 static void dump_filter(int16_t *filter)
662 {
663     int i, ph;
664
665     for(ph=0;ph<NB_PHASES;ph++) {
666         printf("%2d: ", ph);
667         for(i=0;i<NB_TAPS;i++) {
668             printf(" %5.2f", filter[ph * NB_TAPS + i] / 256.0);
669         }
670         printf("\n");
671     }
672 }
673
674 #ifdef HAVE_MMX
675 int mm_flags;
676 #endif
677
678 int main(int argc, char **argv)
679 {
680     int x, y, v, i, xsize, ysize;
681     ImgReSampleContext *s;
682     float fact, factors[] = { 1/2.0, 3.0/4.0, 1.0, 4.0/3.0, 16.0/9.0, 2.0 };
683     char buf[256];
684
685     /* build test image */
686     for(y=0;y<YSIZE;y++) {
687         for(x=0;x<XSIZE;x++) {
688             if (x < XSIZE/2 && y < YSIZE/2) {
689                 if (x < XSIZE/4 && y < YSIZE/4) {
690                     if ((x % 10) <= 6 &&
691                         (y % 10) <= 6)
692                         v = 0xff;
693                     else
694                         v = 0x00;
695                 } else if (x < XSIZE/4) {
696                     if (x & 1) 
697                         v = 0xff;
698                     else 
699                         v = 0;
700                 } else if (y < XSIZE/4) {
701                     if (y & 1) 
702                         v = 0xff;
703                     else 
704                         v = 0;
705                 } else {
706                     if (y < YSIZE*3/8) {
707                         if ((y+x) & 1) 
708                             v = 0xff;
709                         else 
710                             v = 0;
711                     } else {
712                         if (((x+3) % 4) <= 1 &&
713                             ((y+3) % 4) <= 1)
714                             v = 0xff;
715                         else
716                             v = 0x00;
717                     }
718                 }
719             } else if (x < XSIZE/2) {
720                 v = ((x - (XSIZE/2)) * 255) / (XSIZE/2);
721             } else if (y < XSIZE/2) {
722                 v = ((y - (XSIZE/2)) * 255) / (XSIZE/2);
723             } else {
724                 v = ((x + y - XSIZE) * 255) / XSIZE;
725             }
726             img[(YSIZE - y) * XSIZE + (XSIZE - x)] = v;
727         }
728     }
729     save_pgm("/tmp/in.pgm", img, XSIZE, YSIZE);
730     for(i=0;i<sizeof(factors)/sizeof(float);i++) {
731         fact = factors[i];
732         xsize = (int)(XSIZE * fact);
733         ysize = (int)((YSIZE - 100) * fact);
734         s = img_resample_full_init(xsize, ysize, XSIZE, YSIZE, 50 ,50, 0, 0);
735         printf("Factor=%0.2f\n", fact);
736         dump_filter(&s->h_filters[0][0]);
737         component_resample(s, img1, xsize, xsize, ysize,
738                            img + 50 * XSIZE, XSIZE, XSIZE, YSIZE - 100);
739         img_resample_close(s);
740
741         sprintf(buf, "/tmp/out%d.pgm", i);
742         save_pgm(buf, img1, xsize, ysize);
743     }
744
745     /* mmx test */
746 #ifdef HAVE_MMX
747     printf("MMX test\n");
748     fact = 0.72;
749     xsize = (int)(XSIZE * fact);
750     ysize = (int)(YSIZE * fact);
751     mm_flags = MM_MMX;
752     s = img_resample_init(xsize, ysize, XSIZE, YSIZE);
753     component_resample(s, img1, xsize, xsize, ysize,
754                        img, XSIZE, XSIZE, YSIZE);
755
756     mm_flags = 0;
757     s = img_resample_init(xsize, ysize, XSIZE, YSIZE);
758     component_resample(s, img2, xsize, xsize, ysize,
759                        img, XSIZE, XSIZE, YSIZE);
760     if (memcmp(img1, img2, xsize * ysize) != 0) {
761         fprintf(stderr, "mmx error\n");
762         exit(1);
763     }
764     printf("MMX OK\n");
765 #endif
766     return 0;
767 }
768
769 #endif