/***************************************************************************
                                ksdelunay.cpp
                             -------------------                                         
    begin                :
    copyright            : (C) 2001 by Kamil Dobkowski
    email                : kamildobk@poczta.onet.pl                                     
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.                                   * 
 *                                                                         *
 ***************************************************************************/

#include"mpdelunay.h"
#include<math.h>
#include<algo.h>
#include<qobject.h>

//---------------------------------------------------------------------------------------------//


MPDelunay::MPDelunay( MPSymbolList *args, int columnFrom, int columnTo, const char *id )
: MPSymbol( args, columnFrom, columnTo, id )
 {
  m_triangle_buff = 0;
  m_triangle_buff_len = 0;
  m_point_data = NULL;
  m_triangle_data = NULL;
  m_edge_data = NULL;
  m_edges = 0;
  m_edge_buffer_capacity = 0;
  m_point_data = NULL;
  m_points = 0;
  m_triangles = 0;
  m_active_triangles = 0;
  m_p1 = NULL;
  m_p2 = NULL;
  m_p3 = NULL;
  m_result = NULL;
  m_complete = NULL;
  m_error = QString::null;
 }

//---------------------------------------------------------------------------------------------//

MPDelunay::~MPDelunay()
 {
  delete[] m_triangle_buff;
 }

//---------------------------------------------------------------------------------------------//


void MPDelunay::checkArgs( MPError& error )
 {
  stop();
  delete m_triangle_buff; m_triangle_buff = 0;
  m_rows = 0;
  m_cols = 0;
  m_error = QString::null;

  for ( int i=0; i<m_args->count(); i++ ) m_args->at(i)->checkArgs(error);
  MPSymbol::checkArgs( error );
  if ( m_args->count() != 2 ) { error.setWrongNumberOfArguments( m_args->count(), 2, this ); return; }

  triangulation( m_args->at(0), m_args->at(1), true );

  if ( !m_error.isEmpty() )  {
	error.setError( m_error,this );
	m_rows = 0;
	m_cols = 0;
	m_error = QString::null;
	delete m_triangle_buff; m_triangle_buff = 0;
	return;
	}
 }

//---------------------------------------------------------------------------------------------//

void MPDelunay::stop()
 {
  delete[] m_complete;
  m_complete = NULL;

  delete[] m_point_data;
  m_point_data = NULL;
  m_points = 0;

  delete[] m_triangle_data;
  m_triangle_data = NULL;
  m_active_triangles = 0;
  m_triangles = 0;

  delete[] m_edge_data;
  m_edge_data = NULL;
  m_edge_buffer_capacity = 0;
  m_edges = 0;

  delete[] m_result;
  m_result = NULL;

  m_p1 = NULL;
  m_p2 = NULL;
  m_p3 = NULL;
 }

//-------------------------------------------------------------------------//

double MPDelunay::value( int row , int col )
 {
  return m_triangle_buff[row+col*m_triangle_buff_len];
 }


//---------------------------------------------------------------------------------------------//

void MPDelunay::set_error( const QString& message )
 {
  m_error += message + "\n";
 }

//---------------------------------------------------------------------------------------------//

