otbImageKeywordlist.cxx 15.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
/*=========================================================================

  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.


13 14
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15 16 17 18 19
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#include "otbImageKeywordlist.h"

20 21
#include <cassert>

22 23
#include "otbMacro.h"

24
#include "gdal_priv.h"
25

Julien Malik's avatar
Julien Malik committed
26 27
#include "ossim/base/ossimKeywordlist.h"
#include "ossim/base/ossimString.h"
28
#include "ossim/ossimPluginProjectionFactory.h"
Julien Malik's avatar
Julien Malik committed
29
#include "ossim/imaging/ossimImageHandlerRegistry.h"
30
#include "ossim/ossimTileMapModel.h"
Julien Malik's avatar
Julien Malik committed
31
#include "ossim/projection/ossimProjectionFactoryRegistry.h"
32 33
#include "ossim/projection/ossimRpcModel.h"

Julien Malik's avatar
Julien Malik committed
34
#include "otbSensorModelAdapter.h"
35

36 37 38 39
namespace otb
{

ImageKeywordlist
40
::ImageKeywordlist()
41 42 43 44
{
}

ImageKeywordlist
45
::ImageKeywordlist(const Self& p) : m_Keywordlist(p.m_Keywordlist)
46 47 48 49 50 51 52 53
{
}

ImageKeywordlist
::~ImageKeywordlist()
{
}

54
void
55
ImageKeywordlist::
OTB Bot's avatar
STYLE  
OTB Bot committed
56
operator =(const Self& p)
57
{
58
  m_Keywordlist = p.m_Keywordlist;
59 60
}

61
bool
62 63
ImageKeywordlist
::operator ==(const Self& p) const
64 65 66 67
{
  return m_Keywordlist == p.m_Keywordlist;
}

68 69 70 71
void
ImageKeywordlist::
SetKeywordlist(const ossimKeywordlist& kwl)
{
72 73 74 75 76
  m_Keywordlist.clear();
  for (ossimKeywordlist::KeywordMap::const_iterator it = kwl.getMap().begin();
       it != kwl.getMap().end();
       ++it)
    {
77 78
    std::string first(it->first);
    std::string second(it->second);
79 80
    m_Keywordlist[first] = second;
    }
81 82
}

83
const std::string&
84
ImageKeywordlist::
85
GetMetadataByKey(const std::string& key) const
86 87 88
{
  // Search for the key in the output map
  KeywordlistMap::const_iterator it = m_Keywordlist.find(key);
89

90
  // If the key can not be found, throw an exception
OTB Bot's avatar
STYLE  
OTB Bot committed
91
  if (it == m_Keywordlist.end())
92
    {
93
    itkGenericExceptionMacro(<< "Keywordlist has no output with key " << key);
94 95 96 97 98
    }

  // Then if everything is ok, return the ossinString
  return it->second;
}
99

100 101 102 103 104 105 106 107
bool
ImageKeywordlist::
HasKey(const std::string& key) const
{
  KeywordlistMap::const_iterator it = m_Keywordlist.find(key);
  return (it != m_Keywordlist.end());
}

OTB Bot's avatar
STYLE  
OTB Bot committed
108
void
109 110 111 112 113 114
ImageKeywordlist::
ClearMetadataByKey(const std::string& key)
{
  m_Keywordlist[key] = "";
}

115 116 117 118 119 120 121
void
ImageKeywordlist::
AddKey(const std::string& key, const std::string& value)
{
  m_Keywordlist[key] = value;
}

122
void
123
ImageKeywordlist::
124
convertToOSSIMKeywordlist(ossimKeywordlist& kwl) const
125
{
126 127 128 129 130
  ossimKeywordlist::KeywordMap ossimMap;
  for(KeywordlistMap::const_iterator it = m_Keywordlist.begin();
      it != m_Keywordlist.end();
      ++it)
    {
131
    ossimMap[it->first] = it->second;
132 133
    }
  kwl.getMap() = ossimMap;
134 135
}

136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
bool
ImageKeywordlist::
convertToGDALRPC(GDALRPCInfo &rpc) const
{
  ossimKeywordlist geom_kwl;
  this->convertToOSSIMKeywordlist(geom_kwl);
  
  ossimRefPtr<ossimRpcModel> rpcModel = new ossimRpcModel;
  if (rpcModel->loadState(geom_kwl))
    {
    ossimRpcModel::rpcModelStruct ossimRpcStruct;
    rpcModel->getRpcParameters(ossimRpcStruct);
    
    if (ossimRpcStruct.type == 'B')
      {
      rpc.dfSAMP_OFF = ossimRpcStruct.sampOffset;
      rpc.dfLINE_OFF = ossimRpcStruct.lineOffset;
      rpc.dfSAMP_SCALE = ossimRpcStruct.sampScale;
      rpc.dfLINE_SCALE = ossimRpcStruct.lineScale;
      rpc.dfLAT_OFF = ossimRpcStruct.latOffset;
      rpc.dfLONG_OFF = ossimRpcStruct.lonOffset;
      rpc.dfHEIGHT_OFF = ossimRpcStruct.hgtOffset;
      rpc.dfLAT_SCALE = ossimRpcStruct.latScale;
      rpc.dfLONG_SCALE = ossimRpcStruct.lonScale;
      rpc.dfHEIGHT_SCALE = ossimRpcStruct.hgtScale;
      
      memcpy(rpc.adfLINE_NUM_COEFF, ossimRpcStruct.lineNumCoef, sizeof(double) * 20);
      memcpy(rpc.adfLINE_DEN_COEFF, ossimRpcStruct.lineDenCoef, sizeof(double) * 20);
      memcpy(rpc.adfSAMP_NUM_COEFF, ossimRpcStruct.sampNumCoef, sizeof(double) * 20);
      memcpy(rpc.adfSAMP_DEN_COEFF, ossimRpcStruct.sampDenCoef, sizeof(double) * 20);
      
      return true;
      }
    }
  
OTB Bot's avatar
STYLE  
OTB Bot committed
171
  return false;
172 173
}

174
void
175
ImageKeywordlist::
176
Print(std::ostream& os, itk::Indent indent) const
177
{
178 179 180 181 182 183 184 185 186
  this->PrintSelf(os, indent.GetNextIndent());
}

void
ImageKeywordlist::
PrintSelf(std::ostream& os, itk::Indent indent) const
{
  ossimKeywordlist kwl;
  convertToOSSIMKeywordlist(kwl);
OTB Bot's avatar
STYLE  
OTB Bot committed
187
  os << indent << " Ossim Keyword list:" << std::endl;
Emmanuel Christophe's avatar
Emmanuel Christophe committed
188
  os << indent << kwl;
189 190
}

191
std::ostream &
OTB Bot's avatar
STYLE  
OTB Bot committed
192
operator <<(std::ostream& os, const ImageKeywordlist& kwl)
193 194 195
{
  kwl.Print(os);
  return os;
196 197
}

198 199
ImageKeywordlist
ReadGeometry(const std::string& filename)
200 201 202 203 204
{
  return ReadGeometryFromImage(filename);
}

ImageKeywordlist
205
ReadGeometryFromImage(const std::string& filename, bool checkRpcTag)
206 207 208 209
{
  // Trying to read ossim MetaData
  bool             hasMetaData = false;
  ossimKeywordlist geom_kwl; // = new ossimKeywordlist();
210
  ImageKeywordlist otb_kwl;
211

212 213 214
  /****************************************************/
  /* First try : test the OSSIM plugins factory       */
  /****************************************************/
215 216 217 218 219 220
  /** Before, the pluginfactory was tested if the ossim one returned false.
      But in the case TSX, the images tif were considered as ossimQuickbirdTiffTileSource
      thus a TSX tif image wasn't read with TSX Model. We don't use the ossimRegisteryFactory
      because the default include factory contains ossimQuickbirdTiffTileSource. */
  ossimProjection * projection = ossimplugins::ossimPluginProjectionFactory::instance()
                                 ->createProjection(ossimFilename(filename.c_str()), 0);
221
  if (projection)
222
    {
223
    otbMsgDevMacro(<< "OSSIM plugin projection instantiated ! ");
224

225
    hasMetaData = projection->saveState(geom_kwl);
226
    otb_kwl.SetKeywordlist(geom_kwl);
227 228 229

    // Free memory
    delete projection;
230
    projection = 0;
231 232
    }

233 234 235
  /***********************************************/
  /* Second try : the OSSIM projection factory   */
  /***********************************************/
236 237 238 239
  if (!hasMetaData)
    {
    ossimImageHandler* handler = ossimImageHandlerRegistry::instance()
                                 ->open(ossimFilename(filename.c_str()));
240
    if (handler)
241 242
      {
      otbMsgDevMacro(<< "OSSIM Open Image SUCCESS ! ");
243

244 245 246
      // Add ossimPlugins model
      ossimProjectionFactoryRegistry::instance()->registerFactory(ossimplugins::ossimPluginProjectionFactory::instance());

247 248
      ossimRefPtr<ossimImageGeometry> geom = handler->getImageGeometry();
      if (geom.valid())
249
        {
250
        projection = geom->getProjection();
251 252 253 254
        if (projection)
          {
          hasMetaData = projection->saveState(geom_kwl);
          }
255
        }
256
      
257 258
      // if the handler has found a sensor model, copy the tags found
      if (hasMetaData && dynamic_cast<ossimSensorModel*>(projection))
259
        {
260 261
        otbMsgDevMacro(<<"OSSIM sensor projection instantiated ! ");
        otb_kwl.SetKeywordlist(geom_kwl);
262
        }
263
      else
264
        {
265 266 267 268 269 270 271 272 273 274 275 276
        hasMetaData = false;
        }
      // Free memory
      delete handler;
      }
    }
  
  /**********************************************************/
  /* Third try : look for external geom file and RPC tags   */
  /**********************************************************/
  if (!hasMetaData)
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
277
    // If still no metadata, try the ".geom" file
278 279 280 281
    ossimFilename ossimGeomFile = ossimFilename(filename).setExtension(".geom");
    otb_kwl = ReadGeometryFromGEOMFile(ossimGeomFile);
    
    // also check any RPC tags
282 283 284 285 286
    ImageKeywordlist rpc_kwl;
    if (checkRpcTag)
      {
      rpc_kwl = ReadGeometryFromRPCTag(filename);
      }
287 288 289 290 291 292 293 294 295 296 297 298 299 300
     
    if (otb_kwl.HasKey("type"))
      {
      // external geom has a "type" keyword
      std::string geomType(otb_kwl.GetMetadataByKey("type"));
      ossimProjection *testProj = ossimProjectionFactoryRegistry::instance()
        ->createProjection(ossimString(geomType));
      if (dynamic_cast<ossimSensorModel*>(testProj))
        {
        // "type" keyword corresponds to a sensor model : don't erase it
        if (rpc_kwl.GetSize() > 0)
          {
          rpc_kwl.ClearMetadataByKey("type");
          }
301
        
302 303 304 305 306
        ossimKeywordlist ossim_test_kwl;
        otb_kwl.convertToOSSIMKeywordlist(ossim_test_kwl);
        testProj = ossimProjectionFactoryRegistry::instance()
          ->createProjection(ossim_test_kwl);
        if (testProj)
307
          {
308 309
          // external geom contains a valid sensor geometry
          hasMetaData = true;
310 311
          }
        }
312 313 314 315 316 317 318 319 320
      }
      
    // copy keywords found in RPC tags if the external geom is not valid
    if (!hasMetaData && rpc_kwl.GetSize() > 0)
      {
      const ImageKeywordlist::KeywordlistMap &kwlMap = rpc_kwl.GetKeywordlist();
      for (ImageKeywordlist::KeywordlistMap::const_iterator it = kwlMap.begin();
          it != kwlMap.end();
          ++it)
321
        {
322 323 324 325
        if (it->second != "")
          {
          otb_kwl.AddKey(it->first , it->second);
          }
326
        }
327
      hasMetaData = true;
328 329 330
      }
    }

