976 lines
24 KiB
C#
976 lines
24 KiB
C#
/*
|
|
* Created by SharpDevelop.
|
|
* User: Burhan
|
|
* Date: 17/06/2014
|
|
* Time: 11:30 م
|
|
*
|
|
* To change this template use Tools | Options | Coding | Edit Standard Headers.
|
|
*/
|
|
|
|
/*
|
|
* The author of this software is Steven Fortune. Copyright (c) 1994 by AT&T
|
|
* Bell Laboratories.
|
|
* Permission to use, copy, modify, and distribute this software for any
|
|
* purpose without fee is hereby granted, provided that this entire notice
|
|
* is included in all copies of any software which is or includes a copy
|
|
* or modification of this software and in all copies of the supporting
|
|
* documentation for such software.
|
|
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
|
* WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
|
|
* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
|
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
|
*/
|
|
|
|
/*
|
|
* This code was originally written by Stephan Fortune in C code. I, Shane O'Sullivan,
|
|
* have since modified it, encapsulating it in a C++ class and, fixing memory leaks and
|
|
* adding accessors to the Voronoi Edges.
|
|
* Permission to use, copy, modify, and distribute this software for any
|
|
* purpose without fee is hereby granted, provided that this entire notice
|
|
* is included in all copies of any software which is or includes a copy
|
|
* or modification of this software and in all copies of the supporting
|
|
* documentation for such software.
|
|
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
|
* WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
|
|
* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
|
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
|
*/
|
|
|
|
/*
|
|
* Java Version by Zhenyu Pan
|
|
* Permission to use, copy, modify, and distribute this software for any
|
|
* purpose without fee is hereby granted, provided that this entire notice
|
|
* is included in all copies of any software which is or includes a copy
|
|
* or modification of this software and in all copies of the supporting
|
|
* documentation for such software.
|
|
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
|
* WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
|
|
* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
|
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
|
*/
|
|
|
|
/*
|
|
* C# Version by Burhan Joukhadar
|
|
*
|
|
* Permission to use, copy, modify, and distribute this software for any
|
|
* purpose without fee is hereby granted, provided that this entire notice
|
|
* is included in all copies of any software which is or includes a copy
|
|
* or modification of this software and in all copies of the supporting
|
|
* documentation for such software.
|
|
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
|
* WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
|
|
* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
|
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
|
*/
|
|
|
|
using Microsoft.Xna.Framework;
|
|
using System;
|
|
using System.Collections.Generic;
|
|
|
|
namespace Voronoi2
|
|
{
|
|
/// <summary>
|
|
/// Description of Voronoi.
|
|
/// </summary>
|
|
public class Voronoi
|
|
{
|
|
// ************* Private members ******************
|
|
double borderMinX, borderMaxX, borderMinY, borderMaxY;
|
|
int siteidx;
|
|
double xmin, xmax, ymin, ymax, deltax, deltay;
|
|
int nvertices;
|
|
int nedges;
|
|
int nsites;
|
|
Site[] sites;
|
|
Site bottomsite;
|
|
int sqrt_nsites;
|
|
double minDistanceBetweenSites;
|
|
int PQcount;
|
|
int PQmin;
|
|
int PQhashsize;
|
|
Halfedge[] PQhash;
|
|
|
|
const int LE = 0;
|
|
const int RE = 1;
|
|
|
|
int ELhashsize;
|
|
Halfedge[] ELhash;
|
|
Halfedge ELleftend, ELrightend;
|
|
List<GraphEdge> allEdges;
|
|
|
|
|
|
// ************* Public methods ******************
|
|
// ******************************************
|
|
|
|
// constructor
|
|
public Voronoi ( double minDistanceBetweenSites )
|
|
{
|
|
siteidx = 0;
|
|
sites = null;
|
|
|
|
allEdges = null;
|
|
this.minDistanceBetweenSites = minDistanceBetweenSites;
|
|
}
|
|
|
|
/**
|
|
*
|
|
* @param xValuesIn Array of X values for each site.
|
|
* @param yValuesIn Array of Y values for each site. Must be identical length to yValuesIn
|
|
* @param minX The minimum X of the bounding box around the voronoi
|
|
* @param maxX The maximum X of the bounding box around the voronoi
|
|
* @param minY The minimum Y of the bounding box around the voronoi
|
|
* @param maxY The maximum Y of the bounding box around the voronoi
|
|
* @return
|
|
*/
|
|
// تستدعى هذه العملية لإنشاء مخطط فورونوي
|
|
public List<GraphEdge> generateVoronoi ( double[] xValuesIn, double[] yValuesIn, double minX, double maxX, double minY, double maxY )
|
|
{
|
|
sort(xValuesIn, yValuesIn, xValuesIn.Length);
|
|
|
|
// Check bounding box inputs - if mins are bigger than maxes, swap them
|
|
double temp = 0;
|
|
if ( minX > maxX )
|
|
{
|
|
temp = minX;
|
|
minX = maxX;
|
|
maxX = temp;
|
|
}
|
|
if ( minY > maxY )
|
|
{
|
|
temp = minY;
|
|
minY = maxY;
|
|
maxY = temp;
|
|
}
|
|
|
|
borderMinX = minX;
|
|
borderMinY = minY;
|
|
borderMaxX = maxX;
|
|
borderMaxY = maxY;
|
|
|
|
siteidx = 0;
|
|
voronoi_bd ();
|
|
return allEdges;
|
|
}
|
|
|
|
|
|
/*********************************************************
|
|
* Private methods - implementation details
|
|
********************************************************/
|
|
|
|
private void sort ( double[] xValuesIn, double[] yValuesIn, int count )
|
|
{
|
|
sites = null;
|
|
allEdges = new List<GraphEdge>();
|
|
|
|
nsites = count;
|
|
nvertices = 0;
|
|
nedges = 0;
|
|
|
|
double sn = (double)nsites + 4;
|
|
sqrt_nsites = (int) Math.Sqrt ( sn );
|
|
|
|
// Copy the inputs so we don't modify the originals
|
|
double[] xValues = new double[count];
|
|
double[] yValues = new double[count];
|
|
for (int i = 0; i < count; i++)
|
|
{
|
|
xValues[i] = xValuesIn[i];
|
|
yValues[i] = yValuesIn[i];
|
|
}
|
|
sortNode ( xValues, yValues, count );
|
|
}
|
|
|
|
private void qsort ( Site[] sites )
|
|
{
|
|
List<Site> listSites = new List<Site>( sites.Length );
|
|
for ( int i = 0; i < sites.Length; i++ )
|
|
{
|
|
listSites.Add ( sites[i] );
|
|
}
|
|
|
|
listSites.Sort ( new SiteSorterYX () );
|
|
|
|
// Copy back into the array
|
|
for (int i=0; i < sites.Length; i++)
|
|
{
|
|
sites[i] = listSites[i];
|
|
}
|
|
}
|
|
|
|
private void sortNode ( double[] xValues, double[] yValues, int numPoints )
|
|
{
|
|
nsites = numPoints;
|
|
sites = new Site[nsites];
|
|
xmin = xValues[0];
|
|
ymin = yValues[0];
|
|
xmax = xValues[0];
|
|
ymax = yValues[0];
|
|
|
|
for ( int i = 0; i < nsites; i++ )
|
|
{
|
|
sites[i] = new Site();
|
|
sites[i].coord.setPoint ( xValues[i], yValues[i] );
|
|
sites[i].sitenbr = i;
|
|
|
|
if ( xValues[i] < xmin )
|
|
xmin = xValues[i];
|
|
else if ( xValues[i] > xmax )
|
|
xmax = xValues[i];
|
|
|
|
if ( yValues[i] < ymin )
|
|
ymin = yValues[i];
|
|
else if ( yValues[i] > ymax )
|
|
ymax = yValues[i];
|
|
}
|
|
|
|
qsort ( sites );
|
|
deltax = xmax - xmin;
|
|
deltay = ymax - ymin;
|
|
}
|
|
|
|
private Site nextone ()
|
|
{
|
|
Site s;
|
|
if ( siteidx < nsites )
|
|
{
|
|
s = sites[siteidx];
|
|
siteidx++;
|
|
return s;
|
|
}
|
|
return null;
|
|
}
|
|
|
|
private Edge bisect ( Site s1, Site s2 )
|
|
{
|
|
double dx, dy, adx, ady;
|
|
Edge newedge;
|
|
|
|
newedge = new Edge();
|
|
|
|
newedge.reg[0] = s1;
|
|
newedge.reg[1] = s2;
|
|
|
|
newedge.ep [0] = null;
|
|
newedge.ep[1] = null;
|
|
|
|
dx = s2.coord.x - s1.coord.x;
|
|
dy = s2.coord.y - s1.coord.y;
|
|
|
|
adx = dx > 0 ? dx : -dx;
|
|
ady = dy > 0 ? dy : -dy;
|
|
newedge.c = (double)(s1.coord.x * dx + s1.coord.y * dy + (dx * dx + dy* dy) * 0.5);
|
|
|
|
if ( adx > ady )
|
|
{
|
|
newedge.a = 1.0;
|
|
newedge.b = dy / dx;
|
|
newedge.c /= dx;
|
|
}
|
|
else
|
|
{
|
|
newedge.a = dx / dy;
|
|
newedge.b = 1.0;
|
|
newedge.c /= dy;
|
|
}
|
|
|
|
newedge.edgenbr = nedges;
|
|
nedges++;
|
|
|
|
return newedge;
|
|
}
|
|
|
|
private void makevertex ( Site v )
|
|
{
|
|
v.sitenbr = nvertices;
|
|
nvertices++;
|
|
}
|
|
|
|
private bool PQinitialize ()
|
|
{
|
|
PQcount = 0;
|
|
PQmin = 0;
|
|
PQhashsize = 4 * sqrt_nsites;
|
|
PQhash = new Halfedge[ PQhashsize ];
|
|
|
|
for ( int i = 0; i < PQhashsize; i++ )
|
|
{
|
|
PQhash [i] = new Halfedge();
|
|
}
|
|
return true;
|
|
}
|
|
|
|
private int PQbucket ( Halfedge he )
|
|
{
|
|
int bucket;
|
|
|
|
bucket = (int) ((he.ystar - ymin) / deltay * PQhashsize);
|
|
if ( bucket < 0 )
|
|
bucket = 0;
|
|
if ( bucket >= PQhashsize )
|
|
bucket = PQhashsize - 1;
|
|
if ( bucket < PQmin )
|
|
PQmin = bucket;
|
|
|
|
return bucket;
|
|
}
|
|
|
|
// push the HalfEdge into the ordered linked list of vertices
|
|
private void PQinsert ( Halfedge he, Site v, double offset )
|
|
{
|
|
Halfedge last, next;
|
|
|
|
he.vertex = v;
|
|
he.ystar = (double)(v.coord.y + offset);
|
|
last = PQhash [ PQbucket (he) ];
|
|
|
|
while
|
|
(
|
|
(next = last.PQnext) != null
|
|
&&
|
|
(he.ystar > next.ystar || (he.ystar == next.ystar && v.coord.x > next.vertex.coord.x))
|
|
)
|
|
{
|
|
last = next;
|
|
}
|
|
|
|
he.PQnext = last.PQnext;
|
|
last.PQnext = he;
|
|
PQcount++;
|
|
}
|
|
|
|
// remove the HalfEdge from the list of vertices
|
|
private void PQdelete ( Halfedge he )
|
|
{
|
|
Halfedge last;
|
|
|
|
if (he.vertex != null)
|
|
{
|
|
last = PQhash [ PQbucket (he) ];
|
|
while ( last.PQnext != he )
|
|
{
|
|
last = last.PQnext;
|
|
}
|
|
|
|
last.PQnext = he.PQnext;
|
|
PQcount--;
|
|
he.vertex = null;
|
|
}
|
|
}
|
|
|
|
private bool PQempty ()
|
|
{
|
|
return ( PQcount == 0 );
|
|
}
|
|
|
|
private Point PQ_min ()
|
|
{
|
|
Point answer = new Point ();
|
|
|
|
while ( PQhash[PQmin].PQnext == null )
|
|
{
|
|
PQmin++;
|
|
}
|
|
|
|
answer.x = PQhash[PQmin].PQnext.vertex.coord.x;
|
|
answer.y = PQhash[PQmin].PQnext.ystar;
|
|
return answer;
|
|
}
|
|
|
|
private Halfedge PQextractmin ()
|
|
{
|
|
Halfedge curr;
|
|
|
|
curr = PQhash[PQmin].PQnext;
|
|
PQhash[PQmin].PQnext = curr.PQnext;
|
|
PQcount--;
|
|
|
|
return curr;
|
|
}
|
|
|
|
private Halfedge HEcreate(Edge e, int pm)
|
|
{
|
|
Halfedge answer = new Halfedge();
|
|
answer.ELedge = e;
|
|
answer.ELpm = pm;
|
|
answer.PQnext = null;
|
|
answer.vertex = null;
|
|
|
|
return answer;
|
|
}
|
|
|
|
private bool ELinitialize()
|
|
{
|
|
ELhashsize = 2 * sqrt_nsites;
|
|
ELhash = new Halfedge[ELhashsize];
|
|
|
|
for (int i = 0; i < ELhashsize; i++)
|
|
{
|
|
ELhash[i] = null;
|
|
}
|
|
|
|
ELleftend = HEcreate ( null, 0 );
|
|
ELrightend = HEcreate ( null, 0 );
|
|
ELleftend.ELleft = null;
|
|
ELleftend.ELright = ELrightend;
|
|
ELrightend.ELleft = ELleftend;
|
|
ELrightend.ELright = null;
|
|
ELhash[0] = ELleftend;
|
|
ELhash[ELhashsize - 1] = ELrightend;
|
|
|
|
return true;
|
|
}
|
|
|
|
private Halfedge ELright( Halfedge he )
|
|
{
|
|
return he.ELright;
|
|
}
|
|
|
|
private Halfedge ELleft( Halfedge he )
|
|
{
|
|
return he.ELleft;
|
|
}
|
|
|
|
private Site leftreg( Halfedge he )
|
|
{
|
|
if (he.ELedge == null)
|
|
{
|
|
return bottomsite;
|
|
}
|
|
return (he.ELpm == LE ? he.ELedge.reg[LE] : he.ELedge.reg[RE]);
|
|
}
|
|
|
|
private void ELinsert( Halfedge lb, Halfedge newHe )
|
|
{
|
|
newHe.ELleft = lb;
|
|
newHe.ELright = lb.ELright;
|
|
(lb.ELright).ELleft = newHe;
|
|
lb.ELright = newHe;
|
|
}
|
|
|
|
/*
|
|
* This delete routine can't reclaim node, since pointers from hash table
|
|
* may be present.
|
|
*/
|
|
private void ELdelete( Halfedge he )
|
|
{
|
|
(he.ELleft).ELright = he.ELright;
|
|
(he.ELright).ELleft = he.ELleft;
|
|
he.deleted = true;
|
|
}
|
|
|
|
/* Get entry from hash table, pruning any deleted nodes */
|
|
private Halfedge ELgethash( int b )
|
|
{
|
|
Halfedge he;
|
|
if (b < 0 || b >= ELhashsize)
|
|
return null;
|
|
|
|
he = ELhash[b];
|
|
if (he == null || !he.deleted )
|
|
return he;
|
|
|
|
/* Hash table points to deleted half edge. Patch as necessary. */
|
|
ELhash[b] = null;
|
|
return null;
|
|
}
|
|
|
|
private Halfedge ELleftbnd( Point p )
|
|
{
|
|
int bucket;
|
|
Halfedge he;
|
|
|
|
/* Use hash table to get close to desired halfedge */
|
|
// use the hash function to find the place in the hash map that this
|
|
// HalfEdge should be
|
|
bucket = (int) ((p.x - xmin) / deltax * ELhashsize);
|
|
|
|
// make sure that the bucket position is within the range of the hash
|
|
// array
|
|
if ( bucket < 0 ) bucket = 0;
|
|
if ( bucket >= ELhashsize ) bucket = ELhashsize - 1;
|
|
|
|
he = ELgethash ( bucket );
|
|
|
|
// if the HE isn't found, search backwards and forwards in the hash map
|
|
// for the first non-null entry
|
|
if ( he == null )
|
|
{
|
|
for ( int i = 1; i < ELhashsize; i++ )
|
|
{
|
|
if ( (he = ELgethash ( bucket - i ) ) != null )
|
|
break;
|
|
if ( (he = ELgethash ( bucket + i ) ) != null )
|
|
break;
|
|
}
|
|
}
|
|
|
|
/* Now search linear list of halfedges for the correct one */
|
|
if ( he == ELleftend || ( he != ELrightend && right_of (he, p) ) )
|
|
{
|
|
// keep going right on the list until either the end is reached, or
|
|
// you find the 1st edge which the point isn't to the right of
|
|
do
|
|
{
|
|
he = he.ELright;
|
|
}
|
|
while ( he != ELrightend && right_of(he, p) );
|
|
he = he.ELleft;
|
|
}
|
|
else
|
|
// if the point is to the left of the HalfEdge, then search left for
|
|
// the HE just to the left of the point
|
|
{
|
|
do
|
|
{
|
|
he = he.ELleft;
|
|
}
|
|
while ( he != ELleftend && !right_of(he, p) );
|
|
}
|
|
|
|
/* Update hash table and reference counts */
|
|
if ( bucket > 0 && bucket < ELhashsize - 1)
|
|
{
|
|
ELhash[bucket] = he;
|
|
}
|
|
|
|
return he;
|
|
}
|
|
|
|
private void pushGraphEdge( Site leftSite, Site rightSite, Vector2 point1, Vector2 point2 )
|
|
{
|
|
GraphEdge newEdge = new GraphEdge ();
|
|
allEdges.Add ( newEdge );
|
|
newEdge.point1 = point1;
|
|
newEdge.point2 = point2;
|
|
|
|
newEdge.site1 = leftSite;
|
|
newEdge.site2 = rightSite;
|
|
}
|
|
|
|
private void clip_line( Edge e )
|
|
{
|
|
double pxmin, pxmax, pymin, pymax;
|
|
Site s1, s2;
|
|
|
|
double x1 = e.reg[0].coord.x;
|
|
double y1 = e.reg[0].coord.y;
|
|
double x2 = e.reg[1].coord.x;
|
|
double y2 = e.reg[1].coord.y;
|
|
double x = x2- x1;
|
|
double y = y2 - y1;
|
|
|
|
// if the distance between the two points this line was created from is
|
|
// less than the square root of 2 عن جد؟, then ignore it
|
|
if ( Math.Sqrt ( (x*x) + (y*y) ) < minDistanceBetweenSites )
|
|
{
|
|
return;
|
|
}
|
|
pxmin = borderMinX;
|
|
pymin = borderMinY;
|
|
pxmax = borderMaxX;
|
|
pymax = borderMaxY;
|
|
|
|
if ( e.a == 1.0 && e.b >= 0.0 )
|
|
{
|
|
s1 = e.ep[1];
|
|
s2 = e.ep[0];
|
|
}
|
|
else
|
|
{
|
|
s1 = e.ep[0];
|
|
s2 = e.ep[1];
|
|
}
|
|
|
|
if ( e.a == 1.0 )
|
|
{
|
|
y1 = pymin;
|
|
|
|
if ( s1 != null && s1.coord.y > pymin )
|
|
y1 = s1.coord.y;
|
|
if ( y1 > pymax )
|
|
y1 = pymax;
|
|
x1 = e.c - e.b * y1;
|
|
y2 = pymax;
|
|
|
|
if ( s2 != null && s2.coord.y < pymax )
|
|
y2 = s2.coord.y;
|
|
if ( y2 < pymin )
|
|
y2 = pymin;
|
|
x2 = e.c - e.b * y2;
|
|
if ( ( (x1 > pxmax) & (x2 > pxmax) ) | ( (x1 < pxmin) & (x2 < pxmin) ) )
|
|
return;
|
|
|
|
if ( x1 > pxmax )
|
|
{
|
|
x1 = pxmax;
|
|
y1 = ( e.c - x1 ) / e.b;
|
|
}
|
|
if ( x1 < pxmin )
|
|
{
|
|
x1 = pxmin;
|
|
y1 = ( e.c - x1 ) / e.b;
|
|
}
|
|
if ( x2 > pxmax )
|
|
{
|
|
x2 = pxmax;
|
|
y2 = ( e.c - x2 ) / e.b;
|
|
}
|
|
if ( x2 < pxmin )
|
|
{
|
|
x2 = pxmin;
|
|
y2 = ( e.c - x2 ) / e.b;
|
|
}
|
|
|
|
}
|
|
else
|
|
{
|
|
x1 = pxmin;
|
|
if ( s1 != null && s1.coord.x > pxmin )
|
|
x1 = s1.coord.x;
|
|
if ( x1 > pxmax )
|
|
x1 = pxmax;
|
|
y1 = e.c - e.a * x1;
|
|
|
|
x2 = pxmax;
|
|
if ( s2 != null && s2.coord.x < pxmax )
|
|
x2 = s2.coord.x;
|
|
if ( x2 < pxmin )
|
|
x2 = pxmin;
|
|
y2 = e.c - e.a * x2;
|
|
|
|
if (((y1 > pymax) & (y2 > pymax)) | ((y1 < pymin) & (y2 < pymin)))
|
|
return;
|
|
|
|
if ( y1 > pymax )
|
|
{
|
|
y1 = pymax;
|
|
x1 = ( e.c - y1 ) / e.a;
|
|
}
|
|
if ( y1 < pymin )
|
|
{
|
|
y1 = pymin;
|
|
x1 = ( e.c - y1 ) / e.a;
|
|
}
|
|
if ( y2 > pymax )
|
|
{
|
|
y2 = pymax;
|
|
x2 = ( e.c - y2 ) / e.a;
|
|
}
|
|
if ( y2 < pymin )
|
|
{
|
|
y2 = pymin;
|
|
x2 = ( e.c - y2 ) / e.a;
|
|
}
|
|
}
|
|
|
|
pushGraphEdge(e.reg[0], e.reg[1], new Vector2((float)x1, (float)y1), new Vector2((float)x2, (float)y2));
|
|
}
|
|
|
|
private void endpoint( Edge e, int lr, Site s )
|
|
{
|
|
e.ep[lr] = s;
|
|
if ( e.ep[RE - lr] == null )
|
|
return;
|
|
clip_line ( e );
|
|
}
|
|
|
|
/* returns true if p is to right of halfedge e */
|
|
private bool right_of(Halfedge el, Point p)
|
|
{
|
|
Edge e;
|
|
Site topsite;
|
|
bool right_of_site;
|
|
bool above, fast;
|
|
double dxp, dyp, dxs, t1, t2, t3, yl;
|
|
|
|
e = el.ELedge;
|
|
topsite = e.reg[1];
|
|
|
|
if ( p.x > topsite.coord.x )
|
|
right_of_site = true;
|
|
else
|
|
right_of_site = false;
|
|
|
|
if ( right_of_site && el.ELpm == LE )
|
|
return true;
|
|
if (!right_of_site && el.ELpm == RE )
|
|
return false;
|
|
|
|
if ( e.a == 1.0 )
|
|
{
|
|
dxp = p.x - topsite.coord.x;
|
|
dyp = p.y - topsite.coord.y;
|
|
fast = false;
|
|
|
|
if ( (!right_of_site & (e.b < 0.0)) | (right_of_site & (e.b >= 0.0)) )
|
|
{
|
|
above = dyp >= e.b * dxp;
|
|
fast = above;
|
|
}
|
|
else
|
|
{
|
|
above = p.x + p.y * e.b > e.c;
|
|
if ( e.b < 0.0 )
|
|
above = !above;
|
|
if ( !above )
|
|
fast = true;
|
|
}
|
|
if ( !fast )
|
|
{
|
|
dxs = topsite.coord.x - ( e.reg[0] ).coord.x;
|
|
above = e.b * (dxp * dxp - dyp * dyp)
|
|
< dxs * dyp * (1.0 + 2.0 * dxp / dxs + e.b * e.b);
|
|
|
|
if ( e.b < 0 )
|
|
above = !above;
|
|
}
|
|
}
|
|
else // e.b == 1.0
|
|
{
|
|
yl = e.c - e.a * p.x;
|
|
t1 = p.y - yl;
|
|
t2 = p.x - topsite.coord.x;
|
|
t3 = yl - topsite.coord.y;
|
|
above = t1 * t1 > t2 * t2 + t3 * t3;
|
|
}
|
|
return ( el.ELpm == LE ? above : !above );
|
|
}
|
|
|
|
private Site rightreg(Halfedge he)
|
|
{
|
|
if (he.ELedge == (Edge) null)
|
|
// if this halfedge has no edge, return the bottom site (whatever
|
|
// that is)
|
|
{
|
|
return (bottomsite);
|
|
}
|
|
|
|
// if the ELpm field is zero, return the site 0 that this edge bisects,
|
|
// otherwise return site number 1
|
|
return (he.ELpm == LE ? he.ELedge.reg[RE] : he.ELedge.reg[LE]);
|
|
}
|
|
|
|
private double dist( Site s, Site t )
|
|
{
|
|
double dx, dy;
|
|
dx = s.coord.x - t.coord.x;
|
|
dy = s.coord.y - t.coord.y;
|
|
return Math.Sqrt ( dx * dx + dy * dy );
|
|
}
|
|
|
|
// create a new site where the HalfEdges el1 and el2 intersect - note that
|
|
// the Point in the argument list is not used, don't know why it's there
|
|
private Site intersect( Halfedge el1, Halfedge el2 )
|
|
{
|
|
Edge e1, e2, e;
|
|
Halfedge el;
|
|
double d, xint, yint;
|
|
bool right_of_site;
|
|
Site v; // vertex
|
|
|
|
e1 = el1.ELedge;
|
|
e2 = el2.ELedge;
|
|
|
|
if ( e1 == null || e2 == null )
|
|
return null;
|
|
|
|
// if the two edges bisect the same parent, return null
|
|
if ( e1.reg[1] == e2.reg[1] )
|
|
return null;
|
|
|
|
d = e1.a * e2.b - e1.b * e2.a;
|
|
if ( -1.0e-10 < d && d < 1.0e-10 )
|
|
return null;
|
|
|
|
xint = ( e1.c * e2.b - e2.c * e1.b ) / d;
|
|
yint = ( e2.c * e1.a - e1.c * e2.a ) / d;
|
|
|
|
if ( (e1.reg[1].coord.y < e2.reg[1].coord.y)
|
|
|| (e1.reg[1].coord.y == e2.reg[1].coord.y && e1.reg[1].coord.x < e2.reg[1].coord.x) )
|
|
{
|
|
el = el1;
|
|
e = e1;
|
|
}
|
|
else
|
|
{
|
|
el = el2;
|
|
e = e2;
|
|
}
|
|
|
|
right_of_site = xint >= e.reg[1].coord.x;
|
|
if ((right_of_site && el.ELpm == LE)
|
|
|| (!right_of_site && el.ELpm == RE))
|
|
return null;
|
|
|
|
// create a new site at the point of intersection - this is a new vector
|
|
// event waiting to happen
|
|
v = new Site();
|
|
v.coord.x = xint;
|
|
v.coord.y = yint;
|
|
return v;
|
|
}
|
|
|
|
/*
|
|
* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax, deltax,
|
|
* deltay (can all be estimates). Performance suffers if they are wrong;
|
|
* better to make nsites, deltax, and deltay too big than too small. (?)
|
|
*/
|
|
private bool voronoi_bd()
|
|
{
|
|
Site newsite, bot, top, temp, p;
|
|
Site v;
|
|
Point newintstar = null;
|
|
int pm;
|
|
Halfedge lbnd, rbnd, llbnd, rrbnd, bisector;
|
|
Edge e;
|
|
|
|
PQinitialize();
|
|
ELinitialize();
|
|
|
|
bottomsite = nextone();
|
|
newsite = nextone();
|
|
while (true)
|
|
{
|
|
if (!PQempty())
|
|
{
|
|
newintstar = PQ_min();
|
|
}
|
|
// if the lowest site has a smaller y value than the lowest vector
|
|
// intersection,
|
|
// process the site otherwise process the vector intersection
|
|
|
|
if (newsite != null && (PQempty()
|
|
|| newsite.coord.y < newintstar.y
|
|
|| (newsite.coord.y == newintstar.y
|
|
&& newsite.coord.x < newintstar.x)))
|
|
{
|
|
/* new site is smallest -this is a site event */
|
|
// get the first HalfEdge to the LEFT of the new site
|
|
lbnd = ELleftbnd((newsite.coord));
|
|
// get the first HalfEdge to the RIGHT of the new site
|
|
rbnd = ELright(lbnd);
|
|
// if this halfedge has no edge,bot =bottom site (whatever that
|
|
// is)
|
|
bot = rightreg(lbnd);
|
|
// create a new edge that bisects
|
|
e = bisect(bot, newsite);
|
|
|
|
// create a new HalfEdge, setting its ELpm field to 0
|
|
bisector = HEcreate(e, LE);
|
|
// insert this new bisector edge between the left and right
|
|
// vectors in a linked list
|
|
ELinsert(lbnd, bisector);
|
|
|
|
// if the new bisector intersects with the left edge,
|
|
// remove the left edge's vertex, and put in the new one
|
|
if ((p = intersect(lbnd, bisector)) != null)
|
|
{
|
|
PQdelete(lbnd);
|
|
PQinsert(lbnd, p, dist(p, newsite));
|
|
}
|
|
lbnd = bisector;
|
|
// create a new HalfEdge, setting its ELpm field to 1
|
|
bisector = HEcreate(e, RE);
|
|
// insert the new HE to the right of the original bisector
|
|
// earlier in the IF stmt
|
|
ELinsert(lbnd, bisector);
|
|
|
|
// if this new bisector intersects with the new HalfEdge
|
|
if ((p = intersect(bisector, rbnd)) != null)
|
|
{
|
|
// push the HE into the ordered linked list of vertices
|
|
PQinsert(bisector, p, dist(p, newsite));
|
|
}
|
|
newsite = nextone();
|
|
} else if (!PQempty())
|
|
/* intersection is smallest - this is a vector event */
|
|
{
|
|
// pop the HalfEdge with the lowest vector off the ordered list
|
|
// of vectors
|
|
lbnd = PQextractmin();
|
|
// get the HalfEdge to the left of the above HE
|
|
llbnd = ELleft(lbnd);
|
|
// get the HalfEdge to the right of the above HE
|
|
rbnd = ELright(lbnd);
|
|
// get the HalfEdge to the right of the HE to the right of the
|
|
// lowest HE
|
|
rrbnd = ELright(rbnd);
|
|
// get the Site to the left of the left HE which it bisects
|
|
bot = leftreg(lbnd);
|
|
// get the Site to the right of the right HE which it bisects
|
|
top = rightreg(rbnd);
|
|
|
|
v = lbnd.vertex; // get the vertex that caused this event
|
|
makevertex(v); // set the vertex number - couldn't do this
|
|
// earlier since we didn't know when it would be processed
|
|
endpoint(lbnd.ELedge, lbnd.ELpm, v);
|
|
// set the endpoint of
|
|
// the left HalfEdge to be this vector
|
|
endpoint(rbnd.ELedge, rbnd.ELpm, v);
|
|
// set the endpoint of the right HalfEdge to
|
|
// be this vector
|
|
ELdelete(lbnd); // mark the lowest HE for
|
|
// deletion - can't delete yet because there might be pointers
|
|
// to it in Hash Map
|
|
PQdelete(rbnd);
|
|
// remove all vertex events to do with the right HE
|
|
ELdelete(rbnd); // mark the right HE for
|
|
// deletion - can't delete yet because there might be pointers
|
|
// to it in Hash Map
|
|
pm = LE; // set the pm variable to zero
|
|
|
|
if (bot.coord.y > top.coord.y)
|
|
// if the site to the left of the event is higher than the
|
|
// Site
|
|
{ // to the right of it, then swap them and set the 'pm'
|
|
// variable to 1
|
|
temp = bot;
|
|
bot = top;
|
|
top = temp;
|
|
pm = RE;
|
|
}
|
|
e = bisect(bot, top); // create an Edge (or line)
|
|
// that is between the two Sites. This creates the formula of
|
|
// the line, and assigns a line number to it
|
|
bisector = HEcreate(e, pm); // create a HE from the Edge 'e',
|
|
// and make it point to that edge
|
|
// with its ELedge field
|
|
ELinsert(llbnd, bisector); // insert the new bisector to the
|
|
// right of the left HE
|
|
endpoint(e, RE - pm, v); // set one endpoint to the new edge
|
|
// to be the vector point 'v'.
|
|
// If the site to the left of this bisector is higher than the
|
|
// right Site, then this endpoint
|
|
// is put in position 0; otherwise in pos 1
|
|
|
|
// if left HE and the new bisector intersect, then delete
|
|
// the left HE, and reinsert it
|
|
if ((p = intersect(llbnd, bisector)) != null)
|
|
{
|
|
PQdelete(llbnd);
|
|
PQinsert(llbnd, p, dist(p, bot));
|
|
}
|
|
|
|
// if right HE and the new bisector intersect, then
|
|
// reinsert it
|
|
if ((p = intersect(bisector, rrbnd)) != null)
|
|
{
|
|
PQinsert(bisector, p, dist(p, bot));
|
|
}
|
|
} else
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd))
|
|
{
|
|
e = lbnd.ELedge;
|
|
clip_line(e);
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
} // Voronoi Class End
|
|
} // namespace Voronoi2 End |