source: subversion/applications/utils/coastcheck/closeshp.c @ 34318

Last change on this file since 34318 was 13971, checked in by jonb, 11 years ago

coastcheck: Fix objects in the overlapping area between tiles.

File size: 48.7 KB
Line 
1/* closeshp - Takes two shape files and bounding box and returns a new shape
2 * file with closed polygon as determined by the bounding box.
3 * By Martijn van Oosterhout <kleptog@svana.org> Copyright 2008
4 * Licence: GPL
5 *
6 * This is a similar idea to close-areas.pl for osmarender, except
7 * considerably easier due to the fact that all the segments are in order
8 * and connected already. There are two shapefiles: one where coast2shp has
9 * managed to close the way, these are type POLYGON. The second is where the
10 * closing didn't work, these are type ARC.
11 *
12 * The main differences between the two are:
13 * - If a closed polygon completely surrounds (but does not touch the edge
14 *   of) the bounding box it produce something, where as an ARC does not.
15 * - If a polygon intersects the box you always get something whereas if a
16 *   arc ends in the box, it gets dropped.
17 */
18
19#define _GNU_SOURCE
20#include <stdio.h>
21#include <unistd.h>
22#include <stdlib.h>
23#include <string.h>
24#include <assert.h>
25#include <getopt.h>
26#include <errno.h>
27#include <math.h>
28
29#include <shapefil.h>
30
31#define VERBOSE 0
32#define TEST 0
33#define UNUSED __attribute__((unused))
34
35#define MERC_MAX (20037508.34f)
36/* UserParam: DIVISIONS: Indicates how many blocks to divide the world into,
37 * in each direction. The default is 400x400, which works out to be 100
38 * mercator km */
39#define DIVISIONS 400
40#define MERC_BLOCK (2*MERC_MAX/DIVISIONS)
41/* UserParam: RESOLUTION: The simplification level. Essentially, any nodes
42 * closer to each other than this many units will be simplified away. When
43 * generating files for lower zooms, setting this produces simplified shapes.
44 * A good rule of thumb is 40,000,000/(2^(max_zoom+8)). In other words, 1 is
45 * good for upto zoom 17. Zero disables this simplification.
46 * Side Effect: Non-zero causes all shapes to be loaded in memory */
47#define RESOLUTION  0
48/* UserParam: TILE_OVERLAP: Number of mercator metres the tiles overlap, so
49 * the antialising doesn't cause wierd effects. As a rule of thumb
50 * RESOLUTION*2^num_levels where num_levels is the number of zoom levels you
51 * want the shapes to be good for. */
52#define TILE_OVERLAP  150
53
54/* No further user edittable parameters */
55#define INIT_MAX_NODES (1<<19)
56int MAX_NODES = INIT_MAX_NODES;
57#define MAX_SEGS  100
58int MAX_SUBAREAS;   /* Was a define, not anymore. Auto grown array, starting at... */
59#define INIT_MAX_SUBAREAS 1024
60#define MAX_NODES_PER_ARC 10000
61
62static double *v_x, *v_y;
63static int *Parts, *Parttypes;
64
65static char used_bitmap[DIVISIONS * DIVISIONS];
66
67static SHPHandle shp_out;
68static DBFHandle dbf_out;
69struct source;
70
71/* Column of error flag */
72#define DBF_OUT_ERROR   0
73#define DBF_OUT_TILE_X   1
74#define DBF_OUT_TILE_Y   2
75struct segment
76{
77  struct source *src;   /* Shapefile of this segment */
78  int index;       /* Index of object in shapefile */
79  double sx, sy;   /* Where the segment enters the box */
80  double ex, ey;   /* Where the segment leaves the box */
81  int snode, enode; /* Numbers of nodes in polygon belonging to this segment */
82  double sa, ea;   /* Angle (in radians) of start/endpoint from centre of box */
83  int next;        /* In joining stage, have we used this one already */
84};
85
86/* This is used to track areas that fall completely within a box. After we
87 * have dealt will all the areas that intersect, we go through these to
88 * determine if they are contained. If a negative subarea is contained with
89 * a positive one, we append it to the shape. If a positive has no
90 * surrounding we just output it. And a positive within a positive is an
91 * error. */
92struct subarea
93{
94  struct source *src;   /* Shapefile of this segment */
95  int index;       /* Index of this object in shapefile */
96  double x, y;     /* representative point on shape */
97  double areasize; /* Positive area is land, negative is lake */
98  int used;        /* Have we output this subarea yet */
99};
100
101struct state
102{
103  int x,y;
104 
105  double lb[2];
106  double rt[2];
107
108//  SHPHandle shp;
109//  DBFHandle dbf;
110//  SHPTree  *shx;
111   
112  struct segment seg_list[MAX_SEGS];
113  int seg_count;
114 
115  struct subarea *sub_areas;
116  int subarea_count;
117  int subarea_nodecount;
118 
119  int enclosed;
120};
121
122struct source
123{
124  SHPHandle shp;
125  DBFHandle dbf;
126  SHPTree *shx;
127  int shape_count;
128  int in_memory;
129  SHPObject **shapes;
130};
131
132void OutputSegs( struct state *state );
133void Process( struct state *state, struct source *src, int polygon );
134static void InitBitmap( struct source *src, struct source *temp );
135static double CalcArea( const SHPObject *obj );
136static void SplitCoastlines( struct source *poly, struct source *arc, char *out_filename );
137static int contains( double x, double y, double *v_x, double *v_y, int vertices );
138static int CopyShapeToArray( struct source *src, int index, int snode, int enode, int node_count );
139static void ResizeSubareas( struct state *state, int count );
140void ProcessSubareas(struct state *state, int *pc, int *nc);
141static int SourceOpen( struct source *src, char *filename, int shapetype );
142
143static inline void mark_tile( int x, int y )
144{
145  used_bitmap[ x + y * DIVISIONS ] = 1;
146}
147static inline int check_tile2( int x, int y )
148{
149  if (x >= 0 && x < DIVISIONS && y >= 0 && y < DIVISIONS)
150    return used_bitmap[ x + y * DIVISIONS ];
151  else
152    return 0;
153}
154static inline int check_tile( int x, int y )
155{
156  /* The generated tiles overlap one another so it is possible that a
157   * shape near the edge of an adjacent tile may be relevant too
158   */
159  return 
160    check_tile2(x-1, y-1) || check_tile2(x, y-1) || check_tile2(x+1, y-1) ||
161    check_tile2(x-1, y)   || check_tile2(x, y)   || check_tile2(x+1, y)   ||
162    check_tile2(x-1, y+1) || check_tile2(x, y+1) || check_tile2(x+1, y+1);
163}
164
165static inline SHPObject *source_get_shape( struct source *src, int i )
166{
167  if( src->in_memory )
168  {
169    if( !src->shapes[i] )
170      src->shapes[i] = SHPReadObject( src->shp, i );
171    return src->shapes[i];
172  }
173  return SHPReadObject( src->shp, i );
174}
175
176static inline void source_release_shape( struct source *src, SHPObject *obj )
177{
178  if( !src->in_memory )
179    SHPDestroyObject(obj);
180}
181
182int main( int argc, char *argv[] )
183{
184  char out_filename[256];
185 
186  if( argc != 4 )
187  {
188    fprintf( stderr, "closeshp poly_shape arc_shape output_shape\n" );
189    return 1;
190  }
191  char *poly_file = argv[1];
192  char *arc_file = argv[2];
193  char *out_file = argv[3];
194
195  struct source poly, arc, simplified;
196  memset( &poly, 0, sizeof(poly) );
197  memset( &arc, 0, sizeof(arc) );
198  memset( &simplified, 0, sizeof(simplified) );
199  poly.in_memory = arc.in_memory = 0;
200  simplified.in_memory = 1;
201
202  /* open shapefiles and dbf files */
203  if( SourceOpen( &poly, poly_file, SHPT_POLYGON ) )
204    return 1;
205  if( SourceOpen( &arc, arc_file, SHPT_ARC ) )
206    return 1;
207
208  sprintf( out_filename, "%s_p", out_file );
209  shp_out = SHPCreate( out_filename, SHPT_POLYGON );
210  if( !shp_out )
211  {
212    fprintf( stderr, "Couldn't create shapefile '%s': %s\n", out_file, strerror(errno));
213    return 1;
214  }
215 
216  dbf_out = DBFCreate( out_filename );
217  if( !dbf_out )
218  {
219    fprintf( stderr, "Couldn't create DBF '%s': %s\n", out_file, strerror(errno));
220    return 1;
221  }
222//  DBFAddField( dbf_out, "way_id", FTInteger, 11, 0 );
223//  DBFAddField( dbf_out, "orientation", FTInteger, 1, 0 );
224  DBFAddField( dbf_out, "error", FTInteger, 1, 0 );
225  DBFAddField( dbf_out, "tile_x", FTInteger, 4, 0 );
226  DBFAddField( dbf_out, "tile_y", FTInteger, 4, 0 );
227
228  /* Initialise node arrays */
229  v_x = malloc( MAX_NODES * sizeof(double) );
230  v_y = malloc( MAX_NODES * sizeof(double) );
231 
232  if( !v_x || !v_y )
233  {
234    fprintf( stderr, "Couldn't allocate memory for nodes\n" );
235    return 1;
236  }
237  /* Setup state */
238  struct state state;
239  memset( &state, 0, sizeof(state) );
240  ResizeSubareas(&state, INIT_MAX_SUBAREAS);
241 
242  // Temporary file for simplified shapes
243  {
244      // Make the file and open it
245      char filename[32];
246      sprintf( filename, "tmp_coastline_%d", getpid() );
247      simplified.shp = SHPCreate( filename, SHPT_ARC );
248      if( !simplified.shp )
249      {
250          fprintf( stderr, "Couldn't open temporary shapefile: %s\n", strerror(errno) );
251          return 1;
252      }
253      // And now unlink them so they're cleaned up when we die
254      int len = strlen( filename );
255      strcpy( filename+len, ".shp" );
256      unlink(filename);
257      strcpy( filename+len, ".shx" );
258      unlink(filename);
259     
260      // Setup the bitmap while creating the simplified polygons
261      InitBitmap( &arc, &simplified );
262      InitBitmap( &poly, &simplified );
263     
264      // And index the simplified polygons
265      simplified.shx = SHPCreateTree( simplified.shp, 2, 10, NULL, NULL );
266      if( !simplified.shx )
267      {
268          fprintf( stderr, "Couldn't build temporary shapetree\n" );
269          return 1;
270      }
271      SHPGetInfo( simplified.shp, &simplified.shape_count, NULL, NULL, NULL );
272      simplified.shapes = calloc( simplified.shape_count, sizeof(SHPObject*) );
273  }
274  /* Split coastlines into arcs no longer than MAX_NODES_PER_ARC long */
275  sprintf( out_filename, "%s_i", out_file );
276  SplitCoastlines( &poly, &arc, out_filename );
277
278#if !TEST 
279  for( int i=0; i<DIVISIONS; i++ )
280    for( int j=0; j<DIVISIONS; j++ )  //Divide the world into mercator blocks approx 100km x 100km
281#else
282  for( int i=204; i<=204; i++ )
283    for( int j=248; j<=248; j++ )  //Divide the world into mercator blocks approx 100km x 100km
284#endif
285    {
286      state.x = i;
287      state.y = j;
288     
289      double left   = -MERC_MAX + (i*MERC_BLOCK) - TILE_OVERLAP;
290      double right  = -MERC_MAX + ((i+1)*MERC_BLOCK) + TILE_OVERLAP;
291      double bottom = -MERC_MAX + (j*MERC_BLOCK) - TILE_OVERLAP;
292      double top    = -MERC_MAX + ((j+1)*MERC_BLOCK) + TILE_OVERLAP;
293     
294      if( left  < -MERC_MAX ) left  = -MERC_MAX;
295      if( right > +MERC_MAX ) right = +MERC_MAX;
296     
297      state.lb[0] = left;
298      state.lb[1] = bottom;
299      state.rt[0] = right;
300      state.rt[1] = top;
301     
302      if(isatty(STDERR_FILENO))
303//        fprintf( stderr, "\rProcessing (%d,%d)  (%.2f,%.2f)-(%.2f,%.2f)   ", i, j, left, bottom, right, top );
304        fprintf( stderr, "\rProcessing (%d,%d)  (%.2f,%.2f)-(%.2f,%.2f)   ", i, j, state.lb[0], state.lb[1], state.rt[0], state.rt[1] );
305      state.seg_count = 0;
306      state.subarea_count = 0;
307      state.subarea_nodecount = 0;
308      state.enclosed = 0;
309     
310      // Optimisation: if we have determined nothing enters the tile, we can used the simplified tiles
311      // Basically, we only need to test for enclosure
312      if( check_tile( state.x, state.y ) )
313      {
314        Process( &state, &poly, 1 );
315        Process( &state, &arc,  0 );
316      }
317      else
318        Process( &state, &simplified, 0 );
319
320      OutputSegs( &state );
321    }
322   
323  free( state.sub_areas );
324  free( simplified.shapes );
325  SHPDestroyTree( poly.shx );
326  SHPDestroyTree( arc.shx );
327  SHPDestroyTree( simplified.shx );
328  DBFClose( poly.dbf );
329  DBFClose( arc.dbf );
330  DBFClose( dbf_out );
331  SHPClose( poly.shp );
332  SHPClose( arc.shp );
333  SHPClose( simplified.shp );
334  SHPClose( shp_out );
335 
336  printf("\n");
337  return 0;
338}
339
340static void ResizeSubareas( struct state *state, int count )
341{
342  if( count < MAX_SUBAREAS )
343  {
344    fprintf( stderr, "Tried to resize smaller??? (%d < %d)\n", count, MAX_SUBAREAS );
345    exit(1);
346  }
347  fprintf( stderr, "Resizing subarea array to %d\n", count );
348 
349  struct subarea *new_sa = malloc( count * sizeof(struct subarea) );
350  free(Parts);
351  Parts = malloc( 2 * count * sizeof(int) );  // parts and Parttypes are allocated in one chunk
352 
353  if( !new_sa || !Parts )
354  {
355    fprintf( stderr, "Out of memory resizing subarea array (count=%d)\n", count );
356    exit(1);
357  }
358 
359  MAX_SUBAREAS = count;
360  memcpy( new_sa, state->sub_areas, state->subarea_count * sizeof(struct subarea) );
361  free(state->sub_areas);
362  state->sub_areas = new_sa;
363 
364  Parttypes = Parts + MAX_SUBAREAS;
365}
366
367static void SplitCoastlines2( int show, struct source *arc, SHPHandle shp_arc_out, DBFHandle dbf_arc_out )
368{
369  int count = arc->shape_count;
370  for( int i=0; i<count; i++ )
371  {
372    SHPObject *obj = source_get_shape( arc, i );
373    int way_id = DBFReadIntegerAttribute( arc->dbf, i, 2 );
374   
375    if( obj->nVertices <= MAX_NODES_PER_ARC )
376    {
377      int new_id = SHPWriteObject( shp_arc_out, -1, obj );
378      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
379      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 0, way_id );
380      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 1, show );
381      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 2, obj->nVertices < 4 );  /* Flag not real objects */
382      source_release_shape(arc, obj);
383      continue;
384    }
385    int arcs = (obj->nVertices / MAX_NODES_PER_ARC) + 1;
386    int len = (obj->nVertices / arcs) + 1;
387//    printf( "Splitting object with %d vertices, len=%d, arcs=%d\n", obj->nVertices, len, arcs );
388   
389    for( int j=0; j<arcs; j++ )
390    {
391      int this_len = (j==arcs-1)? obj->nVertices - (j*len): len+1;
392//      printf( "Subobject start=%d, length=%d\n", j*len, this_len );
393      SHPObject *new_obj = SHPCreateSimpleObject( SHPT_ARC, this_len, &obj->padfX[j*len], &obj->padfY[j*len], &obj->padfZ[j*len] );
394      int new_id = SHPWriteObject( shp_arc_out, -1, new_obj );
395      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
396      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 0, way_id );
397      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 1, show );
398      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 2, 0 );
399      SHPDestroyObject(new_obj);
400    }
401    source_release_shape(arc, obj);
402  }
403}
404
405/* The first param is currently unused, but if people ever want to get
406 * access to the completed bits of coastline, this is where to change it */
407static void SplitCoastlines( struct source *poly UNUSED, struct source *arc, char *out_filename )
408{
409  SHPHandle shp_arc_out = SHPCreate( out_filename, SHPT_ARC );
410  if( !shp_arc_out )
411  {
412    fprintf( stderr, "Couldn't create shapefile '%s': %s\n", out_filename, strerror(errno));
413    return;
414  }
415 
416  DBFHandle dbf_arc_out = DBFCreate( out_filename );
417  if( !dbf_arc_out )
418  {
419    fprintf( stderr, "Couldn't create DBF '%s': %s\n", out_filename, strerror(errno));
420    return;
421  }
422  DBFAddField( dbf_arc_out, "way_id",   FTInteger, 11, 0 );
423  DBFAddField( dbf_arc_out, "complete", FTInteger, 11, 0 );
424  DBFAddField( dbf_arc_out, "error",    FTInteger, 11, 0 );
425
426//  SplitCoastlines2( shp_poly, dbf_poly, shp_arc_out, dbf_arc_out );
427  SplitCoastlines2( 0, arc, shp_arc_out, dbf_arc_out );
428 
429  SHPClose( shp_arc_out );
430  DBFClose( dbf_arc_out );
431}
432
433static void InitBitmap( struct source *src, struct source *temp )
434{
435  int count = src->shape_count, res_count = 0;
436  int type1 = 0, type2 = 0;
437  size_t bytes_used = 0;
438 
439  if( RESOLUTION )
440    src->shapes = calloc( count, sizeof( SHPObject* ) );
441 
442  for( int i=0; i<count; i++ )
443  {
444    SHPObject *obj = source_get_shape( src, i );
445   
446    int orig = 1;  /* Does obj refer to the original shape */
447#if RESOLUTION != 0
448    /* The first loop: This one takes into account the resolution to simplify the shape */
449    if( RESOLUTION )
450    {
451      int k = 0;
452      for( int j=0; j<obj->nVertices; j++ )
453      {
454        /* If this is not the first or the last node, check that it moved a
455         * significant distance before we accept it. There is an extra
456         * buffer here to avoid lines walking an edge switching between two
457         * cells continuously */
458        if( k > 0 && j != (obj->nVertices-1) )
459        {
460          if( ( abs( v_x[k-1] - obj->padfX[j] ) < 0.75*RESOLUTION ) &&
461              ( abs( v_y[k-1] - obj->padfY[j] ) < 0.75*RESOLUTION ) )
462          {
463            type1++;  /* Type 1 optimisation */
464            continue;
465          }
466        }
467       
468        /* Calculate the position of this node */
469        double new_x = ( floor( obj->padfX[j] / RESOLUTION ) + 0.5 ) * RESOLUTION;
470        double new_y = ( floor( obj->padfY[j] / RESOLUTION ) + 0.5 ) * RESOLUTION;
471       
472        /* This optimisation is: if the step first the last node to this one
473         * is the same direction as from the last node to the one before, we
474         * can coalesce them into one step. */
475        if( k > 2 )
476        {
477          double x1 = new_x - v_x[k-1];
478          double y1 = new_y - v_y[k-1];
479          double x2 = v_x[k-1] - v_x[k-2];
480          double y2 = v_y[k-1] - v_y[k-2];
481         
482          double t1 = x1*y2;
483          double t2 = x2*y1;
484         
485          /* Check they point the same direction. The first two tests are
486           * needed other opposite directions would also match */
487          if( (x1*x2) >= 0.0f && (y1*y2) >= 0.0f && fabs(t1-t2) < 1e-6 )
488          {
489            v_x[k-1] = new_x;  /* Overwrite last node */
490            v_y[k-1] = new_y;
491            type2++;           /* Type 2 optimisation */
492            continue;
493          }
494        }
495        v_x[k] = new_x;
496        v_y[k] = new_y;
497        k++;
498       
499        if( k >= MAX_NODES-2 )
500        {
501          MAX_NODES = MAX_NODES*2;
502          v_x = realloc( v_x, MAX_NODES * sizeof(v_x[0]) );
503          v_y = realloc( v_y, MAX_NODES * sizeof(v_y[0]) );
504          fprintf( stderr, "MAX_NODES raised to %d (%ldMB)\n", MAX_NODES, ((long)MAX_NODES*2*sizeof(v_x[0]))>>20 );
505        }
506      }
507      /* If we didn't move far enough then we get skipped altogether (think small islands, low zoom) */
508      if( k > 3 )
509      {
510        /* Note we switch obj to using the simplfied shape */
511        SHPObject *new_obj = SHPCreateSimpleObject( obj->nSHPType, k, v_x, v_y, NULL );
512        bytes_used += sizeof(SHPObject) + k*sizeof(v_x[0])*2;
513        src->shapes[i] = new_obj;
514        source_release_shape(src, obj);
515        obj = new_obj;
516        orig = 0;   /* So we don't accedently free this shape at the end of the loop */
517      }
518    }
519#endif
520    /* The second loop: sets up the bitmap and makes the super simplified shapes that are only used for enclosure testing */
521    int old_x = -1, old_y = -1;
522    int k = 0;
523
524    if( obj->nVertices < 4 )
525    {
526      source_release_shape(src, obj);
527      continue;
528    }
529    for( int j=0; j<obj->nVertices; j++ )
530    {
531      /* Calculate the tile we're in */
532      int new_x = floor( (obj->padfX[j] / MERC_BLOCK) + (DIVISIONS/2.0) );
533      int new_y = floor( (obj->padfY[j] / MERC_BLOCK) + (DIVISIONS/2.0) );
534     
535      if( new_x < 0 ) new_x = 0;
536      if( new_x >= DIVISIONS ) new_x = DIVISIONS-1;
537      if( new_y < 0 ) new_y = 0;
538      if( new_y >= DIVISIONS ) new_y = DIVISIONS-1;
539     
540      if( new_x == old_x && new_y == old_y )
541          continue;
542         
543      mark_tile( new_x, new_y );
544     
545      v_x[k] = obj->padfX[j];
546      v_y[k] = obj->padfY[j];
547      k++;
548     
549      if( k > MAX_NODES-5 )
550      {
551          fprintf( stderr, "Overflow in simplification: %d > %d\n", k, MAX_NODES-5 );
552          exit(1);
553      }
554     
555      old_x = new_x;
556      old_y = new_y;
557    }
558    if( k > 4 )
559    {
560      v_x[k] = obj->padfX[obj->nVertices-1];
561      v_y[k] = obj->padfY[obj->nVertices-1];
562      k++;
563     
564      SHPObject *new_obj = SHPCreateSimpleObject( SHPT_ARC, k, v_x, v_y, NULL );
565      bytes_used += sizeof(SHPObject) + k*sizeof(v_x[0])*2; /* Note: memory not used now, but later */
566     
567//      fprintf( stderr, "Compressed %d nodes to %d\n", obj->nVertices, k );
568      int new_id = SHPWriteObject( temp->shp, -1, new_obj );
569      if( new_id < 0 ) { fprintf( stderr, "Temp output failure: %m\n"); exit(1); }
570     
571      SHPDestroyObject( new_obj );
572     
573      res_count++;
574    }
575    if( orig )
576      source_release_shape(src, obj);
577  }
578 
579  if( RESOLUTION )
580    src->in_memory = 1;
581 
582  fprintf( stderr, "%s: Found %d large objects out of %d, removed %d+%d nodes, approx memory %.1fMB\n", SHPTypeName(src->shp->nShapeType), res_count, count, type1, type2, bytes_used/(1024.0*1024.0) );
583}
584
585static double CalcArea( const SHPObject *obj )
586{
587  int i;
588  double base_x = obj->dfXMin;
589  double base_y = obj->dfYMin;
590  double area = 0;
591  int n = (obj->nParts <= 1) ? obj->nVertices : obj->panPartStart[1];
592//  if(VERBOSE) printf( "CalcArea: n=%d\n", n );
593  for( i=0; i<n; i++ )
594  {
595    int p = ((i==0)?n:i)-1;
596    double x1 = obj->padfX[p] - base_x;
597    double x2 = obj->padfX[i] - base_x;
598    double y1 = obj->padfY[p] - base_y;
599    double y2 = obj->padfY[i] - base_y;
600   
601//    if(VERBOSE) printf( "i=%d (%f,%f)-(%f,%f) %f\n", i, x1,y1,x2,y2,x1*y2 - y1*x2);
602    area += x1*y2 - y1*x2;
603  }
604  return area;
605}
606
607static const int table[3][3] = { {6, 5, 4}, {7, -1, 3}, { 0, 1, 2 } };
608
609/* Determines the quadrant the given point is in. -1 means it's inside the box */
610/* The rule is: left and bottom edges are in, right and top edges are out */
611static inline int GetPosition( double *lb, double *rt, double X, double Y )
612{
613  int x, y;
614
615  x = (X >= rt[0]) - ( X < lb[0] ) + 1;
616  y = (Y >= rt[1]) - ( Y < lb[1] ) + 1;
617
618  return table[x][y];
619}
620
621struct intersect
622{
623  double x, y;  /* Coordinates */
624  double t;     /* 0..1, where the intersection happens */
625};
626
627/* The rule is: left and bottom edges are in, right and top edges are out */
628/* All corners are out, except the lower left one */
629static int CalculateIntersections( double x1, double y1, double x2, double y2, 
630                                   double *lb, double *rt, struct intersect *intersections )
631{
632  int count = 0;
633  double x, y, t;
634 
635  /* Left side */
636  if( (x1 < lb[0]) != (x2 < lb[0]) )
637  {
638    /* Determine intersection */
639    x = lb[0];
640    t = (x1-lb[0]) / (x1-x2);
641    y = y1 + t*(y2-y1);
642   
643    if( lb[1] <= y && y < rt[1] ) /* Include only if in range */
644      intersections[count++] = (struct intersect){ x: x, y: y, t: t };
645  }
646  /* Right side */
647  if( (x1 >= rt[0]) != (x2 >= rt[0]) )
648  {
649    /* Determine intersection */
650    x = rt[0];
651    t = (x1-rt[0]) / (x1-x2);
652    y = y1 + t*(y2-y1);
653   
654    if( lb[1] <= y && y < rt[1] ) /* Include only if in range */
655      intersections[count++] = (struct intersect){ x: x, y: y, t: t };
656  }
657  /* Top side */
658  if( (y1 >= rt[1]) != (y2 >= rt[1]) )
659  {
660    /* Determine intersection */
661    y = rt[1];
662    t = (y1-rt[1]) / (y1-y2);
663    x = x1 + t*(x2-x1);
664   
665    if( lb[0] <= x && x < rt[0] ) /* Include only if in range */
666      intersections[count++] = (struct intersect){ x: x, y: y, t: t };
667  }
668  /* Bottom side */
669  if( (y1 < lb[1]) != (y2 < lb[1]) )
670  {
671    /* Determine intersection */
672    y = lb[1];
673    t = (y1-lb[1]) / (y1-y2);
674    x = x1 + t*(x2-x1);
675   
676    if( lb[0] <= x && x < rt[0] ) /* Include only if in range */
677      intersections[count++] = (struct intersect){ x: x, y: y, t: t };
678  }
679
680  /* Check the count, if we went over we killed the caller's stack... */
681  if( count > 2 )
682  {
683    fprintf( stderr, "Too many intersections (%d)\n", count );
684    exit(1);
685  }
686 
687  if( count == 2 )
688  {
689    /* If there's two intersections, reorder them to match the intersection order */
690    if( intersections[0].t > intersections[1].t )
691    {
692      struct intersect i;
693      i = intersections[1];
694      intersections[1] = intersections[0];
695      intersections[0] = i;
696    }
697  }
698#if 0
699  if(count == 0)
700  {
701    printf( "\nCalculate intersections: (%.2f,%.2f)-(%.2f,%.2f) hit %d\n", x1, y1, x2, y2, count );
702    for( int i=0; i<count; i++ )
703      printf( "   (%.2f,%.2f) t=%.6f\n", intersections[i].x, intersections[i].y, intersections[i].t );
704  } 
705#endif
706  return count;
707}
708
709static int seg_compare( const void *a, const void *b )
710{
711  const struct segment *aa = (struct segment*)a;
712  const struct segment *bb = (struct segment*)b;
713 
714  if( aa->sa < bb->sa )
715    return -1;
716  if( aa->sa > bb->sa )
717    return 1;
718  return 0;
719}
720
721/* We currently don't use anything from the source DBF file, but the cabability is there */
722void Process( struct state *state, struct source *src, int polygon )
723{
724  int count;
725  int *list = SHPTreeFindLikelyShapes( src->shx, state->lb, state->rt, &count );
726  int poly_start;
727 
728  for( int poly = 0; poly < count; poly++ )
729  {
730    /* Here we track parts that have gone across the box */
731    int vertex;
732    int intersected = 0;
733   
734    int id = list[poly];
735
736    /* Now we have a candidate object, we need to process it */
737    SHPObject *obj = source_get_shape( src, id );
738    if( !obj ) /* May be NULL if this object was simplfied away */
739      continue;
740   
741    /* If it's got less than 4 vertices it's not a real object */
742    /* No need to mark it as error here, done in SplitCoastlines */
743    if( obj->nVertices < 4 )
744    {
745      source_release_shape(src, obj );
746      continue;
747    }
748
749    if( polygon &&
750        state->lb[0] < obj->dfXMin && obj->dfXMax < state->rt[0] &&
751        state->lb[1] < obj->dfYMin && obj->dfYMax < state->rt[1] )
752    {
753      int sa = state->subarea_count;
754      if( sa > MAX_SUBAREAS-5 )
755        ResizeSubareas( state, 2*MAX_SUBAREAS );
756
757      state->sub_areas[sa].src = src;
758      state->sub_areas[sa].index = id;
759      state->sub_areas[sa].x = obj->padfX[0];
760      state->sub_areas[sa].y = obj->padfY[0];
761      state->sub_areas[sa].areasize = CalcArea( obj );
762      state->sub_areas[sa].used = 0;
763
764      state->subarea_count++;
765      state->subarea_nodecount += obj->nVertices;
766      source_release_shape(src, obj );
767      continue;
768    }
769 
770    if(VERBOSE) fprintf( stderr, "\nProcessing object %d (%d vertices)\n", poly, obj->nVertices );
771    /* First we need to move along the object until we leave the box. For polygons we
772     * know it will eventually, since we determined already this object is larger
773     * than the box */
774
775    for( vertex=0; vertex < obj->nVertices && GetPosition( state->lb, state->rt, obj->padfX[vertex], obj->padfY[vertex] ) == -1; vertex++ )
776      ;
777
778    if( vertex == obj->nVertices )
779    {
780      if( polygon )
781        fprintf( stderr, "Object %d did not leave box (%d vertices, polygon:%d) (%.2f,%.2f-%.2f,%.2f)\n", id, 
782                obj->nVertices, polygon, obj->dfXMin, obj->dfYMin, obj->dfXMax, obj->dfYMax );
783      source_release_shape(src, obj );
784      continue;
785    }
786    /* We need to mark this point, so when we loop back we know where to stop */
787    poly_start = vertex+1;
788   
789    /* This tracks the current sector. */
790    int curr_sect = GetPosition( state->lb, state->rt, obj->padfX[vertex], obj->padfY[vertex] );
791   
792    int winding = 0, max_winding = 0, min_winding = 0;
793   
794    /* We need this flag so when we loop round a polygon we know if we're at the first node or the last */
795    int started = 0;
796
797    /* The basic trick is to analyse each new point and determine the
798     * sector. As long as the sector doesn't change, we're cool. They're
799     * numbered in such a way that any line has to go through the sectors in
800     * numerical order (modulo 8) *or* it intersects the box. The
801     * reason we do this is that if the line never intersects the box, we need
802     * to determine the "winding number" determine if we're inside or outside. */
803
804    if( poly_start == obj->nVertices )
805      poly_start = 1;
806//    fprintf( stderr, "poly_start=%d, vertex=%d\n", poly_start, vertex );
807    for(;;)
808    {
809      vertex++;
810      /* First we need to handle the step to the next node */
811      if( polygon )
812      {
813        /* For polygons we loop around to the start point */
814        if( vertex == obj->nVertices )
815        {
816//          fprintf( stderr, "\nLooping..." );
817//          sleep(1);
818          vertex = 1;
819        }
820        if( vertex == poly_start && started )
821          break;
822      }
823      else
824      {
825        /* Else we just stop at the end */
826        if( vertex == obj->nVertices )
827          break;
828      }
829      started = 1;
830     
831      if( vertex >= obj->nVertices )  /* Shouldn't happen */
832      {
833        fprintf( stderr, "Somehow %d >= %d\n", vertex, obj->nVertices );
834        break;
835      }
836     
837      int sect = GetPosition( state->lb, state->rt, obj->padfX[vertex], obj->padfY[vertex] );
838      if( sect == curr_sect )
839        continue;
840      if(VERBOSE)
841      printf("Moved from %d to %d\n", curr_sect, sect );
842      /* Now we know we've moved to another sector, so we need to know the intersection points */
843      struct intersect intersections[2];
844      int int_count;
845
846//      if( id == 133912 ) fprintf(stderr, "\nVertex %d: Moved from sector %d to %d", vertex, curr_sect, sect );
847     
848      /* If we moved to adjacent positive sector, we don't need to check for intersections... */
849      if( sect >= 0 && curr_sect >= 0 && ( ((curr_sect-sect)&7) == 1 || ((sect-curr_sect)&7) == 1 ) )
850        int_count = 0;
851      else
852        int_count = CalculateIntersections( obj->padfX[vertex-1], obj->padfY[vertex-1], 
853                                            obj->padfX[vertex],   obj->padfY[vertex],
854                                            state->lb, state->rt,
855                                            intersections );
856
857      /* There are corner cases with the calculations of intersections, if
858       * you move exactly on to the edge. With floating point numbers you'd
859       * think it was possibly, but, well, it is. Espescially around the 0
860       * meridian. In this case we get zero intersections even though we
861       * changed inside/outside. What we do is basically pretend we havn't
862       * changed sector at all and check the next point */
863      if( int_count == 0 && (sect == -1 || curr_sect == -1) )
864      {
865        fprintf( stderr, "Went from %d to %d without intersection.\n"
866                      "line (%f,%f)-(%f,%f), box (%f,%f)-(%f,%f)\n",
867                      curr_sect, sect,
868                      obj->padfX[vertex-1], obj->padfY[vertex-1],
869                      obj->padfX[vertex],   obj->padfY[vertex],
870                      state->lb[0], state->lb[1], state->rt[0], state->rt[1] );
871      }
872      /* Another possibility is that we went from positive to positive with
873       * only one intersection. Not good. Have not yet thought of a
874       * heuristic to deal with this case */
875      if( int_count == 1 && ( (sect != -1) == (curr_sect != -1) ) )
876      {
877        fprintf( stderr, "Went from %d to %d with 1 intersection.\n"
878                      "line (%f,%f)-(%f,%f), box (%f,%f)-(%f,%f)\n",
879                      curr_sect, sect,
880                      obj->padfX[vertex-1], obj->padfY[vertex-1],
881                      obj->padfX[vertex],   obj->padfY[vertex],
882                      state->lb[0], state->lb[1], state->rt[0], state->rt[1] );
883      }
884     
885      /* finally, if we got two intersections, we mave have passed straight through the box */
886      if( int_count == 2 && (sect == -1 || curr_sect == -1) )
887      {
888        fprintf( stderr, "Went from %d to %d with 2 intersections.\n"
889                      "line (%f,%f)-(%f,%f), box (%f,%f)-(%f,%f)\n",
890                      curr_sect, sect,
891                      obj->padfX[vertex-1], obj->padfY[vertex-1],
892                      obj->padfX[vertex],   obj->padfY[vertex],
893                      state->lb[0], state->lb[1], state->rt[0], state->rt[1] );
894      }
895     
896      struct segment *seg_ptr = &state->seg_list[state->seg_count];
897     
898      /* Now we know the number of intersections. */
899      if( int_count == 0 )
900      {
901        /* If we have no intersections, we were outside and we're still
902         * outside. Then we need to track the winding number. With no
903         * intersections we can move a maximum of three sections, so we can
904         * tell if we moved clockwise or anti-clockwise. */
905        if( !intersected )
906        {
907          int diff = (sect - curr_sect) & 7;  /* Like mod 8, but always positive */
908          if( diff > 4 )
909            diff -= 8;
910           
911          winding += diff;
912         
913          if( winding < min_winding )
914            min_winding = winding;
915          if( winding > max_winding )
916            max_winding = winding;
917        }
918      }
919      else if( int_count == 1 )
920      {
921        /* If we have one intersection, we went from inside to out, or vice
922         * versa. So we need to add a section or finish an old one. */
923        intersected = 1;
924        if( curr_sect != -1 )
925        {
926          /* Going in... */
927          seg_ptr->src = src;
928          seg_ptr->index = id;
929          seg_ptr->sx = intersections[0].x;
930          seg_ptr->sy = intersections[0].y;
931          seg_ptr->snode = vertex;
932        }
933        else
934        {
935          /* Going out... */
936          seg_ptr->ex = intersections[0].x;
937          seg_ptr->ey = intersections[0].y;
938          seg_ptr->enode = vertex-1;
939         
940          (state->seg_count)++;
941        }
942      }
943      else if( int_count == 2 )
944      {
945        /* If we have two intersections, we went straight through. So we need
946         * to make a segment with no internal nodes */
947        intersected = 1;
948       
949        /* Going in and out in one go... */
950        seg_ptr->src = NULL;
951        seg_ptr->index = -2;
952        seg_ptr->sx = intersections[0].x;
953        seg_ptr->sy = intersections[0].y;
954        seg_ptr->snode = -1;
955        seg_ptr->ex = intersections[1].x;
956        seg_ptr->ey = intersections[1].y;
957        seg_ptr->enode = -1;
958       
959        (state->seg_count)++;
960      }
961      // We need 6 spare for possible enclosure test later
962      if( state->seg_count > (MAX_SEGS-6) )
963      {
964          fprintf( stderr, "MAX_SEGS overflow!!!, aborting\n" );
965          exit(1);
966      }
967      curr_sect = sect;
968    }
969    /* Generally, if we have a high winding number we consider ourselves
970     * inside. However, in the case of bits that stick out, while the ends
971     * may show to not have a high winding number, the min/max winding
972     * numbers can show that the shape made a loop around the box, even
973     * though the endpoints may not look like it. So we also trigger if
974     * the end winding number >3 and the greatest difference at least 6 */
975   
976    /* However, it's more complicated than that. We could be enclosed by a
977     * continental landmass, but also by a large lake. So we must take the
978     * *smallest* enclosing area to determine whether we're enclosed or not.
979     * To facilitate this the enclosed fields tracks the area of the
980     * smallest encloser and the sign determines te positive or negativeness
981     * of the area. */
982   
983    if( !intersected )
984    {
985      double enclosed = 0;
986      if ((winding >= 3) || (winding >= 2 && (max_winding-min_winding) > 6))
987      {
988        if(VERBOSE)
989          printf( "Decided enclosed: intersected=%d, winding=%d, max_winding=%d, min_winding=%d\n", intersected, winding, max_winding, min_winding );
990        enclosed = +1;
991      }
992      if (winding < -4) 
993      {
994        if(VERBOSE)
995          printf( "Decided unenclosed: intersected=%d, winding=%d, max_winding=%d, min_winding=%d\n", intersected, winding, max_winding, min_winding );
996        enclosed = -1;
997      }
998
999      if( enclosed )
1000      {
1001        // Here we don't need to worry too much about small areas, since
1002        // they're unlikely to enclose anything
1003        double size = CalcArea(obj);
1004        int intsize = copysign( ceil(10*log2(1+fabs(size))), enclosed ); // Scale the size down to a number that will fit in an integer
1005        if( state->enclosed == 0 )
1006          state->enclosed = intsize;
1007        else if( abs(state->enclosed) > abs(intsize) )
1008          state->enclosed = intsize;
1009         
1010        if(VERBOSE)
1011          printf( "(%d,%d) New state->enclosed: %d, size=%f\n", state->x, state->y, state->enclosed, size );
1012      }
1013    }
1014   
1015    source_release_shape(src, obj);
1016  }
1017  free(list);
1018//  printf( "segcount: %d\n", state->seg_count );
1019}
1020
1021/* Compare function to sort subareas in *descending* order */
1022static int subarea_compare( const void *a, const void *b )
1023{
1024  struct subarea *sa = (struct subarea *)a;
1025  struct subarea *sb = (struct subarea *)b;
1026 
1027  double diff = fabs(sb->areasize) - fabs(sa->areasize);
1028 
1029  if( diff > 0 ) return +1;
1030  if( diff < 0 ) return -1;
1031  return 0;
1032}
1033
1034void OutputSegs( struct state *state )
1035{
1036  /* At this point we've processed the whole object and have a list of
1037   * segments which cross our box. We now compute the angles they make to
1038   * the centre, sort them on that basis and try to make closed areas. There
1039   * are as many starts as finishes so it will always finish, but a senseble
1040   * end result depends on the original polygon having been sane (i.e.
1041   * simple) */
1042   
1043  double centre_x = (state->lb[0] + state->rt[0])/2;
1044  double centre_y = (state->lb[1] + state->rt[1])/2;
1045
1046  struct segment *seg_list = state->seg_list;
1047  int seg_count = state->seg_count;
1048 
1049//  fprintf( stderr, "%d sub %d  ", state->subarea_count, state->subarea_nodecount );
1050  // First we must sort the subareas by decreasing size
1051  qsort( state->sub_areas, state->subarea_count, sizeof( state->sub_areas[0] ), subarea_compare );
1052 
1053  if(VERBOSE)
1054    printf( "enclosed=%d, seg_count=%d, subarea_count=%d\n", state->enclosed, seg_count, state->subarea_count );
1055   
1056  if( seg_count == 0 )
1057  {
1058    /* No intersections at all, so we check the winding number. With the
1059     * water-on-the-right rule we're looking for a positive winding. */
1060//      if( abs(winding) > 4 )
1061//        printf( "\nNot intersected, winding = %d, min_winding = %d, max_winding = %d\n", winding, min_winding, max_winding );
1062     
1063    if( state->enclosed > 0 )  // Enclosed by a negative area does not count
1064    {
1065      seg_list[0].sx = seg_list[0].ex = state->lb[0];
1066      seg_list[0].sy = seg_list[0].ey = centre_y;
1067      seg_list[0].snode = seg_list[0].enode = -1;
1068      seg_list[0].sa = +M_PI;
1069      seg_list[0].ea = -M_PI;
1070      seg_list[0].next = -1;
1071      seg_count++;
1072    }
1073  }
1074  else
1075  {
1076    /* Work out the angles */
1077    for( int k=0; k<seg_count; k++ )
1078    {
1079      seg_list[k].sa = atan2( seg_list[k].sy - centre_y, seg_list[k].sx - centre_x );
1080      seg_list[k].ea = atan2( seg_list[k].ey - centre_y, seg_list[k].ex - centre_x );
1081      seg_list[k].next = -1;
1082    }
1083  }
1084  if( seg_count > 0 )
1085  {
1086    /* Create the helper nodes for the corners */
1087    for( int k=0; k<4; k++, seg_count++ )
1088    {
1089      seg_list[seg_count].sx = seg_list[seg_count].ex = (k<2) ? state->lb[0] : state->rt[0];
1090      seg_list[seg_count].sy = seg_list[seg_count].ey = (((k+1)&3)<2) ? state->lb[1] : state->rt[1];
1091      seg_list[seg_count].sa = seg_list[seg_count].ea = atan2( seg_list[seg_count].sy - centre_y, seg_list[seg_count].sx - centre_x );
1092      seg_list[seg_count].snode = seg_list[seg_count].enode = -2;
1093      seg_list[seg_count].next = -2;
1094    }
1095   
1096    /* Sort the nodes by increasing angle */
1097    qsort( seg_list, seg_count, sizeof(seg_list[0]), seg_compare );
1098   
1099    for(;;)
1100    {
1101      int part_count = 1;
1102      Parts[0] = 0;
1103      Parttypes[0] = SHPP_RING;
1104     
1105      /* First we need to find an unused segment */
1106      int curr;
1107      for( curr=0; curr < seg_count && seg_list[curr].next != -1; curr++ )
1108        ;
1109      if( curr == seg_count )
1110        break;
1111      int node_count = 0;
1112      for(;;)
1113      {
1114        if(VERBOSE) printf( "Part %d: ndc=%d, curr=%d (%d-%d)\n", part_count, node_count, curr, seg_list[curr].snode, seg_list[curr].enode );
1115        v_x[node_count] = seg_list[curr].sx;
1116        v_y[node_count] = seg_list[curr].sy;
1117        node_count++;
1118
1119        if( seg_list[curr].snode >= 0 )
1120        {
1121          node_count = CopyShapeToArray( seg_list[curr].src, seg_list[curr].index, seg_list[curr].snode, seg_list[curr].enode, node_count );
1122          v_x[node_count] = seg_list[curr].ex;
1123          v_y[node_count] = seg_list[curr].ey;
1124          node_count++;
1125        }
1126       
1127        double angle = seg_list[curr].ea;
1128        /* Determine the first unused segment with a start angle greater than the current angle, with wrapping */
1129        int next;
1130        for( next = 0; next < seg_count && seg_list[next].sa <= angle; next++ )
1131          ;
1132        if( next == seg_count )
1133          next = 0;
1134         
1135        seg_list[curr].next = next;
1136        /* If we come to an already used segment, we're done */
1137        if( seg_list[next].next >= 0 )
1138          break;
1139         
1140        curr = next;
1141      }
1142      v_x[node_count] = v_x[0];
1143      v_y[node_count] = v_y[0];
1144      node_count++;
1145     
1146      ProcessSubareas( state, &part_count, &node_count );
1147     
1148      if( part_count > MAX_SUBAREAS - 2 )
1149        fprintf( stderr, "(%d,%d) Subarea overflow: %d > %d\n", state->x, state->y, part_count, MAX_SUBAREAS-2 );
1150      if( node_count > MAX_NODES - 100 )
1151        fprintf( stderr, "(%d,%d) Node overflow: %d > %d\n", state->x, state->y, node_count, MAX_NODES - 100 );
1152//        fprintf( stderr, "Created object: %d verticies\n", node_count );
1153//      SHPObject *shape = SHPCreateSimpleObject( SHPT_POLYGON, node_count, v_x, v_y, NULL );
1154      SHPObject *shape = SHPCreateObject( SHPT_POLYGON, -1, 
1155                                            part_count, Parts, Parttypes, 
1156                                            node_count, v_x, v_y, NULL, NULL );
1157      // If a wrongly oriented shape crosses a boundary, sometimes we can see that...
1158      // Must do this prior to rewinding object
1159      int inverted = (CalcArea( shape ) < 0);
1160      if(VERBOSE) printf( "Created shape orientation: %f,%s\n", CalcArea( shape ), !inverted?"good":"bad" );
1161      SHPRewindObject( NULL, shape );
1162      int new_id = SHPWriteObject( shp_out, -1, shape );
1163      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
1164      SHPDestroyObject( shape );
1165     
1166      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_ERROR, inverted ? 1 : 0 );
1167      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_X, state->x );
1168      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_Y, state->y );
1169    }
1170  }
1171  /* Check for any remaining sub areas... we just output them, but mark negative ones as errors */
1172  for( int k=0; k<state->subarea_count; k++ )
1173  {
1174    struct subarea *sub_area = &state->sub_areas[k];
1175    if( sub_area->used )
1176      continue;
1177     
1178    if(VERBOSE) printf( "Remaining subarea %d (area=%f)\n", k, sub_area->areasize );
1179    int node_count = CopyShapeToArray( sub_area->src, sub_area->index, 0, -1, 0 );
1180    int part_count = 1;
1181    Parttypes[0] = SHPP_RING;
1182   
1183    sub_area->used = 1;  // Need to mark as used first, or it will match itself
1184
1185    ProcessSubareas( state, &part_count, &node_count );
1186   
1187//      SHPObject *obj = SHPReadObject( sub_area->src, sub_area->index );
1188    SHPObject *obj = SHPCreateObject( SHPT_POLYGON, -1, 
1189                                      part_count, Parts, Parttypes, 
1190                                      node_count, v_x, v_y, NULL, NULL );
1191    int new_id = SHPWriteObject( shp_out, -1, obj );
1192    if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
1193    SHPDestroyObject( obj );
1194    DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_ERROR, (sub_area->areasize > 0)?0:1 );
1195    DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_X, state->x );
1196    DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_Y, state->y );
1197  }
1198//    printf( "\nMatching %s with %d vertices, id=%d\n", polygon ? "polygon" : "arc", obj->nVertices, id );
1199  if(VERBOSE)
1200  for( int k=0; k<seg_count; k++ )
1201  {
1202    printf( "%2d %6d (%11.2f,%11.2f)-(%11.2f,%11.2f) (%5d-%5d) (%4.2f-%4.2f) => %2d\n", 
1203        k, seg_list[k].index, seg_list[k].sx, seg_list[k].sy, seg_list[k].ex, seg_list[k].ey, 
1204        seg_list[k].snode, seg_list[k].enode, seg_list[k].sa, seg_list[k].ea,
1205        seg_list[k].next );
1206  }
1207}
1208
1209static int contains( double x, double y, double *v_x, double *v_y, int vertices )
1210{
1211#if 1
1212  /* Referenced from Wikipedia */
1213  /* http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html */
1214    {
1215      int i, j, c = 0;
1216      for (i = 0, j = vertices-1; i < vertices; j = i++) {
1217        if ((((v_y[i]<=y) && (y<v_y[j])) ||
1218             ((v_y[j]<=y) && (y<v_y[i]))) &&
1219            (x < (v_x[j] - v_x[i]) * (y - v_y[i]) / (v_y[j] - v_y[i]) + v_x[i]))
1220
1221          c = !c;
1222      }
1223      return c;
1224    }
1225#else
1226  /* Home cooked algorithm, with a bug */
1227  double last_angle = atan2( v_y[0] - y, v_x[0] - x );
1228  double angle = 0;
1229 
1230#if TEST
1231  int debug = (round(x) == round(6538707.962709)) && (round(y) == round(2365505.004066));
1232#endif
1233
1234  double last_x = v_x[0];
1235  double last_y = v_y[0];
1236 
1237  for( int i=0; i<vertices; i++ )
1238  {
1239    /* Only bother checking if we've changed quadrants */
1240    if( i != vertices-1 &&
1241        (v_x[i] < x) == (last_x < x) &&
1242        (v_y[i] < y) == (last_y < y) )
1243      continue;
1244    last_x = v_x[i];
1245    last_y = v_y[i];
1246   
1247    double new_angle = atan2( v_y[i] - y, v_x[i] - x );
1248   
1249#if TEST
1250    if(debug)
1251      printf( "i=%d (%f,%f)-(%f,%f) last_angle=%.2f new_angle=%.2f angle=%.2f offset=%.2f\n", i, v_x[i-1], v_y[i-1], last_x, last_y, last_angle, new_angle, angle, new_angle-angle );
1252#endif
1253    /* What we want to do is set angle to the same angle as new_angle module
1254     * 2pi in such a way that it is within pi of the current angle. */
1255    double diff = fmod(fmod(new_angle-last_angle,2*M_PI) + 3*M_PI,2*M_PI) - M_PI;
1256    last_angle = new_angle;
1257    angle += diff;
1258   
1259  }
1260  double windings = round( angle / (2*M_PI) );
1261
1262#if TEST
1263  if(debug)
1264    printf( "Determined for (%.2f,%.2f) windings=%f\n", x, y, windings );
1265#endif
1266  return windings > 0;
1267#endif
1268}
1269
1270static int CopyShapeToArray( struct source *src, int index, int snode, int enode, int node_count )
1271{
1272  SHPObject *obj = source_get_shape( src, index );
1273  int k;
1274 
1275  if( enode == -1 )   /* This mean to copy the whole object */
1276    enode = obj->nVertices - 1;
1277  //          printf( "Copying %d - %d (max %d)\n", seg_list[curr].snode, seg_list[curr].enode, obj->nVertices );
1278 
1279  /* The slightly odd construction of the loop is because it need to handle sections of polygon which may wrap around the end of the shape */
1280  if( snode != enode && enode == 0 )
1281    enode = obj->nVertices-1;
1282  // nodes[0] == node[nVertices-1]
1283  // nodes[1] == node[nVertices]
1284  for( k = snode; k != enode; k++ )
1285  {
1286    if( node_count > MAX_NODES - 100 )
1287    {
1288//      fprintf( stderr, "Node overflow...\n" );
1289      break;
1290    }
1291    if( k == obj->nVertices )
1292    {
1293      k = 1;
1294      if( k == enode )
1295        break;
1296    }
1297     
1298    v_x[node_count] = obj->padfX[k];
1299    v_y[node_count] = obj->padfY[k];
1300    node_count++;
1301  }
1302  /* k = enode now... */
1303  v_x[node_count] = obj->padfX[k];
1304  v_y[node_count] = obj->padfY[k];
1305  node_count++;
1306  source_release_shape(src, obj);
1307 
1308  return node_count;
1309}
1310
1311void ProcessSubareas(struct state *state, int *pc, int *nc)
1312{
1313  int part_count = *pc;
1314  int node_count = *nc;
1315 
1316  if(VERBOSE) printf( "ProcessSubareas(pc=%d,nc=%d)\n", *pc, *nc );
1317 
1318  // Setting this simplifies the later code as it can be assumed that
1319  // every ring ends at the beginning of the next part
1320  Parts[part_count] = node_count;
1321  // Now we have the outer ring we need to find if any of the subareas are contained in it
1322  for( int k=0; k<state->subarea_count; k++ )
1323  {
1324    struct subarea *sub_area = &state->sub_areas[k];
1325    if( sub_area->used )  // Already used, skip
1326      continue;
1327    if( !contains( sub_area->x, sub_area->y, v_x, v_y, Parts[1] ) ) // If not contained skip
1328      continue;
1329    // Now we have to verify that this new subarea is not contained in any of the existing ones
1330    if(VERBOSE) printf( "subarea %d is contained (size=%f,%s)\n", k, sub_area->areasize, (sub_area->areasize < 0)?"good":"bad");
1331    int contained = 1;
1332    for( int m=1; m < part_count; m++ )
1333    {
1334      if( contains( sub_area->x, sub_area->y, v_x+Parts[m], v_y+Parts[m], Parts[m+1]-Parts[m] ) )
1335      {
1336        contained = 0;
1337        if(VERBOSE) printf( "subarea %d is excluded\n", k);
1338        break;
1339      }
1340    }
1341    if( !contained )
1342      continue;
1343     
1344    // We know the surrounding object has positive area
1345    if( sub_area->areasize < 0 )
1346    {
1347//          printf( "Appending subshape\n" );
1348      Parttypes[0] = SHPP_OUTERRING; /* Multipart object now */
1349      // Append to object
1350      if( part_count >= MAX_SUBAREAS-2 )
1351      {
1352       // ResizeSubareas( state, MAX_SUBAREAS*2 )
1353        fprintf( stderr, "(%d,%d) Parts array overflow (%d)\n", state->x, state->y, MAX_SUBAREAS );
1354      }
1355      else
1356      {
1357        Parttypes[part_count] = SHPP_INNERRING;
1358        Parts[part_count] = node_count;
1359        node_count = CopyShapeToArray( sub_area->src, sub_area->index, 0, -1, node_count );
1360        Parts[part_count+1] = node_count;
1361      }
1362     
1363      part_count++;
1364    }
1365    else
1366    {
1367      // Error, copy object
1368      SHPObject *obj = source_get_shape( sub_area->src, sub_area->index );
1369      if(VERBOSE) printf( "Outputting error shape (x,y)=(%f,%f), n=%d\n", obj->padfX[0], obj->padfY[0], obj->nVertices );
1370      int new_id = SHPWriteObject( shp_out, -1, obj );
1371      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
1372      source_release_shape(sub_area->src, obj );
1373     
1374      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_ERROR, 1 );
1375      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_X, state->x );
1376      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_Y, state->y );
1377    }
1378     
1379    sub_area->used = 1;
1380  }
1381 
1382  *nc = node_count;
1383  *pc = part_count;
1384}
1385
1386static int SourceOpen( struct source *src, char *filename, int shapetype )
1387{
1388  src->shp = SHPOpen( filename, "rb" );
1389  if( !src->shp )
1390  {
1391    fprintf( stderr, "Couldn't open '%s': %s\n", filename, strerror(errno));
1392    return 1;
1393  }
1394
1395  src->dbf = DBFOpen( filename, "rb" );
1396  if( !src->dbf )
1397  {
1398    fprintf( stderr, "Couldn't open DBF file for '%s'\n", filename );
1399    return 1;
1400  }
1401  if( DBFGetFieldIndex( src->dbf, "way_id" ) != 2 )
1402  {
1403    fprintf( stderr, "Unexpected format DBF file '%s'\n", filename );
1404    return 1;
1405  }
1406
1407  /* Check shapefiles are the right type and initialise length */
1408  int type;
1409  SHPGetInfo( src->shp, &src->shape_count, &type, NULL, NULL );
1410  if( type != shapetype )
1411  {
1412    fprintf( stderr, "'%s' is not a %s shapefile\n", filename, SHPTypeName(shapetype) );
1413    return 1;
1414  }
1415  /* Build indexes on files, we need them... */
1416  src->shx = SHPCreateTree( src->shp, 2, 10, NULL, NULL );
1417  if( !src->shx )
1418  {
1419    fprintf( stderr, "Couldn't create shape indexes\n" );
1420    return 1;
1421  }
1422  return 0;
1423}
Note: See TracBrowser for help on using the repository browser.