diff --git a/include/otbSARCartesianMeanFunctor.h b/include/otbSARCartesianMeanFunctor.h index 8e1eedc87316019189f1cc36bc5d6591c5434b7b..a6a3488c3f8bd804223ba075b5b21d57f47c5074 100644 --- a/include/otbSARCartesianMeanFunctor.h +++ b/include/otbSARCartesianMeanFunctor.h @@ -150,9 +150,10 @@ public: nbCol_into_outputRegion = this->GetNbColSAR() - firstCol_into_outputRegion + m_originC; } + // Because doubles are use internally, doubles are likelly to be mixed + // with TInputPixel==float double cE, cS; double dc; - int col1, col2; // Elevation angle (for shadow) double tgElevationMax = 0.; double tgElevation = 0.; @@ -221,8 +222,8 @@ public: dc = cS - cE; // Colunm into // TODO: he current maille - col1 = (int) (std::min (cE, cS))+1; // +1 ? - col2 = (int) (std::max (cE, cS)); + int const col1 = (int) (std::min (cE, cS))+1; // +1 ? + int const col2 = (int) (std::max (cE, cS)); // tgElevationMax per mailles (not sure) if (isFirstMaille) @@ -236,7 +237,8 @@ public: } // Loop on colunm - for (int k = std::max(col1, firstCol_into_outputRegion); k<= std::min(col2, nbCol_into_outputRegion+firstCol_into_outputRegion-1); k++) + auto const K = std::min(col2, nbCol_into_outputRegion+firstCol_into_outputRegion-1); + for (int k = std::max(col1, firstCol_into_outputRegion); k<= K; k++) { a = (cS - k) / dc; b = 1.0 - a; diff --git a/include/otbSARDEMPolygonsAnalysisImageFilter.h b/include/otbSARDEMPolygonsAnalysisImageFilter.h index f3e273cb39d687e3a416ed6bf7a1e19ff725e8e3..b31da3b5460a34f4c20c5a33613c8dcf6a518f31 100644 --- a/include/otbSARDEMPolygonsAnalysisImageFilter.h +++ b/include/otbSARDEMPolygonsAnalysisImageFilter.h @@ -88,12 +88,7 @@ public: using ImageInPointType = typename ImageInType::PointType; using ImageInValueType = typename ImageInPixelType::ValueType; #define OLD_ITER -#if 0 && defined(OLD_ITER) - using FastPixelType = ImageInPixelType; -#else - // using FastPixelType = otb::Span<ImageInValueType const>; using FastPixelType = ImageInValueType const*; -#endif /** Typedef to describe the image index, size types and spacing for inout image. */ using ImageInIndexType = typename ImageInType::IndexType; using ImageInIndexValueType = typename ImageInType::IndexValueType; @@ -141,7 +136,7 @@ public: // Define Point2DType and Point3DType using Point2DType = itk::Point<double,2>; - using Point3DType = itk::Point<double,3>; + // using Point3DType = itk::Point<double,3>; // Define RSTRansform For Localisation inverse (SAR to DEM) using RSTransformType2D = otb::GenericRSTransform<double,2,2>; @@ -220,10 +215,10 @@ protected: * Select for the current SAR Line (ind_LineSAR), the polygons that intersect the current SAR line. * This selection is always made with the SAR Line (even in ML geometry for the main output). */ - bool PolygonAndSarLineIntersection(const int ind_LineSAR, const double ligUL, const double colUL, - const double ligUR, const double colUR, - const double ligLR, const double colLR, const double ligLL, - const double colLL, + bool PolygonAndSarLineIntersection(const int ind_LineSAR, const ImageInValueType ligUL, const ImageInValueType colUL, + const ImageInValueType ligUR, const ImageInValueType colUR, + const ImageInValueType ligLR, const ImageInValueType colLR, const ImageInValueType ligLL, + const ImageInValueType colLL, int & i_InSideUp, int & i_InSideDown, int & i_OutSideUp, int & i_OutSideDown); diff --git a/include/otbSARDEMPolygonsAnalysisImageFilter.txx b/include/otbSARDEMPolygonsAnalysisImageFilter.txx index 21f76e9703ae012eb19c11c088f9f23d47f0d5e6..58b0adcc2402fd22d4ba336a1aaff833a3e98991 100644 --- a/include/otbSARDEMPolygonsAnalysisImageFilter.txx +++ b/include/otbSARDEMPolygonsAnalysisImageFilter.txx @@ -513,23 +513,25 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction> bool SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction > -::PolygonAndSarLineIntersection(const int ind_LineSAR, const double ligUL, const double colUL, - const double ligUR, const double colUR, - const double ligLR, const double colLR, const double ligLL, const double colLL, +::PolygonAndSarLineIntersection(const int ind_LineSAR, + const ImageInValueType ligUL, const ImageInValueType colUL, + const ImageInValueType ligUR, const ImageInValueType colUR, + const ImageInValueType ligLR, const ImageInValueType colLR, + const ImageInValueType ligLL, const ImageInValueType colLL, int & i_InSideUp, int & i_InSideDown, int & i_OutSideUp, int & i_OutSideDown) { - double const Lmax_ulur = std::max(ligUL, ligUR); - double const Lmin_ulur = std::min(ligUL, ligUR); + auto const Lmax_ulur = std::max(ligUL, ligUR); + auto const Lmin_ulur = std::min(ligUL, ligUR); - double const Lmax_urlr = std::max(ligLR, ligUR); - double const Lmin_urlr = std::min(ligLR, ligUR); + auto const Lmax_urlr = std::max(ligLR, ligUR); + auto const Lmin_urlr = std::min(ligLR, ligUR); - double const Lmax_lllr = std::max(ligLL, ligLR); - double const Lmin_lllr = std::min(ligLL, ligLR); + auto const Lmax_lllr = std::max(ligLL, ligLR); + auto const Lmin_lllr = std::min(ligLL, ligLR); - double const Lmax_ulll = std::max(ligUL, ligLL); - double const Lmin_ulll = std::min(ligUL, ligLL); + auto const Lmax_ulll = std::max(ligUL, ligLL); + auto const Lmin_ulll = std::min(ligUL, ligLL); // Selection for ind_LineSAR intersection auto const intersect_ulur = Lmax_ulur >= ind_LineSAR && Lmin_ulur <= ind_LineSAR; @@ -547,8 +549,8 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF } int countSides = 0; - double CLZY_Up [4][2]; - double CLZY_Down [4][2]; + ImageInValueType CLZY_Up [4][2]; + ImageInValueType CLZY_Down [4][2]; int i_Up [4]; int i_Down [4]; @@ -661,12 +663,12 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF else // if (countSides == 3) { // Select extermes colunms among the three CLZY points to cover the biggest part of SAR Image - double const Cmax_0 = std::max(CLZY_Up[0][0], CLZY_Down[0][0]); - double const Cmin_0 = std::min(CLZY_Up[0][0], CLZY_Down[0][0]); - double const Cmax_1 = std::max(CLZY_Up[1][0], CLZY_Down[1][0]); - double const Cmin_1 = std::min(CLZY_Up[1][0], CLZY_Down[1][0]); - double const Cmax_2 = std::max(CLZY_Up[2][0], CLZY_Down[2][0]); - double const Cmin_2 = std::min(CLZY_Up[2][0], CLZY_Down[2][0]); + auto const Cmax_0 = std::max(CLZY_Up[0][0], CLZY_Down[0][0]); + auto const Cmin_0 = std::min(CLZY_Up[0][0], CLZY_Down[0][0]); + auto const Cmax_1 = std::max(CLZY_Up[1][0], CLZY_Down[1][0]); + auto const Cmin_1 = std::min(CLZY_Up[1][0], CLZY_Down[1][0]); + auto const Cmax_2 = std::max(CLZY_Up[2][0], CLZY_Down[2][0]); + auto const Cmin_2 = std::min(CLZY_Up[2][0], CLZY_Down[2][0]); // Max colunm if (Cmax_0 > Cmax_1 && Cmax_0 > Cmax_2) @@ -793,18 +795,25 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF int const nbLinesIntoDEM_withSelectedPolygon = CLZY_InSideUp_Vec.size(); + // Memorize the "external" member variables. + // * -1 index descreases, +1 index increases + // * Ideally it should be more efficient to have 4 pairs of loop + // It should be easier to acheive w/ 2 lambdas, but PolygonsScan() isn't a + // bottleneck and doesn't deserve obfuscation. Let's bet on branch + // prediction. + auto const shallIncreaseDirectionL = m_DEMdirL == 1; + auto const shallIncreaseDirectionC = m_DEMdirC == 1; + int ind_LDEM, ind_CDEM; // Loop on polygons (for lines into DEM with al least one selected polygon) // the loop direction depends on m_DEMdirL for (int i = 0 ; i < nbLinesIntoDEM_withSelectedPolygon; i++) { - // If 1, the lines increase - if (m_DEMdirL == 1) + if (shallIncreaseDirectionL) { ind_LDEM = i; } - // If -1, the lines decrease else { ind_LDEM = nbLinesIntoDEM_withSelectedPolygon - (i+1); @@ -814,12 +823,10 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF int const nbColIntoDEM_withSelectedPolygon = CLZY_InSideUp_Vec[ind_LDEM]->size(); for (int j = 0; j < nbColIntoDEM_withSelectedPolygon; j++) { - // If 1, the columns increase - if (m_DEMdirC == 1) + if (shallIncreaseDirectionC) { ind_CDEM = j; } - // If -1, the colunms decrease else { ind_CDEM = nbColIntoDEM_withSelectedPolygon - (j+1); @@ -829,7 +836,6 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF // estimate the wanted value (for example amplitude) auto as_vlv = [nbcompo](FastPixelType p) { return ImageInPixelType(p, nbcompo, false); - // return ImageInPixelType(p.data(), p.size(), false); }; m_FunctionOnPolygon->estimate(ind_LineSAR, as_vlv((*CLZY_InSideUp_Vec[ind_LDEM])[ind_CDEM]), @@ -926,7 +932,7 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF auto const margin = m_Margin; - auto const nodata = m_NoData; + ImageInValueType const nodata = m_NoData; // Pixel of otb::VectorImage (input) : 4 Components : C, L, Z and Y @@ -938,11 +944,11 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF tabInIt[3] = &InIt_LL; #else // Make sure they are all proxies => move assign will continue to move - FastPixelType tabIn[4] ; - auto& UL = tabIn[0]; - auto& UR = tabIn[1]; - auto& LR = tabIn[2]; - auto& LL = tabIn[3]; + // FastPixelType tabIn[4] ; + // auto& UL = tabIn[0]; + // auto& UR = tabIn[1]; + // auto& LR = tabIn[2]; + // auto& LL = tabIn[3]; #endif // For each line of output : construct polygon after polygon, determine if polygon intersects the current @@ -969,6 +975,8 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF for (int ind_Line_Into_SARGeo = ind_Line_start; ind_Line_Into_SARGeo < ind_Line_end; ind_Line_Into_SARGeo++) { + ImageInValueType const lo_margin = ind_Line_Into_SARGeo - margin; + ImageInValueType const hi_margin = ind_Line_Into_SARGeo + margin; ///////////////////////////// Polygon Construction and Check intersection ///////////////////////// #if defined(OLD_ITER) InIt_UL.GoToBegin(); @@ -1013,8 +1021,10 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF InIt_Lx.GoToBeginOfLine(); // Put each Iterator to this start colunm // +1 colunm for UR and LR - UL = as_span(InIt_Ux.Get()); // expect the line to not be empty... - LL = as_span(InIt_Lx.Get()); + // UL = as_ptr(InIt_Ux.Get()); // expect the line to not be empty... + // LL = as_ptr(InIt_Lx.Get()); + auto UL = as_ptr(InIt_Ux.Get()); // expect the line to not be empty... + auto LL = as_ptr(InIt_Lx.Get()); ++InIt_Ux; ++InIt_Lx; #endif @@ -1023,7 +1033,7 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF #if defined(OLD_ITER) while (!InIt_UR.IsAtEndOfLine() && !InIt_LR.IsAtEndOfLine()) #else - while (!InIt_Ux.IsAtEndOfLine() && !InIt_Lx.IsAtEndOfLine()) + while (!InIt_Ux.IsAtEndOfLine()) #endif { @@ -1036,29 +1046,32 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF assert(UR.IsAProxy()); assert(LL.IsAProxy()); assert(LR.IsAProxy()); - double const UL_line = UL[1]; + auto const UL_line = UL[1]; #else - UR = as_span(InIt_Ux.Get()); - LR = as_span(InIt_Lx.Get()); - double const UL_line = UL[1]; + assert(!InIt_Lx.IsAtEndOfLine() && "Iterators are defined on the same region and are on the same column"); + // UR = as_ptr(InIt_Ux.Get()); + // LR = as_ptr(InIt_Lx.Get()); + auto const UR = as_ptr(InIt_Ux.Get()); + auto const LR = as_ptr(InIt_Lx.Get()); + auto const UL_line = UL[1]; #endif // Protect Memory access and input access - if (UL_line > ind_Line_Into_SARGeo - margin && UL_line < ind_Line_Into_SARGeo + margin) + if (UL_line > lo_margin && UL_line < hi_margin) { - double const UL_col = UL[0]; - double const UR_col = UR[0]; - double const LR_col = LR[0]; - double const LL_col = LL[0]; + auto const UL_col = UL[0]; + auto const UR_col = UR[0]; + auto const LR_col = LR[0]; + auto const LL_col = LL[0]; // Check if No Data if (!(UL_col == nodata || UR_col == nodata || LR_col == nodata || LL_col == nodata)) { // Determine if intersection - double const UR_line = UR[1]; - double const LR_line = LR[1]; - double const LL_line = LL[1]; + auto const UR_line = UR[1]; + auto const LR_line = LR[1]; + auto const LL_line = LL[1]; // With lig and col of UL, UR, LR and LL select iterators for InSide and Outside of // each polygon (if intersection) @@ -1084,17 +1097,19 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF Polygon_SideOutUp_VecLine [countDEMLine].push_back(as_ptr(tabInIt[i_OutSideUp ]->Get())); Polygon_SideOutDown_VecLine[countDEMLine].push_back(as_ptr(tabInIt[i_OutSideDown]->Get())); #else - auto& inSideUp = tabIn[i_InSideUp ]; - auto& inSideDown = tabIn[i_InSideDown ]; - auto& outSideUp = tabIn[i_OutSideUp ]; - auto& outSideDown = tabIn[i_OutSideDown]; + // TODO: make tab at the last moment as intersections are a + // rare case + FastPixelType tabIn [4] = {UL, UR, LR, LL}; + auto const inSideUp = tabIn[i_InSideUp ]; + auto const inSideDown = tabIn[i_InSideDown ]; + auto const outSideUp = tabIn[i_OutSideUp ]; + auto const outSideDown = tabIn[i_OutSideDown]; Polygon_SideInUp_VecLine [countDEMLine].push_back(inSideUp); Polygon_SideInDown_VecLine [countDEMLine].push_back(inSideDown); Polygon_SideOutUp_VecLine [countDEMLine].push_back(outSideUp); Polygon_SideOutDown_VecLine[countDEMLine].push_back(outSideDown); #endif - } } } @@ -1106,9 +1121,9 @@ SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TF ++InIt_LR; ++InIt_LL; #else - LL = std::move(LR); + LL = LR; ++InIt_Lx; - UL = std::move(UR); + UL = UR; ++InIt_Ux; #endif