src/WeatherCM/FGVoronoi.cpp

   1 /*****************************************************************************
   2
   3  Module:       FGVoronoi.cpp
   4  Author:       Christian Mayer
   5  Date started: 28.05.99
   6
   7  -------- Copyright (C) 1999 Christian Mayer (fgfs@christianmayer.de) --------
   8
   9  This program is free software; you can redistribute it and/or modify it under
  10  the terms of the GNU General Public License as published by the Free Software
  11  Foundation; either version 2 of the License, or (at your option) any later
  12  version.
  13
  14  This program is distributed in the hope that it will be useful, but WITHOUT
  15  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  16  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  17  details.
  18
  19  You should have received a copy of the GNU General Public License along with
  20  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  21  Place - Suite 330, Boston, MA  02111-1307, USA.
  22
  23  Further information about the GNU General Public License can also be found on
  24  the world wide web at http://www.gnu.org.
  25
  26 FUNCTIONAL DESCRIPTION
  27 ------------------------------------------------------------------------------
  28 library for Voronoi Diagram calculation based on Steven Fortune 'Sweep2'
  29 FGVoronoi is the wraper to feed the voronoi calulation with a vetor of points
  30 and any class you want, as it uses templates
  31 NOTE: Sweep2 didn't free *any* memory. So I'm doing what I can, but that's not
  32       good enough...
  33
  34 HISTORY
  35 ------------------------------------------------------------------------------
  36 30.05.99   Christian Mayer      Created
  37 16.06.99   Durk Talsma          Portability for Linux
  38 11.10.1999 Christian Mayer      changed set<> to map<> on Bernie Bright's
  39                                 suggestion
  40 19.10.1999 Christian Mayer      change to use PLIB's sg instead of Point[2/3]D
  41                                 and lots of wee code cleaning
  42 *****************************************************************************/
  43
  44 /****************************************************************************/
  45 /* INCLUDES                                                                 */
  46 /****************************************************************************/
  47 #include "FGVoronoi.h"
  48 #include <Voronoi/voronoi.h>
  49 #include <Voronoi/my_memory.h>
  50
  51 #include <Voronoi/defs.h>
  52
  53 extern "C" {
  54     //forward definitions
  55     void voronoi(int triangulate, struct Site *(*nextsite)());
  56     void geominit();
  57     void freeinit(struct Freelist *fl, int size);
  58     struct Site *nextone();
  59     bool readsites(PointList input);
  60 };
  61
  62 /****************************************************************************/
  63 /********************************** CODE ************************************/
  64 /****************************************************************************/
  65 FGVoronoiOutputList Voronoiate(const FGVoronoiInputList& input)
  66 {
  67     FGVoronoiOutputList ret_list;
  68
  69     PointList p2ds;
  70
  71     FGVoronoiInputList::const_iterator it1;
  72
  73     //get the points
  74     for (it1 = input.begin(); it1 != input.end(); it1++)
  75     {
  76         p2ds.push_back(VoronoiPoint(it1->position[0], it1->position[1]));
  77     }
  78
  79     cl.clear(); //make sure that it's empty
  80
  81     if (readsites(p2ds) == false)
  82         return ret_list;
  83
  84     freeinit(&sfl, sizeof *sites);
  85     siteidx = 0;
  86     geominit();
  87     voronoi(triangulate, nextone);
  88
  89     _my_free();
  90     /*free(sites);    //to prevent a *big* memory leak...
  91     free(ELhash);   //Sweep2 didn't free *any* memory...
  92     free(PQhash);*/
  93
  94     for (cellList::iterator it2 = cl.begin(); it2 != cl.end(); it2++)
  95     {
  96         it2->sort();
  97         it2->strip();
  98
  99         //uncomment for debugging
 100         //cout << *it2;
 101
 102         Point2DList boundary;
 103
 104         //copy points
 105         PointList::iterator it3 = it2->boundary.begin();
 106
 107         if (it3->infinity == false)
 108             boundary.push_back( sgVec2Wrap(it3->p) );
 109         else
 110         {
 111             sgVec2 direction_vector;
 112             sgCopyVec2(direction_vector, it3->p);
 113
 114             it3++;
 115             sgAddVec2(direction_vector, it3->p);
 116
 117             boundary.push_back( direction_vector );
 118         }
 119
 120         for (; it3 != it2->boundary.end(); it3++)
 121         {
 122             boundary.push_back( sgVec2Wrap(it3->p) );
 123
 124         }
 125
 126         it3--;
 127         if (it3->infinity == true)
 128         {
 129             sgVec2 direction_vector;
 130             sgCopyVec2(direction_vector, it3->p);
 131
 132             it3--;
 133             sgAddVec2(direction_vector, it3->p);
 134
 135             boundary.pop_back();
 136             boundary.push_back(direction_vector);
 137         }
 138
 139         ret_list.push_back(FGVoronoiOutput(boundary, input[it2->ID].value));
 140
 141     }
 142
 143     return ret_list;
 144
 145 }
 146
 147 extern "C"
 148 {
 149
 150 /* return a single in-storage site */
 151 struct Site *nextone(void)
 152 {
 153     struct Site *s;
 154     if(siteidx < nsites)
 155     {
 156         s = &sites[siteidx];
 157         siteidx ++;
 158         return(s);
 159     }
 160     else
 161         return( (struct Site *)NULL);
 162 }
 163
 164
 165 /* sort sites on y, then x, coord */
 166 //int scomp(const struct Point *s1, const struct Point *s2)
 167 int scomp(const void *s1, const void *s2)
 168 {
 169     if(((Point*)s1) -> y < ((Point*)s2) -> y) return(-1);
 170     if(((Point*)s1) -> y > ((Point*)s2) -> y) return(1);
 171     if(((Point*)s1) -> x < ((Point*)s2) -> x) return(-1);
 172     if(((Point*)s1) -> x > ((Point*)s2) -> x) return(1);
 173     return(0);
 174 }
 175
 176
 177
 178 /* read all sites, sort, and compute xmin, xmax, ymin, ymax */
 179 bool readsites(PointList input)
 180 {
 181     int i;
 182
 183     if (input.size() == 0)
 184         return false;   //empty array
 185
 186     PointList::iterator It = input.begin();
 187
 188     nsites=0;
 189     sites = (struct Site *) myalloc(4000*sizeof(*sites));
 190
 191     while(It != input.end())
 192     {
 193         sites[nsites].coord.x = It->p[0];
 194         sites[nsites].coord.y = It->p[1];
 195
 196         sites[nsites].sitenbr = nsites;
 197         sites[nsites].refcnt = 0;
 198         nsites ++;
 199         It++;
 200
 201         if (nsites % 4000 == 0)
 202             sites = (struct Site *) my_realloc(sites,(nsites+4000)*sizeof(*sites));
 203     };
 204
 205     qsort(sites, nsites, sizeof (*sites), scomp);
 206     xmin=sites[0].coord.x;
 207     xmax=sites[0].coord.x;
 208
 209     for(i=1; i<nsites; i+=1)
 210     {
 211         if(sites[i].coord.x < xmin)
 212             xmin = sites[i].coord.x;
 213         if(sites[i].coord.x > xmax)
 214             xmax = sites[i].coord.x;
 215     };
 216
 217     ymin = sites[0].coord.y;
 218     ymax = sites[nsites-1].coord.y;
 219
 220     return true;
 221 }
 222
 223
 224 }
 225
 226