version update
[goodguy/cinelerra.git] / cinelerra-5.1 / db / utils / x.C
1 #include<stdio.h>
2 #include<stdarg.h>
3 #include<stdint.h>
4 #include<stdlib.h>
5 #include<string.h>
6 #include<unistd.h>
7 #include<math.h>
8
9 void write_pbm(uint8_t *tp, int w, int h, const char *fmt, ...)
10 {
11   va_list ap;    va_start(ap, fmt);
12   char fn[256];  vsnprintf(fn, sizeof(fn), fmt, ap);
13   va_end(ap);
14   FILE *fp = !strcmp(fn,"-") ? stdout : fopen(fn,"w");
15   if( fp ) {
16     fprintf(fp,"P5\n%d %d\n255\n",w,h);
17     fwrite(tp,w,h,fp);
18     fclose(fp);
19   }
20 }
21
22 class Scale {
23   int iw, ih;
24   double rx, ry, rw, rh;
25   int ix0, ix1, iy0, iy1, rsh;
26   double xf0, xf1, yf0, yf1, pw;
27   uint8_t *idat, *odat;
28   uint8_t *row(int y) { return idat + (iw<<rsh) * (y>>rsh); }
29
30   inline double scalex(uint8_t *rp, double v, double yw) {
31     int x = ix0;  rp += x;
32     v += yw*xf0 * *rp++;
33     while( ++x < ix1 ) v += yw * *rp++;
34     v += yw*xf1 * *rp++;
35     return v;
36   }
37
38   inline double scalexy() {
39     double yw, v = 0;
40     int y = iy0;
41     if( (yw = pw * yf0) > 0 ) v = scalex(row(y), v, yw);
42     while( ++y < iy1 ) v = scalex(row(y), v, pw);
43     if( (yw = pw * yf1) > 0 ) v = scalex(row(y), v, yw);
44     return v;
45   }
46 public:
47   Scale(uint8_t *idata, int type, int iw, int ih,
48                 double rx, double ry, double rw, double rh) {
49     rsh = type ? 1 : 0;
50     this->idat = idata;
51     this->iw = iw;  this->ih = ih;
52     this->rx = rx;  this->ry = ry;
53     this->rw = rw;  this->rh = rh;
54   }
55   void scale(uint8_t *odata, int ow, int oh, int sx, int sy, int sw, int sh);
56 };
57
58 void Scale::
59 scale(uint8_t *odata, int ow, int oh, int sx, int sy, int sw, int sh)
60 {
61   pw = (double)(ow * oh) / (rw * rh);
62   odat = odata + sy*ow + sx;
63   double ly = 0.0;
64
65   for( int dy=1; dy<=sh; ++dy, odat+=ow ) {
66     iy0 = (int) ly;  yf0 = 1.0 - (ly-iy0);
67     double ny = (dy * rh) / oh + ry;
68     iy1 = (int) ny;  yf1 = ny-iy1;
69     uint8_t *bp = odat;
70     double lx = rx;
71
72     for( int dx=1; dx<=sw; ++dx ) {
73       ix0 = (int)lx;  xf0 = 1.0 - (lx-ix0);
74       double nx = (dx*rw) / ow + rx;
75       ix1 = (int)nx;  xf1 = nx-ix1;
76       *bp++ = scalexy();
77       xf0 = 1.0-xf1;  lx = nx;
78     }
79     yf0 = 1.0 - yf1;  ly = ny;
80   }
81 }
82
83
84 double mean(uint8_t *dat, int n, int stride)
85 {
86   int s = 0;
87   for( int i=0; i<n; ++i, dat+=stride ) s += *dat;
88   return (double)s / n;
89 }
90
91 double variance(uint8_t *dat, double m, int n, int stride)
92 {
93   double ss = 0, s = 0;
94   for( int i=0; i<n; ++i, dat+=stride ) {
95     double dx = *dat-m;
96     s += dx;  ss += dx*dx;
97   }
98   return (ss - s*s / n) / (n-1);
99 }
100
101 double std_dev(uint8_t *dat, double m, int n, int stride)
102 {
103         return sqrt(variance(dat, m, n, stride));
104 }
105
106 double centroid(uint8_t *dat, int n, int stride)
107 {
108   int s = 0, ss = 0;
109   for( int i=0; i<n; dat+=stride ) { s += *dat;  ss += ++i * *dat; }
110   return s > 0 ? (double)ss / s : n / 2.;
111 }
112
113 void deviation(uint8_t *a, uint8_t *b, int sz, int margin, int stride, double &dev, int &ofs)
114 {
115   double best = 1e100;  int ii = -1;
116   for( int i=-margin; i<=margin; ++i ) {
117     int aofs = i < 0 ? 0 : i*stride;
118     int bofs = i < 0 ? -i*stride : 0;
119     uint8_t *ap = a + aofs, *bp = b + bofs;
120     int ierr = 0; int n = sz - abs(i);
121     for( int j=n; --j>=0; ap+=stride,bp+=stride ) ierr += abs(*ap - *bp);
122     double err = (double)ierr / n;
123     if( err < best ) { best = err;  ii = i; }
124   }
125   dev = best;
126   ofs = ii;
127 }
128
129 void diff(uint8_t *a, uint8_t *b, int w, int h, int dx, int dy, int bias, double &dev)
130 {
131   int axofs = dx < 0 ? 0 : dx;
132   int ayofs = w * (dy < 0 ? 0 : dy);
133   int aofs = ayofs + axofs;
134   int bxofs = dx < 0 ? -dx : 0;
135   int byofs = w * (dy < 0 ? -dy : 0);
136   int bofs = byofs + bxofs;
137   int sz = w * h;
138   uint8_t *ap = a + aofs, *bp = b + bofs;
139   int ierr = 0, ww = w-abs(dx), hh = h-abs(dy), ssz = ww * hh;
140 uint8_t dif[ssz], *dp = dif;
141   for( int y=hh; --y>=0; ap+=w, bp+=w ) {
142     a = ap;  b = bp;
143     for( int x=ww; --x>=0; ++a, ++b ) { ierr += abs(*a-bias - *b); *dp++ = *a-bias-*b+128; }
144   }
145 write_pbm(dif, ww, hh, "/tmp/x");
146   dev = (double)ierr / ssz;
147 }
148
149
150 int main(int ac, char **av)
151 {
152   char line[120];
153   FILE *afp = fopen(av[1],"r");
154   FILE *bfp = fopen(av[2],"r");
155   fgets(line,sizeof(line),afp);
156   fgets(line,sizeof(line),afp);
157   int aw, ah;  sscanf(line,"%d %d",&aw, &ah);
158   fgets(line,sizeof(line),afp);
159   fgets(line,sizeof(line),bfp);
160   fgets(line,sizeof(line),bfp);
161   int bw, bh;  sscanf(line,"%d %d",&bw, &bh);
162   if( aw != bw || ah != bh ) exit(1);
163   int w = aw, h = ah, sz = w * h;
164   fgets(line,sizeof(line),bfp);
165   uint8_t a[sz], b[sz];
166   uint8_t *ap = a, *bp = b;
167   int ach, bch, n = 0;
168   while( (ach=getc(afp))>=0 && (bch=getc(bfp))>=0 ) {
169     *ap++ = ach;  *bp++ = bch;  if( ++n >= sz ) break;
170   }
171   double dev;
172   for( int dy=-2; dy<=2; ++dy ) {
173     printf(" dy=%d: ", dy);
174     for( int dx=-2; dx<=2; ++dx ) {
175       diff(a, b, w, h, dx, dy, 0, dev),
176       printf(" dx=%d/dev %f ", dx, dev);
177     }
178     printf("\n");
179   }
180   uint8_t *ax = a + sz/2, *ay = a + w/2;
181   uint8_t *bx = b + sz/2, *by = b + w/2;
182   double xam = mean(ax, w, 1), xbm = mean(bx, w, 1);
183   double yam = mean(ay, h, w), ybm = mean(by, h, w);
184   double am = (xam + yam) / 2, bm = (xbm + ybm) / 2;
185   int bias = am-bm;
186   printf("xam/yam %f/%f, xbm/ybm %f/%f, am/bm %f/%f, bias %d\n",xam,yam,xbm,ybm,am,bm,bias);
187   double xasd = std_dev(ax, xam, w, 1), xbsd = std_dev(bx, xbm, w, 1);
188   double yasd = std_dev(ay, yam, h, w), ybsd = std_dev(by, ybm, h, w);
189   double asd = (xasd + yasd) / 2, bsd = (xbsd + ybsd) / 2;
190   printf("xasd/yasd %f/%f, xbsd/ybsd %f/%f, asd/bsd %f/%f\n",xasd,yasd,xbsd,ybsd,asd,bsd);
191   double xdev, ydev;  int xofs, yofs;
192   deviation(ax, bx, w, 5, 1, xdev, xofs);
193   deviation(ay, by, h, 5, w, ydev, yofs);
194   diff(a, b, w, h, xofs, yofs, bias, dev);
195   printf("xdev/ofs %f/%d, ydev/ofs %f/%d, dev %f\n",xdev,xofs,ydev,yofs,dev);
196
197   double cax = centroid(ax, w, 1),  cay = centroid(ay, h, w);
198   double cbx = centroid(bx, w, 1),  cby = centroid(by, h, w);
199   double dx = cax - cbx, dy = cay - cby;
200   printf("cax/cay %f/%f,  cbx/cby %f/%f, dx/dy %f/%f", cax, cay, cbx, cby, dx, dy);
201   int ww = w+7, hh = h+7;  uint8_t aa[ww*hh];  memset(aa,0,ww*hh);
202   ap = a;  bp = aa + 3*ww + 3;  uint8_t c[sz];
203   for( int y=0; y<h; ++y, bp+=7 ) for( int x=0; x<w; ++x ) *bp++ = *ap++;
204   Scale(aa,0, ww,hh, 3+dx,3+dy,w,h).scale(a, w,h, 0,0,w,h);
205 write_pbm(a,w,h,"/tmp/a");
206   diff(a, b, w, h, 0, 0, 0, dev);
207   printf(" dev %f\n", dev);
208   return 0;
209 }
210
211