/**************************************************************************
 * obj.c -- routines to handle WaveFront .obj format files.
 *
 * Written by Curtis Olson, started October 1997.
 *
 * Copyright (C) 1997  Curtis L. Olson  - curt@infoplane.com
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 * $Id$
 * (Log is kept at end of this file)
 **************************************************************************/


#include <stdio.h>
#include <string.h>

#include "obj.h"

#include <Math/mat3.h>


/* what do ya' know, here's some global variables */
static double nodes[MAXNODES][3];
static double normals[MAXNODES][3];
static int faces[MAXNODES][3];
int ncount, vncount, fcount;

static int ccw_list[MAXNODES];
int ccw_list_ptr;

static int cw_list[MAXNODES];
int cw_list_ptr;

FILE *in, *out;

double refx, refy, refz;


/* some simple list routines */

/* reset the list */
void list_init(int *list_ptr) {
    *list_ptr = 0;
}


/* add to list */
void list_add(int *list, int *list_ptr, int node) {
    if ( *list_ptr >= MAXNODES ) {
	printf("ERROR: list overflow in list_add()\n");
	exit(-1);
    }

    list[*list_ptr] = node;
    *list_ptr += 1;

    /* printf("list pointer = %d  adding %d\n", *list_ptr, node); */
}


/* fix the cw list and append to ccw_list */
void fix_cw_list(int *list, int list_ptr) {
    int i, j, len;

    if ( list_ptr < 3 ) {
	printf("List is empty ... skipping\n");
	return;
    }

    printf("Fixing cw list, size = %d\n", list_ptr);

    i = 0;
    while ( i < list_ptr ) { 
	/* do next strip */

	/* find length */
	len = 0;
	/* scan rest of strip (until -1) */
	while ( ((i+len) < list_ptr) && (list[i+len] != -1) ) { 
	    // printf("len = %d item = %d\n", len, list[i+len] );
	    len++;
	}
	// printf("          Final length = %d\n", len);

	if ( (len % 2) != 0 ) {
	    /* if length is odd, just reverse order of nodes to reverse
	       winding */
	    if ( ccw_list_ptr ) {
		list_add(ccw_list, &ccw_list_ptr, -1);
	    }
	    for ( j = i + len - 1; j >= i; j-- ) {
		// printf(" odd -> item = %d\n", list[j] );
		list_add(ccw_list, &ccw_list_ptr, list[j]);
	    }
	} else {
	    /* if length is even, reverse order of (n-1) nodes to
	       reverse winding, and create an orphan triangle for the
	       last "nth" node */
	    if ( ccw_list_ptr ) {
		list_add(ccw_list, &ccw_list_ptr, -1);
	    }
	    for ( j = i + len - 2; j >= i; j-- ) {
		// printf(" even -> item = %d\n", list[j] );
		list_add(ccw_list, &ccw_list_ptr, list[j]);
	    }

	    // printf(" even bonus -> item = %d\n", list[i + len - 1] );
	    // printf(" even bonus -> item = %d\n", list[i + len - 2] );
	    // printf(" even bonus -> item = %d\n", list[i + len - 3] );
	    list_add(ccw_list, &ccw_list_ptr, -1);
	    list_add(ccw_list, &ccw_list_ptr, list[i + len - 3]);
	    list_add(ccw_list, &ccw_list_ptr, list[i + len - 2]);
	    list_add(ccw_list, &ccw_list_ptr, list[i + len - 1]);
	}

	i += len + 1;
    }
}


// Calculate distance between (0,0,0) and the specified point
static double calc_dist(double x, double y, double z) {
    return ( sqrt(x*x + y*y + z*z) );
}


void dump_global_bounds( void ) {
    double dist, radius;
    int i;

    radius = 0.0;

    fprintf(out, "\n");

    for ( i = 1; i < ncount; i++ ) {

	dist = calc_dist(nodes[i][0] - refx, nodes[i][1] - refy, 
			 nodes[i][2] - refz);
	// printf("node = %.2f %.2f %.2f dist = %.2f\n", 
        //        nodes[i][0], nodes[i][1], nodes[i][2],
	//        dist);

	if ( dist > radius ) {
	    radius = dist;
	}

    }

    fprintf(out, "gbs %.5f %.5f %.5f %.2f\n", refx, refy, refz, radius);
}


