/******************************************************************************
* An implementation of the marching cube algorithm (test file).		      *
*									      *
*						Gershon Elber, Dec 1992.      *
******************************************************************************/

#include <stdio.h>
#include "irit_sm.h"
#include "mrchcube.h"
#include "iritprsr.h"
#include "allocate.h"

#define INPUT_ASCII	1
#define INPUT_INTEGER	2
#define INPUT_LONG	3
#define INPUT_BYTE	4
#define INPUT_FLOAT	5
#define INPUT_DOUBLE	6


static RealType
    CubeWidth = 0.1;
static int
    DumpBinary = FALSE,
    DumpOneLayer = 0,
    DumpPolySkip = 1,
    SkipInputData = 1,
    DumpNormals = 1,
    DataWidth = 100,
    DataDepth = 100,
    DataHeight = 100,
    InputFormat = INPUT_ASCII;

static RealType GetOneScalar(void);
static MCCubeCornerScalarStruct *GetCube(CagdBType FirstTime);
static void EstimateGradient(MCCubeCornerScalarStruct *CCS);
static void DumpPolys(MCPolygonStruct *Polys);
static void DumpPolysBinary(MCPolygonStruct *Polys);

/******************************************************************************
* Test routine of Marching Cubes.					      *
* -t Threshold : the level (real number) to compute iso surface at.	      *
* -d Width Depth Height : dimension of sampled data.			      *
* -f [ascii | int | long | byte | float | double] : single datum format.      *
* -s n : dump out only every n polygons. 1 means dump all polygons.           *
* -S n : skip n rows/cols/planes in input data. 1 means skip nothing.         *
* -o n : process only layer number n.					      *
* -n : do not do normals.						      *
******************************************************************************/
void main(int argc, char **argv)
{
    MCCubeCornerScalarStruct *CCS;
    RealType Threshold = 10.5;

    while (argc > 1) {
	if (argv[1][0] == '-') {
	    switch (argv[1][1]) {
		case 't':
		    sscanf(argv[2], "%lf", &Threshold);
		    fprintf(stderr, "Threshold set at %lf\n", Threshold);
		    argc -= 2;
		    argv += 2;
		    break;
		case 'd':
		    sscanf(argv[2], "%d", &DataWidth);
		    sscanf(argv[3], "%d", &DataDepth);
		    sscanf(argv[4], "%d", &DataHeight);
		    fprintf(stderr, "Dimensions %dx%dx%d\n",
			    DataWidth, DataDepth, DataHeight);
		    argc -= 4;
		    argv += 4;
		    break;
		case 'f':
		    if (strncmp(argv[2], "ascii", 3) == 0)
			InputFormat = INPUT_ASCII;
		    else if (strncmp(argv[2], "int", 3) == 0)
			InputFormat = INPUT_INTEGER;
		    else if (strncmp(argv[2], "long", 3) == 0)
			InputFormat = INPUT_LONG;
		    else if (strncmp(argv[2], "byte", 3) == 0)
			InputFormat = INPUT_BYTE;
		    else if (strncmp(argv[2], "float", 3) == 0)
			InputFormat = INPUT_FLOAT;
		    else if (strncmp(argv[2], "double", 3) == 0)
			InputFormat = INPUT_DOUBLE;
		    else
		    {
			fprintf(stderr, "Unknown input format \"%s\"\n",
				argv[2]);
			exit(1);
		    }
		    fprintf(stderr, "Input format is set to \"%s\"\n",
			    argv[2]);
		    argc -= 2;
		    argv += 2;
		    break;
		case 's':
		    sscanf(argv[2], "%d", &DumpPolySkip);
		    fprintf(stderr, "DumpPolySkip is %d\n", DumpPolySkip);
		    argc -= 2;
		    argv += 2;
		    break;
		case 'S':
		    sscanf(argv[2], "%d", &SkipInputData);
		    fprintf(stderr, "SkipInputData is %d\n", SkipInputData);
		    CubeWidth *= SkipInputData;
		    argc -= 2;
		    argv += 2;
		    break;
		case 'n':
		    DumpNormals = 0;
		    fprintf(stderr, "Normal are not dumped\n");
		    argc--;
		    argv++;
		    break;
		case 'b':
		    DumpBinary = TRUE;
		    fprintf(stderr, "Binary data is dumped out\n");
		    argc--;
		    argv++;
		    break;
		case 'o':
		    sscanf(argv[2], "%d", &DumpOneLayer);
		    fprintf(stderr, "DumpOneLayer is %d\n", DumpOneLayer);
		    argc -= 2;
		    argv += 2;
		    break;
		default:
		    fprintf(stderr, "Unknown option %c\n", argv[1][1]);
		    exit(1);
	    }
	}
	else {
	    fprintf(stderr, "Unknown command line %s\n", argv[1]);
	    exit(1);
	}
    }

    GetCube(TRUE);

    if (DumpBinary) {
	MCPolygonStruct
	    *AllPolys = NULL;

	while (CCS = GetCube(FALSE)) {
	    MCPolygonStruct *PolyTmp,
	        *Polys = MCThresholdCube(CCS, Threshold);

	    if (Polys != NULL) {
		for (PolyTmp = Polys;
		     PolyTmp -> Pnext != NULL;
		     PolyTmp = PolyTmp -> Pnext);
		PolyTmp -> Pnext = AllPolys;
		AllPolys = Polys;
	    }
	}

	DumpPolysBinary(AllPolys);
    }
    else {
	printf("[OBJECT MC\n");
	while (CCS = GetCube(FALSE)) {
	    MCPolygonStruct
	        *Polys = MCThresholdCube(CCS, Threshold);

	    DumpPolys(Polys);
	}
	printf("]\n");
    }

    exit(0);
}

