/*
	Computation of Form factor for Full Radiosity Solution.
*/
#include <stdio.h>
#include <math.h>
#ifdef __WATCOMC__
#include <string.h>
#else
#include <memory.h>
#endif
#include <malloc.h>

#include "GraphicsGems.h"
#include "data_structure.h"
#include "objects.h"
#include "vox.h"

#include "raddecl.h"

extern int hemicube_resolution;
extern int verbose_flag;
int ProjectionByRaytracing=0;

static double **TopFace_delta_form_factor;
static double **SideFace_delta_form_factor;
struct Hemicube{
        short int *topface;
        short int *rightface;
        short int *leftface;
        short int *frontface;
        short int *backface;
};
static double **F; /* Form Factors */
static double **A; /* N x N Matrix */
static double **I;  /* Computed Intensity */
static double **deltaI;  /* Delta Intensity for progressive radiosity */
static double *B;
static struct Hemicube hemicube;

#define Iter_Threshold 0.0001
#define VeryHigh 1.e5

static double **alloc_matrix(rows,cols)
int rows,cols;
{
	char *s="Cannot Allocate Matrix.";
	double **m;
	int i;
	m = (double **)malloc(rows*sizeof(double *));
	if (m==NULL) error(s);
	for (i=0;i<rows;i++) {
		m[i] = (double *)malloc(cols*sizeof(double));
		if (m[i] == NULL) error(s);
	}
	return(m);
}
static void allocate_hemicube()
{
	TopFace_delta_form_factor=
		alloc_matrix(hemicube_resolution,hemicube_resolution);
	SideFace_delta_form_factor=
		alloc_matrix(hemicube_resolution/2,hemicube_resolution);
	hemicube.topface=(short int *)malloc(
		hemicube_resolution*hemicube_resolution
		*sizeof(short int));
	hemicube.rightface=(short int *)malloc(
                hemicube_resolution/2*hemicube_resolution
                *sizeof(short int));
	hemicube.leftface=(short int *)malloc(
                hemicube_resolution/2*hemicube_resolution
                *sizeof(short int));
	hemicube.frontface=(short int *)malloc(
                hemicube_resolution/2*hemicube_resolution
                *sizeof(short int));
	hemicube.backface=(short int *)malloc(
                hemicube_resolution/2*hemicube_resolution
                *sizeof(short int));
}
static void free_matrix(m,rows,cols)
double **m;
int rows,cols;
{
	int i;
	for (i=0;i<rows;i++) free(m[i]);
	free(m);
}
static void deallocate_hemicube()
{
	free_matrix(TopFace_delta_form_factor,
		hemicube_resolution,hemicube_resolution);
	free_matrix(SideFace_delta_form_factor,
		hemicube_resolution/2,hemicube_resolution);
	free(hemicube.topface);
	free(hemicube.rightface);
	free(hemicube.leftface);
	free(hemicube.frontface);
	free(hemicube.backface);
}

/*
        compute_delta_form_factors : For the faces of the hemicube.
                Side faces are symmetric, so only one computation
                is enough.
*/
static compute_delta_form_factors()
{
        int i,j;
        double x, y, y2, delta;
        double factor;
        double z,zfactor,z2;
        double halfdelta;
        double delta_area_by_pi;

	delta  = 2.0/hemicube_resolution;
	halfdelta =delta/2.0;
	delta_area_by_pi=delta*delta/PI;

        for(i=0,y=(-1.0+halfdelta); i < hemicube_resolution; i++,y+=delta){
                y2 = y*y;
                for (j=0,x=(-1.0+halfdelta); j < hemicube_resolution; j++,x+=delta){
                        factor=1.0/(x*x+y2+1.0);
                        TopFace_delta_form_factor[i][j]=
                                delta_area_by_pi*factor*factor;
                }
        }
        for (i=0,z=halfdelta; i < hemicube_resolution/2; i++,z+=delta){
                z2 = z*z; zfactor = z*delta_area_by_pi;
                for(j=0,y=(-1.0+halfdelta); j < hemicube_resolution; j++,y+=delta){
                        factor=1.0/(y*y+z2+1.0);
                        SideFace_delta_form_factor[i][j]=zfactor*factor*factor;
                }
        }
}

