/*
 *	tiffcine -- recognize and crop out frames from
 *		a scan of a strip of movie film
 *
 *	Copyright 2002 Richard J Kinch
 *
 *	RJK March 2002
 */

#include <stdio.h>
#include <math.h>
#if FFTW
#include <fftw.h>
#endif
#include "tiffio.h"

#define mag(c) sqrt((c).re*(c).re+(c).im*(c).im)
#define ang(c) (fabs((c).re)>=1e-6?atan((c).im/(c).re):0)
#define deg(r) ((r)*180.0/3.1415926)

/* Greyscale value for a given RGB raster word				*/
#define greyscale(v) (TIFFGetR(v)+TIFFGetG(v)+TIFFGetB(v))
#define threshold(v) ((greyscale(v)>382))
#define greyscale_white 255*3
#define greyscale_black 0

/*
  SOME DIMENSIONS AT 600 PIXELS/INCH FOR 16 MM FILM:

  Frame vertical interval is 180 pixels
  Film total width is 377 pixels
  Clear image width between sprockets 248 pixels
    Center-to-center of sprockets 291 pixels
  Sprocket height 0.050in +/- 0.0004in, i.e., 30 pixels +/- 0.24 pixel
  16 mm full frame at 600 dpi is 377 x 180 pixels
  Image059.tif sample image frames at X=62, with lower-left sprocket centered at Y=39
  Image011.tif sample image frames at X=57, with lower-left sprocket centered at Y=39

  #define START_OF_FIRST_FRAME 57
 
*/

/* 16mm */
#define FILM_FRAME_HEIGHT 180
#define FILM_PERFORATION_HEIGHT 30	/* Must be < FILM_FRAME_HEIGHT */
#define CROPPED_FRAME_HEIGHT 180
#define CROPPED_FRAME_WIDTH 291
#define CROPPED_FRAME_LEFT 39
#define LEFT_SPROCKET_CENTERLINE 40	/* image058.tif			*/
#define RIGHT_SPROCKET_CENTERLINE 331	/* image058.tif			*/

/* Number of text columns to use to plot greyscale values		*/
#define PLOT_COLUMNS 75

void
crop(uint32 *raster,uint32 w,uint32 h,
	float xres, float yres, uint16 resunits,
	uint32 x0, uint32 y0, uint32 xN, uint32 yN, char *fname,
	int rotate) {
    /* Given an input RGB raster of dimensions (w,h), write a new TIFF	*/
    /* file named fname, the new image being cropped from the input	*/
    /* raster starting at (x0,y0) and having size (xN,yN).  If "rotate"	*/
    /* we apply a 90 degree clockwise rotation to the output file.	*/
    TIFF *out;
    unsigned char *buf, *p;
    uint32 *q;
    uint32 row, i;
    uint32 cw, ch;		/* Cropped width/height after rotation	*/
    tsize_t scanlinesize;
    if (w==0 || h==0) return;
    if (x0>=w || y0>=h || (x0+xN)>w || (y0+yN)>h) return;
    if (fopen(fname,"rb")) {
	fprintf(stderr,"File \"%s\" already exists\n",fname);
	return;
	}
    out = TIFFOpen(fname,"w");
    if (out==NULL) return;
    if (rotate) { cw = yN ; ch = xN ; } else { cw = xN ; ch = yN ; }
    TIFFSetField(out,TIFFTAG_IMAGEWIDTH,cw);
    TIFFSetField(out,TIFFTAG_IMAGELENGTH,ch);
    TIFFSetField(out,TIFFTAG_PHOTOMETRIC,((uint16)PHOTOMETRIC_RGB));
    TIFFSetField(out,TIFFTAG_COMPRESSION,((uint16)COMPRESSION_NONE));
    TIFFSetField(out,TIFFTAG_SAMPLESPERPIXEL,((uint16)3));
    TIFFSetField(out,TIFFTAG_PLANARCONFIG,((uint16)PLANARCONFIG_CONTIG));
    TIFFSetField(out,TIFFTAG_BITSPERSAMPLE,((uint16)8));
    TIFFSetField(out,TIFFTAG_RESOLUTIONUNIT,(uint16)resunits);
    TIFFSetField(out,TIFFTAG_XRESOLUTION,(float)xres);
    TIFFSetField(out,TIFFTAG_YRESOLUTION,(float)yres);
    scanlinesize =  TIFFScanlineSize(out);
    if (scanlinesize<=0) return;
    buf = (unsigned char *) _TIFFmalloc(scanlinesize);
    if (buf==NULL) return;
    if (rotate) for (row=0; row<ch; row++) {
	p = buf;
	/* q = raster+(x0+ch-1-row)+(y0*w); */
	q = raster+(x0+row)+(y0*w);
	for (i=0; i<cw; i++) {
	    *p++ = TIFFGetR(*q);
	    *p++ = TIFFGetG(*q);
	    *p++ = TIFFGetB(*q);
	    q += w;
	    }
	TIFFWriteScanline(out,buf,row,0);
	}
    else for (row=0; row<ch; row++) {
	p = buf;
	q = raster+w*(y0+ch-1-row)+x0;
	for (i=0; i<cw; i++,q++) {
	    *p++ = TIFFGetR(*q);
	    *p++ = TIFFGetG(*q);
	    *p++ = TIFFGetB(*q);
	    }
	TIFFWriteScanline(out,buf,row,0);
	}
    TIFFClose(out);
    }

#if FFTW
void
plot(fftw_complex *Z,int N,double scale) {
    int i,j;
    double v;
    char *s;
    printf("\nN=%d\n\n",N);
    for (i=0; i<N; i++) if (mag(Z[i])) {
	v = mag(Z[i])*scale;
	printf("%4d %5.1f %6.2f ",i,v,deg(ang(Z[i])));
	if (v>PLOT_COLUMNS) { v = PLOT_COLUMNS; s = "...\n"; }
	else s = "\n";
	for (j=0; j<v; j++) putchar('*');
	printf(s);
	}
    }

/* 16mm film ideal perforation scanline, and its 1D DFT			*/
fftw_complex hPerf[FILM_FRAME_HEIGHT], HPerf[FILM_FRAME_HEIGHT];
/* Raster scanline workspace and its 1D DFT				*/
fftw_complex hScan[FILM_FRAME_HEIGHT], HScan[FILM_FRAME_HEIGHT];
/* Correlation workspace and its 1D DFT					*/
fftw_complex hCorr[FILM_FRAME_HEIGHT], HCorr[FILM_FRAME_HEIGHT];

/* FFT and IFFT plans for analyzine a full frame scanline		*/
fftw_plan FFT, IFFT;

void
fft_correlate(uint32 *raster,uint32 w,uint32 h,int N,int Y,fftw_complex *H) {
    /* Using a frame height of N, find offset that best correlates	*/
    /* a scanline of the raster to an ideal pattern, where the ideal	*/
    /* is given as H, namely the complex conjugate of the DFT of the	*/
    /* ideal.  Raster has width w and height h in pixels.  The		*/
    /* is found by correlating raster scanline Y to the ideal		*/
    /* scanline, which is performed by vector multiplying the DFT of	*/
    /* raster scanline Y by H, then taking the inverse DFT to get the	*/
    /* correlation response as a function of offset, and detecting the	*/
    /* peak of the resulting correlation response.			*/
    int i;
    uint32 v;
    double rr, ii;
    /* Load scanline							*/
    for (i=0; i<N; i++) {
	v = *(raster+Y*w+i);
	hScan[i].re = greyscale(v);
	}
    /* Transform scanline						*/
    fftw_one(FFT,hScan,HScan);
    /* Correlate scanline						*/
    for (i=0; i<N; i++) {
	rr = HScan[i].re * H[i].re - HScan[i].im * H[i].im;
	ii = HScan[i].im * H[i].re + HScan[i].re * H[i].im;
	HCorr[i].re = rr;
	HCorr[i].im = ii;
	}
    /* Inverse transform correlation					*/
    fftw_one(IFFT,HCorr,hCorr);
#if 0
    /* Display results							*/
    plot(hCorr,N,1.0/N/N/100000*PLOT_COLUMNS);
#endif
    }
#endif


void
plot_scanline(uint32 *raster,uint32 w,uint32 Y) {
    /* Plot greyscale values of scanline Y having width w in raster.	*/
    uint32 i;
    uint32 j,n;
    uint32 *p;
    printf("# scanline=%ld\n",(long)Y);
    p = raster + w*Y;
    for (i=0; i<w; i++,p++) {
	n = greyscale(*p)*PLOT_COLUMNS/(256*3)+1;
	printf("%3d ",greyscale(*p));
	for (j=0; j<n; j++) putchar('*');
	putchar('\n');
	}
    }


#if 0
    /* Cut out frames based on manual coordinates			*/
    for (i=0; i<24; i++) {
	char fname[32];
	sprintf(fname,"_%02d.tif",i);
	crop(raster,w,h,xres,yres,resunits,
		START_OF_FIRST_FRAME+i*FILM_FRAME_HEIGHT,CROPPED_FRAME_LEFT,CROPPED_FRAME_HEIGHT,CROPPED_FRAME_WIDTH,fname,0);
	}
#elif 0
    /* Output a text plot of the sprocket centerline scanline		*/
    plot_scanline(raster,w,RIGHT_SPROCKET_CENTERLINE);
#endif


int
xor_correlate(uint32 *raster,uint32 w,uint32 h,int X0,int X,int Y) {
    /* Given an input raster of size (w,h), we examine scanline Y at	*/
    /* offset X0 for the sequence of X pixels to attempt recognition of	*/
    /* the perforation by XOR correlation.  We presume to thus return	*/
    /* the offset of the leading edge of the perforation, such that the	*/
    /* caller should add half the perforation height to yield the	*/
    /* correct framing.							*/
    int weight;
    int offset;
    int max_weight=0, max_weight_offset=0;
    int x;
    uint32 *scanline=raster+Y*w+X0;
    for (offset=0; offset<X; offset++) {
	weight = 0;
	/* Pixels in scanline where we expect white			*/
	for (x=offset; x<offset+FILM_PERFORATION_HEIGHT; x++) {
	    if (threshold(scanline[x%X])) weight++;
	    }
	/* Pixels in scanline where we expect black			*/
	for (x=offset+FILM_PERFORATION_HEIGHT; x<offset+X; x++) {
	    if (!threshold(scanline[x%X])) weight++;
	    }
	if (weight>max_weight) {
	    max_weight = weight;
	    max_weight_offset = offset;
	    }
	}
    return(max_weight_offset+X0);
    }


