/**************************************************************************/
/*                                                                        */
/*	COPYRIGHT NOTICE:  This software was written and/or developed     */
/*      by John A. Hipp, Ph.D. at the Orthopaedic Biomechanics Lab        */
/*      of Beth Israel Hospital and Harvard Medical School.               */
/*                                                                        */
/**************************************************************************/
/*                                                                        */
/*      It may NOT BE USED without the EXPRESS consent of John Hipp.      */
/*                                                                        */
/**************************************************************************/
/* generates a 3D lattice for use with verifying 3D morphology soft */

/* takes as input:

	3d_grid nx sx ny sy nz sz >! raw.byt

		where 	nx is the number of equally spaced struts long the x axis
				sx is the width in pixels of each strut
				etc
   program also calulates the predicted and actual volume fraction of 
	the solid phase as well as the surface area
*/


#include <stdio.h>

#define asize 128

FILE *fpo;
char *calloc();


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

double dx,dy,dist, avfract, pvfract, sarea;
long tsize, scnt;
int nstruts, astrut;
int i, nx, sx, ny, sy, nz, sz, xx, yy, zz;
int xgap, ygap, zgap, xp, yp, zp, cur_pt, ix, iy, iz;
int zdim, ydim, xdim;
unsigned char background, grid;
unsigned char *array;

/* check for proper command line paramaters */
if(argc !=7) {
	fprintf(stderr,"ERROR: wrong command line format\n\n");
	fprintf(stderr,"use: 3d_grid nx sx ny sy nz sz >! outfile\n");
	exit(-1);
	}

nx  = atoi(argv[1]);
sx = atoi(argv[2]);
ny = atoi(argv[3]);
sy = atoi(argv[4]);
nz = atoi(argv[5]);
sz = atoi(argv[6]);

tsize = asize * asize * asize;
array = (unsigned char *)calloc(tsize, sizeof(unsigned char));
if(array == NULL) {
	fprintf(stderr,"ERROR: could not allocate memory for array \n\n");
	exit(-1);
	}

fpo = stdout;

xgap = (asize - nx * sx) / (nx + 1);
if(xgap < 0) {
	fprintf(stderr,"ERROR: parameters give no gap between struts alng x-axis\n");
	exit(-1);
	}
ygap = (asize - ny * sy) / (ny + 1);
if(xgap < 0) {
	fprintf(stderr,"ERROR: parameters give no gap between struts alng y-axis\n");
	exit(-1);
	}
zgap = (asize - nz * sz) / (nz + 1);
if(xgap < 0) {
	fprintf(stderr,"ERROR: parameters give no gap between struts alng z-axis\n");
	exit(-1);
	}

background = 200;
grid = 20;

/* initialize array to background */

for(zdim = 0; zdim < asize; zdim++) {
	for(ydim = 0; ydim < asize; ydim++) {
		for(xdim = 0; xdim < asize; xdim++) {
			*(array + (zdim * asize + ydim) * asize + xdim) = background;
		}
	}
}

/* draw the struts perpendicular to the xy plane */

for(yy=0; yy < ny; yy++) {   /* loop over each strut along the y axis */
	for(iy=0; iy < sy; iy++) { /* loop over the pixels accross each strut */
		yp = (yy + 1) * ygap + yy * sy + iy;
	
		for(xx=0; xx < nx; xx++) { /* loop over each strut along the x-axis */
			for(ix=0; ix < sx; ix++) { /* loop over the pixels accross each strut */
				xp = (xx + 1) * xgap + xx * sx + ix;

				for(zp=0;zp<asize;zp++) { /* struts are continuous in z-direction */
					
					cur_pt = (zp * asize + yp) * asize + xp;
					*(array + cur_pt) = grid;

					}
	
				} /* end loop over the pixels accross each strut */
			} /* end loop over each strut along the x-axis */

		} /* end loop over the pixels accross each strut */
	} /* end loop over each strut along the y axis */


/* draw the struts perpendicular to the xz plane */

