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