Spatial: диаграмма Вороного (Java)

Всем привет!
Сегодня хочу описать решение по генерации диаграммы Вороного с помощью java.

Диаграмма Вороного конечного множества точек S на плоскости представляет такое разбиение плоскости, при котором каждая область этого разбиения образует множество точек, более близких к одному из элементов множества S, чем к любому другому элементу множества.


У диаграммы очень много областей применения, особенно при решении spatial-задач.
/*Например, я с помощью диаграммы Вороного решал задачу кластеризации на карте Санкт-Петербурга */
К сожалению, в Oracle Spatial нет эффективного механизма (я ошибаюсь?) получения данной диаграммы. Да и к тому же Spatial - платная опция...
Поэтому, возможно, приведенное далее решение будет вам интересно!


В решение используется алгоритм от Steven Fortune - "метод заметающей прямой".


Само решение большей частью взято отсюда - http://stackoverflow.com/questions/2346148/fastest-way-to-get-the-set-of-convex-polygons-formed-by-voronoi-line-segments

И доработано таким образом, что при входных данных:
1) набор гео-точек (например, дома города);
2) ограничивающий точки полигон (например, административная граница города)

- возвращается набор полигонов, представляющих собой ячейки диаграммы Вороного.
Причем ячейки Вороного "обрезаны" ограничивающим полигоном.

В качестве примера рассмотрим формирования ячеек Вороного для Василеостровского района Санкт-Петербурга.

1) Скачаем свежие shape-файлы по Санкт-Петербургу

2) Скачиваем - если вы еще не обзавелись - Oracle Mapbuilder

3) Запускаем Mapbuilder
(указываем явно кодировку, чтобы не "потерять" кириллические символы)

java -Dfile.encoding=UTF8 -XX:MaxPermSize=512m -Xmx1024m -jar mapbuilder.jar

4) Создаем новую схему в БД для spatial-тестов

CREATE USER GEO IDENTIFIED BY GEO DEFAULT TABLESPACE USERS QUOTA UNLIMITED ON USERS;

GRANT CONNECT TO GEO;
GRANT RESOURCE TO GEO;

5) Импортируем 2 shape-файла в БД (аналогично этому) - boundary-polygon.shp и poi-point.shp в таблицы D_ZONE_POLY_8307 и D_POI_POINT_8307 соответственно.

6) Меняем у загруженных гео-данных в сферическую систему координат (возможно, вы в дальнейшем будете выводить эти данные в Mapviewer/OBIEE11g)

create table D_ZONE_POLY_3785 as select * from D_ZONE_POLY_8307;
update D_ZONE_POLY_3785 set geometry = sdo_cs.transform(geometry,'USE_SPHERICAL', 3785);
commit;

delete from user_sdo_geom_metadata where table_name = 'D_ZONE_POLY_3785';
insert into user_sdo_geom_metadata
values('D_ZONE_POLY_3785', 'GEOMETRY', sdo_dim_array(sdo_dim_element('X', -20037508.3427, 20037508.3427, 0.05), sdo_dim_element('Y', -20037508.3427, 20037508.3427, 0.05)), 3785);

create index D_ZONE_POLY_3785_sidx on D_ZONE_POLY_3785(GEOMETRY) indextype is mdsys.spatial_index;

create table D_POI_POINT_3785 as select * from D_POI_POINT_8307;
update D_POI_POINT_3785 set geometry = sdo_cs.transform(geometry,'USE_SPHERICAL', 3785);
commit;

delete from user_sdo_geom_metadata where table_name = 'D_POI_POINT_3785';
insert into user_sdo_geom_metadata
values('D_POI_POINT_3785', 'GEOMETRY', sdo_dim_array(sdo_dim_element('X', -20037508.3427, 20037508.3427, 0.05), sdo_dim_element('Y', -20037508.3427, 20037508.3427, 0.05)), 3785);

create index D_POI_POINT_3785_sidx on D_POI_POINT_3785(GEOMETRY) indextype is mdsys.spatial_index;

7) Вычленяем из загруженных данных часть для проведения теста.
Тестовая выборка будет содержать полигон административной границы Василеостровского района Санкт-Петербурга и его дома.

create table D_ZONE_POLY_TEST as
select * from D_ZONE_POLY_3785 t
where name = 'Василеостровский район';

delete from user_sdo_geom_metadata where table_name = 'D_ZONE_POLY_TEST';
insert into user_sdo_geom_metadata
values('D_ZONE_POLY_TEST', 'GEOMETRY', sdo_dim_array(sdo_dim_element('X', -20037508.3427, 20037508.3427, 0.05), sdo_dim_element('Y', -20037508.3427, 20037508.3427, 0.05)), 3785);

