// ................................................................. 3dhisto.cc
//   30.JUL.2007

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <corso.h>

extern XVideoData xvd;

main()
{
	int i,icol,iPrism,ii,iWin,ix,ix0,iy0,iy,j,l,k,kk,nComm,nPar,nX,nY,nTot,xc,yc;
	int iind[5];
	int *ind,*ParFont;
	int **face;
	unsigned long fcol[5],color[5];
	double alpha,alphadeg,beta,betadeg,DeltaX,DeltaY,dist,dx,dy,gamma,gammadeg,
	x,x0,xLen,y,y0,yLen,zEnlarge,zHigh,zLen,zScale;
	double fintens[5],light[3],rnew[3],rold[3],yyy[5];
	const double deg2rad=M_PI/180.0;
	double *yy;
	double **CosDir,**DataMat,**NewXax,**NewYax,**NewZax,**OldXax,**OldYax,
	**OldZax,**RotCosDir,**rotmat;
	double ***NewPrism,***OldPrism;
	char **comandi,**ParName,**ParVal;
	clock_t musec1,musec2;
	timespec ts1,ts2;
	XColor col;
	XPoint ptswap;
	XPoint pt2d[5];
	
	// ...........................................................................
	iWin=0;
	// ............................................................. modo grafico
	StartXWindow(800,600,0,0,"3D Plot",2);
	XSelectInput(xvd.display,xvd.win[0],ButtonPressMask|ButtonReleaseMask|
			PointerMotionMask|KeyPressMask|KeyReleaseMask|ExposureMask|
					StructureNotifyMask);
	xc=xvd.width[0]/2;
	yc=xvd.height[0]/2;
	// .......................................................... numero di punti
	nX=7;
	nY=7;
	nTot=nX*nY;
	xLen=400.0;
	yLen=400.0;
	zLen=400.0;
	dist=200000.0;
	DeltaX=xLen/double(nX);
	DeltaY=yLen/double(nY);
	dx=0.5*DeltaX;
	dy=0.5*DeltaY;
	zEnlarge=1.0;
	// ...................................................... allocazione memoria
	ind=new int [nTot];
	face=new int * [5];
	for (i=0;i<5;i++) face[i]=new int [4];
	yy=new double [nTot];
	CosDir=new double * [5];
	for (i=0;i<5;i++) CosDir[i]=new double [3];
	RotCosDir=new double * [5];
	for (i=0;i<5;i++) RotCosDir[i]=new double [3];
	DataMat=new double * [7];
	for (i=0;i<7;i++) DataMat[i]=new double [7];
	rotmat=new double * [3];
	for (i=0;i<3;i++) rotmat[i]=new double[3];
	NewXax=new double * [nX];
	OldXax=new double * [nX];
	NewYax=new double * [nY];
	OldYax=new double * [nY];
	for (i=0;i<nX;i++)
	{
		NewXax[i]=new double [3];
		OldXax[i]=new double [3];
	}
	for (i=0;i<nY;i++)
	{
		NewYax[i]=new double [3];
		OldYax[i]=new double [3];
	}
	
	NewPrism=new double ** [nTot];
	OldPrism=new double ** [nTot];
	for (i=0;i<nTot;i++)
	{
		NewPrism[i]=new double * [9];
		OldPrism[i]=new double * [9];
		for (j=0;j<9;j++)
		{
			NewPrism[i][j]=new double [3];
			OldPrism[i][j]=new double [3];
		}
	}
	// .......................................................... facce e spigoli
	for (i=0;i<4;i++)
	{
		face[i][0]=i;
		face[i][1]=(i+1)%4;
		face[i][2]=face[i][1]+4;
		face[i][3]=face[i][0]+4;
	}
	for (i=0;i<4;i++) face[4][i]=4+i;
	CosDir[0][0]=CosDir[0][2]=0.0;
	CosDir[0][1]=-1.0;
	CosDir[1][1]=CosDir[1][2]=0.0;
	CosDir[1][0]=1.0;
	CosDir[2][0]=CosDir[2][2]=0.0;
	CosDir[2][1]=1.0;
	CosDir[3][1]=CosDir[3][2]=0.0;
	CosDir[3][0]=-1.0;
	CosDir[4][0]=CosDir[4][1]=0.0;
	CosDir[4][2]=1.0;
	// ........................................................... direzione luce
	light[0]=light[1]=-sqrt(3.0);
	light[2]=sqrt(3.0);
	// ............................................................. matrice dati
	for (i=0;i<7;i++) for (j=0;j<7;j++)
	{
		x=(double(i)-3.0)*(double(i)-3.0)+(double(j)-3.0)*(double(j)-3.0);
		DataMat[i][j]=exp(-x/4.0);
	}
	// ......................................................... massimo e minimo
	for (i=0;i<nX;i++) for (j=0;j<nY;j++)
	{
		if (i==0 && j==0) zHigh=DataMat[0][0];
		else if (zHigh<DataMat[i][j]) zHigh=DataMat[i][j];
	}
	zScale=double(zLen)/zHigh;
	// ............................................................ tavola colori
	for (i=0;i<5;i++)
	{
		col.red=int(65535.0*double(8-i)/8.0);
		col.green=int(65535.0*double(8-i)/8.0);
		col.blue=int(65535.0*double(8-i)/8.0);
		XAllocColor(xvd.display,xvd.colormap,&col);
		color[i]=col.pixel;
	}
	// ............................................. allocazione memoria per menu'
	nPar=7;
	ParName=new char * [nPar];
	ParVal=new char * [nPar];
	for (i=0;i<nPar;i++) ParVal[i]=new char [24];
	ParName[0]="Origine x";
	ParName[1]="Origine y";
	ParName[2]="a";
	ParName[3]="b";
	ParName[4]="g";
	ParName[5]="Distanza";
	ParName[6]="Ingrandimento z";
	nComm=1;
	comandi=new char * [nComm];
	comandi[0]="Fine";
	ParFont=new int [nPar];
	for (i=0;i<nPar;i++) ParFont[i]=1;
	ParFont[2]=ParFont[3]=ParFont[4]=5;
	// ................................................. origine delle coordinate
	ix0=100;
	iy0=500;
	// ................................................................... angoli
	alphadeg=0.0;
	betadeg=0.0;
	gammadeg=0.0;
	// ....................................................... parametri per menu
	sprintf(ParVal[0],"%.d",ix0);
	sprintf(ParVal[1],"%.d",xvd.height[0]-iy0);
	sprintf(ParVal[2],"%.4lf",alphadeg);
	sprintf(ParVal[3],"%.4lf",betadeg);
	sprintf(ParVal[4],"%.4lf",gammadeg);
	sprintf(ParVal[5],"%.0lf",dist);
	sprintf(ParVal[6],"%.4lf",zEnlarge);
	// .......................................................... disegna il menu
	XWinPermMenu(1,"3D Histo",ParFont,ParName,ParVal,nPar,16,comandi,nComm,0);
	// .............................................................. sincronizza
	XSync(xvd.display,TRUE);
	// .................................................................... ciclo
	ts1.tv_sec=0;
	for (;;)
	{
		XSync(xvd.display,FALSE);
		musec1=clock();
		// .......................................................... leggi il menu
		ii=XWinPermMenu(1,"3D Histo",ParFont,ParName,ParVal,nPar,16,comandi,nComm,1);
		if (ii==0) break;
		//if (ii==nComm)
		gammadeg+=0.1;
		if (TRUE)
		{ // ......................................................... leggi i dati
			ix0=atoi(ParVal[0]);
			iy0=xvd.height[0]-atoi(ParVal[1]);
			alphadeg=atof(ParVal[2]);
			betadeg=atof(ParVal[3]);
			gammadeg=atof(ParVal[4]);
			dist=atof(ParVal[5]);
			zEnlarge=atof(ParVal[6]);
			// .................................................... gradi -> radianti
			alpha=alphadeg*deg2rad;
			beta=betadeg*deg2rad;
			gamma=gammadeg*deg2rad;
			// ................................................. matrice di rotazione
			rotmat[0][0]=cos(gamma)*cos(beta);
			rotmat[0][1]=-sin(gamma)*cos(beta);
			rotmat[0][2]=sin(beta);
			rotmat[1][0]=cos(gamma)*sin(beta)*sin(alpha)+sin(gamma)*cos(alpha);
			rotmat[1][1]=cos(gamma)*cos(alpha)-sin(gamma)*sin(beta)*sin(alpha);
			rotmat[1][2]=-cos(beta)*sin(alpha);
			rotmat[2][0]=sin(gamma)*sin(alpha)-cos(gamma)*sin(beta)*cos(alpha);
			rotmat[2][1]=sin(gamma)*sin(beta)*cos(alpha)+sin(alpha)*cos(gamma);
			rotmat[2][2]=cos(beta)*cos(alpha);
			// .............................. coordinate tridimensionali prerotazione
			iPrism=0;
			for (i=0;i<nX;i++)
			{
				x0=double(i)*DeltaX+dx;
				for (j=0;j<nY;j++)
				{
					y0=double(j)*DeltaY+dy;
					OldPrism[iPrism][8][0]=x0;
					OldPrism[iPrism][8][1]=y0;
					OldPrism[iPrism][8][2]=0.0;
					for (k=0;k<4;k++)
					{
						if (k==0 || k==3) OldPrism[iPrism][k][0]=x0-dx;
						else OldPrism[iPrism][k][0]=x0+dx;
						OldPrism[iPrism][k+4][0]=OldPrism[iPrism][k][0];
						if (k<2) OldPrism[iPrism][k][1]=y0-dy;
						else OldPrism[iPrism][k][1]=y0+dy;
						OldPrism[iPrism][k+4][1]=OldPrism[iPrism][k][1];
						OldPrism[iPrism][k][2]=0.0;
						OldPrism[iPrism][k+4][2]=DataMat[i][j]*zScale*zEnlarge;
					}
					iPrism++;
				}
			}
			// .................................................... assi prerotazione
			for (i=0;i<nX;i++)
			{
				OldXax[i][0]=double(i)*DeltaX+dx;
				OldXax[i][1]=OldXax[i][2]=0.0;
			}
			for (i=0;i<nY;i++)
			{
				OldYax[i][0]=xLen;
				OldYax[i][1]=double(i)*DeltaY+dy;
				OldYax[i][2]=0.0;
			}
			// ................................................. rotazione istogramma
			for (i=0;i<nTot;i++)
			{
				for (j=0;j<9;j++) MatVectMult(rotmat,OldPrism[i][j],NewPrism[i][j],3);
				yy[i]=NewPrism[i][8][1];
			}
			// ....................................................... rotazione assi
			for (i=0;i<nX;i++)  MatVectMult(rotmat,OldXax[i],NewXax[i],3);
			for (i=0;i<nY;i++)  MatVectMult(rotmat,OldYax[i],NewYax[i],3);
			// ........................................... rotazione coseni direttori
			for (i=0;i<5;i++) MatVectMult(rotmat,CosDir[i],RotCosDir[i],3);
			// ............................................ illuminazione delle facce
			for (i=0;i<5;i++)
			{
				fintens[i]=0;
				for (j=0;j<3;j++) fintens[i]+=RotCosDir[i][j]*light[j];
				iind[i]=i;
			}
			DblPtDnSrt(fintens,iind,5);
			for (i=0;i<5;i++) fcol[iind[i]]=color[i];
			// ............................................... ordina i punti ruotati
			for (i=0;i<nTot;i++) ind[i]=i;
			DblPtDnSrt(yy,ind,nTot);
			// .................................................. pulisci la finestra
			XSetForeground(xvd.display,xvd.gc,xvd.white);
			XFillRectangle(xvd.display,xvd.backpix[0],xvd.gc,0,0,800,600);
			// .............................................................. disegno
			for (ii=0;ii<nTot;ii++)
			{
				iPrism=ind[ii];
				// ................................. profondita' dei centri delle facce
				for (j=0;j<5;j++)
				{
					yyy[j]=0.0;
					for (k=0;k<4;k++) yyy[j]+=NewPrism[iPrism][face[j][k]][1];
					yyy[j]/=4.0;
				}
				// .................................................... ordina le facce
				for (j=0;j<5;j++) iind[j]=j;
				DblPtDnSrt(yyy,iind,5);
				// ................................................... disegna le facce
				for (i=0;i<5;i++)
				{
					j=iind[i];
					// ................................... proietta i punti sullo schermo
					for (kk=0;kk<4;kk++)
					{
						pt2d[kk].x=xc+int((NewPrism[iPrism][face[j][kk]][0]+double(ix0-xc))
								*dist/(dist+NewPrism[iPrism][face[j][kk]][1])+0.5);
						pt2d[kk].y=yc-int((NewPrism[iPrism][face[j][kk]][2]+double(yc-iy0))
								*dist/(dist+NewPrism[iPrism][face[j][kk]][1])+0.5);
					}
					pt2d[4]=pt2d[0];
					// .............................................. pulisci il poligono
					/*
					if (j==0 || j==4) icol=color[0];
					else icol=color[50];
					icol=color[90-i*15];
					*/
					XSetForeground(xvd.display,xvd.gc,fcol[j]);
					XFillPolygon(xvd.display,xvd.backpix[0],xvd.gc,pt2d,4,Convex,
											 CoordModeOrigin);
					// .............................................. disegna il reticolo
					XSetForeground(xvd.display,xvd.gc,xvd.black);
					XDrawLines(xvd.display,xvd.backpix[0],xvd.gc,pt2d,5,
										 CoordModeOrigin);
				}
			}
			// ..................................................... disegna gli assi
			XSetForeground(xvd.display,xvd.gc,xvd.blue[2]);
			// ............................................................... asse x
			for (i=0;i<nX;i++)
			{
				ix=xc+int((NewXax[i][0]+double(ix0-xc))*dist/(dist+NewXax[i][1])+0.5);
				iy=yc-int((NewXax[i][2]+double(yc-iy0))*dist/(dist+NewXax[i][1])+0.5)+5;
				XCPrintf(iWin,ix,iy,1,TRUE,"%d",i);
			}
			// ............................................................... asse y
			for (i=0;i<nY;i++)
			{
				ix=xc+int((NewYax[i][0]+double(ix0-xc))*dist/(dist+NewYax[i][1])+0.5);
				iy=yc-int((NewYax[i][2]+double(yc-iy0))*dist/(dist+NewYax[i][1])+0.5)+5;
				XCPrintf(iWin,ix,iy,1,TRUE,"%d",i);
			}
		}
		// ...................................................... copia sul monitor
		XCopyArea(xvd.display,xvd.backpix[0],xvd.win[0],xvd.gc,0,0,800,600,0,0);
		// ................................................................ aspetta
		musec2=clock();
		ts1.tv_nsec=10000000-(musec2-musec1)*100000; 
		if (ts1.tv_nsec>0) nanosleep(&ts1,&ts2);
	}
	// ............................................................... erase menu
	XWinPermMenu(1,"3D Plot",ParFont,ParName,ParVal,nPar,16,comandi,nComm,2);
	// ............................................................. clear memory
	delete [] ParFont;
	delete [] comandi;
	for (i=nPar-1;i>=0;i--) delete ParVal[i];
	delete [] ParVal;
	delete [] ParName;
	for (i=nTot-1;i>=0;i--)
	{
		for (j=8;j>=0;j--)
		{
			delete [] OldPrism[i][j];
			delete [] NewPrism[i][j];
		}
		delete OldPrism[i];
		delete NewPrism[i];
	}
	delete [] OldPrism;
	delete [] NewPrism;
	
	for (i=nX-1;i>=0;i--)
	{
		delete [] OldXax[i];
		delete [] NewXax[i];
	}
	for (i=nY-1;i>=0;i--)
	{
		delete [] OldYax[i];
		delete [] NewYax[i];
	}
	delete [] OldYax;
	delete [] NewYax;
	delete [] OldXax;
	delete [] NewXax;
	for (i=2;i>=0;i--) delete [] rotmat[i];
	delete [] rotmat;
	for (i=6;i>=0;i--) delete [] DataMat[i];
	delete [] DataMat;
	for (i=4;i>=0;i--) delete [] RotCosDir[i];
	delete [] RotCosDir;
	for (i=4;i>=0;i--) delete [] CosDir[i];
	delete [] CosDir;
	delete [] yy;
	for (i=4;i>=0;i--) delete face[i];
	delete [] face;	
	delete [] ind;
	// ....................................................... chiudi le finestre
	CloseXWindow();
}