static void increment_F(rows,raster,f,inc)
int rows;
short int *raster;
double f[];
double **inc;
{
	int i,j,k,patchnum;
	for(i=k=0;i<rows;i++)
		for(j=0;j<hemicube_resolution;j++,k++)
			if (patchnum=raster[k])
				f[patchnum-1] += inc[i][j];
}

static void extract_F_i(i)
int i;
{
	int nn;
       	for(nn=0;nn<total_surface_grid_points;nn++)F[i][nn]=0.;
	increment_F(hemicube_resolution,hemicube.topface,
		F[i],TopFace_delta_form_factor);
	increment_F(hemicube_resolution/2,hemicube.rightface,
		F[i],SideFace_delta_form_factor);
	increment_F(hemicube_resolution/2,hemicube.leftface,
		F[i],SideFace_delta_form_factor);
	increment_F(hemicube_resolution/2,hemicube.frontface,
		F[i],SideFace_delta_form_factor);
	increment_F(hemicube_resolution/2,hemicube.backface,
		F[i],SideFace_delta_form_factor);
}

static init_fullmatrix_radiosity()
{
	/* Allocate space for form factor */
	F = alloc_matrix(total_surface_grid_points,total_surface_grid_points);
	A = alloc_matrix(total_surface_grid_points,total_surface_grid_points);
	I = alloc_matrix(color_channels,total_surface_grid_points);
	B = (double *)malloc(total_surface_grid_points*sizeof(double));
	allocate_hemicube();
	compute_delta_form_factors();
}
static deinit_fullmatrix_radiosity()
{
	/* Free space allocated for FormFactor */
	free_matrix(F,total_surface_grid_points,total_surface_grid_points);
	free_matrix(A,total_surface_grid_points,total_surface_grid_points);
	free_matrix(I,color_channels,total_surface_grid_points);
	free(B);
	deallocate_hemicube();
}

static init_progressive_radiosity()
{
	/* Allocate space for form factor */
	F = alloc_matrix(1,total_surface_grid_points);
	I = alloc_matrix(color_channels,total_surface_grid_points);
	deltaI = alloc_matrix(color_channels,total_surface_grid_points);
	allocate_hemicube();
	compute_delta_form_factors();
}
static deinit_progressive_radiosity()
{
	free_matrix(F,1,total_surface_grid_points);
	free_matrix(I,color_channels,total_surface_grid_points);
	free_matrix(deltaI,color_channels,total_surface_grid_points);
	deallocate_hemicube();
}

static void from_specular_surface(ray,objectnum,t,u,v,n)
Ray *ray;
int objectnum;
double t,u,v;
int *n;
{
	Vector3 N;
        /* .... set the new ray */
        point_on_line(ray->start,ray->direction,t,ray->start);
        N = ofunc[object[objectnum].surface_geometry_type].get_surface_normal
                (object[objectnum].object_specific_structure,u,v);
        mirror_reflection(&N,&(ray->direction),&(object[objectnum]));
        (ray->number)++;
	*n=get_diffuse_surface(ray);
}
static void from_diffuse_surface(ray,objectnum,t,u,v,n)
Ray *ray;
int objectnum;
double t,u,v;
int *n;
{
	int uindex=grid_u(objectnum,u);
	int vindex=grid_v(objectnum,v);
	*n=(object[objectnum].start_grid_index+
		gridindex(objectnum,uindex,vindex)+1);
}
static void from_other_surface(ray,objectnum,t,u,v,n)
Ray *ray;
int objectnum;
double t,u,v;
int *n;
{
        error("Not Implemented.");
}
static void (*get_patchnum[MAX_REFLECTION_TYPE])()={
        from_diffuse_surface,
        from_specular_surface,
        from_other_surface
};

