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