create index D_ZONE_POLY_TEST_sidx on D_ZONE_POLY_TEST(GEOMETRY) indextype is mdsys.spatial_index;

В таблицу D_POI_POINT_TEST сохраняем лишь дома, принадлежащие выбранному району города.

create table D_POI_POINT_TEST as
select p.*
  from D_POI_POINT_3785 p, D_ZONE_POLY_TEST z
 where SDO_CONTAINS(z.geometry, p.geometry) = 'TRUE';
 
delete from user_sdo_geom_metadata where table_name = 'D_POI_POINT_TEST';
insert into user_sdo_geom_metadata
values('D_POI_POINT_TEST', 'GEOMETRY', sdo_dim_array(sdo_dim_element('X', -20037508.3427, 20037508.3427, 0.05), sdo_dim_element('Y', -20037508.3427, 20037508.3427, 0.05)), 3785);

create index D_POI_POINT_TEST_sidx on D_POI_POINT_TEST(GEOMETRY) indextype is mdsys.spatial_index; 

8) Создадим таблицу для хранения результатов теста

create table D_VORONOI_POLY(GEOMETRY MDSYS.SDO_GEOMETRY, OSM_ID NUMBER, NEIGHBORS VARCHAR2(4000));

delete from user_sdo_geom_metadata where table_name = 'D_VORONOI_POLY';
insert into user_sdo_geom_metadata
values('D_VORONOI_POLY', 'GEOMETRY', sdo_dim_array(sdo_dim_element('X', -20037508.3427, 20037508.3427, 0.05), sdo_dim_element('Y', -20037508.3427, 20037508.3427, 0.05)), 3785);

create index D_VORONOI_POLY_sidx on D_VORONOI_POLY(GEOMETRY) indextype is mdsys.spatial_index; 

9) Скачаем архив с java-проектом генерации диаграммы Вороного.
В главном файле проекта - src\ru\servplus\geom\Voronoi.java - следует изменить настройки JDBC-доступа к тестовой схеме БД.

package ru.servplus.geom;

import java.sql.Connection;
import java.sql.DriverManager;
import java.sql.PreparedStatement;
import java.sql.ResultSet;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.ListIterator;
import java.util.PriorityQueue;
import java.util.TreeSet;

import oracle.jdbc.*;
import oracle.spatial.geometry.JGeometry;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.geom.LinearRing;



public class Voronoi {
 boolean debug = false;
 
 Geometry zonePoly;
 GeometryFactory geometryFactory;
 
 double xmin;
 double ymin;
 double xmax;
 double ymax;
 double width;
 double height;
 
 Point topleft;
    Point topright;
    Point bottomleft;
    Point bottomright;  
 
 Connection conn = null;
 
