/*
 * prende una foto e la riflette su una sfera
 * programma "esher" versione 4
 *
 *  Sergio Steffe' - Laboratorio Sperimentale di Matematica Computazionale
 *  Dipartimento di Matematica - Universita' di Pisa - AA 2004-2005
 *  ultima modifica: 2005-06-03
 *
 * le immagini sono in formato file.pnm ASCII
 *
 * l'occhio e' posto nell' origine
 *
 * la sfera ha centro (disf,0,0) e raggio rasf, 0 < rasf < disf
 *    (flags -d e -r, defaults disf=2 rasf=1). 
 *
 * l'immagine da calcolare e' quadrata lc x lc pixel (flag -l default 512), 
 *     perpendicolare all' asse x e centrata in modo che il cerchio intercettato sia 
 *     massimo.
 *           
 * l'immagine data e' posta sul piano posto dietro l'occhio (flag -b  default 10) 
 *    sul piano  y-z, centrata, di dimensione calcolata sapendo la grandezza di ogni 
 *    pixel (flag -s default 0.1 per pixel), ed e' spostabile verticalmente (flag -v)
 *
 * calcolo dell'immagine:
 * per ogni py, pz ( tra -lx*h/2 e +lx*h/2 e -lz*h/2, lz*h/2) si prende il 
 * raggio da occhio fino a sfera; se fuori sfera si usa un fondo
 * si calcola il punto di intersezione sulla superfice della sfera
 * (xs,yx,zs) e la normale (nx,ny,nz) esterna
 * si calcola il raggio uscente (xs,ys,zs) + t * (ux,uy,uz) , t>0
 * se non ha intersezione col piano della foto data, si mette il valore gray ,
 * altrimenti il valore del pixel dell immagine data che e' piu' vicino (opzionale 
 *     flag -p se nel piano ma fuori dall' immagine)
 *
 */
#define _GNU_SOURCE  /* per getline in stdio.h */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>  /* per exit e return */
#include <unistd.h>  /* per getopt  versione normale */
#include <ctype.h>   /* per isprint */


/* prototipi C99 */
double round(double);
/* fine prototipi */

/* usage */
static void usage(char *command)
{
        fprintf(stderr,
"Usage: %s OUTPUTIMAGE.pnm [OPTIONS] \n"
"-h,            help\n"
"-r radius	radius of the sphere (default 1)\n"
"-d distance	distance of the sphere (default 2)\n"
"-l pixel	pixel of the (square) output image (default 512)\n"
"-b back	distance of the backward input image (default 10)\n"
"-s scale	pixel with for input image (default 0.1)\n"
"-v vertical    vertical position of the center of input image (default 0)\n"
"-i inputimage.pnm  must be pnm ASCII color image (header P3)\n"
"-f frontimage.pnm  optional - ASCII color image (header P3)\n"
"-p             prolong image in the plain (default no)"
"     mirrors the input image in a sphere\n"
         ,command);
}
/* end usage */

int main(argc,argv)
int argc;
char *argv[];

