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 =========================================