Commit 4d1a5660 authored by Emmanuel Christophe's avatar Emmanuel Christophe

ENH: improve TileMapTransform interface and usage

parent d57a7a74
......@@ -18,17 +18,10 @@
#ifndef __otbTileMapTransform_h
#define __otbTileMapTransform_h
#include <iostream>
#include <sstream>
#include <stdio.h>
#include "itkTransform.h"
#include "itkExceptionObject.h"
#include "itkMacro.h"
#include "otbMapProjection.h"
namespace ossimplugins {
class ossimTileMapModel;
}
// Only for the enum definition
#include "otbGenericMapProjection.h"
namespace otb
{
......@@ -54,7 +47,6 @@ public:
typedef itk::SmartPointer<const Self> ConstPointer;
typedef typename Superclass::ScalarType ScalarType;
typedef ossimplugins::ossimTileMapModel OssimTileMapTransformType;
typedef itk::Point<ScalarType, NInputDimensions> InputPointType;
typedef itk::Point<ScalarType, NOutputDimensions> OutputPointType;
......@@ -71,11 +63,13 @@ public:
itkStaticConstMacro(SpaceDimension, unsigned int, NInputDimensions);
itkStaticConstMacro(ParametersDimension, unsigned int, NInputDimensions * (NInputDimensions + 1));
itkGetConstMacro(Depth, int);
itkSetMacro(Depth, int);
void SetLevel(unsigned int level);
unsigned int GetLevel() const;
OutputPointType TransformPoint(const InputPointType& point) const;
virtual InputPointType Origin();
virtual void PrintMap() const;
......@@ -88,11 +82,12 @@ public:
protected:
TileMapTransform();
virtual ~TileMapTransform();
OssimTileMapTransformType* m_TileMapTransform;
private:
TileMapTransform(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
int m_Depth;
};
} // namespace otb
......
......@@ -21,28 +21,20 @@
#include "otbTileMapTransform.h"
#include "otbMacro.h"
#include "ossimTileMapModel.h"
#include "base/ossimGpt.h"
#include "base/ossimDpt.h"
namespace otb
{
template<TransformDirection::TransformationDirection TTransformDirection, class TScalarType,
unsigned int NInputDimensions, unsigned int NOutputDimensions>
TileMapTransform<TTransformDirection, TScalarType, NInputDimensions, NOutputDimensions>
::TileMapTransform() : Superclass(SpaceDimension, ParametersDimension)
{
m_TileMapTransform = new OssimTileMapTransformType();
}
::TileMapTransform() : Superclass(SpaceDimension, ParametersDimension), m_Depth(0)
{}
template<TransformDirection::TransformationDirection TTransformDirection, class TScalarType,
unsigned int NInputDimensions, unsigned int NOutputDimensions>
TileMapTransform<TTransformDirection, TScalarType, NInputDimensions, NOutputDimensions>
::~TileMapTransform()
{
delete m_TileMapTransform;
}
{}
template<TransformDirection::TransformationDirection TTransformDirection, class TScalarType,
unsigned int NInputDimensions, unsigned int NOutputDimensions>
......@@ -54,73 +46,38 @@ TileMapTransform<TTransformDirection, TScalarType, NInputDimensions, NOutputDime
if (DirectionOfMapping == TransformDirection::INVERSE)
{
// otbMsgDevMacro(<< "Cartographic coordinates: (" << point[0] << "," << point[1] << ")");
//from "itk::point" to "ossim::ossimDpt"
ossimDpt ossimDPoint(point[0], point[1]);
//map projection
ossimGpt ossimGPoint;
// ossimGPoint=m_TileMapTransform->inverse(ossimDPoint);
m_TileMapTransform->lineSampleToWorld(ossimDPoint, ossimGPoint);
// otbGenericMsgDebugMacro(<< "Inverse : " << std::endl << m_TileMapTransform->print(std::cout));
outputPoint[0] = ossimGPoint.lon;
outputPoint[1] = ossimGPoint.lat;
// otbMsgDevMacro(<< "Geographic coordinates (long/lat) : (" << outputPoint[0] << "," << outputPoint[1] << ")");
outputPoint[0] = static_cast<double>(point[0])/(1 << m_Depth)/256 *360.0-180.0;
double y = static_cast<double>(point[1])/(1 << m_Depth)/256;
double ex = exp(4*M_PI*(y-0.5));
outputPoint[1] = -180.0/M_PI*asin((ex-1)/(ex+1));
}
if (DirectionOfMapping == TransformDirection::FORWARD)
{
// otbMsgDevMacro(<< "Geographic coordinates (long/lat) : (" << point[1] << "," << point[0] << ")");
//from "itk::point" to "ossim::ossimGpt"
ossimGpt ossimGPoint(point[1], point[0]);
//map projection
ossimDpt ossimDPoint;
// ossimDPoint=m_TileMapTransform->forward(ossimGPoint);
m_TileMapTransform->worldToLineSample(ossimGPoint, ossimDPoint);
// otbGenericMsgDebugMacro(<< "Forward : ========================= \n"
// << m_TileMapTransform->print(std::cout));
outputPoint[0] = ossimDPoint.x;
outputPoint[1] = ossimDPoint.y;
// otbMsgDevMacro(<< "Cartographic coordinates: (" << outputPoint[0] << "," << outputPoint[1] << ")");
double x = (180.0 + point[0]) / 360.0;
double y = - point[1] * M_PI / 180; // convert to radians
y = 0.5 * log((1+sin(y)) / (1 - sin(y)));
y *= 1.0/(2 * M_PI); // scale factor from radians to normalized
y += 0.5; // and make y range from 0 - 1
outputPoint[0] = floor(x*pow(2.,static_cast<double>(m_Depth))*256);
outputPoint[1] = floor(y*pow(2.,static_cast<double>(m_Depth))*256);
}
return outputPoint;
}
///\return The geographic point corresponding to (0, 0)
template<TransformDirection::TransformationDirection TTransformDirection, class TScalarType,
unsigned int NInputDimensions, unsigned int NOutputDimensions>
typename TileMapTransform<TTransformDirection, TScalarType, NInputDimensions, NOutputDimensions>::InputPointType
TileMapTransform<TTransformDirection, TScalarType, NInputDimensions, NOutputDimensions>
::Origin()
{
ossimGpt ossimOrigin = m_TileMapTransform->origin();
InputPointType otbOrigin;
otbOrigin[0] = ossimOrigin.lat;
otbOrigin[1] = ossimOrigin.lon;
return otbOrigin;
}
template<TransformDirection::TransformationDirection TTransformDirection, class TScalarType,
unsigned int NInputDimensions, unsigned int NOutputDimensions>
void
TileMapTransform<TTransformDirection, TScalarType, NInputDimensions, NOutputDimensions>
::PrintMap() const
{
std::cout << m_TileMapTransform->print(std::cout);
}
{}
template<TransformDirection::TransformationDirection TTransformDirection, class TScalarType,
unsigned int NInputDimensions, unsigned int NOutputDimensions>
void TileMapTransform<TTransformDirection, TScalarType, NInputDimensions, NOutputDimensions>
::SetLevel(unsigned int level)
{
m_TileMapTransform->setDepth(level);
this->SetDepth(level);
}
template<TransformDirection::TransformationDirection TTransformDirection, class TScalarType,
......@@ -128,7 +85,7 @@ template<TransformDirection::TransformationDirection TTransformDirection, class
unsigned int TileMapTransform<TTransformDirection, TScalarType, NInputDimensions, NOutputDimensions>
::GetLevel() const
{
return m_TileMapTransform->getDepth();
return this->GetDepth();
}
} // namespace otb
......
......@@ -59,7 +59,7 @@
#include "otbForwardSensorModel.h"
#include "otbExtractROI.h"
#include "otbImageFileWriter.h"
#include "ossimTileMapModel.h"
#include "otbTileMapTransform.h"
#include "otbWorldFile.h"
// Software Guide : EndCodeSnippet
......@@ -137,16 +137,10 @@ int main(int argc, char* argv[])
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef otb::InverseSensorModel<double> ModelType;
ModelType::Pointer model = ModelType::New();
bool resInvModel = model->SetImageGeometry(readerTile->GetOutput()->GetImageKeywordlist());
dynamic_cast<ossimplugins::ossimTileMapModel*>(model->GetOssimModel())->setDepth(depth);
if (!resInvModel)
{
std::cerr << "Unable to create a model" << std::endl;
return 1;
}
typedef otb::TileMapTransform<otb::TransformDirection::FORWARD> TransformType;
TransformType::Pointer transform = TransformType::New();
transform->SetDepth(depth);
typedef itk::Point <double, 2> PointType;
PointType lonLatPoint;
......@@ -154,7 +148,7 @@ int main(int argc, char* argv[])
lonLatPoint[1] = lat;
PointType tilePoint;
tilePoint = model->TransformPoint(lonLatPoint);
tilePoint = transform->TransformPoint(lonLatPoint);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
......@@ -204,33 +198,28 @@ int main(int argc, char* argv[])
// new image in a GIS system. For this, we need to compute the coordinates
// of the top left corner and the spacing in latitude and longitude.
//
// For that, we use the forward model to convert the corner coordinates into
// For that, we use the inverse transform to convert the corner coordinates into
// latitude and longitude.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef otb::ForwardSensorModel<double> ForwardModelType;
ForwardModelType::Pointer modelForward = ForwardModelType::New();
bool resFwdModel = modelForward->SetImageGeometry(readerTile->GetOutput()->GetImageKeywordlist());
dynamic_cast<ossimplugins::ossimTileMapModel*>(modelForward->GetOssimModel())->setDepth(
depth);
if (!resFwdModel)
{
std::cerr << "Unable to create a forward model" << std::endl;
return 1;
}
typedef otb::TileMapTransform<otb::TransformDirection::INVERSE> InverseTransformType;
InverseTransformType::Pointer transformInverse = InverseTransformType::New();
transformInverse->SetDepth(depth);
double lonUL, latUL, lonSpacing, latSpacing;
tilePoint[0] = startX - sizeX / 2;
tilePoint[1] = startY - sizeY / 2;
lonLatPoint = modelForward->TransformPoint(tilePoint);
lonLatPoint = transformInverse->TransformPoint(tilePoint);
lonUL = lonLatPoint[0];
latUL = lonLatPoint[1];
tilePoint[0] = startX + sizeX / 2;
tilePoint[1] = startY + sizeY / 2;
lonLatPoint = modelForward->TransformPoint(tilePoint);
lonLatPoint = transformInverse->TransformPoint(tilePoint);
lonSpacing = (lonLatPoint[0] - lonUL) / (sizeX - 1);
latSpacing = (lonLatPoint[1] - latUL) / (sizeY - 1);
// Software Guide : EndCodeSnippet
......
......@@ -16,11 +16,12 @@
=========================================================================*/
#include <fstream>
#include <iomanip>
#include "otbTileMapTransform.h"
int otbTileMapTransform(int argc, char* argv[])
{
const char * outFileName = argv[1];
std::ofstream file;
file.open(outFileName);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment