1 /****************************************************************************** 2 * 3 * Project: GDAL 4 * Purpose: Implements Geolocation array based transformer. 5 * Author: Frank Warmerdam, warmerdam@pobox.com 6 * 7 ****************************************************************************** 8 * Copyright (c) 2006, Frank Warmerdam <warmerdam@pobox.com> 9 * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.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 "cpl_port.h" 31 #include "gdal_alg.h" 32 33 #include <climits> 34 #include <cmath> 35 #include <cstdlib> 36 #include <cstring> 37 38 #include <algorithm> 39 #include <limits> 40 41 #include "cpl_conv.h" 42 #include "cpl_error.h" 43 #include "cpl_minixml.h" 44 #include "cpl_string.h" 45 #include "cpl_vsi.h" 46 #include "gdal.h" 47 #include "gdal_priv.h" 48 49 CPL_CVSID("$Id: gdalgeoloc.cpp 040f61f730ba200425e9791d8cf2511ba978751b 2020-02-27 23:24:20 +0100 Even Rouault $") 50 51 CPL_C_START 52 CPLXMLNode *GDALSerializeGeoLocTransformer( void *pTransformArg ); 53 void *GDALDeserializeGeoLocTransformer( CPLXMLNode *psTree ); 54 CPL_C_END 55 56 /************************************************************************/ 57 /* ==================================================================== */ 58 /* GDALGeoLocTransformer */ 59 /* ==================================================================== */ 60 /************************************************************************/ 61 62 63 //Constants to track down systematic shifts 64 const double FSHIFT = 0.5; 65 const double ISHIFT = 0.5; 66 const double OVERSAMPLE_FACTOR=1.3; 67 68 typedef struct { 69 GDALTransformerInfo sTI; 70 71 bool bReversed; 72 73 // Map from target georef coordinates back to geolocation array 74 // pixel line coordinates. Built only if needed. 75 size_t nBackMapWidth; 76 size_t nBackMapHeight; 77 double adfBackMapGeoTransform[6]; // Maps georef to pixel/line. 78 float *pafBackMapX; 79 float *pafBackMapY; 80 81 // Geolocation bands. 82 GDALDatasetH hDS_X; 83 GDALRasterBandH hBand_X; 84 GDALDatasetH hDS_Y; 85 GDALRasterBandH hBand_Y; 86 int bSwapXY; 87 88 // Located geolocation data. 89 size_t nGeoLocXSize; 90 size_t nGeoLocYSize; 91 double *padfGeoLocX; 92 double *padfGeoLocY; 93 94 int bHasNoData; 95 double dfNoDataX; 96 97 // Geolocation <-> base image mapping. 98 double dfPIXEL_OFFSET; 99 double dfPIXEL_STEP; 100 double dfLINE_OFFSET; 101 double dfLINE_STEP; 102 103 char ** papszGeolocationInfo; 104 105 } GDALGeoLocTransformInfo; 106 107 /************************************************************************/ 108 /* GeoLocLoadFullData() */ 109 /************************************************************************/ 110 111 static bool GeoLocLoadFullData( GDALGeoLocTransformInfo *psTransform ) 112 113 { 114 const int nXSize_XBand = GDALGetRasterXSize( psTransform->hDS_X ); 115 const int nYSize_XBand = GDALGetRasterYSize( psTransform->hDS_X ); 116 const int nXSize_YBand = GDALGetRasterXSize( psTransform->hDS_Y ); 117 const int nYSize_YBand = GDALGetRasterYSize( psTransform->hDS_Y ); 118 119 // Is it a regular grid ? That is: 120 // The XBAND contains the x coordinates for all lines. 121 // The YBAND contains the y coordinates for all columns. 122 const bool bIsRegularGrid = ( nYSize_XBand == 1 && nYSize_YBand == 1 ); 123 124 const int nXSize = nXSize_XBand; 125 int nYSize = 0; 126 if( bIsRegularGrid ) 127 { 128 nYSize = nXSize_YBand; 129 } 130 else 131 { 132 nYSize = nYSize_XBand; 133 } 134 135 psTransform->nGeoLocXSize = static_cast<size_t>(nXSize); 136 psTransform->nGeoLocYSize = static_cast<size_t>(nYSize); 137 138 psTransform->padfGeoLocY = static_cast<double *>( 139 VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize)); 140 psTransform->padfGeoLocX = static_cast<double *>( 141 VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize)); 142 143 if( psTransform->padfGeoLocX == nullptr || 144 psTransform->padfGeoLocY == nullptr ) 145 { 146 return false; 147 } 148 149 if( bIsRegularGrid ) 150 { 151 // Case of regular grid. 152 // The XBAND contains the x coordinates for all lines. 153 // The YBAND contains the y coordinates for all columns. 154 155 double* padfTempX = static_cast<double *>( 156 VSI_MALLOC2_VERBOSE(nXSize, sizeof(double))); 157 double* padfTempY = static_cast<double *>( 158 VSI_MALLOC2_VERBOSE(nYSize, sizeof(double))); 159 if( padfTempX == nullptr || padfTempY == nullptr ) 160 { 161 CPLFree(padfTempX); 162 CPLFree(padfTempY); 163 return false; 164 } 165 166 CPLErr eErr = 167 GDALRasterIO( psTransform->hBand_X, GF_Read, 168 0, 0, nXSize, 1, 169 padfTempX, nXSize, 1, 170 GDT_Float64, 0, 0 ); 171 172 for( size_t j = 0; j < static_cast<size_t>(nYSize); j++ ) 173 { 174 memcpy( psTransform->padfGeoLocX + j * nXSize, 175 padfTempX, 176 nXSize * sizeof(double) ); 177 } 178 179 if( eErr == CE_None ) 180 { 181 eErr = GDALRasterIO( psTransform->hBand_Y, GF_Read, 182 0, 0, nYSize, 1, 183 padfTempY, nYSize, 1, 184 GDT_Float64, 0, 0 ); 185 186 for( size_t j = 0; j < static_cast<size_t>(nYSize); j++ ) 187 { 188 for( size_t i = 0; i < static_cast<size_t>(nXSize); i++ ) 189 { 190 psTransform->padfGeoLocY[j * nXSize + i] = padfTempY[j]; 191 } 192 } 193 } 194 195 CPLFree(padfTempX); 196 CPLFree(padfTempY); 197 198 if( eErr != CE_None ) 199 return false; 200 } 201 else 202 { 203 if( GDALRasterIO( psTransform->hBand_X, GF_Read, 204 0, 0, nXSize, nYSize, 205 psTransform->padfGeoLocX, nXSize, nYSize, 206 GDT_Float64, 0, 0 ) != CE_None 207 || GDALRasterIO( psTransform->hBand_Y, GF_Read, 208 0, 0, nXSize, nYSize, 209 psTransform->padfGeoLocY, nXSize, nYSize, 210 GDT_Float64, 0, 0 ) != CE_None ) 211 return false; 212 } 213 214 psTransform->dfNoDataX = 215 GDALGetRasterNoDataValue( psTransform->hBand_X, 216 &(psTransform->bHasNoData) ); 217 218 return true; 219 } 220 221 /************************************************************************/ 222 /* GeoLocGenerateBackMap() */ 223 /************************************************************************/ 224 225 static bool GeoLocGenerateBackMap( GDALGeoLocTransformInfo *psTransform ) 226 227 { 228 const size_t nXSize = psTransform->nGeoLocXSize; 229 const size_t nYSize = psTransform->nGeoLocYSize; 230 const size_t nXYCount = nXSize * nYSize; 231 const int nMaxIter = 3; 232 233 /* -------------------------------------------------------------------- */ 234 /* Scan forward map for lat/long extents. */ 235 /* -------------------------------------------------------------------- */ 236 double dfMinX = 0.0; 237 double dfMaxX = 0.0; 238 double dfMinY = 0.0; 239 double dfMaxY = 0.0; 240 bool bInit = false; 241 242 for( size_t i = 0; i < nXYCount; i++ ) 243 { 244 if( !psTransform->bHasNoData || 245 psTransform->padfGeoLocX[i] != psTransform->dfNoDataX ) 246 { 247 if( bInit ) 248 { 249 dfMinX = std::min(dfMinX, psTransform->padfGeoLocX[i]); 250 dfMaxX = std::max(dfMaxX, psTransform->padfGeoLocX[i]); 251 dfMinY = std::min(dfMinY, psTransform->padfGeoLocY[i]); 252 dfMaxY = std::max(dfMaxY, psTransform->padfGeoLocY[i]); 253 } 254 else 255 { 256 bInit = true; 257 dfMinX = psTransform->padfGeoLocX[i]; 258 dfMaxX = psTransform->padfGeoLocX[i]; 259 dfMinY = psTransform->padfGeoLocY[i]; 260 dfMaxY = psTransform->padfGeoLocY[i]; 261 } 262 } 263 } 264 265 /* -------------------------------------------------------------------- */ 266 /* Decide on resolution for backmap. We aim for slightly */ 267 /* higher resolution than the source but we can't easily */ 268 /* establish how much dead space there is in the backmap, so it */ 269 /* is approximate. */ 270 /* -------------------------------------------------------------------- */ 271 const double dfTargetPixels = (static_cast<double>(nXSize) * nYSize * OVERSAMPLE_FACTOR); 272 const double dfPixelSize = sqrt((dfMaxX - dfMinX) * (dfMaxY - dfMinY) 273 / dfTargetPixels); 274 if( dfPixelSize == 0.0 ) 275 { 276 CPLError(CE_Failure, CPLE_AppDefined, "Invalid pixel size for backmap"); 277 return false; 278 } 279 280 const double dfBMXSize = (dfMaxX - dfMinX) / dfPixelSize + 1; 281 const double dfBMYSize = (dfMaxY - dfMinY) / dfPixelSize + 1; 282 283 if( !(dfBMXSize > 0 && dfBMXSize < INT_MAX) || 284 !(dfBMYSize > 0 && dfBMYSize < INT_MAX) ) 285 { 286 CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %f x %f", 287 dfBMXSize, dfBMYSize); 288 return false; 289 } 290 291 const size_t nBMXSize = static_cast<size_t>(dfBMXSize); 292 const size_t nBMYSize = static_cast<size_t>(dfBMYSize); 293 294 if( nBMYSize > std::numeric_limits<size_t>::max() / nBMXSize ) 295 { 296 CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %f x %f", 297 dfBMXSize, dfBMYSize); 298 return false; 299 } 300 301 psTransform->nBackMapWidth = nBMXSize; 302 psTransform->nBackMapHeight = nBMYSize; 303 304 dfMinX -= dfPixelSize / 2.0; 305 dfMaxY += dfPixelSize / 2.0; 306 307 308 psTransform->adfBackMapGeoTransform[0] = dfMinX; 309 psTransform->adfBackMapGeoTransform[1] = dfPixelSize; 310 psTransform->adfBackMapGeoTransform[2] = 0.0; 311 psTransform->adfBackMapGeoTransform[3] = dfMaxY; 312 psTransform->adfBackMapGeoTransform[4] = 0.0; 313 psTransform->adfBackMapGeoTransform[5] = -dfPixelSize; 314 315 /* -------------------------------------------------------------------- */ 316 /* Allocate backmap, and initialize to nodata value (-1.0). */ 317 /* -------------------------------------------------------------------- */ 318 GByte *pabyValidFlag = static_cast<GByte *>( 319 VSI_CALLOC_VERBOSE(nBMXSize, nBMYSize)); 320 321 psTransform->pafBackMapX = static_cast<float *>( 322 VSI_MALLOC3_VERBOSE(nBMXSize, nBMYSize, sizeof(float))); 323 psTransform->pafBackMapY = static_cast<float *>( 324 VSI_MALLOC3_VERBOSE(nBMXSize, nBMYSize, sizeof(float))); 325 326 float *wgtsBackMap = static_cast<float *>( 327 VSI_MALLOC3_VERBOSE(nBMXSize, nBMYSize, sizeof(float))); 328 329 if( pabyValidFlag == nullptr || 330 psTransform->pafBackMapX == nullptr || 331 psTransform->pafBackMapY == nullptr || 332 wgtsBackMap == nullptr) 333 { 334 CPLFree( pabyValidFlag ); 335 CPLFree( wgtsBackMap ); 336 return false; 337 } 338 339 const size_t nBMXYCount = nBMXSize * nBMYSize; 340 for( size_t i = 0; i < nBMXYCount; i++ ) 341 { 342 psTransform->pafBackMapX[i] = 0.0; 343 psTransform->pafBackMapY[i] = 0.0; 344 wgtsBackMap[i] = 0.0; 345 pabyValidFlag[i] = 0; 346 } 347 348 /* -------------------------------------------------------------------- */ 349 /* Run through the whole geoloc array forward projecting and */ 350 /* pushing into the backmap. */ 351 /* Initialize to the nMaxIter+1 value so we can spot genuinely */ 352 /* valid pixels in the hole-filling loop. */ 353 /* -------------------------------------------------------------------- */ 354 355 for( size_t iY = 0; iY < nYSize; iY++ ) 356 { 357 for( size_t iX = 0; iX < nXSize; iX++ ) 358 { 359 if( psTransform->bHasNoData && 360 psTransform->padfGeoLocX[iX + iY * nXSize] 361 == psTransform->dfNoDataX ) 362 continue; 363 364 const size_t i = iX + iY * nXSize; 365 366 const double dBMX = static_cast<double>( 367 (psTransform->padfGeoLocX[i] - dfMinX) / dfPixelSize) - FSHIFT; 368 369 const double dBMY = static_cast<double>( 370 (dfMaxY - psTransform->padfGeoLocY[i]) / dfPixelSize) - FSHIFT; 371 372 373 //Get top left index by truncation 374 const int iBMX = static_cast<int>(dBMX); 375 const int iBMY = static_cast<int>(dBMY); 376 const double fracBMX = dBMX - iBMX; 377 const double fracBMY = dBMY - iBMY; 378 379 //Check if the center is in range 380 if( iBMX < -1 || iBMY < -1 || 381 (iBMX > 0 && static_cast<size_t>(iBMX) > nBMXSize) || 382 (iBMY > 0 && static_cast<size_t>(iBMY) > nBMYSize) ) 383 continue; 384 385 //Check logic for top left pixel 386 if ((iBMX >= 0) && (iBMY >= 0) && 387 (static_cast<size_t>(iBMX) < nBMXSize) && 388 (static_cast<size_t>(iBMY) < nBMYSize)) 389 { 390 const double tempwt = (1.0 - fracBMX) * (1.0 - fracBMY); 391 psTransform->pafBackMapX[iBMX + iBMY * nBMXSize] += 392 static_cast<float>( tempwt * ( 393 (iX + FSHIFT) * psTransform->dfPIXEL_STEP + 394 psTransform->dfPIXEL_OFFSET)); 395 396 psTransform->pafBackMapY[iBMX + iBMY * nBMXSize] += 397 static_cast<float>( tempwt * ( 398 (iY + FSHIFT) * psTransform->dfLINE_STEP + 399 psTransform->dfLINE_OFFSET)); 400 wgtsBackMap[iBMX + iBMY * nBMXSize] += static_cast<float>(tempwt); 401 402 //For backward compatibility 403 pabyValidFlag[iBMX + iBMY * nBMXSize] = static_cast<GByte>(nMaxIter+1); 404 } 405 406 //Check logic for top right pixel 407 if ((iBMY >= 0) && 408 (static_cast<size_t>(iBMX+1) < nBMXSize) && 409 (static_cast<size_t>(iBMY) < nBMYSize)) 410 { 411 const double tempwt = fracBMX * (1.0 - fracBMY); 412 413 psTransform->pafBackMapX[iBMX + 1 + iBMY * nBMXSize] += 414 static_cast<float>( tempwt * ( 415 (iX + FSHIFT) * psTransform->dfPIXEL_STEP + 416 psTransform->dfPIXEL_OFFSET)); 417 418 psTransform->pafBackMapY[iBMX + 1 + iBMY * nBMXSize] += 419 static_cast<float>( tempwt * ( 420 (iY + FSHIFT)* psTransform->dfLINE_STEP + 421 psTransform->dfLINE_OFFSET)); 422 wgtsBackMap[iBMX + 1 + iBMY * nBMXSize] += static_cast<float>(tempwt); 423 424 //For backward compatibility 425 pabyValidFlag[iBMX + 1 + iBMY * nBMXSize] = static_cast<GByte>(nMaxIter+1); 426 } 427 428 //Check logic for bottom right pixel 429 if ((static_cast<size_t>(iBMX+1) < nBMXSize) && 430 (static_cast<size_t>(iBMY+1) < nBMYSize)) 431 { 432 const double tempwt = fracBMX * fracBMY; 433 psTransform->pafBackMapX[iBMX + 1 + (iBMY+1) * nBMXSize] += 434 static_cast<float>( tempwt * ( 435 (iX + FSHIFT) * psTransform->dfPIXEL_STEP + 436 psTransform->dfPIXEL_OFFSET)); 437 438 psTransform->pafBackMapY[iBMX + 1 + (iBMY+1) * nBMXSize] += 439 static_cast<float>( tempwt * ( 440 (iY + FSHIFT) * psTransform->dfLINE_STEP + 441 psTransform->dfLINE_OFFSET)); 442 wgtsBackMap[iBMX + 1 + (iBMY+1) * nBMXSize] += static_cast<float>(tempwt); 443 444 //For backward compatibility 445 pabyValidFlag[iBMX + 1 + (iBMY+1) * nBMXSize] = static_cast<GByte>(nMaxIter+1); 446 } 447 448 //Check logic for bottom left pixel 449 if ((iBMX >= 0) && 450 (static_cast<size_t>(iBMX) < nBMXSize) && 451 (static_cast<size_t>(iBMY+1) < nBMYSize)) 452 { 453 const double tempwt = (1.0 - fracBMX) * fracBMY; 454 psTransform->pafBackMapX[iBMX + (iBMY+1) * nBMXSize] += 455 static_cast<float>( tempwt * ( 456 (iX + FSHIFT) * psTransform->dfPIXEL_STEP + 457 psTransform->dfPIXEL_OFFSET)); 458 459 psTransform->pafBackMapY[iBMX + (iBMY+1) * nBMXSize] += 460 static_cast<float>(tempwt * ( 461 (iY + FSHIFT) * psTransform->dfLINE_STEP + 462 psTransform->dfLINE_OFFSET)); 463 wgtsBackMap[iBMX + (iBMY+1) * nBMXSize] += static_cast<float>(tempwt); 464 465 //For backward compatibility 466 pabyValidFlag[iBMX + (iBMY+1) * nBMXSize] = static_cast<GByte>(nMaxIter+1); 467 } 468 469 } 470 } 471 472 473 //Each pixel in the backmap may have multiple entries. 474 //We now go in average it out using the weights 475 for( size_t i = 0; i < nBMXYCount; i++ ) 476 { 477 //Setting these to -1 for backward compatibility 478 if (pabyValidFlag[i] == 0) 479 { 480 psTransform->pafBackMapX[i] = -1.0; 481 psTransform->pafBackMapY[i] = -1.0; 482 } 483 else 484 { 485 //Check if pixel was only touch during neighbor scan 486 //But no real weight was added as source point matched 487 //backmap grid node 488 if (wgtsBackMap[i] > 0) 489 { 490 psTransform->pafBackMapX[i] /= wgtsBackMap[i]; 491 psTransform->pafBackMapY[i] /= wgtsBackMap[i]; 492 pabyValidFlag[i] = static_cast<GByte>(nMaxIter+1); 493 } 494 else 495 { 496 psTransform->pafBackMapX[i] = -1.0; 497 psTransform->pafBackMapY[i] = -1.0; 498 pabyValidFlag[i] = 0; 499 } 500 } 501 } 502 503 /* -------------------------------------------------------------------- */ 504 /* Now, loop over the backmap trying to fill in holes with */ 505 /* nearby values. */ 506 /* -------------------------------------------------------------------- */ 507 for( int iIter = 0; iIter < nMaxIter; iIter++ ) 508 { 509 size_t nNumValid = 0; 510 for( size_t iBMY = 0; iBMY < nBMYSize; iBMY++ ) 511 { 512 for( size_t iBMX = 0; iBMX < nBMXSize; iBMX++ ) 513 { 514 // If this point is already set, ignore it. 515 if( pabyValidFlag[iBMX + iBMY*nBMXSize] ) 516 { 517 nNumValid++; 518 continue; 519 } 520 521 int nCount = 0; 522 double dfXSum = 0.0; 523 double dfYSum = 0.0; 524 const int nMarkedAsGood = nMaxIter - iIter; 525 526 // Left? 527 if( iBMX > 0 && 528 pabyValidFlag[iBMX-1+iBMY*nBMXSize] > nMarkedAsGood ) 529 { 530 dfXSum += psTransform->pafBackMapX[iBMX-1+iBMY*nBMXSize]; 531 dfYSum += psTransform->pafBackMapY[iBMX-1+iBMY*nBMXSize]; 532 nCount++; 533 } 534 // Right? 535 if( iBMX + 1 < nBMXSize && 536 pabyValidFlag[iBMX+1+iBMY*nBMXSize] > nMarkedAsGood ) 537 { 538 dfXSum += psTransform->pafBackMapX[iBMX+1+iBMY*nBMXSize]; 539 dfYSum += psTransform->pafBackMapY[iBMX+1+iBMY*nBMXSize]; 540 nCount++; 541 } 542 // Top? 543 if( iBMY > 0 && 544 pabyValidFlag[iBMX+(iBMY-1)*nBMXSize] > nMarkedAsGood ) 545 { 546 dfXSum += psTransform->pafBackMapX[iBMX+(iBMY-1)*nBMXSize]; 547 dfYSum += psTransform->pafBackMapY[iBMX+(iBMY-1)*nBMXSize]; 548 nCount++; 549 } 550 // Bottom? 551 if( iBMY + 1 < nBMYSize && 552 pabyValidFlag[iBMX+(iBMY+1)*nBMXSize] > nMarkedAsGood ) 553 { 554 dfXSum += psTransform->pafBackMapX[iBMX+(iBMY+1)*nBMXSize]; 555 dfYSum += psTransform->pafBackMapY[iBMX+(iBMY+1)*nBMXSize]; 556 nCount++; 557 } 558 // Top-left? 559 if( iBMX > 0 && iBMY > 0 && 560 pabyValidFlag[iBMX-1+(iBMY-1)*nBMXSize] > nMarkedAsGood ) 561 { 562 dfXSum += 563 psTransform->pafBackMapX[iBMX-1+(iBMY-1)*nBMXSize]; 564 dfYSum += 565 psTransform->pafBackMapY[iBMX-1+(iBMY-1)*nBMXSize]; 566 nCount++; 567 } 568 // Top-right? 569 if( iBMX + 1 < nBMXSize && iBMY > 0 && 570 pabyValidFlag[iBMX+1+(iBMY-1)*nBMXSize] > nMarkedAsGood ) 571 { 572 dfXSum += 573 psTransform->pafBackMapX[iBMX+1+(iBMY-1)*nBMXSize]; 574 dfYSum += 575 psTransform->pafBackMapY[iBMX+1+(iBMY-1)*nBMXSize]; 576 nCount++; 577 } 578 // Bottom-left? 579 if( iBMX > 0 && iBMY + 1 < nBMYSize && 580 pabyValidFlag[iBMX-1+(iBMY+1)*nBMXSize] > nMarkedAsGood ) 581 { 582 dfXSum += 583 psTransform->pafBackMapX[iBMX-1+(iBMY+1)*nBMXSize]; 584 dfYSum += 585 psTransform->pafBackMapY[iBMX-1+(iBMY+1)*nBMXSize]; 586 nCount++; 587 } 588 // Bottom-right? 589 if( iBMX + 1 < nBMXSize && iBMY + 1 < nBMYSize && 590 pabyValidFlag[iBMX+1+(iBMY+1)*nBMXSize] > nMarkedAsGood ) 591 { 592 dfXSum += 593 psTransform->pafBackMapX[iBMX+1+(iBMY+1)*nBMXSize]; 594 dfYSum += 595 psTransform->pafBackMapY[iBMX+1+(iBMY+1)*nBMXSize]; 596 nCount++; 597 } 598 599 if( nCount > 0 ) 600 { 601 psTransform->pafBackMapX[iBMX + iBMY * nBMXSize] = 602 static_cast<float>(dfXSum/nCount); 603 psTransform->pafBackMapY[iBMX + iBMY * nBMXSize] = 604 static_cast<float>(dfYSum/nCount); 605 // Genuinely valid points will have value iMaxIter + 1. 606 // On each iteration mark newly valid points with a 607 // descending value so that it will not be used on the 608 // current iteration only on subsequent ones. 609 pabyValidFlag[iBMX+iBMY*nBMXSize] = 610 static_cast<GByte>(nMaxIter - iIter); 611 } 612 } 613 } 614 if( nNumValid == nBMXSize * nBMYSize ) 615 break; 616 } 617 618 CPLFree( wgtsBackMap ); 619 CPLFree( pabyValidFlag ); 620 621 return true; 622 } 623 624 /************************************************************************/ 625 /* GDALGeoLocRescale() */ 626 /************************************************************************/ 627 628 static void GDALGeoLocRescale( char**& papszMD, const char* pszItem, 629 double dfRatio, double dfDefaultVal ) 630 { 631 const double dfVal = 632 dfRatio * 633 CPLAtofM(CSLFetchNameValueDef(papszMD, pszItem, 634 CPLSPrintf("%.18g", dfDefaultVal))); 635 636 papszMD = CSLSetNameValue(papszMD, pszItem, CPLSPrintf("%.18g", dfVal)); 637 638 } 639 640 /************************************************************************/ 641 /* GDALCreateSimilarGeoLocTransformer() */ 642 /************************************************************************/ 643 644 static 645 void* GDALCreateSimilarGeoLocTransformer( void *hTransformArg, 646 double dfRatioX, double dfRatioY ) 647 { 648 VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarGeoLocTransformer", 649 nullptr); 650 651 GDALGeoLocTransformInfo *psInfo = 652 static_cast<GDALGeoLocTransformInfo *>(hTransformArg); 653 654 char** papszGeolocationInfo = CSLDuplicate(psInfo->papszGeolocationInfo); 655 656 if( dfRatioX != 1.0 || dfRatioY != 1.0 ) 657 { 658 GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_OFFSET", dfRatioX, 0.0); 659 GDALGeoLocRescale(papszGeolocationInfo, "LINE_OFFSET", dfRatioY, 0.0); 660 GDALGeoLocRescale( 661 papszGeolocationInfo, "PIXEL_STEP", 1.0 / dfRatioX, 1.0); 662 GDALGeoLocRescale(papszGeolocationInfo, 663 "LINE_STEP", 1.0 / dfRatioY, 1.0); 664 } 665 666 psInfo = static_cast<GDALGeoLocTransformInfo*>( 667 GDALCreateGeoLocTransformer( 668 nullptr, papszGeolocationInfo, psInfo->bReversed)); 669 670 CSLDestroy(papszGeolocationInfo); 671 672 return psInfo; 673 } 674 675 /************************************************************************/ 676 /* GDALCreateGeoLocTransformer() */ 677 /************************************************************************/ 678 679 /** Create GeoLocation transformer */ 680 void *GDALCreateGeoLocTransformer( GDALDatasetH hBaseDS, 681 char **papszGeolocationInfo, 682 int bReversed ) 683 684 { 685 686 if( CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET") == nullptr 687 || CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET") == nullptr 688 || CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP") == nullptr 689 || CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP") == nullptr 690 || CSLFetchNameValue(papszGeolocationInfo, "X_BAND") == nullptr 691 || CSLFetchNameValue(papszGeolocationInfo, "Y_BAND") == nullptr ) 692 { 693 CPLError( CE_Failure, CPLE_AppDefined, 694 "Missing some geolocation fields in " 695 "GDALCreateGeoLocTransformer()" ); 696 return nullptr; 697 } 698 699 /* -------------------------------------------------------------------- */ 700 /* Initialize core info. */ 701 /* -------------------------------------------------------------------- */ 702 GDALGeoLocTransformInfo *psTransform = 703 static_cast<GDALGeoLocTransformInfo *>( 704 CPLCalloc(sizeof(GDALGeoLocTransformInfo), 1)); 705 706 psTransform->bReversed = CPL_TO_BOOL(bReversed); 707 708 memcpy( psTransform->sTI.abySignature, 709 GDAL_GTI2_SIGNATURE, 710 strlen(GDAL_GTI2_SIGNATURE) ); 711 psTransform->sTI.pszClassName = "GDALGeoLocTransformer"; 712 psTransform->sTI.pfnTransform = GDALGeoLocTransform; 713 psTransform->sTI.pfnCleanup = GDALDestroyGeoLocTransformer; 714 psTransform->sTI.pfnSerialize = GDALSerializeGeoLocTransformer; 715 psTransform->sTI.pfnCreateSimilar = GDALCreateSimilarGeoLocTransformer; 716 717 psTransform->papszGeolocationInfo = CSLDuplicate( papszGeolocationInfo ); 718 719 /* -------------------------------------------------------------------- */ 720 /* Pull geolocation info from the options/metadata. */ 721 /* -------------------------------------------------------------------- */ 722 psTransform->dfPIXEL_OFFSET = 723 CPLAtof(CSLFetchNameValue( papszGeolocationInfo, "PIXEL_OFFSET" )); 724 psTransform->dfLINE_OFFSET = 725 CPLAtof(CSLFetchNameValue( papszGeolocationInfo, "LINE_OFFSET" )); 726 psTransform->dfPIXEL_STEP = 727 CPLAtof(CSLFetchNameValue( papszGeolocationInfo, "PIXEL_STEP" )); 728 psTransform->dfLINE_STEP = 729 CPLAtof(CSLFetchNameValue( papszGeolocationInfo, "LINE_STEP" )); 730 731 /* -------------------------------------------------------------------- */ 732 /* Establish access to geolocation dataset(s). */ 733 /* -------------------------------------------------------------------- */ 734 const char *pszDSName = CSLFetchNameValue( papszGeolocationInfo, 735 "X_DATASET" ); 736 if( pszDSName != nullptr ) 737 { 738 CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true); 739 psTransform->hDS_X = GDALOpenShared( pszDSName, GA_ReadOnly ); 740 } 741 else 742 { 743 psTransform->hDS_X = hBaseDS; 744 if( hBaseDS ) 745 { 746 GDALReferenceDataset( psTransform->hDS_X ); 747 psTransform->papszGeolocationInfo = 748 CSLSetNameValue( psTransform->papszGeolocationInfo, 749 "X_DATASET", 750 GDALGetDescription( hBaseDS ) ); 751 } 752 } 753 754 pszDSName = CSLFetchNameValue( papszGeolocationInfo, "Y_DATASET" ); 755 if( pszDSName != nullptr ) 756 { 757 CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true); 758 psTransform->hDS_Y = GDALOpenShared( pszDSName, GA_ReadOnly ); 759 } 760 else 761 { 762 psTransform->hDS_Y = hBaseDS; 763 if( hBaseDS ) 764 { 765 GDALReferenceDataset( psTransform->hDS_Y ); 766 psTransform->papszGeolocationInfo = 767 CSLSetNameValue( psTransform->papszGeolocationInfo, 768 "Y_DATASET", 769 GDALGetDescription( hBaseDS ) ); 770 } 771 } 772 773 if( psTransform->hDS_X == nullptr || 774 psTransform->hDS_Y == nullptr ) 775 { 776 GDALDestroyGeoLocTransformer( psTransform ); 777 return nullptr; 778 } 779 780 /* -------------------------------------------------------------------- */ 781 /* Get the band handles. */ 782 /* -------------------------------------------------------------------- */ 783 const int nXBand = 784 std::max(1, atoi(CSLFetchNameValue( papszGeolocationInfo, "X_BAND" ))); 785 psTransform->hBand_X = GDALGetRasterBand( psTransform->hDS_X, nXBand ); 786 787 const int nYBand = 788 std::max(1, atoi(CSLFetchNameValue( papszGeolocationInfo, "Y_BAND" ))); 789 psTransform->hBand_Y = GDALGetRasterBand( psTransform->hDS_Y, nYBand ); 790 791 if( psTransform->hBand_X == nullptr || 792 psTransform->hBand_Y == nullptr ) 793 { 794 GDALDestroyGeoLocTransformer( psTransform ); 795 return nullptr; 796 } 797 798 psTransform->bSwapXY = CPLTestBool(CSLFetchNameValueDef( 799 papszGeolocationInfo, "SWAP_XY", "NO")); 800 801 /* -------------------------------------------------------------------- */ 802 /* Check that X and Y bands have the same dimensions */ 803 /* -------------------------------------------------------------------- */ 804 const int nXSize_XBand = GDALGetRasterXSize( psTransform->hDS_X ); 805 const int nYSize_XBand = GDALGetRasterYSize( psTransform->hDS_X ); 806 const int nXSize_YBand = GDALGetRasterXSize( psTransform->hDS_Y ); 807 const int nYSize_YBand = GDALGetRasterYSize( psTransform->hDS_Y ); 808 if( nYSize_XBand == 1 || nYSize_YBand == 1 ) 809 { 810 if( nYSize_XBand != 1 || nYSize_YBand != 1 ) 811 { 812 CPLError(CE_Failure, CPLE_AppDefined, 813 "X_BAND and Y_BAND should have both nYSize == 1"); 814 GDALDestroyGeoLocTransformer( psTransform ); 815 return nullptr; 816 } 817 } 818 else if( nXSize_XBand != nXSize_YBand || 819 nYSize_XBand != nYSize_YBand ) 820 { 821 CPLError(CE_Failure, CPLE_AppDefined, 822 "X_BAND and Y_BAND do not have the same dimensions"); 823 GDALDestroyGeoLocTransformer( psTransform ); 824 return nullptr; 825 } 826 827 if( static_cast<size_t>(nXSize_XBand) > 828 std::numeric_limits<size_t>::max() / nYSize_XBand ) 829 { 830 CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d", 831 nXSize_XBand, nYSize_XBand); 832 GDALDestroyGeoLocTransformer( psTransform ); 833 return nullptr; 834 } 835 836 /* -------------------------------------------------------------------- */ 837 /* Load the geolocation array. */ 838 /* -------------------------------------------------------------------- */ 839 if( !GeoLocLoadFullData( psTransform ) 840 || !GeoLocGenerateBackMap( psTransform ) ) 841 { 842 GDALDestroyGeoLocTransformer( psTransform ); 843 return nullptr; 844 } 845 846 return psTransform; 847 } 848 849 /************************************************************************/ 850 /* GDALDestroyGeoLocTransformer() */ 851 /************************************************************************/ 852 853 /** Destroy GeoLocation transformer */ 854 void GDALDestroyGeoLocTransformer( void *pTransformAlg ) 855 856 { 857 if( pTransformAlg == nullptr ) 858 return; 859 860 GDALGeoLocTransformInfo *psTransform = 861 static_cast<GDALGeoLocTransformInfo *>(pTransformAlg); 862 863 CPLFree( psTransform->pafBackMapX ); 864 CPLFree( psTransform->pafBackMapY ); 865 CSLDestroy( psTransform->papszGeolocationInfo ); 866 CPLFree( psTransform->padfGeoLocX ); 867 CPLFree( psTransform->padfGeoLocY ); 868 869 if( psTransform->hDS_X != nullptr 870 && GDALDereferenceDataset( psTransform->hDS_X ) == 0 ) 871 GDALClose( psTransform->hDS_X ); 872 873 if( psTransform->hDS_Y != nullptr 874 && GDALDereferenceDataset( psTransform->hDS_Y ) == 0 ) 875 GDALClose( psTransform->hDS_Y ); 876 877 CPLFree( pTransformAlg ); 878 } 879 880 /************************************************************************/ 881 /* GDALGeoLocTransform() */ 882 /************************************************************************/ 883 884 /** Use GeoLocation transformer */ 885 int GDALGeoLocTransform( void *pTransformArg, 886 int bDstToSrc, 887 int nPointCount, 888 double *padfX, double *padfY, 889 CPL_UNUSED double *padfZ, 890 int *panSuccess ) 891 { 892 GDALGeoLocTransformInfo *psTransform = 893 static_cast<GDALGeoLocTransformInfo *>(pTransformArg); 894 895 if( psTransform->bReversed ) 896 bDstToSrc = !bDstToSrc; 897 898 /* -------------------------------------------------------------------- */ 899 /* Do original pixel line to target geox/geoy. */ 900 /* -------------------------------------------------------------------- */ 901 if( !bDstToSrc ) 902 { 903 const size_t nXSize = psTransform->nGeoLocXSize; 904 905 for( int i = 0; i < nPointCount; i++ ) 906 { 907 if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL ) 908 { 909 panSuccess[i] = FALSE; 910 continue; 911 } 912 913 const double dfGeoLocPixel = 914 (padfX[i] - psTransform->dfPIXEL_OFFSET) 915 / psTransform->dfPIXEL_STEP; 916 const double dfGeoLocLine = 917 (padfY[i] - psTransform->dfLINE_OFFSET) 918 / psTransform->dfLINE_STEP; 919 920 size_t iX = static_cast<size_t>(std::max(0.0, dfGeoLocPixel)); 921 iX = std::min(iX, psTransform->nGeoLocXSize-1); 922 size_t iY = static_cast<size_t>(std::max(0.0, dfGeoLocLine)); 923 iY = std::min(iY, psTransform->nGeoLocYSize-1); 924 925 const double *padfGLX = psTransform->padfGeoLocX + iX + iY * nXSize; 926 const double *padfGLY = psTransform->padfGeoLocY + iX + iY * nXSize; 927 928 if( psTransform->bHasNoData && 929 padfGLX[0] == psTransform->dfNoDataX ) 930 { 931 panSuccess[i] = FALSE; 932 padfX[i] = HUGE_VAL; 933 padfY[i] = HUGE_VAL; 934 continue; 935 } 936 937 // This assumes infinite extension beyond borders of available 938 // data based on closest grid square. 939 940 if( iX + 1 < psTransform->nGeoLocXSize && 941 iY + 1 < psTransform->nGeoLocYSize && 942 (!psTransform->bHasNoData || 943 (padfGLX[1] != psTransform->dfNoDataX && 944 padfGLX[nXSize] != psTransform->dfNoDataX && 945 padfGLX[nXSize + 1] != psTransform->dfNoDataX) )) 946 { 947 padfX[i] = 948 (1 - (dfGeoLocLine -iY)) 949 * (padfGLX[0] + 950 (dfGeoLocPixel-iX) * (padfGLX[1] - padfGLX[0])) 951 + (dfGeoLocLine -iY) 952 * (padfGLX[nXSize] + (dfGeoLocPixel-iX) * 953 (padfGLX[nXSize+1] - padfGLX[nXSize])); 954 padfY[i] = 955 (1 - (dfGeoLocLine -iY)) 956 * (padfGLY[0] + 957 (dfGeoLocPixel-iX) * (padfGLY[1] - padfGLY[0])) 958 + (dfGeoLocLine -iY) 959 * (padfGLY[nXSize] + (dfGeoLocPixel-iX) * 960 (padfGLY[nXSize+1] - padfGLY[nXSize])); 961 } 962 else if( iX + 1 < psTransform->nGeoLocXSize && 963 (!psTransform->bHasNoData || 964 padfGLX[1] != psTransform->dfNoDataX) ) 965 { 966 padfX[i] = 967 padfGLX[0] + (dfGeoLocPixel-iX) * (padfGLX[1] - padfGLX[0]); 968 padfY[i] = 969 padfGLY[0] + (dfGeoLocPixel-iX) * (padfGLY[1] - padfGLY[0]); 970 } 971 else if( iY + 1 < psTransform->nGeoLocYSize && 972 (!psTransform->bHasNoData || 973 padfGLX[nXSize] != psTransform->dfNoDataX) ) 974 { 975 padfX[i] = padfGLX[0] 976 + (dfGeoLocLine -iY) * (padfGLX[nXSize] - padfGLX[0]); 977 padfY[i] = padfGLY[0] 978 + (dfGeoLocLine -iY) * (padfGLY[nXSize] - padfGLY[0]); 979 } 980 else 981 { 982 padfX[i] = padfGLX[0]; 983 padfY[i] = padfGLY[0]; 984 } 985 986 if( psTransform->bSwapXY ) 987 { 988 std::swap(padfX[i], padfY[i]); 989 } 990 991 panSuccess[i] = TRUE; 992 } 993 } 994 995 /* -------------------------------------------------------------------- */ 996 /* geox/geoy to pixel/line using backmap. */ 997 /* -------------------------------------------------------------------- */ 998 else 999 { 1000 for( int i = 0; i < nPointCount; i++ ) 1001 { 1002 if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL ) 1003 { 1004 panSuccess[i] = FALSE; 1005 continue; 1006 } 1007 1008 if( psTransform->bSwapXY ) 1009 { 1010 std::swap(padfX[i], padfY[i]); 1011 } 1012 1013 const double dfBMX = 1014 ((padfX[i] - psTransform->adfBackMapGeoTransform[0]) 1015 / psTransform->adfBackMapGeoTransform[1]) - ISHIFT; 1016 const double dfBMY = 1017 ((padfY[i] - psTransform->adfBackMapGeoTransform[3]) 1018 / psTransform->adfBackMapGeoTransform[5]) - ISHIFT; 1019 1020 // FIXME: in the case of ]-1,0[, dfBMX-iBMX will be wrong 1021 // We should likely error out if values are < 0 ==> affects a few 1022 // autotest results 1023 if( !(dfBMX > -1 && dfBMY > -1 && 1024 dfBMX < psTransform->nBackMapWidth && 1025 dfBMY < psTransform->nBackMapHeight) ) 1026 { 1027 panSuccess[i] = FALSE; 1028 padfX[i] = HUGE_VAL; 1029 padfY[i] = HUGE_VAL; 1030 continue; 1031 } 1032 1033 const int iBMX = static_cast<int>(dfBMX); 1034 const int iBMY = static_cast<int>(dfBMY); 1035 1036 const size_t iBM = iBMX + iBMY * psTransform->nBackMapWidth; 1037 if( psTransform->pafBackMapX[iBM] < 0 ) 1038 { 1039 panSuccess[i] = FALSE; 1040 padfX[i] = HUGE_VAL; 1041 padfY[i] = HUGE_VAL; 1042 continue; 1043 } 1044 1045 const float* pafBMX = psTransform->pafBackMapX + iBM; 1046 const float* pafBMY = psTransform->pafBackMapY + iBM; 1047 1048 if( static_cast<size_t>(iBMX + 1) < psTransform->nBackMapWidth && 1049 static_cast<size_t>(iBMY + 1) < psTransform->nBackMapHeight && 1050 pafBMX[1] >=0 && pafBMX[psTransform->nBackMapWidth] >= 0 && 1051 pafBMX[psTransform->nBackMapWidth+1] >= 0) 1052 { 1053 padfX[i] = 1054 (1-(dfBMY - iBMY)) 1055 * (pafBMX[0] + (dfBMX - iBMX) * (pafBMX[1] - pafBMX[0])) 1056 + (dfBMY - iBMY) 1057 * (pafBMX[psTransform->nBackMapWidth] + 1058 (dfBMX - iBMX) * (pafBMX[psTransform->nBackMapWidth+1] - 1059 pafBMX[psTransform->nBackMapWidth])); 1060 padfY[i] = 1061 (1-(dfBMY - iBMY)) 1062 * (pafBMY[0] + (dfBMX - iBMX) * (pafBMY[1] - pafBMY[0])) 1063 + (dfBMY - iBMY) 1064 * (pafBMY[psTransform->nBackMapWidth] + 1065 (dfBMX - iBMX) * (pafBMY[psTransform->nBackMapWidth+1] - 1066 pafBMY[psTransform->nBackMapWidth])); 1067 } 1068 else if( static_cast<size_t>(iBMX + 1) < psTransform->nBackMapWidth && 1069 pafBMX[1] >=0) 1070 { 1071 padfX[i] = pafBMX[0] + 1072 (dfBMX - iBMX) * (pafBMX[1] - pafBMX[0]); 1073 padfY[i] = pafBMY[0] + 1074 (dfBMX - iBMX) * (pafBMY[1] - pafBMY[0]); 1075 } 1076 else if( static_cast<size_t>(iBMY + 1) < psTransform->nBackMapHeight && 1077 pafBMX[psTransform->nBackMapWidth] >= 0 ) 1078 { 1079 padfX[i] = 1080 pafBMX[0] + 1081 (dfBMY - iBMY) * (pafBMX[psTransform->nBackMapWidth] - 1082 pafBMX[0]); 1083 padfY[i] = 1084 pafBMY[0] + 1085 (dfBMY - iBMY) * (pafBMY[psTransform->nBackMapWidth] - 1086 pafBMY[0]); 1087 } 1088 else 1089 { 1090 padfX[i] = pafBMX[0]; 1091 padfY[i] = pafBMY[0]; 1092 } 1093 panSuccess[i] = TRUE; 1094 } 1095 } 1096 1097 return TRUE; 1098 } 1099 1100 /************************************************************************/ 1101 /* GDALSerializeGeoLocTransformer() */ 1102 /************************************************************************/ 1103 1104 CPLXMLNode *GDALSerializeGeoLocTransformer( void *pTransformArg ) 1105 1106 { 1107 VALIDATE_POINTER1( pTransformArg, "GDALSerializeGeoLocTransformer", nullptr ); 1108 1109 GDALGeoLocTransformInfo *psInfo = 1110 static_cast<GDALGeoLocTransformInfo *>(pTransformArg); 1111 1112 CPLXMLNode *psTree = 1113 CPLCreateXMLNode( nullptr, CXT_Element, "GeoLocTransformer" ); 1114 1115 /* -------------------------------------------------------------------- */ 1116 /* Serialize bReversed. */ 1117 /* -------------------------------------------------------------------- */ 1118 CPLCreateXMLElementAndValue( 1119 psTree, "Reversed", 1120 CPLString().Printf( "%d", static_cast<int>(psInfo->bReversed) ) ); 1121 1122 /* -------------------------------------------------------------------- */ 1123 /* geoloc metadata. */ 1124 /* -------------------------------------------------------------------- */ 1125 char **papszMD = psInfo->papszGeolocationInfo; 1126 CPLXMLNode *psMD= CPLCreateXMLNode( psTree, CXT_Element, "Metadata" ); 1127 1128 for( int i = 0; papszMD != nullptr && papszMD[i] != nullptr; i++ ) 1129 { 1130 char *pszKey = nullptr; 1131 const char *pszRawValue = CPLParseNameValue( papszMD[i], &pszKey ); 1132 1133 CPLXMLNode *psMDI = CPLCreateXMLNode( psMD, CXT_Element, "MDI" ); 1134 CPLSetXMLValue( psMDI, "#key", pszKey ); 1135 CPLCreateXMLNode( psMDI, CXT_Text, pszRawValue ); 1136 1137 CPLFree( pszKey ); 1138 } 1139 1140 return psTree; 1141 } 1142 1143 /************************************************************************/ 1144 /* GDALDeserializeGeoLocTransformer() */ 1145 /************************************************************************/ 1146 1147 void *GDALDeserializeGeoLocTransformer( CPLXMLNode *psTree ) 1148 1149 { 1150 /* -------------------------------------------------------------------- */ 1151 /* Collect metadata. */ 1152 /* -------------------------------------------------------------------- */ 1153 CPLXMLNode *psMetadata = CPLGetXMLNode( psTree, "Metadata" ); 1154 1155 if( psMetadata == nullptr || 1156 psMetadata->eType != CXT_Element 1157 || !EQUAL(psMetadata->pszValue, "Metadata") ) 1158 return nullptr; 1159 1160 char **papszMD = nullptr; 1161 1162 for( CPLXMLNode *psMDI = psMetadata->psChild; 1163 psMDI != nullptr; 1164 psMDI = psMDI->psNext ) 1165 { 1166 if( !EQUAL(psMDI->pszValue, "MDI") 1167 || psMDI->eType != CXT_Element 1168 || psMDI->psChild == nullptr 1169 || psMDI->psChild->psNext == nullptr 1170 || psMDI->psChild->eType != CXT_Attribute 1171 || psMDI->psChild->psChild == nullptr ) 1172 continue; 1173 1174 papszMD = 1175 CSLSetNameValue( papszMD, 1176 psMDI->psChild->psChild->pszValue, 1177 psMDI->psChild->psNext->pszValue ); 1178 } 1179 1180 /* -------------------------------------------------------------------- */ 1181 /* Get other flags. */ 1182 /* -------------------------------------------------------------------- */ 1183 const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0")); 1184 1185 /* -------------------------------------------------------------------- */ 1186 /* Generate transformation. */ 1187 /* -------------------------------------------------------------------- */ 1188 void *pResult = GDALCreateGeoLocTransformer( nullptr, papszMD, bReversed ); 1189 1190 /* -------------------------------------------------------------------- */ 1191 /* Cleanup GCP copy. */ 1192 /* -------------------------------------------------------------------- */ 1193 CSLDestroy( papszMD ); 1194 1195 return pResult; 1196 } 1197