/******************************************************************************
* Read stdin and returns one real at a time.				      *
******************************************************************************/
static RealType GetOneScalar(void)
{
    short i;
    long l;
    float f;
    double d;
    RealType r;

    switch (InputFormat) {
	case INPUT_ASCII:
	    if (scanf("%lf", &r) != 1)
		return INFINITY;
	    break;
	case INPUT_INTEGER:
	    i = getchar();
	    r = getchar() * 256 + i;
	    break;
	case INPUT_LONG:
	    if (read(0, &l, 4) != 4)
		return INFINITY;
	    r = l;
	    break;
	case INPUT_BYTE:
	    if ((r = getchar()) == EOF)
		return INFINITY;
	    break;
	case INPUT_FLOAT:
	    if (read(0, &f, 4) != 4)
		return INFINITY;
	    r = f;
	    break;
	case INPUT_DOUBLE:
	    if (read(0, &d, 8) != 8)
		return INFINITY;
	    r = d;
	    break;
	default:
	    fprintf(stderr, "Input format requested not supported.\n");
	    exit(1);
	    break;
    }

    return r;
}

/******************************************************************************
* Read stdin and returns one cube at a time.				      *
******************************************************************************/
static MCCubeCornerScalarStruct *GetCube(CagdBType FirstTime)
{
    static MCCubeCornerScalarStruct
	CCS;
    static int
	ProcessedOneLayer = FALSE,
	LayerCountX = -1,
	LayerCountY = 0,
	LayerNumber = -1;
    static RealType
	*LayerOne = NULL,
	*LayerTwo = NULL;
    int i, j;
    RealType *p;

    if (FirstTime) {
	/* Initialize the CCS constant data */
	CCS.CubeDim[0] = CubeWidth;
	CCS.CubeDim[1] = CubeWidth;
	CCS.CubeDim[2] = CubeWidth;
	LayerNumber = -SkipInputData;
	LayerCountX = -1;
	LayerCountY = 0;

	if (LayerOne == NULL) {
	    LayerOne = (RealType *) IritMalloc(sizeof(RealType) * DataWidth *
					       DataDepth);
	    LayerTwo = (RealType *) IritMalloc(sizeof(RealType) * DataWidth *
					       DataDepth);
	}

	for (p = LayerTwo, i = 0; i < DataWidth * DataDepth; i++)
	    if ((*p++ = GetOneScalar()) == INFINITY)
		return NULL;
	return NULL;
    }

    if (LayerCountX == -1) {		    /* Read the next layer. */
	do {
	    LayerNumber += SkipInputData;
	    if ((DumpOneLayer > 0 && ProcessedOneLayer) ||
		LayerNumber >= DataHeight - 1) {
		IritFree((VoidPtr) LayerOne);
		IritFree((VoidPtr) LayerTwo);
		return NULL;
	    }

	    p = LayerOne;
	    LayerOne = LayerTwo;
	    LayerTwo = p;

	    for (j = 0; j < SkipInputData; j++) {
		for (p = LayerTwo, i = 0; i < DataWidth * DataDepth; i++)
		    if ((*p++ = GetOneScalar()) == INFINITY)
		        return NULL;
	    }

	    LayerCountX = 0;
	    LayerCountY = 0;

	    fprintf(stderr, "Doing layer %d\n", LayerNumber);
	}
	while (DumpOneLayer > LayerNumber);

	ProcessedOneLayer = TRUE;
    }
    

    CCS.Vrtx0Lctn[0] = LayerCountX * CubeWidth / SkipInputData;
    CCS.Vrtx0Lctn[1] = LayerCountY * CubeWidth / SkipInputData;
    CCS.Vrtx0Lctn[2] = LayerNumber * CubeWidth / SkipInputData;
    CCS.Corners[0] = LayerOne[LayerCountY * DataWidth + LayerCountX];
    CCS.Corners[1] = LayerOne[LayerCountY * DataWidth +
			      LayerCountX + SkipInputData];
    CCS.Corners[2] = LayerOne[(LayerCountY + SkipInputData) * DataWidth +
			      LayerCountX + SkipInputData];
    CCS.Corners[3] = LayerOne[(LayerCountY + SkipInputData) * DataWidth +
			      LayerCountX];
    CCS.Corners[4] = LayerTwo[LayerCountY * DataWidth + LayerCountX];
    CCS.Corners[5] = LayerTwo[LayerCountY * DataWidth +
			      LayerCountX + SkipInputData];
    CCS.Corners[6] = LayerTwo[(LayerCountY + SkipInputData) * DataWidth +
			      LayerCountX + SkipInputData];
    CCS.Corners[7] = LayerTwo[(LayerCountY + SkipInputData) * DataWidth +
			      LayerCountX];

    EstimateGradient(&CCS);

    LayerCountX += SkipInputData;
    if (LayerCountX >= DataWidth - SkipInputData) {
	LayerCountY += SkipInputData;
	LayerCountX = 0;
	if (LayerCountY >= DataDepth - SkipInputData)
	    LayerCountX = -1; /* No more data */
    }

    return &CCS;
}