    // The set of points that control the centers of the cells
    private LinkedList<Point> pts;
    // A list of line segments that defines where the cells are divided
    private LinkedList<Edge> output;
    // Special list of points that are redundant and need to be removed in post processing
    private LinkedList<Point> mechs;
    // The sites that have not yet been processed, in acending order of X coordinate
    private PriorityQueue<Point> sites;
    // Possible upcoming cirlce events in acending order of X coordinate
    private PriorityQueue<CircleEvent> events;
    // The root of the binary search tree of the parabolic wave front
    private Arc root;
       
    
 public static void main(String[] args) {  
  System.out.println("start!");
  try {
   Voronoi vor = new Voronoi();
   vor.execute();
  }
     catch(Exception e)
     {
      e.printStackTrace();
     }

  System.out.println("end!");        
 }    
    
    
 public Voronoi() throws Exception{
        
  //Init the lists
        pts = new LinkedList<Point>();
        output = new LinkedList<Edge>();
        mechs = new LinkedList<Point>();
        sites = new PriorityQueue<Point>(100);
        events = new PriorityQueue<CircleEvent>(100);
        
  JGeometry geom = null;
  PreparedStatement pstmt = null;
  ResultSet rs = null;
  oracle.sql.STRUCT st = null; 
  Point p;
  int key;
  
     String jdbcUrl = "jdbc:oracle:thin:@localhost:1521:ORCL";
     String jdbcUsr = "GEO";
     String jdbcPwd = "GEO";
      
     Class.forName ("oracle.jdbc.OracleDriver").newInstance();
     this.conn = DriverManager.getConnection(jdbcUrl, jdbcUsr, jdbcPwd);
  
     if(conn == null)
     {
       System.out.print("cannot get connection!");
     }
     
     //calculate zone MBR
     pstmt = conn.prepareStatement(
       "select max(decode(t.id, 1, t.x, null)) x_min, "+
          "       max(decode(t.id, 2, t.x, null)) x_max, "+
          "       max(decode(t.id, 1, t.y, null)) y_min, "+
          "       max(decode(t.id, 2, t.y, null)) y_max  "+
          "    from TABLE (select sdo_util.getvertices(SDO_GEOM.SDO_MBR(geometry)) from D_ZONE_POLY_TEST) t");
     
     rs = pstmt.executeQuery(); 
     while (rs.next()){
      xmin = rs.getDouble("X_MIN") - 1000;
      ymin = rs.getDouble("Y_MIN") - 1000;
      xmax = rs.getDouble("X_MAX") + 1000;
      ymax = rs.getDouble("Y_MAX") + 1000;
     }          
     height = ymax - ymin;
  width = xmax - xmin;     
  pstmt.close();  
  
     topleft = new Point(xmin,ymin);
        topright = new Point(xmax,ymin);
        bottomleft = new Point(xmin,ymax);
        bottomright = new Point(xmax,ymax);  
  
  //get zone Points
     pstmt = conn.prepareStatement(
       "select OSM_ID, GEOMETRY " +
       "from D_POI_POINT_TEST " +
       "where 1=1 ");
     rs = pstmt.executeQuery(); 
          
     while (rs.next())
     {
      st = (oracle.sql.STRUCT) rs.getObject("GEOMETRY");
      key = rs.getInt("OSM_ID");
      geom = JGeometry.load(st);
      double[] coords = geom.getPoint(); 
      p = new Point(coords[0],coords[1]);
      p.key = key;            
            pts.add(p);  
     }         
     pstmt.close();
        
     ListIterator<Point> points = pts.listIterator(0);
    
  //get zone polygon
  geometryFactory = new GeometryFactory();
  
     pstmt = conn.prepareStatement("select GEOMETRY from D_ZONE_POLY_TEST ");
     rs = pstmt.executeQuery(); 

     zonePoly = null;
     LinearRing zoneLr;
     while (rs.next())
     {
      st = (oracle.sql.STRUCT) rs.getObject("GEOMETRY");
      geom = JGeometry.load(st);
      Coordinate[] coords = new Coordinate[geom.getNumPoints()];

      for (int i=0; i<geom.getNumPoints(); i++){
       coords[i] = new Coordinate(geom.getOrdinatesArray()[i*2],geom.getOrdinatesArray()[i*2 + 1]);
      } 
      
      zoneLr = geometryFactory.createLinearRing(coords);
      zonePoly = geometryFactory.createPolygon(zoneLr, null);
     }   
     pstmt.close();
 }

 
 public void execute() throws Exception{

  JGeometry geom = null;
  PreparedStatement pstmt = null;
  //ResultSet rs = null;
  oracle.sql.STRUCT st = null;  
  Point p;
  
  if(conn == null)
     {
         System.out.print("cannot get connection!");
     } 

        // run fortune's algorithm to get new line segments
        runFortune(pts);            
        
        APoint v,ov,frst;
        Point n;
        StringBuffer neighborsConcat;
        String neighborsConcatStr = null;
        ListIterator<Point> siteit = pts.listIterator(0);
        
        pstmt = conn.prepareStatement("DELETE FROM D_VORONOI_POLY");
        pstmt.execute(); 
        
        pstmt = conn.prepareStatement("INSERT INTO D_VORONOI_POLY(GEOMETRY, OSM_ID, NEIGHBORS) VALUES(?,?,?)");
        
        ((OraclePreparedStatement)pstmt).setExecuteBatch (pts.size());
        
        while (siteit.hasNext()) {
            p = siteit.next();
            
            if (p.sverts != null){
             Iterator<APoint> vert = p.sverts.iterator();
             if (vert.hasNext()){
              v = vert.next();
              
              double[] ring = new double[(p.sverts.size() + 1) * 2];
              int f = 0;
              frst = v;
              
              ring[f++] = frst.x;
           ring[f++] = frst.y;
              
              while (vert.hasNext()) {
               ov = v;
               v = vert.next();
               
               ring[f++] = v.x;
               ring[f++] = v.y;
              }
              
              ring[f++] = frst.x;
           ring[f++] = frst.y;
              
              Coordinate[] coords = new Coordinate[ring.length/2];

              for (int i=0; i<ring.length/2; i++){
               coords[i] = new Coordinate(ring[i*2],ring[i*2 + 1]);
              } 
              
              LinearRing cellLr = geometryFactory.createLinearRing(coords);
              Geometry cellPoly = geometryFactory.createPolygon(cellLr, null);           
           
              Geometry cellIntersect = zonePoly.intersection(cellPoly);

              if (cellIntersect != null) {
                               
               Object[] intersectRings = new Object[cellIntersect.getNumGeometries()];
               for(int zz=0; zz<cellIntersect.getNumGeometries(); zz++) {
                
                Coordinate[] insersectCoord = cellIntersect.getGeometryN(zz).getCoordinates();
                   double[] intersectRing = new double[insersectCoord.length * 2];
                   for (int i=0; i< insersectCoord.length;i++){
                    intersectRing[i*2] = insersectCoord[i].x;
                    intersectRing[i*2+1] = insersectCoord[i].y;
                   }
                
                   intersectRings[zz] = intersectRing;
               }               
               
               try {
       geom = JGeometry.createLinearPolygon(intersectRings, 2, 3785);
              
       if ( intersectRings.length > 1 && geom.getType() == 3 ) {
        geom.setType( JGeometry.GTYPE_MULTIPOLYGON );
        //System.out.println("GTYPE_MULTIPOLYGON - "+p.key); 
       }
       
       st = JGeometry.store(geom, conn);
       
       pstmt.setObject(1, st);
       pstmt.setObject(2, p.key);
       
       
      } catch (RuntimeException e) {
       e.printStackTrace();
       System.out.println(p.key);       
      }
                     
                     neighborsConcat = new StringBuffer();                     
                     
                     Iterator<Point> neighbors = p.neighbors.iterator();

                     while (neighbors.hasNext()){
                      n = neighbors.next();
                      neighborsConcat.append(n.key+",");
                     }
                     if (neighborsConcat.length() > 0)
                      neighborsConcatStr = neighborsConcat.substring(0, neighborsConcat.length()-1);
                     pstmt.setObject(3, neighborsConcatStr);

                     pstmt.executeUpdate();                                 
              }               
             }
            }            
        }    
        
        ((OraclePreparedStatement)pstmt).sendBatch(); // JDBC sends the queued request

        conn.commit();

        pstmt.close();             
 }


