--- /dev/null
+#include<stdio.h>
+#include<stdarg.h>
+#include<stdint.h>
+#include<stdlib.h>
+#include<string.h>
+#include<unistd.h>
+#include<math.h>
+
+void write_pbm(uint8_t *tp, int w, int h, const char *fmt, ...)
+{
+ va_list ap; va_start(ap, fmt);
+ char fn[256]; vsnprintf(fn, sizeof(fn), fmt, ap);
+ va_end(ap);
+ FILE *fp = !strcmp(fn,"-") ? stdout : fopen(fn,"w");
+ if( fp ) {
+ fprintf(fp,"P5\n%d %d\n255\n",w,h);
+ fwrite(tp,w,h,fp);
+ fclose(fp);
+ }
+}
+
+class Scale {
+ int iw, ih;
+ double rx, ry, rw, rh;
+ int ix0, ix1, iy0, iy1, rsh;
+ double xf0, xf1, yf0, yf1, pw;
+ uint8_t *idat, *odat;
+ uint8_t *row(int y) { return idat + (iw<<rsh) * (y>>rsh); }
+
+ inline double scalex(uint8_t *rp, double v, double yw) {
+ int x = ix0; rp += x;
+ v += yw*xf0 * *rp++;
+ while( ++x < ix1 ) v += yw * *rp++;
+ v += yw*xf1 * *rp++;
+ return v;
+ }
+
+ inline double scalexy() {
+ double yw, v = 0;
+ int y = iy0;
+ if( (yw = pw * yf0) > 0 ) v = scalex(row(y), v, yw);
+ while( ++y < iy1 ) v = scalex(row(y), v, pw);
+ if( (yw = pw * yf1) > 0 ) v = scalex(row(y), v, yw);
+ return v;
+ }
+public:
+ Scale(uint8_t *idata, int type, int iw, int ih,
+ double rx, double ry, double rw, double rh) {
+ rsh = type ? 1 : 0;
+ this->idat = idata;
+ this->iw = iw; this->ih = ih;
+ this->rx = rx; this->ry = ry;
+ this->rw = rw; this->rh = rh;
+ }
+ void scale(uint8_t *odata, int ow, int oh, int sx, int sy, int sw, int sh);
+};
+
+void Scale::
+scale(uint8_t *odata, int ow, int oh, int sx, int sy, int sw, int sh)
+{
+ pw = (double)(ow * oh) / (rw * rh);
+ odat = odata + sy*ow + sx;
+ double ly = 0.0;
+
+ for( int dy=1; dy<=sh; ++dy, odat+=ow ) {
+ iy0 = (int) ly; yf0 = 1.0 - (ly-iy0);
+ double ny = (dy * rh) / oh + ry;
+ iy1 = (int) ny; yf1 = ny-iy1;
+ uint8_t *bp = odat;
+ double lx = rx;
+
+ for( int dx=1; dx<=sw; ++dx ) {
+ ix0 = (int)lx; xf0 = 1.0 - (lx-ix0);
+ double nx = (dx*rw) / ow + rx;
+ ix1 = (int)nx; xf1 = nx-ix1;
+ *bp++ = scalexy();
+ xf0 = 1.0-xf1; lx = nx;
+ }
+ yf0 = 1.0 - yf1; ly = ny;
+ }
+}
+
+
+double mean(uint8_t *dat, int n, int stride)
+{
+ int s = 0;
+ for( int i=0; i<n; ++i, dat+=stride ) s += *dat;
+ return (double)s / n;
+}
+
+double variance(uint8_t *dat, double m, int n, int stride)
+{
+ double ss = 0, s = 0;
+ for( int i=0; i<n; ++i, dat+=stride ) {
+ double dx = *dat-m;
+ s += dx; ss += dx*dx;
+ }
+ return (ss - s*s / n) / (n-1);
+}
+
+double std_dev(uint8_t *dat, double m, int n, int stride)
+{
+ return sqrt(variance(dat, m, n, stride));
+}
+
+double centroid(uint8_t *dat, int n, int stride)
+{
+ int s = 0, ss = 0;
+ for( int i=0; i<n; dat+=stride ) { s += *dat; ss += ++i * *dat; }
+ return s > 0 ? (double)ss / s : n / 2.;
+}
+
+void deviation(uint8_t *a, uint8_t *b, int sz, int margin, int stride, double &dev, int &ofs)
+{
+ double best = 1e100; int ii = -1;
+ for( int i=-margin; i<=margin; ++i ) {
+ int aofs = i < 0 ? 0 : i*stride;
+ int bofs = i < 0 ? -i*stride : 0;
+ uint8_t *ap = a + aofs, *bp = b + bofs;
+ int ierr = 0; int n = sz - abs(i);
+ for( int j=n; --j>=0; ap+=stride,bp+=stride ) ierr += abs(*ap - *bp);
+ double err = (double)ierr / n;
+ if( err < best ) { best = err; ii = i; }
+ }
+ dev = best;
+ ofs = ii;
+}
+
+void diff(uint8_t *a, uint8_t *b, int w, int h, int dx, int dy, int bias, double &dev)
+{
+ int axofs = dx < 0 ? 0 : dx;
+ int ayofs = w * (dy < 0 ? 0 : dy);
+ int aofs = ayofs + axofs;
+ int bxofs = dx < 0 ? -dx : 0;
+ int byofs = w * (dy < 0 ? -dy : 0);
+ int bofs = byofs + bxofs;
+ int sz = w * h;
+ uint8_t *ap = a + aofs, *bp = b + bofs;
+ int ierr = 0, ww = w-abs(dx), hh = h-abs(dy), ssz = ww * hh;
+uint8_t dif[ssz], *dp = dif;
+ for( int y=hh; --y>=0; ap+=w, bp+=w ) {
+ a = ap; b = bp;
+ for( int x=ww; --x>=0; ++a, ++b ) { ierr += abs(*a-bias - *b); *dp++ = *a-bias-*b+128; }
+ }
+write_pbm(dif, ww, hh, "/tmp/x");
+ dev = (double)ierr / ssz;
+}
+
+
+int main(int ac, char **av)
+{
+ char line[120];
+ FILE *afp = fopen(av[1],"r");
+ FILE *bfp = fopen(av[2],"r");
+ fgets(line,sizeof(line),afp);
+ fgets(line,sizeof(line),afp);
+ int aw, ah; sscanf(line,"%d %d",&aw, &ah);
+ fgets(line,sizeof(line),afp);
+ fgets(line,sizeof(line),bfp);
+ fgets(line,sizeof(line),bfp);
+ int bw, bh; sscanf(line,"%d %d",&bw, &bh);
+ if( aw != bw || ah != bh ) exit(1);
+ int w = aw, h = ah, sz = w * h;
+ fgets(line,sizeof(line),bfp);
+ uint8_t a[sz], b[sz];
+ uint8_t *ap = a, *bp = b;
+ int ach, bch, n = 0;
+ while( (ach=getc(afp))>=0 && (bch=getc(bfp))>=0 ) {
+ *ap++ = ach; *bp++ = bch; if( ++n >= sz ) break;
+ }
+ double dev;
+ for( int dy=-2; dy<=2; ++dy ) {
+ printf(" dy=%d: ", dy);
+ for( int dx=-2; dx<=2; ++dx ) {
+ diff(a, b, w, h, dx, dy, 0, dev),
+ printf(" dx=%d/dev %f ", dx, dev);
+ }
+ printf("\n");
+ }
+ uint8_t *ax = a + sz/2, *ay = a + w/2;
+ uint8_t *bx = b + sz/2, *by = b + w/2;
+ double xam = mean(ax, w, 1), xbm = mean(bx, w, 1);
+ double yam = mean(ay, h, w), ybm = mean(by, h, w);
+ double am = (xam + yam) / 2, bm = (xbm + ybm) / 2;
+ int bias = am-bm;
+ printf("xam/yam %f/%f, xbm/ybm %f/%f, am/bm %f/%f, bias %d\n",xam,yam,xbm,ybm,am,bm,bias);
+ double xasd = std_dev(ax, xam, w, 1), xbsd = std_dev(bx, xbm, w, 1);
+ double yasd = std_dev(ay, yam, h, w), ybsd = std_dev(by, ybm, h, w);
+ double asd = (xasd + yasd) / 2, bsd = (xbsd + ybsd) / 2;
+ printf("xasd/yasd %f/%f, xbsd/ybsd %f/%f, asd/bsd %f/%f\n",xasd,yasd,xbsd,ybsd,asd,bsd);
+ double xdev, ydev; int xofs, yofs;
+ deviation(ax, bx, w, 5, 1, xdev, xofs);
+ deviation(ay, by, h, 5, w, ydev, yofs);
+ diff(a, b, w, h, xofs, yofs, bias, dev);
+ printf("xdev/ofs %f/%d, ydev/ofs %f/%d, dev %f\n",xdev,xofs,ydev,yofs,dev);
+
+ double cax = centroid(ax, w, 1), cay = centroid(ay, h, w);
+ double cbx = centroid(bx, w, 1), cby = centroid(by, h, w);
+ double dx = cax - cbx, dy = cay - cby;
+ printf("cax/cay %f/%f, cbx/cby %f/%f, dx/dy %f/%f", cax, cay, cbx, cby, dx, dy);
+ int ww = w+7, hh = h+7; uint8_t aa[ww*hh]; memset(aa,0,ww*hh);
+ ap = a; bp = aa + 3*ww + 3; uint8_t c[sz];
+ for( int y=0; y<h; ++y, bp+=7 ) for( int x=0; x<w; ++x ) *bp++ = *ap++;
+ Scale(aa,0, ww,hh, 3+dx,3+dy,w,h).scale(a, w,h, 0,0,w,h);
+write_pbm(a,w,h,"/tmp/a");
+ diff(a, b, w, h, 0, 0, 0, dev);
+ printf(" dev %f\n", dev);
+ return 0;
+}
+
+
projects / goodguy / history.git / blobdiff