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