void MPDelunay::triangulation( MPSymbol *xColumn, MPSymbol *yColumn, bool useLargeBuffers )
 {
  m_error = QString::null;
  if ( xColumn->rows() != yColumn->rows() ) {
	set_error( QObject::tr("X and Y must have the same number of rows !") );
	stop();
	return;
	}
  if ( xColumn->cols() != 1 ||
       yColumn->cols() != 1  ) {
	set_error( QObject::tr("X and Y must be column vectors !") );
	stop();
	return;
	};

  if ( xColumn->rows() < 3 ) {
	set_error( QObject::tr("There must be at least three points !") );
	stop();
	return;
	}

  /////////////////////////
  // prepare input data  //
  /////////////////////////

  // copy all points to the internal buffer
  // internal buffer will hold also 3 vertices of bounding super triangle
  m_points = xColumn->rows();
  m_point_data = new point_data_t[m_points+3];
  for( int i=0; i<m_points; i++ ) {
	 m_point_data[i].index = i;
	 m_point_data[i].x = xColumn->value( i, 0 );
	 m_point_data[i].y = yColumn->value( i, 0 );;
	}

  // presort values using x coordinates - simple bubble sort
  sort( m_point_data, m_point_data+m_points );

  point_data_t *new_end = unique( m_point_data, m_point_data+m_points );

  if ( new_end != m_point_data+m_points ) {
	m_points = new_end - m_point_data;
	}

  if ( m_points < 3 ) {
	set_error( QObject::tr("There must be at least three different points !") );
	stop();
	return;
	}

  // calculate the bounding box
  double x_min = m_point_data[0].x;
  double y_min = m_point_data[0].y;
  double x_max = m_point_data[0].x;
  double y_max = m_point_data[0].y;
  for( int i=1; i<m_points; i++ ) {
	 if ( x_min > m_point_data[i].x ) x_min = m_point_data[i].x;
	 if ( y_min > m_point_data[i].y ) y_min = m_point_data[i].y;
	 if ( x_max < m_point_data[i].x ) x_max = m_point_data[i].x;
	 if ( y_max < m_point_data[i].y ) y_max = m_point_data[i].y;	
	}

  // construct the bounding super triangle
  // anti-clockwise(?) order
  // this will be the default order of all triangles
  m_point_data[m_points+0].x =  x_min - 1.0;
  m_point_data[m_points+0].y =  y_min - 1.0;
  m_point_data[m_points+1].x =  x_min - 1.0;
  m_point_data[m_points+1].y =  y_max + ( y_max-y_min ) + 1.0;
  m_point_data[m_points+2].x =  x_max + ( x_max-x_min ) + 1.0;
  m_point_data[m_points+2].y =  y_min - 1.0;
  m_point_data[m_points+0].index = -1;
  m_point_data[m_points+1].index = -1;
  m_point_data[m_points+2].index = -1;

  // allocate nessesary space for the output
  // max. number of triangles is 2*(m_points+3)
  m_triangles = 0;
  int triangles_max = (m_points+3)*2;
  m_result = new long[3*triangles_max]; // buffer for output vertices
  m_p1 = m_result;
  m_p2 = m_p1 + triangles_max;
  m_p3 = m_p2 + triangles_max;
  m_complete = new bool[triangles_max];

  if ( useLargeBuffers ) {
	m_triangle_data = new triangle_data_t[triangles_max];
	}

  // add the first triangle
  add_triangle( m_points, m_points+1, m_points+2 );

  /////////////////////////
  // start triangulation //
  /////////////////////////

  // for each point
  for( int point_nr=0; point_nr<m_points; point_nr++ ) {
	double curr_x = m_point_data[point_nr].x;
        double curr_y = m_point_data[point_nr].y;
	//////////////////////////////////////////////
	// complete edge buffer for the given point //
	//////////////////////////////////////////////
	// for each triangle
	for( int triangle_nr=0; triangle_nr<m_triangles; triangle_nr++ ) {
	    // all completed triangles are put at the end, so if
	    // completed triangle appeared -  all remaining triangles
            // should be not considered
	    if ( m_complete[triangle_nr] ) break;
	    // test if the current point is inside triangle
	    if ( point_in_circum_circle(curr_x,curr_y,triangle_nr) ) {
		// add edges of the current triangle to the buffer
		int p1 = m_p1[triangle_nr];
		int p2 = m_p2[triangle_nr];
		int p3 = m_p3[triangle_nr];
                add_edge( p1, p2 );
		add_edge( p2, p3 );
		add_edge( p3, p1 );
		// remove the current triangle from the triangle list
                del_triangle( triangle_nr );
		// after removal there is a different triangle at the current index
		// so we have to test this index once again
		// do not increment triangle_nr in the next iteration
		triangle_nr--;
		}
	     // triangle was marked as complete - put it on the end
	     if ( m_complete[triangle_nr] ) {
		 put_at_end( triangle_nr );
 		 // now at the current index is a different triangle
		 // so we have to test this index once again
		 triangle_nr--;
		}	
	   }

	/////////////////////////////////////////////////////////
	// make a new triangles using edges in the edge buffer //
        /////////////////////////////////////////////////////////
        // tag multiple edges
        for ( int j=0; j<m_edges-1; j++ )
		for ( int k=j+1; k<m_edges; k++ ) {
			edge_data_t *edge1 = &m_edge_data[j];
			edge_data_t *edge2 = &m_edge_data[k];
			if ( (edge1->p1 == edge2->p2 && edge1->p2 == edge2->p1) ||
			     (edge1->p1 == edge2->p1 && edge1->p2 == edge2->p2)  ) {
				edge1->p1 = -1;	
				edge1->p2 = -1;	
				edge2->p1 = -1;
				edge2->p2 = -1;
		   		}
			}

	// form a new triangle
	for ( int j=0; j<m_edges; j++) {		
		edge_data_t *edge = &m_edge_data[j];
		if ( edge->p1 >= 0 && edge->p2 >= 0 ) add_triangle( edge->p1, edge->p2, point_nr );			
		}

	// clear the edge buffer
	m_edges = 0;

	// emit sig progress
	//if ( (point_nr%progress_step) == 0 ) {
	//	//emit sigProgress( 30+point_nr*70/m_points, &cancel );
	//	}
	}

  // remove triangles with the supertriangle vertices
  // ( supertriangle vertices have index greater then m_points-1 )
  m_active_triangles = m_triangles;
  for( int triangle_nr=0; triangle_nr<m_triangles; triangle_nr++ ) {
      if ( m_p1[triangle_nr] >= m_points ||
	   m_p2[triangle_nr] >= m_points ||
           m_p3[triangle_nr] >= m_points  ) {
		del_triangle( triangle_nr );
		triangle_nr --;
		}
     };

  // change index to point buffer with an index to the apriopriate row.
  for( int i=0; i<m_triangles; i++ ) {
	m_p1[i] = m_point_data[m_p1[i]].index;
	m_p2[i] = m_point_data[m_p2[i]].index;
	m_p3[i] = m_point_data[m_p3[i]].index;
	}

  // free all buffers but m_result
  delete m_triangle_buff; m_triangle_buff = NULL;
  m_rows = m_triangles;
  m_cols = 3;
  m_triangle_buff_len = triangles_max;
  m_triangle_buff = m_result;
  m_result = NULL;
  stop();
 }