    void runFortune(LinkedList<Point> pts) {

     sites.clear();
     events.clear();
     output.clear();
     mechs.clear();
     root = null;

     ListIterator<Point> i = pts.listIterator(0);
     Point meh;
     while (i.hasNext()) {
      meh = i.next();
      sites.offer(meh);
      if (meh.sverts == null)
       meh.sverts = new TreeSet<APoint>();    
      meh.sverts.clear();
      if (meh.neighbors == null)
       meh.neighbors = new TreeSet<Point>();    
      meh.neighbors.clear();
     }

     // Process the queues; select the top element with smaller x coordinate.
     while (sites.size() > 0) {
      if ((events.size() > 0) && ((events.peek().xpos) <= (sites.peek().x))) {
       processCircleEvent(events.poll());
      } else {
       //process a site event by adding a curve to the parabolic front
       frontInsert(sites.poll());
      }
     }

     // After all points are processed, do the remaining circle events.
     while (events.size() > 0) {
      processCircleEvent(events.poll());
     }

     // Clean up dangling edges.
     finishEdges();

     //make polygons from edges
     polyFinish();
    }

    private void processCircleEvent(CircleEvent event) {
        if (event.valid) {
         Arc parc = event.a;
         
            //start a new edge
            Edge edgy = new Edge(event.p, parc.prev.focus, parc.next.focus);
            
            // Remove the associated arc from the front;
            if (parc.prev != null) {
                parc.prev.next = parc.next;
                parc.prev.edge1 = edgy;
            }
            if (parc.next != null) {
                parc.next.prev = parc.prev;
                parc.next.edge0 = edgy;
            }
            
            // Finish the edges before and after this arc.
            if (parc.edge0 != null) {
             parc.edge0.finish(event.p);
            }
            if (parc.edge1 != null) {
             parc.edge1.finish(event.p);
            }

            // Recheck circle events on either side of p:
            if (parc.prev != null)
                checkCircleEvent(parc.prev, event.xpos);
            if (parc.next != null)
                checkCircleEvent(parc.next, event.xpos);

        }
    }