static int get_diffuse_surface(ray)
Ray *ray;
{
	double t_entry,t_exit;
	int HitBoundingBox();
        if (HitBoundingBox(&(volume_grid.extent.min),&(volume_grid.extent.max),
			&(ray->start),&(ray->direction),&t_entry,&t_exit)){
                int get_nearest_object_in_voxel();
                Vlist *create_vlist();
                Vlist *vlist,*templist;
                double t,u,v;
        	int nearest_object= UNDEFINED;
                vlist=templist=create_vlist(
                        &(ray->start),&(ray->direction),t_entry,t_exit
                );
                while (templist!=NULL){
                        Voxel *vox = volume_grid.voxels+templist->voxnum;
                        if(vox->nobjects)
				nearest_object=get_nearest_object_in_voxel
				(ray,templist->t_near,templist->t_far,vox,&t,&u,&v);
                        if (nearest_object!= UNDEFINED) break;
                        templist=templist->next;
                }
                purge_vlist(vlist);
                if (nearest_object!= UNDEFINED){
			int n;
			get_patchnum[object[nearest_object].surface_reflection_type]
			(ray,nearest_object,t,u,v,&n);
			return(n);
		}
	}
	return(0);
}
/*
int dump_buffer(rows,cols,b)
int rows,cols;
char *b;
{
	int i,j,k;
	fprintf(stderr,"\nDumping buffer\n");
	for(i=k=0;i<rows;i++){
		fprintf(stderr,"%02d.....",i);
		for(j=0;j<cols;j++,k++)
			fprintf(stderr,"%1d",b[k]);
		fprintf(stderr,"\n");
	}
}
int dump_raster(rows,cols,b)
int rows,cols;
short int *b;
{
	int i,j,k;
	fprintf(stderr,"\nDumping raster\n");
	for(i=k=0;i<rows;i++){
		fprintf(stderr,"%02d.....",i);
		for(j=0;j<cols;j++,k++)
			fprintf(stderr,"%1d",b[k]);
		fprintf(stderr,"\n");
	}
}
*/
static long project(C,N,U,V,rows,raster)
Point3 *C;
Vector3 *N,*U,*V;
int rows;
short int *raster;
{
	static long RayNumber=0;
	Ray ray;
	Vector3 Dir,delta_U,delta_V,half_delta_U,half_delta_V;
	int i,j,k;

	half_delta_U=delta_U = *U;
	V3Scale(&delta_U,2./hemicube_resolution);
	V3Scale(&half_delta_U,1./hemicube_resolution);
	half_delta_V=delta_V = *V;
	V3Scale(&delta_V,2./hemicube_resolution);
	V3Scale(&half_delta_V,1./hemicube_resolution);
	V3Add(N,V,&Dir); V3Sub(&Dir,U,&Dir);
	V3Add(&Dir,&half_delta_U,&Dir); V3Sub(&Dir,&half_delta_V,&Dir);
	ray.number=total_rays+RayNumber;
#if defined (DEBUG)
	ray.intersections=0;
#endif
	ray.path_length= -1;
	if (ProjectionByRaytracing)
	for(i=k=0;i<rows;i++){
		Vector3 dir;
		dir = Dir;
		for(j=0;j<hemicube_resolution;j++,k++){
			ray.start = *C;
			ray.direction=dir;
			V3Normalize(&(ray.direction));
			raster[k]=get_diffuse_surface(&ray);
			(ray.number)++;
			V3Add(&dir,&delta_U,&dir);
		}
		V3Sub(&Dir,&delta_V,&Dir);
	}
	else {
	double *zbuffer;
	char *surface_buffer;
	Matrix4 view_matrix;
	int n=rows*hemicube_resolution;
	extern void (*geometry_specific_scan_convert[MAX_GEOMETRY_TYPE])();
	/* Clean up the raster */
	memset((char *)raster,0,n*sizeof(short int));
	/* Set up a Z-buffer of equal dimension */
	zbuffer=(double *)malloc(n*sizeof(double));
	if (zbuffer == NULL)
	  error("Cannot allocate Z-buffer memory for scan conversion.");
	for(i=0;i<n;i++)zbuffer[i]=VeryHigh;
	surface_buffer=(char *)calloc(n,1);
	if (surface_buffer == NULL)
	  error("Cannot allocate surface-buffer memory for scan conversion.");
	/* Compute View transform matrix */
	view_transform(C,U,V,N,rows,hemicube_resolution,&view_matrix);
	for(i=0;i<number_objects;i++){
		/* Project each object on the raster */
		geometry_specific_scan_convert[
			object[i].surface_geometry_type
		](i,N,&view_matrix,rows,hemicube_resolution,
			zbuffer,surface_buffer,raster);
	}
	free((char *)zbuffer);
/*
dump_raster(rows,hemicube_resolution,raster);
dump_buffer(rows,hemicube_resolution,surface_buffer);
*/
	for(i=k=0;i<rows;i++)
	for(j=0;j<hemicube_resolution;j++,k++)
		if ((surface_buffer[k]-1)>DIFFUSE){
			ray.start = *C;
			ray.direction.x=Dir.x+j*delta_U.x-i*delta_V.x;
			ray.direction.y=Dir.y+j*delta_U.y-i*delta_V.y;
			ray.direction.z=Dir.z+j*delta_U.z-i*delta_V.z;
			V3Normalize(&(ray.direction));
			raster[k]=get_diffuse_surface(&ray);
			(ray.number)++;
		}
	free(surface_buffer);
	}
	return(RayNumber = (ray.number-total_rays));
}