331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347
  // In some corner cases, ossim is unable to create a ossimImageHandler
  // instance for the provided file.
  // This can happen for TIFF+geom files where the ossimTiffTileSource
  // fails when trying to open the TIFF file (for example
  // because ossim is compiled with a libtiff without BigTIFF support
  // and when the actual TIFF file exceeds the Large File limit).
  // In such cases, we don't actually need ossim to interpret the TIFF file,
  // and just want a valid ossimKeywordlist corresponding to a sensor model.

  // Here is implemented a shortcut used as a fallback scenario,
  // where we bypass the ossimImageHandlerRegistry
  // and directly create a ossimKeywordlist from the ".geom" file.

  // We then verify it is a valid sensor model by using otb::SensorModelAdapter
  // which uses ossimSensorModelFactory and ossimPluginProjectionFactory internally,
  // thus by-passing the need for a valid ossimImageHandler.
  
348 349 350 351 352 353 354
  if (!hasMetaData)
    {
    otbMsgDevMacro(<< "OSSIM MetaData not present ! ");
    }
  else
    {
    otbMsgDevMacro(<< "OSSIM MetaData present ! ");
355
    //otbMsgDevMacro(<< geom_kwl);
356 357 358 359 360
    }

  return otb_kwl;
}

361 362 363 364
ImageKeywordlist
ReadGeometryFromGEOMFile(const std::string& filename)
{
  ossimKeywordlist geom_kwl;
365
  ImageKeywordlist otb_kwl;
366
  ossimFilename ossimGeomFile = ossimFilename(filename);
367 368 369 370 371 372 373 374 375

  if (ossimGeomFile.exists() && ossimGeomFile.isFile())
    {
    // Interpret the geom file as a KWL
    ossimKeywordlist kwl(ossimGeomFile);

    // Check that the geom file results in a valid ossimKeywordlist
    if (kwl.getErrorStatus() == ossimErrorCodes::OSSIM_OK)
      {
376 377 378 379 380 381
//       // Be sure there is a corresponding instance of ossimSensorModel
//       // which understands this kwl
//       SensorModelAdapter::Pointer sensorModel = SensorModelAdapter::New();
//       ImageKeywordlist otbkwl;
//       otbkwl.SetKeywordlist(kwl);
//       sensorModel->CreateProjection(otbkwl);
OTB Bot's avatar
STYLE  
OTB Bot committed
382
//
383 384 385 386 387
//       if (sensorModel->IsValidSensorModel())
//         {
//         geom_kwl = kwl;
//         }
      geom_kwl = kwl;
388 389 390 391 392 393 394 395
      }
    }

  otb_kwl.SetKeywordlist(geom_kwl);

  return otb_kwl;
}