    void frontInsert(Point focus) {
        if (root == null) {
            root = new Arc(focus);
            return;
        }

        Arc parc = root;
        while (parc != null) {
            CircleResultPack rez = intersect(focus, parc);
            if (rez.valid) {
                // New parabola intersects parc.  If necessary, duplicate parc.
             
                if (parc.next != null) {
                    CircleResultPack rezz = intersect(focus, parc.next);
                    if (!rezz.valid){
                        Arc bla = new Arc(parc.focus);
                        bla.prev = parc;
                        bla.next = parc.next;
                        parc.next.prev = bla;
                        parc.next = bla;
                    }
                
                } 
                else {
                    parc.next = new Arc(parc.focus);
                    parc.next.prev = parc;
                }
                parc.next.edge1 = parc.edge1;

                // Add new arc between parc and parc.next.
                Arc bla = new Arc(focus);
                bla.prev = parc;
                bla.next = parc.next;
                parc.next.prev = bla;
                parc.next = bla;

                parc = parc.next; // Now parc points to the new arc.

                // Add new half-edges connected to parc's endpoints.
                parc.edge0 = new Edge(rez.center, parc.prev.focus, parc.focus);
                parc.prev.edge1 = parc.edge0;
                parc.edge1 = new Edge(rez.center, parc.focus, parc.next.focus);
                parc.next.edge0 = parc.edge1;
                
                //save information for removing colinear edges
                rez.center.mech = true;
                rez.center.mech0 = parc.edge0;
                rez.center.mech1 = parc.edge1;
                mechs.add(rez.center);

                // Check for new circle events around the new arc:
                checkCircleEvent(parc, focus.x);
                checkCircleEvent(parc.prev, focus.x);
                checkCircleEvent(parc.next, focus.x);

                return;
            }

            //proceed to next arc
            parc = parc.next;
        }

        // Special case: If p never intersects an arc, append it to the list.        
        parc = root;
        while (parc.next != null)
            parc = parc.next; // Find the last node.
        parc.next = new Arc(focus);
        parc.next.prev = parc;
        Point start = new Point(0, (parc.next.focus.y + parc.focus.y) / 2);
        parc.next.edge0 = new Edge(start, parc.focus, parc.next.focus);
        parc.edge1 = parc.next.edge0;
        
        //focus.verts.addFirst(start);

    }

    void checkCircleEvent(Arc parc, double xpos) {
        // Invalidate any old event.
        if ((parc.event != null) && (parc.event.xpos != xpos))
            parc.event.valid = false;
        parc.event = null;

        if ((parc.prev == null) || (parc.next == null)) return;

        CircleResultPack result = circle(parc.prev.focus, parc.focus, parc.next.focus);
        if (result.valid && result.rightmostX > xpos) {
            // Create new event.
            parc.event = new CircleEvent(result.rightmostX, result.center, parc);
            events.offer(parc.event);
        }

    }

    // Find the rightmost point on the circle through a,b,c.
    CircleResultPack circle(Point a, Point b, Point c) {
        CircleResultPack result = new CircleResultPack();

        // Check that bc is a "right turn" from ab.
        if ((b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y) > 0) {
            result.valid = false;
            return result;
        }

        // Algorithm from O'Rourke 2ed p. 189.
        double A = b.x - a.x;
        double B = b.y - a.y;
        double C = c.x - a.x;
        double D = c.y - a.y;
        double E = A * (a.x + b.x) + B * (a.y + b.y);
        double F = C * (a.x + c.x) + D * (a.y + c.y);
        double G = 2 * (A * (c.y - b.y) - B * (c.x - b.x));

        if (G == 0) { // Points are co-linear.
            result.valid = false;
            return result;
        }

        // centerpoint of the circle.
        Point o = new Point((D * E - B * F) / G, (A * F - C * E) / G);
        result.center = o;

        // o.x plus radius equals max x coordinate.
        result.rightmostX = o.x + Math.sqrt(Math.pow(a.x - o.x, 2.0) + Math.pow(a.y - o.y, 2.0));

        result.valid = true;
        return result;
    }

    // Will a new parabola at point p intersect with arc i?
    CircleResultPack intersect(Point p, Arc i) {
        CircleResultPack res = new CircleResultPack();
        res.valid = false;
        if (i.focus.x == p.x) return res;

        double a = 0.0;
        double b = 0.0;
        if (i.prev != null) // Get the intersection of i->prev, i.
            a = intersection(i.prev.focus, i.focus, p.x).y;
        if (i.next != null) // Get the intersection of i->next, i.
            b = intersection(i.focus, i.next.focus, p.x).y;

        if ((i.prev == null || a <= p.y) && (i.next == null || p.y <= b)) {
            res.center = new Point(xmin, p.y);

            // Plug it back into the parabola equation to get the x coordinate
            res.center.x = (i.focus.x * i.focus.x + (i.focus.y - res.center.y) * (i.focus.y - res.center.y) - p.x * p.x) / (2 * i.focus.x - 2 * p.x);

            res.valid = true;
            return res;
        }
        return res;
    }

