1 : /******************************************************************************
2 : * $Id: contour.cpp 21236 2010-12-11 17:28:32Z rouault $
3 : *
4 : * Project: Contour Generation
5 : * Purpose: Core algorithm implementation for contour line generation.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
10 : * Copyright (c) 2003, Applied Coherent Technology Corporation, www.actgate.com
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "gdal_priv.h"
32 : #include "gdal_alg.h"
33 : #include "ogr_api.h"
34 :
35 : CPL_CVSID("$Id: contour.cpp 21236 2010-12-11 17:28:32Z rouault $");
36 :
37 : #ifdef OGR_ENABLED
38 :
39 : // The amount of a contour interval that pixels should be fudged by if they
40 : // match a contour level exactly.
41 :
42 : #define FUDGE_EXACT 0.001
43 :
44 : // The amount of a pixel that line ends need to be within to be considered to
45 : // match for joining purposes.
46 :
47 : #define JOIN_DIST 0.0001
48 :
49 : /************************************************************************/
50 : /* GDALContourItem */
51 : /************************************************************************/
52 : class GDALContourItem
53 : {
54 : public:
55 : int bRecentlyAccessed;
56 : double dfLevel;
57 :
58 : int nPoints;
59 : int nMaxPoints;
60 : double *padfX;
61 : double *padfY;
62 :
63 : int bLeftIsHigh;
64 :
65 : double dfTailX;
66 :
67 : GDALContourItem( double dfLevel );
68 : ~GDALContourItem();
69 :
70 : int AddSegment( double dfXStart, double dfYStart,
71 : double dfXEnd, double dfYEnd, int bLeftHigh );
72 : void MakeRoomFor( int );
73 : int Merge( GDALContourItem * );
74 : void PrepareEjection();
75 : };
76 :
77 : /************************************************************************/
78 : /* GDALContourLevel */
79 : /************************************************************************/
80 : class GDALContourLevel
81 : {
82 : double dfLevel;
83 :
84 : int nEntryMax;
85 : int nEntryCount;
86 : GDALContourItem **papoEntries;
87 :
88 : public:
89 : GDALContourLevel( double );
90 : ~GDALContourLevel();
91 :
92 287936 : double GetLevel() { return dfLevel; }
93 200244 : int GetContourCount() { return nEntryCount; }
94 164845 : GDALContourItem *GetContour( int i) { return papoEntries[i]; }
95 : void AdjustContour( int );
96 : void RemoveContour( int );
97 : int FindContour( double dfX, double dfY );
98 : int InsertContour( GDALContourItem * );
99 : };
100 :
101 : /************************************************************************/
102 : /* GDALContourGenerator */
103 : /************************************************************************/
104 : class GDALContourGenerator
105 : {
106 : int nWidth;
107 : int nHeight;
108 : int iLine;
109 :
110 : double *padfLastLine;
111 : double *padfThisLine;
112 :
113 : int nLevelMax;
114 : int nLevelCount;
115 : GDALContourLevel **papoLevels;
116 :
117 : int iLastLevel;
118 :
119 : int bNoDataActive;
120 : double dfNoDataValue;
121 :
122 : int bFixedLevels;
123 : double dfContourInterval;
124 : double dfContourOffset;
125 :
126 : CPLErr AddSegment( double dfLevel,
127 : double dfXStart, double dfYStart,
128 : double dfXEnd, double dfYEnd, int bLeftHigh );
129 :
130 : CPLErr ProcessPixel( int iPixel );
131 : CPLErr ProcessRect( double, double, double,
132 : double, double, double,
133 : double, double, double,
134 : double, double, double );
135 :
136 : void Intersect( double, double, double,
137 : double, double, double,
138 : double, double, int *, double *, double * );
139 :
140 : GDALContourLevel *FindLevel( double dfLevel );
141 :
142 : public:
143 : GDALContourWriter pfnWriter;
144 : void *pWriterCBData;
145 :
146 : GDALContourGenerator( int nWidth, int nHeight,
147 : GDALContourWriter pfnWriter, void *pWriterCBData );
148 : ~GDALContourGenerator();
149 :
150 : void SetNoData( double dfNoDataValue );
151 6 : void SetContourLevels( double dfContourInterval,
152 : double dfContourOffset = 0.0 )
153 6 : { this->dfContourInterval = dfContourInterval;
154 6 : this->dfContourOffset = dfContourOffset; }
155 :
156 : void SetFixedLevels( int, double * );
157 : CPLErr FeedLine( double *padfScanline );
158 : CPLErr EjectContours( int bOnlyUnused = FALSE );
159 :
160 : };
161 :
162 : /************************************************************************/
163 : /* GDAL_CG_Create() */
164 : /************************************************************************/
165 :
166 : GDALContourGeneratorH
167 0 : GDAL_CG_Create( int nWidth, int nHeight, int bNoDataSet, double dfNoDataValue,
168 : double dfContourInterval, double dfContourBase,
169 : GDALContourWriter pfnWriter, void *pCBData )
170 :
171 : {
172 : GDALContourGenerator *poCG = new GDALContourGenerator( nWidth, nHeight,
173 0 : pfnWriter, pCBData );
174 :
175 0 : if( bNoDataSet )
176 0 : poCG->SetNoData( dfNoDataValue );
177 :
178 0 : poCG->SetContourLevels( dfContourInterval, dfContourBase );
179 0 : return (GDALContourGeneratorH) poCG;
180 : }
181 :
182 : /************************************************************************/
183 : /* GDAL_CG_FeedLine() */
184 : /************************************************************************/
185 :
186 0 : CPLErr GDAL_CG_FeedLine( GDALContourGeneratorH hCG, double *padfScanline )
187 :
188 : {
189 0 : VALIDATE_POINTER1( hCG, "GDAL_CG_FeedLine", CE_Failure );
190 :
191 0 : return ((GDALContourGenerator *) hCG)->FeedLine( padfScanline );
192 : }
193 :
194 : /************************************************************************/
195 : /* GDAL_CG_Destroy() */
196 : /************************************************************************/
197 :
198 0 : void GDAL_CG_Destroy( GDALContourGeneratorH hCG )
199 :
200 : {
201 0 : delete ((GDALContourGenerator *) hCG);
202 0 : }
203 :
204 : /************************************************************************/
205 : /* ==================================================================== */
206 : /* GDALContourGenerator */
207 : /* ==================================================================== */
208 : /************************************************************************/
209 :
210 : /************************************************************************/
211 : /* GDALContourGenerator() */
212 : /************************************************************************/
213 :
214 8 : GDALContourGenerator::GDALContourGenerator( int nWidthIn, int nHeightIn,
215 : GDALContourWriter pfnWriterIn,
216 : void *pWriterCBDataIn )
217 : {
218 8 : nWidth = nWidthIn;
219 8 : nHeight = nHeightIn;
220 :
221 8 : padfLastLine = (double *) CPLCalloc(sizeof(double),nWidth);
222 8 : padfThisLine = (double *) CPLCalloc(sizeof(double),nWidth);
223 :
224 8 : pfnWriter = pfnWriterIn;
225 8 : pWriterCBData = pWriterCBDataIn;
226 :
227 8 : iLine = -1;
228 :
229 8 : bNoDataActive = FALSE;
230 8 : dfNoDataValue = -1000000.0;
231 8 : dfContourInterval = 10.0;
232 8 : dfContourOffset = 0.0;
233 :
234 8 : nLevelMax = 0;
235 8 : nLevelCount = 0;
236 8 : papoLevels = NULL;
237 8 : bFixedLevels = FALSE;
238 8 : }
239 :
240 : /************************************************************************/
241 : /* ~GDALContourGenerator() */
242 : /************************************************************************/
243 :
244 8 : GDALContourGenerator::~GDALContourGenerator()
245 :
246 : {
247 : int i;
248 :
249 68 : for( i = 0; i < nLevelCount; i++ )
250 60 : delete papoLevels[i];
251 8 : CPLFree( papoLevels );
252 :
253 8 : CPLFree( padfLastLine );
254 8 : CPLFree( padfThisLine );
255 8 : }
256 :
257 : /************************************************************************/
258 : /* SetFixedLevels() */
259 : /************************************************************************/
260 :
261 2 : void GDALContourGenerator::SetFixedLevels( int nFixedLevelCount,
262 : double *padfFixedLevels )
263 :
264 : {
265 2 : bFixedLevels = TRUE;
266 8 : for( int i = 0; i < nFixedLevelCount; i++ )
267 6 : FindLevel( padfFixedLevels[i] );
268 2 : }
269 :
270 : /************************************************************************/
271 : /* SetNoData() */
272 : /************************************************************************/
273 :
274 3 : void GDALContourGenerator::SetNoData( double dfNewValue )
275 :
276 : {
277 3 : bNoDataActive = TRUE;
278 3 : dfNoDataValue = dfNewValue;
279 3 : }
280 :
281 : /************************************************************************/
282 : /* ProcessPixel() */
283 : /************************************************************************/
284 :
285 159403 : CPLErr GDALContourGenerator::ProcessPixel( int iPixel )
286 :
287 : {
288 : double dfUpLeft, dfUpRight, dfLoLeft, dfLoRight;
289 159403 : int bSubdivide = FALSE;
290 :
291 : /* -------------------------------------------------------------------- */
292 : /* Collect the four corner pixel values. Value left or right */
293 : /* of the scanline are taken from the nearest pixel on the */
294 : /* scanline itself. */
295 : /* -------------------------------------------------------------------- */
296 159403 : dfUpLeft = padfLastLine[MAX(0,iPixel-1)];
297 159403 : dfUpRight = padfLastLine[MIN(nWidth-1,iPixel)];
298 :
299 159403 : dfLoLeft = padfThisLine[MAX(0,iPixel-1)];
300 159403 : dfLoRight = padfThisLine[MIN(nWidth-1,iPixel)];
301 :
302 : /* -------------------------------------------------------------------- */
303 : /* Check if we have any nodata values. */
304 : /* -------------------------------------------------------------------- */
305 159403 : if( bNoDataActive
306 : && ( dfUpLeft == dfNoDataValue
307 : || dfLoLeft == dfNoDataValue
308 : || dfLoRight == dfNoDataValue
309 : || dfUpRight == dfNoDataValue ) )
310 0 : bSubdivide = TRUE;
311 :
312 : /* -------------------------------------------------------------------- */
313 : /* Check if we have any nodata, if so, go to a special case of */
314 : /* code. */
315 : /* -------------------------------------------------------------------- */
316 159403 : if( iPixel > 0 && iPixel < nWidth
317 : && iLine > 0 && iLine < nHeight && !bSubdivide )
318 : {
319 : return ProcessRect( dfUpLeft, iPixel - 0.5, iLine - 0.5,
320 : dfLoLeft, iPixel - 0.5, iLine + 0.5,
321 : dfLoRight, iPixel + 0.5, iLine + 0.5,
322 155217 : dfUpRight, iPixel + 0.5, iLine - 0.5 );
323 : }
324 :
325 : /* -------------------------------------------------------------------- */
326 : /* Prepare subdivisions. */
327 : /* -------------------------------------------------------------------- */
328 4186 : int nGoodCount = 0;
329 4186 : double dfASum = 0.0;
330 4186 : double dfCenter, dfTop=0.0, dfRight=0.0, dfLeft=0.0, dfBottom=0.0;
331 :
332 4186 : if( dfUpLeft != dfNoDataValue )
333 : {
334 4186 : dfASum += dfUpLeft;
335 4186 : nGoodCount++;
336 : }
337 :
338 4186 : if( dfLoLeft != dfNoDataValue )
339 : {
340 4186 : dfASum += dfLoLeft;
341 4186 : nGoodCount++;
342 : }
343 :
344 4186 : if( dfLoRight != dfNoDataValue )
345 : {
346 4186 : dfASum += dfLoRight;
347 4186 : nGoodCount++;
348 : }
349 :
350 4186 : if( dfUpRight != dfNoDataValue )
351 : {
352 4186 : dfASum += dfUpRight;
353 4186 : nGoodCount++;
354 : }
355 :
356 4186 : if( nGoodCount == 0.0 )
357 0 : return CE_None;
358 :
359 4186 : dfCenter = dfASum / nGoodCount;
360 :
361 4186 : if( dfUpLeft != dfNoDataValue )
362 : {
363 4186 : if( dfUpRight != dfNoDataValue )
364 4186 : dfTop = (dfUpLeft + dfUpRight) / 2.0;
365 : else
366 0 : dfTop = dfUpLeft;
367 :
368 4186 : if( dfLoLeft != dfNoDataValue )
369 4186 : dfLeft = (dfUpLeft + dfLoLeft) / 2.0;
370 : else
371 0 : dfLeft = dfUpLeft;
372 : }
373 : else
374 : {
375 0 : dfTop = dfUpRight;
376 0 : dfLeft = dfLoLeft;
377 : }
378 :
379 4186 : if( dfLoRight != dfNoDataValue )
380 : {
381 4186 : if( dfUpRight != dfNoDataValue )
382 4186 : dfRight = (dfLoRight + dfUpRight) / 2.0;
383 : else
384 0 : dfRight = dfLoRight;
385 :
386 4186 : if( dfLoLeft != dfNoDataValue )
387 4186 : dfBottom = (dfLoRight + dfLoLeft) / 2.0;
388 : else
389 0 : dfBottom = dfLoRight;
390 : }
391 : else
392 : {
393 0 : dfBottom = dfLoLeft;;
394 0 : dfRight = dfUpRight;
395 : }
396 :
397 : /* -------------------------------------------------------------------- */
398 : /* Process any quadrants that aren't "nodata" anchored. */
399 : /* -------------------------------------------------------------------- */
400 4186 : CPLErr eErr = CE_None;
401 :
402 4186 : if( dfUpLeft != dfNoDataValue && iPixel > 0 && iLine > 0 )
403 : {
404 : eErr = ProcessRect( dfUpLeft, iPixel - 0.5, iLine - 0.5,
405 : dfLeft, iPixel - 0.5, iLine,
406 : dfCenter, iPixel, iLine,
407 2085 : dfTop, iPixel, iLine - 0.5 );
408 : }
409 :
410 4186 : if( dfLoLeft != dfNoDataValue && eErr == CE_None
411 : && iPixel > 0 && iLine < nHeight )
412 : {
413 : eErr = ProcessRect( dfLeft, iPixel - 0.5, iLine,
414 : dfLoLeft, iPixel - 0.5, iLine + 0.5,
415 : dfBottom, iPixel, iLine + 0.5,
416 2085 : dfCenter, iPixel, iLine );
417 : }
418 :
419 4186 : if( dfLoRight != dfNoDataValue && iPixel < nWidth && iLine < nHeight )
420 : {
421 : eErr = ProcessRect( dfCenter, iPixel, iLine,
422 : dfBottom, iPixel, iLine + 0.5,
423 : dfLoRight, iPixel + 0.5, iLine + 0.5,
424 2085 : dfRight, iPixel + 0.5, iLine );
425 : }
426 :
427 4186 : if( dfUpRight != dfNoDataValue && iPixel < nWidth && iLine > 0 )
428 : {
429 : eErr = ProcessRect( dfTop, iPixel, iLine - 0.5,
430 : dfCenter, iPixel, iLine,
431 : dfRight, iPixel + 0.5, iLine,
432 2085 : dfUpRight, iPixel + 0.5, iLine - 0.5 );
433 : }
434 :
435 4186 : return eErr;
436 : }
437 :
438 : /************************************************************************/
439 : /* ProcessRect() */
440 : /************************************************************************/
441 :
442 163557 : CPLErr GDALContourGenerator::ProcessRect(
443 : double dfUpLeft, double dfUpLeftX, double dfUpLeftY,
444 : double dfLoLeft, double dfLoLeftX, double dfLoLeftY,
445 : double dfLoRight, double dfLoRightX, double dfLoRightY,
446 : double dfUpRight, double dfUpRightX, double dfUpRightY )
447 :
448 : {
449 : /* -------------------------------------------------------------------- */
450 : /* Identify the range of elevations over this rect. */
451 : /* -------------------------------------------------------------------- */
452 : int iStartLevel, iEndLevel;
453 :
454 163557 : double dfMin = MIN(MIN(dfUpLeft,dfUpRight),MIN(dfLoLeft,dfLoRight));
455 163557 : double dfMax = MAX(MAX(dfUpLeft,dfUpRight),MAX(dfLoLeft,dfLoRight));
456 :
457 :
458 : /* -------------------------------------------------------------------- */
459 : /* Compute the set of levels to compute contours for. */
460 : /* -------------------------------------------------------------------- */
461 :
462 : /*
463 : ** If we are using fixed levels, then find the min/max in the levels
464 : ** table.
465 : */
466 163557 : if( bFixedLevels )
467 : {
468 53114 : int nStart=0, nEnd=nLevelCount-1, nMiddle;
469 :
470 53114 : iStartLevel = -1;
471 211734 : while( nStart <= nEnd )
472 : {
473 106228 : nMiddle = (nEnd + nStart) / 2;
474 :
475 106228 : double dfMiddleLevel = papoLevels[nMiddle]->GetLevel();
476 :
477 106228 : if( dfMiddleLevel < dfMin )
478 12482 : nStart = nMiddle + 1;
479 93746 : else if( dfMiddleLevel > dfMin )
480 93024 : nEnd = nMiddle - 1;
481 : else
482 : {
483 722 : iStartLevel = nMiddle;
484 722 : break;
485 : }
486 : }
487 :
488 53114 : if( iStartLevel == -1 )
489 52392 : iStartLevel = nEnd + 1;
490 :
491 53114 : iEndLevel = iStartLevel;
492 156300 : while( iEndLevel < nLevelCount-1
493 50072 : && papoLevels[iEndLevel+1]->GetLevel() < dfMax )
494 0 : iEndLevel++;
495 :
496 53114 : if( iStartLevel >= nLevelCount )
497 0 : return CE_None;
498 :
499 53114 : CPLAssert( iStartLevel >= 0 && iStartLevel < nLevelCount );
500 53114 : CPLAssert( iEndLevel >= 0 && iEndLevel < nLevelCount );
501 : }
502 :
503 : /*
504 : ** Otherwise figure out the start and end using the base and offset.
505 : */
506 : else
507 : {
508 : iStartLevel = (int)
509 110443 : ceil((dfMin - dfContourOffset) / dfContourInterval);
510 : iEndLevel = (int)
511 110443 : floor((dfMax - dfContourOffset) / dfContourInterval);
512 : }
513 :
514 163557 : if( iStartLevel > iEndLevel )
515 98333 : return CE_None;
516 :
517 : /* -------------------------------------------------------------------- */
518 : /* Loop over them. */
519 : /* -------------------------------------------------------------------- */
520 : int iLevel;
521 :
522 137639 : for( iLevel = iStartLevel; iLevel <= iEndLevel; iLevel++ )
523 : {
524 : double dfLevel;
525 :
526 72415 : if( bFixedLevels )
527 53114 : dfLevel = papoLevels[iLevel]->GetLevel();
528 : else
529 19301 : dfLevel = iLevel * dfContourInterval + dfContourOffset;
530 :
531 72415 : int nPoints = 0;
532 : double adfX[4], adfY[4];
533 72415 : CPLErr eErr = CE_None;
534 :
535 : /* Logs how many points we have af left + bottom,
536 : ** and left + bottom + right.
537 : */
538 72415 : int nPoints1 = 0, nPoints2 = 0, nPoints3 = 0;
539 :
540 :
541 : Intersect( dfUpLeft, dfUpLeftX, dfUpLeftY,
542 : dfLoLeft, dfLoLeftX, dfLoLeftY,
543 72415 : dfLoRight, dfLevel, &nPoints, adfX, adfY );
544 72415 : nPoints1 = nPoints;
545 : Intersect( dfLoLeft, dfLoLeftX, dfLoLeftY,
546 : dfLoRight, dfLoRightX, dfLoRightY,
547 72415 : dfUpRight, dfLevel, &nPoints, adfX, adfY );
548 72415 : nPoints2 = nPoints;
549 : Intersect( dfLoRight, dfLoRightX, dfLoRightY,
550 : dfUpRight, dfUpRightX, dfUpRightY,
551 72415 : dfUpLeft, dfLevel, &nPoints, adfX, adfY );
552 72415 : nPoints3 = nPoints;
553 : Intersect( dfUpRight, dfUpRightX, dfUpRightY,
554 : dfUpLeft, dfUpLeftX, dfUpLeftY,
555 72415 : dfLoLeft, dfLevel, &nPoints, adfX, adfY );
556 :
557 72415 : if( nPoints == 1 || nPoints == 3 )
558 8 : CPLDebug( "CONTOUR", "Got nPoints = %d", nPoints );
559 :
560 72415 : if( nPoints >= 2 )
561 : {
562 23530 : if ( nPoints1 == 1 && nPoints2 == 2) // left + bottom
563 : {
564 : eErr = AddSegment( dfLevel,
565 : adfX[0], adfY[0], adfX[1], adfY[1],
566 3117 : dfUpRight > dfLoLeft );
567 : }
568 22881 : else if ( nPoints1 == 1 && nPoints3 == 2 ) // left + right
569 : {
570 : eErr = AddSegment( dfLevel,
571 : adfX[0], adfY[0], adfX[1], adfY[1],
572 5585 : dfUpLeft > dfLoRight );
573 : }
574 14887 : else if ( nPoints1 == 1 && nPoints == 2 ) // left + top
575 : { // Do not do vertical contours on the left, due to symmetry
576 3176 : if ( !(dfUpLeft == dfLevel && dfLoLeft == dfLevel) )
577 : eErr = AddSegment( dfLevel,
578 : adfX[0], adfY[0], adfX[1], adfY[1],
579 3115 : dfUpLeft > dfLoRight );
580 : }
581 11611 : else if( nPoints2 == 1 && nPoints3 == 2) // bottom + right
582 : {
583 : eErr = AddSegment( dfLevel,
584 : adfX[0], adfY[0], adfX[1], adfY[1],
585 3076 : dfUpLeft > dfLoRight );
586 : }
587 8277 : else if ( nPoints2 == 1 && nPoints == 2 ) // bottom + top
588 : {
589 : eErr = AddSegment( dfLevel,
590 : adfX[0], adfY[0], adfX[1], adfY[1],
591 2818 : dfLoLeft > dfUpRight );
592 : }
593 5282 : else if ( nPoints3 == 1 && nPoints == 2 ) // right + top
594 : { // Do not do horizontal contours on upside, due to symmetry
595 2641 : if ( !(dfUpRight == dfLevel && dfUpLeft == dfLevel) )
596 : eErr = AddSegment( dfLevel,
597 : adfX[0], adfY[0], adfX[1], adfY[1],
598 2586 : dfLoLeft > dfUpRight );
599 : }
600 : else
601 : {
602 : // If we get here it is a serious error!
603 0 : CPLDebug( "CONTOUR", "Contour state not implemented!");
604 : }
605 :
606 20413 : if( eErr != CE_None )
607 0 : return eErr;
608 : }
609 :
610 72415 : if( nPoints == 4 )
611 : {
612 : // Do not do horizontal contours on upside, due to symmetry
613 461 : if ( !(dfUpRight == dfLevel && dfUpLeft == dfLevel) )
614 : {
615 : /* -------------------------------------------------------------------- */
616 : /* If we get here, we know the first was left+bottom, */
617 : /* so we are at right+top, therefore "left is high" */
618 : /* if low-left is larger than up-right. */
619 : /* We do not do a diagonal check here as we are dealing with */
620 : /* a saddle point. */
621 : /* -------------------------------------------------------------------- */
622 : eErr = AddSegment( dfLevel,
623 : adfX[2], adfY[2], adfX[3], adfY[3],
624 461 : ( dfLoRight > dfUpRight) );
625 461 : if( eErr != CE_None )
626 0 : return eErr;
627 : }
628 : }
629 : }
630 :
631 65224 : return CE_None;
632 : }
633 :
634 : /************************************************************************/
635 : /* Intersect() */
636 : /************************************************************************/
637 :
638 289660 : void GDALContourGenerator::Intersect( double dfVal1, double dfX1, double dfY1,
639 : double dfVal2, double dfX2, double dfY2,
640 : double dfNext,
641 : double dfLevel, int *pnPoints,
642 : double *padfX, double *padfY )
643 :
644 : {
645 310502 : if( dfVal1 < dfLevel && dfVal2 >= dfLevel )
646 : {
647 20842 : double dfRatio = (dfLevel - dfVal1) / (dfVal2 - dfVal1);
648 :
649 20842 : padfX[*pnPoints] = dfX1 * (1.0 - dfRatio) + dfX2 * dfRatio;
650 20842 : padfY[*pnPoints] = dfY1 * (1.0 - dfRatio) + dfY2 * dfRatio;
651 20842 : (*pnPoints)++;
652 : }
653 289500 : else if( dfVal1 > dfLevel && dfVal2 <= dfLevel )
654 : {
655 20682 : double dfRatio = (dfLevel - dfVal2) / (dfVal1 - dfVal2);
656 :
657 20682 : padfX[*pnPoints] = dfX2 * (1.0 - dfRatio) + dfX1 * dfRatio;
658 20682 : padfY[*pnPoints] = dfY2 * (1.0 - dfRatio) + dfY1 * dfRatio;
659 20682 : (*pnPoints)++;
660 : }
661 248136 : else if( dfVal1 == dfLevel && dfVal2 == dfLevel && dfNext != dfLevel )
662 : {
663 232 : padfX[*pnPoints] = dfX2;
664 232 : padfY[*pnPoints] = dfY2;
665 232 : (*pnPoints)++;
666 : }
667 289660 : }
668 :
669 : /************************************************************************/
670 : /* AddSegment() */
671 : /************************************************************************/
672 :
673 20758 : CPLErr GDALContourGenerator::AddSegment( double dfLevel,
674 : double dfX1, double dfY1,
675 : double dfX2, double dfY2,
676 : int bLeftHigh)
677 :
678 : {
679 20758 : GDALContourLevel *poLevel = FindLevel( dfLevel );
680 : GDALContourItem *poTarget;
681 : int iTarget;
682 :
683 : /* -------------------------------------------------------------------- */
684 : /* Check all active contours for any that this might attach */
685 : /* to. Eventually this should be recoded to find the contours */
686 : /* of the correct level more efficiently. */
687 : /* -------------------------------------------------------------------- */
688 :
689 20758 : if( dfY1 < dfY2 )
690 5207 : iTarget = poLevel->FindContour( dfX1, dfY1 );
691 : else
692 15551 : iTarget = poLevel->FindContour( dfX2, dfY2 );
693 :
694 20758 : if( iTarget != -1 )
695 : {
696 148 : poTarget = poLevel->GetContour( iTarget );
697 :
698 148 : poTarget->AddSegment( dfX1, dfY1, dfX2, dfY2, bLeftHigh );
699 :
700 148 : poLevel->AdjustContour( iTarget );
701 :
702 148 : return CE_None;
703 : }
704 :
705 : /* -------------------------------------------------------------------- */
706 : /* No existing contour found, lets create a new one. */
707 : /* -------------------------------------------------------------------- */
708 20610 : poTarget = new GDALContourItem( dfLevel );
709 :
710 20610 : poTarget->AddSegment( dfX1, dfY1, dfX2, dfY2, bLeftHigh );
711 :
712 20610 : poLevel->InsertContour( poTarget );
713 :
714 20610 : return CE_None;
715 : }
716 :
717 : /************************************************************************/
718 : /* FeedLine() */
719 : /************************************************************************/
720 :
721 1055 : CPLErr GDALContourGenerator::FeedLine( double *padfScanline )
722 :
723 : {
724 : /* -------------------------------------------------------------------- */
725 : /* Switch current line to "lastline" slot, and copy new data */
726 : /* into new "this line". */
727 : /* -------------------------------------------------------------------- */
728 1055 : double *padfTempLine = padfLastLine;
729 1055 : padfLastLine = padfThisLine;
730 1055 : padfThisLine = padfTempLine;
731 :
732 : /* -------------------------------------------------------------------- */
733 : /* If this is the end of the lines (NULL passed in), copy the */
734 : /* last line. */
735 : /* -------------------------------------------------------------------- */
736 1055 : if( padfScanline == NULL )
737 : {
738 8 : memcpy( padfThisLine, padfLastLine, sizeof(double) * nWidth );
739 : }
740 : else
741 : {
742 1047 : memcpy( padfThisLine, padfScanline, sizeof(double) * nWidth );
743 : }
744 :
745 : /* -------------------------------------------------------------------- */
746 : /* Perturb any values that occur exactly on level boundaries. */
747 : /* -------------------------------------------------------------------- */
748 : int iPixel;
749 :
750 159403 : for( iPixel = 0; iPixel < nWidth; iPixel++ )
751 : {
752 158348 : if( bNoDataActive && padfThisLine[iPixel] == dfNoDataValue )
753 0 : continue;
754 :
755 158348 : double dfLevel = (padfThisLine[iPixel] - dfContourOffset)
756 158348 : / dfContourInterval;
757 :
758 158348 : if( dfLevel - (int) dfLevel == 0.0 )
759 : {
760 102024 : padfThisLine[iPixel] += dfContourInterval * FUDGE_EXACT;
761 : }
762 : }
763 :
764 : /* -------------------------------------------------------------------- */
765 : /* If this is the first line we need to initialize the previous */
766 : /* line from the first line of data. */
767 : /* -------------------------------------------------------------------- */
768 1055 : if( iLine == -1 )
769 : {
770 8 : memcpy( padfLastLine, padfThisLine, sizeof(double) * nWidth );
771 8 : iLine = 0;
772 : }
773 :
774 : /* -------------------------------------------------------------------- */
775 : /* Clear the recently used flags on the contours so we can */
776 : /* check later which ones were touched for this scanline. */
777 : /* -------------------------------------------------------------------- */
778 : int iLevel, iContour;
779 :
780 8118 : for( iLevel = 0; iLevel < nLevelCount; iLevel++ )
781 : {
782 7063 : GDALContourLevel *poLevel = papoLevels[iLevel];
783 :
784 35559 : for( iContour = 0; iContour < poLevel->GetContourCount(); iContour++ )
785 28496 : poLevel->GetContour( iContour )->bRecentlyAccessed = FALSE;
786 : }
787 :
788 : /* -------------------------------------------------------------------- */
789 : /* Process each pixel. */
790 : /* -------------------------------------------------------------------- */
791 160458 : for( iPixel = 0; iPixel < nWidth+1; iPixel++ )
792 : {
793 159403 : CPLErr eErr = ProcessPixel( iPixel );
794 159403 : if( eErr != CE_None )
795 0 : return eErr;
796 : }
797 :
798 : /* -------------------------------------------------------------------- */
799 : /* eject any pending contours. */
800 : /* -------------------------------------------------------------------- */
801 1055 : CPLErr eErr = EjectContours( padfScanline != NULL );
802 :
803 1055 : iLine++;
804 :
805 1055 : if( iLine == nHeight && eErr == CE_None )
806 8 : return FeedLine( NULL );
807 : else
808 1047 : return eErr;
809 : }
810 :
811 : /************************************************************************/
812 : /* EjectContours() */
813 : /************************************************************************/
814 :
815 1055 : CPLErr GDALContourGenerator::EjectContours( int bOnlyUnused )
816 :
817 : {
818 : int iLevel;
819 1055 : CPLErr eErr = CE_None;
820 :
821 : /* -------------------------------------------------------------------- */
822 : /* Process all contours of all levels that match our criteria */
823 : /* -------------------------------------------------------------------- */
824 8172 : for( iLevel = 0; iLevel < nLevelCount && eErr == CE_None; iLevel++ )
825 : {
826 7117 : GDALContourLevel *poLevel = papoLevels[iLevel];
827 : int iContour;
828 :
829 63340 : for( iContour = 0;
830 : iContour < poLevel->GetContourCount() && eErr == CE_None;
831 : /* increment in loop if we don't consume it. */ )
832 : {
833 : int iC2;
834 49106 : GDALContourItem *poTarget = poLevel->GetContour( iContour );
835 :
836 49106 : if( bOnlyUnused && poTarget->bRecentlyAccessed )
837 : {
838 28496 : iContour++;
839 28496 : continue;
840 : }
841 :
842 20610 : poLevel->RemoveContour( iContour );
843 :
844 : // Try to find another contour we can merge with in this level.
845 :
846 87852 : for( iC2 = 0; iC2 < poLevel->GetContourCount(); iC2++ )
847 : {
848 87095 : GDALContourItem *poOther = poLevel->GetContour( iC2 );
849 :
850 87095 : if( poOther->Merge( poTarget ) )
851 19853 : break;
852 : }
853 :
854 : // If we didn't merge it, then eject (write) it out.
855 20610 : if( iC2 == poLevel->GetContourCount() )
856 : {
857 757 : if( pfnWriter != NULL )
858 : {
859 : // If direction is wrong, then reverse before ejecting.
860 757 : poTarget->PrepareEjection();
861 :
862 : eErr = pfnWriter( poTarget->dfLevel, poTarget->nPoints,
863 : poTarget->padfX, poTarget->padfY,
864 757 : pWriterCBData );
865 : }
866 : }
867 :
868 20610 : delete poTarget;
869 : }
870 : }
871 :
872 1055 : return eErr;
873 : }
874 :
875 : /************************************************************************/
876 : /* FindLevel() */
877 : /************************************************************************/
878 :
879 20764 : GDALContourLevel *GDALContourGenerator::FindLevel( double dfLevel )
880 :
881 : {
882 20764 : int nStart=0, nEnd=nLevelCount-1, nMiddle;
883 :
884 : /* -------------------------------------------------------------------- */
885 : /* Binary search to find the requested level. */
886 : /* -------------------------------------------------------------------- */
887 99346 : while( nStart <= nEnd )
888 : {
889 78522 : nMiddle = (nEnd + nStart) / 2;
890 :
891 78522 : double dfMiddleLevel = papoLevels[nMiddle]->GetLevel();
892 :
893 78522 : if( dfMiddleLevel < dfLevel )
894 23975 : nStart = nMiddle + 1;
895 54547 : else if( dfMiddleLevel > dfLevel )
896 33843 : nEnd = nMiddle - 1;
897 : else
898 20704 : return papoLevels[nMiddle];
899 : }
900 :
901 : /* -------------------------------------------------------------------- */
902 : /* Didn't find the level, create a new one and insert it in */
903 : /* order. */
904 : /* -------------------------------------------------------------------- */
905 60 : GDALContourLevel *poLevel = new GDALContourLevel( dfLevel );
906 :
907 60 : if( nLevelMax == nLevelCount )
908 : {
909 10 : nLevelMax = nLevelMax * 2 + 10;
910 : papoLevels = (GDALContourLevel **)
911 10 : CPLRealloc( papoLevels, sizeof(void*) * nLevelMax );
912 : }
913 :
914 60 : if( nLevelCount - nEnd - 1 > 0 )
915 : memmove( papoLevels + nEnd + 2, papoLevels + nEnd + 1,
916 28 : (nLevelCount - nEnd - 1) * sizeof(void*) );
917 60 : papoLevels[nEnd+1] = poLevel;
918 60 : nLevelCount++;
919 :
920 60 : return poLevel;
921 : }
922 :
923 : /************************************************************************/
924 : /* ==================================================================== */
925 : /* GDALContourLevel */
926 : /* ==================================================================== */
927 : /************************************************************************/
928 :
929 : /************************************************************************/
930 : /* GDALContourLevel() */
931 : /************************************************************************/
932 :
933 60 : GDALContourLevel::GDALContourLevel( double dfLevelIn )
934 :
935 : {
936 60 : dfLevel = dfLevelIn;
937 60 : nEntryMax = 0;
938 60 : nEntryCount = 0;
939 60 : papoEntries = NULL;
940 60 : }
941 :
942 : /************************************************************************/
943 : /* ~GDALContourLevel() */
944 : /************************************************************************/
945 :
946 60 : GDALContourLevel::~GDALContourLevel()
947 :
948 : {
949 60 : CPLAssert( nEntryCount == 0 );
950 60 : CPLFree( papoEntries );
951 60 : }
952 :
953 : /************************************************************************/
954 : /* AdjustContour() */
955 : /* */
956 : /* Assume the indicated contour's tail may have changed, and */
957 : /* adjust it up or down in the list of contours to re-establish */
958 : /* proper ordering. */
959 : /************************************************************************/
960 :
961 148 : void GDALContourLevel::AdjustContour( int iChanged )
962 :
963 : {
964 296 : while( iChanged > 0
965 0 : && papoEntries[iChanged]->dfTailX < papoEntries[iChanged-1]->dfTailX )
966 : {
967 0 : GDALContourItem *poTemp = papoEntries[iChanged];
968 0 : papoEntries[iChanged] = papoEntries[iChanged-1];
969 0 : papoEntries[iChanged-1] = poTemp;
970 0 : iChanged--;
971 : }
972 :
973 411 : while( iChanged < nEntryCount-1
974 106 : && papoEntries[iChanged]->dfTailX > papoEntries[iChanged+1]->dfTailX )
975 : {
976 9 : GDALContourItem *poTemp = papoEntries[iChanged];
977 9 : papoEntries[iChanged] = papoEntries[iChanged+1];
978 9 : papoEntries[iChanged+1] = poTemp;
979 9 : iChanged++;
980 : }
981 148 : }
982 :
983 : /************************************************************************/
984 : /* RemoveContour() */
985 : /************************************************************************/
986 :
987 20610 : void GDALContourLevel::RemoveContour( int iTarget )
988 :
989 : {
990 20610 : if( iTarget < nEntryCount )
991 : memmove( papoEntries + iTarget, papoEntries + iTarget + 1,
992 20610 : (nEntryCount - iTarget - 1) * sizeof(void*) );
993 20610 : nEntryCount--;
994 20610 : }
995 :
996 : /************************************************************************/
997 : /* FindContour() */
998 : /* */
999 : /* Perform a binary search to find the requested "tail" */
1000 : /* location. If not available return -1. In theory there can */
1001 : /* be more than one contour with the same tail X and different */
1002 : /* Y tails ... ensure we check against them all. */
1003 : /************************************************************************/
1004 :
1005 20758 : int GDALContourLevel::FindContour( double dfX, double dfY )
1006 :
1007 : {
1008 20758 : int nStart = 0, nEnd = nEntryCount-1, nMiddle;
1009 :
1010 114454 : while( nEnd >= nStart )
1011 : {
1012 77566 : nMiddle = (nEnd + nStart) / 2;
1013 :
1014 77566 : double dfMiddleX = papoEntries[nMiddle]->dfTailX;
1015 :
1016 77566 : if( dfMiddleX < dfX )
1017 43829 : nStart = nMiddle + 1;
1018 33737 : else if( dfMiddleX > dfX )
1019 29109 : nEnd = nMiddle - 1;
1020 : else
1021 : {
1022 21073 : while( nMiddle > 0
1023 7832 : && fabs(papoEntries[nMiddle]->dfTailX-dfX) < JOIN_DIST )
1024 3985 : nMiddle--;
1025 :
1026 15366 : while( nMiddle < nEntryCount
1027 5133 : && fabs(papoEntries[nMiddle]->dfTailX-dfX) < JOIN_DIST )
1028 : {
1029 1125 : if( fabs(papoEntries[nMiddle]->padfY[papoEntries[nMiddle]->nPoints-1] - dfY) < JOIN_DIST )
1030 148 : return nMiddle;
1031 977 : nMiddle++;
1032 : }
1033 :
1034 4480 : return -1;
1035 : }
1036 : }
1037 :
1038 16130 : return -1;
1039 : }
1040 :
1041 : /************************************************************************/
1042 : /* InsertContour() */
1043 : /* */
1044 : /* Ensure the newly added contour is placed in order according */
1045 : /* to the X value relative to the other contours. */
1046 : /************************************************************************/
1047 :
1048 20610 : int GDALContourLevel::InsertContour( GDALContourItem *poNewContour )
1049 :
1050 : {
1051 : /* -------------------------------------------------------------------- */
1052 : /* Find where to insert by binary search. */
1053 : /* -------------------------------------------------------------------- */
1054 20610 : int nStart = 0, nEnd = nEntryCount-1, nMiddle;
1055 :
1056 127431 : while( nEnd >= nStart )
1057 : {
1058 86211 : nMiddle = (nEnd + nStart) / 2;
1059 :
1060 86211 : double dfMiddleX = papoEntries[nMiddle]->dfTailX;
1061 :
1062 86211 : if( dfMiddleX < poNewContour->dfLevel )
1063 77279 : nStart = nMiddle + 1;
1064 8932 : else if( dfMiddleX > poNewContour->dfLevel )
1065 8932 : nEnd = nMiddle - 1;
1066 : else
1067 : {
1068 0 : nEnd = nMiddle - 1;
1069 0 : break;
1070 : }
1071 : }
1072 :
1073 : /* -------------------------------------------------------------------- */
1074 : /* Do we need to grow the array? */
1075 : /* -------------------------------------------------------------------- */
1076 20610 : if( nEntryMax == nEntryCount )
1077 : {
1078 160 : nEntryMax = nEntryMax * 2 + 10;
1079 : papoEntries = (GDALContourItem **)
1080 160 : CPLRealloc( papoEntries, sizeof(void*) * nEntryMax );
1081 : }
1082 :
1083 : /* -------------------------------------------------------------------- */
1084 : /* Insert the new contour at the appropriate location. */
1085 : /* -------------------------------------------------------------------- */
1086 20610 : if( nEntryCount - nEnd - 1 > 0 )
1087 : memmove( papoEntries + nEnd + 2, papoEntries + nEnd + 1,
1088 3718 : (nEntryCount - nEnd - 1) * sizeof(void*) );
1089 20610 : papoEntries[nEnd+1] = poNewContour;
1090 20610 : nEntryCount++;
1091 :
1092 20610 : return nEnd+1;
1093 : }
1094 :
1095 :
1096 : /************************************************************************/
1097 : /* ==================================================================== */
1098 : /* GDALContourItem */
1099 : /* ==================================================================== */
1100 : /************************************************************************/
1101 :
1102 : /************************************************************************/
1103 : /* GDALContourItem() */
1104 : /************************************************************************/
1105 :
1106 20610 : GDALContourItem::GDALContourItem( double dfLevelIn )
1107 :
1108 : {
1109 20610 : dfLevel = dfLevelIn;
1110 20610 : bRecentlyAccessed = FALSE;
1111 20610 : nPoints = 0;
1112 20610 : nMaxPoints = 0;
1113 20610 : padfX = NULL;
1114 20610 : padfY = NULL;
1115 :
1116 20610 : bLeftIsHigh = FALSE;
1117 :
1118 20610 : dfTailX = 0.0;
1119 20610 : }
1120 :
1121 : /************************************************************************/
1122 : /* ~GDALContourItem() */
1123 : /************************************************************************/
1124 :
1125 20610 : GDALContourItem::~GDALContourItem()
1126 :
1127 : {
1128 20610 : CPLFree( padfX );
1129 20610 : CPLFree( padfY );
1130 20610 : }
1131 :
1132 : /************************************************************************/
1133 : /* AddSegment() */
1134 : /************************************************************************/
1135 :
1136 20758 : int GDALContourItem::AddSegment( double dfXStart, double dfYStart,
1137 : double dfXEnd, double dfYEnd,
1138 : int bLeftHigh)
1139 :
1140 : {
1141 20758 : MakeRoomFor( nPoints + 1 );
1142 :
1143 : /* -------------------------------------------------------------------- */
1144 : /* If there are no segments, just add now. */
1145 : /* -------------------------------------------------------------------- */
1146 20758 : if( nPoints == 0 )
1147 : {
1148 20610 : nPoints = 2;
1149 :
1150 20610 : padfX[0] = dfXStart;
1151 20610 : padfY[0] = dfYStart;
1152 20610 : padfX[1] = dfXEnd;
1153 20610 : padfY[1] = dfYEnd;
1154 20610 : bRecentlyAccessed = TRUE;
1155 :
1156 20610 : dfTailX = padfX[1];
1157 :
1158 : // Here we know that the left of this vector is the high side
1159 20610 : bLeftIsHigh = bLeftHigh;
1160 :
1161 20610 : return TRUE;
1162 : }
1163 :
1164 : /* -------------------------------------------------------------------- */
1165 : /* Try to matching up with one of the ends, and insert. */
1166 : /* -------------------------------------------------------------------- */
1167 245 : if( fabs(padfX[nPoints-1]-dfXStart) < JOIN_DIST
1168 97 : && fabs(padfY[nPoints-1]-dfYStart) < JOIN_DIST )
1169 : {
1170 97 : padfX[nPoints] = dfXEnd;
1171 97 : padfY[nPoints] = dfYEnd;
1172 97 : nPoints++;
1173 :
1174 97 : bRecentlyAccessed = TRUE;
1175 :
1176 97 : dfTailX = dfXEnd;
1177 :
1178 97 : return TRUE;
1179 : }
1180 102 : else if( fabs(padfX[nPoints-1]-dfXEnd) < JOIN_DIST
1181 51 : && fabs(padfY[nPoints-1]-dfYEnd) < JOIN_DIST )
1182 : {
1183 51 : padfX[nPoints] = dfXStart;
1184 51 : padfY[nPoints] = dfYStart;
1185 51 : nPoints++;
1186 :
1187 51 : bRecentlyAccessed = TRUE;
1188 :
1189 51 : dfTailX = dfXStart;
1190 :
1191 51 : return TRUE;
1192 : }
1193 : else
1194 0 : return FALSE;
1195 : }
1196 :
1197 : /************************************************************************/
1198 : /* Merge() */
1199 : /************************************************************************/
1200 :
1201 87095 : int GDALContourItem::Merge( GDALContourItem *poOther )
1202 :
1203 : {
1204 87095 : if( poOther->dfLevel != dfLevel )
1205 0 : return FALSE;
1206 :
1207 : /* -------------------------------------------------------------------- */
1208 : /* Try to matching up with one of the ends, and insert. */
1209 : /* -------------------------------------------------------------------- */
1210 96681 : if( fabs(padfX[nPoints-1]-poOther->padfX[0]) < JOIN_DIST
1211 9586 : && fabs(padfY[nPoints-1]-poOther->padfY[0]) < JOIN_DIST )
1212 : {
1213 8867 : MakeRoomFor( nPoints + poOther->nPoints - 1 );
1214 :
1215 : memcpy( padfX + nPoints, poOther->padfX + 1,
1216 8867 : sizeof(double) * (poOther->nPoints-1) );
1217 : memcpy( padfY + nPoints, poOther->padfY + 1,
1218 8867 : sizeof(double) * (poOther->nPoints-1) );
1219 8867 : nPoints += poOther->nPoints - 1;
1220 :
1221 8867 : bRecentlyAccessed = TRUE;
1222 :
1223 8867 : dfTailX = padfX[nPoints-1];
1224 :
1225 8867 : return TRUE;
1226 : }
1227 83602 : else if( fabs(padfX[0]-poOther->padfX[poOther->nPoints-1]) < JOIN_DIST
1228 5374 : && fabs(padfY[0]-poOther->padfY[poOther->nPoints-1]) < JOIN_DIST )
1229 : {
1230 5035 : MakeRoomFor( nPoints + poOther->nPoints - 1 );
1231 :
1232 : memmove( padfX + poOther->nPoints - 1, padfX,
1233 5035 : sizeof(double) * nPoints );
1234 : memmove( padfY + poOther->nPoints - 1, padfY,
1235 5035 : sizeof(double) * nPoints );
1236 : memcpy( padfX, poOther->padfX,
1237 5035 : sizeof(double) * (poOther->nPoints-1) );
1238 : memcpy( padfY, poOther->padfY,
1239 5035 : sizeof(double) * (poOther->nPoints-1) );
1240 5035 : nPoints += poOther->nPoints - 1;
1241 :
1242 5035 : bRecentlyAccessed = TRUE;
1243 :
1244 5035 : dfTailX = padfX[nPoints-1];
1245 :
1246 5035 : return TRUE;
1247 : }
1248 77117 : else if( fabs(padfX[nPoints-1]-poOther->padfX[poOther->nPoints-1]) < JOIN_DIST
1249 3924 : && fabs(padfY[nPoints-1]-poOther->padfY[poOther->nPoints-1]) < JOIN_DIST )
1250 : {
1251 : int i;
1252 :
1253 3607 : MakeRoomFor( nPoints + poOther->nPoints - 1 );
1254 :
1255 19445 : for( i = 0; i < poOther->nPoints-1; i++ )
1256 : {
1257 15838 : padfX[i+nPoints] = poOther->padfX[poOther->nPoints-i-2];
1258 15838 : padfY[i+nPoints] = poOther->padfY[poOther->nPoints-i-2];
1259 : }
1260 :
1261 3607 : nPoints += poOther->nPoints - 1;
1262 :
1263 3607 : bRecentlyAccessed = TRUE;
1264 :
1265 3607 : dfTailX = padfX[nPoints-1];
1266 :
1267 3607 : return TRUE;
1268 : }
1269 72219 : else if( fabs(padfX[0]-poOther->padfX[0]) < JOIN_DIST
1270 2633 : && fabs(padfY[0]-poOther->padfY[0]) < JOIN_DIST )
1271 : {
1272 : int i;
1273 :
1274 2344 : MakeRoomFor( nPoints + poOther->nPoints - 1 );
1275 :
1276 : memmove( padfX + poOther->nPoints - 1, padfX,
1277 2344 : sizeof(double) * nPoints );
1278 : memmove( padfY + poOther->nPoints - 1, padfY,
1279 2344 : sizeof(double) * nPoints );
1280 :
1281 17190 : for( i = 0; i < poOther->nPoints-1; i++ )
1282 : {
1283 14846 : padfX[i] = poOther->padfX[poOther->nPoints - i - 1];
1284 14846 : padfY[i] = poOther->padfY[poOther->nPoints - i - 1];
1285 : }
1286 :
1287 2344 : nPoints += poOther->nPoints - 1;
1288 :
1289 2344 : bRecentlyAccessed = TRUE;
1290 :
1291 2344 : dfTailX = padfX[nPoints-1];
1292 :
1293 2344 : return TRUE;
1294 : }
1295 : else
1296 67242 : return FALSE;
1297 : }
1298 :
1299 : /************************************************************************/
1300 : /* MakeRoomFor() */
1301 : /************************************************************************/
1302 :
1303 40611 : void GDALContourItem::MakeRoomFor( int nNewPoints )
1304 :
1305 : {
1306 40611 : if( nNewPoints > nMaxPoints )
1307 : {
1308 23223 : nMaxPoints = nNewPoints * 2 + 50;
1309 23223 : padfX = (double *) CPLRealloc(padfX,sizeof(double) * nMaxPoints);
1310 23223 : padfY = (double *) CPLRealloc(padfY,sizeof(double) * nMaxPoints);
1311 : }
1312 40611 : }
1313 :
1314 : /************************************************************************/
1315 : /* PrepareEjection() */
1316 : /************************************************************************/
1317 :
1318 757 : void GDALContourItem::PrepareEjection()
1319 :
1320 : {
1321 : /* If left side is the high side, then reverse to get curve normal
1322 : ** pointing downwards
1323 : */
1324 757 : if( bLeftIsHigh )
1325 : {
1326 : int i;
1327 :
1328 : // Reverse the arrays
1329 8088 : for( i = 0; i < nPoints / 2; i++ )
1330 : {
1331 : double dfTemp;
1332 7570 : dfTemp = padfX[i];
1333 7570 : padfX[i] = padfX[ nPoints - i - 1];
1334 7570 : padfX[ nPoints - i - 1] = dfTemp;
1335 :
1336 7570 : dfTemp = padfY[i];
1337 7570 : padfY[i] = padfY[ nPoints - i - 1];
1338 7570 : padfY[ nPoints - i - 1] = dfTemp;
1339 : }
1340 : }
1341 757 : }
1342 :
1343 :
1344 : /************************************************************************/
1345 : /* ==================================================================== */
1346 : /* Additional C Callable Functions */
1347 : /* ==================================================================== */
1348 : /************************************************************************/
1349 :
1350 : /************************************************************************/
1351 : /* OGRContourWriter() */
1352 : /************************************************************************/
1353 :
1354 757 : CPLErr OGRContourWriter( double dfLevel,
1355 : int nPoints, double *padfX, double *padfY,
1356 : void *pInfo )
1357 :
1358 : {
1359 757 : OGRContourWriterInfo *poInfo = (OGRContourWriterInfo *) pInfo;
1360 : OGRFeatureH hFeat;
1361 : OGRGeometryH hGeom;
1362 : int iPoint;
1363 :
1364 757 : hFeat = OGR_F_Create( OGR_L_GetLayerDefn( (OGRLayerH) poInfo->hLayer ) );
1365 :
1366 757 : if( poInfo->nIDField != -1 )
1367 757 : OGR_F_SetFieldInteger( hFeat, poInfo->nIDField, poInfo->nNextID++ );
1368 :
1369 757 : if( poInfo->nElevField != -1 )
1370 139 : OGR_F_SetFieldDouble( hFeat, poInfo->nElevField, dfLevel );
1371 :
1372 757 : hGeom = OGR_G_CreateGeometry( wkbLineString );
1373 :
1374 22272 : for( iPoint = nPoints-1; iPoint >= 0; iPoint-- )
1375 : {
1376 : OGR_G_SetPoint( hGeom, iPoint,
1377 21515 : poInfo->adfGeoTransform[0]
1378 21515 : + poInfo->adfGeoTransform[1] * padfX[iPoint]
1379 21515 : + poInfo->adfGeoTransform[2] * padfY[iPoint],
1380 21515 : poInfo->adfGeoTransform[3]
1381 21515 : + poInfo->adfGeoTransform[4] * padfX[iPoint]
1382 21515 : + poInfo->adfGeoTransform[5] * padfY[iPoint],
1383 86060 : dfLevel );
1384 : }
1385 :
1386 757 : OGR_F_SetGeometryDirectly( hFeat, hGeom );
1387 :
1388 757 : OGR_L_CreateFeature( (OGRLayerH) poInfo->hLayer, hFeat );
1389 757 : OGR_F_Destroy( hFeat );
1390 :
1391 757 : return CE_None;
1392 : }
1393 : #endif // OGR_ENABLED
1394 :
1395 : /************************************************************************/
1396 : /* GDALContourGenerate() */
1397 : /************************************************************************/
1398 :
1399 : /**
1400 : * Create vector contours from raster DEM.
1401 : *
1402 : * This algorithm will generate contours vectors for the input raster band
1403 : * on the requested set of contour levels. The vector contours are written
1404 : * to the passed in OGR vector layer. Also, a NODATA value may be specified
1405 : * to identify pixels that should not be considered in contour line generation.
1406 : *
1407 : * The gdal/apps/gdal_contour.cpp mainline can be used as an example of
1408 : * how to use this function.
1409 : *
1410 : * ALGORITHM RULES
1411 :
1412 : For contouring purposes raster pixel values are assumed to represent a point
1413 : value at the center of the corresponding pixel region. For the purpose of
1414 : contour generation we virtually connect each pixel center to the values to
1415 : the left, right, top and bottom. We assume that the pixel value is linearly
1416 : interpolated between the pixel centers along each line, and determine where
1417 : (if any) contour lines will appear onlong these line segements. Then the
1418 : contour crossings are connected.
1419 :
1420 : This means that contour lines nodes won't actually be on pixel edges, but
1421 : rather along vertical and horizontal lines connecting the pixel centers.
1422 :
1423 : \verbatim
1424 : General Case:
1425 :
1426 : 5 | | 3
1427 : -- + ---------------- + --
1428 : | |
1429 : | |
1430 : | |
1431 : | |
1432 : 10 + |
1433 : |\ |
1434 : | \ |
1435 : -- + -+-------------- + --
1436 : 12 | 10 | 1
1437 :
1438 :
1439 : Saddle Point:
1440 :
1441 : 5 | | 12
1442 : -- + -------------+-- + --
1443 : | \ |
1444 : | \|
1445 : | +
1446 : | |
1447 : + |
1448 : |\ |
1449 : | \ |
1450 : -- + -+-------------- + --
1451 : 12 | | 1
1452 :
1453 : or:
1454 :
1455 : 5 | | 12
1456 : -- + -------------+-- + --
1457 : | __/ |
1458 : | ___/ |
1459 : | ___/ __+
1460 : | / __/ |
1461 : +' __/ |
1462 : | __/ |
1463 : | ,__/ |
1464 : -- + -+-------------- + --
1465 : 12 | | 1
1466 : \endverbatim
1467 :
1468 : Nodata:
1469 :
1470 : In the "nodata" case we treat the whole nodata pixel as a no-mans land.
1471 : We extend the corner pixels near the nodata out to half way and then
1472 : construct extra lines from those points to the center which is assigned
1473 : an averaged value from the two nearby points (in this case (12+3+5)/3).
1474 :
1475 : \verbatim
1476 : 5 | | 3
1477 : -- + ---------------- + --
1478 : | |
1479 : | |
1480 : | 6.7 |
1481 : | +---------+ 3
1482 : 10 +___ |
1483 : | \____+ 10
1484 : | |
1485 : -- + -------+ +
1486 : 12 | 12 (nodata)
1487 :
1488 : \endverbatim
1489 :
1490 : *
1491 : * @param hBand The band to read raster data from. The whole band will be
1492 : * processed.
1493 : *
1494 : * @param dfContourInterval The elevation interval between contours generated.
1495 : *
1496 : * @param dfContourBase The "base" relative to which contour intervals are
1497 : * applied. This is normally zero, but could be different. To generate 10m
1498 : * contours at 5, 15, 25, ... the ContourBase would be 5.
1499 : *
1500 : * @param nFixedLevelCount The number of fixed levels. If this is greater than
1501 : * zero, then fixed levels will be used, and ContourInterval and ContourBase
1502 : * are ignored.
1503 : *
1504 : * @param padfFixedLevels The list of fixed contour levels at which contours
1505 : * should be generated. It will contain FixedLevelCount entries, and may be
1506 : * NULL if fixed levels are disabled (FixedLevelCount = 0).
1507 : *
1508 : * @param bUseNoData If TRUE the dfNoDataValue will be used.
1509 : *
1510 : * @param dfNoDataValue the value to use as a "nodata" value. That is, a
1511 : * pixel value which should be ignored in generating contours as if the value
1512 : * of the pixel were not known.
1513 : *
1514 : * @param hLayer the layer to which new contour vectors will be written.
1515 : * Each contour will have a LINESTRING geometry attached to it. This
1516 : * is really of type OGRLayerH, but void * is used to avoid pulling the
1517 : * ogr_api.h file in here.
1518 : *
1519 : * @param iIDField if not -1 this will be used as a field index to indicate
1520 : * where a unique id should be written for each feature (contour) written.
1521 : *
1522 : * @param iElevField if not -1 this will be used as a field index to indicate
1523 : * where the elevation value of the contour should be written.
1524 : *
1525 : * @param pfnProgress a GDALProgressFunc that may be used to report progress
1526 : * to the user, or to interrupt the algorithm. May be NULL if not required.
1527 : *
1528 : * @param pProgressArg the callback data for the pfnProgress function.
1529 : *
1530 : * @return CE_None on success or CE_Failure if an error occurs.
1531 : */
1532 :
1533 8 : CPLErr GDALContourGenerate( GDALRasterBandH hBand,
1534 : double dfContourInterval, double dfContourBase,
1535 : int nFixedLevelCount, double *padfFixedLevels,
1536 : int bUseNoData, double dfNoDataValue,
1537 : void *hLayer, int iIDField, int iElevField,
1538 : GDALProgressFunc pfnProgress, void *pProgressArg )
1539 :
1540 : {
1541 : #ifndef OGR_ENABLED
1542 : CPLError(CE_Failure, CPLE_NotSupported, "GDALContourGenerate() unimplemented in a non OGR build");
1543 : return CE_Failure;
1544 : #else
1545 8 : VALIDATE_POINTER1( hBand, "GDALContourGenerate", CE_Failure );
1546 :
1547 : OGRContourWriterInfo oCWI;
1548 :
1549 8 : if( pfnProgress == NULL )
1550 2 : pfnProgress = GDALDummyProgress;
1551 :
1552 8 : if( !pfnProgress( 0.0, "", pProgressArg ) )
1553 : {
1554 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1555 0 : return CE_Failure;
1556 : }
1557 :
1558 : /* -------------------------------------------------------------------- */
1559 : /* Setup contour writer information. */
1560 : /* -------------------------------------------------------------------- */
1561 : GDALDatasetH hSrcDS;
1562 :
1563 8 : oCWI.hLayer = (OGRLayerH) hLayer;
1564 :
1565 8 : oCWI.nElevField = iElevField;
1566 8 : oCWI.nIDField = iIDField;
1567 :
1568 8 : hSrcDS = GDALGetBandDataset( hBand );
1569 8 : GDALGetGeoTransform( hSrcDS, oCWI.adfGeoTransform );
1570 8 : oCWI.nNextID = 0;
1571 :
1572 : /* -------------------------------------------------------------------- */
1573 : /* Setup contour generator. */
1574 : /* -------------------------------------------------------------------- */
1575 8 : int nXSize = GDALGetRasterBandXSize( hBand );
1576 8 : int nYSize = GDALGetRasterBandYSize( hBand );
1577 :
1578 8 : GDALContourGenerator oCG( nXSize, nYSize, OGRContourWriter, &oCWI );
1579 :
1580 8 : if( nFixedLevelCount > 0 )
1581 2 : oCG.SetFixedLevels( nFixedLevelCount, padfFixedLevels );
1582 : else
1583 6 : oCG.SetContourLevels( dfContourInterval, dfContourBase );
1584 :
1585 8 : if( bUseNoData )
1586 3 : oCG.SetNoData( dfNoDataValue );
1587 :
1588 : /* -------------------------------------------------------------------- */
1589 : /* Feed the data into the contour generator. */
1590 : /* -------------------------------------------------------------------- */
1591 : int iLine;
1592 : double *padfScanline;
1593 8 : CPLErr eErr = CE_None;
1594 :
1595 8 : padfScanline = (double *) VSIMalloc(sizeof(double) * nXSize);
1596 8 : if (padfScanline == NULL)
1597 : {
1598 : CPLError( CE_Failure, CPLE_OutOfMemory,
1599 0 : "VSIMalloc(): Out of memory in GDALContourGenerate" );
1600 0 : return CE_Failure;
1601 : }
1602 :
1603 1055 : for( iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ )
1604 : {
1605 : GDALRasterIO( hBand, GF_Read, 0, iLine, nXSize, 1,
1606 1047 : padfScanline, nXSize, 1, GDT_Float64, 0, 0 );
1607 1047 : eErr = oCG.FeedLine( padfScanline );
1608 :
1609 1047 : if( eErr == CE_None
1610 : && !pfnProgress( (iLine+1) / (double) nYSize, "", pProgressArg ) )
1611 : {
1612 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1613 0 : eErr = CE_Failure;
1614 : }
1615 : }
1616 :
1617 8 : CPLFree( padfScanline );
1618 :
1619 8 : return eErr;
1620 : #endif // OGR_ENABLED
1621 : }
|
