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