Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
N
NORMLIM_sigma0
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Package Registry
Container Registry
Model registry
Operate
Environments
Terraform modules
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
s1-tiling
NORMLIM_sigma0
Merge requests
!4
Resolve "Migrate code to OTB 8.x"
Code
Review changes
Check out branch
Download
Patches
Plain diff
Merged
Resolve "Migrate code to OTB 8.x"
5-migrate-code-to-otb-8-x
into
master
Overview
0
Commits
6
Pipelines
0
Changes
1
Merged
Luc Hermitte
requested to merge
5-migrate-code-to-otb-8-x
into
master
2 years ago
Overview
0
Commits
6
Pipelines
0
Changes
1
Expand
Closes
#5 (closed)
0
0
Merge request reports
Viewing commit
3d892f9c
Prev
Next
Show latest version
1 file
+
15
−
1
Inline
Compare changes
Side-by-side
Inline
Show whitespace changes
Show one file at a time
3d892f9c
ENH: Trace performances benchmarks
· 3d892f9c
Luc Hermitte
authored
1 year ago
include/otbSARDEMProjectionImageFilter2.hxx
0 → 100644
+
436
−
0
Options
/*
* Copyright (C) 2005-2023 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*
* https://www.orfeo-toolbox.org/
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
// This file is a fork of DiapOTB SARDEMProjection adapted to OTB 8 source code
// for the needs of S1Tiling
#ifndef otbSARDEMProjectionImageFilter2_hxx
#define otbSARDEMProjectionImageFilter2_hxx
#include
"otbSARDEMProjectionImageFilter2.h"
#include
"itkImageScanlineConstIterator.h"
#include
"itkImageScanlineIterator.h"
#include
"itkProgressReporter.h"
#include
"itkNumericTraitsPointPixel.h"
#if OTB_VERSION_MAJOR >= 8
# include "otbGeocentricTransform.h"
#else
# if defined(__GNUC__) || defined(__clang__)
# pragma GCC diagnostic push
# pragma GCC diagnostic ignored "-Wunused-parameter"
# include "ossim/base/ossimFilename.h"
# include "ossim/base/ossimGpt.h"
# pragma GCC diagnostic pop
# else
# include "ossim/base/ossimFilename.h"
# include "ossim/base/ossimGpt.h"
# endif
#endif
#include
"otbConfigurationManager.h"
#include
"otbDEMHandler.h"
#include
<memory>
#include
<boost/optional.hpp>
#include
<mutex>
#include
<thread>
#include
<cassert>
namespace
otb
{
/**
* Constructor with default initialization
*/
template
<
class
TImageIn
,
class
TImageOut
>
SARDEMProjectionImageFilter2
<
TImageIn
,
TImageOut
>::
SARDEMProjectionImageFilter2
(
std
::
string
const
&
geoid_filename
)
{
#if OTB_VERSION_MAJOR >= 8
auto
&
demHandler
=
DEMHandler
::
GetInstance
();
assert
(
demHandler
.
GetGeoidFile
()
==
geoid_filename
);
#else
auto
&
demHandler
=
*
DEMHandler
::
Instance
();
if
(
!
geoid_filename
.
empty
())
{
const
std
::
string
isEmg96
=
"egm96"
;
std
::
size_t
found
=
geoid_filename
.
find
(
isEmg96
);
if
(
found
!=
std
::
string
::
npos
)
{
m_geoidEmg96
=
std
::
make_unique
<
ossimGeoidEgm96
>
(
ossimFilename
(
geoid_filename
));
}
}
#endif
}
/**
* Set Sar Image keyWordList
*/
#if OTB_VERSION_MAJOR >= 8
template
<
class
TImageIn
,
class
TImageOut
>
void
SARDEMProjectionImageFilter2
<
TImageIn
,
TImageOut
>
::
SetSARImageMetadata
(
ImageMetadata
sarImageMetadata
)
{
m_SarImageMetadata
=
sarImageMetadata
;
}
#else
template
<
class
TImageIn
,
class
TImageOut
>
void
SARDEMProjectionImageFilter2
<
TImageIn
,
TImageOut
>
::
SetSARImageKeyWorList
(
ImageKeywordlist
sarImageKWL
)
{
m_SarImageKwl
=
sarImageKWL
;
}
#endif
/**
* Method GenerateOutputInformaton()
*/
template
<
class
TImageIn
,
class
TImageOut
>
void
SARDEMProjectionImageFilter2
<
TImageIn
,
TImageOut
>
::
GenerateOutputInformation
()
{
// Call superclass implementation
Superclass
::
GenerateOutputInformation
();
// Get pointers to the input and output
ImageInConstPointer
inputPtr
=
this
->
GetInput
();
ImageOutPointer
outputPtr
=
this
->
GetOutput
();
// Image KeyWordList
#if OTB_VERSION_MAJOR >= 8
auto
inputKwl
=
inputPtr
->
GetImageMetadata
();
#else
ImageKeywordlist
inputKwl
=
inputPtr
->
GetImageKeywordlist
();
#endif
// Set the number of Components for the Output VectorImage
// At Least 4 Components :
// _ C : Colunm of DEM into SAR Image
// _ L : Line of DEM into SAR Image
// _ Z : Coordonates Z
// _ Y : Coordonates Y
// (Z,Y) coordinates, defined into ossimSarSensorMode (if otb version < 8) or SarSensorMode as follows:
// Let n = |sensorPos|,
// ps2 = scalar_product(sensorPos,worldPoint)
// d = distance(sensorPos,worldPoint)
//
// Z = n - ps2/n
// Y = sqrt(d*d - Z*Z)
//
// m_withH = true => 1 component added H : high
// m_withXYZ = true => 3 Three components added Xcart Ycart Zcart : cartesien coordonates
// m_withSatPos = true => 3 Three components added : Satellite Position into cartesian
#if OTB_VERSION_MAJOR >= 8
inputKwl
.
Add
(
std
::
string
(
"band.L"
),
std
::
to_string
(
1
));
inputKwl
.
Add
(
std
::
string
(
"band.C"
),
std
::
to_string
(
0
));
inputKwl
.
Add
(
std
::
string
(
"band.Z"
),
std
::
to_string
(
2
));
inputKwl
.
Add
(
std
::
string
(
"band.Y"
),
std
::
to_string
(
3
));
#else
// Fill KeyWordList
inputKwl
.
AddKey
(
"support_data.band.L"
,
"1"
);
inputKwl
.
AddKey
(
"support_data.band.C"
,
"0"
);
inputKwl
.
AddKey
(
"support_data.band.Z"
,
"2"
);
inputKwl
.
AddKey
(
"support_data.band.Y"
,
"3"
);
#endif
m_NbComponents
=
4
;
if
(
m_withXYZ
)
{
m_NbComponents
+=
3
;
#if OTB_VERSION_MAJOR >= 8
inputKwl
.
Add
(
"band.XCart"
,
"4"
);
inputKwl
.
Add
(
"band.YCart"
,
"5"
);
inputKwl
.
Add
(
"band.ZCart"
,
"6"
);
#else
inputKwl
.
AddKey
(
"support_data.band.XCart"
,
"4"
);
inputKwl
.
AddKey
(
"support_data.band.YCart"
,
"5"
);
inputKwl
.
AddKey
(
"support_data.band.ZCart"
,
"6"
);
#endif
}
if
(
m_withH
)
{
m_indH
=
m_NbComponents
;
m_NbComponents
+=
1
;
#if OTB_VERSION_MAJOR >= 8
inputKwl
.
Add
(
"band.H"
,
std
::
to_string
(
m_indH
));
#else
inputKwl
.
AddKey
(
"support_data.band.H"
,
std
::
to_string
(
m_indH
));
#endif
}
if
(
m_withSatPos
)
{
m_indSatPos
=
m_NbComponents
;
m_NbComponents
+=
3
;
#if OTB_VERSION_MAJOR >= 8
inputKwl
.
Add
(
"band.XSatPos"
,
std
::
to_string
(
m_indSatPos
));
inputKwl
.
Add
(
"band.YSatPos"
,
std
::
to_string
(
m_indSatPos
+
1
));
inputKwl
.
Add
(
"band.ZSatPos"
,
std
::
to_string
(
m_indSatPos
+
2
));
#else
inputKwl
.
AddKey
(
"support_data.band.XSatPos"
,
std
::
to_string
(
m_indSatPos
));
inputKwl
.
AddKey
(
"support_data.band.YSatPos"
,
std
::
to_string
(
m_indSatPos
+
1
));
inputKwl
.
AddKey
(
"support_data.band.ZSatPos"
,
std
::
to_string
(
m_indSatPos
+
2
));
#endif
}
// Add gain and direction into kwl
#if OTB_VERSION_MAJOR >= 8
inputKwl
.
Add
(
"band.DirectionToScanDEMC"
,
std
::
to_string
(
m_directionToScanDEMC
));
inputKwl
.
Add
(
"band.DirectionToScanDEML"
,
std
::
to_string
(
m_directionToScanDEML
));
inputKwl
.
Add
(
"band.Gain"
,
std
::
to_string
(
m_gain
));
#else
inputKwl
.
AddKey
(
"support_data.band.DirectionToScanDEMC"
,
std
::
to_string
(
m_directionToScanDEMC
));
inputKwl
.
AddKey
(
"support_data.band.DirectionToScanDEML"
,
std
::
to_string
(
m_directionToScanDEML
));
inputKwl
.
AddKey
(
"support_data.band.Gain"
,
std
::
to_string
(
m_gain
));
#endif
outputPtr
->
SetNumberOfComponentsPerPixel
(
m_NbComponents
);
#if OTB_VERSION_MAJOR >= 8
// Create and Initilaze the SarSensorModel
m_SarSensorModel
=
std
::
make_unique
<
SarSensorModel
>
(
m_SarImageMetadata
);
#else
// Create and Initilaze the SarSensorModelAdapter
m_SarSensorModel
=
SarSensorModelAdapter
::
New
();
bool
loadOk
=
m_SarSensorModel
->
LoadState
(
m_SarImageKwl
);
if
(
!
loadOk
||
!
m_SarSensorModel
->
IsValidSensorModel
())
{
itkExceptionMacro
(
<<
"SAR image does not contain a valid SAR sensor model."
);
}
#endif
// Compute the output image size, spacing and origin (same as the input)
const
ImageInSizeType
&
inputSize
=
inputPtr
->
GetLargestPossibleRegion
().
GetSize
();
ImageInPointType
origin
=
inputPtr
->
GetOrigin
();
const
ImageInSpacingType
&
inSP
=
inputPtr
->
GetSpacing
();
ImageOutRegionType
outputLargestPossibleRegion
=
inputPtr
->
GetLargestPossibleRegion
();
outputLargestPossibleRegion
.
SetSize
(
static_cast
<
ImageOutSizeType
>
(
inputSize
));
outputPtr
->
SetLargestPossibleRegion
(
outputLargestPossibleRegion
);
outputPtr
->
SetOrigin
(
static_cast
<
ImageOutPointType
>
(
origin
));
outputPtr
->
SetSpacing
(
static_cast
<
ImageOutSpacingType
>
(
inSP
));
// Set Image KeyWordList
#if OTB_VERSION_MAJOR >= 8
outputPtr
->
SetImageMetadata
(
inputKwl
);
#else
outputPtr
->
SetImageKeywordList
(
inputKwl
);
#endif
}
/**
* Method GenerateInputRequestedRegion
*/
template
<
class
TImageIn
,
class
TImageOut
>
void
SARDEMProjectionImageFilter2
<
TImageIn
,
TImageOut
>
::
GenerateInputRequestedRegion
()
{
ImageOutRegionType
outputRequestedRegion
=
this
->
GetOutput
()
->
GetRequestedRegion
();
ImageInRegionType
inputRequestedRegion
=
static_cast
<
ImageInRegionType
>
(
outputRequestedRegion
);
ImageInPointer
inputPtr
=
const_cast
<
ImageInType
*
>
(
this
->
GetInput
()
);
inputPtr
->
SetRequestedRegion
(
inputRequestedRegion
);
}
/**
* Method ThreadedGenerateData
*/
template
<
class
TImageIn
,
class
TImageOut
>
void
SARDEMProjectionImageFilter2
<
TImageIn
,
TImageOut
>
::
ThreadedGenerateData
(
ImageOutRegionType
const
&
outputRegionForThread
,
itk
::
ThreadIdType
threadId
)
{
// Compute corresponding input region (same as output here)
ImageInRegionType
inputRegionForThread
=
static_cast
<
ImageInRegionType
>
(
outputRegionForThread
);
// Define/declare an iterator that will walk the input region for this
// thread.
typedef
itk
::
ImageScanlineConstIterator
<
ImageInType
>
InputIterator
;
InputIterator
InIt
(
this
->
GetInput
(),
inputRegionForThread
);
// Define/declare an iterator that will walk the output region for this
// thread. NO VECTOR IMAGE INTO ITERATOR TEMPLATE
typedef
itk
::
ImageScanlineIterator
<
ImageOutType
>
OutputIterator
;
OutputIterator
OutIt
(
this
->
GetOutput
(),
outputRegionForThread
);
// Support progress methods/callbacks
itk
::
ProgressReporter
progress
(
this
,
threadId
,
outputRegionForThread
.
GetNumberOfPixels
()
);
// Transform all points
InIt
.
GoToBegin
();
OutIt
.
GoToBegin
();
Point3DType
demGeoPoint
;
Point2DType
demLatLonPoint
;
Point2DType
col_row
(
0
);
Point2DType
y_z
(
0
);
Point3DType
xyzCart
;
Point3DType
satPos
;
Point3DType
satVel
;
assert
(
m_NbComponents
==
4
+
(
m_withXYZ
?
3
:
0
)
+
(
m_withH
?
1
:
0
)
+
(
m_withSatPos
?
3
:
0
));
ImageOutPixelType
pixelOutput
;
pixelOutput
.
Reserve
(
m_NbComponents
);
#if OTB_VERSION_MAJOR >= 8
auto
&
demHandler
=
DEMHandler
::
GetInstance
();
// demHandler.ClearElevationParameters(); // NO!!!!!
auto
&
demHandlerTLS
=
demHandler
.
GetHandlerForCurrentThread
();
#else
auto
&
demHandler
=
*
DEMHandler
::
Instance
();
#endif
// For each Line in MNT
while
(
!
InIt
.
IsAtEnd
())
{
InIt
.
GoToBeginOfLine
();
OutIt
.
GoToBeginOfLine
();
// For each column in MNT
while
(
!
InIt
.
IsAtEndOfLine
())
{
// Check if NODATA
if
(
InIt
.
Get
()
!=
m_NoData
)
{
// Trnasform index to Lat/Lon Point
this
->
GetInput
()
->
TransformIndexToPhysicalPoint
(
InIt
.
GetIndex
(),
demLatLonPoint
);
demGeoPoint
[
0
]
=
demLatLonPoint
[
0
];
demGeoPoint
[
1
]
=
demLatLonPoint
[
1
];
// Get elevation from earth geoid thanks to DEMHandler or ossimgeoidEmg96
double
h
;
// Elevation from earth geoid
#if OTB_VERSION_MAJOR >= 8
// TODO: GetGeoidHeight has very bad performances (*): GDAL
// isn't meant to be used for one coordinate at the time.
// Instead: have a multi-pass algorithm
// 1. get in a SINGLE call the geoid height for all points in
// the buffer
// 2. Compute the rest
//
// (*)
// - 96.3% of execution time spent in ThreadedGenerateData
// - 61.3% in GetGeoidHeight
// - 31.3% GDALRasterBand::RasterIO
// - 11.7% in mutex loc+unlock
// - 29.2% OCRCoordinateTransform
// - 20.2% in get_pid
//
h
=
GetGeoidHeight
(
demHandlerTLS
,
demLatLonPoint
);
#else
if
(
m_geoidEmg96
)
{
ossimGpt
gptPt
;
gptPt
.
lon
=
demLatLonPoint
[
0
];
gptPt
.
lat
=
demLatLonPoint
[
1
];
h
=
m_geoidEmg96
->
offsetFromEllipsoid
(
gptPt
);
}
else
{
h
=
demHandler
.
GetHeightAboveEllipsoid
(
demGeoPoint
[
0
],
demGeoPoint
[
1
]);
}
#endif
// Correct the height with earth geoid
demGeoPoint
[
2
]
=
InIt
.
Get
()
+
h
;
// Call the method of sarSensorModel
m_SarSensorModel
->
WorldToLineSampleYZ
(
demGeoPoint
,
col_row
,
y_z
);
// Fill the four channels for output
pixelOutput
[
0
]
=
col_row
[
0
];
pixelOutput
[
1
]
=
col_row
[
1
];
pixelOutput
[
2
]
=
y_z
[
1
];
pixelOutput
[
3
]
=
y_z
[
0
];
// Add channels if required
if
(
m_withXYZ
)
{
assert
(
pixelOutput
.
Size
()
>=
6
);
#if OTB_VERSION_MAJOR >= 8
xyzCart
=
Projection
::
WorldToEcef
(
demGeoPoint
);
#else
SarSensorModelAdapter
::
WorldToCartesian
(
demGeoPoint
,
xyzCart
);
#endif
pixelOutput
[
4
]
=
xyzCart
[
0
];
pixelOutput
[
5
]
=
xyzCart
[
1
];
pixelOutput
[
6
]
=
xyzCart
[
2
];
}
if
(
m_withH
)
{
assert
(
pixelOutput
.
Size
()
>
m_indH
);
pixelOutput
[
m_indH
]
=
demGeoPoint
[
2
];
}
if
(
m_withSatPos
)
{
assert
(
pixelOutput
.
Size
()
>=
m_indSatPos
+
3
);
m_SarSensorModel
->
WorldToSatPositionAndVelocity
(
demGeoPoint
,
satPos
,
satVel
);
pixelOutput
[
m_indSatPos
]
=
satPos
[
0
];
pixelOutput
[
m_indSatPos
+
1
]
=
satPos
[
1
];
pixelOutput
[
m_indSatPos
+
2
]
=
satPos
[
2
];
}
// Set the output (without iterator because no VectorImage into iterator template)
OutIt
.
Set
(
pixelOutput
);
progress
.
CompletedPixel
();
}
else
{
// Fill all channels with NoData
for
(
int
id
=
0
;
id
<
m_NbComponents
;
id
++
)
{
pixelOutput
[
id
]
=
m_NoData
;
}
// Set the output (without itertor because no VectorImage into iterator template)
OutIt
.
Set
(
pixelOutput
);
progress
.
CompletedPixel
();
}
// Next Col
++
InIt
;
++
OutIt
;
}
// Next Line
InIt
.
NextLine
();
OutIt
.
NextLine
();
}
}
}
/*namespace otb*/
#endif
Loading