main(int argc, char **argv) {
    TIFF *tif;
    uint32 w, h;
    size_t npixels;
    uint32 *raster;
    float xres, yres;
    uint16 resunits;
    int i;

    if (argc!=3) {
	fprintf(stderr,"Use: tiffcine input-TIFF-file starting-frame-number\n");
	return(-1);
	}

    /* Open input TIFF file						*/
    tif = TIFFOpen(argv[1],"r");
    if (tif==NULL) {
	fprintf(stderr,"Cannot open input TIFF file\n");
	return(-1);
	}

    /* Load and characterize TIFF image as an RGB raster		*/
    TIFFGetField(tif,TIFFTAG_IMAGEWIDTH,&w);
    TIFFGetField(tif,TIFFTAG_IMAGELENGTH,&h);
    TIFFGetField(tif,TIFFTAG_RESOLUTIONUNIT,&resunits);
    TIFFGetField(tif,TIFFTAG_XRESOLUTION,&xres);
    TIFFGetField(tif,TIFFTAG_YRESOLUTION,&yres);
    npixels = w * h;
    printf("# file=%s w=%ld h=%ld npixels=%ld\n",argv[1],
	    (long)w,(long)h,(long)npixels);
    raster = (uint32*) _TIFFmalloc(npixels * sizeof(uint32));
    if (raster==NULL) { fprintf(stderr,"Out of memory!\n"); return(-1); }
    if (TIFFReadRGBAImage(tif,w,h,raster,0)==0) {
	fprintf(stderr,"Cannot read image from TIFF file!\n");
	return(-1);
	}

#if 0
    /* FFT correlation methods						*/
    /* Initialize FFT plans						*/
    FFT = fftw_create_plan(FILM_FRAME_HEIGHT,FFTW_FORWARD,FFTW_ESTIMATE);
    IFFT = fftw_create_plan(FILM_FRAME_HEIGHT,FFTW_BACKWARD,FFTW_ESTIMATE);

    /* Initialize ideal perforation scanline and its DFT conjugate	*/
    for (i=0; i<FILM_PERFORATION_HEIGHT; i++) hPerf[i].re = 1.0 /* greyscale_white */;
    for (i=FILM_PERFORATION_HEIGHT; i<FILM_FRAME_HEIGHT; i++) hPerf[i].re = 0.0 /* greyscale_black */;
    fftw_one(FFT,hPerf,HPerf);
    for (i=0; i<FILM_FRAME_HEIGHT; i++) hPerf[i].im = -hPerf[i].im;

#if 0
    fft_correlate(raster,w,h,FILM_FRAME_HEIGHT,RIGHT_SPROCKET_CENTERLINE,HPerf);
#else
    {
    /* Make correlation map image of dimensions (FILM_FRAME_HEIGHT,h)	*/
    uint32 *plotimg, v;
    int x, y;
    double max=0.0;
    plotimg = (uint32 *) _TIFFmalloc(FILM_FRAME_HEIGHT*h*sizeof(uint32));
    if (plotimg==NULL) { fprintf(stderr,"Out of memory!\n"); return(-1); }
    for (y=0; y<h; y++) {
	fft_correlate(raster,w,h,FILM_FRAME_HEIGHT,y,HPerf);
	for (x=0; x<FILM_FRAME_HEIGHT; x++) {
	    if (max<hCorr[x].re) max = hCorr[x].re;
	    v = (uint32)(hCorr[x].re/FILM_FRAME_HEIGHT/FILM_FRAME_HEIGHT*4);
	    v &= 0xff;
	    plotimg[y*FILM_FRAME_HEIGHT+x] = (v | (v<<8) | (v<<16));
	    }
	}
    crop(plotimg,FILM_FRAME_HEIGHT,h,xres,yres,resunits,0,0,FILM_FRAME_HEIGHT,h,"_response.tif",0);
    printf("max(hCorr) = %f\n",max);
    _TIFFfree(plotimg); plotimg = NULL;
    }
#endif
#endif

    /* XOR correlation methods						*/
    {
    int frames = w / FILM_FRAME_HEIGHT;
    int offset;
    int frame;
    int filenumber=atoi(argv[2]);
    char fname[32];
    printf("%d frames(s) in strip\n",frames);
    /* Frame 0								*/
    offset = xor_correlate(raster,w,h,0,FILM_FRAME_HEIGHT,RIGHT_SPROCKET_CENTERLINE);
    for (offset=-(FILM_FRAME_HEIGHT/2),frame=0; frame<frames; frame++) {
	offset = xor_correlate(raster,w,h,offset+FILM_FRAME_HEIGHT/2,FILM_FRAME_HEIGHT,RIGHT_SPROCKET_CENTERLINE);
	offset += FILM_PERFORATION_HEIGHT/2;
	sprintf(fname,"_frm%04d.tif",filenumber);
	filenumber++;
	printf("%d %d %s\n",frame,offset,fname);
	crop(raster,w,h,xres,yres,resunits,offset,CROPPED_FRAME_LEFT,CROPPED_FRAME_HEIGHT,CROPPED_FRAME_WIDTH,fname,1);
	}
    }

    /* Discard image and close TIFF file.				*/
    _TIFFfree(raster); raster = NULL;
    TIFFClose(tif); tif = NULL;

#if FFTW
    /* Clean up FFTW plans						*/
    if (FFT) { fftw_destroy_plan(FFT); FFT = NULL; }
    if (IFFT) { fftw_destroy_plan(IFFT); IFFT = NULL; }
#endif

    return(0);
    }