/* dump nodes */
void dump_nodes( void ) {
    int i;

    fprintf(out, "\n");
    for ( i = 1; i < ncount; i++ ) {
	fprintf(out, "v %.5f %.5f %.5f\n",
		nodes[i][0] - refx, nodes[i][1] - refy, nodes[i][2] - refz);
    }
}


/* dump normals */
void dump_normals( void ) {
    int i;

    fprintf(out, "\n");
    for ( i = 1; i < vncount; i++ ) {
	fprintf(out, "vn %.5f %.5f %.5f\n", 
		normals[i][0], normals[i][1], normals[i][2]);
    }
}


/* dump faces */
void dump_faces( void ) {
    int i, n1, n2, n3;
    double x, y, z, xmax, xmin, ymax, ymin, zmax, zmin, dist, radius;

    fprintf(out, "\n");
    for ( i = 1; i < fcount; i++ ) {
	n1 = faces[i][0];
	n2 = faces[i][1];
	n3 = faces[i][2];

	/* calc center of face */
	xmin = xmax = nodes[n1][0];
	ymin = ymax = nodes[n1][1];
	zmin = zmax = nodes[n1][2];

	if ( nodes[n2][0] < xmin ) { xmin = nodes[n2][0]; }
	if ( nodes[n2][0] > xmax ) { xmax = nodes[n2][0]; }
	if ( nodes[n2][1] < ymin ) { ymin = nodes[n2][1]; }
	if ( nodes[n2][1] > ymax ) { ymax = nodes[n2][1]; }
	if ( nodes[n2][2] < zmin ) { zmin = nodes[n2][2]; }
	if ( nodes[n2][2] > zmax ) { zmax = nodes[n2][2]; }

	if ( nodes[n3][0] < xmin ) { xmin = nodes[n3][0]; }
	if ( nodes[n3][0] > xmax ) { xmax = nodes[n3][0]; }
	if ( nodes[n3][1] < ymin ) { ymin = nodes[n3][1]; }
	if ( nodes[n3][1] > ymax ) { ymax = nodes[n3][1]; }
	if ( nodes[n3][2] < zmin ) { zmin = nodes[n3][2]; }
	if ( nodes[n3][2] > zmax ) { zmax = nodes[n3][2]; }

	x = (xmin + xmax) / 2.0;
	y = (ymin + ymax) / 2.0;
	z = (zmin + zmax) / 2.0;

	/* calc bounding radius */
	radius = calc_dist(nodes[n1][0] - x, nodes[n1][1] - y, 
			   nodes[n1][2] - z);

	dist = calc_dist(nodes[n2][0] - x, nodes[n2][1] - y, nodes[n2][2] - z);
	if ( dist > radius ) { radius = dist; }

	dist = calc_dist(nodes[n3][0] - x, nodes[n3][1] - y, nodes[n3][2] - z);
	if ( dist > radius ) { radius = dist; }
	
	/* output data */
	fprintf(out, "bs %.2f %.2f %.2f %.2f\n", x, y, z, radius);
	fprintf(out, "f %d %d %d\n", n1, n2, n3);
    }
}