    // Where do two parabolas intersect?
    Point intersection(Point p0, Point p1, double l) {
        Point res = new Point(xmin, ymin);
        Point p = p0;

        if (p0.x == p1.x) {
            res.y = (p0.y + p1.y) / 2;
        } else if (p1.x == l) {
            res.y = p1.y;
        } else if (p0.x == l) {
            res.y = p0.y;
            p = p1;
        } else {
            // Use the quadratic formula.
            double z0 = 2 * (p0.x - l);
            double z1 = 2 * (p1.x - l);

            double a = 1 / z0 - 1 / z1;
            double b = -2 * (p0.y / z0 - p1.y / z1);
            double c = (p0.y * p0.y + p0.x * p0.x - l * l) / z0 - (p1.y * p1.y + p1.x * p1.x - l * l) / z1;

            res.y = (-b - Math.sqrt((b * b - 4 * a * c))) / (2 * a);
        }
        // Plug back into one of the parabola equations.
        res.x = (p.x * p.x + (p.y - res.y) * (p.y - res.y) - l * l) / (2 * p.x - 2 * l);
        return res;
    }

    void finishEdges() {
        // Advance the sweep line so no parabolas can cross the bounding box.
        double l = width * 2 + height;
        
        // Extend each remaining segment to the new parabola intersections.
     
        Arc i = root;
        while (i != null) {         
            if (i.edge1 != null) {
             Point fp = intersection(i.focus, i.next.focus, l * 2);
                i.edge1.finish(fp);
            }
            i = i.next;
        }
       
    }
    
