Skip to content
Snippets Groups Projects
Commit d7e3333f authored by Luc Hermitte's avatar Luc Hermitte
Browse files

PERF/WIP: First draft for 8-10x perf improvments 

Instead of iterating over all DEM lines again and again, a dichotomic search is
done as long as we know DEM lines progression is monotonic to the exploration.

TODO:

- test the case where lines are sorted in ascending order
- support the case where there is no monotonic progression
parent ee275279
No related branches found
No related tags found
1 merge request!1Resolve "Adapt DiapOTB SARCartesianMeanEstimation to keep shadows"
Hide whitespace changes
Inline Side-by-side
include/otbSARDEMPolygonsAnalysisImageFilter2.h
+ 6
− 7
View file @ d7e3333f
......@@ -87,6 +87,8 @@ public:
using ImageInPointType = typename ImageInType::PointType;
using ImageInSubPixelType = typename ImageInType::InternalPixelType;
static_assert(std::is_floating_point<ImageInSubPixelType>::value, "float or double!");
using FastPixelType = ImageInSubPixelType const*;
/** Typedef to describe the image index, size types and spacing for inout image. */
using ImageInIndexType = typename ImageInType::IndexType;
using ImageInIndexValueType = typename ImageInType::IndexValueType;
......@@ -222,10 +224,10 @@ protected:
* Scan all selected polygons and apply on it, the choosen functor.
*/
void PolygonsScan(const int ind_LineSAR,
std::vector<std::vector<ImageInPixelType> *> const& CLZY_InSideUp_Vec,
std::vector<std::vector<ImageInPixelType> *> const& CLZY_InSideDown_Vec,
std::vector<std::vector<ImageInPixelType> *> const& CLZY_OutSideUp_Vec,
std::vector<std::vector<ImageInPixelType> *> const& CLZY_OutSideDown_Vec,
std::vector<std::vector<FastPixelType> const*> const& CLZY_InSideUp_Vec,
std::vector<std::vector<FastPixelType> const*> const& CLZY_InSideDown_Vec,
std::vector<std::vector<FastPixelType> const*> const& CLZY_OutSideUp_Vec,
std::vector<std::vector<FastPixelType> const*> const& CLZY_OutSideDown_Vec,
int firstCol_into_outputRegion,
otb::Span<ImageOutPixelType> outValueTab, int threadId);
......@@ -271,9 +273,6 @@ private:
int m_NbLinesSAR;
int m_NbColSAR;
// Store optionnalImage
ImageOutPixelType * m_OptionnalResults;
// RSTransform (for inverse localisation)
RSTransformType2D::Pointer m_RSTransform;
};
......
include/otbSARDEMPolygonsAnalysisImageFilter2.hxx
+ 251
− 74
View file @ d7e3333f
......@@ -23,6 +23,7 @@
#include "otbSARDEMPolygonsAnalysisImageFilter2.h"
#include "equal_range_interval.h"
#include "otbVLVPointIterator.h"
#include "otbStringUtilities.h"
#include "otbVectorImage.h"
......@@ -36,9 +37,36 @@
#include <cmath>
#include <type_traits>
#include <cassert>
#include "smart-assert.h"
namespace otb
{
template <typename Cmp = std::less<void>>
struct CmpVLVComponent {
explicit CmpVLVComponent(unsigned comp) noexcept : m_comp(comp) {}
template <typename U>
bool operator()(U d, itk::VariableLengthVector<U> const& v) const noexcept
{ return Cmp{}(d, v[m_comp]); }
template <typename U>
bool operator()(itk::VariableLengthVector<U> const& v, U d) const noexcept
{ return Cmp{}(v[m_comp], d); }
template <typename U>
bool operator()(itk::VariableLengthVector<U> const& v1, itk::VariableLengthVector<U> const& v2) const noexcept
{ return Cmp{}(v1[m_comp], v2[m_comp]); }
private:
unsigned m_comp;
};
template <typename T>
inline
auto as_ptr(itk::VariableLengthVector<T> const& vlv)
{
return &vlv[0];
}
/**
* Set Sar Image keyWordList
*/
......@@ -497,8 +525,8 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF
const ImageInSubPixelType ligUR, const ImageInSubPixelType colUR,
const ImageInSubPixelType ligLR, const ImageInSubPixelType colLR,
const ImageInSubPixelType ligLL, const ImageInSubPixelType colLL,
int & i_InSideUp, int & i_InSideDown,
int & i_OutSideUp, int & i_OutSideDown)
int & i_InSideUp, int & i_InSideDown,
int & i_OutSideUp, int & i_OutSideDown)
{
auto const Lmax_ulur = std::max(ligUL, ligUR);
auto const Lmin_ulur = std::min(ligUL, ligUR);
......@@ -750,13 +778,15 @@ template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, clas
void
SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
::PolygonsScan(const int ind_LineSAR,
std::vector<std::vector<ImageInPixelType> *> const& CLZY_InSideUp_Vec,
std::vector<std::vector<ImageInPixelType> *> const& CLZY_InSideDown_Vec,
std::vector<std::vector<ImageInPixelType> *> const& CLZY_OutSideUp_Vec,
std::vector<std::vector<ImageInPixelType> *> const& CLZY_OutSideDown_Vec,
std::vector<std::vector<FastPixelType> const*> const& CLZY_InSideUp_Vec,
std::vector<std::vector<FastPixelType> const*> const& CLZY_InSideDown_Vec,
std::vector<std::vector<FastPixelType> const*> const& CLZY_OutSideUp_Vec,
std::vector<std::vector<FastPixelType> const*> const& CLZY_OutSideDown_Vec,
int firstCol_into_outputRegion, otb::Span<ImageOutPixelType> outValueTab,
int threadId)
{
ImageInConstPointer inputPtr = this->GetInput();
auto const nbcompo = inputPtr->GetNumberOfComponentsPerPixel();
bool isFirstPolygon = true;
// Elevation angle max (for shadows)
......@@ -803,11 +833,14 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF
// Select the CLY_InSideUp, CLY_InSideDown, CLY_OutSideUp and CLY_OutSideDown and
// estimate the wanted value (for example amplitude)
auto as_vlv = [nbcompo](FastPixelType p) {
return ImageInPixelType(p, nbcompo, false);
};
m_FunctionOnPolygon->estimate(ind_LineSAR,
(*CLZY_InSideUp_Vec [ind_LDEM])[ind_CDEM],
(*CLZY_InSideDown_Vec [ind_LDEM])[ind_CDEM],
(*CLZY_OutSideUp_Vec [ind_LDEM])[ind_CDEM],
(*CLZY_OutSideDown_Vec[ind_LDEM])[ind_CDEM],
as_vlv((*CLZY_InSideUp_Vec [ind_LDEM])[ind_CDEM]),
as_vlv((*CLZY_InSideDown_Vec [ind_LDEM])[ind_CDEM]),
as_vlv((*CLZY_OutSideUp_Vec [ind_LDEM])[ind_CDEM]),
as_vlv((*CLZY_OutSideDown_Vec[ind_LDEM])[ind_CDEM]),
firstCol_into_outputRegion, outValueTab, isFirstPolygon,
tgElevationMaxForCurrentLine, threadId);
......@@ -848,10 +881,10 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF
std::vector<ImageOutPixelType> outValueTab(nbCol);
// Allocate Vectors for polygon storage (vector of vectors to separate line and colunm of polygon into DEM)
std::vector<std::vector<ImageInPixelType> * > Polygon_SideInUp_Vec;
std::vector<std::vector<ImageInPixelType> * > Polygon_SideInDown_Vec;
std::vector<std::vector<ImageInPixelType> * > Polygon_SideOutUp_Vec;
std::vector<std::vector<ImageInPixelType> * > Polygon_SideOutDown_Vec;
std::vector<std::vector<FastPixelType> const*> Polygon_SideInUp_Vec;
std::vector<std::vector<FastPixelType> const*> Polygon_SideInDown_Vec;
std::vector<std::vector<FastPixelType> const*> Polygon_SideOutUp_Vec;
std::vector<std::vector<FastPixelType> const*> Polygon_SideOutDown_Vec;
// A polygon is composed of four lines
// => four sides (UL UR LR LL) : UL = current_point
......@@ -866,26 +899,30 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF
auto const& in_index = inputRegionForThread.GetIndex();
auto const in_sizeX = in_size[0];
auto const in_sizeY = in_size[1];
auto const in_startX = in_index[0];
auto const in_startY = in_index[1];
itk::IndexValueType const in_endX = in_startX + in_sizeX;
itk::IndexValueType const in_endY = in_startY + in_sizeY;
auto const in_full_line = in_sizeX * nBands;
// auto const in_startX = in_index[0];
// auto const in_startY = in_index[1];
// itk::IndexValueType const in_endX = in_startX + in_sizeX;
// itk::IndexValueType const in_endY = in_startY + in_sizeY;
// auto const in_full_line = in_sizeX * nBands;
ImageInSubPixelType const *in_data = input->GetBufferPointer();
unsigned long const in_line_off = nBands * input->GetBufferedRegion().GetSize(0);
in_data += nBands * input->ComputeOffset(in_index);
// Allocate Vectors for polygon storage for this line of DEM
std::vector<std::vector<ImageInPixelType>> Polygon_SideInUp_VecLine(nbDEMLines);
std::vector<std::vector<ImageInPixelType>> Polygon_SideInDown_VecLine(nbDEMLines);
std::vector<std::vector<ImageInPixelType>> Polygon_SideOutUp_VecLine(nbDEMLines);
std::vector<std::vector<ImageInPixelType>> Polygon_SideOutDown_VecLine(nbDEMLines);
std::vector<std::vector<FastPixelType>> Polygon_SideInUp_VecLine (nbDEMLines);
std::vector<std::vector<FastPixelType>> Polygon_SideInDown_VecLine (nbDEMLines);
std::vector<std::vector<FastPixelType>> Polygon_SideOutUp_VecLine (nbDEMLines);
std::vector<std::vector<FastPixelType>> Polygon_SideOutDown_VecLine(nbDEMLines);
// Memorize member variables so they can be optimized in loops
auto const margin = m_Margin;
auto const margin = m_Margin;
ImageInSubPixelType const nodata = m_NoData;
std::size_t cpt_in_window = 0;
std::size_t cpt_out_window = 0;
std::size_t cpt_lines_with_no_match = 0;
// For each line of output : construct polygon after polygon, determine if polygon intersects the current
// line and estimate the amplitude with the contribution of all selected polygon
for (OutIt.GoToBeginOfLine(); !OutIt.IsAtEnd() ; OutIt.NextLine())
......@@ -912,6 +949,12 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF
ImageInSubPixelType const lo_margin = ind_Line_Into_SARGeo - margin;
ImageInSubPixelType const hi_margin = ind_Line_Into_SARGeo + margin;
///////////////////////////// Polygon Construction and Check intersection /////////////////////////
// Using the input region associated to the output line doesn't make a
// big difference...
ImageOutRegionType outputRegionForLine(
OutIt.GetIndex(),
{outputRegionForThread.GetSize()[0], 1});
ImageInRegionType inputRegionForWindow = OutputRegionToInputRegion(outputRegionForLine);
int countDEMLine = 0;
......@@ -922,66 +965,196 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF
return VLVPointIterator<ImageInSubPixelType const>(in_data + nBands*dx + dy*in_line_off, nBands);
};
bool has_match = false;
int nb_up = 0;
int nb_dw = 0;
auto prev = (*point(0, dy))[1];
auto register_pixel = [&](std::size_t dx, std::size_t dy, auto const& pt_UL)
{
++cpt_in_window;
auto const pt_UR = *point(dx+1, dy );
auto const pt_LL = *point(dx , dy+1);
auto const pt_LR = *point(dx+1, dy+1);
auto const UL_col = pt_UL[0];
auto const UR_col = pt_UR[0];
auto const LR_col = pt_LR[0];
auto const LL_col = pt_LL[0];
// For each column
for (unsigned dx = 0 ; dx < in_sizeX -1 ; ++dx)
// Check if No Data
if (!(UL_col == nodata || UR_col == nodata || LR_col == nodata
|| LL_col == nodata))
{
// Determine if intersection
auto const UL_line = pt_UL[1];
auto const UR_line = pt_UR[1];
auto const LR_line = pt_LR[1];
auto const LL_line = pt_LL[1];
// Pixel of otb::VectorImage (input) : 4 Components : C, L, Z and Y
int i_InSideUp, i_InSideDown, i_OutSideUp, i_OutSideDown;
// With lig and col of UL, UR, LR and LL select iterators for InSide and Outside of
// each polygon (if intersection)
bool const hasIntersection = this->PolygonAndSarLineIntersection(
ind_Line_Into_SARGeo,
UL_line, UL_col,
UR_line, UR_col,
LR_line, LR_col,
LL_line, LL_col,
i_InSideUp,
i_InSideDown,
i_OutSideUp,
i_OutSideDown);
// If intersection == true : Retrive the intersected sides (two intersected among
// the four sides of the polygon)
if (hasIntersection)
{
FastPixelType tabInPt[] = { as_ptr(pt_UL), as_ptr(pt_UR), as_ptr(pt_LR), as_ptr(pt_LL)};
// Store the polygons
Polygon_SideInUp_VecLine [countDEMLine].push_back(tabInPt[i_InSideUp]);
Polygon_SideInDown_VecLine [countDEMLine].push_back(tabInPt[i_InSideDown]);
Polygon_SideOutUp_VecLine [countDEMLine].push_back(tabInPt[i_OutSideUp]);
Polygon_SideOutDown_VecLine[countDEMLine].push_back(tabInPt[i_OutSideDown]);
}
}
}; // ------------------------------[ register_pixel
// assume for now progression is always monotonic
// TODO: precheck the input window to fall back to the old input method
// otherwise
auto is_in_window = [lo_margin, hi_margin](ImageInSubPixelType line)
{
static_assert(std::is_same<typename std::remove_cv<decltype(lo_margin)>::type, ImageInSubPixelType>::value, "Types shall match!");
return lo_margin < line && line < hi_margin;
}; // ------------------------------[ is_in_window
auto line = [](VLVPointIterator<ImageInSubPixelType const> it)
{
auto pt_UL = point(dx, dy);
auto const UL_line = pt_UL[1];
return (*it)[1];
}; // ------------------------------[ line
// Protect Memory access and input access
static_assert(std::is_same<decltype(UL_line), decltype(lo_margin)>::value, "no convertion!!!");
if (UL_line > lo_margin && UL_line < hi_margin) // true 2% of the time...
bool is_line_monotonic = true;
if (is_line_monotonic)
{
auto find_lo_hi_window = [lo_margin, hi_margin, line, is_in_window](auto line_start, auto line_end)
{
auto const pt_UR = point(dx+1, dy );
auto const pt_LL = point(dx , dy+1);
auto const pt_LR = point(dx+1, dy+1);
VLVPointIterator<ImageInSubPixelType const> tabInPt[] = { pt_UL, pt_UR, pt_LR, pt_LL};
auto const UL_col = pt_UL[0];
auto const UR_col = pt_UR[0];
auto const LR_col = pt_LR[0];
auto const LL_col = pt_LL[0];
// Check if No Data
if (!(UL_col == nodata || UR_col == nodata || LR_col == nodata
|| LL_col == nodata))
if (line(line_start) < line(line_end-1))
{
// Determine if intersection
auto const UR_line = pt_UR[1];
auto const LR_line = pt_LR[1];
auto const LL_line = pt_LL[1];
// Pixel of otb::VectorImage (input) : 4 Components : C, L, Z and Y
int i_InSideUp, i_InSideDown, i_OutSideUp, i_OutSideDown;
// With lig and col of UL, UR, LR and LL select iterators for InSide and Outside of
// each polygon (if intersection)
bool const hasIntersection = this->PolygonAndSarLineIntersection(ind_Line_Into_SARGeo,
UL_line, UL_col,
UR_line, UR_col,
LR_line, LR_col,
LL_line, LL_col,
i_InSideUp,
i_InSideDown,
i_OutSideUp,
i_OutSideDown);
// If intersection == true : Retrive the intersected sides (two intersected among
// the four sides of the polygon)
if (hasIntersection)
auto lo_hi = equal_range_interval(
line_start, line_end,
lo_margin, hi_margin, CmpVLVComponent<>(1));
while (lo_hi.first != lo_hi.second && lo_hi.first != line_end && line(lo_hi.first) == lo_margin)
{
// Store the polygons
Polygon_SideInUp_VecLine [countDEMLine].push_back(*tabInPt[i_InSideUp]);
Polygon_SideInDown_VecLine [countDEMLine].push_back(*tabInPt[i_InSideDown]);
Polygon_SideOutUp_VecLine [countDEMLine].push_back(*tabInPt[i_OutSideUp]);
Polygon_SideOutDown_VecLine[countDEMLine].push_back(*tabInPt[i_OutSideDown]);
// DiapOTB only keep line in ]lo_margin, hi_margin[ but
// equal_range/lower_bound returns [lo, hi[
++lo_hi.first;
}
return lo_hi;
}
else
{
auto const lo_hi = equal_range_interval(
std::make_reverse_iterator(line_end),
std::make_reverse_iterator(line_start),
lo_margin, hi_margin, CmpVLVComponent<>(1));
auto lo = lo_hi.first;
auto hi = lo_hi.second;
auto lob = lo.base();
auto hib = hi.base();
// DiapOTB only keep line in ]lo_margin, hi_margin[ but
// equal_range/lower_bound returns [lo, hi[
// Also in the end after reversing, we need to return range that
// exclude hi, and keep lo
while (lob > hib && !is_in_window(line(lob-1)))
{
// std::cout << "L("<<line(lob)<<")";
--lob;
++lo; // 4 logs
}
//
while (hib < lob && !is_in_window(line(hib)))
{
// std::cout << "H("<<line(hib)<<")";
++hib;
--hi; // 4 logs
}
// std::cout
// << /*&(*lo)[0] <<*/ "(" << (*lo)[1] << '/' << line(lob) << ") -- "
// << /*&(*hi)[0] <<*/ "(" << (*hi)[1] << '/' << line(hib) << ")\t";
return std::make_pair(hib, lob);
}
}; // ------------------------------[ find_lo_hi_window
auto const line_start = point(0, dy);
auto const line_end = point(in_sizeX, dy);
auto const window = find_lo_hi_window(line_start, line_end);
auto const first = window.first;
auto const last = window.second;
// std::cout << "#y="<<(dy+in_startY) << "\t; start[0]="<<(*point(0, dy))[1]
// << "\t; end[" << (in_sizeX-1) << "]= " << (*point(in_sizeX-1, dy))[1]
// << "\t VS ["<<lo_margin<<", " << hi_margin << "]"
// << "\t -> " << (first - line_start) << " -- " << (last - line_start)
// << "\n";
has_match = first.data() != last.data();
if (has_match)
{ // a window has been found!
if (first.data() > line_start.data())
{ // assert before is outside
ASSERT(!is_in_window(line(first-1)), "")(in_data)(line(first-1))(dy)(lo_margin)(hi_margin)(line(line_start))(line(line_end-1));
}
if (last.data() < line_end.data())
{ // assert after is outside
ASSERT(!is_in_window(line(last)), "")(in_data)(line(last))(dy)(lo_margin)(hi_margin)(line(line_start))(line(line_end-1));
}
// after within is inside
ASSERT(is_in_window(line(first)), "")(in_data)(line(first))(dy)(lo_margin)(hi_margin)(line(line_start))(line(line_end-1));
if (last.data() > line_start.data())
ASSERT(is_in_window(line(last-1)), "")(in_data)(line(last-1))(line(last))(dy)(lo_margin)(hi_margin)(line(line_start))(line(line_end-1));
for (auto dx = (first-line_start); dx < (last-line_start) ; ++dx)
{
auto const pt_UL = *point(dx, dy);
auto const UL_line = pt_UL[1];
assert (is_in_window(UL_line));
register_pixel(dx, dy, pt_UL);
}
} // has_match
// cpt_out_window += (hi - lo);
}
else // DEM line is not monotonic in respect of progression of UL_line value
{
// For each column
for (unsigned dx = 0 ; dx < in_sizeX -1 ; ++dx)
{
auto const pt_UL = *point(dx, dy);
auto const UL_line = pt_UL[1];
if (UL_line > prev) ++nb_up;
else ++nb_dw;
prev = UL_line;
has_match = true;
// Protect Memory access and input access
static_assert(std::is_same<decltype(UL_line), decltype(lo_margin)>::value, "no convertion!!!");
if (is_in_window(UL_line)) // true 2% of the time...
{
register_pixel(dx, dy, pt_UL);
}
else
{
++cpt_out_window;
}
}
} // End colunms (input)
} // else: not monotonic
assert((nb_up == 0 || nb_dw == 0) && "Lines have a monotonic progression");
} // End colunms (input)
cpt_lines_with_no_match += has_match ? 0 : 1;
// if vector of for the current Line of DEM is not empty => store this vector into our vectors
......@@ -1033,6 +1206,10 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF
} // End lines (Main output)
std::cout << "IN VS OUT: " << cpt_in_window << " VS " << cpt_out_window
<< " / " << (cpt_out_window + cpt_in_window)
<< " ; lines w/o match: " << cpt_lines_with_no_match << " / " << nbDEMLines
<< std::endl;
}
} /*namespace otb*/
......
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