396 397 398 399 400 401 402 403
ImageKeywordlist
ReadGeometryFromRPCTag(const std::string& filename)
{
  ossimKeywordlist geom_kwl;
  ImageKeywordlist otb_kwl;
  
  //  try to use GeoTiff RPC tag if present.
  // Warning : RPC in subdatasets are not supported
404
  GDALDriverH identifyDriverH = GDALIdentifyDriver(filename.c_str(), NULL);
405 406 407 408 409 410
  if(identifyDriverH == NULL)
    {
    // If no driver has identified the dataset, don't try to open it and exit
    return otb_kwl;
    }
  
411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500
  GDALDatasetH datasetH = GDALOpen(filename.c_str(), GA_ReadOnly);
  if (datasetH != NULL)
    {
    GDALDataset* dataset = static_cast<GDALDataset*>(datasetH);
    GDALRPCInfo rpcStruct;
    if (GDALExtractRPCInfo(dataset->GetMetadata("RPC"),&rpcStruct))
      {
      otbMsgDevMacro(<<"RPC Coeff found");
      std::vector<double> lineNumCoefs;
      std::vector<double> lineDenCoefs;
      std::vector<double> sampNumCoefs;
      std::vector<double> sampDenCoefs;
      
      for (unsigned int k=0; k<20; ++k)
        {
        lineNumCoefs.push_back(rpcStruct.adfLINE_NUM_COEFF[k]);
        lineDenCoefs.push_back(rpcStruct.adfLINE_DEN_COEFF[k]);
        sampNumCoefs.push_back(rpcStruct.adfSAMP_NUM_COEFF[k]);
        sampDenCoefs.push_back(rpcStruct.adfSAMP_DEN_COEFF[k]);
        }
      
      ossimRefPtr<ossimRpcModel> rpcModel = new ossimRpcModel;
      rpcModel->setAttributes( rpcStruct.dfSAMP_OFF,
                              rpcStruct.dfLINE_OFF,
                              rpcStruct.dfSAMP_SCALE,
                              rpcStruct.dfLINE_SCALE,
                              rpcStruct.dfLAT_OFF,
                              rpcStruct.dfLONG_OFF,
                              rpcStruct.dfHEIGHT_OFF,
                              rpcStruct.dfLAT_SCALE,
                              rpcStruct.dfLONG_SCALE,
                              rpcStruct.dfHEIGHT_SCALE,
                              sampNumCoefs,
                              sampDenCoefs,
                              lineNumCoefs,
                              lineDenCoefs);
      
      double errorBias = 0.0;
      double errorRand = 0.0;
      
      // setup other metadata
      rpcModel->setPositionError(errorBias,errorRand,true);
      ossimDrect rectangle(0.0,
                          0.0,
                          static_cast<double>(dataset->GetRasterXSize()-1),
                          static_cast<double>(dataset->GetRasterYSize()-1));
      rpcModel->setImageRect(rectangle);
      
      ossimDpt size;
      size.line = rectangle.height();
      size.samp = rectangle.width();
      rpcModel->setImageSize(size);
      
      // Compute 4 corners and reference point
      rpcModel->updateModel();
      double heightOffset = rpcStruct.dfHEIGHT_OFF;
      ossimGpt ulGpt, urGpt, lrGpt, llGpt;
      ossimGpt refGndPt;
      
      rpcModel->lineSampleHeightToWorld(rectangle.ul(), heightOffset, ulGpt);
      rpcModel->lineSampleHeightToWorld(rectangle.ur(), heightOffset, urGpt);
      rpcModel->lineSampleHeightToWorld(rectangle.lr(), heightOffset, lrGpt);
      rpcModel->lineSampleHeightToWorld(rectangle.ll(), heightOffset, llGpt);
      rpcModel->setGroundRect(ulGpt,urGpt,lrGpt,llGpt);
      
      rpcModel->lineSampleHeightToWorld(rectangle.midPoint(), heightOffset, refGndPt);
      rpcModel->setRefGndPt(refGndPt);
      
      // compute ground sampling distance
      try
        {
        // Method can throw ossimException.
        rpcModel->computeGsd();
        }
      catch (const ossimException& e)
        {
        otbMsgDevMacro(<< "OSSIM Compute ground sampling distance FAILED ! ");
        }
      
      if (rpcModel->saveState(geom_kwl))
        {
        otb_kwl.SetKeywordlist(geom_kwl);
        }
      }
    GDALClose(datasetH);
    }
  
  return otb_kwl;
}

501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517
void
WriteGeometry(const ImageKeywordlist& otb_kwl, const std::string& filename)
{
  // Write the image keyword list if any
  ossimKeywordlist geom_kwl;
  otb_kwl.convertToOSSIMKeywordlist(geom_kwl);

  if (geom_kwl.getSize() > 0)
    {
    otbMsgDevMacro(<< "Exporting keywordlist ...");
    ossimFilename geomFileName(filename);
    geomFileName.setExtension(".geom");
    geom_kwl.write(geomFileName.chars());
    }
}


518
} //namespace otb