    void polyFinish(){
     //clean up mechs
     //These are single lines formed by 3 colinear points
     //they are a normal result of fortune's algorithm
     //We only want 2 points
        Point unpoint;
        ListIterator<Point> mek = mechs.listIterator(0);
        while (mek.hasNext()) {
         unpoint = mek.next();
         unpoint.mech0.start = unpoint.mech1.end;
         unpoint.mech1.done = false;
        }
     
        
        //create a TreeSet<APoint> for every site
        //APoints will be ordered by thier angle from the site.
        Point pk;
        ListIterator<Point> eeek = pts.listIterator(0);
        while (eeek.hasNext()) {
         pk = eeek.next();
            pk.sortSet = new TreeSet<Point>();
            pk.sverts = new TreeSet<APoint>();
            pk.neighbors = new TreeSet<Point>();
        }
        
        //define the bounding box
        Point[] xing = new Point[4];
        
        //the corners of the bounding box
        topleft = new Point(xmin,ymin);
        topright = new Point(xmax,ymin);
        bottomleft = new Point(xmin,ymax);
        bottomright = new Point(xmax,ymax);
        
        //four sorted sets of Points that lie along each side of the bounding box
        TreeSet[] sides = new TreeSet[4];
        sides[0] = new TreeSet<Point>();//top
        sides[1] = new TreeSet<Point>();//right
        sides[2] = new TreeSet<Point>();//bottom
        sides[3] = new TreeSet<Point>();//left
        
        //add the corners to the lists of sidepoints        
        /*
        sides[0].add(topleft);
        sides[0].add(topright);
        sides[1].add(topright);
        sides[1].add(bottomright);
        sides[2].add(bottomleft);
        sides[2].add(bottomright);
        sides[3].add(topleft);
        sides[3].add(bottomleft);
        */
        
        //the edge in question
        Edge e;
        ListIterator<Edge> j;
        
        /*
         * loop over the edges
         * 
         * trim edges that intersect the bounding box,
         * and add edges between those points of intersection,
         * that close all the polygons that were opened by the trimming.
         * 
         */
        j = output.listIterator(0);
        
        while (j.hasNext()) {
            e = j.next();
            
            if (e.done){             
             //check if the edge intersects with the bounds
             xing[0] = lineIntersect( topleft, topright, e.start, e.end );
             xing[1] = lineIntersect( topright, bottomright, e.start, e.end );
             xing[2] = lineIntersect( bottomleft, bottomright, e.start, e.end );
             xing[3] = lineIntersect( topleft, bottomleft, e.start, e.end );
             
             for (int ii=0; ii<4; ii++){
              //top, right, bottom, left
              //inside function catchs double precision mistakes  
              if (xing[ii] != null && inside(xing[ii]) ){
               // e intersects side ii
               if (ii==0)
                xing[ii].y = ymin;
               else if (ii==1)
                xing[ii].x = xmax;
               else if (ii==2)
                xing[ii].y = ymax;
               else if (ii==3)
                xing[ii].x = xmin;
                
               sides[ii].add(xing[ii]);
               xing[ii].assoc = e;
               
               //shorten e
               if (outside(e.start)){
                e.start = xing[ii];               
               }
               else {
                e.end = xing[ii];                
               }
              }
             }
            }
        }
  
        /*
         * loop over side points
         * all intersections have been found, turn lists of side points into edges with proper
         * site associations
         */
        
        Iterator<Point> jj;
        Point sp,osp; //sidepoint, old sidepoint
        
        for (int ii=0; ii<4; ii++){ // for each side
         if (sides[ii].size() >= 2){ // if there was anything more than just the corner points          
          jj = sides[ii].iterator();
          sp = jj.next();

          while (jj.hasNext()){
           osp = sp;
           sp = jj.next();

           //must determine which site this edge should border
           //will be the site that both edges connected to 
           //sp and osp have in common
           //note this does not handle corner cases
           Point commonSite = null;
           if (osp.assoc != null && sp.assoc != null){
            if (osp.assoc.siteA == sp.assoc.siteA || osp.assoc.siteA == sp.assoc.siteB) {
             commonSite = osp.assoc.siteA;
            } else if (osp.assoc.siteB == sp.assoc.siteA || osp.assoc.siteB == sp.assoc.siteB) {
             commonSite = osp.assoc.siteB;
            }
           }
           
           Edge xo = new Edge(osp, commonSite, null);
           xo.finish(sp);
          }
         }
        }
        
        //associate corners with sites
        Point commonSite = null;
        Point corner = null;
        // standing at the corner, looking in, sp will be the sidepoint on your right and osp will be to your left.
        osp = null;
        sp = null;
        
        for (int ii=0; ii<4; ii++) {
         if (ii == 0){
          corner = topleft;
          sp = (Point)(sides[3].first());
          osp = (Point)(sides[0].first());  
          
         } else if (ii == 1){
          corner = topright;
          sp = (Point)(sides[0].last());
          osp = (Point)(sides[1].first());
         
         } else if (ii == 2){
          corner = bottomright;
          sp = (Point)(sides[1].last());
          osp = (Point)(sides[2].last());         
         } else if (ii == 3){
          corner = bottomleft;
          sp = (Point)(sides[2].first());
          osp = (Point)(sides[3].last());
         }
         
         if (osp != null && sp != null) {
             if ( (osp.assoc.siteA == sp.assoc.siteA || osp.assoc.siteA == sp.assoc.siteB) && osp.assoc.siteA != null ) {
              commonSite = osp.assoc.siteA;
             } else if ( (osp.assoc.siteB == sp.assoc.siteA || osp.assoc.siteB == sp.assoc.siteB) && osp.assoc.siteB != null ) {
              commonSite = osp.assoc.siteB;
             }
             
             Edge xo = new Edge( osp, commonSite, null );
             xo.finish(corner);
             Edge po = new Edge( sp, commonSite, null );             
             po.finish(corner); 
         }
         
        }
  
       
        
        /*
         * loop over edges again
         * add the endpoints of the edge to a sorted set of points,
         * the Points are sorted based on thier angle from the site
         * associated with both of the sites that are next to this edge
         * there are only ever two, but for the boundary cases, one of them will be null,
         * ignore it.
         * 
         */
        j = output.listIterator(0);
        while (j.hasNext()) {
            e = j.next();
            if (e.done){
             //add the edge's endpoints to the sorted set of verts for the sites on each side.
             if (e.siteA != null){
              if (!outside(e.start))
               e.siteA.sortSet.add( e.start );
              if (!outside(e.end))
               e.siteA.sortSet.add( e.end );
              
             }
             if (e.siteB != null){
              if (!outside(e.start))
               e.siteB.sortSet.add( e.start );
              if (!outside(e.end))
               e.siteB.sortSet.add( e.end );
             }
             if ((e.siteA != null) && (e.siteB != null)){
              
              Coordinate[] coords = new Coordinate[2];
              coords[0] = new Coordinate(e.start.x, e.start.y);
              coords[1] = new Coordinate(e.end.x, e.end.y);
              LineString ls = geometryFactory.createLineString(coords);
              
              //boolean isCover = zonePoly.covers(ls);
              boolean isIntersect = zonePoly.intersects(ls);
              //geometryFactory
              
              if (isIntersect /*isCover*/){
                  e.siteA.neighbors.add(e.siteB);
                  e.siteB.neighbors.add(e.siteA);               
              }
              
             }
             
            }
        }
        
        /*
         * loop over the sites
         * for each site, take the sorted set of Points (in acending x value)
         * and add them all to the sorted set of APoints (in acending angle from site)
         */
        Iterator<Point> uu;
        Point eko;
        ListIterator<Point> oo = pts.listIterator(0);
        Point site;
        while (oo.hasNext()) {
            site = oo.next();
            
            if (site.sortSet.size() >= 3){
             uu = site.sortSet.iterator();
             while (uu.hasNext()){
              eko = uu.next();              
              site.sverts.add( new APoint(eko, radFromPoints(site,eko)) );                             
             }
            }
        }
        
        
    }

    boolean outside(Point q){
     //is a Point outside the bounding box?
     return ( (q.x < xmin) || 
        (q.x > xmax) || 
        (q.y < ymin) || 
        (q.y > ymax) );
    }
    