{
/* per la lettura del comando */
char *nomecomando;    /* nome del comando */
int c;                /* per il getopt */
char *cvalue = NULL ; /* pre le stringhe delle opzioni */
double rasf,disf;     /* raggio e distanza della sfera */
int lc;	              /* pixel immagine prodotta   */
double ba;            /* distanza backward dell'immagine da specchiare*/
double sca;	      /* scala immagine da specchiare*/
double ver;	      /* spostamento verticale del centro immagine da specchiare*/
char nomein[256], nomefondo[256], nomeout[256]; /* nomi dei files */
int fondosi, inputsi, prolongsi;/* presenza fondo  input  */
/* inizio lettura files dei dati */
FILE *imgin, *imgout, *imgfondo, *fopen(); /* immagine input output */
char imgtype[256], stringaletta[256];
char * line = NULL; /* linea  in immagine input */
size_t len = 0;  /* linea  in immagine input */
ssize_t read;  /* linea  in immagine input */
int zonadati; /* 0 ingresso - 1 lette dimensioni - 2 letto maxvalue  */
int ir,ig,ib; /* contatori colori */
int  xi,yi,toti,vmaxi ; /* immagine ingresso xi x yi profondita vmax 3 colori */
int  *mater, *mateg, *mateb; /* immagine ingresso valori r g b */ 
int xf,yf,zf,totf,vmaxf ; /* immagine fondo xi x yi profondita vmax 3 colori */
int *fondr,*fondg,*fondb;  /* immagine fondo valori r g b */
/* dati per il file di uscita */
int vmax, tot,  *sferar, *sferag, *sferab;
int skyr,skyg,skyb, grayr,grayg,grayb, metalr,metalg,metalb;
/* variabili per i calcoli */
int i,j,ii,jj;        /* general purpose interi */
double aa,bb,cc,nnn,phi,phizero;  /* general purpose double */
double nx,ny,nz,sx,sy,sz,ux,uy,uz; /* normale, punto della sfera e versore riflesso */
double t,t1,t2,tcr,vcr; /* calcolo intersezione raggi */
double px,py,pz; /* coordinate pixel partenza px=disf-rasf */
double xx,yy,zz; /* punto di intersezione */
double hc,hd,hf; /* ampiezza pixels */
double pigreco;
/*fine dichiarazioni */

/* variabili de default */
pigreco=atan(1.0)*4.0;
rasf=1.0; /* raggio sfera */
disf=2.0; /* distanza sfera */
lc=512;   /* pixel immagine output */
ba=-10;   /* distanza back immagine */
sca=0.1;  /* scala immagine */
ver=0.0;  /* shift verticale */
fondosi=0;
inputsi=0;
prolongsi=0; /* prolunga immagine data nel piano in cui giace */

/* getopt */
nomecomando=argv[0];
while ((c = getopt (argc, argv, "hr:d:l:b:s:v:i:f:p")) != -1)
switch (c)
           {
            case 'h':
                usage(nomecomando);
                exit(EXIT_FAILURE);
                break;
            case 'r':
                cvalue = optarg;
                if( 1 != sscanf(cvalue,"%lg\n",&rasf) ) {printf("parametro -r %s non valido\n",cvalue); exit(EXIT_FAILURE);}
                break;
            case 'd':
                cvalue = optarg;
                if( 1 != sscanf(cvalue,"%lg\n",&disf) ) {printf("parametro -d %s non valido\n",cvalue); exit(EXIT_FAILURE);}
                break;
            case 'l':
                cvalue = optarg;
                if( 1 != sscanf(cvalue,"%d\n",&lc)) {printf("parametro -l %s non valido\n",cvalue); exit(EXIT_FAILURE);}
                break;
            case 'b':
                cvalue = optarg;
                if( 1 != sscanf(cvalue,"%lg",&ba)) {printf("parametro -b %s non valido\n",cvalue); exit(EXIT_FAILURE);}
                break;
           case 's':
                cvalue = optarg;
                if( 1 != sscanf(cvalue,"%lg",&sca)) {printf("parametro -s %s non valido\n",cvalue); exit(EXIT_FAILURE);}
                break;
           case 'v':
                cvalue = optarg;
                if( 1 != sscanf(cvalue,"%lg",&ver)) {printf("parametro -v %s non valido\n",cvalue); exit(EXIT_FAILURE);}
                break;
	   case 'i':
                cvalue = optarg;
                if( 1 != sscanf(cvalue,"%s",nomein)) {printf("parametro -i %s non valido\n",cvalue); exit(EXIT_FAILURE);}
	        inputsi=1;
                break;
           case 'f':
                cvalue = optarg;
                if( 1 != sscanf(cvalue,"%s",nomefondo)) {printf("parametro -f %s non valido\n",cvalue); exit(EXIT_FAILURE);}
                fondosi=1;
                break;
	   case 'p':
		prolongsi=1;
		break;
           case '?':
                if (isprint (optopt))
               fprintf (stderr, "Unknown option `-%c'.\n", optopt);
             else
                fprintf (stderr, "Unknown option character `\\x%x'.\n", optopt);
                exit(EXIT_FAILURE);
           default:
             abort ();

          }
if (optind < argc) { 
/* dovrebbe esserci solo il nome del file di output */
 sscanf(argv[optind],"%s",nomeout);
	}
else {
/* manca nome del file di output */
	printf("Manca il nome del file di output !\n"); 
	exit(EXIT_FAILURE);
	}
/* ottenute le opzioni si controlla che ci sia tutto e vada bene */

/* raggio positivo ? */
if (rasf < 0) {
		printf("il raggio deve essere positivo !\n");
		exit(EXIT_FAILURE);
	      }
/* distanza sfera da occhio positiva e superiore al raggio */
if (disf < 0) {
	  	printf("la distanza della sfera deve essere positiva !\n");
                exit(EXIT_FAILURE);
              }
if (disf < rasf) {
                printf("la distanza della sfera deve essere maggiore del raggio !\n");
                exit(EXIT_FAILURE);
              }
/* numero pixel tra 100 e 10000 */
if (lc < 100)  {
                printf("suvvia, almeno 100 pixel per lato!\n");
                exit(EXIT_FAILURE);
	      }
if (lc > 10000)  {
                printf("al massimo 10000 pixel\n");
                exit(EXIT_FAILURE);
              }
/* ba numero negativo */
if (ba > 0) ba=-ba; 
/* scala numero positivo */
if (sca < 0) sca=-sca;
/* input file ? */
if (inputsi == 0)  {
                printf("manca il file di input !\n");
                exit(EXIT_FAILURE);
              }
/*  quasi a posto ... */


/* apri file input in lettura  */
if ((imgin = fopen(nomein,"r")) == NULL) {
        printf("\n non posso leggere %s \n",nomein);
        exit(EXIT_FAILURE);
        };

/* apri file output in scrittura  */
if ((imgout = fopen(nomeout,"w")) == NULL) {
        printf("\n non posso scrivere %s \n",nomeout);
        exit(EXIT_FAILURE);
        };
/* se presente, apri il file del fondo */
if (fondosi == 1) if ((imgfondo = fopen(nomefondo,"r")) == NULL) {
        printf("\n non posso leggere %s \n",nomefondo);
        exit(EXIT_FAILURE);
        };


/* lettura file in ingresso e controllo sommario tipo file */
zonadati = 0;
ir=0;
ig=0;
ib=0;
/* controlla tipo immagine (P3) */
if ((read = getline(&line, &len, imgin)) != -1) {
	sscanf(line,"%s",imgtype);
        if ( strcmp("P3",imgtype) != 0 ) {
        printf("\nTipo di immagine errato: %s\n",nomein);
        exit(EXIT_FAILURE);
	}
 }
/* leggi il contenuto */
while ((read = getline(&line, &len, imgin)) != -1) {
/* scarta eventuale riga di commento */
	sscanf(line,"%s",stringaletta);
	if (strcmp("#",stringaletta) == 0) {	
/*		printf("scarto commento di  %zu  caratteri:\n", read);
                 printf("%s", line); */
			}
	else  if (zonadati == 0){
/* leggi dimensioni immagine */
	sscanf(line,"%i %i\n",&xi,&yi);
	toti=xi*yi;
/* alloca spazio per immagin ingresso mate e uscita sfera */
	mater =  calloc(toti, sizeof(toti));
	mateg =  calloc(toti, sizeof(toti));
	mateb =  calloc(toti, sizeof(toti));
	if ((mater == NULL)||(mateg == NULL)||(mateb == NULL)) {
		printf(" manca memoria per lettura immagine (3x %d bytes)\n",toti);
		exit(EXIT_FAILURE);	
		}
	zonadati = 1;
	}
	else  if (zonadati == 1){
        sscanf(line,"%i \n",&vmaxi);
	zonadati = 2;
		}
/* leggi  dati mate */
	else if (zonadati == 2) {
		sscanf(line,"%i\n", &(mater[ir]));
		ir++;
		zonadati=3;	
		}
	else if (zonadati == 3) {
		sscanf(line,"%i\n", &(mateg[ig]));
                ig++;
                zonadati=4;
		}
	else if (zonadati == 4) {	
		sscanf(line,"%i\n", &(mateb[ib]));
                ib++;
                zonadati=2;
		}
 }
/* controllo consistenza dati */
if ((ir != ig) || (ir != ib)) {
        printf("Immagine input con inconsistenti: red=%d green=%d blu=%d\n ",ir,ig,ib);
	exit(EXIT_FAILURE);
	}
fclose(imgin);
phizero=yi;
phizero=atan(phizero/xi);
/* a questo punto mater, mateg, mateb matrici xi x yi = toti con immagine input */

/* se presente, leggi immagine di fondo  fai controllo sommario tipo file */
if (fondosi == 1) {
  zonadati = 0;
  ir=0;
  ig=0;
  ib=0;
/* controlla tipo immagine (P3) */
  if ((read = getline(&line, &len, imgfondo)) != -1) {
        sscanf(line,"%s",imgtype);
        if ( strcmp("P3",imgtype) != 0 ) {
        printf("\nTipo di immagine errato: %s\n",nomefondo);
        exit(EXIT_FAILURE);
        }
  }
/* leggi il contenuto */
  while ((read = getline(&line, &len, imgfondo)) != -1) {
/* scarta eventuale riga di commento */
        sscanf(line,"%s",stringaletta);
        if (strcmp("#",stringaletta) == 0) {
/*                printf("scarto commento di  %zu  caratteri:\n", read);
                 printf("%s", line); */
                        }
        else  if (zonadati == 0){
/* leggi dimensioni immagine */
        sscanf(line,"%i %i\n",&xf,&yf);
        totf=xf*yf;
/* alloca spazio per immagine del fondo */
        fondr =  calloc(totf, sizeof(totf));
        fondg =  calloc(totf, sizeof(totf));
        fondb =  calloc(totf, sizeof(totf));
        if ((fondr == NULL)||(fondg == NULL)||(fondb == NULL)) {
                printf(" manca memoria per lettura immagine (3x %d bytes)\n",totf);
                exit(EXIT_FAILURE);
                }
        zonadati = 1;
        }
        else  if (zonadati == 1){
        sscanf(line,"%i \n",&vmaxf);
        zonadati = 2;
                }
/* leggi  dati fondo */
        else if (zonadati == 2) {
                sscanf(line,"%i\n", &(fondr[ir]));
                ir++;
                zonadati=3;
                }
        else if (zonadati == 3) {
                sscanf(line,"%i\n", &(fondg[ig]));
                ig++;
                zonadati=4;
                }
        else if (zonadati == 4) {
                sscanf(line,"%i\n", &(fondb[ib]));
                ib++;
                zonadati=2;
                }
 }
/* controllo consistenza dati */
  if ((ir != ig) || (ir != ib)) {
        printf("Immagine fondo con dati inconsistenti: red=%d green=%d blu=%d\n ",ir,ig,ib);
        exit(EXIT_FAILURE);
        }
fclose(imgfondo);
/* a questo punto fondr fondg fondb  matrici xf x yf = toti con immagine input */
}

/* default per calcolo  sfera */

/* colori del fondo (blu, grigio, metal) da usare se non presente immagine di fondo */
skyr=12;
skyg=12;
skyb=228;
grayr=128;
grayg=128;
grayb=128;
metalr=180;
metalg=180;
metalb=180;
vmax=256;  /* 8 bit per colore */
tot=lc*lc; /* pixel immagine sfera */
/* alloca spazio per immagin uscita sfera */
sferar =  calloc(tot, sizeof(vmax));
sferag =  calloc(tot, sizeof(vmax));
sferab =  calloc(tot, sizeof(vmax));
/* l'immagine da calcolare sta su un piano parallelo a y-z e distante da occhio */
/* disf - rasf  e quindi prendiamo come ampiezza dei pixel hc */
hc = 2.0*rasf*disf/(sqrt(disf*disf-rasf*rasf)*lc);
/* in modo che il cerchio immagine occupi quasi tutto il quadrato */

/* default immagine da specchiare */

hd=sca; /* grandezza pixel: immagine occupa (xi*hd) x (yi*hd) spazio */
xx=ba;

/* default del fondo */
/* viene usata solo la parte centrale del fondo */
if (fondosi == 1) {
   if (xf > yf) zf = yf; else zf =xf;
/* si scala la parte centrale di fondo sulla immagine da costruire */
/* gli zf pixel per hf grandezza pixel = lc * hc  */
hf= (lc*hc)/zf;

}

/* inizio calcolo immagine riflessa */

/* prendo tutti i possibili pixel della immagine da costruire */
   px=disf;
for (i=0;i<lc;i++) for (j=0;j<lc;j++) {
   py = - hc* lc/2.0 + i* hc;
   pz =  hc *lc/2.0 -j * hc;
/* il raggio dall'origine a (px,py,pz) incontra la sfera ? */
/* al variare di t, si tratta di vedere se una parabola ha valori negativi */
/* tcr e' il valore critico della parabola vcr il valore dell parabola */
aa=(px*px+py*py+pz*pz);
bb=px*disf;
tcr= bb/aa;
cc=disf*disf-rasf*rasf;
vcr= tcr*tcr*aa - 2*tcr*bb +cc;
if(vcr > 0) {
       if (fondosi == 0) {
/* sky color */
		sferar[j*lc+i]= skyr;
		sferag[j*lc+i]= skyg;
		sferab[j*lc+i]=	skyb;
		}
	else    {
/* fondo */
		ii=yf - round(yf/2.0 + pz/hf);
		jj=xf - round(xf/2.0 - py/hf);	
		sferar[j*lc+i]=fondr[ii*xf+jj] ;
                sferag[j*lc+i]=fondg[ii*xf+jj] ;
                sferab[j*lc+i]=fondb[ii*xf+jj] ;
		}
	}
else {
/* intanto metti bianco poi si vedra' */
	sferar[j*lc+i]=metalr;
        sferag[j*lc+i]=metalg;
        sferab[j*lc+i]=metalb;
/*                                    */
/* calcolo prima intersezione linea - sfera con t tra 0 e tcr */
t1=(bb-sqrt(bb*bb-aa*cc))/aa;
t2=(bb+sqrt(bb*bb-aa*cc))/aa;
if ( t1 > 0) t=t1; else t=t2;
/* calcolo punto sulla sfera */
sx=t*px;
sy=t*py;
sz=t*pz;
/* calcolo normale nel punto */
nx=sx-disf;
ny=sy;
nz=sz;
nnn=sqrt(nx*nx+ny*ny+nz*nz);
nx=nx/nnn;
ny=ny/nnn;
nz=nz/nnn;
/* calcolo vettore riflesso u = 2*(s*n) n - s */
nnn = sx*nx+sy*ny+sz*nz;
ux = 2*nnn*nx - sx;
uy = 2*nnn*ny - sy;
uz = 2*nnn*nz - sz;
/* ora per t>0 il vettore s + u * t incontra immagine dipartimento ? */
/* intanto incontrera il piano x=xx ? occorre ux non zero */
if (ux == 0) {
/* prova con il bianco inutile  */
        sferar[j*lc+i]=vmax;
        sferag[j*lc+i]=vmax;
        sferab[j*lc+i]=vmax;
	}
	else {
/* calcolo t  */
	t= (sx - xx )/ux;	
	if (t > 0) {
/* prova con il bianco  
        sferar[j*lc+i]=vmax;
        sferag[j*lc+i]=vmax;
        sferab[j*lc+i]=vmax;
	*/
/* calcolo coordinate riflesse su xx(=ba) */
	yy= sy+uy*t;
	zz= sz+uz*t + ver; 
/* calcolo pixel vicini immagine dipartimento */
	ii=round(yy/hd + xi/2.0);
	jj=round(zz/hd + yi/2.0);
	if((ii > 0)&&(ii < xi)&&(jj > 0)&&(jj < yi) ){
	sferar[j*lc+i]=mater[jj*xi+ii];
        sferag[j*lc+i]=mateg[jj*xi+ii];
        sferab[j*lc+i]=mateb[jj*xi+ii];
		}
	/* pixel fuori immagine, si prolunga il colore */
	else if(prolongsi == 1) {
                phi= atan2(zz,yy); /* tra -PI a PI incluso */
		if (phi > pigreco-phizero) {
				ii=0;
				jj=round(yi/2.0-tan(phi)*yi/2.0);
				if(jj > yi-1) jj=yi-1;
				if(jj < 0) jj=0;
				sferar[j*lc+i]=mater[jj*xi+ii];
       				sferag[j*lc+i]=mateg[jj*xi+ii];
        			sferab[j*lc+i]=mateb[jj*xi+ii];
				}
                else if (phi > phizero) {
				ii=round(xi/2.0+xi/(tan(phi)*2.0));
				jj=yi-1;
				if (ii < 0) ii=0;
				if (ii > xi-1) ii=xi-1;
				sferar[j*lc+i]=mater[jj*xi+ii];
        			sferag[j*lc+i]=mateg[jj*xi+ii];
        			sferab[j*lc+i]=mateb[jj*xi+ii];
					}
		else if (phi > -1.0*phizero) {
				ii=xi-1;
				jj=round(yi/2.0+tan(phi)*yi/2.0);
				if(jj > yi-1) jj=yi-1;
                                if(jj < 0) jj=0;
				sferar[j*lc+i]=mater[jj*xi+ii];
        			sferag[j*lc+i]=mateg[jj*xi+ii];
        			sferab[j*lc+i]=mateb[jj*xi+ii];
				}
	        else if (phi > -1.0*pigreco+phizero) {
				ii=round(xi/2.0-xi/(tan(phi)*2.0));
				jj=0;
				if (ii < 0) ii=0;
                                if (ii > xi-1) ii=xi-1;
				sferar[j*lc+i]=mater[jj*xi+ii];
        			sferag[j*lc+i]=mateg[jj*xi+ii];
        			sferab[j*lc+i]=mateb[jj*xi+ii];
				}
		else {
				ii=0;
				jj=round(yi/2.0 -tan(phi)*yi/2.0);
				if(jj > yi-1) jj=yi-1;
                                if(jj < 0) jj=0;
				sferar[j*lc+i]=mater[jj*xi+ii];
        			sferag[j*lc+i]=mateg[jj*xi+ii];
        			sferab[j*lc+i]=mateb[jj*xi+ii];
				}
		}
	    else {
				sferar[j*lc+i]=grayr;
                                sferag[j*lc+i]=grayg;
                                sferab[j*lc+i]=grayb;


		}
	   }
	}

     }
  }

/* scrivi tipo immagine, commento, dimensioni immagine e sfera */

fprintf(imgout,"%s\n",imgtype);
fprintf(imgout,"%i %i\n",lc,lc);
fprintf(imgout,"%i\n",vmax);
for (i=0;i<tot;i++) fprintf(imgout,"%d\n%d\n%d\n",sferar[i],sferag[i],sferab[i]);
fclose(imgout);
	
return EXIT_SUCCESS;
}