1 : /******************************************************************************
2 : * $Id: contour.cpp 19881 2010-06-17 18:45:14Z chaitanya $
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 19881 2010-06-17 18:45:14Z chaitanya $");
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 287936 : double GetLevel() { return dfLevel; }
91 200244 : int GetContourCount() { return nEntryCount; }
92 164845 : 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 : void SetContourLevels( double dfContourInterval,
150 6 : 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 : GDAL_CG_Create( int nWidth, int nHeight, int bNoDataSet, double dfNoDataValue,
166 : double dfContourInterval, double dfContourBase,
167 0 : 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 : GDALContourGenerator::GDALContourGenerator( int nWidthIn, int nHeightIn,
213 : GDALContourWriter pfnWriterIn,
214 8 : 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 : void GDALContourGenerator::SetFixedLevels( int nFixedLevelCount,
260 2 : 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 : 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 163557 : 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 106228 : while( iEndLevel < nLevelCount-1
491 : && 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 : 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 289660 : 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 : CPLErr GDALContourGenerator::AddSegment( double dfLevel,
672 : double dfX1, double dfY1,
673 : double dfX2, double dfY2,
674 20758 : 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 : 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 : && 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 305 : while( iChanged < nEntryCount-1
972 : && 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 13241 : while( nMiddle > 0
1021 : && fabs(papoEntries[nMiddle]->dfTailX-dfX) < JOIN_DIST )
1022 3985 : nMiddle--;
1023 :
1024 10233 : while( nMiddle < nEntryCount
1025 : && 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 : int GDALContourItem::AddSegment( double dfXStart, double dfYStart,
1135 : double dfXEnd, double dfYEnd,
1136 20758 : 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 148 : if( fabs(padfX[nPoints-1]-dfXStart) < JOIN_DIST
1166 : && 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 51 : else if( fabs(padfX[nPoints-1]-dfXEnd) < JOIN_DIST
1179 : && 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 87095 : if( fabs(padfX[nPoints-1]-poOther->padfX[0]) < JOIN_DIST
1209 : && 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 78228 : else if( fabs(padfX[0]-poOther->padfX[poOther->nPoints-1]) < JOIN_DIST
1226 : && 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 73193 : else if( fabs(padfX[nPoints-1]-poOther->padfX[poOther->nPoints-1]) < JOIN_DIST
1247 : && 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 69586 : else if( fabs(padfX[0]-poOther->padfX[0]) < JOIN_DIST
1268 : && 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 : CPLErr OGRContourWriter( double dfLevel,
1353 : int nPoints, double *padfX, double *padfY,
1354 757 : 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 : poInfo->adfGeoTransform[0]
1376 : + poInfo->adfGeoTransform[1] * padfX[iPoint]
1377 : + poInfo->adfGeoTransform[2] * padfY[iPoint],
1378 : poInfo->adfGeoTransform[3]
1379 : + poInfo->adfGeoTransform[4] * padfX[iPoint]
1380 : + poInfo->adfGeoTransform[5] * padfY[iPoint],
1381 21515 : 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 :
1392 : /************************************************************************/
1393 : /* GDALContourGenerate() */
1394 : /************************************************************************/
1395 :
1396 : /**
1397 : * Create vector contours from raster DEM.
1398 : *
1399 : * This algorithm will generate contours vectors for the input raster band
1400 : * on the requested set of contour levels. The vector contours are written
1401 : * to the passed in OGR vector layer. Also, a NODATA value may be specified
1402 : * to identify pixels that should not be considered in contour line generation.
1403 : *
1404 : * The gdal/apps/gdal_contour.cpp mainline can be used as an example of
1405 : * how to use this function.
1406 : *
1407 : * ALGORITHM RULES
1408 :
1409 : For contouring purposes raster pixel values are assumed to represent a point
1410 : value at the center of the corresponding pixel region. For the purpose of
1411 : contour generation we virtually connect each pixel center to the values to
1412 : the left, right, top and bottom. We assume that the pixel value is linearly
1413 : interpolated between the pixel centers along each line, and determine where
1414 : (if any) contour lines will appear onlong these line segements. Then the
1415 : contour crossings are connected.
1416 :
1417 : This means that contour lines nodes won't actually be on pixel edges, but
1418 : rather along vertical and horizontal lines connecting the pixel centers.
1419 :
1420 : \verbatim
1421 : General Case:
1422 :
1423 : 5 | | 3
1424 : -- + ---------------- + --
1425 : | |
1426 : | |
1427 : | |
1428 : | |
1429 : 10 + |
1430 : |\ |
1431 : | \ |
1432 : -- + -+-------------- + --
1433 : 12 | 10 | 1
1434 :
1435 :
1436 : Saddle Point:
1437 :
1438 : 5 | | 12
1439 : -- + -------------+-- + --
1440 : | \ |
1441 : | \|
1442 : | +
1443 : | |
1444 : + |
1445 : |\ |
1446 : | \ |
1447 : -- + -+-------------- + --
1448 : 12 | | 1
1449 :
1450 : or:
1451 :
1452 : 5 | | 12
1453 : -- + -------------+-- + --
1454 : | __/ |
1455 : | ___/ |
1456 : | ___/ __+
1457 : | / __/ |
1458 : +' __/ |
1459 : | __/ |
1460 : | ,__/ |
1461 : -- + -+-------------- + --
1462 : 12 | | 1
1463 : \endverbatim
1464 :
1465 : Nodata:
1466 :
1467 : In the "nodata" case we treat the whole nodata pixel as a no-mans land.
1468 : We extend the corner pixels near the nodata out to half way and then
1469 : construct extra lines from those points to the center which is assigned
1470 : an averaged value from the two nearby points (in this case (12+3+5)/3).
1471 :
1472 : \verbatim
1473 : 5 | | 3
1474 : -- + ---------------- + --
1475 : | |
1476 : | |
1477 : | 6.7 |
1478 : | +---------+ 3
1479 : 10 +___ |
1480 : | \____+ 10
1481 : | |
1482 : -- + -------+ +
1483 : 12 | 12 (nodata)
1484 :
1485 : \endverbatim
1486 :
1487 : *
1488 : * @param hBand The band to read raster data from. The whole band will be
1489 : * processed.
1490 : *
1491 : * @param dfContourInterval The elevation interval between contours generated.
1492 : *
1493 : * @param dfContourBase The "base" relative to which contour intervals are
1494 : * applied. This is normally zero, but could be different. To generate 10m
1495 : * contours at 5, 15, 25, ... the ContourBase would be 5.
1496 : *
1497 : * @param nFixedLevelCount The number of fixed levels. If this is greater than
1498 : * zero, then fixed levels will be used, and ContourInterval and ContourBase
1499 : * are ignored.
1500 : *
1501 : * @param padfFixedLevels The list of fixed contour levels at which contours
1502 : * should be generated. It will contain FixedLevelCount entries, and may be
1503 : * NULL if fixed levels are disabled (FixedLevelCount = 0).
1504 : *
1505 : * @param bUseNoData If TRUE the dfNoDataValue will be used.
1506 : *
1507 : * @param dfNoDataValue the value to use as a "nodata" value. That is, a
1508 : * pixel value which should be ignored in generating contours as if the value
1509 : * of the pixel were not known.
1510 : *
1511 : * @param hLayer the layer to which new contour vectors will be written.
1512 : * Each contour will have a LINESTRING geometry attached to it. This
1513 : * is really of type OGRLayerH, but void * is used to avoid pulling the
1514 : * ogr_api.h file in here.
1515 : *
1516 : * @param iIDField if not -1 this will be used as a field index to indicate
1517 : * where a unique id should be written for each feature (contour) written.
1518 : *
1519 : * @param iElevField if not -1 this will be used as a field index to indicate
1520 : * where the elevation value of the contour should be written.
1521 : *
1522 : * @param pfnProgress a GDALProgressFunc that may be used to report progress
1523 : * to the user, or to interrupt the algorithm. May be NULL if not required.
1524 : *
1525 : * @param pProgressArg the callback data for the pfnProgress function.
1526 : *
1527 : * @return CE_None on success or CE_Failure if an error occurs.
1528 : */
1529 :
1530 : CPLErr GDALContourGenerate( GDALRasterBandH hBand,
1531 : double dfContourInterval, double dfContourBase,
1532 : int nFixedLevelCount, double *padfFixedLevels,
1533 : int bUseNoData, double dfNoDataValue,
1534 : void *hLayer, int iIDField, int iElevField,
1535 8 : GDALProgressFunc pfnProgress, void *pProgressArg )
1536 :
1537 : {
1538 8 : VALIDATE_POINTER1( hBand, "GDALContourGenerate", CE_Failure );
1539 :
1540 : OGRContourWriterInfo oCWI;
1541 :
1542 8 : if( pfnProgress == NULL )
1543 2 : pfnProgress = GDALDummyProgress;
1544 :
1545 8 : if( !pfnProgress( 0.0, "", pProgressArg ) )
1546 : {
1547 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1548 0 : return CE_Failure;
1549 : }
1550 :
1551 : /* -------------------------------------------------------------------- */
1552 : /* Setup contour writer information. */
1553 : /* -------------------------------------------------------------------- */
1554 : GDALDatasetH hSrcDS;
1555 :
1556 8 : oCWI.hLayer = (OGRLayerH) hLayer;
1557 :
1558 8 : oCWI.nElevField = iElevField;
1559 8 : oCWI.nIDField = iIDField;
1560 :
1561 8 : hSrcDS = GDALGetBandDataset( hBand );
1562 8 : GDALGetGeoTransform( hSrcDS, oCWI.adfGeoTransform );
1563 8 : oCWI.nNextID = 0;
1564 :
1565 : /* -------------------------------------------------------------------- */
1566 : /* Setup contour generator. */
1567 : /* -------------------------------------------------------------------- */
1568 8 : int nXSize = GDALGetRasterBandXSize( hBand );
1569 8 : int nYSize = GDALGetRasterBandYSize( hBand );
1570 :
1571 8 : GDALContourGenerator oCG( nXSize, nYSize, OGRContourWriter, &oCWI );
1572 :
1573 8 : if( nFixedLevelCount > 0 )
1574 2 : oCG.SetFixedLevels( nFixedLevelCount, padfFixedLevels );
1575 : else
1576 6 : oCG.SetContourLevels( dfContourInterval, dfContourBase );
1577 :
1578 8 : if( bUseNoData )
1579 3 : oCG.SetNoData( dfNoDataValue );
1580 :
1581 : /* -------------------------------------------------------------------- */
1582 : /* Feed the data into the contour generator. */
1583 : /* -------------------------------------------------------------------- */
1584 : int iLine;
1585 : double *padfScanline;
1586 8 : CPLErr eErr = CE_None;
1587 :
1588 8 : padfScanline = (double *) VSIMalloc(sizeof(double) * nXSize);
1589 8 : if (padfScanline == NULL)
1590 : {
1591 : CPLError( CE_Failure, CPLE_OutOfMemory,
1592 0 : "VSIMalloc(): Out of memory in GDALContourGenerate" );
1593 8 : return CE_Failure;
1594 : }
1595 :
1596 1055 : for( iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ )
1597 : {
1598 : GDALRasterIO( hBand, GF_Read, 0, iLine, nXSize, 1,
1599 1047 : padfScanline, nXSize, 1, GDT_Float64, 0, 0 );
1600 1047 : eErr = oCG.FeedLine( padfScanline );
1601 :
1602 1047 : if( eErr == CE_None
1603 : && !pfnProgress( (iLine+1) / (double) nYSize, "", pProgressArg ) )
1604 : {
1605 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1606 0 : eErr = CE_Failure;
1607 : }
1608 : }
1609 :
1610 8 : CPLFree( padfScanline );
1611 :
1612 8 : return eErr;
1613 : }
|
