1 : /******************************************************************************
2 : * $Id: gdalgeoloc.cpp 17017 2009-05-15 18:54:38Z rouault $
3 : *
4 : * Project: GDAL
5 : * Purpose: Implements Geolocation array based transformer.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2006, Frank Warmerdam <warmerdam@pobox.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "gdal_priv.h"
31 : #include "gdal_alg.h"
32 :
33 : #ifdef SHAPE_DEBUG
34 : #include "/u/pkg/shapelib/shapefil.h"
35 :
36 : SHPHandle hSHP = NULL;
37 : DBFHandle hDBF = NULL;
38 : #endif
39 :
40 : CPL_CVSID("$Id: gdalgeoloc.cpp 17017 2009-05-15 18:54:38Z rouault $");
41 :
42 : CPL_C_START
43 : CPLXMLNode *GDALSerializeGeoLocTransformer( void *pTransformArg );
44 : void *GDALDeserializeGeoLocTransformer( CPLXMLNode *psTree );
45 : CPL_C_END
46 :
47 : /************************************************************************/
48 : /* ==================================================================== */
49 : /* GDALGeoLocTransformer */
50 : /* ==================================================================== */
51 : /************************************************************************/
52 :
53 : typedef struct {
54 :
55 : GDALTransformerInfo sTI;
56 :
57 : int bReversed;
58 :
59 : // Map from target georef coordinates back to geolocation array
60 : // pixel line coordinates. Built only if needed.
61 :
62 : int nBackMapWidth;
63 : int nBackMapHeight;
64 : double adfBackMapGeoTransform[6]; // maps georef to pixel/line.
65 : float *pafBackMapX;
66 : float *pafBackMapY;
67 :
68 : // geolocation bands.
69 :
70 : GDALDatasetH hDS_X;
71 : GDALRasterBandH hBand_X;
72 : GDALDatasetH hDS_Y;
73 : GDALRasterBandH hBand_Y;
74 :
75 : // Located geolocation data.
76 : int nGeoLocXSize;
77 : int nGeoLocYSize;
78 : double *padfGeoLocX;
79 : double *padfGeoLocY;
80 :
81 : double dfNoDataX;
82 : double dfNoDataY;
83 :
84 : // geolocation <-> base image mapping.
85 : double dfPIXEL_OFFSET;
86 : double dfPIXEL_STEP;
87 : double dfLINE_OFFSET;
88 : double dfLINE_STEP;
89 :
90 : char ** papszGeolocationInfo;
91 :
92 : } GDALGeoLocTransformInfo;
93 :
94 : /************************************************************************/
95 : /* GeoLocLoadFullData() */
96 : /************************************************************************/
97 :
98 2 : static int GeoLocLoadFullData( GDALGeoLocTransformInfo *psTransform )
99 :
100 : {
101 2 : int nXSize = GDALGetRasterXSize( psTransform->hDS_X );
102 2 : int nYSize = GDALGetRasterYSize( psTransform->hDS_X );
103 :
104 2 : psTransform->nGeoLocXSize = nXSize;
105 2 : psTransform->nGeoLocYSize = nYSize;
106 :
107 : psTransform->padfGeoLocY = (double *)
108 2 : VSIMalloc3(sizeof(double), nXSize, nYSize);
109 : psTransform->padfGeoLocX = (double *)
110 2 : VSIMalloc3(sizeof(double), nXSize, nYSize);
111 :
112 2 : if( psTransform->padfGeoLocX == NULL ||
113 : psTransform->padfGeoLocY == NULL )
114 : {
115 : CPLError(CE_Failure, CPLE_OutOfMemory,
116 0 : "GeoLocLoadFullData : Out of memory");
117 0 : return FALSE;
118 : }
119 :
120 2 : if( GDALRasterIO( psTransform->hBand_X, GF_Read,
121 : 0, 0, nXSize, nYSize,
122 : psTransform->padfGeoLocX, nXSize, nYSize,
123 : GDT_Float64, 0, 0 ) != CE_None
124 : || GDALRasterIO( psTransform->hBand_Y, GF_Read,
125 : 0, 0, nXSize, nYSize,
126 : psTransform->padfGeoLocY, nXSize, nYSize,
127 : GDT_Float64, 0, 0 ) != CE_None )
128 0 : return FALSE;
129 :
130 : psTransform->dfNoDataX = GDALGetRasterNoDataValue( psTransform->hBand_X,
131 2 : NULL );
132 : psTransform->dfNoDataY = GDALGetRasterNoDataValue( psTransform->hBand_Y,
133 2 : NULL );
134 :
135 2 : return TRUE;
136 : }
137 :
138 : /************************************************************************/
139 : /* GeoLocGenerateBackMap() */
140 : /************************************************************************/
141 :
142 2 : static int GeoLocGenerateBackMap( GDALGeoLocTransformInfo *psTransform )
143 :
144 : {
145 2 : int nXSize = GDALGetRasterXSize( psTransform->hDS_X );
146 2 : int nYSize = GDALGetRasterYSize( psTransform->hDS_X );
147 2 : int nMaxIter = 3;
148 :
149 : /* -------------------------------------------------------------------- */
150 : /* Scan forward map for lat/long extents. */
151 : /* -------------------------------------------------------------------- */
152 2 : double dfMinX=0, dfMaxX=0, dfMinY=0, dfMaxY=0;
153 2 : int i, bInit = FALSE;
154 :
155 4682 : for( i = nXSize * nYSize - 1; i >= 0; i-- )
156 : {
157 4680 : if( psTransform->padfGeoLocX[i] != psTransform->dfNoDataX )
158 : {
159 4680 : if( bInit )
160 : {
161 4678 : dfMinX = MIN(dfMinX,psTransform->padfGeoLocX[i]);
162 4678 : dfMaxX = MAX(dfMaxX,psTransform->padfGeoLocX[i]);
163 4678 : dfMinY = MIN(dfMinY,psTransform->padfGeoLocY[i]);
164 4678 : dfMaxY = MAX(dfMaxY,psTransform->padfGeoLocY[i]);
165 : }
166 : else
167 : {
168 2 : bInit = TRUE;
169 2 : dfMinX = dfMaxX = psTransform->padfGeoLocX[i];
170 2 : dfMinY = dfMaxY = psTransform->padfGeoLocY[i];
171 : }
172 : }
173 : }
174 :
175 : /* -------------------------------------------------------------------- */
176 : /* Decide on resolution for backmap. We aim for slightly */
177 : /* higher resolution than the source but we can't easily */
178 : /* establish how much dead space there is in the backmap, so it */
179 : /* is approximate. */
180 : /* -------------------------------------------------------------------- */
181 2 : double dfTargetPixels = (nXSize * nYSize * 1.3);
182 : double dfPixelSize = sqrt((dfMaxX - dfMinX) * (dfMaxY - dfMinY)
183 2 : / dfTargetPixels);
184 : int nBMXSize, nBMYSize;
185 :
186 : nBMYSize = psTransform->nBackMapHeight =
187 2 : (int) ((dfMaxY - dfMinY) / dfPixelSize + 1);
188 : nBMXSize= psTransform->nBackMapWidth =
189 2 : (int) ((dfMaxX - dfMinX) / dfPixelSize + 1);
190 :
191 2 : if (nBMXSize > INT_MAX / nBMYSize)
192 : {
193 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d",
194 0 : nBMXSize, nBMYSize);
195 0 : return FALSE;
196 : }
197 :
198 2 : dfMinX -= dfPixelSize/2.0;
199 2 : dfMaxY += dfPixelSize/2.0;
200 :
201 2 : psTransform->adfBackMapGeoTransform[0] = dfMinX;
202 2 : psTransform->adfBackMapGeoTransform[1] = dfPixelSize;
203 2 : psTransform->adfBackMapGeoTransform[2] = 0.0;
204 2 : psTransform->adfBackMapGeoTransform[3] = dfMaxY;
205 2 : psTransform->adfBackMapGeoTransform[4] = 0.0;
206 2 : psTransform->adfBackMapGeoTransform[5] = -dfPixelSize;
207 :
208 : /* -------------------------------------------------------------------- */
209 : /* Allocate backmap, and initialize to nodata value (-1.0). */
210 : /* -------------------------------------------------------------------- */
211 : GByte *pabyValidFlag;
212 :
213 : pabyValidFlag = (GByte *)
214 2 : VSICalloc(nBMXSize, nBMYSize);
215 :
216 : psTransform->pafBackMapX = (float *)
217 2 : VSIMalloc3(nBMXSize, nBMYSize, sizeof(float));
218 : psTransform->pafBackMapY = (float *)
219 2 : VSIMalloc3(nBMXSize, nBMYSize, sizeof(float));
220 :
221 2 : if( pabyValidFlag == NULL ||
222 : psTransform->pafBackMapX == NULL ||
223 : psTransform->pafBackMapY == NULL )
224 : {
225 : CPLError( CE_Failure, CPLE_OutOfMemory,
226 : "Unable to allocate %dx%d back-map for geolocation array transformer.",
227 0 : nBMXSize, nBMYSize );
228 0 : CPLFree( pabyValidFlag );
229 0 : return FALSE;
230 : }
231 :
232 6258 : for( i = nBMXSize * nBMYSize - 1; i >= 0; i-- )
233 : {
234 6256 : psTransform->pafBackMapX[i] = -1.0;
235 6256 : psTransform->pafBackMapY[i] = -1.0;
236 : }
237 :
238 : /* -------------------------------------------------------------------- */
239 : /* Run through the whole geoloc array forward projecting and */
240 : /* pushing into the backmap. */
241 : /* Initialise to the nMaxIter+1 value so we can spot genuinely */
242 : /* valid pixels in the hole-filling loop. */
243 : /* -------------------------------------------------------------------- */
244 : int iBMX, iBMY;
245 : int iX, iY;
246 :
247 80 : for( iY = 0; iY < nYSize; iY++ )
248 : {
249 4758 : for( iX = 0; iX < nXSize; iX++ )
250 : {
251 4680 : if( psTransform->padfGeoLocX[iX + iY * nXSize]
252 : == psTransform->dfNoDataX )
253 0 : continue;
254 :
255 4680 : i = iX + iY * nXSize;
256 :
257 4680 : iBMX = (int) ((psTransform->padfGeoLocX[i] - dfMinX) / dfPixelSize);
258 4680 : iBMY = (int) ((dfMaxY - psTransform->padfGeoLocY[i]) / dfPixelSize);
259 :
260 4680 : if( iBMX < 0 || iBMY < 0 || iBMX >= nBMXSize || iBMY >= nBMYSize )
261 2 : continue;
262 :
263 4678 : psTransform->pafBackMapX[iBMX + iBMY * nBMXSize] =
264 4678 : (float)(iX * psTransform->dfPIXEL_STEP + psTransform->dfPIXEL_OFFSET);
265 4678 : psTransform->pafBackMapY[iBMX + iBMY * nBMXSize] =
266 4678 : (float)(iY * psTransform->dfLINE_STEP + psTransform->dfLINE_OFFSET);
267 :
268 4678 : pabyValidFlag[iBMX + iBMY * nBMXSize] = nMaxIter+1;
269 :
270 : }
271 : }
272 :
273 : /* -------------------------------------------------------------------- */
274 : /* Now, loop over the backmap trying to fill in holes with */
275 : /* nearby values. */
276 : /* -------------------------------------------------------------------- */
277 : int iIter;
278 : int nNumValid, nExpectedValid;
279 :
280 8 : for( iIter = 0; iIter < nMaxIter; iIter++ )
281 : {
282 6 : nNumValid = 0;
283 6 : nExpectedValid = (nBMYSize - 2) * (nBMXSize - 2);
284 270 : for( iBMY = 1; iBMY < nBMYSize-1; iBMY++ )
285 : {
286 17688 : for( iBMX = 1; iBMX < nBMXSize-1; iBMX++ )
287 : {
288 : // if this point is already set, ignore it.
289 17424 : if( pabyValidFlag[iBMX + iBMY*nBMXSize] )
290 : {
291 12860 : nNumValid++;
292 12860 : continue;
293 : }
294 :
295 4564 : int nCount = 0;
296 4564 : double dfXSum = 0.0, dfYSum = 0.0;
297 4564 : int nMarkedAsGood = nMaxIter - iIter;
298 :
299 : // left?
300 4564 : if( pabyValidFlag[iBMX-1+iBMY*nBMXSize] > nMarkedAsGood )
301 : {
302 788 : dfXSum += psTransform->pafBackMapX[iBMX-1+iBMY*nBMXSize];
303 788 : dfYSum += psTransform->pafBackMapY[iBMX-1+iBMY*nBMXSize];
304 788 : nCount++;
305 : }
306 : // right?
307 4564 : if( pabyValidFlag[iBMX+1+iBMY*nBMXSize] > nMarkedAsGood )
308 : {
309 812 : dfXSum += psTransform->pafBackMapX[iBMX+1+iBMY*nBMXSize];
310 812 : dfYSum += psTransform->pafBackMapY[iBMX+1+iBMY*nBMXSize];
311 812 : nCount++;
312 : }
313 : // top?
314 4564 : if( pabyValidFlag[iBMX+(iBMY-1)*nBMXSize] > nMarkedAsGood )
315 : {
316 756 : dfXSum += psTransform->pafBackMapX[iBMX+(iBMY-1)*nBMXSize];
317 756 : dfYSum += psTransform->pafBackMapY[iBMX+(iBMY-1)*nBMXSize];
318 756 : nCount++;
319 : }
320 : // bottom?
321 4564 : if( pabyValidFlag[iBMX+(iBMY+1)*nBMXSize] > nMarkedAsGood )
322 : {
323 736 : dfXSum += psTransform->pafBackMapX[iBMX+(iBMY+1)*nBMXSize];
324 736 : dfYSum += psTransform->pafBackMapY[iBMX+(iBMY+1)*nBMXSize];
325 736 : nCount++;
326 : }
327 : // top-left?
328 4564 : if( pabyValidFlag[iBMX-1+(iBMY-1)*nBMXSize] > nMarkedAsGood )
329 : {
330 908 : dfXSum += psTransform->pafBackMapX[iBMX-1+(iBMY-1)*nBMXSize];
331 908 : dfYSum += psTransform->pafBackMapY[iBMX-1+(iBMY-1)*nBMXSize];
332 908 : nCount++;
333 : }
334 : // top-right?
335 4564 : if( pabyValidFlag[iBMX+1+(iBMY-1)*nBMXSize] > nMarkedAsGood )
336 : {
337 778 : dfXSum += psTransform->pafBackMapX[iBMX+1+(iBMY-1)*nBMXSize];
338 778 : dfYSum += psTransform->pafBackMapY[iBMX+1+(iBMY-1)*nBMXSize];
339 778 : nCount++;
340 : }
341 : // bottom-left?
342 4564 : if( pabyValidFlag[iBMX-1+(iBMY+1)*nBMXSize] > nMarkedAsGood )
343 : {
344 734 : dfXSum += psTransform->pafBackMapX[iBMX-1+(iBMY+1)*nBMXSize];
345 734 : dfYSum += psTransform->pafBackMapY[iBMX-1+(iBMY+1)*nBMXSize];
346 734 : nCount++;
347 : }
348 : // bottom-right?
349 4564 : if( pabyValidFlag[iBMX+1+(iBMY+1)*nBMXSize] > nMarkedAsGood )
350 : {
351 916 : dfXSum += psTransform->pafBackMapX[iBMX+1+(iBMY+1)*nBMXSize];
352 916 : dfYSum += psTransform->pafBackMapY[iBMX+1+(iBMY+1)*nBMXSize];
353 916 : nCount++;
354 : }
355 :
356 4564 : if( nCount > 0 )
357 : {
358 1702 : psTransform->pafBackMapX[iBMX + iBMY * nBMXSize] = (float)(dfXSum/nCount);
359 1702 : psTransform->pafBackMapY[iBMX + iBMY * nBMXSize] = (float)(dfYSum/nCount);
360 : // genuinely valid points will have value iMaxIter+1
361 : // On each iteration mark newly valid points with a
362 : // descending value so that it will not be used on the
363 : // current iteration only on subsequent ones.
364 1702 : pabyValidFlag[iBMX+iBMY*nBMXSize] = nMaxIter - iIter;
365 : }
366 : }
367 : }
368 6 : if (nNumValid == nExpectedValid)
369 0 : break;
370 : }
371 :
372 : #ifdef notdef
373 : GDALDatasetH hBMDS = GDALCreate( GDALGetDriverByName( "GTiff" ),
374 : "backmap.tif", nBMXSize, nBMYSize, 2,
375 : GDT_Float32, NULL );
376 : GDALSetGeoTransform( hBMDS, psTransform->adfBackMapGeoTransform );
377 : GDALRasterIO( GDALGetRasterBand(hBMDS,1), GF_Write,
378 : 0, 0, nBMXSize, nBMYSize,
379 : psTransform->pafBackMapX, nBMXSize, nBMYSize,
380 : GDT_Float32, 0, 0 );
381 : GDALRasterIO( GDALGetRasterBand(hBMDS,2), GF_Write,
382 : 0, 0, nBMXSize, nBMYSize,
383 : psTransform->pafBackMapY, nBMXSize, nBMYSize,
384 : GDT_Float32, 0, 0 );
385 : GDALClose( hBMDS );
386 : #endif
387 :
388 2 : CPLFree( pabyValidFlag );
389 :
390 2 : return TRUE;
391 : }
392 :
393 : /************************************************************************/
394 : /* FindGeoLocPosition() */
395 : /************************************************************************/
396 :
397 : #ifdef notdef
398 :
399 : /*
400 : This searching approach has been abandoned because it is too sensitive
401 : to discontinuities in the data. Left in case it might be revived in
402 : the future.
403 : */
404 :
405 : static int FindGeoLocPosition( GDALGeoLocTransformInfo *psTransform,
406 : double dfGeoX, double dfGeoY,
407 : int nStartX, int nStartY,
408 : double *pdfFoundX, double *pdfFoundY )
409 :
410 : {
411 : double adfPathX[5000], adfPathY[5000];
412 :
413 : if( psTransform->padfGeoLocX == NULL )
414 : return FALSE;
415 :
416 : int nXSize = psTransform->nGeoLocXSize;
417 : int nYSize = psTransform->nGeoLocYSize;
418 : int nStepCount = 0;
419 :
420 : // Start in center if we don't have any provided info.
421 : if( nStartX < 0 || nStartY < 0
422 : || nStartX >= nXSize || nStartY >= nYSize )
423 : {
424 : nStartX = nXSize / 2;
425 : nStartY = nYSize / 2;
426 : }
427 :
428 : nStartX = MIN(nStartX,nXSize-2);
429 : nStartY = MIN(nStartY,nYSize-2);
430 :
431 : int iX = nStartX, iY = nStartY;
432 : int iLastX = -1, iLastY = -1;
433 : int iSecondLastX = -1, iSecondLastY = -1;
434 :
435 : while( nStepCount < MAX(nXSize,nYSize) )
436 : {
437 : int iXNext = -1, iYNext = -1;
438 : double dfDeltaXRight, dfDeltaYRight, dfDeltaXDown, dfDeltaYDown;
439 :
440 : double *padfThisX = psTransform->padfGeoLocX + iX + iY * nXSize;
441 : double *padfThisY = psTransform->padfGeoLocY + iX + iY * nXSize;
442 :
443 : double dfDeltaX = dfGeoX - *padfThisX;
444 : double dfDeltaY = dfGeoY - *padfThisY;
445 :
446 : if( iX == nXSize-1 )
447 : {
448 : dfDeltaXRight = *(padfThisX) - *(padfThisX-1);
449 : dfDeltaYRight = *(padfThisY) - *(padfThisY-1);
450 : }
451 : else
452 : {
453 : dfDeltaXRight = *(padfThisX+1) - *padfThisX;
454 : dfDeltaYRight = *(padfThisY+1) - *padfThisY;
455 : }
456 :
457 : if( iY == nYSize - 1 )
458 : {
459 : dfDeltaXDown = *(padfThisX) - *(padfThisX-nXSize);
460 : dfDeltaYDown = *(padfThisY) - *(padfThisY-nXSize);
461 : }
462 : else
463 : {
464 : dfDeltaXDown = *(padfThisX+nXSize) - *padfThisX;
465 : dfDeltaYDown = *(padfThisY+nXSize) - *padfThisY;
466 : }
467 :
468 : double dfRightProjection =
469 : (dfDeltaXRight * dfDeltaX + dfDeltaYRight * dfDeltaY)
470 : / (dfDeltaXRight*dfDeltaXRight + dfDeltaYRight*dfDeltaYRight);
471 :
472 : double dfDownProjection =
473 : (dfDeltaXDown * dfDeltaX + dfDeltaYDown * dfDeltaY)
474 : / (dfDeltaXDown*dfDeltaXDown + dfDeltaYDown*dfDeltaYDown);
475 :
476 : // Are we in our target cell?
477 : if( dfRightProjection >= 0.0 && dfRightProjection < 1.0
478 : && dfDownProjection >= 0.0 && dfDownProjection < 1.0 )
479 : {
480 : *pdfFoundX = iX + dfRightProjection;
481 : *pdfFoundY = iY + dfDownProjection;
482 :
483 : return TRUE;
484 : }
485 :
486 : if( ABS(dfRightProjection) > ABS(dfDownProjection) )
487 : {
488 : // Do we want to move right?
489 : if( dfRightProjection > 1.0 && iX < nXSize-1 )
490 : {
491 : iXNext = iX + MAX(1,(int)(dfRightProjection - nStepCount)/2);
492 : iYNext = iY;
493 : }
494 :
495 : // Do we want to move left?
496 : else if( dfRightProjection < 0.0 && iX > 0 )
497 : {
498 : iXNext = iX - MAX(1,(int)(ABS(dfRightProjection) - nStepCount)/2);
499 : iYNext = iY;
500 : }
501 :
502 : // Do we want to move down.
503 : else if( dfDownProjection > 1.0 && iY < nYSize-1 )
504 : {
505 : iXNext = iX;
506 : iYNext = iY + MAX(1,(int)(dfDownProjection - nStepCount)/2);
507 : }
508 :
509 : // Do we want to move up?
510 : else if( dfDownProjection < 0.0 && iY > 0 )
511 : {
512 : iXNext = iX;
513 : iYNext = iY - MAX(1,(int)(ABS(dfDownProjection) - nStepCount)/2);
514 : }
515 :
516 : // We aren't there, and we have no where to go
517 : else
518 : {
519 : return FALSE;
520 : }
521 : }
522 : else
523 : {
524 : // Do we want to move down.
525 : if( dfDownProjection > 1.0 && iY < nYSize-1 )
526 : {
527 : iXNext = iX;
528 : iYNext = iY + MAX(1,(int)(dfDownProjection - nStepCount)/2);
529 : }
530 :
531 : // Do we want to move up?
532 : else if( dfDownProjection < 0.0 && iY > 0 )
533 : {
534 : iXNext = iX;
535 : iYNext = iY - MAX(1,(int)(ABS(dfDownProjection) - nStepCount)/2);
536 : }
537 :
538 : // Do we want to move right?
539 : else if( dfRightProjection > 1.0 && iX < nXSize-1 )
540 : {
541 : iXNext = iX + MAX(1,(int)(dfRightProjection - nStepCount)/2);
542 : iYNext = iY;
543 : }
544 :
545 : // Do we want to move left?
546 : else if( dfRightProjection < 0.0 && iX > 0 )
547 : {
548 : iXNext = iX - MAX(1,(int)(ABS(dfRightProjection) - nStepCount)/2);
549 : iYNext = iY;
550 : }
551 :
552 : // We aren't there, and we have no where to go
553 : else
554 : {
555 : return FALSE;
556 : }
557 : }
558 : adfPathX[nStepCount] = iX;
559 : adfPathY[nStepCount] = iY;
560 :
561 : nStepCount++;
562 : iX = MAX(0,MIN(iXNext,nXSize-1));
563 : iY = MAX(0,MIN(iYNext,nYSize-1));
564 :
565 : if( iX == iSecondLastX && iY == iSecondLastY )
566 : {
567 : // Are we *near* our target cell?
568 : if( dfRightProjection >= -1.0 && dfRightProjection < 2.0
569 : && dfDownProjection >= -1.0 && dfDownProjection < 2.0 )
570 : {
571 : *pdfFoundX = iX + dfRightProjection;
572 : *pdfFoundY = iY + dfDownProjection;
573 :
574 : return TRUE;
575 : }
576 :
577 : #ifdef SHAPE_DEBUG
578 : if( hSHP != NULL )
579 : {
580 : SHPObject *hObj;
581 :
582 : hObj = SHPCreateSimpleObject( SHPT_ARC, nStepCount,
583 : adfPathX, adfPathY, NULL );
584 : SHPWriteObject( hSHP, -1, hObj );
585 : SHPDestroyObject( hObj );
586 :
587 : int iShape = DBFGetRecordCount( hDBF );
588 : DBFWriteDoubleAttribute( hDBF, iShape, 0, dfGeoX );
589 : DBFWriteDoubleAttribute( hDBF, iShape, 1, dfGeoY );
590 : }
591 : #endif
592 : //CPLDebug( "GeoL", "Looping at step (%d) on search for %g,%g.",
593 : // nStepCount, dfGeoX, dfGeoY );
594 : return FALSE;
595 : }
596 :
597 : iSecondLastX = iLastX;
598 : iSecondLastY = iLastY;
599 :
600 : iLastX = iX;
601 : iLastY = iY;
602 :
603 : }
604 :
605 : //CPLDebug( "GeoL", "Exceeded step count max (%d) on search for %g,%g.",
606 : // MAX(nXSize,nYSize),
607 : // dfGeoX, dfGeoY );
608 :
609 : #ifdef SHAPE_DEBUG
610 : if( hSHP != NULL )
611 : {
612 : SHPObject *hObj;
613 :
614 : hObj = SHPCreateSimpleObject( SHPT_ARC, nStepCount,
615 : adfPathX, adfPathY, NULL );
616 : SHPWriteObject( hSHP, -1, hObj );
617 : SHPDestroyObject( hObj );
618 :
619 : int iShape = DBFGetRecordCount( hDBF );
620 : DBFWriteDoubleAttribute( hDBF, iShape, 0, dfGeoX );
621 : DBFWriteDoubleAttribute( hDBF, iShape, 1, dfGeoY );
622 : }
623 : #endif
624 :
625 : return FALSE;
626 : }
627 : #endif /* def notdef */
628 :
629 :
630 :
631 : /************************************************************************/
632 : /* GDALCreateGeoLocTransformer() */
633 : /************************************************************************/
634 :
635 2 : void *GDALCreateGeoLocTransformer( GDALDatasetH hBaseDS,
636 : char **papszGeolocationInfo,
637 : int bReversed )
638 :
639 : {
640 : GDALGeoLocTransformInfo *psTransform;
641 :
642 2 : if( CSLFetchNameValue(papszGeolocationInfo,"PIXEL_OFFSET") == NULL
643 : || CSLFetchNameValue(papszGeolocationInfo,"LINE_OFFSET") == NULL
644 : || CSLFetchNameValue(papszGeolocationInfo,"PIXEL_STEP") == NULL
645 : || CSLFetchNameValue(papszGeolocationInfo,"LINE_STEP") == NULL
646 : || CSLFetchNameValue(papszGeolocationInfo,"X_BAND") == NULL
647 : || CSLFetchNameValue(papszGeolocationInfo,"Y_BAND") == NULL )
648 : {
649 : CPLError( CE_Failure, CPLE_AppDefined,
650 0 : "Missing some geolocation fields in GDALCreateGeoLocTransformer()" );
651 0 : return NULL;
652 : }
653 :
654 : /* -------------------------------------------------------------------- */
655 : /* Initialize core info. */
656 : /* -------------------------------------------------------------------- */
657 : psTransform = (GDALGeoLocTransformInfo *)
658 2 : CPLCalloc(sizeof(GDALGeoLocTransformInfo),1);
659 :
660 2 : psTransform->bReversed = bReversed;
661 :
662 2 : strcpy( psTransform->sTI.szSignature, "GTI" );
663 2 : psTransform->sTI.pszClassName = "GDALGeoLocTransformer";
664 2 : psTransform->sTI.pfnTransform = GDALGeoLocTransform;
665 2 : psTransform->sTI.pfnCleanup = GDALDestroyGeoLocTransformer;
666 2 : psTransform->sTI.pfnSerialize = GDALSerializeGeoLocTransformer;
667 2 : psTransform->papszGeolocationInfo = CSLDuplicate( papszGeolocationInfo );
668 :
669 : /* -------------------------------------------------------------------- */
670 : /* Pull geolocation info from the options/metadata. */
671 : /* -------------------------------------------------------------------- */
672 : psTransform->dfPIXEL_OFFSET = atof(CSLFetchNameValue( papszGeolocationInfo,
673 2 : "PIXEL_OFFSET" ));
674 : psTransform->dfLINE_OFFSET = atof(CSLFetchNameValue( papszGeolocationInfo,
675 2 : "LINE_OFFSET" ));
676 : psTransform->dfPIXEL_STEP = atof(CSLFetchNameValue( papszGeolocationInfo,
677 2 : "PIXEL_STEP" ));
678 : psTransform->dfLINE_STEP = atof(CSLFetchNameValue( papszGeolocationInfo,
679 2 : "LINE_STEP" ));
680 :
681 : /* -------------------------------------------------------------------- */
682 : /* Establish access to geolocation dataset(s). */
683 : /* -------------------------------------------------------------------- */
684 : const char *pszDSName = CSLFetchNameValue( papszGeolocationInfo,
685 2 : "X_DATASET" );
686 2 : if( pszDSName != NULL )
687 : {
688 2 : psTransform->hDS_X = GDALOpenShared( pszDSName, GA_ReadOnly );
689 : }
690 : else
691 : {
692 0 : psTransform->hDS_X = hBaseDS;
693 0 : GDALReferenceDataset( psTransform->hDS_X );
694 : psTransform->papszGeolocationInfo =
695 : CSLSetNameValue( psTransform->papszGeolocationInfo,
696 : "X_DATASET",
697 0 : GDALGetDescription( hBaseDS ) );
698 : }
699 :
700 2 : pszDSName = CSLFetchNameValue( papszGeolocationInfo, "Y_DATASET" );
701 2 : if( pszDSName != NULL )
702 : {
703 2 : psTransform->hDS_Y = GDALOpenShared( pszDSName, GA_ReadOnly );
704 : }
705 : else
706 : {
707 0 : psTransform->hDS_Y = hBaseDS;
708 0 : GDALReferenceDataset( psTransform->hDS_Y );
709 : psTransform->papszGeolocationInfo =
710 : CSLSetNameValue( psTransform->papszGeolocationInfo,
711 : "Y_DATASET",
712 0 : GDALGetDescription( hBaseDS ) );
713 : }
714 :
715 2 : if (psTransform->hDS_X == NULL ||
716 : psTransform->hDS_Y == NULL)
717 : {
718 0 : GDALDestroyGeoLocTransformer( psTransform );
719 0 : return NULL;
720 : }
721 :
722 : /* -------------------------------------------------------------------- */
723 : /* Get the band handles. */
724 : /* -------------------------------------------------------------------- */
725 : int nBand;
726 :
727 2 : nBand = MAX(1,atoi(CSLFetchNameValue( papszGeolocationInfo, "X_BAND" )));
728 2 : psTransform->hBand_X = GDALGetRasterBand( psTransform->hDS_X, nBand );
729 :
730 2 : nBand = MAX(1,atoi(CSLFetchNameValue( papszGeolocationInfo, "Y_BAND" )));
731 2 : psTransform->hBand_Y = GDALGetRasterBand( psTransform->hDS_Y, nBand );
732 :
733 2 : if (psTransform->hBand_X == NULL ||
734 : psTransform->hBand_Y == NULL)
735 : {
736 0 : GDALDestroyGeoLocTransformer( psTransform );
737 0 : return NULL;
738 : }
739 :
740 : /* -------------------------------------------------------------------- */
741 : /* Check that X and Y bands have the same dimensions */
742 : /* -------------------------------------------------------------------- */
743 2 : int nXSize_XBand = GDALGetRasterXSize( psTransform->hDS_X );
744 2 : int nYSize_XBand = GDALGetRasterYSize( psTransform->hDS_X );
745 2 : int nXSize_YBand = GDALGetRasterXSize( psTransform->hDS_Y );
746 2 : int nYSize_YBand = GDALGetRasterYSize( psTransform->hDS_Y );
747 2 : if (nXSize_XBand != nXSize_YBand ||
748 : nYSize_XBand != nYSize_YBand )
749 : {
750 : CPLError(CE_Failure, CPLE_AppDefined,
751 0 : "X_BAND and Y_BAND do not have the same dimensions");
752 0 : GDALDestroyGeoLocTransformer( psTransform );
753 0 : return NULL;
754 : }
755 :
756 2 : if (nXSize_XBand > INT_MAX / nYSize_XBand)
757 : {
758 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d",
759 0 : nXSize_XBand, nYSize_XBand);
760 0 : GDALDestroyGeoLocTransformer( psTransform );
761 0 : return NULL;
762 : }
763 :
764 : /* -------------------------------------------------------------------- */
765 : /* Load the geolocation array. */
766 : /* -------------------------------------------------------------------- */
767 2 : if( !GeoLocLoadFullData( psTransform )
768 : || !GeoLocGenerateBackMap( psTransform ) )
769 : {
770 0 : GDALDestroyGeoLocTransformer( psTransform );
771 0 : return NULL;
772 : }
773 :
774 2 : return psTransform;
775 : }
776 :
777 : /************************************************************************/
778 : /* GDALDestroyGeoLocTransformer() */
779 : /************************************************************************/
780 :
781 2 : void GDALDestroyGeoLocTransformer( void *pTransformAlg )
782 :
783 : {
784 : GDALGeoLocTransformInfo *psTransform =
785 2 : (GDALGeoLocTransformInfo *) pTransformAlg;
786 :
787 2 : CPLFree( psTransform->pafBackMapX );
788 2 : CPLFree( psTransform->pafBackMapY );
789 2 : CSLDestroy( psTransform->papszGeolocationInfo );
790 2 : CPLFree( psTransform->padfGeoLocX );
791 2 : CPLFree( psTransform->padfGeoLocY );
792 :
793 2 : if( psTransform->hDS_X != NULL
794 : && GDALDereferenceDataset( psTransform->hDS_X ) == 0 )
795 0 : GDALClose( psTransform->hDS_X );
796 :
797 2 : if( psTransform->hDS_Y != NULL
798 : && GDALDereferenceDataset( psTransform->hDS_Y ) == 0 )
799 2 : GDALClose( psTransform->hDS_Y );
800 :
801 2 : CPLFree( pTransformAlg );
802 2 : }
803 :
804 : /************************************************************************/
805 : /* GDALGeoLocTransform() */
806 : /************************************************************************/
807 :
808 260 : int GDALGeoLocTransform( void *pTransformArg, int bDstToSrc,
809 : int nPointCount,
810 : double *padfX, double *padfY, double *padfZ,
811 : int *panSuccess )
812 :
813 : {
814 : GDALGeoLocTransformInfo *psTransform =
815 260 : (GDALGeoLocTransformInfo *) pTransformArg;
816 :
817 260 : if( psTransform->bReversed )
818 0 : bDstToSrc = !bDstToSrc;
819 :
820 : /* -------------------------------------------------------------------- */
821 : /* Do original pixel line to target geox/geoy. */
822 : /* -------------------------------------------------------------------- */
823 260 : if( !bDstToSrc )
824 : {
825 1 : int i, nXSize = psTransform->nGeoLocXSize;
826 :
827 2 : for( i = 0; i < nPointCount; i++ )
828 : {
829 1 : if( !panSuccess[i] )
830 0 : continue;
831 :
832 1 : if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL )
833 : {
834 0 : panSuccess[i] = FALSE;
835 0 : continue;
836 : }
837 :
838 1 : double dfGeoLocPixel = (padfX[i] - psTransform->dfPIXEL_OFFSET)
839 1 : / psTransform->dfPIXEL_STEP;
840 1 : double dfGeoLocLine = (padfY[i] - psTransform->dfLINE_OFFSET)
841 1 : / psTransform->dfLINE_STEP;
842 :
843 : int iX, iY;
844 :
845 1 : iX = MAX(0,(int) dfGeoLocPixel);
846 1 : iX = MIN(iX,psTransform->nGeoLocXSize-2);
847 1 : iY = MAX(0,(int) dfGeoLocLine);
848 1 : iY = MIN(iY,psTransform->nGeoLocYSize-2);
849 :
850 1 : double *padfGLX = psTransform->padfGeoLocX + iX + iY * nXSize;
851 1 : double *padfGLY = psTransform->padfGeoLocY + iX + iY * nXSize;
852 :
853 : // This assumes infinite extension beyond borders of available
854 : // data based on closest grid square.
855 :
856 1 : padfX[i] = padfGLX[0]
857 1 : + (dfGeoLocPixel-iX) * (padfGLX[1] - padfGLX[0])
858 2 : + (dfGeoLocLine -iY) * (padfGLX[nXSize] - padfGLX[0]);
859 1 : padfY[i] = padfGLY[0]
860 1 : + (dfGeoLocPixel-iX) * (padfGLY[1] - padfGLY[0])
861 2 : + (dfGeoLocLine -iY) * (padfGLY[nXSize] - padfGLY[0]);
862 :
863 1 : panSuccess[i] = TRUE;
864 : }
865 : }
866 :
867 : /* -------------------------------------------------------------------- */
868 : /* geox/geoy to pixel/line using backmap. */
869 : /* -------------------------------------------------------------------- */
870 : else
871 : {
872 : int i;
873 :
874 66705 : for( i = 0; i < nPointCount; i++ )
875 : {
876 66446 : if( !panSuccess[i] )
877 0 : continue;
878 :
879 66446 : if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL )
880 : {
881 0 : panSuccess[i] = FALSE;
882 0 : continue;
883 : }
884 :
885 : int iBMX, iBMY;
886 :
887 66446 : iBMX = (int) ((padfX[i] - psTransform->adfBackMapGeoTransform[0])
888 66446 : / psTransform->adfBackMapGeoTransform[1]);
889 66446 : iBMY = (int) ((padfY[i] - psTransform->adfBackMapGeoTransform[3])
890 66446 : / psTransform->adfBackMapGeoTransform[5]);
891 :
892 66446 : int iBM = iBMX + iBMY * psTransform->nBackMapWidth;
893 :
894 69480 : if( iBMX < 0 || iBMY < 0
895 : || iBMX >= psTransform->nBackMapWidth
896 : || iBMY >= psTransform->nBackMapHeight
897 3034 : || psTransform->pafBackMapX[iBM] < 0 )
898 : {
899 63969 : panSuccess[i] = FALSE;
900 63969 : padfX[i] = HUGE_VAL;
901 63969 : padfY[i] = HUGE_VAL;
902 63969 : continue;
903 : }
904 :
905 2477 : padfX[i] = psTransform->pafBackMapX[iBM];
906 2477 : padfY[i] = psTransform->pafBackMapY[iBM];
907 2477 : panSuccess[i] = TRUE;
908 : }
909 : }
910 :
911 : /* -------------------------------------------------------------------- */
912 : /* geox/geoy to pixel/line using search algorithm. */
913 : /* -------------------------------------------------------------------- */
914 : #ifdef notdef
915 : else
916 : {
917 : int i;
918 : int nStartX = -1, nStartY = -1;
919 :
920 : #ifdef SHAPE_DEBUG
921 : hSHP = SHPCreate( "tracks.shp", SHPT_ARC );
922 : hDBF = DBFCreate( "tracks.dbf" );
923 : DBFAddField( hDBF, "GEOX", FTDouble, 10, 4 );
924 : DBFAddField( hDBF, "GEOY", FTDouble, 10, 4 );
925 : #endif
926 : for( i = 0; i < nPointCount; i++ )
927 : {
928 : double dfGeoLocX, dfGeoLocY;
929 :
930 : if( !panSuccess[i] )
931 : continue;
932 :
933 : if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL )
934 : {
935 : panSuccess[i] = FALSE;
936 : continue;
937 : }
938 :
939 : if( !FindGeoLocPosition( psTransform, padfX[i], padfY[i],
940 : -1, -1, &dfGeoLocX, &dfGeoLocY ) )
941 : {
942 : padfX[i] = HUGE_VAL;
943 : padfY[i] = HUGE_VAL;
944 :
945 : panSuccess[i] = FALSE;
946 : continue;
947 : }
948 : nStartX = (int) dfGeoLocX;
949 : nStartY = (int) dfGeoLocY;
950 :
951 : padfX[i] = dfGeoLocX * psTransform->dfPIXEL_STEP
952 : + psTransform->dfPIXEL_OFFSET;
953 : padfY[i] = dfGeoLocY * psTransform->dfLINE_STEP
954 : + psTransform->dfLINE_OFFSET;
955 :
956 : panSuccess[i] = TRUE;
957 : }
958 :
959 : #ifdef SHAPE_DEBUG
960 : if( hSHP != NULL )
961 : {
962 : DBFClose( hDBF );
963 : hDBF = NULL;
964 :
965 : SHPClose( hSHP );
966 : hSHP = NULL;
967 : }
968 : #endif
969 : }
970 : #endif
971 :
972 260 : return TRUE;
973 : }
974 :
975 : /************************************************************************/
976 : /* GDALSerializeGeoLocTransformer() */
977 : /************************************************************************/
978 :
979 0 : CPLXMLNode *GDALSerializeGeoLocTransformer( void *pTransformArg )
980 :
981 : {
982 0 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeGeoLocTransformer", NULL );
983 :
984 : CPLXMLNode *psTree;
985 : GDALGeoLocTransformInfo *psInfo =
986 0 : (GDALGeoLocTransformInfo *)(pTransformArg);
987 :
988 0 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "GeoLocTransformer" );
989 :
990 : /* -------------------------------------------------------------------- */
991 : /* Serialize bReversed. */
992 : /* -------------------------------------------------------------------- */
993 : CPLCreateXMLElementAndValue(
994 : psTree, "Reversed",
995 0 : CPLString().Printf( "%d", psInfo->bReversed ) );
996 :
997 : /* -------------------------------------------------------------------- */
998 : /* geoloc metadata. */
999 : /* -------------------------------------------------------------------- */
1000 0 : char **papszMD = psInfo->papszGeolocationInfo;
1001 : CPLXMLNode *psMD= CPLCreateXMLNode( psTree, CXT_Element,
1002 0 : "Metadata" );
1003 :
1004 0 : for( int i = 0; papszMD != NULL && papszMD[i] != NULL; i++ )
1005 : {
1006 : const char *pszRawValue;
1007 : char *pszKey;
1008 : CPLXMLNode *psMDI;
1009 :
1010 0 : pszRawValue = CPLParseNameValue( papszMD[i], &pszKey );
1011 :
1012 0 : psMDI = CPLCreateXMLNode( psMD, CXT_Element, "MDI" );
1013 0 : CPLSetXMLValue( psMDI, "#key", pszKey );
1014 0 : CPLCreateXMLNode( psMDI, CXT_Text, pszRawValue );
1015 :
1016 0 : CPLFree( pszKey );
1017 : }
1018 :
1019 0 : return psTree;
1020 : }
1021 :
1022 : /************************************************************************/
1023 : /* GDALDeserializeGeoLocTransformer() */
1024 : /************************************************************************/
1025 :
1026 1 : void *GDALDeserializeGeoLocTransformer( CPLXMLNode *psTree )
1027 :
1028 : {
1029 : void *pResult;
1030 : int bReversed;
1031 1 : char **papszMD = NULL;
1032 : CPLXMLNode *psMDI, *psMetadata;
1033 :
1034 : /* -------------------------------------------------------------------- */
1035 : /* Collect metadata. */
1036 : /* -------------------------------------------------------------------- */
1037 1 : psMetadata = CPLGetXMLNode( psTree, "Metadata" );
1038 :
1039 1 : if( psMetadata->eType != CXT_Element
1040 : || !EQUAL(psMetadata->pszValue,"Metadata") )
1041 0 : return NULL;
1042 :
1043 10 : for( psMDI = psMetadata->psChild; psMDI != NULL;
1044 : psMDI = psMDI->psNext )
1045 : {
1046 9 : if( !EQUAL(psMDI->pszValue,"MDI")
1047 : || psMDI->eType != CXT_Element
1048 : || psMDI->psChild == NULL
1049 : || psMDI->psChild->psNext == NULL
1050 : || psMDI->psChild->eType != CXT_Attribute
1051 : || psMDI->psChild->psChild == NULL )
1052 0 : continue;
1053 :
1054 : papszMD =
1055 : CSLSetNameValue( papszMD,
1056 : psMDI->psChild->psChild->pszValue,
1057 9 : psMDI->psChild->psNext->pszValue );
1058 : }
1059 :
1060 : /* -------------------------------------------------------------------- */
1061 : /* Get other flags. */
1062 : /* -------------------------------------------------------------------- */
1063 1 : bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
1064 :
1065 : /* -------------------------------------------------------------------- */
1066 : /* Generate transformation. */
1067 : /* -------------------------------------------------------------------- */
1068 1 : pResult = GDALCreateGeoLocTransformer( NULL, papszMD, bReversed );
1069 :
1070 : /* -------------------------------------------------------------------- */
1071 : /* Cleanup GCP copy. */
1072 : /* -------------------------------------------------------------------- */
1073 1 : CSLDestroy( papszMD );
1074 :
1075 1 : return pResult;
1076 : }
|
