From c74aefe94fcf497ee72e806284dfa484efb0cb31 Mon Sep 17 00:00:00 2001 From: StaronC Date: Tue, 30 Jan 2024 17:03:19 +0100 Subject: [PATCH] Included radiometry factor in new form equation, added possibility to use multi-scale approach, initialise values of unknowns with values obtained in previous iteration along with hard and soft freezing unknowns. Added option to display last translations and radiometry unknown values. Added command and equation to specifically compute radiometric deformations between images and another for translations. --- MMVII/include/MMVII_DeclareAllCmd.h | 6 + MMVII/include/MMVII_PhgrDist.h | 10 + MMVII/include/MMVII_TplSymbTriangle.h | 118 +++ MMVII/src/Appli/cSpecMMVII_Appli.cpp | 6 + .../ApplyDelaunayOnRandomGeneratedPoints.cpp | 251 +++++ MMVII/src/MeshDisplacement/SimulDispl.cpp | 313 +++--- .../MeshDisplacement/TriangleDeformation.cpp | 921 ++++++++++++++++++ .../MeshDisplacement/TriangleDeformation.h | 180 ++++ .../TriangleDeformationRad.cpp | 405 ++++++++ .../MeshDisplacement/TriangleDeformationRad.h | 92 ++ .../TriangleDeformationRadiometry.cpp | 457 +++++++++ .../TriangleDeformationRadiometry.h | 109 +++ .../TriangleDeformationTrRad.cpp | 583 +++++++++++ .../TriangleDeformationTrRad.h | 113 +++ .../TriangleDeformationTranslation.cpp | 463 +++++++++ .../TriangleDeformationTranslation.h | 111 +++ .../src/SymbDerGen/Formulas_TrianglesDeform.h | 324 ++++++ MMVII/src/SymbDerGen/GenerateCodes.cpp | 32 + 18 files changed, 4320 insertions(+), 174 deletions(-) create mode 100644 MMVII/include/MMVII_TplSymbTriangle.h create mode 100644 MMVII/src/MeshDisplacement/ApplyDelaunayOnRandomGeneratedPoints.cpp create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformation.cpp create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformation.h create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformationRad.cpp create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformationRad.h create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformationRadiometry.cpp create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformationRadiometry.h create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformationTrRad.cpp create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformationTrRad.h create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformationTranslation.cpp create mode 100644 MMVII/src/MeshDisplacement/TriangleDeformationTranslation.h create mode 100644 MMVII/src/SymbDerGen/Formulas_TrianglesDeform.h diff --git a/MMVII/include/MMVII_DeclareAllCmd.h b/MMVII/include/MMVII_DeclareAllCmd.h index 462d8f05dc..2048c1d161 100755 --- a/MMVII/include/MMVII_DeclareAllCmd.h +++ b/MMVII/include/MMVII_DeclareAllCmd.h @@ -82,6 +82,12 @@ extern cSpecMMVII_Appli TheSpec_ConvertV1V2_GCPIM; extern cSpecMMVII_Appli TheSpec_SpecSerial; extern cSpecMMVII_Appli TheSpec_CGPReport; extern cSpecMMVII_Appli TheSpec_TiePReport; +extern cSpecMMVII_Appli TheSpec_RandomGeneratedDelaunay; +extern cSpecMMVII_Appli TheSpec_ComputeTriangleDeformation; +extern cSpecMMVII_Appli TheSpec_ComputeTriangleDeformationTrRad; +extern cSpecMMVII_Appli TheSpec_ComputeTriangleDeformationTranslation; +extern cSpecMMVII_Appli TheSpec_ComputeTriangleDeformationRadiometry; +extern cSpecMMVII_Appli TheSpec_ComputeTriangleDeformationRad; extern cSpecMMVII_Appli TheSpec_PoseCmpReport; extern cSpecMMVII_Appli TheSpec_BlockCamInit; // RIGIDBLOC extern cSpecMMVII_Appli TheSpec_ClinoInit; diff --git a/MMVII/include/MMVII_PhgrDist.h b/MMVII/include/MMVII_PhgrDist.h index 53ab8925f8..698fb2d8db 100755 --- a/MMVII/include/MMVII_PhgrDist.h +++ b/MMVII/include/MMVII_PhgrDist.h @@ -198,6 +198,16 @@ NS_SymbolicDerivative::cCalculator * EqTopoDZ(bool WithDerive,int aSzBuf NS_SymbolicDerivative::cCalculator * EqDeformImHomotethy(bool WithDerive,int aSzBuf); /// Variant of "EqDeformImHomotethy", case where we use linear approximation NS_SymbolicDerivative::cCalculator * EqDeformImLinearGradHomotethy(bool WithDerive,int aSzBuf); +/// Equation used to optimize translation and radiometry transformation between two images. +NS_SymbolicDerivative::cCalculator *EqDeformTri(bool WithDerive, int aSzBuf); +/// Equation used to optimize translation and radiometry transformation between two images with alternative equation +NS_SymbolicDerivative::cCalculator *EqDeformTriTrRad(bool WithDerive, int aSzBuf); +/// Equation used to optimize translation transformation between two images. +NS_SymbolicDerivative::cCalculator *EqDeformTriTranslation(bool WithDerive, int aSzBuf); +/// Equation used to optimize radiometric transformation between two images. +NS_SymbolicDerivative::cCalculator *EqDeformTriRadiometry(bool WithDerive, int aSzBuf); +/// Equation used to optimize radiometric transformation between two images with alternate equation +NS_SymbolicDerivative::cCalculator *EqDeformTriRad(bool WithDerive, int aSzBuf); // ............. Covariance propagation ............. diff --git a/MMVII/include/MMVII_TplSymbTriangle.h b/MMVII/include/MMVII_TplSymbTriangle.h new file mode 100644 index 0000000000..d13757a48b --- /dev/null +++ b/MMVII/include/MMVII_TplSymbTriangle.h @@ -0,0 +1,118 @@ +#ifndef _MMVII_TplSymbTriangle_ +#define _MMVII_TplSymbTriangle_ + +#include "SymbDer/SymbolicDerivatives.h" +#include "SymbDer/SymbDer_MACRO.h" + +#include "MMVII_Ptxd.h" +#include "MMVII_Geom2D.h" + +#include "../src/MeshDisplacement/TriangleDeformation.h" + +/** + \file MMVII_TplSymbTriangle.h + \brief Contains helpers for triangle as formula +**/ + +using namespace NS_SymbolicDerivative; + +namespace MMVII +{ + + /** + Compute the formula for composition of a function and an image. + + Let I be an image, considered as function R^2 -> R using an interpolation model, this computes + the formula corresponding to the function + + x,y -> I(Fx(x,y),Fy(x,y)) + + We use here the bilinear interpolation. Let X0,Y0 be integers such that: + + X0 <= Fx(x,y) < X0+1 + Y0 <= Fy(x,y) < Y0+1 + Let X1=X0+1 , Y1=Y0+1, X=Fx(x,y),Y=Fy(x,y) + + Let I00= I[X0,Y0] , I10 = I[X1,Y0] ... the bilinear formula is : + + Bil(I,X,Y) = + I00 * (X1-X) * (Y1-Y) + + I10 * (X-X0) * (Y1-Y) + + I01 * (X1-X) * (Y-Y0) + + I11 * (X-X0) * (Y-Y0) + + In this function the value X0, Y0, I00, I10, I01, I11 are communicated as element + of the observation vector aVObs, starting from aKObs0. + **/ + + template + TypeFunc FormalBilinTri_Formula( + const std::vector &aVObs, + int aKObs0, + const TypeFunc &FX, + const TypeFunc &FY) + { + TypeFunc aX0(aVObs.at(aKObs0)); + TypeFunc aY0(aVObs.at(aKObs0 + 1)); + TypeFunc aCst1 = CreateCste(1.0, aX0); // create a symbolic formula for constant 1 + + TypeFunc aWX1 = FX - aX0; // weight for I10 and I11 + TypeFunc aWX0 = aCst1 - aWX1; // weight for I00 and I01 + TypeFunc aWY1 = FY - aY0; // weight for I10 and I11 + TypeFunc aWY0 = aCst1 - aWY1; // weight for I00 and I10 + + return aWX0 * aWY0 * aVObs.at(aKObs0 + 2) // I00 + + aWX1 * aWY0 * aVObs.at(aKObs0 + 3) // I10 + + aWX0 * aWY1 * aVObs.at(aKObs0 + 4) // I01 + + aWX1 * aWY1 * aVObs.at(aKObs0 + 5); // I11 + } + + // standard name for observation + std::vector FormalBilinIm2D_NameObs(const std::string &aPrefix); + + template + void FormalInterpBarycenter_SetObs( + std::vector &aVObs, // vector of observation to fill + const int aK0, // first index where fill the vector + const cPtInsideTriangles &aPixInsideTriangle) + { + // push integer coordinate of point + SetOrPush(aVObs, aK0, Type(aPixInsideTriangle.GetCartesianCoordinates().x())); + SetOrPush(aVObs, aK0 + 1, Type(aPixInsideTriangle.GetCartesianCoordinates().y())); + SetOrPush(aVObs, aK0 + 2, Type(aPixInsideTriangle.GetBarycenterCoordinates().x())); + SetOrPush(aVObs, aK0 + 3, Type(aPixInsideTriangle.GetBarycenterCoordinates().y())); + SetOrPush(aVObs, aK0 + 4, Type(aPixInsideTriangle.GetBarycenterCoordinates().z())); + SetOrPush(aVObs, aK0 + 5, Type(aPixInsideTriangle.GetPixelValue())); + } + + /* + This is the "companion" function of FormalBilinIm2D_Formula, it fills + the vector aVObs with X0,Y0,I00, that will be used in FormalBilinIm2D_Formula. + */ + + template + void FormalBilinTri_SetObs( + std::vector &aVObs, // vector of observation to fill + const int aK0, // first index where fill the vector + const cPt2dr aPtIm, // point in image + const cDataIm2D &aDIm // image + ) + { + // compute coordinate of left-high corner of including pixel + cPt2di aP0 = Pt_round_down(aPtIm); + + // push integer coordinate of point + SetOrPush(aVObs, aK0, Type(aP0.x())); + SetOrPush(aVObs, aK0 + 1, Type(aP0.y())); + + // push values of image at its 4 corners + SetOrPush(aVObs, aK0 + 2, (Type)aDIm.GetV(aP0)); + SetOrPush(aVObs, aK0 + 3, (Type)aDIm.GetV(aP0 + cPt2di(1, 0))); + SetOrPush(aVObs, aK0 + 4, (Type)aDIm.GetV(aP0 + cPt2di(0, 1))); + SetOrPush(aVObs, aK0 + 5, (Type)aDIm.GetV(aP0 + cPt2di(1, 1))); + } + + constexpr int TriangleDisplacement_NbObs = 6; +}; + +#endif // _MMVII_TplSymbTriangle_ \ No newline at end of file diff --git a/MMVII/src/Appli/cSpecMMVII_Appli.cpp b/MMVII/src/Appli/cSpecMMVII_Appli.cpp index 8c014ea618..a1d1cc607d 100755 --- a/MMVII/src/Appli/cSpecMMVII_Appli.cpp +++ b/MMVII/src/Appli/cSpecMMVII_Appli.cpp @@ -239,6 +239,12 @@ std::vector & cSpecMMVII_Appli::InternVecAll() TheVecAll.push_back(&TheSpec_TestProj); TheVecAll.push_back(&TheSpec_ChSysCo); TheVecAll.push_back(&TheSpec_CreateCalib); + TheVecAll.push_back(&TheSpec_RandomGeneratedDelaunay); + TheVecAll.push_back(&TheSpec_ComputeTriangleDeformation); + TheVecAll.push_back(&TheSpec_ComputeTriangleDeformationTrRad); + TheVecAll.push_back(&TheSpec_ComputeTriangleDeformationTranslation); + TheVecAll.push_back(&TheSpec_ComputeTriangleDeformationRadiometry); + TheVecAll.push_back(&TheSpec_ComputeTriangleDeformationRad); TheVecAll.push_back(&TheSpec_ImportTiePMul); TheVecAll.push_back(&TheSpec_ImportMesImGCP); TheVecAll.push_back(&TheSpecImportExtSens); diff --git a/MMVII/src/MeshDisplacement/ApplyDelaunayOnRandomGeneratedPoints.cpp b/MMVII/src/MeshDisplacement/ApplyDelaunayOnRandomGeneratedPoints.cpp new file mode 100644 index 0000000000..63077a9308 --- /dev/null +++ b/MMVII/src/MeshDisplacement/ApplyDelaunayOnRandomGeneratedPoints.cpp @@ -0,0 +1,251 @@ +#include "cMMVII_Appli.h" +#include "MMVII_Geom2D.h" + + +/** + \file ApplyDelaunayOnRandomGeneratedPoints.cpp + + \brief file for generating random points distributed uniformely + and applying 2D Delaunay triangulation. + +**/ + +namespace MMVII +{ + + /* ==================================================== */ + /* */ + /* cAppli_RandomGeneratedDelaunay */ + /* */ + /* ==================================================== */ + + class cAppli_RandomGeneratedDelaunay : public cMMVII_Appli + { + public: + typedef cIm2D tIm; + typedef cDataIm2D tDIm; + typedef cTriangulation2D tTriangule2dr; + + cAppli_RandomGeneratedDelaunay(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec); + + int Exe() override; + cCollecSpecArg2007 &ArgObl(cCollecSpecArg2007 &anArgObl) override; + cCollecSpecArg2007 &ArgOpt(cCollecSpecArg2007 &anArgOpt) override; + + void ApplyAndSaveDelaunayTriangulationOnPoints(); + void GeneratePointsForDelaunay(); + void ConstructUniformRandomVector(); + void ComputeShiftedTrianglesAndPoints(const cTriangle & aTri, std::vector & ShiftedTriangleCoordinates, + const double aUniformRandomRedChannel, const double aUniformRandomGreenChannel, + const double aUniformRandomBlueChannel); + + private: + // == Mandatory args ==== + + std::string mNameInputImage; + std::string mNamePlyFile; + int mNumberPointsToGenerate; // number of generated points + + // == Optionnal args ==== + + int mRandomUniformLawUpperBoundXAxis; // Uniform law generate numbers from [0, mRandomUniformLawUpperBound [ for x-axis + int mRandomUniformLawUpperBoundYAxis; // Uniform law generate numbers from [0, mRandomUniformLawUpperBound [ for y-axis + bool mPlyFileisBinary; + bool mShiftTriangles; + std::string mNameModifedTrianglePlyFile; + + // == Internal variables ==== + + tIm mImIn; ///< memory representation of the image + tDIm *mDImIn; ///< memory representation of the image + cPt2di mSz; // Size of images + std::vector mVectorPts; // A vector containing a set of points + }; + + cAppli_RandomGeneratedDelaunay::cAppli_RandomGeneratedDelaunay(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) : cMMVII_Appli(aVArgs, aSpec), + mRandomUniformLawUpperBoundXAxis(1), + mRandomUniformLawUpperBoundYAxis(1), + mPlyFileisBinary(false), + mShiftTriangles(true), + mImIn(cPt2di(1, 1)), + mDImIn(nullptr) + { + } + + cCollecSpecArg2007 &cAppli_RandomGeneratedDelaunay::ArgObl(cCollecSpecArg2007 &anArgObl) + { + return anArgObl + << Arg2007(mNameInputImage, "Name of input image file.", {{eTA2007::FileImage}, {eTA2007::FileDirProj}}) + << Arg2007(mNamePlyFile, "Name of main triangulation file to save in .ply format.", {{eTA2007::FileCloud}}) + << Arg2007(mNumberPointsToGenerate, "Number of points you want to generate for triangulation."); + } + + cCollecSpecArg2007 &cAppli_RandomGeneratedDelaunay::ArgOpt(cCollecSpecArg2007 &anArgOpt) + { + return anArgOpt + << AOpt2007(mRandomUniformLawUpperBoundXAxis, "RandomUniformLawXAxis", "Maximum value that the uniform law can draw from for x-axis.", {eTA2007::HDV}) + << AOpt2007(mRandomUniformLawUpperBoundYAxis, "RandomUniformLawYAxis", "Maximum value that the uniform law can draw from for y-axis.", {eTA2007::HDV}) + << AOpt2007(mShiftTriangles, "ShiftTriangles", "Whether to shift points of triangles after application of Delaunay triangulation.", {eTA2007::HDV}) + << AOpt2007(mPlyFileisBinary, "PlyFileIsBinary", "Whether to save the .ply file binarised or not.", {eTA2007::HDV}) + << AOpt2007(mNameModifedTrianglePlyFile, "NamePlyFileShiftedTriangles", "Name of .ply file for shifted triangles.", {eTA2007::FileCloud}); + } + + //========================================================= + + void cAppli_RandomGeneratedDelaunay::ComputeShiftedTrianglesAndPoints(const cTriangle & aTri, + std::vector & ShiftedTriangleCoordinates, + const double aUniformRandomRedChannel, + const double aUniformRandomGreenChannel, + const double aUniformRandomBlueChannel) + { + const cTriangle2DCompiled aCompTri(aTri); + + // Compute shifts depending on point of triangle + const cPt2dr PercentDiffA = 0.01 * (aTri.Pt(0) - aTri.Pt(2)); + const cPt2dr PercentDiffB = 0.015 * (aTri.Pt(1) - aTri.Pt(0)); + const cPt2dr PercentDiffC = 0.02 * (aTri.Pt(2) - aTri.Pt(1)); + + ShiftedTriangleCoordinates.push_back(aTri.Pt(0) + PercentDiffA); + ShiftedTriangleCoordinates.push_back(aTri.Pt(1) + PercentDiffB); + ShiftedTriangleCoordinates.push_back(aTri.Pt(2) + PercentDiffC); + + // Get pixels inside each triangle and shift them + std::vector aVectorToFillwithInsidePixels; + aCompTri.PixelsInside(aVectorToFillwithInsidePixels); + for (long unsigned int FilledPixel=0; FilledPixel < aVectorToFillwithInsidePixels.size(); FilledPixel++) + { + if (FilledPixel % 10 == 0) + { + const cPt2dr aFilledPoint(aVectorToFillwithInsidePixels[FilledPixel].x(), aVectorToFillwithInsidePixels[FilledPixel].y()); + const cPt3dr barycenter_coordinates = aCompTri.CoordBarry(aFilledPoint); + const cPt2dr ShiftedInsidePixels = cPt2dr(aFilledPoint.x() + barycenter_coordinates.x() * PercentDiffA.x() + + barycenter_coordinates.y() * PercentDiffB.x() + barycenter_coordinates.z() * PercentDiffC.x(), + aFilledPoint.y() + barycenter_coordinates.x() * PercentDiffA.y() + + barycenter_coordinates.y() * PercentDiffB.y() + barycenter_coordinates.z() * PercentDiffC.y()); + + StdOut() << aFilledPoint.x() << " " << aFilledPoint.y() << " " << aUniformRandomRedChannel + << " " << aUniformRandomGreenChannel << " " << aUniformRandomBlueChannel << std::endl; + StdOut() << ShiftedInsidePixels.x() << " " << ShiftedInsidePixels.y() << " " + << aUniformRandomRedChannel << " " << aUniformRandomGreenChannel << " " << aUniformRandomBlueChannel << std::endl; + } + } + } + + void cAppli_RandomGeneratedDelaunay::ApplyAndSaveDelaunayTriangulationOnPoints() + { + tTriangule2dr aDelTri(mVectorPts); + + aDelTri.MakeDelaunay(); + + std::vector ShiftedTriangleCoordinates; + + // Loop over all triangle + for (size_t aKt = 0; aKt < aDelTri.NbFace(); aKt++) + { + const cTriangle aTri = aDelTri.KthTri(aKt); + + if (mShiftTriangles) + { + // for colouring points in representation + const double aUniformRandomRedChannel = RandUnif_N(256); + const double aUniformRandomGreenChannel = RandUnif_N(256); + const double aUniformRandomBlueChannel = RandUnif_N(256); + ComputeShiftedTrianglesAndPoints(aTri, ShiftedTriangleCoordinates, aUniformRandomRedChannel, + aUniformRandomGreenChannel, aUniformRandomBlueChannel); + } + } + if (mShiftTriangles) + { + tTriangule2dr aModifiedDelTri(ShiftedTriangleCoordinates); + + aModifiedDelTri.MakeDelaunay(); + // Save files to .ply format + aModifiedDelTri.WriteFile(mNameModifedTrianglePlyFile, mPlyFileisBinary); + } + + // Save files to .ply format + aDelTri.WriteFile(mNamePlyFile, mPlyFileisBinary); + } + + void cAppli_RandomGeneratedDelaunay::ConstructUniformRandomVector() + { + // Use current time as seed for random generator + srand(time(0)); + // Generate coordinates from drawing lines and columns of coordinates from a uniform distribution + for (int aNbPt = 0; aNbPt < mNumberPointsToGenerate; aNbPt++) + { + const double aUniformRandomXAxis = RandUnif_N(mRandomUniformLawUpperBoundXAxis); + const double aUniformRandomYAxis = RandUnif_N(mRandomUniformLawUpperBoundYAxis); + const cPt2dr aUniformRandomPt(aUniformRandomXAxis, aUniformRandomYAxis); // cPt2dr format + mVectorPts.push_back(aUniformRandomPt); + } + } + + void cAppli_RandomGeneratedDelaunay::GeneratePointsForDelaunay() + { + // If user hasn't defined another value than the default value, it is changed + if (mRandomUniformLawUpperBoundXAxis == 1 && mRandomUniformLawUpperBoundXAxis == 1) + { + // Maximum value of coordinates are drawn from [0, NumberOfImageLines[ + mRandomUniformLawUpperBoundXAxis = mSz.x(); + mRandomUniformLawUpperBoundYAxis = mSz.y(); + } + else + { + if (mRandomUniformLawUpperBoundXAxis != 1 && mRandomUniformLawUpperBoundYAxis == 1) + mRandomUniformLawUpperBoundYAxis = mSz.y(); + else + { + if (mRandomUniformLawUpperBoundXAxis == 1 && mRandomUniformLawUpperBoundYAxis != 1) + mRandomUniformLawUpperBoundXAxis = mSz.x(); + } + } + + ConstructUniformRandomVector(); + + ApplyAndSaveDelaunayTriangulationOnPoints(); // Apply Delaunay triangulation on generated points. + } + + //---------------------------------------- + + int cAppli_RandomGeneratedDelaunay::Exe() + { + /* + MMVII RandomGeneratedDelaunay pair18_im1_720.png OriginalTriangles.ply 20 50 NamePlyFileShiftedTriangles=ShiftedTriangles.ply > OriginalAndShiftedPoints.xyz + awk NR%2==1 < OriginalAndShiftedPoints.xyz > OriginalPoints.xyz + awk NR%2==0 < OriginalAndShiftedPoints.xyz > ShiftedPoints.xyz + */ + + mImIn = tIm::FromFile(mNameInputImage); + + mDImIn = &mImIn.DIm(); + mSz = mDImIn->Sz(); + + GeneratePointsForDelaunay(); + + return EXIT_SUCCESS; + } + + /* ================================== */ + /* */ + /* MMVII */ + /* */ + /* ================================== */ + + tMMVII_UnikPApli Alloc_RandomGeneratedDelaunay(const std::vector &aVArgs, const cSpecMMVII_Appli &aSpec) + { + return tMMVII_UnikPApli(new cAppli_RandomGeneratedDelaunay(aVArgs, aSpec)); + } + + cSpecMMVII_Appli TheSpec_RandomGeneratedDelaunay( + "RandomGeneratedDelaunay", + Alloc_RandomGeneratedDelaunay, + "Generate random points and apply Delaunay triangulation", + {eApF::ImProc}, // category + {eApDT::Image}, // input + {eApDT::Ply}, // output + __FILE__); + +}; // MMVII diff --git a/MMVII/src/MeshDisplacement/SimulDispl.cpp b/MMVII/src/MeshDisplacement/SimulDispl.cpp index 26ab111481..f63894c74f 100644 --- a/MMVII/src/MeshDisplacement/SimulDispl.cpp +++ b/MMVII/src/MeshDisplacement/SimulDispl.cpp @@ -3,211 +3,176 @@ #include "MMVII_Linear2DFiltering.h" #include "MMVII_Tpl_Images.h" - /** \file SimulDispl.cpp \brief file for generating simulation of smooth displacement - - */ + **/ namespace MMVII { + /* ======================================== */ + /* */ + /* cAppli_SimulDispl */ + /* */ + /* ======================================== */ -/* ==================================================== */ -/* */ -/* cAppli_CalibratedSpaceResection */ -/* */ -/* ==================================================== */ - -class cAppli_SimulDispl : public cMMVII_Appli -{ - public : - typedef cIm2D tIm; - typedef cDataIm2D tDIm; - typedef cIm2D tImDispl; - typedef cDataIm2D tDImDispl; + class cAppli_SimulDispl : public cMMVII_Appli + { + public: + typedef cIm2D tIm; + typedef cDataIm2D tDIm; + typedef cIm2D tImDispl; + typedef cDataIm2D tDImDispl; - cAppli_SimulDispl(const std::vector & aVArgs,const cSpecMMVII_Appli &); + cAppli_SimulDispl(const std::vector &aVArgs, + const cSpecMMVII_Appli &); int Exe() override; - cCollecSpecArg2007 & ArgObl(cCollecSpecArg2007 & anArgObl) override; - cCollecSpecArg2007 & ArgOpt(cCollecSpecArg2007 & anArgOpt) override; - - private : - - tImDispl GenerateSmoothRandDispl(); - - // == Mandatory args ==== - std::string mNameImage; ///< name of the input image to deform - - // == Optionnal args ==== - - tREAL8 mAmplDef; - bool mWithDisc; - - // == Internal variables ==== - tIm mImIn; ///< memory representation of the image - tDIm * mDImIn; ///< memory representation of the image - cPt2di mSz; - tIm mImOut; ///< memory representation of the image - tDIm * mDImOut; ///< memory representation of the image - - -}; - -cAppli_SimulDispl::cAppli_SimulDispl -( - const std::vector & aVArgs, - const cSpecMMVII_Appli & aSpec -) : - cMMVII_Appli (aVArgs,aSpec), - mAmplDef (2.0), - mWithDisc (true), - mImIn (cPt2di(1,1)), - mDImIn (nullptr), - mImOut (cPt2di(1,1)), - mDImOut (nullptr) -{ -} - - - -cCollecSpecArg2007 & cAppli_SimulDispl::ArgObl(cCollecSpecArg2007 & anArgObl) -{ - return anArgObl - << Arg2007(mNameImage,"Name of image to deform",{{eTA2007::FileImage},{eTA2007::FileDirProj}}) - ; -} - -cCollecSpecArg2007 & cAppli_SimulDispl::ArgOpt(cCollecSpecArg2007 & anArgOpt) -{ - - return anArgOpt - << AOpt2007(mAmplDef,"Ampl","Amplitude of deformation",{eTA2007::HDV}) - << AOpt2007(mWithDisc,"WithDisc","Do we add disconinuities",{eTA2007::HDV}) - /* - << AOpt2007(mNbIterBundle,"NbIterBund","Number of bundle iteration, after ransac init",{eTA2007::HDV}) - << AOpt2007(mShowBundle,"ShowBundle","Show detail of bundle results",{eTA2007::HDV}) - << AOpt2007(mThrsReject,"ThrRej","Threshold for rejection of outlayer, in pixel",{eTA2007::HDV}) - << AOpt2007(mMaxErrOK,"MaxErr","Max error acceptable for initial resection",{eTA2007::HDV}) - << mPhProj.DPPointsMeasures().ArgDirOutOpt("DirFiltered","Directory for filtered point") - */ - ; -} - -//================================================ - - -cAppli_SimulDispl::tImDispl cAppli_SimulDispl::GenerateSmoothRandDispl() -{ - tREAL8 aDeZoom = 10.0; - tREAL8 aNbBlob = 10.0; - - cPt2di aSzRed = Pt_round_up(ToR(mSz)/aDeZoom); - - tImDispl aResSsEch(aSzRed); - - for (const auto & aPix : aResSsEch.DIm()) - aResSsEch.DIm().SetV(aPix,RandUnif_C()); - - ExpFilterOfStdDev(aResSsEch.DIm(),5,Norm2(aSzRed)/aNbBlob); - NormalizedAvgDev(aResSsEch.DIm(),1e-10,2.0); - - tImDispl aRes(mSz); - for (const auto & aPix : aRes.DIm()) + cCollecSpecArg2007 &ArgObl(cCollecSpecArg2007 &anArgObl) override; + cCollecSpecArg2007 &ArgOpt(cCollecSpecArg2007 &anArgOpt) override; + + private: + tImDispl GenerateSmoothRandDispl(); + + // == Mandatory args ==== + std::string mNameImage; ///< name of the input image to deform + + // == Optionnal args ==== + + tREAL8 mAmplDef; + bool mWithDisc; + + // == Internal variables ==== + tIm mImIn; ///< memory representation of the image + tDIm *mDImIn; ///< memory representation of the image + cPt2di mSz; + tIm mImOut; ///< memory representation of the image + tDIm *mDImOut; ///< memory representation of the image + }; + + cAppli_SimulDispl::cAppli_SimulDispl( + const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) : cMMVII_Appli(aVArgs, aSpec), + mAmplDef(2.0), + mWithDisc(false), + mImIn(cPt2di(1, 1)), + mDImIn(nullptr), + mImOut(cPt2di(1, 1)), + mDImOut(nullptr) { - cPt2dr aPixSE = ToR(aPix)/aDeZoom; - aRes.DIm().SetV(aPix,aResSsEch.DIm().DefGetVBL(aPixSE,0)); } - - return aRes; -} - - -int cAppli_SimulDispl::Exe() -{ - - mImIn = tIm::FromFile(mNameImage); - cDataFileIm2D aDescFile = cDataFileIm2D::Create(mNameImage,false); - - mDImIn = & mImIn.DIm(); - mSz = mDImIn->Sz(); - - mImOut = tIm(mSz); - mDImOut = & mImOut.DIm(); - - - for (const auto & aPix : *mDImIn) + cCollecSpecArg2007 &cAppli_SimulDispl::ArgObl(cCollecSpecArg2007 &anArgObl) { - mDImOut->SetV(aPix,255-mDImIn->GetV(aPix)); + return anArgObl + << Arg2007(mNameImage, "Name of image to deform", {{eTA2007::FileImage}, + {eTA2007::FileDirProj}}); } - - tImDispl aImDispx = GenerateSmoothRandDispl(); - tImDispl aImDispy = GenerateSmoothRandDispl(); - tImDispl aImRegion = GenerateSmoothRandDispl(); - if (mWithDisc) + cCollecSpecArg2007 &cAppli_SimulDispl::ArgOpt(cCollecSpecArg2007 &anArgOpt) { - for (const auto & aPix : aImRegion.DIm()) - { - aImRegion.DIm().SetV(aPix,aImRegion.DIm().GetV(aPix)>0); - if (aImRegion.DIm().GetV(aPix) ) - { - std::swap(aImDispx.DIm().GetReference_V(aPix),aImDispy.DIm().GetReference_V(aPix)); - } - } - } - aImDispx.DIm().ToFile("DeplX.tif"); - aImDispy.DIm().ToFile("DeplY.tif"); - aImRegion.DIm().ToFile("Region.tif"); + return anArgOpt + << AOpt2007(mAmplDef, "Ampl", "Amplitude of deformation", {eTA2007::HDV}) + << AOpt2007(mWithDisc, "WithDisc", "Do we add disconinuities", {eTA2007::HDV}); + } + //================================================ - for (const auto & aPix : mImOut.DIm()) + cAppli_SimulDispl::tImDispl cAppli_SimulDispl::GenerateSmoothRandDispl() { - tREAL8 aDx = aImDispx.DIm().GetV(aPix); - tREAL8 aDy = aImDispy.DIm().GetV(aPix); - cPt2dr aPixR = ToR(aPix) + cPt2dr(aDx,aDy); + const tREAL8 aDeZoom = 10.0; + const tREAL8 aNbBlob = 10.0; - mDImOut->SetV(aPix,mDImIn->DefGetVBL(aPixR,0)); - } - - mDImOut->ToFile("Reech.tif",aDescFile.Type()); + const cPt2di aSzRed = Pt_round_up(ToR(mSz) / aDeZoom); - StdOut() << "hello , size of image = " << mImIn.DIm().Sz() << " " << mDImIn->Sz() << "\n"; - return EXIT_SUCCESS; -} + tImDispl aResSsEch(aSzRed); -/* ==================================================== */ -/* */ -/* MMVII */ -/* */ -/* ==================================================== */ + for (const cPt2di &aPix : aResSsEch.DIm()) + aResSsEch.DIm().SetV(aPix, RandUnif_C()); + ExpFilterOfStdDev(aResSsEch.DIm(), 5, Norm2(aSzRed) / aNbBlob); + NormalizedAvgDev(aResSsEch.DIm(), 1e-10, mAmplDef); -tMMVII_UnikPApli Alloc_SimulDispl(const std::vector & aVArgs,const cSpecMMVII_Appli & aSpec) -{ - return tMMVII_UnikPApli(new cAppli_SimulDispl(aVArgs,aSpec)); -} + tImDispl aRes(mSz); + for (const cPt2di &aPix : aRes.DIm()) + { + const cPt2dr aPixSE = ToR(aPix) / aDeZoom; + aRes.DIm().SetV(aPix, aResSsEch.DIm().DefGetVBL(aPixSE, 0)); + } -cSpecMMVII_Appli TheSpec_SimulDispl -( - "SimulDispl", - Alloc_SimulDispl, - "Generate smooth displacement and deformed image", - {eApF::ImProc}, - {eApDT::Image}, - {eApDT::Image}, - __FILE__ -); + return aRes; + } + int cAppli_SimulDispl::Exe() + { + mImIn = tIm::FromFile(mNameImage); + cDataFileIm2D aDescFile = cDataFileIm2D::Create(mNameImage, false); + + mDImIn = &mImIn.DIm(); + mSz = mDImIn->Sz(); + + mImOut = tIm(mSz); + mDImOut = &mImOut.DIm(); + + for (const cPt2di &aPix : *mDImIn) + mDImOut->SetV(aPix, 255 - mDImIn->GetV(aPix)); + + tImDispl aImDispx = GenerateSmoothRandDispl(); + tImDispl aImDispy = GenerateSmoothRandDispl(); + tImDispl aImRegion = GenerateSmoothRandDispl(); + + if (mWithDisc) + { + for (const cPt2di &aPix : aImRegion.DIm()) + { + aImRegion.DIm().SetV(aPix, aImRegion.DIm().GetV(aPix) > 0); + if (aImRegion.DIm().GetV(aPix)) + std::swap(aImDispx.DIm().GetReference_V(aPix), + aImDispy.DIm().GetReference_V(aPix)); + } + } + + aImDispx.DIm().ToFile("DeplX.tif"); + aImDispy.DIm().ToFile("DeplY.tif"); + aImRegion.DIm().ToFile("Region.tif"); + + for (const auto &aPix : mImOut.DIm()) + { + const tREAL8 aDx = aImDispx.DIm().GetV(aPix); + const tREAL8 aDy = aImDispy.DIm().GetV(aPix); + const cPt2dr aPixR = ToR(aPix) - cPt2dr(aDx, aDy); + + mDImOut->SetV(aPix, mDImIn->DefGetVBL(aPixR, 0)); + } + + mDImOut->ToFile("image_post.tif", aDescFile.Type()); + + StdOut() << "Size of image = [" << mImIn.DIm().Sz().x() + << ", " << mDImIn->SzY() << "]" << std::endl; + return EXIT_SUCCESS; + } -#if (0) -#endif + /* ====================================== */ + /* */ + /* MMVII */ + /* */ + /* ====================================== */ + tMMVII_UnikPApli Alloc_SimulDispl(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) + { + return tMMVII_UnikPApli(new cAppli_SimulDispl(aVArgs, aSpec)); + } + cSpecMMVII_Appli TheSpec_SimulDispl( + "SimulDispl", + Alloc_SimulDispl, + "Generate smooth displacement and deformed image", + {eApF::ImProc}, + {eApDT::Image}, + {eApDT::Image}, + __FILE__); }; // MMVII - diff --git a/MMVII/src/MeshDisplacement/TriangleDeformation.cpp b/MMVII/src/MeshDisplacement/TriangleDeformation.cpp new file mode 100644 index 0000000000..26c7f3f3c8 --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformation.cpp @@ -0,0 +1,921 @@ +#include "cMMVII_Appli.h" + +#include "MMVII_TplSymbTriangle.h" + +#include "TriangleDeformation.h" + +#include + +/** + \file TriangleDeformation.cpp + + \brief file for computing 2D deformations between 2 images + thanks to triangles. +**/ + +namespace MMVII +{ + /****************************************/ + /* */ + /* cPtInsideTriangles */ + /* */ + /****************************************/ + + cPtInsideTriangles::cPtInsideTriangles(const cTriangle2DCompiled &aCompTri, // a compiled triangle + const std::vector &aVectorFilledwithInsidePixels, // vector containing pixels insisde triangles + const size_t aFilledPixel, // a counter that is looping over pixels in triangles + const cDataIm2D &aDIm) // image + { + mFilledIndices = cPt2dr(aVectorFilledwithInsidePixels[aFilledPixel].x(), aVectorFilledwithInsidePixels[aFilledPixel].y()); + mBarycenterCoordinatesOfPixel = aCompTri.CoordBarry(mFilledIndices); + mValueOfPixel = aDIm.GetV(cPt2di(mFilledIndices.x(), mFilledIndices.y())); + } + + cPt3dr cPtInsideTriangles::GetBarycenterCoordinates() const { return mBarycenterCoordinatesOfPixel; } // Accessor + cPt2dr cPtInsideTriangles::GetCartesianCoordinates() const { return mFilledIndices; } // Accessor + tREAL8 cPtInsideTriangles::GetPixelValue() const { return mValueOfPixel; } // Accessor + + /******************************************/ + /* */ + /* cTriangleDeformation */ + /* */ + /******************************************/ + + cAppli_cTriangleDeformation::cAppli_cTriangleDeformation(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) : cMMVII_Appli(aVArgs, aSpec), + mRandomUniformLawUpperBoundLines(1), + mRandomUniformLawUpperBoundCols(1), + mShow(true), + mComputeAvgMax(false), + mUseMultiScaleApproach(true), + mSigmaGaussFilterStep(1), + mGenerateDisplacementImage(true), + mInitialiseWithPreviousIter(false), + mFreezeTranslationX(false), + mFreezeTranslationY(false), + mFreezeRadTranslation(false), + mFreezeRadScale(false), + mWeightRadTranslation(-1), + mWeightRadScale(-1), + mNumberOfIterGaussFilter(3), + mNumberOfEndIterations(2), + mFolderSaveResult(""), + mDisplayLastTranslationValues(false), + mSzImPre(cPt2di(1, 1)), + mImPre(mSzImPre), + mDImPre(nullptr), + mSzImPost(cPt2di(1, 1)), + mImPost(mSzImPost), + mDImPost(nullptr), + mSzImOut(cPt2di(1, 1)), + mImOut(mSzImOut), + mDImOut(nullptr), + mSzImDiff(cPt2di(1, 1)), + mImDiff(mSzImDiff), + mDImDiff(nullptr), + mSzImDepX(cPt2di(1, 1)), + mImDepX(mSzImDepX), + mDImDepX(nullptr), + mSzImDepY(cPt2di(1, 1)), + mImDepY(mSzImDepY), + mDImDepY(nullptr), + mVectorPts({cPt2dr(0, 0)}), + mDelTri(mVectorPts), + mSys(nullptr), + mEqTriDeform(nullptr) + { + mEqTriDeform = EqDeformTri(true, 1); // true means with derivative, 1 is size of buffer + } + + cAppli_cTriangleDeformation::~cAppli_cTriangleDeformation() + { + delete mSys; + delete mEqTriDeform; + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformation::ArgObl(cCollecSpecArg2007 &anArgObl) + { + return anArgObl + << Arg2007(mNamePreImage, "Name of pre-image file.", {{eTA2007::FileImage}, {eTA2007::FileDirProj}}) + << Arg2007(mNamePostImage, "Name of post-image file.", {eTA2007::FileImage}) + << Arg2007(mNumberPointsToGenerate, "Number of points you want to generate for triangulation.") + << Arg2007(mNumberOfScales, "Total number of scales to run in multi-scale approach or iterations if multi-scale approach is not applied in optimisation process."); + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformation::ArgOpt(cCollecSpecArg2007 &anArgOpt) + { + return anArgOpt + << AOpt2007(mRandomUniformLawUpperBoundCols, "RandomUniformLawUpperBoundXAxis", + "Maximum value that the uniform law can draw from on the x-axis.", {eTA2007::HDV}) + << AOpt2007(mRandomUniformLawUpperBoundLines, "RandomUniformLawUpperBoundYAxis", + "Maximum value that the uniform law can draw from for on the y-axis.", {eTA2007::HDV}) + << AOpt2007(mShow, "Show", "Whether to print minimisation results.", {eTA2007::HDV}) + << AOpt2007(mComputeAvgMax, "ComputeAvgMaxDiffIm", + "Whether to compute the average and maximum pixel value of the difference image between post and pre image or not.", {eTA2007::HDV}) + << AOpt2007(mUseMultiScaleApproach, "UseMultiScaleApproach", "Whether to use multi-scale approach or not.", {eTA2007::HDV}) + << AOpt2007(mSigmaGaussFilterStep, "SigmaGaussFilterStep", "Sigma value to use for Gauss filter in multi-stage approach.", {eTA2007::HDV}) + << AOpt2007(mGenerateDisplacementImage, "GenerateDisplacementImage", + "Whether to generate and save an image having been translated.", {eTA2007::HDV}) + << AOpt2007(mInitialiseWithPreviousIter, "InitialiseWithPreviousIteration", + "Whether to initialise or not the solution with the values obtained at the previous iteration.", {eTA2007::HDV}) + << AOpt2007(mFreezeTranslationX, "FreezeXTranslation", + "Whether to freeze or not x-translation to certain value during computation.", {eTA2007::HDV}) + << AOpt2007(mFreezeTranslationY, "FreezeYTranslation", + "Whether to freeze or not y-translation to certain value during computation.", {eTA2007::HDV}) + << AOpt2007(mFreezeRadTranslation, "FreezeRadTranslation", + "Whether to freeze radiometry translation factor in computation or not.", {eTA2007::HDV}) + << AOpt2007(mFreezeRadScale, "FreezeRadScaling", + "Whether to freeze radiometry scaling factor in computation or not.", {eTA2007::HDV}) + << AOpt2007(mWeightRadTranslation, "WeightRadiometryTranslation", + "A value to weight radiometry translation for soft freezing of coefficient.", {eTA2007::HDV}) + << AOpt2007(mWeightRadScale, "WeightRadiometryScaling", + "A value to weight radiometry scaling for soft freezing of coefficient.", {eTA2007::HDV}) + << AOpt2007(mNumberOfIterGaussFilter, "NumberOfIterationsGaussFilter", + "Number of iterations to run in Gauss filter algorithm.", {eTA2007::HDV}) + << AOpt2007(mNumberOfEndIterations, "NumberOfEndIterations", + "Number of iterations to run on original images in multi-scale approach.", {eTA2007::HDV}) + << AOpt2007(mFolderSaveResult, "FolderToSaveResults", + "Folder name where to store produced results", {eTA2007::HDV}) + << AOpt2007(mDisplayLastTranslationValues, "DisplayLastTranslationsValues", + "Whether to display the final values of unknowns linked to point translation.", {eTA2007::HDV}) + << AOpt2007(mDisplayLastRadiometryValues, "DisplayLastRadiometryValues", + "Whether to display or not the last values of radiometry unknowns after optimisation process.", {eTA2007::HDV}); + } + + void cAppli_cTriangleDeformation::ConstructUniformRandomVectorAndApplyDelaunay(std::vector aVectorPts, const int aNumberOfPointsToGenerate, + const int aRandomUniformLawUpperBoundLines, const int aRandomUniformLawUpperBoundCols, + cTriangulation2D &aDelaunayTri) + { + aVectorPts.pop_back(); // eliminate initialisation values + // Generate coordinates from drawing lines and columns of coordinates from a uniform distribution + for (int aNbPt = 0; aNbPt < aNumberOfPointsToGenerate; aNbPt++) + { + const tREAL8 aUniformRandomLine = RandUnif_N(aRandomUniformLawUpperBoundLines); + const tREAL8 aUniformRandomCol = RandUnif_N(aRandomUniformLawUpperBoundCols); + const cPt2dr aUniformRandomPt(aUniformRandomCol, aUniformRandomLine); // cPt2dr format + aVectorPts.push_back(aUniformRandomPt); + } + aDelaunayTri = aVectorPts; + + aDelaunayTri.MakeDelaunay(); // Delaunay triangulate randomly generated points. + } + + void cAppli_cTriangleDeformation::GeneratePointsForDelaunay(std::vector aVectorPts, const int aNumberOfPointsToGenerate, + int aRandomUniformLawUpperBoundLines, int aRandomUniformLawUpperBoundCols, + cTriangulation2D &aDelaunayTri, const cPt2di &aSzImPre) + { + // If user hasn't defined another value than the default value, it is changed + if (aRandomUniformLawUpperBoundLines == 1 && aRandomUniformLawUpperBoundCols == 1) + { + // Maximum value of coordinates are drawn from [0, NumberOfImageLines[ for lines + aRandomUniformLawUpperBoundLines = aSzImPre.y(); + // Maximum value of coordinates are drawn from [0, NumberOfImageColumns[ for columns + aRandomUniformLawUpperBoundCols = aSzImPre.x(); + } + else + { + if (aRandomUniformLawUpperBoundLines != 1 && aRandomUniformLawUpperBoundCols == 1) + aRandomUniformLawUpperBoundCols = aSzImPre.x(); + else + { + if (aRandomUniformLawUpperBoundLines == 1 && aRandomUniformLawUpperBoundCols != 1) + aRandomUniformLawUpperBoundLines = aSzImPre.y(); + } + } + + ConstructUniformRandomVectorAndApplyDelaunay(aVectorPts, aNumberOfPointsToGenerate, + aRandomUniformLawUpperBoundLines, aRandomUniformLawUpperBoundCols, + aDelaunayTri); + } + + void cAppli_cTriangleDeformation::InitialisationAfterExe(cTriangulation2D &aDelaunayTri, + cResolSysNonLinear *&aSys) + { + const size_t aStartNumberPts = 4 * aDelaunayTri.NbPts(); + tDenseVect aVInit(aStartNumberPts, eModeInitImage::eMIA_Null); + + for (size_t aStartKtNumber = 0; aStartKtNumber < aStartNumberPts; aStartKtNumber++) + { + if (aStartKtNumber % 4 == 3) + aVInit(aStartKtNumber) = 1; + } + + aSys = new cResolSysNonLinear(eModeSSR::eSSR_LsqDense, aVInit); + } + + void cAppli_cTriangleDeformation::InitialisationBeforeIteration(cTriangulation2D &aDelaunayTri, + cResolSysNonLinear *&aSys) + { + const size_t aINumberPts = 4 * aDelaunayTri.NbPts(); + tDenseVect aVInitBeforeIteration(aINumberPts, eModeInitImage::eMIA_Null); + + tDenseVect aIntermediateVCur = aSys->CurGlobSol(); // Get current solution. + + for (size_t aITr = 0; aITr < aDelaunayTri.NbFace(); aITr++) + { + const cPt3di aIntermediateIndicesOfTriKnots = aDelaunayTri.KthFace(aITr); + + const tIntVect aIntermediateVecInd = {4 * aIntermediateIndicesOfTriKnots.x(), 4 * aIntermediateIndicesOfTriKnots.x() + 1, + 4 * aIntermediateIndicesOfTriKnots.y(), 4 * aIntermediateIndicesOfTriKnots.y() + 1, + 4 * aIntermediateIndicesOfTriKnots.z(), 4 * aIntermediateIndicesOfTriKnots.z() + 1}; + + for (size_t aIKnotNumber=0; aIKnotNumber < aIntermediateVecInd.size(); aIKnotNumber += 2) + { + aVInitBeforeIteration(aIntermediateVecInd.at(aIKnotNumber)) = aIntermediateVCur(aIntermediateVecInd.at(aIKnotNumber)); + aVInitBeforeIteration(aIntermediateVecInd.at(aIKnotNumber + 1)) = aIntermediateVCur(aIntermediateVecInd.at(aIKnotNumber + 1)); + } + } + + aSys = new cResolSysNonLinear(eModeSSR::eSSR_LsqDense, aVInitBeforeIteration); + } + + void cAppli_cTriangleDeformation::SubtractPrePostImageAndComputeAvgAndMax() + { + mImDiff = tIm(mSzImPre); + mDImDiff = &mImDiff.DIm(); + + for (const cPt2di &aDiffPix : *mDImDiff) + mDImDiff->SetV(aDiffPix, mDImPre->GetV(aDiffPix) - mDImPost->GetV(aDiffPix)); + const int aNumberOfPixelsInImage = mSzImPre.x() * mSzImPre.y(); + + tREAL8 aSumPixelValuesInDiffImage = 0; + tREAL8 aMaxPixelValuesInDiffImage = 0; + tREAL8 aDiffImPixelValue = 0; + for (const cPt2di &aDiffPix : *mDImDiff) + { + aDiffImPixelValue = mDImDiff->GetV(aDiffPix); + aSumPixelValuesInDiffImage += aDiffImPixelValue; + if (aDiffImPixelValue > aMaxPixelValuesInDiffImage) + aMaxPixelValuesInDiffImage = aDiffImPixelValue; + } + StdOut() << "The average value of the difference image between the Pre and Post images is : " + << aSumPixelValuesInDiffImage / (tREAL8)aNumberOfPixelsInImage << std::endl; + StdOut() << "The maximum value of the difference image between the Pre and Post images is : " + << aMaxPixelValuesInDiffImage << std::endl; + } + + cPt2dr cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaToFilledPixel(const cPt2dr &aCurrentTranslationPointA, + const cPt2dr &aCurrentTranslationPointB, + const cPt2dr &aCurrentTranslationPointC, + const tDoubleVect &aVObs) + { + // apply current barycenter translation formula for x and y on current observations. + const tREAL8 aXTriCoord = aVObs[0] + aVObs[2] * aCurrentTranslationPointA.x() + aVObs[3] * aCurrentTranslationPointB.x() + + aVObs[4] * aCurrentTranslationPointC.x(); + const tREAL8 aYTriCoord = aVObs[1] + aVObs[2] * aCurrentTranslationPointA.y() + aVObs[3] * aCurrentTranslationPointB.y() + + aVObs[4] * aCurrentTranslationPointC.y(); + + const cPt2dr aCurrentTranslatedPixel = cPt2dr(aXTriCoord, aYTriCoord); + + return aCurrentTranslatedPixel; + } + + tREAL8 cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForTranslationRadiometry(const tREAL8 aCurrentRadTranslationPointA, + const tREAL8 aCurrentRadTranslationPointB, + const tREAL8 aCurrentRadTranslationPointC, + const tDoubleVect &aVObs) + { + const tREAL8 aCurentRadTranslation = aVObs[2] * aCurrentRadTranslationPointA + aVObs[3] * aCurrentRadTranslationPointB + + aVObs[4] * aCurrentRadTranslationPointC; + return aCurentRadTranslation; + } + + tREAL8 cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForScalingRadiometry(const tREAL8 aCurrentRadScalingPointA, + const tREAL8 aCurrentRadScalingPointB, + const tREAL8 aCurrentRadScalingPointC, + const tDoubleVect &aVObs) + { + const tREAL8 aCurrentRadScaling = aVObs[2] * aCurrentRadScalingPointA + aVObs[3] * aCurrentRadScalingPointB + + aVObs[4] * aCurrentRadScalingPointC; + return aCurrentRadScaling; + } + + void cAppli_cTriangleDeformation::LoadImageAndData(tIm &aCurIm, tDIm *&aCurDIm, const std::string &aPreOrPostImage, tIm &aImPre, tIm &aImPost) + { + (aPreOrPostImage == "pre") ? aCurIm = aImPre : aCurIm = aImPost; + aCurDIm = &aCurIm.DIm(); + } + + bool cAppli_cTriangleDeformation::ManageDifferentCasesOfEndIterations(const int aIterNumber, const int aNumberOfScales, const int aNumberOfEndIterations, + bool aIsLastIters, tIm &aImPre, tIm &aImPost, tIm aCurPreIm, tDIm *aCurPreDIm, + tIm aCurPostIm, tDIm *aCurPostDIm) + { + switch (aNumberOfEndIterations) + { + case 1: // one last iteration + if (aIterNumber == aNumberOfScales) + { + aIsLastIters = true; + LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", aImPre, aImPost); + LoadImageAndData(aCurPostIm, aCurPostDIm, "post", aImPre, aImPost); + } + break; + case 2: // two last iterations + if ((aIterNumber == aNumberOfScales) || (aIterNumber == aNumberOfScales + aNumberOfEndIterations - 1)) + { + aIsLastIters = true; + LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", aImPre, aImPost); + LoadImageAndData(aCurPostIm, aCurPostDIm, "post", aImPre, aImPost); + } + break; + case 3: // three last iterations + if ((aIterNumber == aNumberOfScales) || (aIterNumber == mNumberOfScales + aNumberOfEndIterations - 2) || + (aIterNumber == aNumberOfScales + aNumberOfEndIterations - 1)) + { + aIsLastIters = true; + LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", aImPre, aImPost); + LoadImageAndData(aCurPostIm, aCurPostDIm, "post", aImPre, aImPost); + } + break; + default: // default is two last iterations + if ((aIterNumber == aNumberOfScales) || (aIterNumber == aNumberOfScales + aNumberOfEndIterations - 1)) + { + aIsLastIters = true; + LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", aImPre, aImPost); + LoadImageAndData(aCurPostIm, aCurPostDIm, "post", aImPre, aImPost); + } + break; + } + return aIsLastIters; + } + + void cAppli_cTriangleDeformation::LoopOverTrianglesAndUpdateParameters(const int aIterNumber, const bool aUserDefinedFolderName) + { + //----------- allocate vec of obs : + tDoubleVect aVObs(12, 0.0); // 6 for ImagePre and 6 for ImagePost + + //----------- extract current parameters + tDenseVect aVCur = mSys->CurGlobSol(); // Get current solution. + + /* + for (int aUnk=0; aUnkToFile(mFolderSaveResult + "/GaussFilteredImPre_iter_" + std::to_string(aIterNumber) + ".tif"); + else + aCurPreDIm->ToFile("GaussFilteredImPre_iter_" + std::to_string(aIterNumber) + ".tif"); + } + } + else if (mUseMultiScaleApproach && mIsLastIters) + { + LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + //----------- declaration of indicators of convergence + tREAL8 aSomDif = 0; // sum of difference between untranslated pixel and translated one. + size_t aNbOut = 0; // number of translated pixels out of image + + // Count number of pixels inside triangles for normalisation + size_t aTotalNumberOfInsidePixels = 0; + + // hard constraint : freeze radiometric coefficients + if (mFreezeRadTranslation || mFreezeRadScale || + mFreezeTranslationX || mFreezeTranslationY) + { + for (size_t aTr = 0; aTr < mDelTri.NbFace(); aTr++) + { + const cPt3di aIndicesOfTriKnots = mDelTri.KthFace(aTr); + + const tIntVect aVecInd = {4 * aIndicesOfTriKnots.x(), 4 * aIndicesOfTriKnots.x() + 1, + 4 * aIndicesOfTriKnots.x() + 2, 4 * aIndicesOfTriKnots.x() + 3, + 4 * aIndicesOfTriKnots.y(), 4 * aIndicesOfTriKnots.y() + 1, + 4 * aIndicesOfTriKnots.y() + 2, 4 * aIndicesOfTriKnots.y() + 3, + 4 * aIndicesOfTriKnots.z(), 4 * aIndicesOfTriKnots.z() + 1, + 4 * aIndicesOfTriKnots.z() + 2, 4 * aIndicesOfTriKnots.z() + 3}; + + if (mFreezeTranslationX) + { + mSys->SetFrozenVar(aVecInd.at(2), aVCur(0)); + mSys->SetFrozenVar(aVecInd.at(6), aVCur(4)); + mSys->SetFrozenVar(aVecInd.at(10), aVCur(8)); + } + if (mFreezeTranslationY) + { + mSys->SetFrozenVar(aVecInd.at(2), aVCur(1)); + mSys->SetFrozenVar(aVecInd.at(6), aVCur(5)); + mSys->SetFrozenVar(aVecInd.at(10), aVCur(9)); + } + if (mFreezeRadTranslation) + { + mSys->SetFrozenVar(aVecInd.at(2), aVCur(2)); + mSys->SetFrozenVar(aVecInd.at(6), aVCur(6)); + mSys->SetFrozenVar(aVecInd.at(10), aVCur(10)); + } + if (mFreezeRadScale) + { + mSys->SetFrozenVar(aVecInd.at(3), aVCur(3)); + mSys->SetFrozenVar(aVecInd.at(7), aVCur(7)); + mSys->SetFrozenVar(aVecInd.at(11), aVCur(11)); + } + } + } + + // Loop over all triangles to add the observations on each point + for (size_t aTr = 0; aTr < mDelTri.NbFace(); aTr++) + { + const tTri2dr aTri = mDelTri.KthTri(aTr); + const cPt3di aIndicesOfTriKnots = mDelTri.KthFace(aTr); + + const cTriangle2DCompiled aCompTri(aTri); + + std::vector aVectorToFillWithInsidePixels; + aCompTri.PixelsInside(aVectorToFillWithInsidePixels); // get pixels inside triangle + + //----------- index of unknown, finds the associated pixels of current triangle + const tIntVect aVecInd = {4 * aIndicesOfTriKnots.x(), 4 * aIndicesOfTriKnots.x() + 1, + 4 * aIndicesOfTriKnots.x() + 2, 4 * aIndicesOfTriKnots.x() + 3, + 4 * aIndicesOfTriKnots.y(), 4 * aIndicesOfTriKnots.y() + 1, + 4 * aIndicesOfTriKnots.y() + 2, 4 * aIndicesOfTriKnots.y() + 3, + 4 * aIndicesOfTriKnots.z(), 4 * aIndicesOfTriKnots.z() + 1, + 4 * aIndicesOfTriKnots.z() + 2, 4 * aIndicesOfTriKnots.z() + 3}; + + const cPt2dr aCurTrPointA = cPt2dr(aVCur(aVecInd.at(0)), + aVCur(aVecInd.at(1))); // current translation 1st point of triangle + const cPt2dr aCurTrPointB = cPt2dr(aVCur(aVecInd.at(4)), + aVCur(aVecInd.at(5))); // current translation 2nd point of triangle + const cPt2dr aCurTrPointC = cPt2dr(aVCur(aVecInd.at(8)), + aVCur(aVecInd.at(9))); // current translation 3rd point of triangle + + const tREAL8 aCurRadTrPointA = aVCur(aVecInd.at(2)); // current translation on radiometry 1st point of triangle + const tREAL8 aCurRadScPointA = aVCur(aVecInd.at(3)); // current scale on radiometry 3rd point of triangle + const tREAL8 aCurRadTrPointB = aVCur(aVecInd.at(6)); // current translation on radiometry 2nd point of triangle + const tREAL8 aCurRadScPointB = aVCur(aVecInd.at(7)); // current scale on radiometry 3rd point of triangle + const tREAL8 aCurRadTrPointC = aVCur(aVecInd.at(10)); // current translation on radiometry 3rd point of triangle + const tREAL8 aCurRadScPointC = aVCur(aVecInd.at(11)); // current scale on radiometry 3rd point of triangle + + // soft constraint radiometric translation + if (!mFreezeRadTranslation) + { + if (mWeightRadTranslation >= 0) + { + const int aSolStep = 4; // adapt step to solution vector configuration + const int aSolStart = 2; + for (size_t aIndCurSol = aSolStart; aIndCurSol < aVecInd.size() - 1; aIndCurSol += aSolStep) + { + const int aIndices = aVecInd.at(aIndCurSol); + mSys->AddEqFixVar(aIndices, aVCur(aIndices), mWeightRadTranslation); + } + } + } + + // soft constraint radiometric scaling + if (!mFreezeRadScale) + { + if (mWeightRadScale >= 0) + { + const int aSolStep = 4; // adapt step to solution vector configuration + const int aSolStart = 3; + for (size_t aIndCurSol = aSolStart; aIndCurSol < aVecInd.size(); aIndCurSol += aSolStep) + { + const int aIndices = aVecInd.at(aIndCurSol); + mSys->AddEqFixVar(aIndices, aVCur(aIndices), mWeightRadScale); + } + } + } + + const size_t aNumberOfInsidePixels = aVectorToFillWithInsidePixels.size(); + + // Loop over all pixels inside triangle + // size_t is necessary as there can be a lot of pixels in triangles + for (size_t aFilledPixel = 0; aFilledPixel < aNumberOfInsidePixels; aFilledPixel++) + { + const cPtInsideTriangles aPixInsideTriangle = cPtInsideTriangles(aCompTri, aVectorToFillWithInsidePixels, + aFilledPixel, *aCurPreDIm); + // prepare for barycenter translation formula by filling aVObs with different coordinates + FormalInterpBarycenter_SetObs(aVObs, 0, aPixInsideTriangle); + + // image of a point in triangle by current translation + const cPt2dr aTranslatedFilledPoint = ApplyBarycenterTranslationFormulaToFilledPixel(aCurTrPointA, aCurTrPointB, + aCurTrPointC, aVObs); + // radiometry translation of pixel by current radiometry translation of triangle knots + const tREAL8 aRadiometryTranslation = ApplyBarycenterTranslationFormulaForTranslationRadiometry(aCurRadTrPointA, + aCurRadTrPointB, + aCurRadTrPointC, + aVObs); + // radiometry translation of pixel by current radiometry scaling of triangle knots + const tREAL8 aRadiometryScaling = ApplyBarycenterTranslationFormulaForScalingRadiometry(aCurRadScPointA, + aCurRadScPointB, + aCurRadScPointC, + aVObs); + + if (aCurPostDIm->InsideBL(aTranslatedFilledPoint)) // avoid errors + { + // prepare for application of bilinear formula + FormalBilinTri_SetObs(aVObs, TriangleDisplacement_NbObs, aTranslatedFilledPoint, *aCurPostDIm); + + // Now add observation + mSys->CalcAndAddObs(mEqTriDeform, aVecInd, aVObs); + + // compute indicators + const tREAL8 aRadiomValueImPre = aRadiometryScaling * aVObs[5] + aRadiometryTranslation; + const tREAL8 aDif = aRadiomValueImPre - aCurPostDIm->GetVBL(aTranslatedFilledPoint); // residual + aSomDif += std::abs(aDif); + } + else + aNbOut++; // Count number of pixels translated outside post image + + aTotalNumberOfInsidePixels += aNumberOfInsidePixels; + } + } + + if (mUseMultiScaleApproach && !mIsLastIters && aIterNumber != 0) + { + const bool aGenerateIntermediateMaps = false; // if a generating intermediate displacement maps is wanted + if (aGenerateIntermediateMaps) + GenerateDisplacementMapsAndOutputImages(aVCur, aIterNumber, aUserDefinedDir); + } + + // Update all parameter taking into account previous observation + mSys->SolveUpdateReset(); + + if (mShow) + StdOut() << aIterNumber + 1 << ", " << aSomDif / aTotalNumberOfInsidePixels + << ", " << aNbOut << std::endl; + } + + void cAppli_cTriangleDeformation::FillDisplacementMapsAndOutputImage(const cPtInsideTriangles &aLastPixInsideTriangle, + const cPt2dr &aLastTranslatedFilledPoint, + const tREAL8 aLastRadiometryTranslation, + const tREAL8 aLastRadiometryScaling) + { + const tREAL8 aLastXCoordinate = aLastPixInsideTriangle.GetCartesianCoordinates().x(); + const tREAL8 aLastYCoordinate = aLastPixInsideTriangle.GetCartesianCoordinates().y(); + const tREAL8 aLastPixelValue = aLastPixInsideTriangle.GetPixelValue(); + + const cPt2di aLastCoordinate = cPt2di(aLastXCoordinate, aLastYCoordinate); + mDImDepX->SetV(aLastCoordinate, + aLastTranslatedFilledPoint.x() - aLastXCoordinate); + mDImDepY->SetV(aLastCoordinate, + aLastTranslatedFilledPoint.y() - aLastYCoordinate); + const tREAL8 aLastXTranslatedCoord = aLastXCoordinate + mDImDepX->GetV(aLastCoordinate); + const tREAL8 aLastYTranslatedCoord = aLastYCoordinate + mDImDepY->GetV(aLastCoordinate); + + const tREAL8 aLastRadiometryValue = aLastRadiometryScaling * aLastPixelValue + + aLastRadiometryTranslation; + + // Build image with intensities displaced + // deal with different cases of pixel being translated out of image + if (aLastXTranslatedCoord < 0 && aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, 0))); + else if (aLastXTranslatedCoord >= mSzImOut.x() && aLastYTranslatedCoord >= mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, mSzImOut.y() - 1))); + else if (aLastXTranslatedCoord < 0 && aLastYTranslatedCoord >= mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, mSzImOut.y() - 1))); + else if (aLastXTranslatedCoord >= mSzImOut.x() && aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, 0))); + else if (aLastXTranslatedCoord >= 0 && aLastXTranslatedCoord < mSzImOut.x() && + aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(aLastXTranslatedCoord, 0))); + else if (aLastXTranslatedCoord >= 0 && aLastXTranslatedCoord < mSzImOut.x() && + aLastYTranslatedCoord > mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(aLastXTranslatedCoord, mSzImOut.y() - 1))); + else if (aLastYTranslatedCoord >= 0 && aLastYTranslatedCoord < mSzImOut.y() && + aLastXTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, aLastYTranslatedCoord))); + else if (aLastYTranslatedCoord >= 0 && aLastYTranslatedCoord < mSzImOut.y() && + aLastXTranslatedCoord > mSzImOut.x()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, aLastYTranslatedCoord))); + else + // at the translated pixel the untranslated pixel value is given computed with the right radiometry values + mDImOut->SetV(cPt2di(aLastXTranslatedCoord, aLastYTranslatedCoord), aLastRadiometryValue); + } + + void cAppli_cTriangleDeformation::GenerateDisplacementMapsAndOutputImages(const tDenseVect &aVFinalSol, const int aIterNumber, + const bool aUserDefinedFolderName) + { + mImOut = tIm(mSzImPre); + mDImOut = &mImOut.DIm(); + mSzImOut = cPt2di(mDImOut->Sz().x(), mDImOut->Sz().y()); + + mImDepX = tIm(mSzImPre, 0, eModeInitImage::eMIA_Null); + mDImDepX = &mImDepX.DIm(); + + mImDepY = tIm(mSzImPre, 0, eModeInitImage::eMIA_Null); + mDImDepY = &mImDepY.DIm(); + + tIm aLastPreIm = tIm(mSzImPre); + tDIm *aLastPreDIm = nullptr; + LoadImageAndData(aLastPreIm, aLastPreDIm, "pre", mImPre, mImPost); + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aLastPreIm = mImPre.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aLastPreDIm = &aLastPreIm.DIm(); + } + + tDoubleVect aLastVObs(12, 0.0); + + for (const cPt2di &aOutPix : *mDImOut) // Initialise output images + mDImOut->SetV(aOutPix, aLastPreDIm->GetV(aOutPix)); + + for (size_t aLTr = 0; aLTr < mDelTri.NbFace(); aLTr++) + { + const tTri2dr aLastTri = mDelTri.KthTri(aLTr); + const cPt3di aLastIndicesOfTriKnots = mDelTri.KthFace(aLTr); + + const cTriangle2DCompiled aLastCompTri(aLastTri); + + std::vector aLastVectorToFillWithInsidePixels; + aLastCompTri.PixelsInside(aLastVectorToFillWithInsidePixels); + + const tIntVect aLastVecInd = { + 4 * aLastIndicesOfTriKnots.x(), + 4 * aLastIndicesOfTriKnots.x() + 1, + 4 * aLastIndicesOfTriKnots.x() + 2, + 4 * aLastIndicesOfTriKnots.x() + 3, + 4 * aLastIndicesOfTriKnots.y(), + 4 * aLastIndicesOfTriKnots.y() + 1, + 4 * aLastIndicesOfTriKnots.y() + 2, + 4 * aLastIndicesOfTriKnots.y() + 3, + 4 * aLastIndicesOfTriKnots.z(), + 4 * aLastIndicesOfTriKnots.z() + 1, + 4 * aLastIndicesOfTriKnots.z() + 2, + 4 * aLastIndicesOfTriKnots.z() + 3, + }; + + const cPt2dr aLastTrPointA = cPt2dr(aVFinalSol(aLastVecInd.at(0)), + aVFinalSol(aLastVecInd.at(1))); // last translation 1st point of triangle + const cPt2dr aLastTrPointB = cPt2dr(aVFinalSol(aLastVecInd.at(4)), + aVFinalSol(aLastVecInd.at(5))); // last translation 2nd point of triangle + const cPt2dr aLastTrPointC = cPt2dr(aVFinalSol(aLastVecInd.at(8)), + aVFinalSol(aLastVecInd.at(9))); // last translation 3rd point of triangle + + const tREAL8 aLastRadTrPointA = aVFinalSol(aLastVecInd.at(2)); // last translation on radiometry 1st point of triangle + const tREAL8 aLastRadScPointA = aVFinalSol(aLastVecInd.at(3)); // last scale on radiometry 3rd point of triangle + const tREAL8 aLastRadTrPointB = aVFinalSol(aLastVecInd.at(6)); // last translation on radiometry 2nd point of triangle + const tREAL8 aLastRadScPointB = aVFinalSol(aLastVecInd.at(7)); // last scale on radiometry 3rd point of triangle + const tREAL8 aLastRadTrPointC = aVFinalSol(aLastVecInd.at(10)); // last translation on radiometry 3rd point of triangle + const tREAL8 aLastRadScPointC = aVFinalSol(aLastVecInd.at(11)); // last scale on radiometry 3rd point of triangle + + const size_t aLastNumberOfInsidePixels = aLastVectorToFillWithInsidePixels.size(); + + for (size_t aLastFilledPixel = 0; aLastFilledPixel < aLastNumberOfInsidePixels; aLastFilledPixel++) + { + const cPtInsideTriangles aLastPixInsideTriangle = cPtInsideTriangles(aLastCompTri, aLastVectorToFillWithInsidePixels, + aLastFilledPixel, *aLastPreDIm); + // prepare for barycenter translation formula by filling aVObs with different coordinates + FormalInterpBarycenter_SetObs(aLastVObs, 0, aLastPixInsideTriangle); + + // image of a point in triangle by current translation + const cPt2dr aLastTranslatedFilledPoint = ApplyBarycenterTranslationFormulaToFilledPixel(aLastTrPointA, aLastTrPointB, + aLastTrPointC, aLastVObs); + + const tREAL8 aLastRadiometryTranslation = ApplyBarycenterTranslationFormulaForTranslationRadiometry(aLastRadTrPointA, + aLastRadTrPointB, + aLastRadTrPointC, + aLastVObs); + + const tREAL8 aLastRadiometryScaling = ApplyBarycenterTranslationFormulaForScalingRadiometry(aLastRadScPointA, + aLastRadScPointB, + aLastRadScPointC, + aLastVObs); + + FillDisplacementMapsAndOutputImage(aLastPixInsideTriangle, aLastTranslatedFilledPoint, + aLastRadiometryTranslation, aLastRadiometryScaling); + } + } + + // save displacement maps in x and y to image files + if (mUseMultiScaleApproach) + { + if (aUserDefinedFolderName) + { + mDImDepX->ToFile(mFolderSaveResult + "/DisplacedPixelsX_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + mDImDepY->ToFile(mFolderSaveResult + "/DisplacedPixelsY_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + } + else + { + mDImDepX->ToFile("DisplacedPixelsX_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + mDImDepY->ToFile("DisplacedPixelsY_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + } + if (aIterNumber == mNumberOfScales + mNumberOfEndIterations - 1) + { + if (aUserDefinedFolderName) + mDImOut->ToFile(mFolderSaveResult + "/DisplacedPixels_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + else + mDImOut->ToFile("DisplacedPixels_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + } + } + else + { + if (aUserDefinedFolderName) + { + mDImDepX->ToFile(mFolderSaveResult + "/DisplacedPixelsX_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + mDImDepY->ToFile(mFolderSaveResult + "/DisplacedPixelsY_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + mDImOut->ToFile(mFolderSaveResult + "/DisplacedPixels_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + } + else + { + mDImDepX->ToFile("DisplacedPixelsX_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + mDImDepY->ToFile("DisplacedPixelsY_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + mDImOut->ToFile("DisplacedPixels_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + } + } + } + + void cAppli_cTriangleDeformation::DisplayLastUnknownValues(const tDenseVect &aVFinalSol, const bool aDisplayLastRadiometryValues, + const bool aDisplayLastTranslationValues) + { + if (aDisplayLastRadiometryValues && aDisplayLastTranslationValues) + { + for (int aFinalUnk = 0; aFinalUnk < aVFinalSol.DIm().Sz(); aFinalUnk++) + { + StdOut() << aVFinalSol(aFinalUnk) << " "; + if (aFinalUnk % 4 == 3 && aFinalUnk != 0) + StdOut() << std::endl; + } + } + if (aDisplayLastRadiometryValues && !aDisplayLastTranslationValues) + { + for (int aFinalUnk = 0; aFinalUnk < aVFinalSol.DIm().Sz(); aFinalUnk++) + { + if (aFinalUnk % 4 == 0 || aFinalUnk % 4 == 1) + StdOut() << aVFinalSol(aFinalUnk) << " "; + if (aFinalUnk % 4 == 3 && aFinalUnk != 0) + StdOut() << std::endl; + } + } + else if (aDisplayLastRadiometryValues && !aDisplayLastTranslationValues) + { + for (int aFinalUnk = 0; aFinalUnk < aVFinalSol.DIm().Sz(); aFinalUnk++) + { + if (aFinalUnk % 4 == 2 || aFinalUnk % 4 == 3) + StdOut() << aVFinalSol(aFinalUnk) << " "; + if (aFinalUnk % 4 == 3 && aFinalUnk != 0) + StdOut() << std::endl; + } + } + } + + void cAppli_cTriangleDeformation::GenerateDisplacementMapsAndDisplayLastValuesUnknowns(const int aIterNumber, const bool aDisplayLastRadiometryValues, + const bool aDisplayLastTranslationValues, const bool aUserDefinedFolderName) + { + tDenseVect aVFinalSol = mSys->CurGlobSol(); + + if (mGenerateDisplacementImage) + GenerateDisplacementMapsAndOutputImages(aVFinalSol, aIterNumber, aUserDefinedFolderName); + + if (aDisplayLastRadiometryValues || aDisplayLastTranslationValues) + DisplayLastUnknownValues(aVFinalSol, aDisplayLastRadiometryValues, aDisplayLastTranslationValues); + } + + void cAppli_cTriangleDeformation::DoOneIteration(const int aIterNumber, const bool aUserDefinedFolderName) + { + LoopOverTrianglesAndUpdateParameters(aIterNumber, aUserDefinedFolderName); // Iterate over triangles and solve system + + // Show final translation results and produce displacement maps + if (mUseMultiScaleApproach) + { + if (aIterNumber == (mNumberOfScales + mNumberOfEndIterations - 1)) + GenerateDisplacementMapsAndDisplayLastValuesUnknowns(aIterNumber, mDisplayLastRadiometryValues, + mDisplayLastTranslationValues, aUserDefinedFolderName); + } + else + { + if (aIterNumber == (mNumberOfScales - 1)) + GenerateDisplacementMapsAndDisplayLastValuesUnknowns(aIterNumber, mDisplayLastRadiometryValues, + mDisplayLastTranslationValues, aUserDefinedFolderName); + } + } + + //----------------------------------------- + + int cAppli_cTriangleDeformation::Exe() + { + // read pre and post images and their sizes + mImPre = tIm::FromFile(mNamePreImage); + mImPost = tIm::FromFile(mNamePostImage); + + mDImPre = &mImPre.DIm(); + mSzImPre = mDImPre->Sz(); + + mDImPost = &mImPost.DIm(); + mSzImPost = mDImPost->Sz(); + + if (mUseMultiScaleApproach) + mSigmaGaussFilter = mNumberOfScales * mSigmaGaussFilterStep; + + if (mComputeAvgMax) + SubtractPrePostImageAndComputeAvgAndMax(); + + if (mShow) + StdOut() << "Iter, " + << "Diff, " + << "NbOut" << std::endl; + + // Generate triangulated knots coordinates + GeneratePointsForDelaunay(mVectorPts, mNumberPointsToGenerate, mRandomUniformLawUpperBoundLines, + mRandomUniformLawUpperBoundCols, mDelTri, mSzImPre); + + // Initialise unknown values of problem + InitialisationAfterExe(mDelTri, mSys); + + bool aUseDefinedFolderName = false; + if (!mFolderSaveResult.empty()) + { + aUseDefinedFolderName = true; + if (!std::filesystem::exists(mFolderSaveResult)) + std::filesystem::create_directory(mFolderSaveResult); + } + + if (mUseMultiScaleApproach) + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales + mNumberOfEndIterations; aIterNumber++) + { + DoOneIteration(aIterNumber, aUseDefinedFolderName); + if (mInitialiseWithPreviousIter) + InitialisationBeforeIteration(mDelTri, mSys); + } + } + else + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales; aIterNumber++) + { + DoOneIteration(aIterNumber, aUseDefinedFolderName); + if (mInitialiseWithPreviousIter) + InitialisationBeforeIteration(mDelTri, mSys); + } + } + + return EXIT_SUCCESS; + } + + /********************************************/ + // ::MMVII // + /********************************************/ + + tMMVII_UnikPApli Alloc_cTriangleDeformation(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) + { + return tMMVII_UnikPApli(new cAppli_cTriangleDeformation(aVArgs, aSpec)); + } + + cSpecMMVII_Appli TheSpec_ComputeTriangleDeformation( + "ComputeTriangleDeformation", + Alloc_cTriangleDeformation, + "Compute 2D deformation of triangles between images using triangles", + {eApF::ImProc}, // category + {eApDT::Image}, // input + {eApDT::Image}, // output + __FILE__); + +}; // namespace MMVII \ No newline at end of file diff --git a/MMVII/src/MeshDisplacement/TriangleDeformation.h b/MMVII/src/MeshDisplacement/TriangleDeformation.h new file mode 100644 index 0000000000..ed48e4f7e8 --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformation.h @@ -0,0 +1,180 @@ +#ifndef _TRIANGLEDEFORMATION_H_ +#define _TRIANGLEDEFORMATION_H_ + +#include "MMVII_Geom2D.h" +#include "MMVII_PhgrDist.h" + +using namespace NS_SymbolicDerivative; + +namespace MMVII +{ + /****************************************/ + /* */ + /* cPtInsideTriangles */ + /* */ + /****************************************/ + + class cPtInsideTriangles + { + public: + cPtInsideTriangles(const cTriangle2DCompiled &aCompTri, // a compiled triangle + const std::vector &aVectorFilledwithInsidePixels, // vector containing pixels inside triangles + const size_t aFilledPixel, // a counter that is looping over pixels in triangles + const cDataIm2D &aDIm); // image + cPt3dr GetBarycenterCoordinates() const; // Accessor for barycenter coordinates + cPt2dr GetCartesianCoordinates() const; // Accessor for cartesian coordinates + tREAL8 GetPixelValue() const; // Accessor for pixel value at coordinates + + private: + cPt3dr mBarycenterCoordinatesOfPixel; // Barycentric coordinates of pixel. + cPt2dr mFilledIndices; // 2D cartesian coordinates of pixel. + tREAL8 mValueOfPixel; // Intensity in image at pixel. + }; + + /******************************************/ + /* */ + /* cTriangleDeformation */ + /* */ + /******************************************/ + + class cAppli_cTriangleDeformation : public cMMVII_Appli + { + public: + typedef cIm2D tIm; + typedef cDataIm2D tDIm; + typedef cTriangle tTri2dr; + typedef cDenseVect tDenseVect; + typedef std::vector tIntVect; + typedef std::vector tDoubleVect; + + cAppli_cTriangleDeformation(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec); + ~cAppli_cTriangleDeformation(); + + int Exe() override; + cCollecSpecArg2007 &ArgObl(cCollecSpecArg2007 &anArgObl) override; + cCollecSpecArg2007 &ArgOpt(cCollecSpecArg2007 &anArgOpt) override; + + // Build uniform vector of coordinates and apply Delaunay triangulation + void ConstructUniformRandomVectorAndApplyDelaunay(std::vector aVectorPts, const int aNumberOfPointsToGenerate, + const int aRandomUniformLawUpperBoundLines, const int aRandomUniformLawUpperBoundCols, + cTriangulation2D &aDelaunayTri); + // Generate coordinates from uniform law for Delaunay triangulation application + void GeneratePointsForDelaunay(std::vector aVectorPts, const int aNumberOfPointsToGenerate, + int aRandomUniformLawUpperBoundLines, int aRandomUniformLawUpperBoundCols, + cTriangulation2D &aDelaunayTri, const cPt2di &aSzImPre); + // Construct difference image and compute average and max pixel value on ths image + void SubtractPrePostImageAndComputeAvgAndMax(); + // Iterate of triangles and inside pixels + virtual void DoOneIteration(const int aIterNumber, const bool aUserDefinedFolderName); + // Loops over all triangles and solves system to update parameters at end of iteration + virtual void LoopOverTrianglesAndUpdateParameters(const int aIterNumber, const bool aUserDefinedFolderName); + // Generate displacement maps of last solution + virtual void GenerateDisplacementMapsAndOutputImages(const tDenseVect &aVFinalSol, const int aIterNumber, + const bool aUserDefinedFolderName); + // Display values of unknowns at last iteration of optimisation process + void DisplayLastUnknownValues(const tDenseVect &aVFinalSol, const bool aDisplayLastRadiometryValues, + const bool aDisplayLastTranslationValues); + // Generates Displacement maps and coordinates of points in triangulation at last iteration + virtual void GenerateDisplacementMapsAndDisplayLastValuesUnknowns(const int aIterNumber, const bool aDisplayLastRadiometryValues, + const bool aDisplayLastTranslationValues, const bool aUserDefinedFolderName); + // Initialise values of unknowns at the beginning of optimsation process after user has input information + void InitialisationAfterExe(cTriangulation2D &aDelaunayTri, + cResolSysNonLinear *&aSys); + // Initialise values of unknowns before current iteration + void InitialisationBeforeIteration(cTriangulation2D &aDelaunayTri, cResolSysNonLinear *&aSys); + // Loads current pre and post images + void LoadImageAndData(tIm &aCurIm, tDIm *&aCurDIm, const std::string &aPreOrPostImage, tIm &aImPre, tIm &aImPost); + // Load image and data according to number of iterations to optimise on original image + bool ManageDifferentCasesOfEndIterations(const int aIterNumber, const int aNumberOfScales, const int aNumberOfEndIterations, + bool aIsLastIters, tIm &aImPre, tIm &aImPost, tIm aCurPreIm, tDIm *aCurPreDIm, + tIm aCurPostIm, tDIm *aCurPostDIm); + // Fill displacement maps and output image + virtual void FillDisplacementMapsAndOutputImage(const cPtInsideTriangles &aLastPixInsideTriangle, + const cPt2dr &aLastTranslatedFilledPoint, + const tREAL8 aLastRadiometryTranslation, + const tREAL8 aLastRadiometryScaling); + // Apply barycentric translation formula to current translation values + cPt2dr ApplyBarycenterTranslationFormulaToFilledPixel(const cPt2dr &aCurrentTranslationPointA, + const cPt2dr &aCurrentTranslationPointB, + const cPt2dr &aCurrentTranslationPointC, + const tDoubleVect &aVObs); + // Apply barycentric translation formula to current radiometric translation values + tREAL8 ApplyBarycenterTranslationFormulaForTranslationRadiometry(const tREAL8 aCurrentRadTranslationPointA, + const tREAL8 aCurrentRadTranslationPointB, + const tREAL8 aCurrentRadTranslationPointC, + const tDoubleVect &aVObs); + // Apply barycentric translation formula to current radiometric scaling values + tREAL8 ApplyBarycenterTranslationFormulaForScalingRadiometry(const tREAL8 aCurrentRadScalingPointA, + const tREAL8 aCurrentRadScalingPointB, + const tREAL8 aCurrentRadScalingPointC, + const tDoubleVect &aVObs); + + private: + // == Mandatory args ==== + + std::string mNamePreImage; // Name of given pre-image + std::string mNamePostImage; // Name of given post-image + int mNumberPointsToGenerate; // number of generated points + int mNumberOfScales; // number of iterations in optimisation process + + // == Optionnal args ==== + + int mRandomUniformLawUpperBoundLines; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundLines [ + int mRandomUniformLawUpperBoundCols; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundCols [ + bool mShow; // Print result, export image ... + bool mComputeAvgMax; // Compute average and maximum pixel value of difference image between pre and post images + bool mUseMultiScaleApproach; // Apply multi-scale approach or not + int mSigmaGaussFilterStep; // Decreasing step of sigma value during iterations + bool mGenerateDisplacementImage; // Generate image with displaced pixels + bool mInitialiseWithPreviousIter; // Initialise unknown value at beginning of every iteration with values obtained at previous iteration + bool mFreezeTranslationX; // Freeze x-translation or not during optimisation + bool mFreezeTranslationY; // Freeze y-translation or not during optimisation + bool mFreezeRadTranslation; // Freeze radiometry translation or not during optimisation + bool mFreezeRadScale; // Freeze radiometry scaling or not during optimisation + double mWeightRadTranslation; // Weight given to radiometry translation if soft freezing is applied (default : negative => not applied) + double mWeightRadScale; // Weight given to radiometry scaling if soft freezing is applied (default : negative => not applied) + int mNumberOfIterGaussFilter; // Number of iterations to be done in Gauss filter algorithm + int mNumberOfEndIterations; // Number of iterations to do while using original image in multi-scale approach + std::string mFolderSaveResult; // Folder name to save results + bool mDisplayLastTranslationValues; // Whether to display the final coordinates of the translated points + bool mDisplayLastRadiometryValues; // Display final values of radiometry unknowns at last iteration of optimisation process + + // == Internal variables ==== + + cPt2di mSzImPre; // size of image + tIm mImPre; // memory representation of the image + tDIm *mDImPre; // memory representation of the image + + cPt2di mSzImPost; // size of image + tIm mImPost; // memory representation of the image + tDIm *mDImPost; // memory representation of the image + + cPt2di mSzImOut; // size of image + tIm mImOut; // memory representation of the image + tDIm *mDImOut; // memory representation of the image + + cPt2di mSzImDiff; // size of image + tIm mImDiff; // memory representation of the image + tDIm *mDImDiff; // memory representation of the image + + cPt2di mSzImDepX; // size of image + tIm mImDepX; // memory representation of the image + tDIm *mDImDepX; // memory representation of the image + + cPt2di mSzImDepY; // size of image + tIm mImDepY; // memory representation of the image + tDIm *mDImDepY; // memory representation of the image + + std::vector mVectorPts; // A vector containing a set of points + cTriangulation2D mDelTri; // A Delaunay triangle + + double mSigmaGaussFilter; // Value of sigma in gauss filter + bool mIsLastIters; // Determines whether optimisation process is at last iters to optimise on original image + + cResolSysNonLinear *mSys; // Non Linear Sys for solving problem + cCalculator *mEqTriDeform; // calculator giving access to values and derivatives + }; +} + +#endif // _TRIANGLEDEFORMATION_H_ diff --git a/MMVII/src/MeshDisplacement/TriangleDeformationRad.cpp b/MMVII/src/MeshDisplacement/TriangleDeformationRad.cpp new file mode 100644 index 0000000000..6936ed57aa --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformationRad.cpp @@ -0,0 +1,405 @@ +#include "MMVII_TplSymbTriangle.h" + +#include "TriangleDeformationRad.h" + +/** + \file TriangleDeformationRad.cpp + + \brief file computing radiometric deformations + between 2 images using triangles. +**/ + +namespace MMVII +{ + + /**********************************************/ + /* */ + /* cTriangleDeformationRad */ + /* */ + /**********************************************/ + + cAppli_cTriangleDeformationRad::cAppli_cTriangleDeformationRad(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) : cAppli_cTriangleDeformationRadiometry(aVArgs, aSpec), + mRandomUniformLawUpperBoundLines(1), + mRandomUniformLawUpperBoundCols(1), + mShow(true), + mGenerateOutputImage(true), + mUseMultiScaleApproach(false), + mWeightRadTranslation(-1), + mWeightRadScale(-1), + mDisplayLastRadiometryValues(false), + mSigmaGaussFilterStep(1), + mNumberOfIterGaussFilter(3), + mNumberOfEndIterations(2), + mSzImPre(cPt2di(1, 1)), + mImPre(mSzImPre), + mDImPre(nullptr), + mSzImPost(cPt2di(1, 1)), + mImPost(mSzImPost), + mDImPost(nullptr), + mSzImOut(cPt2di(1, 1)), + mImOut(mSzImOut), + mDImOut(nullptr), + mVectorPts({cPt2dr(0, 0)}), + mDelTri(mVectorPts), + mSys(nullptr), + mEqRadTri(nullptr) + + { + mEqRadTri = EqDeformTriRad(true, 1); // true means with derivative, 1 is size of buffer + // mEqRadiometryTri->SetDebugEnabled(true); + } + + cAppli_cTriangleDeformationRad::~cAppli_cTriangleDeformationRad() + { + delete mSys; + delete mEqRadTri; + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformationRad::ArgObl(cCollecSpecArg2007 &anArgObl) + { + return anArgObl + << Arg2007(mNamePreImage, "Name of pre-image file.", {{eTA2007::FileImage}, {eTA2007::FileDirProj}}) + << Arg2007(mNamePostImage, "Name of post-image file.", {eTA2007::FileImage}) + << Arg2007(mNumberPointsToGenerate, "Number of points you want to generate for triangulation.") + << Arg2007(mNumberOfScales, "Total number of scales to run in multi-scale approach optimisation process."); + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformationRad::ArgOpt(cCollecSpecArg2007 &anArgOpt) + { + return anArgOpt + << AOpt2007(mRandomUniformLawUpperBoundCols, "RandomUniformLawUpperBoundXAxis", + "Maximum value that the uniform law can draw from on the x-axis.", {eTA2007::HDV}) + << AOpt2007(mRandomUniformLawUpperBoundLines, "RandomUniformLawUpperBoundYAxis", + "Maximum value that the uniform law can draw from for on the y-axis.", {eTA2007::HDV}) + << AOpt2007(mShow, "Show", "Whether to print minimisation results.", {eTA2007::HDV}) + << AOpt2007(mGenerateOutputImage, "GenerateDisplacementImage", + "Whether to generate and save the output image with computed radiometry", {eTA2007::HDV}) + << AOpt2007(mUseMultiScaleApproach, "UseMultiScaleApproach", "Whether to use multi-scale approach or not.", {eTA2007::HDV}) + << AOpt2007(mWeightRadTranslation, "WeightRadiometryTranslation", + "A value to weight radiometry translation for soft freezing of coefficient.", {eTA2007::HDV}) + << AOpt2007(mWeightRadScale, "WeightRadiometryScaling", + "A value to weight radiometry scaling for soft freezing of coefficient.", {eTA2007::HDV}) + << AOpt2007(mDisplayLastRadiometryValues, "DisplayLastRadiometryValues", + "Whether to display or not the last values of radiometry unknowns after optimisation process.", {eTA2007::HDV}) + << AOpt2007(mSigmaGaussFilterStep, "SigmaGaussFilterStep", + "Sigma value to use for Gauss filter in multi-stage approach.", {eTA2007::HDV}) + << AOpt2007(mNumberOfIterGaussFilter, "NumberOfIterationsGaussFilter", + "Number of iterations to run in Gauss filter algorithm.", {eTA2007::HDV}) + << AOpt2007(mNumberOfEndIterations, "NumberOfEndIterations", + "Number of iterations to run on original images in multi-scale approach.", {eTA2007::HDV}); + } + + void cAppli_cTriangleDeformationRad::LoopOverTrianglesAndUpdateParametersRadiometry(const int aIterNumber) + { + //----------- allocate vec of obs : + tDoubleVect aVObs(12, 0.0); // 6 for ImagePre and 6 for ImagePost + + //----------- extract current parameters + tDenseVect aVCur = mSys->CurGlobSol(); // Get current solution. + + tIm aCurPreIm = tIm(mSzImPre); + tDIm *aCurPreDIm = nullptr; + tIm aCurPostIm = tIm(mSzImPost); + tDIm *aCurPostDIm = nullptr; + + mIsLastIters = false; + + if (mUseMultiScaleApproach) + mIsLastIters = cAppli_cTriangleDeformation::ManageDifferentCasesOfEndIterations(aIterNumber, mNumberOfScales, mNumberOfEndIterations, + mIsLastIters, mImPre, mImPost, aCurPreIm, aCurPreDIm, + aCurPostIm, aCurPostDIm); + else + { + cAppli_cTriangleDeformation::LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + cAppli_cTriangleDeformation::LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aCurPreIm = mImPre.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aCurPostIm = mImPost.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + + aCurPreDIm = &aCurPreIm.DIm(); + aCurPostDIm = &aCurPostIm.DIm(); + + mSigmaGaussFilter -= 1; + + const bool aSaveGaussImage = false; + if (aSaveGaussImage) + aCurPreDIm->ToFile("GaussFilteredImPre_iter_" + std::to_string(aIterNumber) + ".tif"); + } + else if (mUseMultiScaleApproach && mIsLastIters) + { + LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + //----------- declaration of indicator of convergence + tREAL8 aSomDif = 0; // sum of difference between untranslated pixel and translated one. + size_t aNbOut = 0; // number of translated pixels out of image + + // Count number of pixels inside triangles for normalisation + size_t aTotalNumberOfInsidePixels = 0; + + // Loop over all triangles to add the observations on each point + for (size_t aTr = 0; aTr < mDelTri.NbFace(); aTr++) + { + const tTri2dr aTri = mDelTri.KthTri(aTr); + const cPt3di aIndicesOfTriKnots = mDelTri.KthFace(aTr); + + const cTriangle2DCompiled aCompTri(aTri); + + std::vector aVectorToFillWithInsidePixels; + aCompTri.PixelsInside(aVectorToFillWithInsidePixels); // get pixels inside triangle + + //----------- index of unknown, finds the associated pixels of current triangle knots + const tIntVect aVecInd = {2 * aIndicesOfTriKnots.x(), 2 * aIndicesOfTriKnots.x() + 1, + 2 * aIndicesOfTriKnots.y(), 2 * aIndicesOfTriKnots.y() + 1, + 2 * aIndicesOfTriKnots.z(), 2 * aIndicesOfTriKnots.z() + 1}; + + const tREAL8 aCurRadTrPointA = aVCur(aVecInd.at(0)); // current translation on radiometry 1st point of triangle + const tREAL8 aCurRadScPointA = aVCur(aVecInd.at(1)); // current scale on radiometry 3rd point of triangle + const tREAL8 aCurRadTrPointB = aVCur(aVecInd.at(2)); // current translation on radiometry 2nd point of triangle + const tREAL8 aCurRadScPointB = aVCur(aVecInd.at(3)); // current scale on radiometry 3rd point of triangle + const tREAL8 aCurRadTrPointC = aVCur(aVecInd.at(4)); // current translation on radiometry 3rd point of triangle + const tREAL8 aCurRadScPointC = aVCur(aVecInd.at(5)); // current scale on radiometry 3rd point of triangle + + // soft constraint radiometric translation + if (mWeightRadTranslation >= 0) + { + const int aSolStep = 2; // adapt step to solution vector configuration + const int aSolStart = 0; + for (size_t aIndCurSol = aSolStart; aIndCurSol < aVecInd.size() - 1; aIndCurSol += aSolStep) + { + const int aIndices = aVecInd.at(aIndCurSol); + mSys->AddEqFixVar(aIndices, aVCur(aIndices), mWeightRadTranslation); + } + } + + // soft constraint radiometric scaling + if (mWeightRadScale >= 0) + { + const int aSolStep = 2; // adapt step to solution vector configuration + const int aSolStart = 1; + for (size_t aIndCurSol = aSolStart; aIndCurSol < aVecInd.size(); aIndCurSol += aSolStep) + { + const int aIndices = aVecInd.at(aIndCurSol); + mSys->AddEqFixVar(aIndices, aVCur(aIndices), mWeightRadScale); + } + } + + const size_t aNumberOfInsidePixels = aVectorToFillWithInsidePixels.size(); + + // Loop over all pixels inside triangle + // size_t is necessary as there can be a lot of pixels in triangles + for (size_t aFilledPixel = 0; aFilledPixel < aNumberOfInsidePixels; aFilledPixel++) + { + const cPtInsideTriangles aPixInsideTriangle = cPtInsideTriangles(aCompTri, aVectorToFillWithInsidePixels, + aFilledPixel, *aCurPreDIm); + // prepare for barycenter translation formula by filling aVObs with different coordinates + FormalInterpBarycenter_SetObs(aVObs, 0, aPixInsideTriangle); + + // radiometry translation of pixel by current radiometry translation of triangle knots + const tREAL8 aRadiometryTranslation = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForTranslationRadiometry(aCurRadTrPointA, + aCurRadTrPointB, + aCurRadTrPointC, + aVObs); + + // radiometry scaling of pixel by current radiometry scaling of triangle knots + const tREAL8 aRadiometryScaling = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForScalingRadiometry(aCurRadScPointA, + aCurRadScPointB, + aCurRadScPointC, + aVObs); + + const cPt2dr aInsideTrianglePoint = aPixInsideTriangle.GetCartesianCoordinates(); + const cPt2di aEastTranslatedPoint = cPt2di(aInsideTrianglePoint.x(), aInsideTrianglePoint.y()) + cPt2di(1, 0); + const cPt2di aSouthTranslatedPoint = cPt2di(aInsideTrianglePoint.x(), aInsideTrianglePoint.y()) + cPt2di(0, 1); + + if (aCurPostDIm->InsideBL(cPt2dr(aEastTranslatedPoint.x(), aEastTranslatedPoint.y()))) // avoid errors + { + if (aCurPostDIm->InsideBL(cPt2dr(aSouthTranslatedPoint.x(), aSouthTranslatedPoint.y()))) // avoid errors + { + // prepare for application of bilinear formula + FormalBilinTri_SetObs(aVObs, TriangleDisplacement_NbObs, aInsideTrianglePoint, *aCurPostDIm); + + // Now add observation + mSys->CalcAndAddObs(mEqRadTri, aVecInd, aVObs); + + // compute indicators + const tREAL8 aBilinearRadiomValue = aRadiometryScaling * aCurPostDIm->GetVBL(aPixInsideTriangle.GetCartesianCoordinates()) + + aRadiometryTranslation; + const tREAL8 aDif = aVObs[5] - aBilinearRadiomValue; // residual : aValueImPre - aBilinearRadiomValue + aSomDif += std::abs(aDif); + } + } + else + aNbOut++; + aTotalNumberOfInsidePixels += aNumberOfInsidePixels; + } + } + + // Update all parameter taking into account previous observation + mSys->SolveUpdateReset(); + + if (mShow) + StdOut() << aIterNumber + 1 << ", " << aSomDif / aTotalNumberOfInsidePixels + << ", " << aNbOut << std::endl; + } + + void cAppli_cTriangleDeformationRad::GenerateOutputImageAndDisplayLastRadiometryValues(const tDenseVect &aVFinalSol, const int aIterNumber) + { + mImOut = tIm(mSzImPre); + mDImOut = &mImOut.DIm(); + mSzImOut = cPt2di(mDImOut->Sz().x(), mDImOut->Sz().y()); + + tIm aLastPostIm = tIm(mSzImPost); + tDIm *aLastPostDIm = nullptr; + cAppli_cTriangleDeformation::LoadImageAndData(aLastPostIm, aLastPostDIm, "post", mImPre, mImPost); + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aLastPostIm = mImPost.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aLastPostDIm = &aLastPostIm.DIm(); + } + + for (const cPt2di &aOutPix : *mDImOut) // Initialise output image + mDImOut->SetV(aOutPix, aLastPostDIm->GetV(aOutPix)); + + for (size_t aLTr = 0; aLTr < mDelTri.NbFace(); aLTr++) + { + const tTri2dr aLastTri = mDelTri.KthTri(aLTr); + const cPt3di aLastIndicesOfTriKnots = mDelTri.KthFace(aLTr); + + const cTriangle2DCompiled aLastCompTri(aLastTri); + + std::vector aLastVectorToFillWithInsidePixels; + aLastCompTri.PixelsInside(aLastVectorToFillWithInsidePixels); + + const tIntVect aLastVecInd = {2 * aLastIndicesOfTriKnots.x(), 2 * aLastIndicesOfTriKnots.x() + 1, + 2 * aLastIndicesOfTriKnots.y(), 2 * aLastIndicesOfTriKnots.y() + 1, + 2 * aLastIndicesOfTriKnots.z(), 2 * aLastIndicesOfTriKnots.z() + 1}; + + const tREAL8 aLastRadTrPointA = aVFinalSol(aLastVecInd.at(0)); // last translation on radiometry 1st point of triangle + const tREAL8 aLastRadScPointA = aVFinalSol(aLastVecInd.at(1)); // last scale on radiometry 1st point of triangle + const tREAL8 aLastRadTrPointB = aVFinalSol(aLastVecInd.at(2)); // last translation on radiometry 2nd point of triangle + const tREAL8 aLastRadScPointB = aVFinalSol(aLastVecInd.at(3)); // last scale on radiometry 2nd point of triangle + const tREAL8 aLastRadTrPointC = aVFinalSol(aLastVecInd.at(4)); // last translation on radiometry 3rd point of triangle + const tREAL8 aLastRadScPointC = aVFinalSol(aLastVecInd.at(5)); // last scale on radiometry 3rd point of triangle + + const size_t aLastNumberOfInsidePixels = aLastVectorToFillWithInsidePixels.size(); + + for (size_t aLastFilledPixel = 0; aLastFilledPixel < aLastNumberOfInsidePixels; aLastFilledPixel++) + { + const cPtInsideTriangles aLastPixInsideTriangle = cPtInsideTriangles(aLastCompTri, aLastVectorToFillWithInsidePixels, + aLastFilledPixel, *aLastPostDIm); + + // radiometry translation of pixel by current radiometry translation + const tREAL8 aLastRadiometryTranslation = cAppli_cTriangleDeformationRadiometry::ApplyLastBarycenterInterpolationFormulaRadiometryTranslation(aLastRadTrPointA, + aLastRadTrPointB, + aLastRadTrPointC, + aLastPixInsideTriangle); + + const tREAL8 aLastRadiometryScaling = cAppli_cTriangleDeformationRadiometry::ApplyLastBarycenterInterpolationFormulaRadiometryScaling(aLastRadScPointA, + aLastRadScPointB, + aLastRadScPointC, + aLastPixInsideTriangle); + + cAppli_cTriangleDeformationRadiometry::FillOutputImage(aLastPixInsideTriangle, aLastRadiometryTranslation, + aLastRadiometryScaling); + } + } + + // save output image with calculated radiometries to image file + mDImOut->ToFile("OutputImage_" + std::to_string(mNumberOfScales) + "_" + std::to_string(mNumberPointsToGenerate) + ".tif"); + + if (mDisplayLastRadiometryValues) + { + for (int aFinalUnk = 0; aFinalUnk < aVFinalSol.DIm().Sz(); aFinalUnk++) + { + StdOut() << aVFinalSol(aFinalUnk) << " "; + if (aFinalUnk % 2 == 1 && aFinalUnk != 0) + StdOut() << std::endl; + } + } + } + + void cAppli_cTriangleDeformationRad::DoOneIterationRadiometry(const int aIterNumber) + { + LoopOverTrianglesAndUpdateParametersRadiometry(aIterNumber); // Iterate over triangles and solve system + + tDenseVect aVFinalSol = mSys->CurGlobSol(); + // Show final translation results and produce displacement maps + if (mUseMultiScaleApproach) + { + if (aIterNumber == (mNumberOfScales + mNumberOfEndIterations - 1)) + GenerateOutputImageAndDisplayLastRadiometryValues(aVFinalSol, aIterNumber); + } + else + { + if (aIterNumber == (mNumberOfScales - 1)) + GenerateOutputImageAndDisplayLastRadiometryValues(aVFinalSol, aIterNumber); + } + } + + //----------------------------------------- + + int cAppli_cTriangleDeformationRad::Exe() + { + // read pre and post images and their sizes + mImPre = tIm::FromFile(mNamePreImage); + mImPost = tIm::FromFile(mNamePostImage); + + mDImPre = &mImPre.DIm(); + mSzImPre = mDImPre->Sz(); + + mDImPost = &mImPost.DIm(); + mSzImPost = mDImPost->Sz(); + + if (mUseMultiScaleApproach) + mSigmaGaussFilter = mNumberOfScales * mSigmaGaussFilterStep; + + if (mShow) + StdOut() << "Iter, " + << "Diff, " + << "NbOut" << std::endl; + + cAppli_cTriangleDeformation::GeneratePointsForDelaunay(mVectorPts, mNumberPointsToGenerate, mRandomUniformLawUpperBoundLines, + mRandomUniformLawUpperBoundCols, mDelTri, mSzImPre); + + cAppli_cTriangleDeformationRadiometry::InitialisationAfterExeRadiometry(mDelTri, mSys); + + if (mUseMultiScaleApproach) + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales + mNumberOfEndIterations; aIterNumber++) + DoOneIterationRadiometry(aIterNumber); + } + else + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales; aIterNumber++) + DoOneIterationRadiometry(aIterNumber); + } + + return EXIT_SUCCESS; + } + + /********************************************/ + // ::MMVII // + /********************************************/ + + tMMVII_UnikPApli Alloc_cTriangleDeformationRad(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) + { + return tMMVII_UnikPApli(new cAppli_cTriangleDeformationRad(aVArgs, aSpec)); + } + + cSpecMMVII_Appli TheSpec_ComputeTriangleDeformationRad( + "ComputeTriangleDeformationRad", + Alloc_cTriangleDeformationRad, + "Compute radiometric deformation between images using triangles and alternative equation", + {eApF::ImProc}, // category + {eApDT::Image}, // input + {eApDT::Image}, // output + __FILE__); + +}; // namespace MMVII \ No newline at end of file diff --git a/MMVII/src/MeshDisplacement/TriangleDeformationRad.h b/MMVII/src/MeshDisplacement/TriangleDeformationRad.h new file mode 100644 index 0000000000..43dd72d74c --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformationRad.h @@ -0,0 +1,92 @@ +#ifndef _TRIANGLEDEFORMATIONRAD_H_ +#define _TRIANGLEDEFORMATIONRAD_H_ + +#include "MMVII_Geom2D.h" +#include "MMVII_PhgrDist.h" +#include "TriangleDeformation.h" +#include "TriangleDeformationRadiometry.h" + +using namespace NS_SymbolicDerivative; + +namespace MMVII +{ + + /******************************************************/ + /* */ + /* cAppli_cTriangleDeformationRad */ + /* */ + /******************************************************/ + + class cAppli_cTriangleDeformationRad : public cAppli_cTriangleDeformationRadiometry + { + public: + typedef cIm2D tIm; + typedef cDataIm2D tDIm; + typedef cTriangle tTri2dr; + typedef cDenseVect tDenseVect; + typedef std::vector tIntVect; + typedef std::vector tDoubleVect; + + cAppli_cTriangleDeformationRad(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec); + ~cAppli_cTriangleDeformationRad(); + + int Exe() override; + cCollecSpecArg2007 &ArgObl(cCollecSpecArg2007 &anArgObl) override; + cCollecSpecArg2007 &ArgOpt(cCollecSpecArg2007 &anArgOpt) override; + + // Iterate of triangles and inside pixels + virtual void DoOneIterationRadiometry(const int aIterNumber); + // Loops over all triangles and solves system to update parameters at end of iteration + virtual void LoopOverTrianglesAndUpdateParametersRadiometry(const int aIterNumber); + // Generate displacement maps of last solution + virtual void GenerateOutputImageAndDisplayLastRadiometryValues(const tDenseVect &aVFinalSol, const int aIterNumber); + + private: + // == Mandatory args ==== + + std::string mNamePreImage; // Name of given pre-image + std::string mNamePostImage; // Name of given post-image + int mNumberPointsToGenerate; // number of generated points + int mNumberOfScales; // number of iterations in optimisation process + + // == Optionnal args ==== + + int mRandomUniformLawUpperBoundLines; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundLines [ + int mRandomUniformLawUpperBoundCols; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundCols [ + bool mShow; // Print result, export image ... + bool mGenerateOutputImage; // Generate image with displaced pixels + bool mUseMultiScaleApproach; // Apply multi-scale approach or not + tREAL8 mWeightRadTranslation; // Weight given to radiometry translation if soft freezing is applied (default : negative => not applied) + tREAL8 mWeightRadScale; // Weight given to radiometry scaling if soft freezing is applied (default : negative => not applied) + bool mDisplayLastRadiometryValues; // Display final values of radiometry unknowns at last iteration of optimisation process + int mSigmaGaussFilterStep; // Decreasing step of sigma value during iterations + int mNumberOfIterGaussFilter; // Number of iterations to be done in Gauss filter algorithm + int mNumberOfEndIterations; // Number of iterations to do while using original image in multi-scale approach + + // == Internal variables ==== + + cPt2di mSzImPre; // size of image + tIm mImPre; // memory representation of the image + tDIm *mDImPre; // memory representation of the image + + cPt2di mSzImPost; // size of image + tIm mImPost; // memory representation of the image + tDIm *mDImPost; // memory representation of the image + + cPt2di mSzImOut; // size of image + tIm mImOut; // memory representation of the image + tDIm *mDImOut; // memory representation of the image + + tREAL8 mSigmaGaussFilter; // Value of sigma in gauss filter + bool mIsLastIters; // Determines whether optimisation process is at last iters to optimise on original image + + std::vector mVectorPts; // A vector containing a set of points + cTriangulation2D mDelTri; // A Delaunay triangle + + cResolSysNonLinear *mSys; // Non Linear Sys for solving problem + cCalculator *mEqRadTri; // calculator giving access to values and derivatives + }; +} + +#endif // _TRIANGLEDEFORMATIONRAD_H_ diff --git a/MMVII/src/MeshDisplacement/TriangleDeformationRadiometry.cpp b/MMVII/src/MeshDisplacement/TriangleDeformationRadiometry.cpp new file mode 100644 index 0000000000..487b64e344 --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformationRadiometry.cpp @@ -0,0 +1,457 @@ +#include "MMVII_TplSymbTriangle.h" + +#include "TriangleDeformationRadiometry.h" + +/** + \file TriangleDeformationRadiometry.cpp + + \brief file computing radiometric deformations + between 2 images using triangles. +**/ + +namespace MMVII +{ + + /**********************************************/ + /* */ + /* cTriangleDeformationRadiometry */ + /* */ + /**********************************************/ + + cAppli_cTriangleDeformationRadiometry::cAppli_cTriangleDeformationRadiometry(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) : cAppli_cTriangleDeformation(aVArgs, aSpec), + mRandomUniformLawUpperBoundLines(1), + mRandomUniformLawUpperBoundCols(1), + mShow(true), + mGenerateOutputImage(true), + mUseMultiScaleApproach(false), + mWeightRadTranslation(-1), + mWeightRadScale(-1), + mDisplayLastRadiometryValues(false), + mSigmaGaussFilterStep(1), + mNumberOfIterGaussFilter(3), + mNumberOfEndIterations(2), + mSzImPre(cPt2di(1, 1)), + mImPre(mSzImPre), + mDImPre(nullptr), + mSzImPost(cPt2di(1, 1)), + mImPost(mSzImPost), + mDImPost(nullptr), + mSzImOut(cPt2di(1, 1)), + mImOut(mSzImOut), + mDImOut(nullptr), + mVectorPts({cPt2dr(0, 0)}), + mDelTri(mVectorPts), + mSys(nullptr), + mEqRadiometryTri(nullptr) + + { + mEqRadiometryTri = EqDeformTriRadiometry(true, 1); // true means with derivative, 1 is size of buffer + // mEqRadiometryTri->SetDebugEnabled(true); + } + + cAppli_cTriangleDeformationRadiometry::~cAppli_cTriangleDeformationRadiometry() + { + delete mSys; + delete mEqRadiometryTri; + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformationRadiometry::ArgObl(cCollecSpecArg2007 &anArgObl) + { + return anArgObl + << Arg2007(mNamePreImage, "Name of pre-image file.", {{eTA2007::FileImage}, {eTA2007::FileDirProj}}) + << Arg2007(mNamePostImage, "Name of post-image file.", {eTA2007::FileImage}) + << Arg2007(mNumberPointsToGenerate, "Number of points you want to generate for triangulation.") + << Arg2007(mNumberOfScales, "Total number of scales to run in multi-scale approach optimisation process."); + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformationRadiometry::ArgOpt(cCollecSpecArg2007 &anArgOpt) + { + return anArgOpt + << AOpt2007(mRandomUniformLawUpperBoundCols, "RandomUniformLawUpperBoundXAxis", + "Maximum value that the uniform law can draw from on the x-axis.", {eTA2007::HDV}) + << AOpt2007(mRandomUniformLawUpperBoundLines, "RandomUniformLawUpperBoundYAxis", + "Maximum value that the uniform law can draw from for on the y-axis.", {eTA2007::HDV}) + << AOpt2007(mShow, "Show", "Whether to print minimisation results.", {eTA2007::HDV}) + << AOpt2007(mGenerateOutputImage, "GenerateDisplacementImage", + "Whether to generate and save the output image with computed radiometry", {eTA2007::HDV}) + << AOpt2007(mUseMultiScaleApproach, "UseMultiScaleApproach", "Whether to use multi-scale approach or not.", {eTA2007::HDV}) + << AOpt2007(mWeightRadTranslation, "WeightRadiometryTranslation", + "A value to weight radiometry translation for soft freezing of coefficient.", {eTA2007::HDV}) + << AOpt2007(mWeightRadScale, "WeightRadiometryScaling", + "A value to weight radiometry scaling for soft freezing of coefficient.", {eTA2007::HDV}) + << AOpt2007(mDisplayLastRadiometryValues, "DisplayLastRadiometryValues", + "Whether to display or not the last values of radiometry unknowns after optimisation process.", {eTA2007::HDV}) + << AOpt2007(mSigmaGaussFilterStep, "SigmaGaussFilterStep", + "Sigma value to use for Gauss filter in multi-stage approach.", {eTA2007::HDV}) + << AOpt2007(mNumberOfIterGaussFilter, "NumberOfIterationsGaussFilter", + "Number of iterations to run in Gauss filter algorithm.", {eTA2007::HDV}) + << AOpt2007(mNumberOfEndIterations, "NumberOfEndIterations", + "Number of iterations to run on original images in multi-scale approach.", {eTA2007::HDV}); + } + + void cAppli_cTriangleDeformationRadiometry::InitialisationAfterExeRadiometry(cTriangulation2D &aDelaunayTri, + cResolSysNonLinear *&aSys) + { + const size_t aNumberPts = 2 * aDelaunayTri.NbPts(); + tDenseVect aVInit(aNumberPts, eModeInitImage::eMIA_Null); // eMIA_V1 + + for (size_t aKtNumber = 0; aKtNumber < aNumberPts; aKtNumber++) + { + if (aKtNumber % 2 == 1) + aVInit(aKtNumber) = 1; + } + + aSys = new cResolSysNonLinear(eModeSSR::eSSR_LsqDense, aVInit); + } + + void cAppli_cTriangleDeformationRadiometry::LoopOverTrianglesAndUpdateParametersRadiometry(const int aIterNumber) + { + //----------- allocate vec of obs : + tDoubleVect aVObs(12, 0.0); // 6 for ImagePre and 6 for ImagePost + + //----------- extract current parameters + tDenseVect aVCur = mSys->CurGlobSol(); // Get current solution. + + tIm aCurPreIm = tIm(mSzImPre); + tDIm *aCurPreDIm = nullptr; + tIm aCurPostIm = tIm(mSzImPost); + tDIm *aCurPostDIm = nullptr; + + mIsLastIters = false; + + if (mUseMultiScaleApproach) + mIsLastIters = cAppli_cTriangleDeformation::ManageDifferentCasesOfEndIterations(aIterNumber, mNumberOfScales, mNumberOfEndIterations, + mIsLastIters, mImPre, mImPost, aCurPreIm, aCurPreDIm, + aCurPostIm, aCurPostDIm); + else + { + cAppli_cTriangleDeformation::LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + cAppli_cTriangleDeformation::LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aCurPreIm = mImPre.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aCurPostIm = mImPost.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + + aCurPreDIm = &aCurPreIm.DIm(); + aCurPostDIm = &aCurPostIm.DIm(); + + mSigmaGaussFilter -= 1; + + const bool aSaveGaussImage = false; + if (aSaveGaussImage) + aCurPreDIm->ToFile("GaussFilteredImPre_iter_" + std::to_string(aIterNumber) + ".tif"); + } + else if (mUseMultiScaleApproach && mIsLastIters) + { + cAppli_cTriangleDeformation::LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + cAppli_cTriangleDeformation::LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + //----------- declaration of indicator of convergence + tREAL8 aSomDif = 0; // sum of difference between untranslated pixel and translated one. + size_t aNbOut = 0; // number of translated pixels out of image + + // Count number of pixels inside triangles for normalisation + size_t aTotalNumberOfInsidePixels = 0; + + // Loop over all triangles to add the observations on each point + for (size_t aTr = 0; aTr < mDelTri.NbFace(); aTr++) + { + const tTri2dr aTri = mDelTri.KthTri(aTr); + const cPt3di aIndicesOfTriKnots = mDelTri.KthFace(aTr); + + const cTriangle2DCompiled aCompTri(aTri); + + std::vector aVectorToFillWithInsidePixels; + aCompTri.PixelsInside(aVectorToFillWithInsidePixels); // get pixels inside triangle + + //----------- index of unknown, finds the associated pixels of current triangle knots + const tIntVect aVecInd = {2 * aIndicesOfTriKnots.x(), 2 * aIndicesOfTriKnots.x() + 1, + 2 * aIndicesOfTriKnots.y(), 2 * aIndicesOfTriKnots.y() + 1, + 2 * aIndicesOfTriKnots.z(), 2 * aIndicesOfTriKnots.z() + 1}; + + const tREAL8 aCurRadTrPointA = aVCur(aVecInd.at(0)); // current translation on radiometry 1st point of triangle + const tREAL8 aCurRadScPointA = aVCur(aVecInd.at(1)); // current scale on radiometry 3rd point of triangle + const tREAL8 aCurRadTrPointB = aVCur(aVecInd.at(2)); // current translation on radiometry 2nd point of triangle + const tREAL8 aCurRadScPointB = aVCur(aVecInd.at(3)); // current scale on radiometry 3rd point of triangle + const tREAL8 aCurRadTrPointC = aVCur(aVecInd.at(4)); // current translation on radiometry 3rd point of triangle + const tREAL8 aCurRadScPointC = aVCur(aVecInd.at(5)); // current scale on radiometry 3rd point of triangle + + // soft constraint radiometric translation + if (mWeightRadTranslation >= 0) + { + const int aSolStep = 2; // adapt step to solution vector configuration + const int aSolStart = 0; + for (size_t aIndCurSol = aSolStart; aIndCurSol < aVecInd.size() - 1; aIndCurSol += aSolStep) + { + const int aIndices = aVecInd.at(aIndCurSol); + mSys->AddEqFixVar(aIndices, aVCur(aIndices), mWeightRadTranslation); + } + } + + // soft constraint radiometric scaling + if (mWeightRadScale >= 0) + { + const int aSolStep = 2; // adapt step to solution vector configuration + const int aSolStart = 1; + for (size_t aIndCurSol = aSolStart; aIndCurSol < aVecInd.size(); aIndCurSol += aSolStep) + { + const int aIndices = aVecInd.at(aIndCurSol); + mSys->AddEqFixVar(aIndices, aVCur(aIndices), mWeightRadScale); + } + } + + const size_t aNumberOfInsidePixels = aVectorToFillWithInsidePixels.size(); + + // Loop over all pixels inside triangle + // size_t is necessary as there can be a lot of pixels in triangles + for (size_t aFilledPixel = 0; aFilledPixel < aNumberOfInsidePixels; aFilledPixel++) + { + const cPtInsideTriangles aPixInsideTriangle = cPtInsideTriangles(aCompTri, aVectorToFillWithInsidePixels, + aFilledPixel, *aCurPreDIm); + // prepare for barycenter translation formula by filling aVObs with different coordinates + FormalInterpBarycenter_SetObs(aVObs, 0, aPixInsideTriangle); + + // radiometry translation of pixel by current radiometry translation of triangle knots + const tREAL8 aRadiometryTranslation = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForTranslationRadiometry(aCurRadTrPointA, + aCurRadTrPointB, + aCurRadTrPointC, + aVObs); + + // radiometry scaling of pixel by current radiometry scaling of triangle knots + const tREAL8 aRadiometryScaling = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForScalingRadiometry(aCurRadScPointA, + aCurRadScPointB, + aCurRadScPointC, + aVObs); + + const cPt2dr aInsideTrianglePoint = aPixInsideTriangle.GetCartesianCoordinates(); + const cPt2di aEastTranslatedPoint = cPt2di(aInsideTrianglePoint.x(), aInsideTrianglePoint.y()) + cPt2di(1, 0); + const cPt2di aSouthTranslatedPoint = cPt2di(aInsideTrianglePoint.x(), aInsideTrianglePoint.y()) + cPt2di(0, 1); + + if (aCurPostDIm->InsideBL(cPt2dr(aEastTranslatedPoint.x(), aEastTranslatedPoint.y()))) // avoid errors + { + if (aCurPostDIm->InsideBL(cPt2dr(aSouthTranslatedPoint.x(), aSouthTranslatedPoint.y()))) // avoid errors + { + // prepare for application of bilinear formula + FormalBilinTri_SetObs(aVObs, TriangleDisplacement_NbObs, aInsideTrianglePoint, *aCurPostDIm); + + // Now add observation + mSys->CalcAndAddObs(mEqRadiometryTri, aVecInd, aVObs); + + // compute indicators + const tREAL8 aRadiomValueImPre = aRadiometryScaling * aVObs[5] + aRadiometryTranslation; + const tREAL8 aDif = aRadiomValueImPre - aCurPostDIm->GetVBL(aInsideTrianglePoint); // residual + aSomDif += std::abs(aDif); + } + } + else + aNbOut++; + aTotalNumberOfInsidePixels += aNumberOfInsidePixels; + } + } + + // Update all parameter taking into account previous observation + mSys->SolveUpdateReset(); + + if (mShow) + StdOut() << aIterNumber + 1 << ", " << aSomDif / aTotalNumberOfInsidePixels + << ", " << aNbOut << std::endl; + } + + void cAppli_cTriangleDeformationRadiometry::FillOutputImage(const cPtInsideTriangles &aLastPixInsideTriangle, + const tREAL8 aLastRadiometryTranslation, + const tREAL8 aLastRadiometryScaling) + { + const tREAL8 aLastXCoordinate = aLastPixInsideTriangle.GetCartesianCoordinates().x(); + const tREAL8 aLastYCoordinate = aLastPixInsideTriangle.GetCartesianCoordinates().y(); + + const cPt2di aLastCoordinate = cPt2di(aLastXCoordinate, aLastYCoordinate); + + const tREAL8 aLastRadiometryValue = aLastRadiometryScaling * aLastPixInsideTriangle.GetPixelValue() + + aLastRadiometryTranslation; + + // Build image with radiometric intensities + mDImOut->SetV(aLastCoordinate, aLastRadiometryValue); + } + + tREAL8 cAppli_cTriangleDeformationRadiometry::ApplyLastBarycenterInterpolationFormulaRadiometryTranslation(const tREAL8 aLastRadTranslationPointA, + const tREAL8 aLastRadTranslationPointB, + const tREAL8 aLastRadTranslationPointC, + const cPtInsideTriangles &aLastPixInsideTriangle) + { + const cPt3dr BarycenterCoordinateOfPoint = aLastPixInsideTriangle.GetBarycenterCoordinates(); + const tREAL8 aCurentRadTranslation = BarycenterCoordinateOfPoint.x() * aLastRadTranslationPointA + BarycenterCoordinateOfPoint.y() * aLastRadTranslationPointB + + BarycenterCoordinateOfPoint.z() * aLastRadTranslationPointC; + return aCurentRadTranslation; + } + + tREAL8 cAppli_cTriangleDeformationRadiometry::ApplyLastBarycenterInterpolationFormulaRadiometryScaling(const tREAL8 aLastRadScalingPointA, + const tREAL8 aLastRadScalingPointB, + const tREAL8 aLastRadScalingPointC, + const cPtInsideTriangles &aLastPixInsideTriangle) + { + const cPt3dr BarycenterCoordinateOfPoint = aLastPixInsideTriangle.GetBarycenterCoordinates(); + const tREAL8 aCurrentRadScaling = BarycenterCoordinateOfPoint.x() * aLastRadScalingPointA + BarycenterCoordinateOfPoint.y() * aLastRadScalingPointB + + BarycenterCoordinateOfPoint.z() * aLastRadScalingPointC; + return aCurrentRadScaling; + } + + void cAppli_cTriangleDeformationRadiometry::GenerateOutputImageAndDisplayLastRadiometryValues(const tDenseVect &aVFinalSol, const int aIterNumber) + { + mImOut = tIm(mSzImPre); + mDImOut = &mImOut.DIm(); + mSzImOut = cPt2di(mDImOut->Sz().x(), mDImOut->Sz().y()); + + tIm aLastPreIm = tIm(mSzImPre); + tDIm *aLastPreDIm = nullptr; + cAppli_cTriangleDeformation::LoadImageAndData(aLastPreIm, aLastPreDIm, "pre", mImPre, mImPost); + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aLastPreIm = mImPre.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aLastPreDIm = &aLastPreIm.DIm(); + } + + for (const cPt2di &aOutPix : *mDImOut) // Initialise output image + mDImOut->SetV(aOutPix, aLastPreDIm->GetV(aOutPix)); + + for (size_t aLTr = 0; aLTr < mDelTri.NbFace(); aLTr++) + { + const tTri2dr aLastTri = mDelTri.KthTri(aLTr); + const cPt3di aLastIndicesOfTriKnots = mDelTri.KthFace(aLTr); + + const cTriangle2DCompiled aLastCompTri(aLastTri); + + std::vector aLastVectorToFillWithInsidePixels; + aLastCompTri.PixelsInside(aLastVectorToFillWithInsidePixels); + + const tIntVect aLastVecInd = {2 * aLastIndicesOfTriKnots.x(), 2 * aLastIndicesOfTriKnots.x() + 1, + 2 * aLastIndicesOfTriKnots.y(), 2 * aLastIndicesOfTriKnots.y() + 1, + 2 * aLastIndicesOfTriKnots.z(), 2 * aLastIndicesOfTriKnots.z() + 1}; + + const tREAL8 aLastRadTrPointA = aVFinalSol(aLastVecInd.at(0)); // last translation on radiometry 1st point of triangle + const tREAL8 aLastRadScPointA = aVFinalSol(aLastVecInd.at(1)); // last scale on radiometry 1st point of triangle + const tREAL8 aLastRadTrPointB = aVFinalSol(aLastVecInd.at(2)); // last translation on radiometry 2nd point of triangle + const tREAL8 aLastRadScPointB = aVFinalSol(aLastVecInd.at(3)); // last scale on radiometry 2nd point of triangle + const tREAL8 aLastRadTrPointC = aVFinalSol(aLastVecInd.at(4)); // last translation on radiometry 3rd point of triangle + const tREAL8 aLastRadScPointC = aVFinalSol(aLastVecInd.at(5)); // last scale on radiometry 3rd point of triangle + + const size_t aLastNumberOfInsidePixels = aLastVectorToFillWithInsidePixels.size(); + + for (size_t aLastFilledPixel = 0; aLastFilledPixel < aLastNumberOfInsidePixels; aLastFilledPixel++) + { + const cPtInsideTriangles aLastPixInsideTriangle = cPtInsideTriangles(aLastCompTri, aLastVectorToFillWithInsidePixels, + aLastFilledPixel, *aLastPreDIm); + + // radiometry translation of pixel by current radiometry translation + const tREAL8 aLastRadiometryTranslation = cAppli_cTriangleDeformationRadiometry::ApplyLastBarycenterInterpolationFormulaRadiometryTranslation(aLastRadTrPointA, + aLastRadTrPointB, + aLastRadTrPointC, + aLastPixInsideTriangle); + + const tREAL8 aLastRadiometryScaling = cAppli_cTriangleDeformationRadiometry::ApplyLastBarycenterInterpolationFormulaRadiometryScaling(aLastRadScPointA, + aLastRadScPointB, + aLastRadScPointC, + aLastPixInsideTriangle); + + FillOutputImage(aLastPixInsideTriangle, aLastRadiometryTranslation, + aLastRadiometryScaling); + } + } + + // save output image with calculated radiometries to image file + mDImOut->ToFile("OutputImage_" + std::to_string(mNumberOfScales) + "_" + std::to_string(mNumberPointsToGenerate) + ".tif"); + + if (mDisplayLastRadiometryValues) + { + for (int aFinalUnk = 0; aFinalUnk < aVFinalSol.DIm().Sz(); aFinalUnk++) + { + StdOut() << aVFinalSol(aFinalUnk) << " "; + if (aFinalUnk % 2 == 1 && aFinalUnk != 0) + StdOut() << std::endl; + } + } + } + + void cAppli_cTriangleDeformationRadiometry::DoOneIterationRadiometry(const int aIterNumber) + { + LoopOverTrianglesAndUpdateParametersRadiometry(aIterNumber); // Iterate over triangles and solve system + + tDenseVect aVFinalSol = mSys->CurGlobSol(); + // Show final translation results and produce displacement maps + if (mUseMultiScaleApproach) + { + if (aIterNumber == (mNumberOfScales + mNumberOfEndIterations - 1)) + GenerateOutputImageAndDisplayLastRadiometryValues(aVFinalSol, aIterNumber); + } + else + { + if (aIterNumber == (mNumberOfScales - 1)) + GenerateOutputImageAndDisplayLastRadiometryValues(aVFinalSol, aIterNumber); + } + } + + //----------------------------------------- + + int cAppli_cTriangleDeformationRadiometry::Exe() + { + // read pre and post images and their sizes + mImPre = tIm::FromFile(mNamePreImage); + mImPost = tIm::FromFile(mNamePostImage); + + mDImPre = &mImPre.DIm(); + mSzImPre = mDImPre->Sz(); + + mDImPost = &mImPost.DIm(); + mSzImPost = mDImPost->Sz(); + + if (mUseMultiScaleApproach) + mSigmaGaussFilter = mNumberOfScales * mSigmaGaussFilterStep; + + if (mShow) + StdOut() << "Iter, " + << "Diff, " + << "NbOut" << std::endl; + + cAppli_cTriangleDeformation::GeneratePointsForDelaunay(mVectorPts, mNumberPointsToGenerate, mRandomUniformLawUpperBoundLines, + mRandomUniformLawUpperBoundCols, mDelTri, mSzImPre); + + cAppli_cTriangleDeformationRadiometry::InitialisationAfterExeRadiometry(mDelTri, mSys); + + if (mUseMultiScaleApproach) + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales + mNumberOfEndIterations; aIterNumber++) + DoOneIterationRadiometry(aIterNumber); + } + else + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales; aIterNumber++) + DoOneIterationRadiometry(aIterNumber); + } + + return EXIT_SUCCESS; + } + + /********************************************/ + // ::MMVII // + /********************************************/ + + tMMVII_UnikPApli Alloc_cTriangleDeformationRadiometry(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) + { + return tMMVII_UnikPApli(new cAppli_cTriangleDeformationRadiometry(aVArgs, aSpec)); + } + + cSpecMMVII_Appli TheSpec_ComputeTriangleDeformationRadiometry( + "ComputeTriangleDeformationRadiometry", + Alloc_cTriangleDeformationRadiometry, + "Compute radiometric deformation between images using triangles", + {eApF::ImProc}, // category + {eApDT::Image}, // input + {eApDT::Image}, // output + __FILE__); + +}; // namespace MMVII \ No newline at end of file diff --git a/MMVII/src/MeshDisplacement/TriangleDeformationRadiometry.h b/MMVII/src/MeshDisplacement/TriangleDeformationRadiometry.h new file mode 100644 index 0000000000..6af4d4f3a9 --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformationRadiometry.h @@ -0,0 +1,109 @@ +#ifndef _TRIANGLEDEFORMATIONRADIOMETRY_H_ +#define _TRIANGLEDEFORMATIONRADIOMETRY_H_ + +#include "MMVII_Geom2D.h" +#include "MMVII_PhgrDist.h" +#include "TriangleDeformation.h" + +using namespace NS_SymbolicDerivative; + +namespace MMVII +{ + + /**************************************************/ + /* */ + /* cAppli_cTriangleDeformationRadiometry */ + /* */ + /**************************************************/ + + class cAppli_cTriangleDeformationRadiometry : public cAppli_cTriangleDeformation + { + public: + typedef cIm2D tIm; + typedef cDataIm2D tDIm; + typedef cTriangle tTri2dr; + typedef cDenseVect tDenseVect; + typedef std::vector tIntVect; + typedef std::vector tDoubleVect; + + cAppli_cTriangleDeformationRadiometry(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec); + ~cAppli_cTriangleDeformationRadiometry(); + + int Exe() override; + cCollecSpecArg2007 &ArgObl(cCollecSpecArg2007 &anArgObl) override; + cCollecSpecArg2007 &ArgOpt(cCollecSpecArg2007 &anArgOpt) override; + + // Iterate of triangles and inside pixels + virtual void DoOneIterationRadiometry(const int aIterNumber); + // Loops over all triangles and solves system to update parameters at end of iteration + virtual void LoopOverTrianglesAndUpdateParametersRadiometry(const int aIterNumber); + // Generate displacement maps of last solution + virtual void GenerateOutputImageAndDisplayLastRadiometryValues(const tDenseVect &aVFinalSol, const int aIterNumber); + // Initialise problem after user has input information + void InitialisationAfterExeRadiometry(cTriangulation2D &aDelaunayTri, + cResolSysNonLinear *&aSys); + // Fill displacement maps and output image + void FillOutputImage(const cPtInsideTriangles &aLastPixInsideTriangle, + const tREAL8 aLastRadiometryTranslation, + const tREAL8 aLastRadiometryScaling); + + // Apply barycentric translation formula to last radiometric translation values + tREAL8 ApplyLastBarycenterInterpolationFormulaRadiometryTranslation(const tREAL8 aLastRadTranslationPointA, + const tREAL8 aLastRadTranslationPointB, + const tREAL8 aLastRadTranslationPointC, + const cPtInsideTriangles &aLastPixInsideTriangle); + // Apply barycentric translation formula to last radiometric scaling values + tREAL8 ApplyLastBarycenterInterpolationFormulaRadiometryScaling(const tREAL8 aLastRadScalingPointA, + const tREAL8 aLastRadScalingPointB, + const tREAL8 aLastRadScalingPointC, + const cPtInsideTriangles &aLastPixInsideTriangle); + + private: + // == Mandatory args ==== + + std::string mNamePreImage; // Name of given pre-image + std::string mNamePostImage; // Name of given post-image + int mNumberPointsToGenerate; // number of generated points + int mNumberOfScales; // number of iterations in optimisation process + + // == Optionnal args ==== + + int mRandomUniformLawUpperBoundLines; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundLines [ + int mRandomUniformLawUpperBoundCols; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundCols [ + bool mShow; // Print result, export image ... + bool mGenerateOutputImage; // Generate image with displaced pixels + bool mUseMultiScaleApproach; // Apply multi-scale approach or not + tREAL8 mWeightRadTranslation; // Weight given to radiometry translation if soft freezing is applied (default : negative => not applied) + tREAL8 mWeightRadScale; // Weight given to radiometry scaling if soft freezing is applied (default : negative => not applied) + bool mDisplayLastRadiometryValues; // Display final values of radiometry unknowns at last iteration of optimisation process + int mSigmaGaussFilterStep; // Decreasing step of sigma value during iterations + int mNumberOfIterGaussFilter; // Number of iterations to be done in Gauss filter algorithm + int mNumberOfEndIterations; // Number of iterations to do while using original image in multi-scale approach + + // == Internal variables ==== + + cPt2di mSzImPre; // size of image + tIm mImPre; // memory representation of the image + tDIm *mDImPre; // memory representation of the image + + cPt2di mSzImPost; // size of image + tIm mImPost; // memory representation of the image + tDIm *mDImPost; // memory representation of the image + + cPt2di mSzImOut; // size of image + tIm mImOut; // memory representation of the image + tDIm *mDImOut; // memory representation of the image + + tREAL8 mSigmaGaussFilter; // Value of sigma in gauss filter + bool mIsLastIters; // Determines whether optimisation process is at last iters to optimise on original image + + std::vector mVectorPts; // A vector containing a set of points + cTriangulation2D mDelTri; // A Delaunay triangle + + cResolSysNonLinear *mSys; // Non Linear Sys for solving problem + cCalculator *mEqRadiometryTri; // calculator giving access to values and derivatives + }; +} + +#endif // _TRIANGLEDEFORMATIONRADIOMETRY_H_ diff --git a/MMVII/src/MeshDisplacement/TriangleDeformationTrRad.cpp b/MMVII/src/MeshDisplacement/TriangleDeformationTrRad.cpp new file mode 100644 index 0000000000..acec5284f8 --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformationTrRad.cpp @@ -0,0 +1,583 @@ +#include "cMMVII_Appli.h" + +#include "MMVII_TplSymbTriangle.h" + +#include "TriangleDeformationTrRad.h" + +/** + \file TriangleDeformationTrRad.cpp + + \brief file for computing 2D deformations (translation and radiometry) + between 2 images thanks to triangles. +**/ + +namespace MMVII +{ + /******************************************/ + /* */ + /* cAppli_cTriangleDeformationTrRad */ + /* */ + /******************************************/ + + cAppli_cTriangleDeformationTrRad::cAppli_cTriangleDeformationTrRad(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) : cAppli_cTriangleDeformation(aVArgs, aSpec), + mRandomUniformLawUpperBoundLines(1), + mRandomUniformLawUpperBoundCols(1), + mShow(true), + mUseMultiScaleApproach(true), + mSigmaGaussFilterStep(1), + mGenerateDisplacementImage(true), + mFreezeTranslationX(false), + mFreezeTranslationY(false), + mFreezeRadTranslation(false), + mFreezeRadScale(false), + mWeightRadTranslation(-1), + mWeightRadScale(-1), + mNumberOfIterGaussFilter(3), + mNumberOfEndIterations(2), + mDisplayLastTranslationValues(false), + mDisplayLastRadiometryValues(false), + mSzImPre(cPt2di(1, 1)), + mImPre(mSzImPre), + mDImPre(nullptr), + mSzImPost(cPt2di(1, 1)), + mImPost(mSzImPost), + mDImPost(nullptr), + mSzImOut(cPt2di(1, 1)), + mImOut(mSzImOut), + mDImOut(nullptr), + mSzImDepX(cPt2di(1, 1)), + mImDepX(mSzImDepX), + mDImDepX(nullptr), + mSzImDepY(cPt2di(1, 1)), + mImDepY(mSzImDepY), + mDImDepY(nullptr), + mVectorPts({cPt2dr(0, 0)}), + mDelTri(mVectorPts), + mSys(nullptr), + mEqTriDeformTrRad(nullptr) + { + mEqTriDeformTrRad = EqDeformTriTrRad(true, 1); // true means with derivative, 1 is size of buffer + } + + cAppli_cTriangleDeformationTrRad::~cAppli_cTriangleDeformationTrRad() + { + delete mSys; + delete mEqTriDeformTrRad; + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformationTrRad::ArgObl(cCollecSpecArg2007 &anArgObl) + { + return anArgObl + << Arg2007(mNamePreImage, "Name of pre-image file.", {{eTA2007::FileImage}, {eTA2007::FileDirProj}}) + << Arg2007(mNamePostImage, "Name of post-image file.", {eTA2007::FileImage}) + << Arg2007(mNumberPointsToGenerate, "Number of points you want to generate for triangulation.") + << Arg2007(mNumberOfScales, "Total number of scales to run in multi-scale approach optimisation process."); + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformationTrRad::ArgOpt(cCollecSpecArg2007 &anArgOpt) + { + return anArgOpt + << AOpt2007(mRandomUniformLawUpperBoundCols, "RandomUniformLawUpperBoundXAxis", + "Maximum value that the uniform law can draw from on the x-axis.", {eTA2007::HDV}) + << AOpt2007(mRandomUniformLawUpperBoundLines, "RandomUniformLawUpperBoundYAxis", + "Maximum value that the uniform law can draw from for on the y-axis.", {eTA2007::HDV}) + << AOpt2007(mShow, "Show", "Whether to print minimisation results.", {eTA2007::HDV}) + << AOpt2007(mUseMultiScaleApproach, "UseMultiScaleApproach", "Whether to use multi-scale approach or not.", {eTA2007::HDV}) + << AOpt2007(mSigmaGaussFilterStep, "SigmaGaussFilterStep", "Sigma value to use for Gauss filter in multi-stage approach.", {eTA2007::HDV}) + << AOpt2007(mGenerateDisplacementImage, "GenerateDisplacementImage", + "Whether to generate and save an image having been translated.", {eTA2007::HDV}) + << AOpt2007(mFreezeTranslationX, "FreezeXTranslation", + "Whether to freeze or not x-translation to certain value during computation.", {eTA2007::HDV}) + << AOpt2007(mFreezeTranslationY, "FreezeYTranslation", + "Whether to freeze or not y-translation to certain value during computation.", {eTA2007::HDV}) + << AOpt2007(mFreezeRadTranslation, "FreezeRadTranslation", + "Whether to freeze radiometry translation factor in computation or not.", {eTA2007::HDV}) + << AOpt2007(mFreezeRadScale, "FreezeRadScaling", + "Whether to freeze radiometry scaling factor in computation or not.", {eTA2007::HDV}) + << AOpt2007(mWeightRadTranslation, "WeightRadiometryTranslation", + "A value to weight radiometry translation for soft freezing of coefficient.", {eTA2007::HDV}) + << AOpt2007(mWeightRadScale, "WeightRadiometryScaling", + "A value to weight radiometry scaling for soft freezing of coefficient.", {eTA2007::HDV}) + << AOpt2007(mNumberOfIterGaussFilter, "NumberOfIterationsGaussFilter", + "Number of iterations to run in Gauss filter algorithm.", {eTA2007::HDV}) + << AOpt2007(mNumberOfEndIterations, "NumberOfEndIterations", + "Number of iterations to run on original images in multi-scale approach.", {eTA2007::HDV}) + << AOpt2007(mDisplayLastTranslationValues, "DisplayLastTranslationsValues", + "Whether to display the final values of unknowns linked to point translation.", {eTA2007::HDV}) + << AOpt2007(mDisplayLastRadiometryValues, "DisplayLastRadiometryValues", + "Whether to display or not the last values of radiometry unknowns after optimisation process.", {eTA2007::HDV}); + } + + void cAppli_cTriangleDeformationTrRad::LoopOverTrianglesAndUpdateParameters(const int aIterNumber) + { + //----------- allocate vec of obs : + tDoubleVect aVObs(12, 0.0); // 6 for ImagePre and 6 for ImagePost + + //----------- extract current parameters + tDenseVect aVCur = mSys->CurGlobSol(); // Get current solution. + + tIm aCurPreIm = tIm(mSzImPre); + tDIm *aCurPreDIm = nullptr; + tIm aCurPostIm = tIm(mSzImPost); + tDIm *aCurPostDIm = nullptr; + + mIsLastIters = false; + + if (mUseMultiScaleApproach) + mIsLastIters = cAppli_cTriangleDeformation::ManageDifferentCasesOfEndIterations(aIterNumber, mNumberOfScales, mNumberOfEndIterations, + mIsLastIters, mImPre, mImPost, aCurPreIm, aCurPreDIm, + aCurPostIm, aCurPostDIm); + else + { + cAppli_cTriangleDeformation::LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + cAppli_cTriangleDeformation::LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aCurPreIm = mImPre.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aCurPostIm = mImPost.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + + aCurPreDIm = &aCurPreIm.DIm(); + aCurPostDIm = &aCurPostIm.DIm(); + + mSigmaGaussFilter -= 1; + + const bool aSaveGaussImage = false; + if (aSaveGaussImage) + aCurPreDIm->ToFile("GaussFilteredImPre_iter_" + std::to_string(aIterNumber) + ".tif"); + } + else if (mUseMultiScaleApproach && mIsLastIters) + { + cAppli_cTriangleDeformation::LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + cAppli_cTriangleDeformation::LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + //----------- declaration of indicator of convergence + tREAL8 aSomDif = 0; // sum of difference between untranslated pixel and translated one. + size_t aNbOut = 0; // number of translated pixels out of image + + // Count number of pixels inside triangles for normalisation + size_t aTotalNumberOfInsidePixels = 0; + + // hard constraint : freeze radiometric coefficients + if (mFreezeRadTranslation || mFreezeRadScale || + mFreezeTranslationX || mFreezeTranslationY) + { + for (size_t aTr = 0; aTr < mDelTri.NbFace(); aTr++) + { + const cPt3di aIndicesOfTriKnots = mDelTri.KthFace(aTr); + + const tIntVect aVecInd = {4 * aIndicesOfTriKnots.x(), 4 * aIndicesOfTriKnots.x() + 1, + 4 * aIndicesOfTriKnots.x() + 2, 4 * aIndicesOfTriKnots.x() + 3, + 4 * aIndicesOfTriKnots.y(), 4 * aIndicesOfTriKnots.y() + 1, + 4 * aIndicesOfTriKnots.y() + 2, 4 * aIndicesOfTriKnots.y() + 3, + 4 * aIndicesOfTriKnots.z(), 4 * aIndicesOfTriKnots.z() + 1, + 4 * aIndicesOfTriKnots.z() + 2, 4 * aIndicesOfTriKnots.z() + 3}; + + if (mFreezeTranslationX) + { + mSys->SetFrozenVar(aVecInd.at(2), aVCur(0)); + mSys->SetFrozenVar(aVecInd.at(6), aVCur(4)); + mSys->SetFrozenVar(aVecInd.at(10), aVCur(8)); + } + if (mFreezeTranslationY) + { + mSys->SetFrozenVar(aVecInd.at(2), aVCur(1)); + mSys->SetFrozenVar(aVecInd.at(6), aVCur(5)); + mSys->SetFrozenVar(aVecInd.at(10), aVCur(9)); + } + if (mFreezeRadTranslation) + { + mSys->SetFrozenVar(aVecInd.at(2), aVCur(2)); + mSys->SetFrozenVar(aVecInd.at(6), aVCur(6)); + mSys->SetFrozenVar(aVecInd.at(10), aVCur(10)); + } + if (mFreezeRadScale) + { + mSys->SetFrozenVar(aVecInd.at(3), aVCur(3)); + mSys->SetFrozenVar(aVecInd.at(7), aVCur(7)); + mSys->SetFrozenVar(aVecInd.at(11), aVCur(11)); + } + } + } + + // Loop over all triangles to add the observations on each point + for (size_t aTr = 0; aTr < mDelTri.NbFace(); aTr++) + { + const tTri2dr aTri = mDelTri.KthTri(aTr); + const cPt3di aIndicesOfTriKnots = mDelTri.KthFace(aTr); + + const cTriangle2DCompiled aCompTri(aTri); + + std::vector aVectorToFillWithInsidePixels; + aCompTri.PixelsInside(aVectorToFillWithInsidePixels); // get pixels inside triangle + + //----------- index of unknown, finds the associated pixels of current triangle + const tIntVect aVecInd = {4 * aIndicesOfTriKnots.x(), 4 * aIndicesOfTriKnots.x() + 1, + 4 * aIndicesOfTriKnots.x() + 2, 4 * aIndicesOfTriKnots.x() + 3, + 4 * aIndicesOfTriKnots.y(), 4 * aIndicesOfTriKnots.y() + 1, + 4 * aIndicesOfTriKnots.y() + 2, 4 * aIndicesOfTriKnots.y() + 3, + 4 * aIndicesOfTriKnots.z(), 4 * aIndicesOfTriKnots.z() + 1, + 4 * aIndicesOfTriKnots.z() + 2, 4 * aIndicesOfTriKnots.z() + 3}; + + const cPt2dr aCurTrPointA = cPt2dr(aVCur(aVecInd.at(0)), + aVCur(aVecInd.at(1))); // current translation 1st point of triangle + const cPt2dr aCurTrPointB = cPt2dr(aVCur(aVecInd.at(4)), + aVCur(aVecInd.at(5))); // current translation 2nd point of triangle + const cPt2dr aCurTrPointC = cPt2dr(aVCur(aVecInd.at(8)), + aVCur(aVecInd.at(9))); // current translation 3rd point of triangle + + const tREAL8 aCurRadTrPointA = aVCur(aVecInd.at(2)); // current translation on radiometry 1st point of triangle + const tREAL8 aCurRadScPointA = aVCur(aVecInd.at(3)); // current scale on radiometry 3rd point of triangle + const tREAL8 aCurRadTrPointB = aVCur(aVecInd.at(6)); // current translation on radiometry 2nd point of triangle + const tREAL8 aCurRadScPointB = aVCur(aVecInd.at(7)); // current scale on radiometry 3rd point of triangle + const tREAL8 aCurRadTrPointC = aVCur(aVecInd.at(10)); // current translation on radiometry 3rd point of triangle + const tREAL8 aCurRadScPointC = aVCur(aVecInd.at(11)); // current scale on radiometry 3rd point of triangle + + // soft constraint radiometric translation + if (!mFreezeRadTranslation) + { + if (mWeightRadTranslation >= 0) + { + const int aSolStep = 4; // adapt step to solution vector configuration + const int aSolStart = 2; + for (size_t aIndCurSol = aSolStart; aIndCurSol < aVecInd.size() - 1; aIndCurSol += aSolStep) + { + const int aIndices = aVecInd.at(aIndCurSol); + mSys->AddEqFixVar(aIndices, aVCur(aIndices), mWeightRadTranslation); + } + } + } + + // soft constraint radiometric scaling + if (!mFreezeRadScale) + { + if (mWeightRadScale >= 0) + { + const int aSolStep = 4; // adapt step to solution vector configuration + const int aSolStart = 3; + for (size_t aIndCurSol = aSolStart; aIndCurSol < aVecInd.size(); aIndCurSol += aSolStep) + { + const int aIndices = aVecInd.at(aIndCurSol); + mSys->AddEqFixVar(aIndices, aVCur(aIndices), mWeightRadScale); + } + } + } + + const size_t aNumberOfInsidePixels = aVectorToFillWithInsidePixels.size(); + + // Loop over all pixels inside triangle + // size_t is necessary as there can be a lot of pixels in triangles + for (size_t aFilledPixel = 0; aFilledPixel < aNumberOfInsidePixels; aFilledPixel++) + { + const cPtInsideTriangles aPixInsideTriangle = cPtInsideTriangles(aCompTri, aVectorToFillWithInsidePixels, + aFilledPixel, *aCurPreDIm); + // prepare for barycenter translation formula by filling aVObs with different coordinates + FormalInterpBarycenter_SetObs(aVObs, 0, aPixInsideTriangle); + + // image of a point in triangle by current translation + const cPt2dr aTranslatedFilledPoint = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaToFilledPixel(aCurTrPointA, aCurTrPointB, + aCurTrPointC, aVObs); + // radiometry translation of pixel by current radiometry translation of triangle knots + const tREAL8 aRadiometryTranslation = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForTranslationRadiometry(aCurRadTrPointA, + aCurRadTrPointB, + aCurRadTrPointC, + aVObs); + // radiometry translation of pixel by current radiometry scaling of triangle knots + const tREAL8 aRadiometryScaling = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForScalingRadiometry(aCurRadScPointA, + aCurRadScPointB, + aCurRadScPointC, + aVObs); + + if (aCurPostDIm->InsideBL(aTranslatedFilledPoint)) // avoid errors + { + // prepare for application of bilinear formula + FormalBilinTri_SetObs(aVObs, TriangleDisplacement_NbObs, aTranslatedFilledPoint, *aCurPostDIm); + + // Now add observation + mSys->CalcAndAddObs(mEqTriDeformTrRad, aVecInd, aVObs); + + // compute indicators + const tREAL8 aBilinearRadiomValue = aRadiometryScaling * aCurPostDIm->GetVBL(aTranslatedFilledPoint) + aRadiometryTranslation; + const tREAL8 aDif = aVObs[5] - aBilinearRadiomValue; // residual : aValueImPre - aBilinearRadiomValue + aSomDif += std::abs(aDif); + } + else + aNbOut++; // Count number of pixels translated outside post image + + aTotalNumberOfInsidePixels += aNumberOfInsidePixels; + } + } + + if (mUseMultiScaleApproach && !mIsLastIters && aIterNumber != 0) + { + const bool aGenerateIntermediateMaps = false; + if (aGenerateIntermediateMaps) + GenerateDisplacementMapsAndOutputImages(aVCur, aIterNumber); + } + + // Update all parameter taking into account previous observation + mSys->SolveUpdateReset(); + + if (mShow) + StdOut() << aIterNumber + 1 << ", " << aSomDif / aTotalNumberOfInsidePixels + << ", " << aNbOut << std::endl; + } + + void cAppli_cTriangleDeformationTrRad::FillDisplacementMapsAndOutputImage(const cPtInsideTriangles &aLastPixInsideTriangle, + const cPt2dr &aLastTranslatedFilledPoint, + const tREAL8 aLastRadiometryTranslation, + const tREAL8 aLastRadiometryScaling) + { + const tREAL8 aLastXCoordinate = aLastPixInsideTriangle.GetCartesianCoordinates().x(); + const tREAL8 aLastYCoordinate = aLastPixInsideTriangle.GetCartesianCoordinates().y(); + + const cPt2di aLastCoordinate = cPt2di(aLastXCoordinate, aLastYCoordinate); + mDImDepX->SetV(aLastCoordinate, + aLastTranslatedFilledPoint.x() - aLastXCoordinate); + mDImDepY->SetV(aLastCoordinate, + aLastTranslatedFilledPoint.y() - aLastYCoordinate); + const tREAL8 aLastXTranslatedCoord = aLastXCoordinate + mDImDepX->GetV(aLastCoordinate); + const tREAL8 aLastYTranslatedCoord = aLastYCoordinate + mDImDepY->GetV(aLastCoordinate); + + const tREAL8 aLastRadiometryValue = aLastRadiometryScaling * aLastPixInsideTriangle.GetPixelValue() + + aLastRadiometryTranslation; + + // Build image with intensities displaced + // deal with different cases of pixel being translated out of image + if (aLastXTranslatedCoord < 0 && aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, 0))); + else if (aLastXTranslatedCoord >= mSzImOut.x() && aLastYTranslatedCoord >= mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, mSzImOut.y() - 1))); + else if (aLastXTranslatedCoord < 0 && aLastYTranslatedCoord >= mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, mSzImOut.y() - 1))); + else if (aLastXTranslatedCoord >= mSzImOut.x() && aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, 0))); + else if (aLastXTranslatedCoord >= 0 && aLastXTranslatedCoord < mSzImOut.x() && + aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(aLastXTranslatedCoord, 0))); + else if (aLastXTranslatedCoord >= 0 && aLastXTranslatedCoord < mSzImOut.x() && + aLastYTranslatedCoord > mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(aLastXTranslatedCoord, mSzImOut.y() - 1))); + else if (aLastYTranslatedCoord >= 0 && aLastYTranslatedCoord < mSzImOut.y() && + aLastXTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, aLastYTranslatedCoord))); + else if (aLastYTranslatedCoord >= 0 && aLastYTranslatedCoord < mSzImOut.y() && + aLastXTranslatedCoord > mSzImOut.x()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, aLastYTranslatedCoord))); + else + // at the translated pixel the untranslated pixel value is given computed with the right radiometry values + mDImOut->SetV(cPt2di(aLastXTranslatedCoord, aLastYTranslatedCoord), aLastRadiometryValue); + } + + void cAppli_cTriangleDeformationTrRad::GenerateDisplacementMapsAndOutputImages(const tDenseVect &aVFinalSol, const int aIterNumber) + { + mImOut = tIm(mSzImPre); + mDImOut = &mImOut.DIm(); + mSzImOut = cPt2di(mDImOut->Sz().x(), mDImOut->Sz().y()); + + mImDepX = tIm(mSzImPre, 0, eModeInitImage::eMIA_Null); + mDImDepX = &mImDepX.DIm(); + + mImDepY = tIm(mSzImPre, 0, eModeInitImage::eMIA_Null); + mDImDepY = &mImDepY.DIm(); + + tIm aLastPostIm = tIm(mSzImPost); + tDIm *aLastPostDIm = nullptr; + cAppli_cTriangleDeformation::LoadImageAndData(aLastPostIm, aLastPostDIm, "post", mImPre, mImPost); + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aLastPostIm = mImPre.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aLastPostDIm = &aLastPostIm.DIm(); + } + + tDoubleVect aLastVObs(12, 0.0); + + for (const cPt2di &aOutPix : *mDImOut) // Initialise output image + mDImOut->SetV(aOutPix, aLastPostDIm->GetV(aOutPix)); + + for (size_t aLTr = 0; aLTr < mDelTri.NbFace(); aLTr++) + { + const tTri2dr aLastTri = mDelTri.KthTri(aLTr); + const cPt3di aLastIndicesOfTriKnots = mDelTri.KthFace(aLTr); + + const cTriangle2DCompiled aLastCompTri(aLastTri); + + std::vector aLastVectorToFillWithInsidePixels; + aLastCompTri.PixelsInside(aLastVectorToFillWithInsidePixels); + + const tIntVect aLastVecInd = { + 4 * aLastIndicesOfTriKnots.x(), + 4 * aLastIndicesOfTriKnots.x() + 1, + 4 * aLastIndicesOfTriKnots.x() + 2, + 4 * aLastIndicesOfTriKnots.x() + 3, + 4 * aLastIndicesOfTriKnots.y(), + 4 * aLastIndicesOfTriKnots.y() + 1, + 4 * aLastIndicesOfTriKnots.y() + 2, + 4 * aLastIndicesOfTriKnots.y() + 3, + 4 * aLastIndicesOfTriKnots.z(), + 4 * aLastIndicesOfTriKnots.z() + 1, + 4 * aLastIndicesOfTriKnots.z() + 2, + 4 * aLastIndicesOfTriKnots.z() + 3, + }; + + const cPt2dr aLastTrPointA = cPt2dr(aVFinalSol(aLastVecInd.at(0)), + aVFinalSol(aLastVecInd.at(1))); // last translation 1st point of triangle + const cPt2dr aLastTrPointB = cPt2dr(aVFinalSol(aLastVecInd.at(4)), + aVFinalSol(aLastVecInd.at(5))); // last translation 2nd point of triangle + const cPt2dr aLastTrPointC = cPt2dr(aVFinalSol(aLastVecInd.at(8)), + aVFinalSol(aLastVecInd.at(9))); // last translation 3rd point of triangle + + const tREAL8 aLastRadTrPointA = aVFinalSol(aLastVecInd.at(2)); // last translation on radiometry 1st point of triangle + const tREAL8 aLastRadScPointA = aVFinalSol(aLastVecInd.at(3)); // last scale on radiometry 3rd point of triangle + const tREAL8 aLastRadTrPointB = aVFinalSol(aLastVecInd.at(6)); // last translation on radiometry 2nd point of triangle + const tREAL8 aLastRadScPointB = aVFinalSol(aLastVecInd.at(7)); // last scale on radiometry 3rd point of triangle + const tREAL8 aLastRadTrPointC = aVFinalSol(aLastVecInd.at(10)); // last translation on radiometry 3rd point of triangle + const tREAL8 aLastRadScPointC = aVFinalSol(aLastVecInd.at(11)); // last scale on radiometry 3rd point of triangle + + const size_t aLastNumberOfInsidePixels = aLastVectorToFillWithInsidePixels.size(); + + for (size_t aLastFilledPixel = 0; aLastFilledPixel < aLastNumberOfInsidePixels; aLastFilledPixel++) + { + const cPtInsideTriangles aLastPixInsideTriangle = cPtInsideTriangles(aLastCompTri, aLastVectorToFillWithInsidePixels, + aLastFilledPixel, *aLastPostDIm); + // prepare for barycenter translation formula by filling aVObs with different coordinates + FormalInterpBarycenter_SetObs(aLastVObs, 0, aLastPixInsideTriangle); + + // image of a point in triangle by current translation + const cPt2dr aLastTranslatedFilledPoint = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaToFilledPixel(aLastTrPointA, aLastTrPointB, + aLastTrPointC, aLastVObs); + + const tREAL8 aLastRadiometryTranslation = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForTranslationRadiometry(aLastRadTrPointA, + aLastRadTrPointB, + aLastRadTrPointC, + aLastVObs); + + const tREAL8 aLastRadiometryScaling = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaForScalingRadiometry(aLastRadScPointA, + aLastRadScPointB, + aLastRadScPointC, + aLastVObs); + + FillDisplacementMapsAndOutputImage(aLastPixInsideTriangle, aLastTranslatedFilledPoint, aLastRadiometryTranslation, + aLastRadiometryScaling); + } + } + + // save displacement maps in x and y to image files + if (mUseMultiScaleApproach) + { + mDImDepX->ToFile("DisplacedPixelsX_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + mDImDepY->ToFile("DisplacedPixelsY_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + if (aIterNumber == mNumberOfScales + mNumberOfEndIterations - 1) + mDImOut->ToFile("DisplacedPixels.tif"); + } + else + { + mDImDepX->ToFile("DisplacedPixelsX_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + mDImDepY->ToFile("DisplacedPixelsY_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + mDImOut->ToFile("DisplacedPixels.tif"); + } + } + + void cAppli_cTriangleDeformationTrRad::GenerateDisplacementMapsAndDisplayLastValuesUnknowns(const int aIterNumber, const bool aDisplayLastRadiometryValues, + const bool aDisplayLastTranslationValues) + { + tDenseVect aVFinalSol = mSys->CurGlobSol(); + + if (mGenerateDisplacementImage) + GenerateDisplacementMapsAndOutputImages(aVFinalSol, aIterNumber); + + if (aDisplayLastRadiometryValues || aDisplayLastTranslationValues) + cAppli_cTriangleDeformation::DisplayLastUnknownValues(aVFinalSol, aDisplayLastRadiometryValues, + aDisplayLastTranslationValues); + } + + void cAppli_cTriangleDeformationTrRad::DoOneIteration(const int aIterNumber) + { + LoopOverTrianglesAndUpdateParameters(aIterNumber); // Iterate over triangles and solve system + + // Show final translation results and produce displacement maps + if (mUseMultiScaleApproach) + { + if (aIterNumber == (mNumberOfScales + mNumberOfEndIterations - 1)) + GenerateDisplacementMapsAndDisplayLastValuesUnknowns(aIterNumber, mDisplayLastRadiometryValues, + mDisplayLastTranslationValues); + } + else + { + if (aIterNumber == (mNumberOfScales - 1)) + GenerateDisplacementMapsAndDisplayLastValuesUnknowns(aIterNumber, mDisplayLastRadiometryValues, + mDisplayLastTranslationValues); + } + } + + //----------------------------------------- + + int cAppli_cTriangleDeformationTrRad::Exe() + { + // read pre and post images and their sizes + mImPre = tIm::FromFile(mNamePreImage); + mImPost = tIm::FromFile(mNamePostImage); + + mDImPre = &mImPre.DIm(); + mSzImPre = mDImPre->Sz(); + + mDImPost = &mImPost.DIm(); + mSzImPost = mDImPost->Sz(); + + if (mUseMultiScaleApproach) + mSigmaGaussFilter = mNumberOfScales * mSigmaGaussFilterStep; + + if (mShow) + StdOut() << "Iter, " + << "Diff, " + << "NbOut" << std::endl; + + cAppli_cTriangleDeformation::GeneratePointsForDelaunay(mVectorPts, mNumberPointsToGenerate, mRandomUniformLawUpperBoundLines, + mRandomUniformLawUpperBoundCols, mDelTri, mSzImPre); + + cAppli_cTriangleDeformation::InitialisationAfterExe(mDelTri, mSys); + + if (mUseMultiScaleApproach) + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales + mNumberOfEndIterations; aIterNumber++) + DoOneIteration(aIterNumber); + } + else + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales; aIterNumber++) + DoOneIteration(aIterNumber); + } + + return EXIT_SUCCESS; + } + + /********************************************/ + // ::MMVII // + /********************************************/ + + tMMVII_UnikPApli Alloc_cTriangleDeformationTrRad(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) + { + return tMMVII_UnikPApli(new cAppli_cTriangleDeformationTrRad(aVArgs, aSpec)); + } + + cSpecMMVII_Appli TheSpec_ComputeTriangleDeformationTrRad( + "ComputeTriangleDeformationTrRad", + Alloc_cTriangleDeformationTrRad, + "Compute 2D deformation of triangles between images using triangles and alternative formula", + {eApF::ImProc}, // category + {eApDT::Image}, // input + {eApDT::Image}, // output + __FILE__); + +}; // namespace MMVII \ No newline at end of file diff --git a/MMVII/src/MeshDisplacement/TriangleDeformationTrRad.h b/MMVII/src/MeshDisplacement/TriangleDeformationTrRad.h new file mode 100644 index 0000000000..5af78c8007 --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformationTrRad.h @@ -0,0 +1,113 @@ +#ifndef _TRIANGLEDEFORMATIONTRAD_H_ +#define _TRIANGLEDEFORMATIONTRAD_H_ + +#include "MMVII_Geom2D.h" +#include "MMVII_PhgrDist.h" +#include "TriangleDeformation.h" + +using namespace NS_SymbolicDerivative; + +namespace MMVII +{ + + /******************************************/ + /* */ + /* cTriangleDeformation */ + /* */ + /******************************************/ + + class cAppli_cTriangleDeformationTrRad : public cAppli_cTriangleDeformation + { + public: + typedef cIm2D tIm; + typedef cDataIm2D tDIm; + typedef cTriangle tTri2dr; + typedef cDenseVect tDenseVect; + typedef std::vector tIntVect; + typedef std::vector tDoubleVect; + + cAppli_cTriangleDeformationTrRad(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec); + ~cAppli_cTriangleDeformationTrRad(); + + int Exe() override; + cCollecSpecArg2007 &ArgObl(cCollecSpecArg2007 &anArgObl) override; + cCollecSpecArg2007 &ArgOpt(cCollecSpecArg2007 &anArgOpt) override; + + // Iterate of triangles and inside pixels + virtual void DoOneIteration(const int aIterNumber); + // Loops over all triangles and solves system to update parameters at end of iteration + virtual void LoopOverTrianglesAndUpdateParameters(const int aIterNumber); + // Generate displacement maps of last solution + virtual void GenerateDisplacementMapsAndOutputImages(const tDenseVect &aVFinalSol, const int aIterNumber); + // Generates Displacement maps and coordinates of points in triangulation at last iteration + virtual void GenerateDisplacementMapsAndDisplayLastValuesUnknowns(const int aIterNumber, const bool aDisplayLastRadiometryValues, + const bool aDisplayLastTranslationValues); + // Fill displacement maps and output image + virtual void FillDisplacementMapsAndOutputImage(const cPtInsideTriangles &aLastPixInsideTriangle, + const cPt2dr &aLastTranslatedFilledPoint, + const tREAL8 aLastRadiometryTranslation, + const tREAL8 aLastRadiometryScaling); + + private: + // == Mandatory args ==== + + std::string mNamePreImage; // Name of given pre-image + std::string mNamePostImage; // Name of given post-image + int mNumberPointsToGenerate; // number of generated points + int mNumberOfScales; // number of iterations in optimisation process + + // == Optionnal args ==== + + int mRandomUniformLawUpperBoundLines; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundLines [ + int mRandomUniformLawUpperBoundCols; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundCols [ + bool mShow; // Print result, export image ... + bool mComputeAvgMax; // Compute average and maximum pixel value of difference image between pre and post images + bool mUseMultiScaleApproach; // Apply multi-scale approach or not + int mSigmaGaussFilterStep; // Decreasing step of sigma value during iterations + bool mGenerateDisplacementImage; // Generate image with displaced pixels + bool mFreezeTranslationX; // Freeze x-translation or not during optimisation + bool mFreezeTranslationY; // Freeze y-translation or not during optimisation + bool mFreezeRadTranslation; // Freeze radiometry translation in computation + bool mFreezeRadScale; // Freeze radiometry scaling in computation + double mWeightRadTranslation; // Weight given to radiometry translation if soft freezing is applied (default : negative => not applied) + double mWeightRadScale; // Weight given to radiometry scaling if soft freezing is applied (default : negative => not applied) + int mNumberOfIterGaussFilter; // Number of iterations to be done in Gauss filter algorithm + int mNumberOfEndIterations; // Number of iterations to do while using original image in multi-scale approach + bool mDisplayLastTranslationValues; // Whether to display the final coordinates of the translated points + bool mDisplayLastRadiometryValues; // Display final values of radiometry unknowns at last iteration of optimisation process + + // == Internal variables ==== + + cPt2di mSzImPre; // size of image + tIm mImPre; // memory representation of the image + tDIm *mDImPre; // memory representation of the image + + cPt2di mSzImPost; // size of image + tIm mImPost; // memory representation of the image + tDIm *mDImPost; // memory representation of the image + + cPt2di mSzImOut; // size of image + tIm mImOut; // memory representation of the image + tDIm *mDImOut; // memory representation of the image + + cPt2di mSzImDepX; // size of image + tIm mImDepX; // memory representation of the image + tDIm *mDImDepX; // memory representation of the image + + cPt2di mSzImDepY; // size of image + tIm mImDepY; // memory representation of the image + tDIm *mDImDepY; // memory representation of the image + + std::vector mVectorPts; // A vector containing a set of points + cTriangulation2D mDelTri; // A Delaunay triangle + + double mSigmaGaussFilter; // Value of sigma in gauss filter + bool mIsLastIters; // Determines whether optimisation process is at last iters to optimise on original image + + cResolSysNonLinear *mSys; // Non Linear Sys for solving problem + cCalculator *mEqTriDeformTrRad; // calculator giving access to values and derivatives + }; +} + +#endif // _TRIANGLEDEFORMATIONTRAD_H_ diff --git a/MMVII/src/MeshDisplacement/TriangleDeformationTranslation.cpp b/MMVII/src/MeshDisplacement/TriangleDeformationTranslation.cpp new file mode 100644 index 0000000000..c3da45f140 --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformationTranslation.cpp @@ -0,0 +1,463 @@ +#include "cMMVII_Appli.h" + +#include "MMVII_TplSymbTriangle.h" + +#include "TriangleDeformationTranslation.h" + +/** + \file TriangleDeformationTranslation.cpp + + \brief file for computing 2D translation between 2 images + thanks to triangles. +**/ + +namespace MMVII +{ + + /******************************************/ + /* */ + /* cTriangleDeformationTranslation */ + /* */ + /******************************************/ + + cAppli_cTriangleDeformationTranslation::cAppli_cTriangleDeformationTranslation(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) : cAppli_cTriangleDeformation(aVArgs, aSpec), + mRandomUniformLawUpperBoundLines(1), + mRandomUniformLawUpperBoundCols(1), + mShow(true), + mUseMultiScaleApproach(true), + mSigmaGaussFilterStep(1), + mGenerateDisplacementImage(true), + mNumberOfIterGaussFilter(3), + mNumberOfEndIterations(2), + mDisplayLastTranslatedPointsCoordinates(false), + mSzImPre(cPt2di(1, 1)), + mImPre(mSzImPre), + mDImPre(nullptr), + mSzImPost(cPt2di(1, 1)), + mImPost(mSzImPost), + mDImPost(nullptr), + mSzImOut(cPt2di(1, 1)), + mImOut(mSzImOut), + mDImOut(nullptr), + mSzImDepX(cPt2di(1, 1)), + mImDepX(mSzImDepX), + mDImDepX(nullptr), + mSzImDepY(cPt2di(1, 1)), + mImDepY(mSzImDepY), + mDImDepY(nullptr), + mVectorPts({cPt2dr(0, 0)}), + mDelTri(mVectorPts), + mSys(nullptr), + mEqTranslationTri(nullptr) + { + mEqTranslationTri = EqDeformTriTranslation(true, 1); // true means with derivative, 1 is size of buffer + // mEqTranslationTri->SetDebugEnabled(true); + } + + cAppli_cTriangleDeformationTranslation::~cAppli_cTriangleDeformationTranslation() + { + delete mSys; + delete mEqTranslationTri; + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformationTranslation::ArgObl(cCollecSpecArg2007 &anArgObl) + { + return anArgObl + << Arg2007(mNamePreImage, "Name of pre-image file.", {{eTA2007::FileImage}, {eTA2007::FileDirProj}}) + << Arg2007(mNamePostImage, "Name of post-image file.", {eTA2007::FileImage}) + << Arg2007(mNumberPointsToGenerate, "Number of points you want to generate for triangulation.") + << Arg2007(mNumberOfScales, "Total number of scales to run in multi-scale approach optimisation process."); + } + + cCollecSpecArg2007 &cAppli_cTriangleDeformationTranslation::ArgOpt(cCollecSpecArg2007 &anArgOpt) + { + return anArgOpt + << AOpt2007(mRandomUniformLawUpperBoundCols, "RandomUniformLawUpperBoundXAxis", + "Maximum value that the uniform law can draw from on the x-axis.", {eTA2007::HDV}) + << AOpt2007(mRandomUniformLawUpperBoundLines, "RandomUniformLawUpperBoundYAxis", + "Maximum value that the uniform law can draw from for on the y-axis.", {eTA2007::HDV}) + << AOpt2007(mShow, "Show", "Whether to print minimisation results.", {eTA2007::HDV}) + << AOpt2007(mUseMultiScaleApproach, "UseMultiScaleApproach", "Whether to use multi-scale approach or not.", {eTA2007::HDV}) + << AOpt2007(mSigmaGaussFilterStep, "SigmaGaussFilterStep", "Sigma value to use for Gauss filter in multi-stage approach.", {eTA2007::HDV}) + << AOpt2007(mGenerateDisplacementImage, "GenerateDisplacementImage", + "Whether to generate and save an image having been translated.", {eTA2007::HDV}) + << AOpt2007(mDisplayLastTranslatedPointsCoordinates, "DisplayLastTranslatedCoordinates", + "Whether to display the final coordinates of the trainslated points.", {eTA2007::HDV}) + << AOpt2007(mNumberOfIterGaussFilter, "NumberOfIterationsGaussFilter", + "Number of iterations to run in Gauss filter algorithm.", {eTA2007::HDV}) + << AOpt2007(mNumberOfEndIterations, "NumberOfEndIterations", + "Number of iterations to run on original images in multi-scale approach.", {eTA2007::HDV}); + } + + void cAppli_cTriangleDeformationTranslation::InitialisationAfterExeTranslation() + { + tDenseVect aVInit(2 * mDelTri.NbPts(), eModeInitImage::eMIA_Null); + + mSys = new cResolSysNonLinear(eModeSSR::eSSR_LsqDense, aVInit); + } + + void cAppli_cTriangleDeformationTranslation::LoopOverTrianglesAndUpdateParametersTranslation(const int aIterNumber) + { + //----------- allocate vec of obs : + tDoubleVect aVObs(12, 0.0); // 6 for ImagePre and 6 for ImagePost + + //----------- extract current parameters + tDenseVect aVCur = mSys->CurGlobSol(); // Get current solution. + + tIm aCurPreIm = tIm(mSzImPre); + tDIm *aCurPreDIm = nullptr; + tIm aCurPostIm = tIm(mSzImPost); + tDIm *aCurPostDIm = nullptr; + + mIsLastIters = false; + + if (mUseMultiScaleApproach) + mIsLastIters = cAppli_cTriangleDeformation::ManageDifferentCasesOfEndIterations(aIterNumber, mNumberOfScales, mNumberOfEndIterations, + mIsLastIters, mImPre, mImPost, aCurPreIm, aCurPreDIm, + aCurPostIm, aCurPostDIm); + else + { + cAppli_cTriangleDeformation::LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + cAppli_cTriangleDeformation::LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aCurPreIm = mImPre.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aCurPostIm = mImPost.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + + aCurPreDIm = &aCurPreIm.DIm(); + aCurPostDIm = &aCurPostIm.DIm(); + + mSigmaGaussFilter -= 1; + + const bool aSaveGaussImage = false; + if (aSaveGaussImage) + aCurPreDIm->ToFile("GaussFilteredImPre_iter_" + std::to_string(aIterNumber) + ".tif"); + } + else if (mUseMultiScaleApproach && mIsLastIters) + { + cAppli_cTriangleDeformation::LoadImageAndData(aCurPreIm, aCurPreDIm, "pre", mImPre, mImPost); + cAppli_cTriangleDeformation::LoadImageAndData(aCurPostIm, aCurPostDIm, "post", mImPre, mImPost); + } + + //----------- declaration of indicator of convergence + tREAL8 aSomDif = 0; // sum of difference between untranslated pixel and translated one. + size_t aNbOut = 0; // number of translated pixels out of image + + // Count number of pixels inside triangles for normalisation + size_t aTotalNumberOfInsidePixels = 0; + + // Loop over all triangles to add the observations on each point + for (size_t aTr = 0; aTr < mDelTri.NbFace(); aTr++) + { + const tTri2dr aTri = mDelTri.KthTri(aTr); + const cPt3di aIndicesOfTriKnots = mDelTri.KthFace(aTr); + + const cTriangle2DCompiled aCompTri(aTri); + + std::vector aVectorToFillWithInsidePixels; + aCompTri.PixelsInside(aVectorToFillWithInsidePixels); // get pixels inside triangle + + //----------- index of unknown, finds the associated pixels of current triangle + const tIntVect aVecInd = {2 * aIndicesOfTriKnots.x(), 2 * aIndicesOfTriKnots.x() + 1, + 2 * aIndicesOfTriKnots.y(), 2 * aIndicesOfTriKnots.y() + 1, + 2 * aIndicesOfTriKnots.z(), 2 * aIndicesOfTriKnots.z() + 1}; + + const cPt2dr aCurTrPointA = cPt2dr(aVCur(aVecInd.at(0)), + aVCur(aVecInd.at(1))); // current translation 1st point of triangle + const cPt2dr aCurTrPointB = cPt2dr(aVCur(aVecInd.at(2)), + aVCur(aVecInd.at(3))); // current translation 2nd point of triangle + const cPt2dr aCurTrPointC = cPt2dr(aVCur(aVecInd.at(4)), + aVCur(aVecInd.at(5))); // current translation 3rd point of triangle + + const size_t aNumberOfInsidePixels = aVectorToFillWithInsidePixels.size(); + + // Loop over all pixels inside triangle + // size_t is necessary as there can be a lot of pixels in triangles + for (size_t aFilledPixel = 0; aFilledPixel < aNumberOfInsidePixels; aFilledPixel++) + { + const cPtInsideTriangles aPixInsideTriangle = cPtInsideTriangles(aCompTri, aVectorToFillWithInsidePixels, + aFilledPixel, *aCurPreDIm); + // prepare for barycenter translation formula by filling aVObs with different coordinates + FormalInterpBarycenter_SetObs(aVObs, 0, aPixInsideTriangle); + + // image of a point in triangle by current translation + const cPt2dr aTranslatedFilledPoint = cAppli_cTriangleDeformation::ApplyBarycenterTranslationFormulaToFilledPixel(aCurTrPointA, aCurTrPointB, + aCurTrPointC, aVObs); + + if (aCurPostDIm->InsideBL(aTranslatedFilledPoint)) // avoid errors + { + // prepare for application of bilinear formula + FormalBilinTri_SetObs(aVObs, TriangleDisplacement_NbObs, aTranslatedFilledPoint, *aCurPostDIm); + + // Now add observation + mSys->CalcAndAddObs(mEqTranslationTri, aVecInd, aVObs); + + // compute indicators + const tREAL8 aDif = aVObs[5] - aCurPostDIm->GetVBL(aTranslatedFilledPoint); // residual - aValueImPre - aCurPostDIm->GetVBL(aTranslatedFilledPoint) + aSomDif += std::abs(aDif); + } + else + aNbOut++; // Count number of pixels translated outside post image + + aTotalNumberOfInsidePixels += aNumberOfInsidePixels; + } + } + + if (mUseMultiScaleApproach && !mIsLastIters && aIterNumber != 0) + { + const bool aGenerateIntermediateMaps = false; + if (aGenerateIntermediateMaps) + GenerateDisplacementMaps(aVCur, aIterNumber); + } + + // Update all parameter taking into account previous observation + mSys->SolveUpdateReset(); + + if (mShow) + StdOut() << aIterNumber + 1 << ", " << aSomDif / aTotalNumberOfInsidePixels + << ", " << aNbOut << std::endl; + } + + cPt2dr cAppli_cTriangleDeformationTranslation::ApplyLastBarycenterTranslationFormulaToInsidePixel(const cPt2dr &aLastTranslationPointA, + const cPt2dr &aLastTranslationPointB, + const cPt2dr &aLastTranslationPointC, + const cPtInsideTriangles &aLastPixInsideTriangle) + { + const cPt2dr aLastCartesianCoordinates = aLastPixInsideTriangle.GetCartesianCoordinates(); + const cPt3dr aLastBarycenterCoordinates = aLastPixInsideTriangle.GetBarycenterCoordinates(); + // apply current barycenter translation formula for x and y on current observations. + const tREAL8 aLastXTriCoord = aLastCartesianCoordinates.x() + aLastBarycenterCoordinates.x() * aLastTranslationPointA.x() + + aLastBarycenterCoordinates.y() * aLastTranslationPointB.x() + aLastBarycenterCoordinates.z() * aLastTranslationPointC.x(); + const tREAL8 aLastYTriCoord = aLastCartesianCoordinates.y() + aLastBarycenterCoordinates.x() * aLastTranslationPointA.y() + + aLastBarycenterCoordinates.y() * aLastTranslationPointB.y() + aLastBarycenterCoordinates.z() * aLastTranslationPointC.y(); + + const cPt2dr aLastTranslatedPixel = cPt2dr(aLastXTriCoord, aLastYTriCoord); + + return aLastTranslatedPixel; + } + + void cAppli_cTriangleDeformationTranslation::FillDisplacementMaps(const cPtInsideTriangles &aLastPixInsideTriangle, + const cPt2dr &aLastTranslatedFilledPoint) + { + const tREAL8 aLastXCoordinate = aLastPixInsideTriangle.GetCartesianCoordinates().x(); + const tREAL8 aLastYCoordinate = aLastPixInsideTriangle.GetCartesianCoordinates().y(); + const tREAL8 aLastPixelValue = aLastPixInsideTriangle.GetPixelValue(); + + const cPt2di aLastCoordinate = cPt2di(aLastXCoordinate, aLastYCoordinate); + + mDImDepX->SetV(aLastCoordinate, + aLastTranslatedFilledPoint.x() - aLastXCoordinate); + mDImDepY->SetV(aLastCoordinate, + aLastTranslatedFilledPoint.y() - aLastYCoordinate); + const tREAL8 aLastXTranslatedCoord = aLastXCoordinate + mDImDepX->GetV(aLastCoordinate); + const tREAL8 aLastYTranslatedCoord = aLastYCoordinate + mDImDepY->GetV(aLastCoordinate); + + // Build image with intensities displaced + // deal with different cases of pixel being translated out of image + if (aLastXTranslatedCoord < 0 && aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, 0))); + else if (aLastXTranslatedCoord >= mSzImOut.x() && aLastYTranslatedCoord >= mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, mSzImOut.y() - 1))); + else if (aLastXTranslatedCoord < 0 && aLastYTranslatedCoord >= mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, mSzImOut.y() - 1))); + else if (aLastXTranslatedCoord >= mSzImOut.x() && aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, 0))); + else if (aLastXTranslatedCoord >= 0 && aLastXTranslatedCoord < mSzImOut.x() && + aLastYTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(aLastXTranslatedCoord, 0))); + else if (aLastXTranslatedCoord >= 0 && aLastXTranslatedCoord < mSzImOut.x() && + aLastYTranslatedCoord > mSzImOut.y()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(aLastXTranslatedCoord, mSzImOut.y() - 1))); + else if (aLastYTranslatedCoord >= 0 && aLastYTranslatedCoord < mSzImOut.y() && + aLastXTranslatedCoord < 0) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(0, aLastYTranslatedCoord))); + else if (aLastYTranslatedCoord >= 0 && aLastYTranslatedCoord < mSzImOut.y() && + aLastXTranslatedCoord > mSzImOut.x()) + mDImOut->SetV(aLastCoordinate, mDImOut->GetV(cPt2di(mSzImOut.x() - 1, aLastYTranslatedCoord))); + else + // at the translated pixel the untranslated pixel value is given + mDImOut->SetV(cPt2di(aLastXTranslatedCoord, aLastYTranslatedCoord), aLastPixelValue); + } + + void cAppli_cTriangleDeformationTranslation::GenerateDisplacementMaps(const tDenseVect &aVFinalSol, const int aIterNumber) + { + mImOut = tIm(mSzImPre); + mDImOut = &mImOut.DIm(); + mSzImOut = cPt2di(mDImOut->Sz().x(), mDImOut->Sz().y()); + + mImDepX = tIm(mSzImPre, 0, eModeInitImage::eMIA_Null); + mDImDepX = &mImDepX.DIm(); + + mImDepY = tIm(mSzImPre, 0, eModeInitImage::eMIA_Null); + mDImDepY = &mImDepY.DIm(); + + tIm aLastPreIm = tIm(mSzImPre); + tDIm *aLastPreDIm = nullptr; + cAppli_cTriangleDeformation::LoadImageAndData(aLastPreIm, aLastPreDIm, "pre", mImPre, mImPost); + + if (mUseMultiScaleApproach && !mIsLastIters) + { + aLastPreIm = mImPre.GaussFilter(mSigmaGaussFilter, mNumberOfIterGaussFilter); + aLastPreDIm = &aLastPreIm.DIm(); + } + + for (const cPt2di &aOutPix : *mDImOut) // Initialise output image + mDImOut->SetV(aOutPix, aLastPreDIm->GetV(aOutPix)); + + for (size_t aLTr = 0; aLTr < mDelTri.NbFace(); aLTr++) + { + const tTri2dr aLastTri = mDelTri.KthTri(aLTr); + const cPt3di aLastIndicesOfTriKnots = mDelTri.KthFace(aLTr); + + const cTriangle2DCompiled aLastCompTri(aLastTri); + + std::vector aLastVectorToFillWithInsidePixels; + aLastCompTri.PixelsInside(aLastVectorToFillWithInsidePixels); + + const tIntVect aLastVecInd = {2 * aLastIndicesOfTriKnots.x(), 2 * aLastIndicesOfTriKnots.x() + 1, + 2 * aLastIndicesOfTriKnots.y(), 2 * aLastIndicesOfTriKnots.y() + 1, + 2 * aLastIndicesOfTriKnots.z(), 2 * aLastIndicesOfTriKnots.z() + 1}; + + const cPt2dr aLastTrPointA = cPt2dr(aVFinalSol(aLastVecInd.at(0)), + aVFinalSol(aLastVecInd.at(1))); // last translation 1st point of triangle + const cPt2dr aLastTrPointB = cPt2dr(aVFinalSol(aLastVecInd.at(2)), + aVFinalSol(aLastVecInd.at(3))); // last translation 2nd point of triangle + const cPt2dr aLastTrPointC = cPt2dr(aVFinalSol(aLastVecInd.at(4)), + aVFinalSol(aLastVecInd.at(5))); // last translation 3rd point of triangle + + const size_t aLastNumberOfInsidePixels = aLastVectorToFillWithInsidePixels.size(); + + for (size_t aLastFilledPixel = 0; aLastFilledPixel < aLastNumberOfInsidePixels; aLastFilledPixel++) + { + const cPtInsideTriangles aLastPixInsideTriangle = cPtInsideTriangles(aLastCompTri, aLastVectorToFillWithInsidePixels, + aLastFilledPixel, *aLastPreDIm); + + // image of a point in triangle by current translation + const cPt2dr aLastTranslatedFilledPoint = ApplyLastBarycenterTranslationFormulaToInsidePixel(aLastTrPointA, aLastTrPointB, + aLastTrPointC, aLastPixInsideTriangle); + + FillDisplacementMaps(aLastPixInsideTriangle, aLastTranslatedFilledPoint); + } + } + + // save displacement maps in x and y to image files + if (mUseMultiScaleApproach) + { + mDImDepX->ToFile("DisplacedPixelsX_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + mDImDepY->ToFile("DisplacedPixelsY_iter_" + std::to_string(aIterNumber) + "_" + + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales + mNumberOfEndIterations) + ".tif"); + if (aIterNumber == mNumberOfScales + mNumberOfEndIterations - 1) + mDImOut->ToFile("DisplacedPixels.tif"); + } + else + { + mDImDepX->ToFile("DisplacedPixelsX_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + mDImDepY->ToFile("DisplacedPixelsY_" + std::to_string(mNumberPointsToGenerate) + "_" + + std::to_string(mNumberOfScales) + ".tif"); + mDImOut->ToFile("DisplacedPixels.tif"); + } + } + + void cAppli_cTriangleDeformationTranslation::GenerateDisplacementMapsAndLastTranslatedPoints(const int aIterNumber) + { + tDenseVect aVFinalSol = mSys->CurGlobSol(); + + if (mGenerateDisplacementImage) + GenerateDisplacementMaps(aVFinalSol, aIterNumber); + + if (mDisplayLastTranslatedPointsCoordinates) + { + for (int aFinalUnk = 0; aFinalUnk < aVFinalSol.DIm().Sz(); aFinalUnk++) + { + StdOut() << aVFinalSol(aFinalUnk) << " "; + if (aFinalUnk % 2 == 1 && aFinalUnk != 0) + StdOut() << std::endl; + } + } + } + + void cAppli_cTriangleDeformationTranslation::DoOneIterationTranslation(const int aIterNumber) + { + LoopOverTrianglesAndUpdateParametersTranslation(aIterNumber); // Iterate over triangles and solve system + + // Show final translation results and produce displacement maps + if (mUseMultiScaleApproach) + { + if (aIterNumber == (mNumberOfScales + mNumberOfEndIterations - 1)) + GenerateDisplacementMapsAndLastTranslatedPoints(aIterNumber); + } + else + { + if (aIterNumber == (mNumberOfScales - 1)) + GenerateDisplacementMapsAndLastTranslatedPoints(aIterNumber); + } + } + + //----------------------------------------- + + int cAppli_cTriangleDeformationTranslation::Exe() + { + // read pre and post images and their sizes + mImPre = tIm::FromFile(mNamePreImage); + mImPost = tIm::FromFile(mNamePostImage); + + mDImPre = &mImPre.DIm(); + mSzImPre = mDImPre->Sz(); + + mDImPost = &mImPost.DIm(); + mSzImPost = mDImPost->Sz(); + + if (mUseMultiScaleApproach) + mSigmaGaussFilter = mNumberOfScales * mSigmaGaussFilterStep; + + if (mShow) + StdOut() << "Iter, " + << "Diff, " + << "NbOut" << std::endl; + + cAppli_cTriangleDeformation::GeneratePointsForDelaunay(mVectorPts, mNumberPointsToGenerate, mRandomUniformLawUpperBoundLines, + mRandomUniformLawUpperBoundCols, mDelTri, mSzImPre); + + InitialisationAfterExeTranslation(); + + if (mUseMultiScaleApproach) + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales + mNumberOfEndIterations; aIterNumber++) + DoOneIterationTranslation(aIterNumber); + } + else + { + for (int aIterNumber = 0; aIterNumber < mNumberOfScales; aIterNumber++) + DoOneIterationTranslation(aIterNumber); + } + + return EXIT_SUCCESS; + } + + /********************************************/ + // ::MMVII // + /********************************************/ + + tMMVII_UnikPApli Alloc_cTriangleDeformationTranslation(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec) + { + return tMMVII_UnikPApli(new cAppli_cTriangleDeformationTranslation(aVArgs, aSpec)); + } + + cSpecMMVII_Appli TheSpec_ComputeTriangleDeformationTranslation( + "ComputeTriangleDeformationTranslation", + Alloc_cTriangleDeformationTranslation, + "Compute 2D translation deformations of triangles between images using triangles", + {eApF::ImProc}, // category + {eApDT::Image}, // input + {eApDT::Image}, // output + __FILE__); + +}; // namespace MMVII \ No newline at end of file diff --git a/MMVII/src/MeshDisplacement/TriangleDeformationTranslation.h b/MMVII/src/MeshDisplacement/TriangleDeformationTranslation.h new file mode 100644 index 0000000000..170687e39c --- /dev/null +++ b/MMVII/src/MeshDisplacement/TriangleDeformationTranslation.h @@ -0,0 +1,111 @@ +#ifndef _TRIANGLEDEFORMATIONTRANSLATION_H_ +#define _TRIANGLEDEFORMATIONTRANSLATION_H_ + +#include "MMVII_Geom2D.h" +#include "MMVII_PhgrDist.h" + +#include "TriangleDeformation.h" + +using namespace NS_SymbolicDerivative; + +namespace MMVII +{ + + /******************************************/ + /* */ + /* cTriangleDeformationTranslation */ + /* */ + /******************************************/ + + class cAppli_cTriangleDeformationTranslation : public cAppli_cTriangleDeformation + { + public: + typedef cIm2D tIm; + typedef cDataIm2D tDIm; + typedef cTriangle tTri2dr; + typedef cDenseVect tDenseVect; + typedef std::vector tIntVect; + typedef std::vector tDoubleVect; + + cAppli_cTriangleDeformationTranslation(const std::vector &aVArgs, + const cSpecMMVII_Appli &aSpec); + ~cAppli_cTriangleDeformationTranslation(); + + int Exe() override; + cCollecSpecArg2007 &ArgObl(cCollecSpecArg2007 &anArgObl) override; + cCollecSpecArg2007 &ArgOpt(cCollecSpecArg2007 &anArgOpt) override; + + // Iterate of triangles and inside pixels + void DoOneIterationTranslation(const int aIterNumber); + // Loops over all triangles and solves system to update parameters at end of iteration + void LoopOverTrianglesAndUpdateParametersTranslation(const int aIterNumber); + // Generate displacement maps of last solution + void GenerateDisplacementMaps(const tDenseVect &aVFinalSol, const int aIterNumber); + // Generates Displacement maps and coordinates of points in triangulation at last iteration + void GenerateDisplacementMapsAndLastTranslatedPoints(const int aIterNumber); + // Initialise problem after user has input information + void InitialisationAfterExeTranslation(); + // Fill displacement maps and output image + void FillDisplacementMaps(const cPtInsideTriangles &aLastPixInsideTriangle, + const cPt2dr &aLastTranslatedFilledPoint); + // Apply barycentration translation formula to last translation values in optimisation process + cPt2dr ApplyLastBarycenterTranslationFormulaToInsidePixel(const cPt2dr &aLastTranslationPointA, + const cPt2dr &aLastTranslationPointB, + const cPt2dr &aLastTranslationPointC, + const cPtInsideTriangles &aLastPixInsideTriangle); + + private: + // == Mandatory args ==== + + std::string mNamePreImage; // Name of given pre-image + std::string mNamePostImage; // Name of given post-image + int mNumberPointsToGenerate; // number of generated points + int mNumberOfScales; // number of iterations in optimisation process + + // == Optionnal args ==== + + int mRandomUniformLawUpperBoundLines; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundLines [ + int mRandomUniformLawUpperBoundCols; // Uniform law generates random coordinates in interval [0, mRandomUniformLawUpperBoundCols [ + bool mShow; // Print result, export image ... + bool mComputeAvgMax; // Compute average and maximum pixel value of difference image between pre and post images + bool mUseMultiScaleApproach; // Apply multi-scale approach or not + int mSigmaGaussFilterStep; // Decreasing step of sigma value during iterations + bool mGenerateDisplacementImage; // Generate image with displaced pixels + int mNumberOfIterGaussFilter; // Number of iterations to be done in Gauss filter algorithm + int mNumberOfEndIterations; // Number of iterations to do while using original image in multi-scale approach + bool mDisplayLastTranslatedPointsCoordinates; // Whether to display the final coordinates of the translated points + + // == Internal variables ==== + + cPt2di mSzImPre; // size of image + tIm mImPre; // memory representation of the image + tDIm *mDImPre; // memory representation of the image + + cPt2di mSzImPost; // size of image + tIm mImPost; // memory representation of the image + tDIm *mDImPost; // memory representation of the image + + cPt2di mSzImOut; // size of image + tIm mImOut; // memory representation of the image + tDIm *mDImOut; // memory representation of the image + + cPt2di mSzImDepX; // size of image + tIm mImDepX; // memory representation of the image + tDIm *mDImDepX; // memory representation of the image + + cPt2di mSzImDepY; // size of image + tIm mImDepY; // memory representation of the image + tDIm *mDImDepY; // memory representation of the image + + double mSigmaGaussFilter; // Value of sigma in gauss filter + bool mIsLastIters; // Determines whether optimisation process is at last iters to optimise on original image + + std::vector mVectorPts; // A vector containing a set of points + cTriangulation2D mDelTri; // A Delaunay triangle + + cResolSysNonLinear *mSys; // Non Linear Sys for solving problem + cCalculator *mEqTranslationTri; // calculator giving access to values and derivatives + }; +} + +#endif // _TRIANGLEDEFORMATIONTRANSLATION_H_ \ No newline at end of file diff --git a/MMVII/src/SymbDerGen/Formulas_TrianglesDeform.h b/MMVII/src/SymbDerGen/Formulas_TrianglesDeform.h new file mode 100644 index 0000000000..053cbad86d --- /dev/null +++ b/MMVII/src/SymbDerGen/Formulas_TrianglesDeform.h @@ -0,0 +1,324 @@ +#ifndef _FORMULAS_TRIANGLES_DEFORM_H_ +#define _FORMULAS_TRIANGLES_DEFORM_H_ + +#include "MMVII_TplSymbTriangle.h" +#include "MMVII_util_tpl.h" + +#include "SymbDer/SymbolicDerivatives.h" +#include "SymbDer/SymbDer_MACRO.h" + +/** + \file Formulas_TrianglesDeform.h + \brief class to generate code for triangle deformations by minimization +**/ + +using namespace NS_SymbolicDerivative; + +namespace MMVII +{ + class cTriangleDeformation + { + public: + cTriangleDeformation() + { + } + + static const std::vector VNamesUnknowns() { return std::vector{"GeomTrXPointA", "GeomTrYPointA", "RadTranslationPointA", "RadScalingPointA", + "GeomTrXPointB", "GeomTrYPointB", "RadTranslationPointB", "RadScalingPointB", + "GeomTrXPointC", "GeomTrYPointC", "RadTranslationPointC", "RadScalingPointC"}; } + static const std::vector VNamesObs() + { + return Append( + std::vector{"PixelCoordinatesX", "PixelCoordinatesY", "AlphaCoordPixel", "BetaCoordPixel", "GammaCoordPixel", "IntensityImPre"}, + FormalBilinIm2D_NameObs("T") // 6 obs for bilinear interpol of Im + ); + } + + std::string FormulaName() const { return "TriangleDeformation"; } + + template + static std::vector formula( + const std::vector &aVUnk, + const std::vector &aVObs) + { + + // extract observation + const auto &aXCoordinate = aVObs[0]; + const auto &aYCoordinate = aVObs[1]; + const auto &aAlphaCoordinate = aVObs[2]; + const auto &aBetaCoordinate = aVObs[3]; + const auto &aGammaCoordinate = aVObs[4]; + const auto &aIntensityImPre = aVObs[5]; + + // extract unknowns + const auto &aGeomTrXPointA = aVUnk[0]; + const auto &aGeomTrYPointA = aVUnk[1]; + const auto &aRadTranslationPointA = aVUnk[2]; + const auto &aRadScalingPointA = aVUnk[3]; + const auto &aGeomTrXPointB = aVUnk[4]; + const auto &aGeomTrYPointB = aVUnk[5]; + const auto &aRadTranslationPointB = aVUnk[6]; + const auto &aRadScalingPointB = aVUnk[7]; + const auto &aGeomTrXPointC = aVUnk[8]; + const auto &aGeomTrYPointC = aVUnk[9]; + const auto &aRadTranslationPointC = aVUnk[10]; + const auto &aRadScalingPointC = aVUnk[11]; + + // Apply barycentric translation formula to coordinates + auto aXTri = aXCoordinate + aAlphaCoordinate * aGeomTrXPointA + aBetaCoordinate * aGeomTrXPointB + aGammaCoordinate * aGeomTrXPointC; + auto aYTri = aYCoordinate + aAlphaCoordinate * aGeomTrYPointA + aBetaCoordinate * aGeomTrYPointB + aGammaCoordinate * aGeomTrYPointC; + + // Apply barycentric interpolation to radiometric factors + auto aRadTranslation = aAlphaCoordinate * aRadTranslationPointA + aBetaCoordinate * aRadTranslationPointB + aGammaCoordinate * aRadTranslationPointC; + auto aRadScaling = aAlphaCoordinate * aRadScalingPointA + aBetaCoordinate * aRadScalingPointB + aGammaCoordinate * aRadScalingPointC; + + // compute formula of bilinear interpolation + auto aBilinearValueTri = FormalBilinTri_Formula(aVObs, TriangleDisplacement_NbObs, aXTri, aYTri); + + // Take into account radiometry in minimisation process + auto aRadiometryValueTri = aRadScaling * aIntensityImPre + aRadTranslation; + + // residual is simply the difference between values in before image and estimated value in new image + return {aRadiometryValueTri - aBilinearValueTri}; + } + }; + + class cTriangleDeformationTrRad + { + public: + cTriangleDeformationTrRad() + { + } + + static const std::vector VNamesUnknowns() { return std::vector{"GeomTrXPointA", "GeomTrYPointA", "RadTranslationPointA", "RadScalingPointA", + "GeomTrXPointB", "GeomTrYPointB", "RadTranslationPointB", "RadScalingPointB", + "GeomTrXPointC", "GeomTrYPointC", "RadTranslationPointC", "RadScalingPointC"}; } + static const std::vector VNamesObs() + { + return Append( + std::vector{"PixelCoordinatesX", "PixelCoordinatesY", "AlphaCoordPixel", "BetaCoordPixel", "GammaCoordPixel", "IntensityImPre"}, + FormalBilinIm2D_NameObs("T") // 6 obs for bilinear interpol of Im + ); + } + + std::string FormulaName() const { return "TriangleDeformationTrRad"; } + + template + static std::vector formula( + const std::vector &aVUnk, + const std::vector &aVObs) + { + + // extract observation + const auto &aXCoordinate = aVObs[0]; + const auto &aYCoordinate = aVObs[1]; + const auto &aAlphaCoordinate = aVObs[2]; + const auto &aBetaCoordinate = aVObs[3]; + const auto &aGammaCoordinate = aVObs[4]; + const auto &aIntensityImPre = aVObs[5]; + + // extract unknowns + const auto &aGeomTrXPointA = aVUnk[0]; + const auto &aGeomTrYPointA = aVUnk[1]; + const auto &aRadTranslationPointA = aVUnk[2]; + const auto &aRadScalingPointA = aVUnk[3]; + const auto &aGeomTrXPointB = aVUnk[4]; + const auto &aGeomTrYPointB = aVUnk[5]; + const auto &aRadTranslationPointB = aVUnk[6]; + const auto &aRadScalingPointB = aVUnk[7]; + const auto &aGeomTrXPointC = aVUnk[8]; + const auto &aGeomTrYPointC = aVUnk[9]; + const auto &aRadTranslationPointC = aVUnk[10]; + const auto &aRadScalingPointC = aVUnk[11]; + + // Apply barycentric translation formula to coordinates + auto aXTri = aXCoordinate + aAlphaCoordinate * aGeomTrXPointA + aBetaCoordinate * aGeomTrXPointB + aGammaCoordinate * aGeomTrXPointC; + auto aYTri = aYCoordinate + aAlphaCoordinate * aGeomTrYPointA + aBetaCoordinate * aGeomTrYPointB + aGammaCoordinate * aGeomTrYPointC; + + // Apply barycentric interpolation to radiometric factors + auto aRadTranslation = aAlphaCoordinate * aRadTranslationPointA + aBetaCoordinate * aRadTranslationPointB + aGammaCoordinate * aRadTranslationPointC; + auto aRadScaling = aAlphaCoordinate * aRadScalingPointA + aBetaCoordinate * aRadScalingPointB + aGammaCoordinate * aRadScalingPointC; + + // compute formula of bilinear interpolation + auto aBilinearValueTri = FormalBilinTri_Formula(aVObs, TriangleDisplacement_NbObs, aXTri, aYTri); + + // Take into account radiometry in minimisation process + auto aRadiometryValueTri = aRadScaling * aBilinearValueTri + aRadTranslation; + + // residual is simply the difference between values in before image and estimated value in new image + return {aIntensityImPre - aRadiometryValueTri}; + } + }; + + class cTriangleDeformationTranslation + { + public: + cTriangleDeformationTranslation() + { + } + + static const std::vector VNamesUnknowns() { return Append(std::vector{"GeomTrXPointA", "GeomTrYPointA"}, + std::vector{"GeomTrXPointB", "GeomTrYPointB"}, + std::vector{"GeomTrXPointC", "GeomTrYPointC"}); } + static const std::vector VNamesObs() + { + return Append( + std::vector{"PixelCoordinatesX", "PixelCoordinatesY", "AlphaCoordPixel", "BetaCoordPixel", "GammaCoordPixel", "IntensityImPre"}, + FormalBilinIm2D_NameObs("T") // 6 obs for bilinear interpol of Image + ); + } + + std::string FormulaName() const { return "TriangleDeformationTranslation"; } + + template + static std::vector formula( + const std::vector &aVUnk, + const std::vector &aVObs) + { + // extract observation + const auto &aXCoordinate = aVObs[0]; + const auto &aYCoordinate = aVObs[1]; + const auto &aAlphaCoordinate = aVObs[2]; + const auto &aBetaCoordinate = aVObs[3]; + const auto &aGammaCoordinate = aVObs[4]; + const auto &aIntensityImPre = aVObs[5]; + + // extract unknowns + const auto &aGeomTrXPointA = aVUnk[0]; + const auto &aGeomTrYPointA = aVUnk[1]; + const auto &aGeomTrXPointB = aVUnk[2]; + const auto &aGeomTrYPointB = aVUnk[3]; + const auto &aGeomTrXPointC = aVUnk[4]; + const auto &aGeomTrYPointC = aVUnk[5]; + + auto aXTri = aXCoordinate + aAlphaCoordinate * aGeomTrXPointA + aBetaCoordinate * aGeomTrXPointB + aGammaCoordinate * aGeomTrXPointC; + auto aYTri = aYCoordinate + aAlphaCoordinate * aGeomTrYPointA + aBetaCoordinate * aGeomTrYPointB + aGammaCoordinate * aGeomTrYPointC; + + // compute formula of bilinear interpolation + auto aEstimatedTranslatedValueTri = FormalBilinTri_Formula(aVObs, TriangleDisplacement_NbObs, aXTri, aYTri); + + // residual is simply the difference between values in before image and estimated value in new image. + return {aIntensityImPre - aEstimatedTranslatedValueTri}; + } + }; + + class cTriangleDeformationRadiometry + { + public: + cTriangleDeformationRadiometry() + { + } + + static const std::vector VNamesUnknowns() { return Append(std::vector{"RadTranslationPointA", "RadScalingPointA"}, + std::vector{"RadTranslationPointB", "RadScalingPointB"}, + std::vector{"RadTranslationPointC", "RadScalingPointC"}); } + + static const std::vector VNamesObs() + { + return Append( + std::vector{"PixelCoordinatesX", "PixelCoordinatesY", "AlphaCoordPixel", "BetaCoordPixel", "GammaCoordPixel", "IntensityImPre"}, + FormalBilinIm2D_NameObs("T") // 6 obs for bilinear interpol of Im + ); + } + + std::string FormulaName() const { return "TriangleDeformationRadiometry"; } + + template + static std::vector formula( + const std::vector &aVUnk, + const std::vector &aVObs) + { + // extract observation + const auto &aXCoordinate = aVObs[0]; + const auto &aYCoordinate = aVObs[1]; + const auto &aAlphaCoordinate = aVObs[2]; + const auto &aBetaCoordinate = aVObs[3]; + const auto &aGammaCoordinate = aVObs[4]; + const auto &aIntensityImPre = aVObs[5]; + + // extract unknowns + const auto &aRadTranslationPointA = aVUnk[0]; + const auto &aRadScalingPointA = aVUnk[1]; + const auto &aRadTranslationPointB = aVUnk[2]; + const auto &aRadScalingPointB = aVUnk[3]; + const auto &aRadTranslationPointC = aVUnk[4]; + const auto &aRadScalingPointC = aVUnk[5]; + + // Apply barycentric interpolation to radiometric factors + auto aRadTranslation = aAlphaCoordinate * aRadTranslationPointA + aBetaCoordinate * aRadTranslationPointB + aGammaCoordinate * aRadTranslationPointC; + auto aRadScaling = aAlphaCoordinate * aRadScalingPointA + aBetaCoordinate * aRadScalingPointB + aGammaCoordinate * aRadScalingPointC; + + // compute formula of bilinear interpolation + auto aBilinearValueTri = FormalBilinTri_Formula(aVObs, TriangleDisplacement_NbObs, aXCoordinate, aYCoordinate); + + // Take into account radiometry in minimisation process + auto aRadiometryValueTri = aRadScaling * aIntensityImPre + aRadTranslation; + + // residual is simply the difference between values in before image and estimated value in new image + return {aRadiometryValueTri - aBilinearValueTri}; + } + }; + +class cTriangleDeformationRad + { + public: + cTriangleDeformationRad() + { + } + + static const std::vector VNamesUnknowns() { return Append(std::vector{"RadTranslationPointA", "RadScalingPointA"}, + std::vector{"RadTranslationPointB", "RadScalingPointB"}, + std::vector{"RadTranslationPointC", "RadScalingPointC"}); } + + static const std::vector VNamesObs() + { + return Append( + std::vector{"PixelCoordinatesX", "PixelCoordinatesY", "AlphaCoordPixel", "BetaCoordPixel", "GammaCoordPixel", "IntensityImPre"}, + FormalBilinIm2D_NameObs("T") // 6 obs for bilinear interpol of Im + ); + } + + std::string FormulaName() const { return "TriangleDeformationRad"; } + + template + static std::vector formula( + const std::vector &aVUnk, + const std::vector &aVObs) + { + // extract observation + const auto &aXCoordinate = aVObs[0]; + const auto &aYCoordinate = aVObs[1]; + const auto &aAlphaCoordinate = aVObs[2]; + const auto &aBetaCoordinate = aVObs[3]; + const auto &aGammaCoordinate = aVObs[4]; + const auto &aIntensityImPre = aVObs[5]; + + // extract unknowns + const auto &aRadTranslationPointA = aVUnk[0]; + const auto &aRadScalingPointA = aVUnk[1]; + const auto &aRadTranslationPointB = aVUnk[2]; + const auto &aRadScalingPointB = aVUnk[3]; + const auto &aRadTranslationPointC = aVUnk[4]; + const auto &aRadScalingPointC = aVUnk[5]; + + // Apply barycentric interpolation to radiometric factors + auto aRadTranslation = aAlphaCoordinate * aRadTranslationPointA + aBetaCoordinate * aRadTranslationPointB + aGammaCoordinate * aRadTranslationPointC; + auto aRadScaling = aAlphaCoordinate * aRadScalingPointA + aBetaCoordinate * aRadScalingPointB + aGammaCoordinate * aRadScalingPointC; + + // SymbPrint(aRadTranslation, "RadTranslation"); + // SymbPrint(aRadScaling, "RadScaling"); + + // compute formula of bilinear interpolation + auto aBilinearValueTri = FormalBilinTri_Formula(aVObs, TriangleDisplacement_NbObs, aXCoordinate, aYCoordinate); + + // Take into account radiometry in minimisation process + auto aRadiometryValueTri = aRadScaling * aBilinearValueTri + aRadTranslation; + + // residual is simply the difference between values in before image and estimated value in new image + return {aIntensityImPre - aRadiometryValueTri}; + } + }; + +}; // MMVII + +#endif // _FORMULAS_TRIANGLES_DEFORM_H_ \ No newline at end of file diff --git a/MMVII/src/SymbDerGen/GenerateCodes.cpp b/MMVII/src/SymbDerGen/GenerateCodes.cpp index dbe20097c3..a44c2dd3f0 100755 --- a/MMVII/src/SymbDerGen/GenerateCodes.cpp +++ b/MMVII/src/SymbDerGen/GenerateCodes.cpp @@ -2,6 +2,7 @@ #include "SymbDer/SymbolicDerivatives.h" #include "SymbDer/SymbDer_GenNameAlloc.h" #include "Formulas_ImagesDeform.h" +#include "Formulas_TrianglesDeform.h" #include "Formulas_CamStenope.h" #include "Formulas_Geom2D.h" #include "Formulas_Radiom.h" @@ -308,6 +309,31 @@ cCalculator * EqDeformImAffinity(bool WithDerive,int aSzBuf) return StdAllocCalc(NameFormula(cDeformImAffinity(),WithDerive),aSzBuf); } +cCalculator *EqDeformTri(bool WithDerive, int aSzBuf) +{ + return StdAllocCalc(NameFormula(cTriangleDeformation(), WithDerive), aSzBuf); +} + +cCalculator *EqDeformTriTrRad(bool WithDerive, int aSzBuf) +{ + return StdAllocCalc(NameFormula(cTriangleDeformationTrRad(), WithDerive), aSzBuf); +} + +cCalculator *EqDeformTriTranslation(bool WithDerive, int aSzBuf) +{ + return StdAllocCalc(NameFormula(cTriangleDeformationTranslation(), WithDerive), aSzBuf); +} + +cCalculator *EqDeformTriRadiometry(bool WithDerive, int aSzBuf) +{ + return StdAllocCalc(NameFormula(cTriangleDeformationRadiometry(), WithDerive), aSzBuf); +} + +cCalculator *EqDeformTriRad(bool WithDerive, int aSzBuf) +{ + return StdAllocCalc(NameFormula(cTriangleDeformationRad(), WithDerive), aSzBuf); +} + // dist3d // Cons distance template cCalculator * TplEqDist3D(bool WithDerive,int aSzBuf) @@ -814,6 +840,12 @@ int cAppliGenCode::Exe() for (const auto IsLinearGrad : {true,false}) GenCodesFormula((tREAL8*)nullptr,cDeformImHomotethy(IsLinearGrad) ,WithDer); + GenCodesFormula((tREAL8 *)nullptr, cTriangleDeformation(), WithDer); + GenCodesFormula((tREAL8 *)nullptr, cTriangleDeformationTrRad(), WithDer); + GenCodesFormula((tREAL8 *)nullptr, cTriangleDeformationTranslation(), WithDer); + GenCodesFormula((tREAL8 *)nullptr, cTriangleDeformationRadiometry(), WithDer); + GenCodesFormula((tREAL8 *)nullptr, cTriangleDeformationRad(), WithDer); + // =============== CODE FOR RADIOMETRY =========================================