//---------------------------------------------------------------------------------------------//

bool MPDelunay::point_in_circum_circle( double xp, double yp, int triangle )
 {
  if ( m_complete[triangle] ) return false;

  // calculate circum circle
  double xc,yc,r2;
  if ( !m_triangle_data || m_triangle_data[triangle].r2 < 0.0 ) {
	 point_data_t *p1 = &m_point_data[m_p1[triangle]];
	 point_data_t *p2 = &m_point_data[m_p2[triangle]];
	 point_data_t *p3 = &m_point_data[m_p3[triangle]];
	 double x1 = p1->x;
	 double y1 = p1->y;
	 double x2 = p2->x;
	 double y2 = p2->y;
	 double x3 = p3->x;
	 double y3 = p3->y;

         static const double EPSILON = 10e-100;
         double m1,m2,mx1,mx2,my1,my2;
   	 if (fabs(y1-y2) < EPSILON && fabs(y2-y3) < EPSILON) {
       		//set_error("Error when calculating circum circle ( points too close ? ) ");
       		return false;
   		}

   	if ( fabs(y2-y1) < EPSILON ) {
       		m2 = - (x3-x2) / (y3-y2);
       		mx2 = (x2 + x3) / 2.0;
       		my2 = (y2 + y3) / 2.0;
       		xc = (x2 + x1) / 2.0;
       		yc = m2 * (xc - mx2) + my2;
    	} else if ( fabs(y3-y2) < EPSILON ) {
       		m1 = - (x2-x1) / (y2-y1);
       		mx1 = (x1 + x2) / 2.0;
       		my1 = (y1 + y2) / 2.0;
       		xc = (x3 + x2) / 2.0;
       		yc = m1 * (xc - mx1) + my1;
    	} else {
       		m1 = - (x2-x1) / (y2-y1);
       		m2 = - (x3-x2) / (y3-y2);
       		mx1 = (x1 + x2) / 2.0;
       		mx2 = (x2 + x3) / 2.0;
       		my1 = (y1 + y2) / 2.0;
       		my2 = (y2 + y3) / 2.0;
       		xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
       		yc = m1 * (xc - mx1) + my1;
    		}

        double dx  = x2 - xc;
        double dy  = y2 - yc;
    	r2 = dx*dx + dy*dy;

   	// remember circum circle
  	if ( m_triangle_data ) {
		triangle_data_t *triangle_data = &m_triangle_data[triangle];
		triangle_data->r2 = r2;
		triangle_data->xc = xc;
		triangle_data->yc = yc;
		}

	}
    // circum circle was calculated befor - use buffered values
    else if ( m_triangle_data ) {
	triangle_data_t *triangle_data = &m_triangle_data[triangle];
	r2 = triangle_data->r2;
	xc = triangle_data->xc;
	yc = triangle_data->yc;	
	}

  // check if this triangle should be considered in further processing
  // we have presorted points, so we are sure that further points will have
  // greater x value than this point, so all further points will be outside
  // of the circum circle if the current point is outside.
  // if ( xc+r < xp ) then won't consider this triangle in further processing
  // xc+r<xp; xc+sqrt(r2)<xp; sqrt(r2)<xp-xc; r2<(xp-xc)^2 && xp>xc
  double dx2 = ( xp - xc ) * ( xp - xc );
  double dy2 = ( yp - yc ) * ( yp - yc );
  if ( xp > xc && r2 < dx2 ) {  m_complete[triangle] = true; return false; }

  // check if point is inside circum circle
  double d2 = dx2 + dy2;

  return ( d2 <= r2 );
 }