/* dump list */
void dump_list(int *list, int list_ptr) {
    double x, y, z, xmax, xmin, ymax, ymin, zmax, zmin, dist, radius;
    int i, j, len, n;

    if ( list_ptr < 3 ) {
	printf("List is empty ... skipping\n");
	return;
    }

    printf("Dumping list, size = %d\n", list_ptr);

    i = 0;
    while ( i < list_ptr ) { 
	/* do next strip */

	if ( (i % 2) == 0 ) {
	    fprintf(out, "\nusemtl desert1\n");
	} else {
	    fprintf(out, "\nusemtl desert2\n");
	}

	/* find length of next tri strip */
	len = 0;
	/* scan rest of strip (until -1) */
	while ( ((i+len) < list_ptr) && (list[i+len] != -1) ) { 
	    // printf("len = %d item = %d\n", len, list[i+len] );
	    len++;
	}
	// printf("strip length = %d\n", len);

	/* calc center of face */
	n = list[i];
	xmin = xmax = nodes[n][0];
	ymin = ymax = nodes[n][1];
	zmin = zmax = nodes[n][2];
	printf("%.2f %.2f %.2f\n", nodes[n][0], nodes[n][1], nodes[n][2]);

	for ( j = i + 1; j < i + len; j++ ) {
	    // printf("j = %d\n", j);
	    n = list[j];
	    if ( nodes[n][0] < xmin ) { xmin = nodes[n][0]; }
	    if ( nodes[n][0] > xmax ) { xmax = nodes[n][0]; }
	    if ( nodes[n][1] < ymin ) { ymin = nodes[n][1]; }
	    if ( nodes[n][1] > ymax ) { ymax = nodes[n][1]; }
	    if ( nodes[n][2] < zmin ) { zmin = nodes[n][2]; }
	    if ( nodes[n][2] > zmax ) { zmax = nodes[n][2]; }
	    printf("%.2f %.2f %.2f\n", nodes[n][0], nodes[n][1], nodes[n][2]);
	}	    
	x = (xmin + xmax) / 2.0;
	y = (ymin + ymax) / 2.0;
	z = (zmin + zmax) / 2.0;
	printf("center = %.2f %.2f %.2f\n", x, y, z);

	/* calc bounding radius */
	n = list[i];
	radius = calc_dist(nodes[n][0] - x, nodes[n][1] - y, nodes[n][2] - z);

	for ( j = i + 1; j < i + len; j++ ) {
	    n = list[j];
	    dist = calc_dist(nodes[n][0] - x, nodes[n][1] - y, 
			     nodes[n][2] - z);
	    if ( dist > radius ) { radius = dist; }
	}
	printf("radius = %.2f\n", radius);

	/* dump bounding sphere and header */
	fprintf(out, "bs %.2f %.2f %.2f %.2f\n", x, y, z, radius);
	fprintf(out, "t %d %d %d\n", list[i], list[i+1], list[i+2]);
	/* printf("t %d %d %d\n", list[i], list[i+1], list[i+2]); */
	i += 3;

	/* dump rest of strip (until -1) */
	while ( (i < list_ptr) && (list[i] != -1) ) { 
	    fprintf(out, "q %d", list[i]);
	    i++;
	    if ( (i < list_ptr) && (list[i] != -1) ) { 
		fprintf(out, " %d", list[i]);
		i++;
	    }
	    fprintf(out, "\n");
	}

	i++;
    }
}


/* Check the direction the current triangle faces, compared to it's
 * pregenerated normal.  Returns the dot product between the target
 * normal and actual normal.  If the dot product is close to 1.0, they
 * nearly match.  If the are close to -1.0, the are nearly
 * opposite. */
double check_cur_face(int n1, int n2, int n3) {
    double v1[3], v2[3], approx_normal[3], dot_prod, temp;

    /* check for the proper rotation by calculating an approximate
     * normal and seeing if it is close to the precalculated normal */
    v1[0] = nodes[n2][0] - nodes[n1][0];
    v1[1] = nodes[n2][1] - nodes[n1][1];
    v1[2] = nodes[n2][2] - nodes[n1][2];
    v2[0] = nodes[n3][0] - nodes[n1][0];
    v2[1] = nodes[n3][1] - nodes[n1][1];
    v2[2] = nodes[n3][2] - nodes[n1][2];

    MAT3cross_product(approx_normal, v1, v2);
    MAT3_NORMALIZE_VEC(approx_normal,temp);
    dot_prod = MAT3_DOT_PRODUCT(normals[n1], approx_normal);

    /* not first triangle */
    /* if ( ((dot_prod < -0.5) && !is_backwards) ||
	 ((dot_prod >  0.5) && is_backwards) ) {
	printf("    Approx normal = %.2f %.2f %.2f\n", approx_normal[0], 
	       approx_normal[1], approx_normal[2]);
	printf("    Dot product = %.4f\n", dot_prod);
    } */
    /* angle = acos(dot_prod); */
    /* printf("Normal ANGLE = %.3f rads.\n", angle); */

    return(dot_prod);
}