for(zz=0; zz < nz; zz++) {   /* loop over each strut along the z axis */
	for(iz=0; iz < sz; iz++) { /* loop over the pixels accross each z strut */
		zp = (zz + 1) * zgap + zz * sz + iz;
	
		for(xx=0; xx < nx; xx++) { /* loop over each strut along the x-axis */
			for(ix=0; ix < sx; ix++) { /* loop over the pixels accross each strut */
				xp = (xx + 1) * xgap + xx * sx + ix;

				for(yp=0;yp<asize;yp++) { /* struts are continuous in y-direction */
					
					cur_pt = (zp * asize + yp) * asize + xp;
					*(array + cur_pt) = grid;
	
					}
	
				} /* end loop over the pixels accross each x strut */
			} /* end loop over each strut along the x-axis */

		} /* end loop over the pixels accross each z strut */
	} /* end loop over each strut along the z axis */

/* draw the struts perpendicular to the yz plane */

for(zz=0; zz < nz; zz++) {   /* loop over each strut along the z axis */
	for(iz=0; iz < sz; iz++) { /* loop over the pixels accross each z strut */
		zp = (zz + 1) * zgap + zz * sz + iz;
	
		for(yy=0; yy < ny; yy++) { /* loop over each strut along the y-axis */
			for(iy=0; iy < sy; iy++) { /* loop over the pixels accross each strut */
				yp = (yy + 1) * ygap + yy * sy + iy;

				for(xp=0;xp<asize;xp++) { /* struts are continuous in x-direction */
					
					cur_pt = (zp * asize + yp) * asize + xp;
					*(array + cur_pt) = grid;

					}
	
				} /* end loop over the pixels accross each y strut */
			} /* end loop over each strut along the y-axis */

		} /* end loop over the pixels accross each z strut */
	} /* end loop over each strut along the z axis */

/* finally! - writ-out the array */

i = fwrite(array,sizeof(unsigned char), tsize,fpo);

fclose(fpo);

/* calculate the actual volume fraction */

scnt = 0;
for(i=0;i<tsize;i++) if(*(array + i) == grid) scnt++;
avfract = (double)scnt / (double)tsize;
fprintf(stderr,"Actual volume fraction for entire cube is %lf\n",avfract);

/* calculate what the volume fraction should be */

/* first consider the struts perpindicular to the xy plane */

nstruts = nx * ny;  /* total # struts */
astrut = sx * sy;   /* x-sect area of each strut */
scnt = nstruts * astrut * asize;

/* now consider the struts perp to the xz plane - don't include intersect with xy struts */

astrut = sx * sz; /* x-sect area of each strut */
nstruts = nx * nz;
scnt += nstruts * astrut * (asize - ny * sy);

/* now consider the struts perp to the yz plane - don't include intersect with xy struts */

astrut = sy * sz; /* x-sect area of each strut */
nstruts = ny * nz;
scnt += nstruts * astrut * (asize - nx * sx);

pvfract = (double)scnt / (double)tsize;
fprintf(stderr,"Predicted volume fraction for entire cube is %lf\n",pvfract);
pvfract -= avfract;
fprintf(stderr,"Predicted - actual is %lf\n",pvfract);

/* now calculate the surface area of the struts */

/* first consider the struts perpindicular to the xy plane */

nstruts = nx * ny;  /* total # struts */
astrut = 2 * sx + 2 * sy;   /* perimeter of each strut */
scnt = nstruts * astrut * (asize - nz * sz);

/* now consider the struts perp to the xz plane  */

astrut = 2 * sx + 2 * sz; /* perimeter of each strut */
nstruts = nx * nz;
scnt += nstruts * astrut * (asize - ny * sy);

/* now consider the struts perp to the yz plane  */

astrut = 2 * sy + 2 * sz; /* perimeter of each strut */
nstruts = ny * nz;
scnt += nstruts * astrut * (asize - nx * sx);

sarea = (double)scnt / (double)tsize;

fprintf(stderr,"Total surface area (pixels**2) = %ld\n",scnt);
fprintf(stderr,"Surface area / Total Volume = %lf\n",sarea);
sarea = (double)scnt / ((double)tsize * avfract);
fprintf(stderr,"Surface area / Volume struts = %lf\n",sarea);

}