/******************************************************************************
* Estimate the gradient at the eight vertices, by first order difference.     *
* Note gradient is estimated from this cube's eight vertices only.	      *
******************************************************************************/
static void EstimateGradient(MCCubeCornerScalarStruct *CCS)
{
    CCS -> GradientX[0] =
	CCS -> GradientX[1] =
	    CCS -> Corners[1] - CCS -> Corners[0];
    CCS -> GradientX[2] =
	CCS -> GradientX[3] =
	    CCS -> Corners[2] - CCS -> Corners[3];
    CCS -> GradientX[4] =
	CCS -> GradientX[5] =
	    CCS -> Corners[5] - CCS -> Corners[4];
    CCS -> GradientX[6] =
	CCS -> GradientX[7] =
	    CCS -> Corners[6] - CCS -> Corners[7];
    
    CCS -> GradientY[0] =
	CCS -> GradientY[3] =
	    CCS -> Corners[3] - CCS -> Corners[0];
    CCS -> GradientY[1] =
	CCS -> GradientY[2] =
	    CCS -> Corners[2] - CCS -> Corners[1];
    CCS -> GradientY[4] =
	CCS -> GradientY[7] =
	    CCS -> Corners[7] - CCS -> Corners[4];
    CCS -> GradientY[5] =
	CCS -> GradientY[6] =
	    CCS -> Corners[6] - CCS -> Corners[5];
    
    CCS -> GradientZ[0] =
	CCS -> GradientZ[4] =
	    CCS -> Corners[4] - CCS -> Corners[0];
    CCS -> GradientZ[1] =
	CCS -> GradientZ[5] =
	    CCS -> Corners[5] - CCS -> Corners[1];
    CCS -> GradientZ[2] =
	CCS -> GradientZ[6] =
	    CCS -> Corners[6] - CCS -> Corners[2];
    CCS -> GradientZ[3] =
	CCS -> GradientZ[7] =
	    CCS -> Corners[7] - CCS -> Corners[3];

    CCS -> HasGradient = TRUE;
}

/******************************************************************************
* Dumps the polygons to stdout (Only triangles).			      *
******************************************************************************/
static void DumpPolys(MCPolygonStruct *Polys)
{
    static int
	DumpPolySkipLcl = 1;
    char *VrtxFrmt;

    if (DumpNormals)
	VrtxFrmt = "\t[[NORMAL %9.6lg %9.6lg %9.6lg] %9.6lg %9.6lg %9.6lg]\n";
    else
	VrtxFrmt = "\t[%9.6lg %9.6lg %9.6lg]\n";

    while (Polys) {
	int i;
	MCPolygonStruct
	    *Poly = Polys;

	Polys = Polys -> Pnext;

	/* Skip polygons in the output (so we can display them - they      */
	/* are just too many of them.					   */
	if (--DumpPolySkipLcl <= 0) {
	    DumpPolySkipLcl = DumpPolySkip;

	    for (i = 2; i < Poly -> NumOfVertices - 1; i++) {
		printf("    [POLYGON 3\n");
		if (DumpNormals) {
		    printf(VrtxFrmt,
			   Poly -> N[0][0], Poly -> N[0][1], Poly -> N[0][2],
			   Poly -> V[0][0], Poly -> V[0][1], Poly -> V[0][2]);
		    printf(VrtxFrmt,
			   Poly -> N[i-1][0], Poly -> N[i-1][1], Poly -> N[i-1][2],
			   Poly -> V[i-1][0], Poly -> V[i-1][1], Poly -> V[i-1][2]);
		    printf(VrtxFrmt,
			   Poly -> N[i][0], Poly -> N[i][1], Poly -> N[i][2],
			   Poly -> V[i][0], Poly -> V[i][1], Poly -> V[i][2]);
		}
		else {
		    printf(VrtxFrmt,
			   Poly -> V[0][0], Poly -> V[0][1], Poly -> V[0][2]);
		    printf(VrtxFrmt,
			   Poly -> V[i-1][0], Poly -> V[i-1][1], Poly -> V[i-1][2]);
		    printf(VrtxFrmt,
			   Poly -> V[i][0], Poly -> V[i][1], Poly -> V[i][2]);
		}
		printf("    ]\n");
	    }
	}

	IritFree((VoidPtr) Poly);
    }
}

/******************************************************************************
* Dumps triangles to stdout as binary file.				      *
******************************************************************************/
static void DumpPolysBinary(MCPolygonStruct *Polys)
{
    static int
	DumpPolySkipLcl = 1;
    int i, j, Handler;
    IPPolygonStruct *Pl,
	*PlHead = NULL;
    IPObjectStruct *PObj;

    while (Polys) {
	IPVertexStruct *V1, *V2, *V3;
	int i;
	MCPolygonStruct
	    *Poly = Polys;

	Polys = Polys -> Pnext;

	/* Skip polygons in the output (so we can display them - they      */
	/* are just too many of them.					   */
	if (--DumpPolySkipLcl <= 0) {
	    DumpPolySkipLcl = DumpPolySkip;

	    for (i = 2; i < Poly -> NumOfVertices - 1; i++) {
		V3 = IPAllocVertex(0, 0, NULL, NULL);
		V2 = IPAllocVertex(0, 0, NULL, V3);
		V1 = IPAllocVertex(0, 0, NULL, V2);

		Pl = IPAllocPolygon(0, 0, V1, PlHead);
		PlHead = Pl;

		for (j = 0; j < 3; j++) {
		    V1 -> Coord[j] = Poly -> V[0][j];
		    V2 -> Coord[j] = Poly -> V[i - 1][j];
		    V3 -> Coord[j] = Poly -> V[i][j];
		}

		if (DumpNormals) {
		    for (j = 0; j < 3; j++) {
			V1 -> Normal[j] = Poly -> N[0][j];
			V2 -> Normal[j] = Poly -> N[i - 1][j];
			V3 -> Normal[j] = Poly -> N[i][j];
		    }
		    IP_SET_NORMAL_VRTX(V1);
		    IP_SET_NORMAL_VRTX(V2);
		    IP_SET_NORMAL_VRTX(V3);
		}
	    }
	}

	IritFree((VoidPtr) Poly);
    }

    PObj = IPAllocObject("MCH", IP_OBJ_POLY, NULL);
    PObj -> U.Pl = PlHead;

    Handler = IritPrsrOpenStreamFromFile(stdout, FALSE, TRUE, FALSE);
    IritPrsrPutObjectToHandler(Handler, PObj);
    IritPrsrCloseStream(Handler, TRUE);

    IPFreeObject(PObj);
}