/* Load a .obj file */
void obj_fix(char *infile, char *outfile) {
    char line[256];
    double dot_prod;
    int first, n1, n2, n3, n4;
    double x, y, z, xmax, xmin, ymax, ymin, zmax, zmin;
    int is_ccw;

    if ( (in = fopen(infile, "r")) == NULL ) {
	printf("Cannot open file: %s\n", infile);
	exit(-1);
    }

    if ( (out = fopen(outfile, "w")) == NULL ) {
	printf("Cannot open file: %s\n", outfile);
	exit(-1);
    }

    list_init(&ccw_list_ptr);
    list_init(&cw_list_ptr);

    /* I start counting at one because that is how the triangle
       program refers to nodes and normals */
    first = 1;
    ncount = 1;
    vncount = 1;
    fcount = 1;

    printf("Reading file:  %s\n", infile);

    while ( fgets(line, 250, in) != NULL ) {
	if ( line[0] == '#' ) {
	    /* pass along the comments verbatim */
	    fprintf(out, "%s", line);
	} else if ( strlen(line) <= 1 ) {
	    /* don't pass along empty lines */
	    // fprintf(out, "%s", line);
	} else if ( strncmp(line, "v ", 2) == 0 ) {
	    /* save vertex to memory and output to file */
            if ( ncount < MAXNODES ) {
                /* printf("vertex = %s", line); */
                sscanf(line, "v %lf %lf %lf\n", &x, &y, &z);
		nodes[ncount][0] = x;
		nodes[ncount][1] = y;
		nodes[ncount][2] = z;

		/* first time through set min's and max'es */
		if ( ncount == 1 ) {
		    xmin = x;
		    xmax = x;
		    ymin = y;
		    ymax = y;
		    zmin = z;
		    zmax = z;
		}
    
		/* keep track of min/max vertex values */
		if ( x < xmin ) xmin = x;
		if ( x > xmax ) xmax = x;
		if ( y < ymin ) ymin = y;
		if ( y > ymax ) ymax = y;
		if ( z < zmin ) zmin = z;
		if ( z > zmax ) zmax = z;		

		// fprintf(out, "v %.2f %.2f %.2f\n", 
		//       nodes[ncount][0], nodes[ncount][1], nodes[ncount][2]);
		ncount++;
            } else {
                printf("Read too many nodes ... dying :-(\n");
                exit(-1);
            }
	} else if ( strncmp(line, "vn ", 3) == 0 ) {
	    /* save vertex normals to memory and output to file */
            if ( vncount < MAXNODES ) {
                /* printf("vertex normal = %s", line); */
                sscanf(line, "vn %lf %lf %lf\n", 
                       &normals[vncount][0], &normals[vncount][1], 
                       &normals[vncount][2]);
		// fprintf(out, "vn %.4f %.4f %.4f\n", normals[vncount][0], 
		//	normals[vncount][1], normals[vncount][2]);
                vncount++;
            } else {
                printf("Read too many vertex normals ... dying :-(\n");
                exit(-1);
            }
	} else if ( line[0] == 't' ) {
	    /* starting a new triangle strip */

	    printf("Starting a new triangle strip\n");

	    n1 = n2 = n3 = n4 = 0;

	    printf("new tri strip = %s", line);
	    sscanf(line, "t %d %d %d %d\n", &n1, &n2, &n3, &n4);

	    /* special cases to handle bugs in our beloved tri striper */
	    if ( (n1 == 4) && (n2 == 2) && (n3 == 2) && (n4 == 1) ) {
		n2 = 3;
	    }
	    if ( (n1 == 3) && (n2 == 1) && (n3 == 1) && (n4 == 0) ) {
		n3 = 4;
	    }

	    dot_prod = check_cur_face(n1, n2, n3);
	    if ( dot_prod < 0.0 ) {
		/* this stripe is backwards (CW) */
		is_ccw = 0;
		printf(" -> Starting a backwards stripe\n");
	    } else {
		/* this stripe is normal (CCW) */
		is_ccw = 1;
	    }

	    if ( is_ccw ) {
		if ( ccw_list_ptr ) {
		    list_add(ccw_list, &ccw_list_ptr, -1);
		}

		list_add(ccw_list, &ccw_list_ptr, n1);
		list_add(ccw_list, &ccw_list_ptr, n2);
		list_add(ccw_list, &ccw_list_ptr, n3);
	    } else {
		if ( cw_list_ptr ) {
		    list_add(cw_list, &cw_list_ptr, -1);
		}

		list_add(cw_list, &cw_list_ptr, n1);
		list_add(cw_list, &cw_list_ptr, n2);
		list_add(cw_list, &cw_list_ptr, n3);
	    }

	    if ( n4 > 0 ) {
		if ( is_ccw ) {
		    list_add(ccw_list, &ccw_list_ptr, n4);
		} else {
		    list_add(cw_list, &cw_list_ptr, n4);
		}
	    }
	} else if ( line[0] == 'f' ) {
	    if ( fcount < MAXNODES ) {
		/* pass along the unoptimized faces verbatim */
		sscanf(line, "f %d %d %d\n", &n1, &n2, &n3);
		faces[fcount][0] = n1;
		faces[fcount][1] = n2;
		faces[fcount][2] = n3;

		fcount++;
	    } else {
		printf("Read too many unoptimized faces ... dying :-(\n");
                exit(-1);
	    }
 
	    // fprintf(out, "%s", line);
	} else if ( line[0] == 'q' ) {
	    /* continue a triangle strip */
	    n1 = n2 = 0;

	    /* printf("continued tri strip = %s ", line); */
	    sscanf(line, "q %d %d\n", &n1, &n2);

	    if ( is_ccw ) {
		list_add(ccw_list, &ccw_list_ptr, n1);
	    } else {
		list_add(cw_list, &cw_list_ptr, n1);
	    }

	    if ( n2 > 0 ) {
		if ( is_ccw ) {
		    list_add(ccw_list, &ccw_list_ptr, n2);
		} else {
		    list_add(cw_list, &cw_list_ptr, n2);
		}
	    }
	} else {
	    printf("Unknown line in %s = %s\n", infile, line);
	}
    }

    /* reference point is the "center" */
    refx = (xmin + xmax) / 2.0;
    refy = (ymin + ymax) / 2.0;
    refz = (zmin + zmax) / 2.0;

    /* convert the cw_list to ccw add append to ccw_list */
    fix_cw_list(cw_list, cw_list_ptr);

    dump_global_bounds();
    dump_nodes();
    dump_normals();
    if ( fcount > 1 ) {
	dump_faces();
    }

    dump_list(ccw_list, ccw_list_ptr);

    fclose(in);
    fclose(out);
}