static long project_on_hemicube(c,m)
Point3 *c;
Matrix4 *m;
{
	Vector3 X,Y,Z,Xneg,Yneg,Zneg;
	long n;

	X.x = 1.; X.y=X.z = 0.;	
	X = transform_vector(&X,m); V3Normalize(&X); Xneg=X; V3Negate(&Xneg);
	Y.y = 1.; Y.x=Y.z = 0.;	
	Y = transform_vector(&Y,m); V3Normalize(&Y); Yneg=Y; V3Negate(&Yneg);
	Z.z = 1.; Z.x=Z.y = 0.;	
	Z = transform_vector(&Z,m); V3Normalize(&Z); Zneg=Z; V3Negate(&Zneg);
	/* Project on Top Face */
	project(c,&Z,&X,&Y,hemicube_resolution,hemicube.topface);
	/* Project on Left Face */
	project(c,&Xneg,&Y,&Z,hemicube_resolution/2,hemicube.leftface);
	/* Project on Right Face */
	project(c,&X,&Yneg,&Z,hemicube_resolution/2,hemicube.rightface);
	/* Project on Front Face */
	project(c,&Yneg,&Xneg,&Z,hemicube_resolution/2,hemicube.frontface);
	/* Project on Back Face */
	n=project(c,&Y,&X,&Z,hemicube_resolution/2,hemicube.backface);
	return(n);
}

static hemicube_method(i,j,k,n,nrays)
int i;	/* Object num */
int j;  /* Vertical Grid */
int k;  /* Horizontal grid */
int n;  /* Patch num */
long *nrays;
{
	Matrix4 m;
	Vector3 N;
	Point3 centre;
	if (object[i].surface_reflection_type==MIRROR)return;
	centre = ofunc[object[i].surface_geometry_type].
		compute_center_point_in_the_grid_cell(i,j,k);
	N=ofunc[object[i].surface_geometry_type].
		get_surface_normal
		(object[i].object_specific_structure,
		(k+0.5)/object[i].grid_h_reso,
		(j+0.5)/object[i].grid_v_reso);
	m  = ofunc[object[i].surface_geometry_type].
		get_local_matrix(
		&N,object[i].object_specific_structure);
	*nrays=project_on_hemicube(&centre,&m);
}
static long compute_F_ij()
{
	int i,j,k,n;
	long nrays;
	if (verbose_flag)fprintf(stderr,"From Factor:");
	for (i=n=0;i<number_objects; i++){
		for(j=0;j<object[i].grid_v_reso;j++)
		for(k=0;k<object[i].grid_h_reso;k++,n++){
			hemicube_method(i,j,k,n,&nrays);
			extract_F_i(n);	
			if (verbose_flag)fprintf(stderr,".");
		}
	}
	if (verbose_flag)fprintf(stderr,"\n");
	return(nrays);
}

