#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "paulslib.h"
#include "bitmaplib.h"

/*
	Proof of concept code to turn lens distorted
   (No really a distortion) into idealis perspective
   projections (aka pinhole camera)
*/

void GiveUsage(char *);

/* 
	Alpha1 and alpha2 are constants for the transformation 
	equations, they will be small positive numbers
	Determined for a particular lens by trial and error
	If the lines are bowed inwards (concave) then increase
   the appropriate constant. If they are bowed out (convex)
   then decrease the corresponding constant.
	eg: for test image seen so far
	Lens    alphax  alphay
	35mm     0.02   0.075
   40mm     0.016  0.059
   50mm     0.006  0.03
*/
double alphax = 0.02, alphay = 0.08;

/* Level of antialiasing */
int aa = 3;

/* Size of the output image */
int nxout,nyout;

int main(int argc,char **argv)
{
	int nx,ny,depth;
	int i,j,i2,j2,ii,jj,aacount;
	double x,y,x2,x3,y2,y3,r;
	COLOUR csum;
	BITMAP4 *image,*output1,*output2;
	BITMAP4 c,grey = {128,128,128,255};
	FILE *fptr;

	/* If there are any arguments assume alphx and alphay */
	if (argc < 2) 
		GiveUsage(argv[0]);

	/* Parse command line */
	for (i=1;i<argc;i++) {
		if (strcmp(argv[i],"-k") == 0) {
			alphax = atof(argv[i+1]);
			alphay = atof(argv[i+2]);
		}
      if (strcmp(argv[i],"-a") == 0) {
         aa = atof(argv[i+1]);
			if (aa < 1)
				aa = 1;
      }
      if (strcmp(argv[i],"-h") == 0) {
         GiveUsage(argv[0]); 
      }
	}

	/* Read TGA file - see bitmaplib.c/h */
	if ((fptr = fopen(argv[argc-1],"rb")) == NULL) {
		fprintf(stderr,"Failed to open image file \"%s\"\n",argv[argc-1]);
		exit(-1);
	}
	TGA_Info(fptr,&nx,&ny,&depth);
	fprintf(stderr,"Image: %d x %d x %d\n",nx,ny,depth);
	if ((image = Create_Bitmap(nx,ny)) == NULL) {
		fprintf(stderr,"Failed to malloc image data\n");
		exit(-1);
	}
	if (TGA_Read(fptr,image,&nx,&ny) != 0) {
		fprintf(stderr,"TGA read failed\n");
		exit(-1);
	}
	fclose(fptr);

	/* Create first output image */
	nxout = nx / 2;
	nyout = ny / 2;
   output1 = Create_Bitmap(nxout,nyout);
   for (i=0;i<nxout;i++) 
      for (j=0;j<nyout;j++) 
         Draw_Pixel(output1,nxout,nyout,i,j,grey);

	/* Forward transform */
   for (i=0;i<nxout;i++) {
      for (j=0;j<nyout;j++) {
			aacount = 0;
			csum.r = 0; csum.g = 0; csum.b = 0;
			for (ii=0;ii<aa;ii++) {
				for (jj=0;jj<aa;jj++) {
					x = (2 * (i+ii/(double)aa) - nxout) / (double)nxout;
					y = (2 * (j+jj/(double)aa) - nyout) / (double)nyout; 
					r = x*x + y*y;
					x3 = x / (1 - alphax * r);
					y3 = y / (1 - alphay * r); 
					x2 = x / (1 - alphax * (x3*x3+y3*y3));
					y2 = y / (1 - alphay * (x3*x3+y3*y3));	
					i2 = (x2 + 1) * nx / 2;
					j2 = (y2 + 1) * ny / 2;
					if (i2 >= 0 && i2 < nx && j2 >= 0 && j2 < ny) {
						c = Get_Pixel(image,nx,ny,i2,j2);
						csum.r += c.r;
						csum.g += c.g;
						csum.b += c.b;
						aacount++;
					}
				}
			}
			if (aacount > 0) {
				c.r = csum.r / aacount;
				c.g = csum.g / aacount;
				c.b = csum.b / aacount;
				c.a = 0;
				Draw_Pixel(output1,nxout,nyout,i,j,c);
			}
      }
   }

	/* Write the first image */
   if ((fptr = fopen("1.tga","wb")) != NULL) {
   	Write_Bitmap(fptr,output1,nxout,nyout,12);
   	fclose(fptr);
	}

	/* Create the secnd image */
	output2 = Create_Bitmap(nxout,nyout);
   for (i=0;i<nxout;i++)
      for (j=0;j<nyout;j++)
         Draw_Pixel(output2,nxout,nyout,i,j,grey);

	/* Transform back */ 
   for (i=0;i<nxout;i++) {
      for (j=0;j<nyout;j++) {
         aacount = 0;
         csum.r = 0; csum.g = 0; csum.b = 0;
         for (ii=0;ii<aa;ii++) {
            for (jj=0;jj<aa;jj++) {
               x = (2 * (i+ii/(double)aa) - nxout) / (double)nxout;
               y = (2 * (j+jj/(double)aa) - nyout) / (double)nyout;
               r = x*x + y*y;
               x2 = x * (1 - alphax * r);
               y2 = y * (1 - alphay * r);
               i2 = (x2 + 1) * nxout / 2;
               j2 = (y2 + 1) * nyout / 2;
               if (i2 >= 0 && i2 < nxout && j2 >= 0 && j2 < nyout) {
                  c = Get_Pixel(output1,nxout,nyout,i2,j2);
                  csum.r += c.r;
                  csum.g += c.g;
                  csum.b += c.b;
                  aacount++;
               }
            }
         }
         if (aacount > 0) {
            c.r = csum.r / aacount;
            c.g = csum.g / aacount;
            c.b = csum.b / aacount;
            c.a = 0;
            Draw_Pixel(output2,nxout,nyout,i,j,c);
         }
      }
   }

   if ((fptr = fopen("2.tga","wb")) != NULL) {
   	Write_Bitmap(fptr,output2,nxout,nyout,12);
   	fclose(fptr);
	}

	exit(0);
}

void GiveUsage(char *s)
{
	fprintf(stderr,"Usage: %s [-k n1 n2] [-a n] [-h] tgafilename\n",s);
	exit(-1);
}