    boolean inside(Point q){
     //is a Point outside the bounding box?
     return ( (q.x >= xmin - 1) && 
        (q.x <= xmax + 1) && 
        (q.y >= ymin - 1) && 
        (q.y <= ymax + 1) );
    }    

    double radFromPoints(Point A, Point B){
        double dx = B.x-A.x;
        double dy = B.y-A.y;
        double rad = Math.atan(dy/dx);
        if(dx<0) rad+=Math.PI;
        return rad;
    }

    class Point implements Comparable<Point> {

        public double x, y;
        public int key;
        
        public TreeSet<Point> sortSet;
        public TreeSet<APoint> sverts;
        public TreeSet<Point> neighbors;
        public Edge assoc;
        
        public boolean mech = false;
        public Edge mech0, mech1;
        
        public int fillcolor;

        public Point(double X, double Y) {
            x = X;
            y = Y;
        }

        public int compareTo(Point foo) {
         //return zero only if this is the same object.
         if (this == foo) {
          return 0;
         } else {
          int rx = ((Double) this.x).compareTo((Double) foo.x);
          if (rx != 0){
           return rx;
          } else {
           int ry = ((Double) this.y).compareTo((Double) foo.y);
           if (ry != 0){
            return ry;
           } else {
            return 1;
           }
          }
         }
        }
        
    }
    
    class APoint implements Comparable<APoint>{
     
     public double x, y;
     public double angle;
     
     public APoint(Point p, double ang){
      x = p.x;
      y = p.y;
      angle = ang;
     }
     
        public int compareTo(APoint foo) {
         //return zero only if this is the same object.
         if (this == foo) {
          return 0;
         } else {
          return ((Double) this.angle).compareTo((Double) foo.angle);
         }
        }
    }

    class CircleEvent implements Comparable<CircleEvent> {

        public double xpos;
        public Point p;
        public Arc a;
        public boolean valid;

        public CircleEvent(double X, Point P, Arc A) {
            xpos = X;
            a = A;
            p = P;
            valid = true;
        }

        public int compareTo(CircleEvent foo) {
            return ((Double) this.xpos).compareTo((Double) foo.xpos);
        }
    }

    class Edge {

        public Point start, end;
        public boolean done;
        public Point siteA, siteB;

        public Edge(Point p, Point sitea, Point siteb) {
            start = p;
            end = new Point(0, 0);
            done = false;
            siteA = sitea;
            siteB = siteb;
            output.add(this);
        }

        public void finish(Point p) {
            if (done) {
                return;
            }
            end = p;
            done = true;
        }
    }

    class Arc {
        //parabolic arc is the set of points eqadistant from a focus point and the beach line

        public Point focus;
        //these object exsit in a linked list
        public Arc next, prev;
        //
        public CircleEvent event;
        //
        public Edge edge0, edge1;

        public Arc(Point p) {
            focus = p;
            next = null;
            prev = null;
            event = null;
            edge0 = null;
            edge1 = null;
        }
    }
    
    class CircleResultPack {
        public boolean valid;
        public Point center;
        public double rightmostX;
    }
    
    Point midpoint(Point a, Point b){
     return new Point( (a.x+b.x)/2, (a.y+b.y)/2 );
    }
    
    Point lineIntersect(Point A, Point B, Point C, Point D){
     /*
      * if line segment AB intersects line segment CD,
      * return the point of intersection
      * else return null
      */
     double ex = B.x - A.x;
     double ey = B.y - A.y;
     double fx = D.x - C.x;
     double fy = D.y - C.y;
     double h = ( (A.x-C.x) * (-ey) + (A.y-C.y) * ex) / (fx * (-ey) + fy * ex);
     if (0 <= h && h <= 1){
      return new Point( C.x + h * fx,
            C.y + h * fy );
     } else return null;
    }     
}

10) После запуска java-программы проверим наличие данных в таблице D_VORONOI_POLY с помощью sqldeveloper (рекомендую скачать его и использовать для работы со Spatial-объектами БД Oracle)


11) Дальше вы сами вправе решать, что делать с получившимися полигонами... склеить их, или использовать в своих эвристических алгоритмах...
Но следует помнить о ряде ограничений данного решения:
- Диаграмма строится по точкам внутри замкнутого полигона, если полигон содержит анклавы - ситуацию придется как-то обрабатывать;
- Вам желательно самим предусмотреть обработку точек-дубликатов (когда несколько РАЗНЫХ домов имеют одинаковые координаты);
- Если ограничивающий полигон не является выпуклым многоугольником, то возможны "некрасивые ситуации" (следует помнить, что сама диаграмма строится внутри MBR - "минимального охватывающего точки прямоугольник" - и уже затем от каждой ячейки Вороного "отрезаются кусочки" за границами охватывающего полигона).


На этом все, надеюсь кому-то это поможет!