You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
413 lines
12 KiB
413 lines
12 KiB
/***************************************************************************
|
|
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 );
|
|
}
|
|
|
|
|
|
|
|
|
|
|