// obj.cxx -- routines to handle WaveFront .obj format files.
//
// Written by Curtis Olson, started October 1997.
//
// Copyright (C) 1997 - 1998  Curtis L. Olson - curt@me.umn.edu
//
// 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 <iostream>
#include <string.h>

#include <vector>
#include "Include/compiler.h"

#ifdef NEEDNAMESPACESTD
using namespace std;
#endif

#include "obj.hxx"

#include <Math/mat3.h>
#include <Math/point3d.hxx>


typedef vector < Point3D > container3;
typedef container3::iterator iterator3;
typedef container3::const_iterator const_iterator3;


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

static int ccw_list[MAXNODES];
int ccw_list_ptr;

static int cw_list[MAXNODES];
int cw_list_ptr;

FILE *in, *out;

Point3D ref;


// 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;
    }
}


void dump_global_bounds( void ) {
    double dist_squared, radius, radius_squared;

    radius = 0.0;

    fprintf(out, "\n");

    
    iterator3 current = nodes.begin();
    iterator3 last = nodes.end();

    // skip first dummy node
    ++current;

    for ( ; current != last; ++current ) {
	dist_squared = ref.distance3Dsquared(*current);
	// cout << "node = " << *current << " dist = " << dist_squared << endl;

	if ( dist_squared > radius_squared ) {
	    radius_squared = dist_squared;
	}
    }

    radius = sqrt(radius_squared);

    fprintf( out, 
	     "gbs %.5f %.5f %.5f %.2f\n", 
	     ref.x(), ref.y(), ref.z(), radius);
}


// dump nodes
void dump_nodes( void ) {
    Point3D p;

    fprintf(out, "\n");

    iterator3 current = nodes.begin();
    iterator3 last = nodes.end();

    // skip first dummy node
    ++current;

    for ( ; current != last; ++current ) {
	p = *current - ref;
	fprintf( out, "v %.5f %.5f %.5f\n", p.x(), p.y(), p.z() );
    }
}


// dump normals
void dump_normals( void ) {
    Point3D p;

    fprintf(out, "\n");

    iterator3 current = normals.begin();
    iterator3 last = normals.end();

    // skip first dummy normal
    ++current;

    for ( ; current != last; ++current ) {
	p = *current;
	fprintf(out, "vn %.5f %.5f %.5f\n", p.x(), p.y(), p.z() );
    }
}


// dump faces
void dump_faces( void ) {
    Point3D p;
    int i, n1, n2, n3;
    double 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].x();
	ymin = ymax = nodes[n1].y();
	zmin = zmax = nodes[n1].z();

	if ( nodes[n2].x() < xmin ) { xmin = nodes[n2].x(); }
	if ( nodes[n2].x() > xmax ) { xmax = nodes[n2].x(); }
	if ( nodes[n2].y() < ymin ) { ymin = nodes[n2].y(); }
	if ( nodes[n2].y() > ymax ) { ymax = nodes[n2].y(); }
	if ( nodes[n2].z() < zmin ) { zmin = nodes[n2].z(); }
	if ( nodes[n2].z() > zmax ) { zmax = nodes[n2].z(); }

	if ( nodes[n3].x() < xmin ) { xmin = nodes[n3].x(); }
	if ( nodes[n3].x() > xmax ) { xmax = nodes[n3].x(); }
	if ( nodes[n3].y() < ymin ) { ymin = nodes[n3].y(); }
	if ( nodes[n3].y() > ymax ) { ymax = nodes[n3].y(); }
	if ( nodes[n3].z() < zmin ) { zmin = nodes[n3].z(); }
	if ( nodes[n3].z() > zmax ) { zmax = nodes[n3].z(); }

	p = Point3D( (xmin + xmax) / 2.0,
		     (ymin + ymax) / 2.0,
		     (zmin + zmax) / 2.0 );

	// calc bounding radius
	radius = p.distance3D(nodes[n1]);

	dist = p.distance3D(nodes[n2]);
	if ( dist > radius ) { radius = dist; }

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


// dump list
void dump_list(int *list, int list_ptr) {
    Point3D p;
    double xmax, xmin, ymax, ymin, zmax, zmin, dist_squared, radius_squared;
    double 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].x();
	ymin = ymax = nodes[n].y();
	zmin = zmax = nodes[n].z();
	// printf("%.2f %.2f %.2f\n", nodes[n].x(), nodes[n].y(), nodes[n].z());

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

	// calc bounding radius
	n = list[i];
	radius_squared = p.distance3Dsquared(nodes[n]);

	for ( j = i + 1; j < i + len; j++ ) {
	    n = list[j];
	    dist_squared = p.distance3Dsquared(nodes[n]);
	    if ( dist_squared > radius_squared ) {
		radius_squared = dist_squared;
	    }
	}
	radius = sqrt(radius_squared);

	// printf("radius = %.2f\n", radius);

	// dump bounding sphere and header
	fprintf(out, "bs %.2f %.2f %.2f %.2f\n", p.x(), p.y(), p.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].x() - nodes[n1].x();
    v1[1] = nodes[n2].y() - nodes[n1].y();
    v1[2] = nodes[n2].z() - nodes[n1].z();
    v2[0] = nodes[n3].x() - nodes[n1].x();
    v2[1] = nodes[n3].y() - nodes[n1].y();
    v2[2] = nodes[n3].z() - nodes[n1].z();

    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) {
    Point3D node, normal;
    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);
    }

    // push dummy records onto the lists since we start counting with "1"
    node = Point3D(0.0, 0.0, 0.0);
    nodes.push_back(node);

    normal = Point3D(0.0, 0.0, 0.0);
    normals.push_back(normal);

    // initialize other lists
    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;
    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
	    // printf("vertex = %s", line);
	    sscanf(line, "v %lf %lf %lf\n", &x, &y, &z);

	    if ( nodes.size() == 1 ) {
		// first time through set min's and max'es
		xmin = x;
		xmax = x;
		ymin = y;
		ymax = y;
		zmin = z;
		zmax = z;
	    } else {
		// update 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;		
	    }

	    node = Point3D(x, y, z);
	    nodes.push_back(node);
	    // fprintf(out, "v %.2f %.2f %.2f\n", 
	    //       node.x(), node.y(), node.z());
	} else if ( strncmp(line, "vn ", 3) == 0 ) {
	    // save vertex normals to memory and output to file
	    // printf("vertex normal = %s", line);
	    sscanf(line, "vn %lf %lf %lf\n", &x, &y, &z);
	    normal = Point3D(x, y, z);
	    normals.push_back(normal);
	} 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"
    ref = Point3D( (xmin + xmax) / 2.0,
		   (ymin + ymax) / 2.0,
		   (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.3  1999/02/01 21:09:40  curt
// Optimizations from Norman Vine.
//
// Revision 1.2  1998/10/21 14:55:55  curt
// Converted to Point3D class.
//
// Revision 1.1  1998/06/08 17:11:46  curt
// Renamed *.[ch] to *.[ch]xx
//
// Revision 1.16  1998/05/27 02:27:22  curt
// Commented out a couple of debugging messages.
//
// 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.
//