static double get_emission_color(onum,channel)
int onum,channel;
{
	int i;
	for(i=0;i<light_sources;i++)
		if (source[i].oindex == onum)
			return(source[i].strength[channel]);
	return(0.);
}

static void solve(channel,color)
int channel;
double color[];
{
	int i,j;
	int m,n;
	int diagonally_dominant();

        for (i=m=0; i < number_objects; i++){
		int grid_points = object[i].grid_v_reso*object[i].grid_h_reso;
		double rho = object[i].reflectance[channel];
		double c = get_emission_color(i,channel);
		for (j=0; j < grid_points; j++,m++){
			for(n=0;n<total_surface_grid_points;n++)
				A[m][n] = -F[m][n]*rho;
			A[m][m] += 1.;
			B[m] = color[m] = c/object[i].grid[j].area;
		}
	}
	if (diagonally_dominant(total_surface_grid_points,A)){
		fprintf(stderr,"Warning : Not Diagonally Dominant.\n");
		fprintf(stderr,
		      "Check for the reflectance of the surfaces.\n");
		fprintf(stderr,"Iteration may not converge.\n");
		fprintf(stderr,"So a predefined number of iterations will be tried.\n");
	}
	gauss_iter(total_surface_grid_points,A,color,B,Iter_Threshold);
}

static set_intensity(diff_flag,fl)
int diff_flag;
FILE *fl;
{
	int i,j,k,m;
	double max_intensity=0.,max_cum_intensity=0.;
        if((fl != NULL)&& verbose_output_flag){
                int rows,cols;
                for (i=m=0; i < number_objects; i++){
                        fprintf(fl,"intensity output for Surface : %d\n",i);
			for(rows=0;rows<object[i].grid_v_reso;rows++)
			for(cols=0;cols<object[i].grid_h_reso;cols++,m++){
                                Point3 centre;
                                centre = ofunc[object[i].
					surface_geometry_type].
				   compute_center_point_in_the_grid_cell
					 (i,rows,cols);
				fprintf(fl,"%g\t%g\t%g",centre.x,
					centre.y,centre.z);
                                for(k=0; k < color_channels; k++)
                                        fprintf(fl,"\t%g", I[k][m]);
                                fprintf(fl,"\n");
                        }
                }
        }
	/* Normalize */
	for(k=0;k<color_channels;k++)
	for(i=0;i<total_surface_grid_points;i++)
		if (I[k][i] > max_intensity) max_intensity=I[k][i];	
	for(k=0;k<color_channels;k++)
	for(i=0;i<total_surface_grid_points;i++)
		I[k][i] /= max_intensity;

	for (i=m=0; i < number_objects; i++){
		double cum_area=0.,cum_intensity[MAXCHANNELS];
		int grid_points = object[i].grid_v_reso*object[i].grid_h_reso;
		Surface_Grid_structure *g=object[i].grid;
		for(k=0; k<color_channels;k++) cum_intensity[k] = 0.;
		for (j=0; j < grid_points; j++,m++,g++){
			cum_area += g->area;
			for(k=0;k<color_channels;k++){
				cum_intensity[k] += I[k][m]*g->area;
				g->normalized_flux_density[k] = ((diff_flag)
					? fabs(I[k][m]-g->normalized_flux_density[k])
					:I[k][m]);
			}
		}
		for(k=0;k<color_channels;k++){
			cum_intensity[k]/=cum_area;
			object[i].normalized_flux_density[k] = cum_intensity[k];
			if (cum_intensity[k] > max_cum_intensity )
				max_cum_intensity = cum_intensity[k];
		}
	}
	for(i=0;i<number_objects; i++)
		for(k=0;k<color_channels;k++)
			object[i].normalized_flux_density[k]/=max_cum_intensity;
}


long fullmatrix_radiosity(diff_flag,fl)
int diff_flag;
FILE *fl;
{
	int i;
	long number=0;
	init_fullmatrix_radiosity();
	number=compute_F_ij();
	for(i=0;i<color_channels;i++)
		solve(i,I[i]);
	set_intensity(diff_flag,fl);
	write_out_intensity(fl);
	deinit_fullmatrix_radiosity();
	return(number);
}

long progressive_radiosity(maxpasses,diff_flag,fl)
int maxpasses;
int diff_flag;
FILE *fl;
{
	int i,j,k,npasses=0;
	long nrays;
	int max_h_grid,max_v_grid,maxpatch,maxobject= UNDEFINED;
	double maxdelta=0.;

	init_progressive_radiosity();
	for(k=0;k<color_channels;k++)
		for(i=0;i<total_surface_grid_points;i++)
			deltaI[k][i]=I[k][i]=0.;
	for(i=0;i<light_sources;i++){
		int index=source[i].oindex,n,nn;
		double cumarea=0.;
		for(j=nn=0;j<object[index].grid_v_reso;j++)
		for(k=0;k<object[index].grid_h_reso;k++,nn++)
			cumarea += object[index].grid[nn].area;
		n=object[index].start_grid_index;
		for(j=nn=0;j<object[index].grid_v_reso;j++)
		for(k=0;k<object[index].grid_h_reso;k++,nn++,n++){
			int l;
			double area=object[index].grid[nn].area;
			for(l=0;l<color_channels;l++){
			   I[l][n]=source[i].strength[l]/cumarea;
			   deltaI[l][n] = I[l][n]*area;
			   if (maxdelta < deltaI[l][n]){
				maxdelta=deltaI[l][n];
				max_h_grid=k;
				max_v_grid=j;
				maxpatch=n;
				maxobject=index;
			   }
			}
		}
	}
	if (maxpasses < 0) maxpasses=total_surface_grid_points;
	/* Distribute */
	while (maxpatch!= UNDEFINED){
		int l,n,nn;
		double area,increment[MAXCHANNELS];
		hemicube_method(maxobject,
			max_v_grid,max_h_grid,
			maxpatch,&nrays);
		extract_F_i(0);	
		area=object[maxobject].grid[max_v_grid*
				object[maxobject].grid_h_reso+
				max_h_grid].area;
		for(l=0;l<color_channels;l++){
		   increment[l]=deltaI[l][maxpatch];
		   deltaI[l][maxpatch]=0.;
		}
		maxobject= UNDEFINED; maxdelta=0.;
		for(i=n=0;i<number_objects;i++){
			for(j=nn=0;j<object[i].grid_v_reso;j++)
			for(k=0;k<object[i].grid_h_reso;k++,n++,nn++)
			if (object[i].surface_reflection_type !=DIFFUSE)
				continue;
			else{
				int l;
				for(l=0;l<color_channels;l++){
				   deltaI[l][n]+=
					object[i].reflectance[l]*
					increment[l]*F[0][n];
				   I[l][n]+=
					object[i].reflectance[l]*
					increment[l]*F[0][n]/
						object[i].grid[nn].area;
				   if(maxdelta < (deltaI[l][n])){
					maxdelta=deltaI[l][n];
					max_h_grid=k;
					max_v_grid=j;
					maxpatch=n;
					maxobject=i;
				   }
				}
			}
		}
		npasses++;
		if (verbose_flag)fprintf(stderr,".");
		if (npasses >= maxpasses) break;
	}
	if (verbose_flag)
		fprintf(stderr,"\n%d Distribution steps.\n",npasses);
	set_intensity(diff_flag,fl);
	write_out_intensity(fl);
	deinit_progressive_radiosity();
	return(nrays);
}