1 /******************************************************************************
2 *
3 * Project: OpenCPN
4 * Purpose: Tesselated Polygon Object
5 * Author: David Register
6 *
7 ***************************************************************************
8 * Copyright (C) 2010 by David S. Register *
9 * *
10 * This program is free software; you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 * This program is distributed in the hope that it will be useful, *
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18 * GNU General Public License for more details. *
19 * *
20 * You should have received a copy of the GNU General Public License *
21 * along with this program; if not, write to the *
22 * Free Software Foundation, Inc., *
23 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. *
24 ***************************************************************************
25 *
26 *
27 */
28 // For compilers that support precompilation, includes "wx.h".
29 #include "wx/wxprec.h"
30
31 #ifndef WX_PRECOMP
32 #include "wx/wx.h"
33 #endif //precompiled headers
34
35 #include "wx/tokenzr.h"
36 #include <wx/mstream.h>
37
38 #include "config.h"
39
40 #include "gdal/ogr_geometry.h"
41
42 #include "cutil.h"
43
44 #include "vector2D.h"
45
46 #include "s52s57.h"
47
48 #include "mygeom.h"
49 #include "georef.h"
50
51 #include "dychart.h"
52
53 #include <stdio.h>
54 #include <stdlib.h>
55 #include <string.h>
56 #include <math.h>
57
58 #include "tesselator.h"
59 #include "Striper.h"
60
61 #ifdef __WXMSW__
62 #include <windows.h>
63 #endif
64
65 static const double CM93_semimajor_axis_meters = 6378388.0; // CM93 semimajor axis
66
67
68
69 // Module Internal Prototypes
70
71
72 #if defined( __UNIX__ ) && !defined(__WXOSX__) // high resolution stopwatch for profiling
73 class OCPNStopWatchTess
74 {
75 public:
OCPNStopWatchTess()76 OCPNStopWatchTess() { Reset(); }
Reset()77 void Reset() { clock_gettime(CLOCK_REALTIME, &tp); }
78
GetTime()79 double GetTime() {
80 timespec tp_end;
81 clock_gettime(CLOCK_REALTIME, &tp_end);
82 return (tp_end.tv_sec - tp.tv_sec) * 1.e3 + (tp_end.tv_nsec - tp.tv_nsec) / 1.e6;
83 }
84
85 private:
86 timespec tp;
87 };
88 #endif
89
90
91 //-----------------------------------------------------------
92 //
93 // Static memory allocators for libtess2
94 //
95 //-----------------------------------------------------------
96 //#define USE_POOL 1
97
stdAlloc(void * userData,unsigned int size)98 void* stdAlloc(void* userData, unsigned int size)
99 {
100 int* allocated = ( int*)userData;
101 TESS_NOTUSED(userData);
102 *allocated += (int)size;
103 return malloc(size);
104 }
105
stdFree(void * userData,void * ptr)106 void stdFree(void* userData, void* ptr)
107 {
108 TESS_NOTUSED(userData);
109 free(ptr);
110 }
111
112 struct MemPool
113 {
114 unsigned char* buf;
115 unsigned int cap;
116 unsigned int size;
117 };
118
poolAlloc(void * userData,unsigned int size)119 void* poolAlloc( void* userData, unsigned int size )
120 {
121 struct MemPool* pool = (struct MemPool*)userData;
122 size = (size+0x7) & ~0x7;
123 if (pool->size + size < pool->cap)
124 {
125 unsigned char* ptr = pool->buf + pool->size;
126 pool->size += size;
127 return ptr;
128 }
129 printf("out of mem: %d < %d!\n", pool->size + size, pool->cap);
130 return 0;
131 }
132
poolFree(void * userData,void * ptr)133 void poolFree( void* userData, void* ptr )
134 {
135 // empty
136 TESS_NOTUSED(userData);
137 TESS_NOTUSED(ptr);
138 }
139
140
141
142
143
144
145
146
147 /**
148 * Returns TRUE if the ring has clockwise winding.
149 *
150 * @return TRUE if clockwise otherwise FALSE.
151 */
152
isRingClockwise(wxPoint2DDouble * pp,int nPointCount)153 bool isRingClockwise(wxPoint2DDouble *pp, int nPointCount)
154
155 {
156 double dfSum = 0.0;
157
158 for( int iVert = 0; iVert < nPointCount-1; iVert++ )
159 {
160 dfSum += pp[iVert].m_x * pp[iVert+1].m_y
161 - pp[iVert].m_y * pp[iVert+1].m_x;
162 }
163
164 dfSum += pp[nPointCount-1].m_x * pp[0].m_y
165 - pp[nPointCount-1].m_y * pp[0].m_x;
166
167 return dfSum < 0.0;
168 }
169
170 //------------------------------------------------------------------------------
171 // Extended_Geometry Implementation
172 //------------------------------------------------------------------------------
Extended_Geometry()173 Extended_Geometry::Extended_Geometry()
174 {
175 vertex_array = NULL;
176 contour_array = NULL;
177 }
178
~Extended_Geometry()179 Extended_Geometry::~Extended_Geometry()
180 {
181 free(vertex_array);
182 free(contour_array);
183 }
184
185
186 //------------------------------------------------------------------------------
187 // PolyTessGeo Implementation
188 //------------------------------------------------------------------------------
PolyTessGeo()189 PolyTessGeo::PolyTessGeo()
190 {
191 m_pxgeom = NULL;
192 m_printStats = false;
193 m_bstripify = false;
194 m_LOD_meters = 0;
195 }
196
197 // Build PolyTessGeo Object from Extended_Geometry
PolyTessGeo(Extended_Geometry * pxGeom)198 PolyTessGeo::PolyTessGeo(Extended_Geometry *pxGeom)
199 {
200 m_ppg_head = NULL;
201 m_bOK = false;
202
203 m_pxgeom = pxGeom;
204 m_printStats = false;
205 m_bstripify = false;
206 m_LOD_meters= 0;
207 }
208
209 // Build PolyTessGeo Object from OGR Polygon
PolyTessGeo(OGRPolygon * poly,bool bSENC_SM,double ref_lat,double ref_lon,double LOD_meters)210 PolyTessGeo::PolyTessGeo(OGRPolygon *poly, bool bSENC_SM, double ref_lat, double ref_lon,
211 double LOD_meters)
212 {
213 m_printStats = false;
214
215 ErrorCode = 0;
216 m_ppg_head = NULL;
217 m_pxgeom = NULL;
218 m_tess_orient = TESS_HORZ;
219
220 m_ref_lat = ref_lat;
221 m_ref_lon = ref_lon;
222 m_LOD_meters = LOD_meters;
223 m_b_senc_sm = bSENC_SM;
224 m_bmerc_transform = false;
225
226 // PolyGeo BBox as lat/lon
227 OGREnvelope Envelope;
228 poly->getEnvelope(&Envelope);
229 xmin = Envelope.MinX;
230 ymin = Envelope.MinY;
231 xmax = Envelope.MaxX;
232 ymax = Envelope.MaxY;
233
234 m_feature_ref_lat = ymin + (ymax - ymin)/2;
235 m_feature_ref_lon = xmin + (xmax - xmin)/2;
236
237 toSM(m_feature_ref_lat, m_feature_ref_lon, ref_lat, ref_lon, &m_feature_easting, &m_feature_northing);
238
239 // Build the array of contour point counts
240 // Get total number of contours
241 m_ncnt = poly->getNumInteriorRings() + 1;
242
243 // build the contour point countarray
244 m_cntr = (int *)malloc(m_ncnt * sizeof(int));
245
246 m_cntr[0] = poly->getExteriorRing()->getNumPoints();
247
248 for(int i = 1; i < m_ncnt ; i++){
249 m_cntr[i] = poly->getInteriorRing( i-1 )->getNumPoints();
250 }
251
252 // Build the point array index table
253
254 OGRPoint p;
255 m_vertexPtrArray = (double **)malloc(m_ncnt * sizeof(double *));
256
257 m_vertexPtrArray[0] = (double *)malloc(m_cntr[0] * sizeof(double) * 2);
258 double *pp = m_vertexPtrArray[0];
259 for(int i = 0 ; i < m_cntr[0] ; i++){
260 poly->getExteriorRing()->getPoint(i, &p);
261
262 // Calculate SM from feature reference point
263 double easting, northing;
264 toSM(p.getY(), p.getX(), m_feature_ref_lat, m_feature_ref_lon, &easting, &northing);
265 *pp++ = easting; // x
266 *pp++ = northing; // y
267 }
268
269 for(int i = 1; i < m_ncnt ; i++){
270 m_vertexPtrArray[i] = (double *)malloc(m_cntr[i] * sizeof(double) * 2);
271 double *pp = m_vertexPtrArray[i];
272 for(int j = 0 ; j < m_cntr[i] ; j++){
273 poly->getInteriorRing(i-1)->getPoint(j, &p);
274
275 double easting, northing;
276 toSM(p.getY(), p.getX(), m_feature_ref_lat, m_feature_ref_lon, &easting, &northing);
277 *pp++ = easting; // x
278 *pp++ = northing; // y
279 }
280 }
281
282
283
284 mx_rate = 1.0;
285 mx_offset = 0.0;
286 my_rate = 1.0;
287 my_offset = 0.0;
288
289 m_bstripify = true;
290
291 earthAxis = CM93_semimajor_axis_meters * mercator_k0;
292
293 m_bcm93 = false;
294
295 BuildTessGLU();
296
297 // Free the working memory
298 for(int i = 0; i < m_ncnt ; i++)
299 free (m_vertexPtrArray[i]);
300
301 free( m_vertexPtrArray );
302 m_vertexPtrArray = nullptr;
303
304 }
305
306
SetExtents(double x_left,double y_bot,double x_right,double y_top)307 void PolyTessGeo::SetExtents(double x_left, double y_bot, double x_right, double y_top)
308 {
309 xmin = x_left;
310 ymin = y_bot;
311 xmax = x_right;
312 ymax = y_top;
313 }
314
315
316
~PolyTessGeo()317 PolyTessGeo::~PolyTessGeo()
318 {
319 delete m_ppg_head;
320 delete m_pxgeom;
321 }
322
BuildDeferredTess(void)323 int PolyTessGeo::BuildDeferredTess(void)
324 {
325 if(m_pxgeom){
326
327 // For cm93
328 m_b_senc_sm = false;
329 m_bmerc_transform = true;
330 m_tess_orient = TESS_HORZ; // prefer horizontal tristrips
331
332 // Get total number of contours
333 m_ncnt = m_pxgeom->n_contours;
334
335 // build the contour point countarray
336 m_cntr = (int *)malloc(m_ncnt * sizeof(int));
337
338 for(int i = 0; i < m_ncnt ; i++){
339 m_cntr[i] = m_pxgeom->contour_array[i];
340 }
341
342
343 // Build the point array index table
344
345 m_vertexPtrArray = (double **)malloc(m_ncnt * sizeof(double *));
346
347 double *vertex_pointer = (double *)m_pxgeom->vertex_array;
348 vertex_pointer += 2; // skip 0
349
350 for(int k = 0 ; k < m_ncnt ; k++){
351 m_vertexPtrArray[k] = vertex_pointer;
352 vertex_pointer += 2 * m_cntr[k];
353 }
354
355 mx_rate = m_pxgeom->x_rate;
356 mx_offset = m_pxgeom->x_offset;
357 my_rate = m_pxgeom->y_rate;
358 my_offset = m_pxgeom->y_offset;
359
360 m_feature_easting = 0;
361 m_feature_northing = 0;
362
363 m_ref_lat = m_pxgeom->ref_lat;
364 m_ref_lon = m_pxgeom->ref_lon;
365
366 // Note Bene... cm93 assumes no flattening constant
367 earthAxis = CM93_semimajor_axis_meters;
368
369
370 m_bcm93 = true;
371
372 // For cm93, we prefer TessGLU, since it produces more TriangleStrips, so renders faster with less storage.
373 int rv = BuildTessGLU();
374
375 // All done with this geometry
376 delete m_pxgeom;
377 m_pxgeom = NULL;
378
379 return rv;
380
381 }
382 else
383 return 0;
384 }
385
386
387
388 void beginCallback(GLenum which, void *polyData);
389 void errorCallback(GLenum errorCode, void *polyData);
390 void endCallback(void *polyData);
391 void vertexCallback(GLvoid *vertex, void *polyData);
392 void combineCallback(GLdouble coords[3],
393 GLdouble *vertex_data[4],
394 GLfloat weight[4], GLdouble **dataOut, void *polyData );
395
396
397 // Build PolyTessGeo Object from OGR Polygon
398 // Using OpenGL/GLU tesselator
BuildTessGLU()399 int PolyTessGeo::BuildTessGLU()
400 {
401
402 unsigned int iir, ip;
403 GLdouble *geoPt;
404
405
406 #if 0
407 // Make a quick sanity check of the polygon coherence
408 bool b_ok = true;
409 OGRLineString *tls = poly->getExteriorRing();
410 if(!tls) {
411 b_ok = false;
412 }
413 else {
414 int tnpta = poly->getExteriorRing()->getNumPoints();
415 if(tnpta < 3 )
416 b_ok = false;
417 }
418
419
420 for( iir=0 ; iir < poly->getNumInteriorRings() ; iir++)
421 {
422 int tnptr = poly->getInteriorRing(iir)->getNumPoints();
423 if( tnptr < 3 )
424 b_ok = false;
425 }
426
427 if( !b_ok )
428 return ERROR_BAD_OGRPOLY;
429
430 #endif
431
432
433
434 // Allocate a work buffer, which will be grown as needed
435 #define NINIT_BUFFER_LEN 10000
436 m_pwork_buf = (GLdouble *)malloc(NINIT_BUFFER_LEN * 2 * sizeof(GLdouble));
437 m_buf_len = NINIT_BUFFER_LEN * 2;
438 m_buf_idx = 0;
439
440 // Create an array to hold pointers to allocated vertices created by "combine" callback,
441 // so that they may be deleted after tesselation.
442 m_pCombineVertexArray = new wxArrayPtrVoid;
443
444 // Create tesselator
445 GLUtessobj = gluNewTess();
446
447 // Register the callbacks
448 gluTessCallback(GLUtessobj, GLU_TESS_BEGIN_DATA, (GLvoid (*) ())&beginCallback);
449 gluTessCallback(GLUtessobj, GLU_TESS_BEGIN_DATA, (GLvoid (*) ())&beginCallback);
450 gluTessCallback(GLUtessobj, GLU_TESS_VERTEX_DATA, (GLvoid (*) ())&vertexCallback);
451 gluTessCallback(GLUtessobj, GLU_TESS_END_DATA, (GLvoid (*) ())&endCallback);
452 gluTessCallback(GLUtessobj, GLU_TESS_COMBINE_DATA, (GLvoid (*) ())&combineCallback);
453 gluTessCallback(GLUtessobj, GLU_TESS_ERROR_DATA, (GLvoid (*) ())&errorCallback);
454
455 gluTessProperty(GLUtessobj, GLU_TESS_WINDING_RULE, GLU_TESS_WINDING_POSITIVE );
456 gluTessNormal(GLUtessobj, 0, 0, 1);
457
458 // Get total number of points(vertices) to build a buffer
459 int npta = 0;
460
461 for( int i=0 ; i < m_ncnt ; i++)
462 npta += m_cntr[i];
463
464
465 geoPt = (GLdouble *)malloc((npta + 6) * 3 * sizeof(GLdouble)); // vertex array
466
467
468
469 // Grow the work buffer if necessary
470
471 if((npta * 4) > m_buf_len)
472 {
473 m_pwork_buf = (GLdouble *)realloc(m_pwork_buf, npta * 4 * sizeof(GLdouble));
474 m_buf_len = npta * 4;
475 }
476
477
478 // Define the polygon
479 gluTessBeginPolygon(GLUtessobj, this);
480
481
482
483 // Check and account for winding direction of ring
484
485
486 double x0, y0, x, y;
487 // OGRPoint p;
488
489 wxPoint2DDouble *pp = (wxPoint2DDouble *)m_vertexPtrArray[0];
490 bool cw = !isRingClockwise(pp, m_cntr[0]);
491
492 //pp++; // skip 0?
493
494 if(cw)
495 {
496 x0 = pp->m_x;
497 y0 = pp->m_y;
498 }
499 else
500 {
501 x0 = pp[m_cntr[0]-1].m_x;
502 y0 = pp[m_cntr[0]-1].m_y;
503 }
504
505
506
507 unsigned int ptValid = m_cntr[0];
508
509 // Transcribe points to vertex array, in proper order with no duplicates
510 // also, accounting for tess_orient
511 GLdouble *ppt = geoPt;
512
513 for(ip = 0 ; ip < (unsigned int)m_cntr[0] ; ip++)
514 {
515 int pidx;
516 if(cw)
517 pidx = m_cntr[0] - ip - 1;
518
519 else
520 pidx = ip;
521
522 x = pp[pidx].m_x;
523 y = pp[pidx].m_y;
524
525 if((fabs(x-x0) > EQUAL_EPS) || (fabs(y-y0) > EQUAL_EPS))
526 {
527 if(m_tess_orient == TESS_VERT)
528 {
529 *ppt++ = x;
530 *ppt++ = y;
531 *ppt++ = 0;
532 }
533 else
534 {
535 *ppt++ = y;
536 *ppt++ = x;
537 *ppt++ = 0;
538 }
539 }
540 else
541 ptValid--;
542
543 x0 = x;
544 y0 = y;
545 }
546
547 // Apply LOD reduction
548
549 if(ptValid > 20 && (m_LOD_meters > .01)) {
550 std::vector<bool> bool_keep(ptValid, false);
551
552 // Keep a few key points
553 bool_keep[0] = true;
554 bool_keep[1] = true;
555 bool_keep[ptValid-1] = true;
556 bool_keep[ptValid-2] = true;
557
558 DouglasPeuckerDI(geoPt, 1, ptValid-2, m_LOD_meters, bool_keep);
559
560 // Create a new buffer
561 double *LOD_result = (double *)malloc((m_cntr[0]) * 3 * sizeof(double));
562 double *plod = LOD_result;
563 int kept_LOD =0;
564
565 for(unsigned int i=0 ; i < ptValid ; i++){
566 if(bool_keep[i]){
567 double x = geoPt[i*3];
568 double y = geoPt[(i*3) + 1];
569 *plod++ = x;
570 *plod++ = y;
571 *plod++ = 0;
572 kept_LOD++;
573 }
574 }
575
576 // Copy the lod points back into the vertex buffer
577 memcpy(geoPt, LOD_result, kept_LOD * 3 * sizeof(double));
578
579 free(LOD_result);
580 ptValid = kept_LOD;
581 }
582
583 // Declare the gluContour and copy the points
584 gluTessBeginContour(GLUtessobj);
585
586 double *DPrun = geoPt;
587 for(ip = 0 ; ip < ptValid ; ip++) {
588 gluTessVertex( GLUtessobj, DPrun, DPrun ) ;
589 DPrun += 3;
590 }
591
592 gluTessEndContour(GLUtessobj);
593
594 // Now the interior contours
595 #if 1
596 int gpIndex = m_cntr[0];
597
598 for(iir=0; iir < (unsigned int)m_ncnt-1; iir++){
599
600 wxPoint2DDouble *pp = (wxPoint2DDouble *)m_vertexPtrArray[iir + 1];
601
602 int npti = m_cntr[iir+1];
603 ptValid = npti;
604
605 double *ppt = &geoPt[gpIndex * 3]; // next available location in geoPT
606 double *DPStart = ppt;
607
608 // Check and account for winding direction of ring
609 bool cw = isRingClockwise(pp, m_cntr[iir+1]);
610
611 if(cw)
612 {
613 x0 = pp[0].m_x;
614 y0 = pp[0].m_y;
615 }
616 else
617 {
618 x0 = pp[npti-1].m_x;
619 y0 = pp[npti-1].m_y;
620 }
621
622 // Transcribe points to vertex array, in proper order with no duplicates
623 // also, accounting for tess_orient
624 for(int ip = 0 ; ip < npti ; ip++)
625 {
626 int pidx;
627 if(cw) // interior contours must be cw
628 pidx = npti - ip - 1;
629 else
630 pidx = ip;
631
632
633 x = pp[pidx].m_x;
634 y = pp[pidx].m_y;
635
636 if((fabs(x-x0) > EQUAL_EPS) || (fabs(y-y0) > EQUAL_EPS))
637 {
638 if(m_tess_orient == TESS_VERT)
639 {
640 *ppt++ = x;
641 *ppt++ = y;
642 *ppt++ = 0;
643 }
644 else
645 {
646 *ppt++ = y;
647 *ppt++ = x;
648 *ppt++ = 0;
649 }
650 }
651 else
652 ptValid--;
653
654 x0 = x;
655 y0 = y;
656
657 }
658
659 // Declare the gluContour and reference the points
660 gluTessBeginContour(GLUtessobj);
661
662 double *DPruni = DPStart;
663 for(ip = 0 ; ip < ptValid ; ip++)
664 {
665 gluTessVertex( GLUtessobj, DPruni, DPruni ) ;
666 DPruni += 3;
667 }
668
669 gluTessEndContour(GLUtessobj);
670
671 gpIndex += m_cntr[iir+1];
672 }
673 #endif
674
675
676
677
678
679 m_pTPG_Last = NULL;
680 m_pTPG_Head = NULL;
681 m_nvmax = 0;
682
683 // Ready to kick off the tesselator
684 gluTessEndPolygon(GLUtessobj); // here it goes
685
686 m_nvertex_max = m_nvmax; // record largest vertex count, updates in callback
687
688
689 // Tesselation all done, so...
690
691 // Create the data structures
692
693 m_ppg_head = new PolyTriGroup;
694 m_ppg_head->m_bSMSENC = m_b_senc_sm;
695
696 m_ppg_head->nContours = m_ncnt;
697
698 m_ppg_head->pn_vertex = m_cntr; // pointer to array of poly vertex counts
699 m_ppg_head->data_type = DATA_TYPE_DOUBLE;
700
701
702 // Transcribe the raw geometry buffer
703 // Converting to float as we go, and
704 // allowing for tess_orient
705 // Also, convert to SM if requested
706
707 // Recalculate the size of the geometry buffer
708 unsigned int nptfinal = m_cntr[0] + 2;
709 // for(int i=0 ; i < nint ; i++)
710 // nptfinal += cntr[i+1] + 2;
711
712 // No longer need the full geometry in the SENC,
713 nptfinal = 1;
714
715 m_nwkb = (nptfinal + 1) * 2 * sizeof(float);
716 m_ppg_head->pgroup_geom = (float *)calloc(sizeof(float), (nptfinal + 1) * 2);
717 float *vro = m_ppg_head->pgroup_geom;
718 ppt = geoPt;
719 float tx,ty;
720
721 for(ip = 0 ; ip < nptfinal ; ip++)
722 {
723 if(TESS_HORZ == m_tess_orient)
724 {
725 ty = *ppt++;
726 tx = *ppt++;
727 }
728 else
729 {
730 tx = *ppt++;
731 ty = *ppt++;
732 }
733
734 if(m_b_senc_sm)
735 {
736 // Calculate SM from chart common reference point
737 double easting, northing;
738 toSM(ty, tx, m_ref_lat, m_ref_lon, &easting, &northing);
739 *vro++ = easting; // x
740 *vro++ = northing; // y
741 }
742 else
743 {
744 *vro++ = tx; // lon
745 *vro++ = ty; // lat
746 }
747
748 ppt++; // skip z
749 }
750
751 m_ppg_head->tri_prim_head = m_pTPG_Head; // head of linked list of TriPrims
752
753
754 // Convert the Triangle vertex arrays into a single memory allocation of floats
755 // to reduce SENC size and enable efficient access later
756
757 // First calculate the total byte size
758 int total_byte_size = 2 * sizeof(float);
759 TriPrim *p_tp = m_ppg_head->tri_prim_head;
760 while( p_tp ) {
761 total_byte_size += p_tp->nVert * 2 * sizeof(float);
762 p_tp = p_tp->p_next; // pick up the next in chain
763 }
764
765 float *vbuf = (float *)malloc(total_byte_size);
766 p_tp = m_ppg_head->tri_prim_head;
767 float *p_run = vbuf;
768 while( p_tp ) {
769 float *pfbuf = p_run;
770 GLdouble *pdouble_buf = (GLdouble *)p_tp->p_vertex;
771
772 for( int i=0 ; i < p_tp->nVert * 2 ; ++i){
773 float x = (float)( *((GLdouble *)pdouble_buf) );
774 pdouble_buf++;
775 *p_run++ = x;
776 }
777
778 free(p_tp->p_vertex);
779 p_tp->p_vertex = (double *)pfbuf;
780 p_tp = p_tp->p_next; // pick up the next in chain
781 }
782 m_ppg_head->bsingle_alloc = true;
783 m_ppg_head->single_buffer = (unsigned char *)vbuf;
784 m_ppg_head->single_buffer_size = total_byte_size;
785 m_ppg_head->data_type = DATA_TYPE_FLOAT;
786
787
788
789
790
791
792 gluDeleteTess(GLUtessobj);
793
794 free( m_pwork_buf );
795 m_pwork_buf = NULL;
796
797 free (geoPt);
798
799 // Free up any "Combine" vertices created
800 for(unsigned int i = 0; i < m_pCombineVertexArray->GetCount() ; i++)
801 free (m_pCombineVertexArray->Item(i));
802 delete m_pCombineVertexArray;
803
804 m_bOK = true;
805
806
807 return 0;
808 }
809
810
811
BuildTess(void)812 int PolyTessGeo::BuildTess(void)
813 {
814 // Setup the tesselator
815
816 TESSalloc ma;
817 TESStesselator* tess = 0;
818 const int nvp = 3;
819 // unsigned char* vflags = 0;
820 #ifdef USE_POOL
821 struct MemPool pool;
822 unsigned char mem[4024*1024];
823 // int nvflags = 0;
824 #else
825 int allocated = 0;
826 #endif
827
828 #ifdef USE_POOL
829
830 pool.size = 0;
831 pool.cap = sizeof(mem);
832 pool.buf = mem;
833 memset(&ma, 0, sizeof(ma));
834 ma.memalloc = poolAlloc;
835 ma.memfree = poolFree;
836 ma.userData = (void*)&pool;
837 ma.extraVertices = 256; // realloc not provided, allow 256 extra vertices.
838
839 #else
840
841 memset(&ma, 0, sizeof(ma));
842 ma.memalloc = stdAlloc;
843 ma.memfree = stdFree;
844 ma.userData = (void*)&allocated;
845 ma.extraVertices = 256; // realloc not provided, allow 256 extra vertices.
846 #endif
847
848 tess = tessNewTess(&ma);
849 if (!tess)
850 return -1;
851
852 //tessSetOption(tess, TESS_CONSTRAINED_DELAUNAY_TRIANGULATION, 1);
853
854
855 // Create the contour vertex arrays
856
857 int iir, ip;
858
859 // Get max number number of points(vertices) in any contour
860 int npta = m_cntr[0];
861 npta += 2; // fluff
862
863 for( iir=0 ; iir < m_ncnt-1 ; iir++){
864 int nptr = m_cntr[iir+1];
865 npta = wxMax(npta, nptr);
866 }
867
868
869 TESSreal *geoPt = (TESSreal *)malloc((npta) * 2 * sizeof(TESSreal)); // tess input vertex array
870 TESSreal *ppt = geoPt;
871
872
873 // Create input structures
874
875 // Exterior Ring
876 int npte = m_cntr[0];
877 unsigned int ptValid = npte;
878
879 // Check and account for winding direction of ring
880
881
882 double x0, y0, x, y;
883 OGRPoint p;
884
885 wxPoint2DDouble *pp = (wxPoint2DDouble *)m_vertexPtrArray[0];
886 bool cw = isRingClockwise(pp, m_cntr[0]);
887
888 //pp++; // skip 0?
889
890 if(cw)
891 {
892 x0 = pp->m_x;
893 y0 = pp->m_y;
894 }
895 else
896 {
897 x0 = pp[npte-1].m_x;
898 y0 = pp[npte-1].m_y;
899 }
900
901
902
903 // Transcribe points to vertex array, in proper order with no duplicates
904 // also, accounting for tess_orient
905 for(ip = 0 ; ip < npte ; ip++)
906 {
907 int pidx;
908 if(cw)
909 pidx = npte - ip - 1;
910
911 else
912 pidx = ip;
913
914 x = pp[pidx].m_x;
915 y = pp[pidx].m_y;
916
917 if((fabs(x-x0) > EQUAL_EPS) || (fabs(y-y0) > EQUAL_EPS))
918 {
919 if(m_tess_orient == TESS_VERT)
920 {
921 *ppt++ = x;
922 *ppt++ = y;
923 }
924 else
925 {
926 *ppt++ = y;
927 *ppt++ = x;
928 }
929 }
930 else
931 ptValid--;
932
933 x0 = x;
934 y0 = y;
935 }
936
937
938 // Apply LOD reduction
939 int beforeLOD = ptValid;
940 int afterLOD = beforeLOD;
941
942 std::vector<bool> bool_keep;
943 if(ptValid > 5 && (m_LOD_meters > .01)){
944
945 for(unsigned int i = 0 ; i < ptValid ; i++)
946 bool_keep.push_back(false);
947
948 // Keep a few key points
949 bool_keep[0] = true;
950 bool_keep[1] = true;
951 bool_keep[ptValid-1] = true;
952 bool_keep[ptValid-2] = true;
953
954 DouglasPeuckerFI(geoPt, 1, ptValid-2, m_LOD_meters, bool_keep);
955
956 // Create a new buffer
957 float *LOD_result = (float *)malloc((npte) * 2 * sizeof(float));
958 float *plod = LOD_result;
959 int kept_LOD =0;
960
961 for(unsigned int i=0 ; i < ptValid ; i++){
962 if(bool_keep[i]){
963 float x = geoPt[i*2];
964 float y = geoPt[(i*2) + 1];
965 *plod++ = x;
966 *plod++ = y;
967 kept_LOD++;
968 }
969 }
970
971 beforeLOD = ptValid;
972 afterLOD = kept_LOD;
973
974 tessAddContour(tess, 2, LOD_result, sizeof(float)*2, kept_LOD);
975
976 free(LOD_result);
977 }
978 else {
979 tessAddContour(tess, 2, geoPt, sizeof(float)*2, ptValid);
980 }
981
982
983
984 #if 1
985 // Now the interior contours
986 for(iir=0; iir < m_ncnt-1; iir++)
987 {
988 ppt = geoPt;
989 wxPoint2DDouble *pp = (wxPoint2DDouble *)m_vertexPtrArray[iir + 1];
990
991 int npti = m_cntr[iir+1];
992 ptValid = npti;
993
994 // Check and account for winding direction of ring
995 bool cw = isRingClockwise(pp, m_cntr[iir+1]);
996
997 if(!cw)
998 {
999 x0 = pp[0].m_x;
1000 y0 = pp[0].m_y;
1001 }
1002 else
1003 {
1004 x0 = pp[npti-1].m_x;
1005 y0 = pp[npti-1].m_y;
1006 }
1007
1008 // Transcribe points to vertex array, in proper order with no duplicates
1009 // also, accounting for tess_orient
1010 for(int ip = 0 ; ip < npti ; ip++)
1011 {
1012 OGRPoint p;
1013 int pidx;
1014 if(!cw) // interior contours must be cw
1015 pidx = npti - ip - 1;
1016 else
1017 pidx = ip;
1018
1019
1020 x = pp[pidx].m_x;
1021 y = pp[pidx].m_y;
1022
1023 if((fabs(x-x0) > EQUAL_EPS) || (fabs(y-y0) > EQUAL_EPS))
1024 {
1025 if(m_tess_orient == TESS_VERT)
1026 {
1027 *ppt++ = x;
1028 *ppt++ = y;
1029 }
1030 else
1031 {
1032 *ppt++ = y;
1033 *ppt++ = x;
1034 }
1035 }
1036 else
1037 ptValid--;
1038
1039 x0 = x;
1040 y0 = y;
1041
1042 }
1043
1044 tessAddContour(tess, 2, geoPt, sizeof(float)*2, ptValid);
1045
1046 }
1047 #endif
1048
1049
1050 // Ready to kick off the tesselator
1051
1052 TriPrim *pTPG_Last = NULL;
1053 TriPrim *pTPG_Head = NULL;
1054
1055 //s_nvmax = 0;
1056
1057 //OCPNStopWatchTess tt0;
1058
1059 if (!tessTesselate(tess, TESS_WINDING_POSITIVE, TESS_POLYGONS, nvp, 2, 0))
1060 return -1;
1061 double tessTime = 0; //tt0.GetTime();
1062
1063
1064 // Tesselation all done, so...
1065 // Create the data structures
1066
1067 // Make a linked list of TriPrims from tess output arrays
1068
1069 const TESSreal* verts = tessGetVertices(tess);
1070 const int* vinds = tessGetVertexIndices(tess);
1071 const int* elems = tessGetElements(tess);
1072 const int nverts = tessGetVertexCount(tess);
1073 int nelems = tessGetElementCount(tess);
1074
1075
1076 bool skip = false;
1077 //skip = true;
1078
1079 double stripTime = 0;
1080 //tt0.Reset();
1081
1082 int bytes_needed_vbo = 0;
1083 float *vbo = 0;
1084
1085 if(m_bstripify && nelems && !skip){
1086
1087
1088 STRIPERCREATE sc;
1089 sc.DFaces = (udword *)elems;
1090 sc.NbFaces = nelems;
1091 sc.AskForWords = false;
1092 sc.ConnectAllStrips = false;
1093 sc.OneSided = false;
1094 sc.SGIAlgorithm = false;
1095
1096
1097 Striper Strip;
1098 Strip.Init(sc);
1099
1100 STRIPERRESULT sr;
1101 Strip.Compute(sr);
1102
1103 stripTime = 0; //tt0.GetTime();
1104
1105 /*
1106 fprintf(stdout, "Number of strips: %d\n", sr.NbStrips);
1107 fprintf(stdout, "Number of points: %d\n", nelems);
1108 uword* Refs = (uword*)sr.StripRuns;
1109 for(udword i=0;i<sr.NbStrips;i++)
1110 {
1111 fprintf(stdout, "Strip %d: ", i);
1112 udword NbRefs = sr.StripLengths[i];
1113 for(udword j=0;j<NbRefs;j++)
1114 {
1115 fprintf(stdout, "%d ", *Refs++);
1116 }
1117 fprintf(stdout, "\n");
1118 }
1119 */
1120 // calculate and allocate the final (float) VBO-like buffer for this entire feature
1121
1122 int *Refs = (int *)sr.StripRuns;
1123 for(unsigned int i=0;i<sr.NbStrips;i++){
1124 unsigned int NbRefs = sr.StripLengths[i]; // vertices per strip
1125 bytes_needed_vbo += NbRefs * 2 * sizeof(float);
1126 }
1127
1128 vbo = (float *)malloc(bytes_needed_vbo);
1129 float *vbo_run = vbo;
1130
1131 for(unsigned int i=0;i<sr.NbStrips;i++){
1132 unsigned int NbRefs = sr.StripLengths[i]; // vertices per strip
1133
1134 if(NbRefs >= 3){ // this is a valid primitive
1135
1136 TriPrim *pTPG = new TriPrim;
1137 if(NULL == pTPG_Last){
1138 pTPG_Head = pTPG;
1139 pTPG_Last = pTPG;
1140 }
1141 else{
1142 pTPG_Last->p_next = pTPG;
1143 pTPG_Last = pTPG;
1144 }
1145
1146 pTPG->p_next = NULL;
1147 pTPG->nVert = NbRefs;
1148
1149 // pTPG->p_vertex = (double *)malloc(NbRefs * 2 * sizeof(double));
1150 // GLdouble *pdd = (GLdouble*)pTPG->p_vertex;
1151
1152 pTPG->p_vertex = (double *)vbo_run;
1153
1154 // Preset LLBox bounding box limits
1155 double sxmax = -1e8;
1156 double sxmin = 1e8;
1157 double symax = -1e8;
1158 double symin = 1e8;
1159
1160 if(NbRefs >3)
1161 pTPG->type = GL_TRIANGLE_STRIP;
1162 else
1163 pTPG->type = GL_TRIANGLES;
1164
1165 // Add the first two points
1166 int vindex[3];
1167 vindex[0] = *Refs++;
1168 vindex[1] = *Refs++;
1169 unsigned int np = 2; // for the first triangle
1170
1171 for(unsigned int i = 2; i < NbRefs ; i++){
1172 vindex[np] = *Refs++;
1173 np++;
1174
1175 for (unsigned int j = 0; j < np; ++j){
1176 double yd = verts[vindex[j]*2];
1177 double xd = verts[vindex[j]*2+1];
1178
1179 // Calculate LLBox bounding box for each Tri-prim
1180 if(m_bmerc_transform){
1181 double valx = ( xd * mx_rate ) + mx_offset;
1182 double valy = ( yd * my_rate ) + my_offset;
1183
1184 // quickly convert to lat/lon
1185 double lat = ( 2.0 * atan ( exp ( valy/CM93_semimajor_axis_meters ) ) - PI/2. ) / DEGREE;
1186 double lon= ( valx / ( DEGREE * CM93_semimajor_axis_meters ) );
1187
1188
1189 sxmax = wxMax(lon, sxmax);
1190 sxmin = wxMin(lon, sxmin);
1191 symax = wxMax(lat, symax);
1192 symin = wxMin(lat, symin);
1193 }
1194 else{
1195 sxmax = wxMax(xd, sxmax);
1196 sxmin = wxMin(xd, sxmin);
1197 symax = wxMax(yd, symax);
1198 symin = wxMin(yd, symin);
1199 }
1200
1201 // Append this point to TriPrim vbo
1202
1203 *vbo_run++ = (xd) + m_feature_easting; // adjust to chart ref coordinates
1204 *vbo_run++ = (yd) + m_feature_northing;
1205 } // For
1206
1207 // Compute the final LLbbox for this TriPrim chain
1208 double minlat, minlon, maxlat, maxlon;
1209 fromSM(sxmin, symin, m_feature_ref_lat, m_feature_ref_lon, &minlat, &minlon);
1210 fromSM(sxmax, symax, m_feature_ref_lat, m_feature_ref_lon, &maxlat, &maxlon);
1211
1212 pTPG->tri_box.Set(minlat, minlon, maxlat, maxlon);
1213
1214 // set for next single point
1215 np = 0;
1216
1217 }
1218 }
1219 else
1220 Refs += sr.StripLengths[i];
1221
1222 } // for strips
1223 }
1224 else{ // not stripified
1225
1226 m_nvertex_max = nverts; // record largest vertex count
1227
1228 bytes_needed_vbo = nelems * nvp * 2 * sizeof(float);
1229 vbo = (float *)malloc(bytes_needed_vbo);
1230 float *vbo_run = vbo;
1231
1232 for (int i = 0; i < nelems; ++i){
1233 const int* p = &elems[i*nvp];
1234
1235 TriPrim *pTPG = new TriPrim;
1236 if(NULL == pTPG_Last){
1237 pTPG_Head = pTPG;
1238 pTPG_Last = pTPG;
1239 }
1240 else{
1241 pTPG_Last->p_next = pTPG;
1242 pTPG_Last = pTPG;
1243 }
1244
1245 pTPG->p_next = NULL;
1246 pTPG->type = GL_TRIANGLES;
1247 pTPG->nVert = nvp;
1248
1249 // pTPG->p_vertex = (double *)malloc(nvp * 2 * sizeof(double));
1250 // GLdouble *pdd = (GLdouble*)pTPG->p_vertex;
1251
1252 pTPG->p_vertex = (double *)vbo_run;
1253
1254 // Preset LLBox bounding box limits
1255 double sxmax = -1e8;
1256 double sxmin = 1e8;
1257 double symax = -1e8;
1258 double symin = 1e8;
1259
1260 for (size_t j = 0; j < nvp && p[j] != TESS_UNDEF; ++j){
1261 double yd = verts[p[j]*2];
1262 double xd = verts[p[j]*2+1];
1263
1264 // Calculate LLBox bounding box for each Tri-prim
1265 if(m_bmerc_transform){
1266 double valx = ( xd * mx_rate ) + mx_offset;
1267 double valy = ( yd * my_rate ) + my_offset;
1268
1269 // quickly convert to lat/lon
1270 double lat = ( 2.0 * atan ( exp ( valy/CM93_semimajor_axis_meters ) ) - PI/2. ) / DEGREE;
1271 double lon= ( valx / ( DEGREE * CM93_semimajor_axis_meters ) );
1272
1273
1274 sxmax = wxMax(lon, sxmax);
1275 sxmin = wxMin(lon, sxmin);
1276 symax = wxMax(lat, symax);
1277 symin = wxMin(lat, symin);
1278 }
1279 else{
1280 sxmax = wxMax(xd, sxmax);
1281 sxmin = wxMin(xd, sxmin);
1282 symax = wxMax(yd, symax);
1283 symin = wxMin(yd, symin);
1284 }
1285
1286 // Append this point to TriPrim, converting to SM if called for
1287
1288 if(m_b_senc_sm){
1289 double easting, northing;
1290 toSM(yd, xd, m_ref_lat, m_ref_lon, &easting, &northing);
1291 *vbo_run++ = easting;
1292 *vbo_run++ = northing;
1293 }
1294 else{
1295 *vbo_run++ = xd;
1296 *vbo_run++ = yd;
1297 }
1298 }
1299
1300 pTPG->tri_box.Set(symin, sxmin, symax, sxmax);
1301 }
1302 } // stripify
1303
1304 if(m_printStats){
1305 int nTri = 0;
1306 int nStrip = 0;
1307
1308 TriPrim *p_tp = pTPG_Head;
1309 while( p_tp ) {
1310 if(p_tp->type == GL_TRIANGLES)
1311 nTri++;
1312 if(p_tp->type == GL_TRIANGLE_STRIP)
1313 nStrip++;
1314
1315 p_tp = p_tp->p_next; // pick up the next in chain
1316 }
1317
1318 if((nTri + nStrip) > 10000){
1319 printf("LOD: %d/%d\n", afterLOD, beforeLOD);
1320
1321 printf("Tess time(ms): %f\n", tessTime);
1322 printf("Strip time(ms): %f\n", stripTime);
1323
1324 printf("Primitives: Tri: %5d Strip: %5d Total: %5d\n", nTri, nStrip, nTri + nStrip);
1325 printf("\n");
1326 }
1327 }
1328
1329
1330
1331 m_ppg_head = new PolyTriGroup;
1332 m_ppg_head->m_bSMSENC = m_b_senc_sm;
1333
1334 m_ppg_head->nContours = m_ncnt;
1335 m_ppg_head->pn_vertex = m_cntr; // pointer to array of poly vertex counts
1336 m_ppg_head->tri_prim_head = pTPG_Head; // head of linked list of TriPrims
1337 //m_ppg_head->data_type = DATA_TYPE_DOUBLE;
1338
1339
1340 // Transcribe the raw geometry buffer
1341 // Converting to float as we go, and
1342 // allowing for tess_orient
1343 // Also, convert to SM if requested
1344
1345 int nptfinal = npta;
1346
1347 // No longer need the full geometry in the SENC,
1348 nptfinal = 1;
1349
1350 m_nwkb = (nptfinal +1) * 2 * sizeof(float);
1351 m_ppg_head->pgroup_geom = (float *)malloc(m_nwkb);
1352 float *vro = m_ppg_head->pgroup_geom;
1353 ppt = geoPt;
1354 float tx,ty;
1355
1356 for(ip = 0 ; ip < nptfinal ; ip++)
1357 {
1358 if(TESS_HORZ == m_tess_orient)
1359 {
1360 ty = *ppt++;
1361 tx = *ppt++;
1362 }
1363 else
1364 {
1365 tx = *ppt++;
1366 ty = *ppt++;
1367 }
1368
1369 if(m_b_senc_sm)
1370 {
1371 // Calculate SM from chart common reference point
1372 double easting, northing;
1373 toSM(ty, tx, 0/*ref_lat*/, 0/*ref_lon*/, &easting, &northing);
1374 *vro++ = easting; // x
1375 *vro++ = northing; // y
1376 }
1377 else
1378 {
1379 *vro++ = tx; // lon
1380 *vro++ = ty; // lat
1381 }
1382
1383 ppt++; // skip z
1384 }
1385
1386
1387 m_ppg_head->bsingle_alloc = true;
1388 m_ppg_head->single_buffer = (unsigned char *)vbo;
1389 m_ppg_head->single_buffer_size = bytes_needed_vbo;
1390 m_ppg_head->data_type = DATA_TYPE_FLOAT;
1391
1392
1393 free (geoPt);
1394
1395 // All allocated buffers are owned now by the m_ppg_head
1396 // And will be freed on dtor of this object
1397
1398 if (tess) tessDeleteTess(tess);
1399
1400 m_bOK = true;
1401
1402
1403 return 0;
1404
1405
1406 }
1407
1408
1409
1410
1411
1412 // GLU tesselation support functions
beginCallback(GLenum which,void * polyData)1413 void beginCallback(GLenum which, void *polyData)
1414 {
1415 PolyTessGeo *pThis = (PolyTessGeo *)polyData;
1416
1417 pThis->m_buf_idx = 0;
1418 pThis->m_nvcall = 0;
1419 pThis->m_gltri_type = which;
1420 }
1421
errorCallback(GLenum errorCode,void * polyData)1422 void errorCallback(GLenum errorCode, void *polyData )
1423 {
1424 const GLubyte *estring;
1425
1426 estring = gluErrorString(errorCode);
1427 printf("Tessellation Error: %s\n", estring);
1428 exit(0);
1429 }
1430
1431
endCallback(void * polyData)1432 void endCallback(void *polyData)
1433 {
1434 PolyTessGeo *pThis = (PolyTessGeo *)polyData;
1435
1436 // Create a TriPrim
1437 if(pThis->m_nvcall > pThis->m_nvmax) // keep track of largest number of triangle vertices
1438 pThis->m_nvmax = pThis->m_nvcall;
1439
1440 switch(pThis->m_gltri_type)
1441 {
1442 case GL_TRIANGLE_FAN:
1443 case GL_TRIANGLE_STRIP:
1444 case GL_TRIANGLES:
1445 {
1446 TriPrim *pTPG = new TriPrim;
1447 if(NULL == pThis->m_pTPG_Last)
1448 {
1449 pThis->m_pTPG_Head = pTPG;
1450 pThis->m_pTPG_Last = pTPG;
1451 }
1452 else
1453 {
1454 pThis->m_pTPG_Last->p_next = pTPG;
1455 pThis->m_pTPG_Last = pTPG;
1456 }
1457
1458 pTPG->p_next = NULL;
1459 pTPG->type = pThis->m_gltri_type;
1460 pTPG->nVert = pThis->m_nvcall;
1461
1462 // Calculate bounding box
1463 double sxmax = -1e8; // this poly BBox
1464 double sxmin = 1e8;
1465 double symax = -1e8;
1466 double symin = 1e8;
1467
1468 GLdouble *pvr = pThis->m_pwork_buf;
1469 for(int iv=0 ; iv < pThis->m_nvcall ; iv++)
1470 {
1471 GLdouble xd, yd;
1472 xd = *pvr++;
1473 yd = *pvr++;
1474
1475 if(pThis->m_bcm93)
1476 {
1477 // cm93 hits here
1478 double valx = ( xd * pThis->mx_rate ) + pThis->mx_offset + pThis->m_feature_easting;
1479 double valy = ( yd * pThis->my_rate ) + pThis->my_offset + pThis->m_feature_northing;
1480
1481 // Convert to lat/lon
1482 double lat = ( 2.0 * atan ( exp ( valy/CM93_semimajor_axis_meters ) ) - PI/2. ) / DEGREE;
1483 double lon = ( valx / ( DEGREE * CM93_semimajor_axis_meters ) );
1484
1485 sxmax = wxMax(lon, sxmax);
1486 sxmin = wxMin(lon, sxmin);
1487 symax = wxMax(lat, symax);
1488 symin = wxMin(lat, symin);
1489 }
1490 else
1491 {
1492 // ENC hits here, values are SM measures from feature reference point
1493 double valx = ( xd * pThis->mx_rate ) + pThis->mx_offset + pThis->m_feature_easting;
1494 double valy = ( yd * pThis->my_rate ) + pThis->my_offset + pThis->m_feature_northing;
1495
1496 sxmax = wxMax(valx, sxmax);
1497 sxmin = wxMin(valx, sxmin);
1498 symax = wxMax(valy, symax);
1499 symin = wxMin(valy, symin);
1500 }
1501 }
1502
1503
1504 // Compute the final LLbbox for this TriPrim chain
1505 if(pThis->m_bcm93)
1506 pTPG->tri_box.Set(symin, sxmin, symax, sxmax);
1507 else{
1508 double minlat, minlon, maxlat, maxlon;
1509 fromSM(sxmin, symin, pThis->m_ref_lat, pThis->m_ref_lon, &minlat, &minlon);
1510 fromSM(sxmax, symax, pThis->m_ref_lat, pThis->m_ref_lon, &maxlat, &maxlon);
1511 pTPG->tri_box.Set(minlat, minlon, maxlat, maxlon);
1512 }
1513
1514
1515 // Transcribe this geometry to TriPrim,
1516 {
1517 GLdouble *pds = pThis->m_pwork_buf;
1518 pTPG->p_vertex = (double *)malloc(pThis->m_nvcall * 2 * sizeof(double));
1519 GLdouble *pdd = (GLdouble*)pTPG->p_vertex;
1520
1521 for(int ip = 0 ; ip < pThis->m_nvcall ; ip++)
1522 {
1523 double dlon = *pds++;
1524 double dlat = *pds++;
1525 *pdd++ = dlon + pThis->m_feature_easting; // adjust to feature ref coordinates
1526 *pdd++ = dlat + pThis->m_feature_northing;
1527 }
1528 }
1529
1530 break;
1531 }
1532 default:
1533 {
1534 // sprintf(buf, "....begin Callback unknown\n");
1535 break;
1536 }
1537 }
1538 }
1539
vertexCallback(GLvoid * vertex,void * polyData)1540 void vertexCallback(GLvoid *vertex, void *polyData)
1541 {
1542 PolyTessGeo *pThis = (PolyTessGeo *)polyData;
1543
1544 GLdouble *pointer;
1545
1546 pointer = (GLdouble *) vertex;
1547
1548 if(pThis->m_buf_idx > pThis->m_buf_len - 4)
1549 {
1550 int new_buf_len = pThis->m_buf_len + 100;
1551 GLdouble * tmp = pThis->m_pwork_buf;
1552
1553 pThis->m_pwork_buf = (GLdouble *)realloc(pThis->m_pwork_buf, new_buf_len * sizeof(GLdouble));
1554 if (NULL == pThis->m_pwork_buf)
1555 {
1556 free(tmp);
1557 tmp = NULL;
1558 }
1559 else
1560 pThis->m_buf_len = new_buf_len;
1561 }
1562
1563 if(pThis->m_tess_orient == TESS_VERT)
1564 {
1565 pThis->m_pwork_buf[pThis->m_buf_idx++] = pointer[0];
1566 pThis->m_pwork_buf[pThis->m_buf_idx++] = pointer[1];
1567 }
1568 else
1569 {
1570 pThis->m_pwork_buf[pThis->m_buf_idx++] = pointer[1];
1571 pThis->m_pwork_buf[pThis->m_buf_idx++] = pointer[0];
1572 }
1573
1574
1575 pThis->m_nvcall++;
1576
1577 }
1578
1579 /* combineCallback is used to create a new vertex when edges
1580 * intersect.
1581 */
combineCallback(GLdouble coords[3],GLdouble * vertex_data[4],GLfloat weight[4],GLdouble ** dataOut,void * polyData)1582 void combineCallback(GLdouble coords[3],
1583 GLdouble *vertex_data[4],
1584 GLfloat weight[4], GLdouble **dataOut, void *polyData )
1585 {
1586 PolyTessGeo *pThis = (PolyTessGeo *)polyData;
1587
1588 GLdouble *vertex = (GLdouble *)malloc(6 * sizeof(GLdouble));
1589
1590 vertex[0] = coords[0];
1591 vertex[1] = coords[1];
1592 vertex[2] = coords[2];
1593 vertex[3] = vertex[4] = vertex[5] = 0. ; //01/13/05 bugfix
1594
1595 *dataOut = vertex;
1596
1597 pThis->m_pCombineVertexArray->Add(vertex);
1598 }
1599
1600
1601
1602 wxStopWatch *s_stwt;
1603
1604
1605 //------------------------------------------------------------------------------
1606 // PolyTriGroup Implementation
1607 //------------------------------------------------------------------------------
PolyTriGroup()1608 PolyTriGroup::PolyTriGroup()
1609 {
1610 pn_vertex = NULL; // pointer to array of poly vertex counts
1611 pgroup_geom = NULL; // pointer to Raw geometry, used for contour line drawing
1612 tri_prim_head = NULL; // head of linked list of TriPrims
1613 m_bSMSENC = false;
1614 bsingle_alloc = false;
1615 single_buffer = NULL;
1616 single_buffer_size = 0;
1617 data_type = DATA_TYPE_DOUBLE;
1618 sfactor = 1.0;
1619
1620
1621 }
1622
~PolyTriGroup()1623 PolyTriGroup::~PolyTriGroup()
1624 {
1625 free(pn_vertex);
1626 free(pgroup_geom);
1627 //Walk the list of TriPrims, deleting as we go
1628 TriPrim *tp_next;
1629 TriPrim *tp = tri_prim_head;
1630
1631 if(bsingle_alloc){
1632 free(single_buffer);
1633 while(tp) {
1634 tp_next = tp->p_next;
1635 delete tp;
1636 tp = tp_next;
1637 }
1638 }
1639 else {
1640 while(tp) {
1641 tp_next = tp->p_next;
1642 tp->FreeMem();
1643 delete tp;
1644 tp = tp_next;
1645 }
1646 }
1647 }
1648
1649 //------------------------------------------------------------------------------
1650 // TriPrim Implementation
1651 //------------------------------------------------------------------------------
TriPrim()1652 TriPrim::TriPrim()
1653 {
1654 }
1655
~TriPrim()1656 TriPrim::~TriPrim()
1657 {
1658 }
1659
FreeMem()1660 void TriPrim::FreeMem()
1661 {
1662 free(p_vertex);
1663 }
1664
1665
1666
1667
1668
1669
1670