/* $Log$
/* Revision 1.15  1998/05/24 02:47:47  curt
/* For each strip, specify a default material property and calculate a center
/* and bounding sphere.
/*
 * Revision 1.14  1998/05/23 15:19:49  curt
 * Output more digits after the decimal place.
 *
 * Revision 1.13  1998/05/20 20:55:19  curt
 * Fixed arbitrary polygon winding problem here so all tristrips are passed
 * to runtime simulator with a consistant counter clockwise winding.
 *
 * Revision 1.12  1998/05/16 13:11:26  curt
 * Fixed an off by one error in node, normal, and face counters.
 *
 * Revision 1.11  1998/04/27 15:59:24  curt
 * Fixed an off by one error.
 *
 * Revision 1.10  1998/04/27 03:33:11  curt
 * Code now calculates a center reference points and outputs everything
 * relative to that.  This is useful in the rendering engine to keep everything
 * close to (0, 0, 0) where we can avoid many GLfloat precision problems.
 *
 * Revision 1.9  1998/04/18 04:01:03  curt
 * Now use libMath rather than having local copies of math routines.
 *
 * Revision 1.8  1998/04/08 23:19:37  curt
 * Adopted Gnu automake/autoconf system.
 *
 * Revision 1.7  1998/03/19 02:51:41  curt
 * Added special case handling to compensate for bugs in our beloved tri striper
 *
 * Revision 1.6  1998/03/03 15:36:12  curt
 * Tweaks for compiling with g++
 *
 * Revision 1.5  1998/03/03 03:37:03  curt
 * Cumulative tweaks.
 *
 * Revision 1.4  1998/01/31 00:41:25  curt
 * Made a few changes converting floats to doubles.
 *
 * Revision 1.3  1998/01/19 19:51:07  curt
 * A couple final pre-release tweaks.
 *
 * Revision 1.2  1998/01/09 23:03:12  curt
 * Restructured to split 1deg x 1deg dem's into 64 subsections.
 *
 * Revision 1.1  1997/12/08 19:28:54  curt
 * Initial revision.
 *
 */