spec: Use %license macro to copy license
[platform/upstream/libtheora.git] / examples / dump_psnr.c
1 /********************************************************************
2  *                                                                  *
3  * THIS FILE IS PART OF THE OggTheora SOFTWARE CODEC SOURCE CODE.   *
4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
7  *                                                                  *
8  * THE Theora SOURCE CODE IS COPYRIGHT (C) 2002-2009                *
9  * by the Xiph.Org Foundation and contributors http://www.xiph.org/ *
10  *                                                                  *
11  ********************************************************************
12
13   function: example dumpvid application; dumps  Theora streams
14   last mod: $Id: dump_video.c 15675 2009-02-06 09:43:27Z tterribe $
15
16  ********************************************************************/
17
18 #if !defined(_GNU_SOURCE)
19 #define _GNU_SOURCE
20 #endif
21 #if !defined(_LARGEFILE_SOURCE)
22 #define _LARGEFILE_SOURCE
23 #endif
24 #if !defined(_LARGEFILE64_SOURCE)
25 #define _LARGEFILE64_SOURCE
26 #endif
27 #if !defined(_FILE_OFFSET_BITS)
28 #define _FILE_OFFSET_BITS 64
29 #endif
30
31 #include <stdio.h>
32 #if !defined(_WIN32)
33 #include <getopt.h>
34 #include <unistd.h>
35 #else
36 #include "getopt.h"
37 #endif
38 #include <stdlib.h>
39 #include <string.h>
40 #include <sys/timeb.h>
41 #include <sys/types.h>
42 #include <sys/stat.h>
43 /*Yes, yes, we're going to hell.*/
44 #if defined(_WIN32)
45 #include <io.h>
46 #endif
47 #include <fcntl.h>
48 #include <math.h>
49 #include <signal.h>
50 #include "theora/theoradec.h"
51
52 const char *optstring = "fsy";
53 struct option options [] = {
54   {"frame-type",no_argument,NULL,'f'},
55   {"summary",no_argument,NULL,'s'},
56   {"luma-only",no_argument,NULL,'y'},
57   {NULL,0,NULL,0}
58 };
59
60 static int show_frame_type;
61 static int summary_only;
62 static int luma_only;
63
64 typedef struct y4m_input y4m_input;
65
66 /*The function used to perform chroma conversion.*/
67 typedef void (*y4m_convert_func)(y4m_input *_y4m,
68  unsigned char *_dst,unsigned char *_aux);
69
70 struct y4m_input{
71   int               frame_w;
72   int               frame_h;
73   int               pic_w;
74   int               pic_h;
75   int               pic_x;
76   int               pic_y;
77   int               fps_n;
78   int               fps_d;
79   int               par_n;
80   int               par_d;
81   char              interlace;
82   int               src_c_dec_h;
83   int               src_c_dec_v;
84   int               dst_c_dec_h;
85   int               dst_c_dec_v;
86   char              chroma_type[16];
87   /*The size of each converted frame buffer.*/
88   size_t            dst_buf_sz;
89   /*The amount to read directly into the converted frame buffer.*/
90   size_t            dst_buf_read_sz;
91   /*The size of the auxilliary buffer.*/
92   size_t            aux_buf_sz;
93   /*The amount to read into the auxilliary buffer.*/
94   size_t            aux_buf_read_sz;
95   y4m_convert_func  convert;
96   unsigned char    *dst_buf;
97   unsigned char    *aux_buf;
98 };
99
100
101 static int y4m_parse_tags(y4m_input *_y4m,char *_tags){
102   int   got_w;
103   int   got_h;
104   int   got_fps;
105   int   got_interlace;
106   int   got_par;
107   int   got_chroma;
108   char *p;
109   char *q;
110   got_w=got_h=got_fps=got_interlace=got_par=got_chroma=0;
111   for(p=_tags;;p=q){
112     /*Skip any leading spaces.*/
113     while(*p==' ')p++;
114     /*If that's all we have, stop.*/
115     if(p[0]=='\0')break;
116     /*Find the end of this tag.*/
117     for(q=p+1;*q!='\0'&&*q!=' ';q++);
118     /*Process the tag.*/
119     switch(p[0]){
120       case 'W':{
121         if(sscanf(p+1,"%d",&_y4m->pic_w)!=1)return -1;
122         got_w=1;
123       }break;
124       case 'H':{
125         if(sscanf(p+1,"%d",&_y4m->pic_h)!=1)return -1;
126         got_h=1;
127       }break;
128       case 'F':{
129         if(sscanf(p+1,"%d:%d",&_y4m->fps_n,&_y4m->fps_d)!=2){
130           return -1;
131         }
132         got_fps=1;
133       }break;
134       case 'I':{
135         _y4m->interlace=p[1];
136         got_interlace=1;
137       }break;
138       case 'A':{
139         if(sscanf(p+1,"%d:%d",&_y4m->par_n,&_y4m->par_d)!=2){
140           return -1;
141         }
142         got_par=1;
143       }break;
144       case 'C':{
145         if(q-p>16)return -1;
146         memcpy(_y4m->chroma_type,p+1,q-p-1);
147         _y4m->chroma_type[q-p-1]='\0';
148         got_chroma=1;
149       }break;
150       /*Ignore unknown tags.*/
151     }
152   }
153   if(!got_w||!got_h||!got_fps||!got_interlace||!got_par)return -1;
154   /*Chroma-type is not specified in older files, e.g., those generated by
155      mplayer.*/
156   if(!got_chroma)strcpy(_y4m->chroma_type,"420");
157   return 0;
158 }
159
160 /*All anti-aliasing filters in the following conversion functions are based on
161    one of two window functions:
162   The 6-tap Lanczos window (for down-sampling and shifts):
163    sinc(\pi*t)*sinc(\pi*t/3), |t|<3  (sinc(t)==sin(t)/t)
164    0,                         |t|>=3
165   The 4-tap Mitchell window (for up-sampling):
166    7|t|^3-12|t|^2+16/3,             |t|<1
167    -(7/3)|x|^3+12|x|^2-20|x|+32/3,  |t|<2
168    0,                               |t|>=2
169   The number of taps is intentionally kept small to reduce computational
170    overhead and limit ringing.
171
172   The taps from these filters are scaled so that their sum is 1, and the result
173    is scaled by 128 and rounded to integers to create a filter whose
174    intermediate values fit inside 16 bits.
175   Coefficients are rounded in such a way as to ensure their sum is still 128,
176    which is usually equivalent to normal rounding.*/
177
178 #define OC_MINI(_a,_b)      ((_a)>(_b)?(_b):(_a))
179 #define OC_MAXI(_a,_b)      ((_a)<(_b)?(_b):(_a))
180 #define OC_CLAMPI(_a,_b,_c) (OC_MAXI(_a,OC_MINI(_b,_c)))
181
182 /*420jpeg chroma samples are sited like:
183   Y-------Y-------Y-------Y-------
184   |       |       |       |
185   |   BR  |       |   BR  |
186   |       |       |       |
187   Y-------Y-------Y-------Y-------
188   |       |       |       |
189   |       |       |       |
190   |       |       |       |
191   Y-------Y-------Y-------Y-------
192   |       |       |       |
193   |   BR  |       |   BR  |
194   |       |       |       |
195   Y-------Y-------Y-------Y-------
196   |       |       |       |
197   |       |       |       |
198   |       |       |       |
199
200   420mpeg2 chroma samples are sited like:
201   Y-------Y-------Y-------Y-------
202   |       |       |       |
203   BR      |       BR      |
204   |       |       |       |
205   Y-------Y-------Y-------Y-------
206   |       |       |       |
207   |       |       |       |
208   |       |       |       |
209   Y-------Y-------Y-------Y-------
210   |       |       |       |
211   BR      |       BR      |
212   |       |       |       |
213   Y-------Y-------Y-------Y-------
214   |       |       |       |
215   |       |       |       |
216   |       |       |       |
217
218   We use a resampling filter to shift the site locations one quarter pixel (at
219    the chroma plane's resolution) to the right.
220   The 4:2:2 modes look exactly the same, except there are twice as many chroma
221    lines, and they are vertically co-sited with the luma samples in both the
222    mpeg2 and jpeg cases (thus requiring no vertical resampling).*/
223 static void y4m_convert_42xmpeg2_42xjpeg(y4m_input *_y4m,unsigned char *_dst,
224  unsigned char *_aux){
225   int c_w;
226   int c_h;
227   int pli;
228   int y;
229   int x;
230   /*Skip past the luma data.*/
231   _dst+=_y4m->pic_w*_y4m->pic_h;
232   /*Compute the size of each chroma plane.*/
233   c_w=(_y4m->pic_w+_y4m->dst_c_dec_h-1)/_y4m->dst_c_dec_h;
234   c_h=(_y4m->pic_h+_y4m->dst_c_dec_v-1)/_y4m->dst_c_dec_v;
235   for(pli=1;pli<3;pli++){
236     for(y=0;y<c_h;y++){
237       /*Filter: [4 -17 114 35 -9 1]/128, derived from a 6-tap Lanczos
238          window.*/
239       for(x=0;x<OC_MINI(c_w,2);x++){
240         _dst[x]=(unsigned char)OC_CLAMPI(0,4*_aux[0]-17*_aux[OC_MAXI(x-1,0)]+
241          114*_aux[x]+35*_aux[OC_MINI(x+1,c_w-1)]-9*_aux[OC_MINI(x+2,c_w-1)]+
242          _aux[OC_MINI(x+3,c_w-1)]+64>>7,255);
243       }
244       for(;x<c_w-3;x++){
245         _dst[x]=(unsigned char)OC_CLAMPI(0,4*_aux[x-2]-17*_aux[x-1]+
246          114*_aux[x]+35*_aux[x+1]-9*_aux[x+2]+_aux[x+3]+64>>7,255);
247       }
248       for(;x<c_w;x++){
249         _dst[x]=(unsigned char)OC_CLAMPI(0,4*_aux[x-2]-17*_aux[x-1]+
250          114*_aux[x]+35*_aux[OC_MINI(x+1,c_w-1)]-9*_aux[OC_MINI(x+2,c_w-1)]+
251          _aux[c_w-1]+64>>7,255);
252       }
253       _dst+=c_w;
254       _aux+=c_w;
255     }
256   }
257 }
258
259 /*This format is only used for interlaced content, but is included for
260    completeness.
261
262   420jpeg chroma samples are sited like:
263   Y-------Y-------Y-------Y-------
264   |       |       |       |
265   |   BR  |       |   BR  |
266   |       |       |       |
267   Y-------Y-------Y-------Y-------
268   |       |       |       |
269   |       |       |       |
270   |       |       |       |
271   Y-------Y-------Y-------Y-------
272   |       |       |       |
273   |   BR  |       |   BR  |
274   |       |       |       |
275   Y-------Y-------Y-------Y-------
276   |       |       |       |
277   |       |       |       |
278   |       |       |       |
279
280   420paldv chroma samples are sited like:
281   YR------Y-------YR------Y-------
282   |       |       |       |
283   |       |       |       |
284   |       |       |       |
285   YB------Y-------YB------Y-------
286   |       |       |       |
287   |       |       |       |
288   |       |       |       |
289   YR------Y-------YR------Y-------
290   |       |       |       |
291   |       |       |       |
292   |       |       |       |
293   YB------Y-------YB------Y-------
294   |       |       |       |
295   |       |       |       |
296   |       |       |       |
297
298   We use a resampling filter to shift the site locations one quarter pixel (at
299    the chroma plane's resolution) to the right.
300   Then we use another filter to move the C_r location down one quarter pixel,
301    and the C_b location up one quarter pixel.*/
302 static void y4m_convert_42xpaldv_42xjpeg(y4m_input *_y4m,unsigned char *_dst,
303  unsigned char *_aux){
304   unsigned char *tmp;
305   int            c_w;
306   int            c_h;
307   int            c_sz;
308   int            pli;
309   int            y;
310   int            x;
311   /*Skip past the luma data.*/
312   _dst+=_y4m->pic_w*_y4m->pic_h;
313   /*Compute the size of each chroma plane.*/
314   c_w=(_y4m->pic_w+1)/2;
315   c_h=(_y4m->pic_h+_y4m->dst_c_dec_h-1)/_y4m->dst_c_dec_h;
316   c_sz=c_w*c_h;
317   /*First do the horizontal re-sampling.
318     This is the same as the mpeg2 case, except that after the horizontal case,
319      we need to apply a second vertical filter.*/
320   tmp=_aux+2*c_sz;
321   for(pli=1;pli<3;pli++){
322     for(y=0;y<c_h;y++){
323       /*Filter: [4 -17 114 35 -9 1]/128, derived from a 6-tap Lanczos
324          window.*/
325       for(x=0;x<OC_MINI(c_w,2);x++){
326         tmp[x]=(unsigned char)OC_CLAMPI(0,4*_aux[0]-17*_aux[OC_MAXI(x-1,0)]+
327          114*_aux[x]+35*_aux[OC_MINI(x+1,c_w-1)]-9*_aux[OC_MINI(x+2,c_w-1)]+
328          _aux[OC_MINI(x+3,c_w-1)]+64>>7,255);
329       }
330       for(;x<c_w-3;x++){
331         tmp[x]=(unsigned char)OC_CLAMPI(0,4*_aux[x-2]-17*_aux[x-1]+
332          114*_aux[x]+35*_aux[x+1]-9*_aux[x+2]+_aux[x+3]+64>>7,255);
333       }
334       for(;x<c_w;x++){
335         tmp[x]=(unsigned char)OC_CLAMPI(0,4*_aux[x-2]-17*_aux[x-1]+
336          114*_aux[x]+35*_aux[OC_MINI(x+1,c_w-1)]-9*_aux[OC_MINI(x+2,c_w-1)]+
337          _aux[c_w-1]+64>>7,255);
338       }
339       tmp+=c_w;
340       _aux+=c_w;
341     }
342     switch(pli){
343       case 1:{
344         tmp-=c_sz;
345         /*Slide C_b up a quarter-pel.
346           This is the same filter used above, but in the other order.*/
347         for(x=0;x<c_w;x++){
348           for(y=0;y<OC_MINI(c_h,3);y++){
349             _dst[y*c_w]=(unsigned char)OC_CLAMPI(0,tmp[0]-
350              9*tmp[OC_MAXI(y-2,0)*c_w]+35*tmp[OC_MAXI(y-1,0)*c_w]+
351              114*tmp[y*c_w]-17*tmp[OC_MINI(y+1,c_h-1)*c_w]+
352              4*tmp[OC_MINI(y+2,c_h-1)*c_w]+64>>7,255);
353           }
354           for(;y<c_h-2;y++){
355             _dst[y*c_w]=(unsigned char)OC_CLAMPI(0,tmp[(y-3)*c_w]-
356              9*tmp[(y-2)*c_w]+35*tmp[(y-1)*c_w]+114*tmp[y*c_w]-
357              17*tmp[(y+1)*c_w]+4*tmp[(y+2)*c_w]+64>>7,255);
358           }
359           for(;y<c_h;y++){
360             _dst[y*c_w]=(unsigned char)OC_CLAMPI(0,tmp[(y-3)*c_w]-
361              9*tmp[(y-2)*c_w]+35*tmp[(y-1)*c_w]+114*tmp[y*c_w]-
362              17*tmp[OC_MINI(y+1,c_h-1)*c_w]+4*tmp[(c_h-1)*c_w]+64>>7,255);
363           }
364           _dst++;
365           tmp++;
366         }
367         _dst+=c_sz-c_w;
368         tmp-=c_w;
369       }break;
370       case 2:{
371         tmp-=c_sz;
372         /*Slide C_r down a quarter-pel.
373           This is the same as the horizontal filter.*/
374         for(x=0;x<c_w;x++){
375           for(y=0;y<OC_MINI(c_h,2);y++){
376             _dst[y*c_w]=(unsigned char)OC_CLAMPI(0,4*tmp[0]-
377              17*tmp[OC_MAXI(y-1,0)*c_w]+114*tmp[y*c_w]+
378              35*tmp[OC_MINI(y+1,c_h-1)*c_w]-9*tmp[OC_MINI(y+2,c_h-1)*c_w]+
379              tmp[OC_MINI(y+3,c_h-1)*c_w]+64>>7,255);
380           }
381           for(;y<c_h-3;y++){
382             _dst[y*c_w]=(unsigned char)OC_CLAMPI(0,4*tmp[(y-2)*c_w]-
383              17*tmp[(y-1)*c_w]+114*tmp[y*c_w]+35*tmp[(y+1)*c_w]-
384              9*tmp[(y+2)*c_w]+tmp[(y+3)*c_w]+64>>7,255);
385           }
386           for(;y<c_h;y++){
387             _dst[y*c_w]=(unsigned char)OC_CLAMPI(0,4*tmp[(y-2)*c_w]-
388              17*tmp[(y-1)*c_w]+114*tmp[y*c_w]+35*tmp[OC_MINI(y+1,c_h-1)*c_w]-
389              9*tmp[OC_MINI(y+2,c_h-1)*c_w]+tmp[(c_h-1)*c_w]+64>>7,255);
390           }
391           _dst++;
392           tmp++;
393         }
394       }break;
395     }
396     /*For actual interlaced material, this would have to be done separately on
397        each field, and the shift amounts would be different.
398       C_r moves down 1/8, C_b up 3/8 in the top field, and C_r moves down 3/8,
399        C_b up 1/8 in the bottom field.
400       The corresponding filters would be:
401        Down 1/8 (reverse order for up): [3 -11 125 15 -4 0]/128
402        Down 3/8 (reverse order for up): [4 -19 98 56 -13 2]/128*/
403   }
404 }
405
406 /*422jpeg chroma samples are sited like:
407   Y---BR--Y-------Y---BR--Y-------
408   |       |       |       |
409   |       |       |       |
410   |       |       |       |
411   Y---BR--Y-------Y---BR--Y-------
412   |       |       |       |
413   |       |       |       |
414   |       |       |       |
415   Y---BR--Y-------Y---BR--Y-------
416   |       |       |       |
417   |       |       |       |
418   |       |       |       |
419   Y---BR--Y-------Y---BR--Y-------
420   |       |       |       |
421   |       |       |       |
422   |       |       |       |
423
424   411 chroma samples are sited like:
425   YBR-----Y-------Y-------Y-------
426   |       |       |       |
427   |       |       |       |
428   |       |       |       |
429   YBR-----Y-------Y-------Y-------
430   |       |       |       |
431   |       |       |       |
432   |       |       |       |
433   YBR-----Y-------Y-------Y-------
434   |       |       |       |
435   |       |       |       |
436   |       |       |       |
437   YBR-----Y-------Y-------Y-------
438   |       |       |       |
439   |       |       |       |
440   |       |       |       |
441
442   We use a filter to resample at site locations one eighth pixel (at the source
443    chroma plane's horizontal resolution) and five eighths of a pixel to the
444    right.*/
445 static void y4m_convert_411_422jpeg(y4m_input *_y4m,unsigned char *_dst,
446  unsigned char *_aux){
447   int c_w;
448   int dst_c_w;
449   int c_h;
450   int pli;
451   int y;
452   int x;
453   /*Skip past the luma data.*/
454   _dst+=_y4m->pic_w*_y4m->pic_h;
455   /*Compute the size of each chroma plane.*/
456   c_w=(_y4m->pic_w+_y4m->src_c_dec_h-1)/_y4m->src_c_dec_h;
457   dst_c_w=(_y4m->pic_w+_y4m->dst_c_dec_h-1)/_y4m->dst_c_dec_h;
458   c_h=(_y4m->pic_h+_y4m->dst_c_dec_v-1)/_y4m->dst_c_dec_v;
459   for(pli=1;pli<3;pli++){
460     for(y=0;y<c_h;y++){
461       /*Filters: [1 110 18 -1]/128 and [-3 50 86 -5]/128, both derived from a
462          4-tap Mitchell window.*/
463       for(x=0;x<OC_MINI(c_w,1);x++){
464         _dst[x<<1]=(unsigned char)OC_CLAMPI(0,111*_aux[0]+
465          18*_aux[OC_MINI(1,c_w-1)]-_aux[OC_MINI(2,c_w-1)]+64>>7,255);
466         _dst[x<<1|1]=(unsigned char)OC_CLAMPI(0,47*_aux[0]+
467          86*_aux[OC_MINI(1,c_w-1)]-5*_aux[OC_MINI(2,c_w-1)]+64>>7,255);
468       }
469       for(;x<c_w-2;x++){
470         _dst[x<<1]=(unsigned char)OC_CLAMPI(0,_aux[x-1]+110*_aux[x]+
471          18*_aux[x+1]-_aux[x+2]+64>>7,255);
472         _dst[x<<1|1]=(unsigned char)OC_CLAMPI(0,-3*_aux[x-1]+50*_aux[x]+
473          86*_aux[x+1]-5*_aux[x+2]+64>>7,255);
474       }
475       for(;x<c_w;x++){
476         _dst[x<<1]=(unsigned char)OC_CLAMPI(0,_aux[x-1]+110*_aux[x]+
477          18*_aux[OC_MINI(x+1,c_w-1)]-_aux[c_w-1]+64>>7,255);
478         if((x<<1|1)<dst_c_w){
479           _dst[x<<1|1]=(unsigned char)OC_CLAMPI(0,-3*_aux[x-1]+50*_aux[x]+
480            86*_aux[OC_MINI(x+1,c_w-1)]-5*_aux[c_w-1]+64>>7,255);
481         }
482       }
483       _dst+=dst_c_w;
484       _aux+=c_w;
485     }
486   }
487 }
488
489 /*The image is padded with empty chroma components at 4:2:0.
490   This costs about 17 bits a frame to code.*/
491 static void y4m_convert_mono_420jpeg(y4m_input *_y4m,unsigned char *_dst,
492  unsigned char *_aux){
493   int c_sz;
494   _dst+=_y4m->pic_w*_y4m->pic_h;
495   c_sz=((_y4m->pic_w+_y4m->dst_c_dec_h-1)/_y4m->dst_c_dec_h)*
496    ((_y4m->pic_h+_y4m->dst_c_dec_v-1)/_y4m->dst_c_dec_v);
497   memset(_dst,128,c_sz*2);
498 }
499
500 #if 0
501 /*Right now just 444 to 420.
502   Not too hard to generalize.*/
503 static void y4m_convert_4xxjpeg_42xjpeg(y4m_input *_y4m,unsigned char *_dst,
504  unsigned char *_aux){
505   unsigned char *tmp;
506   int            c_w;
507   int            c_h;
508   int            pic_sz;
509   int            tmp_sz;
510   int            c_sz;
511   int            pli;
512   int            y;
513   int            x;
514   /*Compute the size of each chroma plane.*/
515   c_w=(_y4m->pic_w+_y4m->dst_c_dec_h-1)/_y4m->dst_c_dec_h;
516   c_h=(_y4m->pic_h+_y4m->dst_c_dec_v-1)/_y4m->dst_c_dec_v;
517   pic_sz=_y4m->pic_w*_y4m->pic_h;
518   tmp_sz=c_w*_y4m->pic_h;
519   c_sz=c_w*c_h;
520   _dst+=pic_sz;
521   for(pli=1;pli<3;pli++){
522     tmp=_aux+pic_sz;
523     /*In reality, the horizontal and vertical steps could be pipelined, for
524        less memory consumption and better cache performance, but we do them
525        separately for simplicity.*/
526     /*First do horizontal filtering (convert to 4:2:2)*/
527     /*Filter: [3 -17 78 78 -17 3]/128, derived from a 6-tap Lanczos window.*/
528     for(y=0;y<_y4m->pic_h;y++){
529       for(x=0;x<OC_MINI(_y4m->pic_w,2);x+=2){
530         tmp[x>>1]=OC_CLAMPI(0,64*_aux[0]+78*_aux[OC_MINI(1,_y4m->pic_w-1)]
531          -17*_aux[OC_MINI(2,_y4m->pic_w-1)]
532          +3*_aux[OC_MINI(3,_y4m->pic_w-1)]+64>>7,255);
533       }
534       for(;x<_y4m->pic_w-3;x+=2){
535         tmp[x>>1]=OC_CLAMPI(0,3*(_aux[x-2]+_aux[x+3])-17*(_aux[x-1]+_aux[x+2])+
536          78*(_aux[x]+_aux[x+1])+64>>7,255);
537       }
538       for(;x<_y4m->pic_w;x+=2){
539         tmp[x>>1]=OC_CLAMPI(0,3*(_aux[x-2]+_aux[_y4m->pic_w-1])-
540          17*(_aux[x-1]+_aux[OC_MINI(x+2,_y4m->pic_w-1)])+
541          78*(_aux[x]+_aux[OC_MINI(x+1,_y4m->pic_w-1)])+64>>7,255);
542       }
543       tmp+=c_w;
544       _aux+=_y4m->pic_w;
545     }
546     _aux-=pic_sz;
547     tmp-=tmp_sz;
548     /*Now do the vertical filtering.*/
549     for(x=0;x<c_w;x++){
550       for(y=0;y<OC_MINI(_y4m->pic_h,2);y+=2){
551         _dst[(y>>1)*c_w]=OC_CLAMPI(0,64*tmp[0]
552          +78*tmp[OC_MINI(1,_y4m->pic_h-1)*c_w]
553          -17*tmp[OC_MINI(2,_y4m->pic_h-1)*c_w]
554          +3*tmp[OC_MINI(3,_y4m->pic_h-1)*c_w]+64>>7,255);
555       }
556       for(;y<_y4m->pic_h-3;y+=2){
557         _dst[(y>>1)*c_w]=OC_CLAMPI(0,3*(tmp[(y-2)*c_w]+tmp[(y+3)*c_w])-
558          17*(tmp[(y-1)*c_w]+tmp[(y+2)*c_w])+78*(tmp[y*c_w]+tmp[(y+1)*c_w])+
559          64>>7,255);
560       }
561       for(;y<_y4m->pic_h;y+=2){
562         _dst[(y>>1)*c_w]=OC_CLAMPI(0,3*(tmp[(y-2)*c_w]
563          +tmp[(_y4m->pic_h-1)*c_w])-17*(tmp[(y-1)*c_w]
564          +tmp[OC_MINI(y+2,_y4m->pic_h-1)*c_w])
565          +78*(tmp[y*c_w]+tmp[OC_MINI(y+1,_y4m->pic_h-1)*c_w])+64>>7,255);
566       }
567       tmp++;
568       _dst++;
569     }
570     _dst-=c_w;
571   }
572 }
573 #endif
574
575 /*No conversion function needed.*/
576 static void y4m_convert_null(y4m_input *_y4m,unsigned char *_dst,
577  unsigned char *_aux){
578 }
579
580 static int y4m_input_open(y4m_input *_y4m,FILE *_fin,char *_skip,int _nskip){
581   char buffer[80];
582   int  ret;
583   int  i;
584   /*Read until newline, or 80 cols, whichever happens first.*/
585   for(i=0;i<79;i++){
586     if(_nskip>0){
587       buffer[i]=*_skip++;
588       _nskip--;
589     }
590     else{
591       ret=fread(buffer+i,1,1,_fin);
592       if(ret<1)return -1;
593     }
594     if(buffer[i]=='\n')break;
595   }
596   /*We skipped too much header data.*/
597   if(_nskip>0)return -1;
598   if(i==79){
599     fprintf(stderr,"Error parsing header; not a YUV2MPEG2 file?\n");
600     return -1;
601   }
602   buffer[i]='\0';
603   if(memcmp(buffer,"YUV4MPEG",8)){
604     fprintf(stderr,"Incomplete magic for YUV4MPEG file.\n");
605     return -1;
606   }
607   if(buffer[8]!='2'){
608     fprintf(stderr,"Incorrect YUV input file version; YUV4MPEG2 required.\n");
609   }
610   ret=y4m_parse_tags(_y4m,buffer+5);
611   if(ret<0){
612     fprintf(stderr,"Error parsing YUV4MPEG2 header.\n");
613     return ret;
614   }
615   if(_y4m->interlace!='p'){
616     fprintf(stderr,"Input video is interlaced; "
617      "Theora only handles progressive scan.\n");
618     return -1;
619   }
620   if(strcmp(_y4m->chroma_type,"420")==0||
621    strcmp(_y4m->chroma_type,"420jpeg")==0){
622     _y4m->src_c_dec_h=_y4m->dst_c_dec_h=_y4m->src_c_dec_v=_y4m->dst_c_dec_v=2;
623     _y4m->dst_buf_read_sz=_y4m->pic_w*_y4m->pic_h
624      +2*((_y4m->pic_w+1)/2)*((_y4m->pic_h+1)/2);
625     /*Natively supported: no conversion required.*/
626     _y4m->aux_buf_sz=_y4m->aux_buf_read_sz=0;
627     _y4m->convert=y4m_convert_null;
628   }
629   else if(strcmp(_y4m->chroma_type,"420mpeg2")==0){
630     _y4m->src_c_dec_h=_y4m->dst_c_dec_h=_y4m->src_c_dec_v=_y4m->dst_c_dec_v=2;
631     _y4m->dst_buf_read_sz=_y4m->pic_w*_y4m->pic_h;
632     /*Chroma filter required: read into the aux buf first.*/
633     _y4m->aux_buf_sz=_y4m->aux_buf_read_sz=
634      2*((_y4m->pic_w+1)/2)*((_y4m->pic_h+1)/2);
635     _y4m->convert=y4m_convert_42xmpeg2_42xjpeg;
636   }
637   else if(strcmp(_y4m->chroma_type,"420paldv")==0){
638     _y4m->src_c_dec_h=_y4m->dst_c_dec_h=_y4m->src_c_dec_v=_y4m->dst_c_dec_v=2;
639     _y4m->dst_buf_read_sz=_y4m->pic_w*_y4m->pic_h;
640     /*Chroma filter required: read into the aux buf first.
641       We need to make two filter passes, so we need some extra space in the
642        aux buffer.*/
643     _y4m->aux_buf_sz=3*((_y4m->pic_w+1)/2)*((_y4m->pic_h+1)/2);
644     _y4m->aux_buf_read_sz=2*((_y4m->pic_w+1)/2)*((_y4m->pic_h+1)/2);
645     _y4m->convert=y4m_convert_42xpaldv_42xjpeg;
646   }
647   else if(strcmp(_y4m->chroma_type,"422")==0){
648     _y4m->src_c_dec_h=_y4m->dst_c_dec_h=2;
649     _y4m->src_c_dec_v=_y4m->dst_c_dec_v=1;
650     _y4m->dst_buf_read_sz=_y4m->pic_w*_y4m->pic_h;
651     /*Chroma filter required: read into the aux buf first.*/
652     _y4m->aux_buf_sz=_y4m->aux_buf_read_sz=2*((_y4m->pic_w+1)/2)*_y4m->pic_h;
653     _y4m->convert=y4m_convert_42xmpeg2_42xjpeg;
654   }
655   else if(strcmp(_y4m->chroma_type,"411")==0){
656     _y4m->src_c_dec_h=4;
657     /*We don't want to introduce any additional sub-sampling, so we
658        promote 4:1:1 material to 4:2:2, as the closest format Theora can
659        handle.*/
660     _y4m->dst_c_dec_h=2;
661     _y4m->src_c_dec_v=_y4m->dst_c_dec_v=1;
662     _y4m->dst_buf_read_sz=_y4m->pic_w*_y4m->pic_h;
663     /*Chroma filter required: read into the aux buf first.*/
664     _y4m->aux_buf_sz=_y4m->aux_buf_read_sz=2*((_y4m->pic_w+3)/4)*_y4m->pic_h;
665     _y4m->convert=y4m_convert_411_422jpeg;
666   }
667   else if(strcmp(_y4m->chroma_type,"444")==0){
668     _y4m->src_c_dec_h=_y4m->dst_c_dec_h=_y4m->src_c_dec_v=_y4m->dst_c_dec_v=1;
669     _y4m->dst_buf_read_sz=_y4m->pic_w*_y4m->pic_h*3;
670     /*Natively supported: no conversion required.*/
671     _y4m->aux_buf_sz=_y4m->aux_buf_read_sz=0;
672     _y4m->convert=y4m_convert_null;
673   }
674   else if(strcmp(_y4m->chroma_type,"444alpha")==0){
675     _y4m->src_c_dec_h=_y4m->dst_c_dec_h=_y4m->src_c_dec_v=_y4m->dst_c_dec_v=1;
676     _y4m->dst_buf_read_sz=_y4m->pic_w*_y4m->pic_h*3;
677     /*Read the extra alpha plane into the aux buf.
678       It will be discarded.*/
679     _y4m->aux_buf_sz=_y4m->aux_buf_read_sz=_y4m->pic_w*_y4m->pic_h;
680     _y4m->convert=y4m_convert_null;
681   }
682   else if(strcmp(_y4m->chroma_type,"mono")==0){
683     _y4m->src_c_dec_h=_y4m->src_c_dec_v=0;
684     _y4m->dst_c_dec_h=_y4m->dst_c_dec_v=2;
685     _y4m->dst_buf_read_sz=_y4m->pic_w*_y4m->pic_h;
686     /*No extra space required, but we need to clear the chroma planes.*/
687     _y4m->aux_buf_sz=_y4m->aux_buf_read_sz=0;
688     _y4m->convert=y4m_convert_mono_420jpeg;
689   }
690   else{
691     fprintf(stderr,"Unknown chroma sampling type: %s\n",_y4m->chroma_type);
692     return -1;
693   }
694   /*The size of the final frame buffers is always computed from the
695      destination chroma decimation type.*/
696   _y4m->dst_buf_sz=_y4m->pic_w*_y4m->pic_h
697    +2*((_y4m->pic_w+_y4m->dst_c_dec_h-1)/_y4m->dst_c_dec_h)*
698    ((_y4m->pic_h+_y4m->dst_c_dec_v-1)/_y4m->dst_c_dec_v);
699   /*Scale the picture size up to a multiple of 16.*/
700   _y4m->frame_w=_y4m->pic_w+15&~0xF;
701   _y4m->frame_h=_y4m->pic_h+15&~0xF;
702   /*Force the offsets to be even so that chroma samples line up like we
703      expect.*/
704   _y4m->pic_x=_y4m->frame_w-_y4m->pic_w>>1&~1;
705   _y4m->pic_y=_y4m->frame_h-_y4m->pic_h>>1&~1;
706   _y4m->dst_buf=(unsigned char *)malloc(_y4m->dst_buf_sz);
707   _y4m->aux_buf=(unsigned char *)malloc(_y4m->aux_buf_sz);
708   return 0;
709 }
710
711 static void y4m_input_get_info(y4m_input *_y4m,th_info *_ti){
712   _ti->frame_width=_y4m->frame_w;
713   _ti->frame_height=_y4m->frame_h;
714   _ti->pic_width=_y4m->pic_w;
715   _ti->pic_height=_y4m->pic_h;
716   _ti->pic_x=_y4m->pic_x;
717   _ti->pic_y=_y4m->pic_y;
718   _ti->fps_numerator=_y4m->fps_n;
719   _ti->fps_denominator=_y4m->fps_d;
720   _ti->aspect_numerator=_y4m->par_n;
721   _ti->aspect_denominator=_y4m->par_d;
722   _ti->pixel_fmt=_y4m->dst_c_dec_h==2?
723    (_y4m->dst_c_dec_v==2?TH_PF_420:TH_PF_422):TH_PF_444;
724 }
725
726 static int y4m_input_fetch_frame(y4m_input *_y4m,FILE *_fin,
727  th_ycbcr_buffer _ycbcr){
728   char frame[6];
729   int  pic_sz;
730   int  frame_c_w;
731   int  frame_c_h;
732   int  c_w;
733   int  c_h;
734   int  c_sz;
735   int  ret;
736   pic_sz=_y4m->pic_w*_y4m->pic_h;
737   frame_c_w=_y4m->frame_w/_y4m->dst_c_dec_h;
738   frame_c_h=_y4m->frame_h/_y4m->dst_c_dec_v;
739   c_w=(_y4m->pic_w+_y4m->dst_c_dec_h-1)/_y4m->dst_c_dec_h;
740   c_h=(_y4m->pic_h+_y4m->dst_c_dec_v-1)/_y4m->dst_c_dec_v;
741   c_sz=c_w*c_h;
742   /*Read and skip the frame header.*/
743   ret=fread(frame,1,6,_fin);
744   if(ret<6)return 0;
745   if(memcmp(frame,"FRAME",5)){
746     fprintf(stderr,"Loss of framing in YUV input data\n");
747     exit(1);
748   }
749   if(frame[5]!='\n'){
750     char c;
751     int  j;
752     for(j=0;j<79&&fread(&c,1,1,_fin)&&c!='\n';j++);
753     if(j==79){
754       fprintf(stderr,"Error parsing YUV frame header\n");
755       return -1;
756     }
757   }
758   /*Read the frame data that needs no conversion.*/
759   if(fread(_y4m->dst_buf,1,_y4m->dst_buf_read_sz,_fin)!=_y4m->dst_buf_read_sz){
760     fprintf(stderr,"Error reading YUV frame data.\n");
761     return -1;
762   }
763   /*Read the frame data that does need conversion.*/
764   if(fread(_y4m->aux_buf,1,_y4m->aux_buf_read_sz,_fin)!=_y4m->aux_buf_read_sz){
765     fprintf(stderr,"Error reading YUV frame data.\n");
766     return -1;
767   }
768   /*Now convert the just read frame.*/
769   (*_y4m->convert)(_y4m,_y4m->dst_buf,_y4m->aux_buf);
770   /*Fill in the frame buffer pointers.*/
771   _ycbcr[0].width=_y4m->frame_w;
772   _ycbcr[0].height=_y4m->frame_h;
773   _ycbcr[0].stride=_y4m->pic_w;
774   _ycbcr[0].data=_y4m->dst_buf-_y4m->pic_x-_y4m->pic_y*_y4m->pic_w;
775   _ycbcr[1].width=frame_c_w;
776   _ycbcr[1].height=frame_c_h;
777   _ycbcr[1].stride=c_w;
778   _ycbcr[1].data=_y4m->dst_buf+pic_sz-(_y4m->pic_x/_y4m->dst_c_dec_h)-
779    (_y4m->pic_y/_y4m->dst_c_dec_v)*c_w;
780   _ycbcr[2].width=frame_c_w;
781   _ycbcr[2].height=frame_c_h;
782   _ycbcr[2].stride=c_w;
783   _ycbcr[2].data=_ycbcr[1].data+c_sz;
784   return 1;
785 }
786
787 static void y4m_input_close(y4m_input *_y4m){
788   free(_y4m->dst_buf);
789   free(_y4m->aux_buf);
790 }
791
792
793
794 typedef struct th_input th_input;
795
796 struct th_input{
797   ogg_sync_state    oy;
798   int               theora_p;
799   ogg_stream_state  to;
800   th_info           ti;
801   th_comment        tc;
802   th_dec_ctx       *td;
803 };
804
805
806
807 /*Grab some more compressed bitstream and sync it for page extraction.*/
808 static int th_input_buffer_data(th_input *_th,FILE *_fin){
809   char *buffer;
810   int bytes;
811   buffer=ogg_sync_buffer(&_th->oy,4096);
812   bytes=fread(buffer,1,4096,_fin);
813   ogg_sync_wrote(&_th->oy,bytes);
814   return bytes;
815 }
816
817 /*Push a page into the appropriate steam.
818   This can be done blindly; a stream won't accept a page that doesn't belong to
819    it.*/
820 static void th_input_queue_page(th_input *_th,ogg_page *_og){
821   if(_th->theora_p)ogg_stream_pagein(&_th->to,_og);
822 }
823
824 static int th_input_open_impl(th_input *_th,th_setup_info **_ts,FILE *_fin,
825  char *_sig,int _nsig){
826   ogg_packet op;
827   ogg_page   og;
828   int        nheaders_left;
829   int        done_headers;
830   ogg_sync_init(&_th->oy);
831   th_info_init(&_th->ti);
832   th_comment_init(&_th->tc);
833   *_ts=NULL;
834   /*Buffer any initial data read for file ID.*/
835   if(_nsig>0){
836     char *buffer;
837     buffer=ogg_sync_buffer(&_th->oy,_nsig);
838     memcpy(buffer,_sig,_nsig);
839     ogg_sync_wrote(&_th->oy,_nsig);
840   }
841   _th->theora_p=0;
842   nheaders_left=0;
843   for(done_headers=0;!done_headers;){
844     if(th_input_buffer_data(_th,_fin)==0)break;
845     while(ogg_sync_pageout(&_th->oy,&og)>0){
846       ogg_stream_state test;
847       /*Is this a mandated initial header?
848         If not, stop parsing.*/
849       if(!ogg_page_bos(&og)){
850         /*Don't leak the page; get it into the appropriate stream.*/
851         th_input_queue_page(_th,&og);
852         done_headers=1;
853         break;
854       }
855       ogg_stream_init(&test,ogg_page_serialno(&og));
856       ogg_stream_pagein(&test,&og);
857       ogg_stream_packetpeek(&test,&op);
858       /*Identify the codec: try Theora.*/
859       if(!_th->theora_p){
860         nheaders_left=th_decode_headerin(&_th->ti,&_th->tc,_ts,&op);
861         if(nheaders_left>=0){
862           /*It is Theora.*/
863           memcpy(&_th->to,&test,sizeof(test));
864           _th->theora_p=1;
865           /*Advance past the successfully processed header.*/
866           if(nheaders_left>0)ogg_stream_packetout(&_th->to,NULL);
867           continue;
868         }
869       }
870       /*Whatever it is, we don't care about it.*/
871       ogg_stream_clear(&test);
872     }
873   }
874   /*We're expecting more header packets.*/
875   while(_th->theora_p&&nheaders_left>0){
876     int ret;
877     while(nheaders_left>0){
878       ret=ogg_stream_packetpeek(&_th->to,&op);
879       if(ret==0)break;
880       if(ret<0)continue;
881       nheaders_left=th_decode_headerin(&_th->ti,&_th->tc,_ts,&op);
882       if(nheaders_left<0){
883         fprintf(stderr,"Error parsing Theora stream headers; "
884          "corrupt stream?\n");
885         return -1;
886       }
887       /*Advance past the successfully processed header.*/
888       else if(nheaders_left>0)ogg_stream_packetout(&_th->to,NULL);
889       _th->theora_p++;
890     }
891     /*Stop now so we don't fail if there aren't enough pages in a short
892        stream.*/
893     if(!(_th->theora_p&&nheaders_left>0))break;
894     /*The header pages/packets will arrive before anything else we care
895        about, or the stream is not obeying spec.*/
896     if(ogg_sync_pageout(&_th->oy,&og)>0)th_input_queue_page(_th,&og);
897     /*We need more data.*/
898     else if(th_input_buffer_data(_th,_fin)==0){
899       fprintf(stderr,"End of file while searching for codec headers.\n");
900       return -1;
901     }
902   }
903   /*And now we have it all.
904     Initialize the decoder.*/
905   if(_th->theora_p){
906     _th->td=th_decode_alloc(&_th->ti,*_ts);
907     if(_th->td!=NULL){
908       fprintf(stderr,"Ogg logical stream %lx is Theora %ix%i %.02f fps video.\n"
909        "Encoded frame content is %ix%i with %ix%i offset.\n",
910        _th->to.serialno,_th->ti.frame_width,_th->ti.frame_height,
911        (double)_th->ti.fps_numerator/_th->ti.fps_denominator,
912        _th->ti.pic_width,_th->ti.pic_height,_th->ti.pic_x,_th->ti.pic_y);
913       return 1;
914     }
915   }
916   return -1;
917 }
918
919 static void th_input_close(th_input *_th){
920   if(_th->theora_p){
921     ogg_stream_clear(&_th->to);
922     th_decode_free(_th->td);
923   }
924   th_comment_clear(&_th->tc);
925   th_info_clear(&_th->ti);
926   ogg_sync_clear(&_th->oy);
927 }
928
929 static int th_input_open(th_input *_th,FILE *_fin,char *_sig,int _nsig){
930   th_input       th;
931   th_setup_info *ts;
932   int            ret;
933   ret=th_input_open_impl(&th,&ts,_fin,_sig,_nsig);
934   th_setup_free(ts);
935   /*Clean up on failure.*/
936   if(ret<0)th_input_close(&th);
937   else memcpy(_th,&th,sizeof(th));
938   return ret;
939 }
940
941 static void th_input_get_info(th_input *_th,th_info *_ti){
942   memcpy(_ti,&_th->ti,sizeof(*_ti));
943 }
944
945 static int th_input_fetch_frame(th_input *_th,FILE *_fin,
946  th_ycbcr_buffer _ycbcr){
947   for(;;){
948     ogg_page   og;
949     ogg_packet op;
950     if(ogg_stream_packetout(&_th->to,&op)>0){
951       if(th_decode_packetin(_th->td,&op,NULL)>=0){
952         th_decode_ycbcr_out(_th->td,_ycbcr);
953         if(!summary_only&&show_frame_type){
954           printf("%c",th_packet_iskeyframe(&op)?'K':'D');
955           if(op.bytes>0)printf("%02i ",op.packet[0]&0x3F);
956           else printf("-- ");
957         }
958         return 1;
959       }
960       else return -1;
961     }
962     while(ogg_sync_pageout(&_th->oy,&og)<=0){
963       if(th_input_buffer_data(_th,_fin)==0)return feof(_fin)?0:-1;
964     }
965     th_input_queue_page(_th,&og);
966   }
967 }
968
969
970
971 typedef struct video_input video_input;
972 typedef void (*video_input_get_info_func)(void *_ctx,th_info *_ti);
973 typedef int (*video_input_fetch_frame_func)(void *_ctx,FILE *_fin,
974  th_ycbcr_buffer _ycbcr);
975 typedef void (*video_input_close_func)(void *_ctx);
976
977 struct video_input{
978   FILE                         *fin;
979   video_input_get_info_func     get_info;
980   video_input_fetch_frame_func  fetch_frame;
981   video_input_close_func        close;
982   union{
983     y4m_input y4m;
984     th_input  th;
985   }ctx;
986 };
987
988 static int video_input_open(video_input *_vid,FILE *_fin){
989   char buffer[4];
990   int  ret;
991   /* look for magic */
992   ret=fread(buffer,1,4,_fin);
993   if(ret<4)fprintf(stderr,"EOF determining file type of file.\n");
994   else{
995     if(!memcmp(buffer,"YUV4",4)){
996       if(y4m_input_open(&_vid->ctx.y4m,_fin,buffer,4)>=0){
997         /*fprintf(stderr,"Original %s is %dx%d %.02f fps %s video.\n",
998          f,_y4m->pic_w,_y4m->pic_h,(double)_y4m->fps_n/_y4m->fps_d,_y4m->chroma_type);*/
999         _vid->fin=_fin;
1000         _vid->get_info=(video_input_get_info_func)y4m_input_get_info;
1001         _vid->fetch_frame=(video_input_fetch_frame_func)y4m_input_fetch_frame;
1002         _vid->close=(video_input_close_func)y4m_input_close;
1003         return 0;
1004       }
1005     }
1006     else if(!memcmp(buffer,"OggS",4)){
1007       if(th_input_open(&_vid->ctx.th,_fin,buffer,4)>=0){
1008         _vid->fin=_fin;
1009         _vid->get_info=(video_input_get_info_func)th_input_get_info;
1010         _vid->fetch_frame=(video_input_fetch_frame_func)th_input_fetch_frame;
1011         _vid->close=(video_input_close_func)th_input_close;
1012         return 0;
1013       }
1014     }
1015     else fprintf(stderr,"Unknown file type.\n");
1016   }
1017   return -1;
1018 }
1019
1020 static void video_input_get_info(video_input *_vid,th_info *_ti){
1021   (*_vid->get_info)(&_vid->ctx,_ti);
1022 }
1023
1024 static int video_input_fetch_frame(video_input *_vid,th_ycbcr_buffer _ycbcr){
1025   return (*_vid->fetch_frame)(&_vid->ctx,_vid->fin,_ycbcr);
1026 }
1027
1028 static void video_input_close(video_input *_vid){
1029   (*_vid->close)(&_vid->ctx);
1030   fclose(_vid->fin);
1031 }
1032
1033
1034
1035 static void usage(char *_argv[]){
1036   fprintf(stderr,"Usage: %s [options] <video1> <video2>\n"
1037    "    <video1> and <video1> may be either YUV4MPEG or Ogg Theora files.\n\n"
1038    "    Options:\n\n"
1039    "      -f --frame-type Show frame type and QI value for each Theora frame.\n"
1040    "      -s --summary    Only output the summary line.\n"
1041    "      -y --luma-only  Only output values for the luma channel.\n",_argv[0]);
1042 }
1043
1044 int main(int _argc,char *_argv[]){
1045   video_input  vid1;
1046   th_info      ti1;
1047   video_input  vid2;
1048   th_info      ti2;
1049   ogg_int64_t  gsqerr;
1050   ogg_int64_t  gnpixels;
1051   ogg_int64_t  gplsqerr[3];
1052   ogg_int64_t  gplnpixels[3];
1053   int          frameno;
1054   FILE        *fin;
1055   int          long_option_index;
1056   int          c;
1057 #ifdef _WIN32
1058   /*We need to set stdin/stdout to binary mode on windows.
1059     Beware the evil ifdef.
1060     We avoid these where we can, but this one we cannot.
1061     Don't add any more, you'll probably go to hell if you do.*/
1062   _setmode(_fileno(stdin),_O_BINARY);
1063 #endif
1064   /*Process option arguments.*/
1065   while((c=getopt_long(_argc,_argv,optstring,options,&long_option_index))!=EOF){
1066     switch(c){
1067       case 'f':show_frame_type=1;break;
1068       case 's':summary_only=1;break;
1069       case 'y':luma_only=1;break;
1070       default:usage(_argv);break;
1071     }
1072   }
1073   if(optind+2!=_argc){
1074     usage(_argv);
1075     exit(1);
1076   }
1077   fin=strcmp(_argv[optind],"-")==0?stdin:fopen(_argv[optind],"rb");
1078   if(fin==NULL){
1079     fprintf(stderr,"Unable to open '%s' for extraction.\n",_argv[optind]);
1080     exit(1);
1081   }
1082   fprintf(stderr,"Opening %s...\n",_argv[optind]);
1083   if(video_input_open(&vid1,fin)<0)exit(1);
1084   video_input_get_info(&vid1,&ti1);
1085   fin=strcmp(_argv[optind+1],"-")==0?stdin:fopen(_argv[optind+1],"rb");
1086   if(fin==NULL){
1087     fprintf(stderr,"Unable to open '%s' for extraction.\n",_argv[optind+1]);
1088     exit(1);
1089   }
1090   fprintf(stderr,"Opening %s...\n",_argv[optind+1]);
1091   if(video_input_open(&vid2,fin)<0)exit(1);
1092   video_input_get_info(&vid2,&ti2);
1093   /*Check to make sure these videos are compatible.*/
1094   if(ti1.pic_width!=ti2.pic_width||ti1.pic_height!=ti2.pic_height){
1095     fprintf(stderr,"Video resolution does not match.\n");
1096     exit(1);
1097   }
1098   if(ti1.pixel_fmt!=ti2.pixel_fmt){
1099     fprintf(stderr,"Pixel formats do not match.\n");
1100     exit(1);
1101   }
1102   if((ti1.pic_x&!(ti1.pixel_fmt&1))!=(ti2.pic_x&!(ti2.pixel_fmt&1))||
1103    (ti1.pic_y&!(ti1.pixel_fmt&2))!=(ti2.pic_y&!(ti2.pixel_fmt&2))){
1104     fprintf(stderr,"Chroma subsampling offsets do not match.\n");
1105     exit(1);
1106   }
1107   if(ti1.fps_numerator*(ogg_int64_t)ti2.fps_denominator!=
1108    ti2.fps_numerator*(ogg_int64_t)ti1.fps_denominator){
1109     fprintf(stderr,"Warning: framerates do not match.\n");
1110   }
1111   if(ti1.aspect_numerator*(ogg_int64_t)ti2.aspect_denominator!=
1112    ti2.aspect_numerator*(ogg_int64_t)ti1.aspect_denominator){
1113     fprintf(stderr,"Warning: aspect ratios do not match.\n");
1114   }
1115   gsqerr=gplsqerr[0]=gplsqerr[1]=gplsqerr[2]=0;
1116   gnpixels=gplnpixels[0]=gplnpixels[1]=gplnpixels[2]=0;
1117   for(frameno=0;;frameno++){
1118     th_ycbcr_buffer f1;
1119     th_ycbcr_buffer f2;
1120     ogg_int64_t     plsqerr[3];
1121     long            plnpixels[3];
1122     ogg_int64_t     sqerr;
1123     long            npixels;
1124     int             ret1;
1125     int             ret2;
1126     int             pli;
1127     ret1=video_input_fetch_frame(&vid1,f1);
1128     ret2=video_input_fetch_frame(&vid2,f2);
1129     if(ret1==0&&ret2==0)break;
1130     else if(ret1<0||ret2<0)break;
1131     else if(ret1==0){
1132       fprintf(stderr,"%s ended before %s.\n",
1133        _argv[optind],_argv[optind+1]);
1134       break;
1135     }
1136     else if(ret2==0){
1137       fprintf(stderr,"%s ended before %s.\n",
1138        _argv[optind+1],_argv[optind]);
1139       break;
1140     }
1141     /*Okay, we got one frame from each.*/
1142     sqerr=0;
1143     npixels=0;
1144     for(pli=0;pli<3;pli++){
1145       int xdec;
1146       int ydec;
1147       int y1;
1148       int y2;
1149       xdec=pli&&!(ti1.pixel_fmt&1);
1150       ydec=pli&&!(ti1.pixel_fmt&2);
1151       plsqerr[pli]=0;
1152       plnpixels[pli]=0;
1153       for(y1=ti1.pic_y>>ydec,y2=ti2.pic_y>>ydec;
1154        y1<ti1.pic_y+ti1.pic_height+ydec>>ydec;y1++,y2++){
1155         int x1;
1156         int x2;
1157         for(x1=ti1.pic_x>>xdec,x2=ti2.pic_x>>xdec;
1158          x1<ti1.pic_x+ti1.pic_width+xdec>>xdec;x1++,x2++){
1159           int d;
1160           d=*(f1[pli].data+y1*f1[pli].stride+x1)-
1161            *(f2[pli].data+y2*f2[pli].stride+x2);
1162           plsqerr[pli]+=d*d;
1163           plnpixels[pli]++;
1164         }
1165       }
1166       sqerr+=plsqerr[pli];
1167       gplsqerr[pli]+=plsqerr[pli];
1168       npixels+=plnpixels[pli];
1169       gplnpixels[pli]+=plnpixels[pli];
1170     }
1171     if(!summary_only){
1172       if(!luma_only){
1173         printf("%08i: %-7lG  (Y': %-7lG  Cb: %-7lG  Cr: %-7lG)\n",frameno,
1174          10*(log10(255*255)+log10(npixels)-log10(sqerr)),
1175          10*(log10(255*255)+log10(plnpixels[0])-log10(plsqerr[0])),
1176          10*(log10(255*255)+log10(plnpixels[1])-log10(plsqerr[1])),
1177          10*(log10(255*255)+log10(plnpixels[2])-log10(plsqerr[2])));
1178       }
1179       else{
1180         printf("%08i: %-7lG\n",frameno,
1181          10*(log10(255*255)+log10(plnpixels[0])-log10(plsqerr[0])));
1182       }
1183     }
1184     gsqerr+=sqerr;
1185     gnpixels+=npixels;
1186   }
1187   if(!luma_only){
1188     printf("Total: %-7lG  (Y': %-7lG  Cb: %-7lG  Cr: %-7lG)\n",
1189      10*(log10(255*255)+log10(gnpixels)-log10(gsqerr)),
1190      10*(log10(255*255)+log10(gplnpixels[0])-log10(gplsqerr[0])),
1191      10*(log10(255*255)+log10(gplnpixels[1])-log10(gplsqerr[1])),
1192      10*(log10(255*255)+log10(gplnpixels[2])-log10(gplsqerr[2])));
1193   }
1194   else{
1195     printf("Total: %-7lG\n",
1196      10*(log10(255*255)+log10(gplnpixels[0])-log10(gplsqerr[0])));
1197   }
1198   video_input_close(&vid1);
1199   video_input_close(&vid2);
1200   return 0;
1201 }