From 9e696680371d3120c321b51f79f2ea52f62c0d9a Mon Sep 17 00:00:00 2001 From: JoLeBelge Date: Wed, 21 Feb 2018 16:44:14 +0100 Subject: [PATCH 1/9] rien de special --- include/general/photogram.h | 3 - include/private/files.h | 2 +- src/CBinaires/mm3d.cpp | 6 + src/RegTIRVIS/Keypoint.cpp | 6 + src/RegTIRVIS/Keypoint.h | 1 + src/RegTIRVIS/RegTIRVIS.cpp | 393 ++++---- src/TpMMPD/DIDRO/cAppuis2homol.cpp | 93 ++ src/TpMMPD/DIDRO/cAppuis2homol.h | 16 + src/TpMMPD/DIDRO/cblendingorthos.cpp | 11 + src/TpMMPD/DIDRO/cblendingorthos.h | 42 + src/TpMMPD/DIDRO/cero_appli.cpp | 6 +- src/TpMMPD/DIDRO/cero_modelonepaire.cpp | 226 ++++- src/TpMMPD/DIDRO/cero_modelonepaire.h | 13 +- src/TpMMPD/DIDRO/cimgeo.cpp | 46 + src/TpMMPD/DIDRO/cimgeo.h | 2 + src/TpMMPD/DIDRO/cimprovehomogtirvis.cpp | 7 + src/TpMMPD/DIDRO/cimprovehomogtirvis.h | 12 + src/TpMMPD/DIDRO/cransac_2dline.cpp | 90 +- src/TpMMPD/DIDRO/cransac_2dline.h | 10 +- src/TpMMPD/DIDRO/divers.cpp | 1146 +++++++++++++++++++++- src/TpMMPD/Sources.cmake | 1 + 21 files changed, 1900 insertions(+), 232 deletions(-) create mode 100644 src/TpMMPD/DIDRO/cAppuis2homol.cpp create mode 100644 src/TpMMPD/DIDRO/cAppuis2homol.h create mode 100644 src/TpMMPD/DIDRO/cblendingorthos.cpp create mode 100644 src/TpMMPD/DIDRO/cblendingorthos.h create mode 100644 src/TpMMPD/DIDRO/cimprovehomogtirvis.cpp create mode 100644 src/TpMMPD/DIDRO/cimprovehomogtirvis.h diff --git a/include/general/photogram.h b/include/general/photogram.h index d8cd2db6c8..87ee0df40f 100644 --- a/include/general/photogram.h +++ b/include/general/photogram.h @@ -863,9 +863,6 @@ class cElHomographie : public cElMap2D virtual void InitFromParams(const std::vector &aSol); virtual std::vector Params() const; - - - bool HasNan() const; Pt2dr Direct (const Pt2dr & aP) const; diff --git a/include/private/files.h b/include/private/files.h index 8529b461b4..241990b6a9 100644 --- a/include/private/files.h +++ b/include/private/files.h @@ -193,7 +193,7 @@ class ELISE_fp ~ELISE_fp(); ELISE_fp(eModeBinTxt ModeBin=eTxtOnPremierLigne); - ELISE_fp(const char *,mode_open,bool svp = false,eModeBinTxt =eTxtOnPremierLigne); + ELISE_fp(const char *,mode_open=READ,bool svp = false,eModeBinTxt =eTxtOnPremierLigne); U_INT1 read_U_INT1(); U_INT2 read_U_INT2(); diff --git a/src/CBinaires/mm3d.cpp b/src/CBinaires/mm3d.cpp index 93b852991f..de758be777 100644 --- a/src/CBinaires/mm3d.cpp +++ b/src/CBinaires/mm3d.cpp @@ -637,10 +637,13 @@ extern int T2V_main(int argc,char ** argv); extern int Tapioca_IDR_main(int argc,char ** argv); extern int resizeImg_main(int argc,char ** argv); extern int resizeHomol_main(int argc,char ** argv); +extern int VarioCamTo8Bits_main(int argc,char ** argv); // test je jo extern int main_test(int argc,char ** argv); +extern int main_test2(int argc,char ** argv); extern int main_ero(int argc,char ** argv); extern int main_ascii2tif(int argc,char ** argv); +extern int GCP2Hom_main(int argc,char ** argv); #if (ELISE_UNIX) extern int DocEx_Introanalyse_main(int,char **); #endif @@ -934,9 +937,12 @@ const std::vector & TestLibAvailableCommands() aRes.push_back(cMMCom("GeoSud",ServiceGeoSud_GeoSud_main,"")); aRes.push_back(cMMCom("Surf",ServiceGeoSud_Surf_main,"")); aRes.push_back(cMMCom("ImageRectification",ImageRectification,"Rectify individual aerial images, ground is assumed to be a plane")); + aRes.push_back(cMMCom("ThermicVarioCamTo8Bits",VarioCamTo8Bits_main,"Convert 16 bits Variocam thermic images to 8 bits")); aRes.push_back(cMMCom("jo_FFH",FilterFileHom_main,"filtrer un fichier de paire d'image")); aRes.push_back(cMMCom("jo_T2V",T2V_main,"appliquer une homographie a un ensemble d'im thermique pour Reg avec images visibles")); aRes.push_back(cMMCom("jo_test",main_test,"test function for didro project")); + aRes.push_back(cMMCom("jo_test2",main_test2,"test function for didro project")); + aRes.push_back(cMMCom("GCP2Hom",GCP2Hom_main,"Convert GCP 2D measures in homol file")); aRes.push_back(cMMCom("TapiocaIDR",Tapioca_IDR_main,"Utiliser Tapioca avec des Images de Résolution Différente (effectue un resample des images)")); aRes.push_back(cMMCom("ResizeImg",resizeImg_main,"Resize image in order to reach a specific image width")); aRes.push_back(cMMCom("ResizeHomol",resizeHomol_main,"Resize Homol pack")); diff --git a/src/RegTIRVIS/Keypoint.cpp b/src/RegTIRVIS/Keypoint.cpp index 8b52b937df..684d69020a 100644 --- a/src/RegTIRVIS/Keypoint.cpp +++ b/src/RegTIRVIS/Keypoint.cpp @@ -23,3 +23,9 @@ KeyPoint::KeyPoint(float x, float y, float size, float angle, float response) KeyPoint::KeyPoint(){} KeyPoint::~KeyPoint(){} + +void KeyPoint::undist(CamStenope * aCalib) +{ + Pt2dr ptUndist=aCalib->DistDirecte(Pt2dr(m_Point.x,m_Point.y)); + m_Point=Pt2df(ptUndist.x,ptUndist.y); +} diff --git a/src/RegTIRVIS/Keypoint.h b/src/RegTIRVIS/Keypoint.h index cc067da288..5597041ec2 100644 --- a/src/RegTIRVIS/Keypoint.h +++ b/src/RegTIRVIS/Keypoint.h @@ -31,6 +31,7 @@ class KeyPoint void setSize(float sz) {m_size=sz;} void setAngle(float angle){m_angle=angle;} void setResponse(float response){m_response=response;} + void undist(CamStenope * aCalib); //getters Pt2df getPoint(){return m_Point;} float getAngle(){return m_angle;} diff --git a/src/RegTIRVIS/RegTIRVIS.cpp b/src/RegTIRVIS/RegTIRVIS.cpp index e4a8c391ff..d517bf4c5a 100644 --- a/src/RegTIRVIS/RegTIRVIS.cpp +++ b/src/RegTIRVIS/RegTIRVIS.cpp @@ -1,17 +1,4 @@ -#include "StdAfx.h" -#include -#include "Image.h" -#include "msd.h" -#include "Keypoint.h" -#include "lab_header.h" -#include "../../uti_image/Digeo/Digeo.h" -#include "DescriptorExtractor.h" -#include "Arbre.h" - -#include "../../uti_image/Digeo/DigeoPoint.h" - -#define distMax 0.75 -#define rayMax 0.5 +#include "cregtirvis.h" /*******************************************************/ @@ -122,12 +109,27 @@ void Resizeim(Im2D & im, Im2D & Out, Pt2dr Newsize) } } } + +void cAppliRegTIRVIS::mkDirPastis(std::string aSH, std::string aImName, std::string aDir) +{ + std::string DirPastis=aDir +"/Homol" + aSH + "/Pastis" + aImName ; + if (!ELISE_fp::exist_file(DirPastis)) + { + //create directory Pastis + ELISE_fp::MkDir(DirPastis); + } +} + + + + //===============================================================================// -// Estimate Homograhy via multiple iterations on keypoints which have not been matched during +// Estimate Homograhy via multiple iterations on keypoints which have not been matched bu Ann during // the correspondence search phase //===============================================================================// -void EnrichKps(std::vector< KeyPoint > Kpsfrom, ArbreKD * Tree, cElHomographie &Homog, int NbIter) + +void cAppliRegTIRVIS::EnrichKps(std::vector< KeyPoint > Kpsfrom, ArbreKD * Tree, cElHomographie &Homog, int NbIter,int ImPairKey) { ElSTDNS set< pair > Voisins ; // where to put nearest neighbours @@ -155,11 +157,26 @@ void EnrichKps(std::vector< KeyPoint > Kpsfrom, ArbreKD * Tree, cElHomographie & std::cout<<" Computed Homography at "<Nb of Tie point used: "<Ecart "<Quality "<If Ok=1: "<Assoc1To2(aKeyAsocHom, ThermalImages[ImPairKey], VisualImages[ImPairKey],true); + HomologousPts.StdPutInFile(aHomolFile); + HomologousPts.SelfSwap(); + aHomolFile= mDirTest+ mICNM->Assoc1To2(aKeyAsocHom, VisualImages[ImPairKey],ThermalImages[ImPairKey] ,true); + HomologousPts.StdPutInFile(aHomolFile); } + + //delete mImgIdx; } @@ -622,104 +639,21 @@ void ParseHomol(string MasterImage, std::vector< cCpleString> ImCpls,std::vector } } -//===============================================================================// -/* OrientedImage class */ -//===============================================================================// -class Orient_Image -{ - public: - Orient_Image - ( - std::string aOriIn, - std::string aName, - cInterfChantierNameManipulateur * aICNM - ); - std::string getName(){return mName;} - CamStenope * getCam(){return mCam;} - std::string getOrifileName(){return mOriFileName;} - CamStenope * mCam; - - protected: - - std::string mName; - std::string mOriFileName; - }; - -Orient_Image::Orient_Image - ( std::string aOriIn, - std::string aName, - cInterfChantierNameManipulateur * aICNM): - mName(aName),mOriFileName(aOriIn+"Orientation-"+mName+".xml") -{ - mCam=CamOrientGenFromFile(mOriFileName,aICNM); -} -//===============================================================================// -/* Main Function */ -//===============================================================================// -int RegTIRVIS_main( int argc, char ** argv ) -{ - - // Call Elise Librairy We will be using because classes that handle - // These transformations are already computed - - std::string Image_Pattern, Oris_VIS_dir,TestDataDIRpat, PlyFileIn; - std::string Option; - bool aTif=false; - - /************************************************************************/ - //Initilialize ElInitArgMain which as i understood captures arguments entered - //by the operator - /**********************************************************************/ - ElInitArgMain - ( - argc,argv, - LArgMain() << EAMC(TestDataDIRpat ,"Set of images to be used to estimate the homography",eSAM_IsPatFile) - << EAMC(Image_Pattern," Images pattern : Dir+ Name: prefixes Thermal: TIR, Visual: VIS", eSAM_IsPatFile) - << EAMC(Oris_VIS_dir," Orientation files for visual images",eSAM_IsExistDirOri) - << EAMC(PlyFileIn, " Ply file that will be used as a 3d model",eSAM_IsExistFile) - ,LArgMain() << EAM(Option,"Option",true," FOR NOW, DON'T DO ANYTHING") - ); - - -/**********************************************************************/ - if (MMVisualMode) return EXIT_SUCCESS; - - - - -/*@@@@@@@@@@@@@@@@@@@@@@@@@ Train DATA IMAGES PROCESSING @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ - - - - - //Take the testdata directory to compute the homography that - //will be used as a predictor for keypoint matching - - // MakeFileDirCompl(TestDataDIRpat); - //std::cout<<"Test data Directory "< Test data directory: "< Pattern image test: "< aSetImTest = *(aICNM0->Get(aPatImTst)); - std::cout<<"==============> Param Chantier manipulateur set \n"; + const std::vector aSetImTest = *(mICNM->Get(mPatImTest)); /*===============================================================================*/ /* Transform all images to 8 bits: Lab and wallis on 8 bits */ /*===============================================================================*/ - std::size_t found = aPatImTst.find_last_of("."); - string ext = aPatImTst.substr(found+1); + found = mPatImTest.find_last_of("."); + ext = mPatImTest.substr(found+1); cout<<"==============> Test images Extension : "< cmd; std::cout< ThermalImages, VisualImages; - std::size_t found_Vis; - std::size_t found_Tir; - std::size_t found_commonPattern; + for (unsigned int j=0; j > Voisins ; // where to put nearest neighbours + + box=Box2dr(Pt2dr(0,0), Pt2dr(VisualSize.x,VisualSize.y) ); + /*********************************************************** */ @@ -836,9 +745,9 @@ int RegTIRVIS_main( int argc, char ** argv ) cElRegex rgxxx("vis_(.*).tif",10); std::string aNameMatch; bool rgxMatch,rgxxMatch,rgxxxMatch; - rgxMatch=rgx.Match(VisualImages[i]); - rgxxMatch=rgxx.Match(VisualImages[i]); - rgxxxMatch=rgxxx.Match(VisualImages[i]); + rgxMatch=rgx.Match(mDirTest+VisualImages[i]); + rgxxMatch=rgxx.Match(mDirTest+VisualImages[i]); + rgxxxMatch=rgxxx.Match(mDirTest+VisualImages[i]); if (rgxMatch||rgxxMatch||rgxxxMatch) { if(rgxMatch) aNameMatch=rgx.KIemeExprPar(1); @@ -851,14 +760,14 @@ int RegTIRVIS_main( int argc, char ** argv ) if (found_commonPattern!=std::string::npos) { //Process the Termal Image: MSD --> Lab --> Wallis --> SIFT Descriptor - Tiff_Im ImageV=Tiff_Im::UnivConvStd(VisualImages[i]); + Tiff_Im ImageV=Tiff_Im::UnivConvStd(mDirTest+VisualImages[i]); /**************************************************************************/ // Check if MSD keypoints have already been computed: No need to do that again std::vector KpsV; - string directory, file; - SplitDirAndFile(directory,file,VisualImages[i]); - string filenameV=Path + file + ".txt"; + //string directory, file; + // SplitDirAndFile(directory,file,VisualImages[i]); + string filenameV=mDirTest + "KpsTEST/" + VisualImages[i] + ".txt"; if (DoesFileExist(filenameV.c_str())) @@ -919,12 +828,12 @@ int RegTIRVIS_main( int argc, char ** argv ) //Process the Termal Image: MSD --> Lab --> Wallis --> SIFT Descriptor - Tiff_Im ImageTh=Tiff_Im::UnivConvStd(ThermalImages[i]); + Tiff_Im ImageTh=Tiff_Im::UnivConvStd(mDirTest+ThermalImages[i]); // Check if MSD keypoints have already been computed: No need to do that again - SplitDirAndFile(directory,file,ThermalImages[i]); + //SplitDirAndFile(directory,file,ThermalImages[i]); std::vector KpsTh; - string filenameTh=Path + file + ".txt"; + string filenameTh=mDirTest + "KpsTEST/"+ ThermalImages[i] + ".txt"; if (DoesFileExist(filenameTh.c_str())) { @@ -959,6 +868,11 @@ int RegTIRVIS_main( int argc, char ** argv ) vector ListTh; aKp=KpsTh.begin(); + + // test jo ; si je corrige la position des Kps avec le modèle de distorion du capteur, est-ce que ça améliore le résultat? + //CamStenope * aCalib; + //if(EAMIsInit(&aCalibName)) aCalib=Std_Cal_From_File(aCalibName); + for (;aKp!=KpsTh.end();++aKp) { DigeoPoint DP; @@ -970,22 +884,51 @@ int RegTIRVIS_main( int argc, char ** argv ) SIFTTh.normalize_and_truncate(descriptor); //std::cout<<"Descp "<getPoint() << "\n"; + /* if(EAMIsInit(&aCalibName)) + {aKp->undist(aCalib); + DP.x=aKp->getPoint().x; + DP.y=aKp->getPoint().y; + } + //std::cout << " Keypoint position after undistort ; " << aKp->getPoint() << "\n"; +*/ ListTh.push_back(DP); } + + + DigeoPoint::writeDigeoFile("DigeoTH.txt",ListTh); + + + std::string aSH="-Init"; + std::string aKeyAsocHom = "NKS-Assoc-CplIm2Hom@"+ aSH +"@txt"; + std::string aHomolFileT2V,aHomolFileV2T, aCmd; + // the success of Ann require that the output directory must exist ; we create them is is do not exist already + mkDirPastis(aSH,VisualImages[i],mDirTest); + mkDirPastis(aSH,ThermalImages[i],mDirTest); + list cmd; - string aCmd=MM3DStr + " Ann "+ std::string("-ratio 0.9") + std::string(" DigeoTH.txt") + std::string(" DigeoV.txt") + std::string(" Matches.txt"); + + aHomolFileV2T= mDirTest+ mICNM->Assoc1To2(aKeyAsocHom, VisualImages[i],ThermalImages[i] ,true); + aCmd=MM3DStr + " Ann "+ std::string("-ratio 0.9") + std::string(" DigeoV.txt") + std::string(" DigeoTH.txt ") + std::string(aHomolFileV2T); + cmd.push_back(aCmd); + if (mDebug) std::cout << aCmd << "\n"; + aHomolFileT2V = mDirTest+ mICNM->Assoc1To2(aKeyAsocHom, ThermalImages[i], VisualImages[i],true); + aCmd=MM3DStr + " Ann "+ std::string("-ratio 0.9") + std::string(" DigeoTH.txt") + std::string(" DigeoV.txt ") + std::string(aHomolFileT2V); + if (mDebug) std::cout << aCmd << "\n"; cmd.push_back(aCmd); cEl_GPAO::DoComInParal(cmd); + //Compute a Robust Homography out of the resulting Matching FIle ElPackHomologue HomologousPts; - bool Exist= ELISE_fp::exist_file("Matches.txt"); + bool Exist= ELISE_fp::exist_file(aHomolFileV2T); if (Exist) { - HomologousPts= ElPackHomologue::FromFile("Matches.txt"); + HomologousPts= ElPackHomologue::FromFile(aHomolFileT2V); } cElComposHomographie Ix(0,0,0); cElComposHomographie Iy(0,0,0); @@ -994,31 +937,120 @@ int RegTIRVIS_main( int argc, char ** argv ) double anEcart,aQuality; bool Ok; - Hout=H2estimate.RobustInit(anEcart,&aQuality,HomologousPts,Ok,50,80.0,2000); + H2estimate=H2estimate.RobustInit(anEcart,&aQuality,HomologousPts,Ok,50,80.0,2000); + mH=&H2estimate; std::cout<< " =================> Initial Homography <=====================\n"; - Hout.Show(); + mH->Show(); + std::cout<<"========= >Nb of Tie point used: "<Ecart "<Quality "<If Ok= 1 " <Inverse().ToXmlGen(),aNameFileHomog); + ArbreKD * ArbreV= new ArbreKD(Pt_of_Point, box, KpsV.size(), 1.0); for (uint i=0; iinsert(pair(i, Pt2dr(KpsV.at(i).getPoint().x, KpsV.at(i).getPoint().y))); } - // Call Enrich Keypoints to use the homography as predictor - EnrichKps(KpsTh,ArbreV,Hout,5); + // Call Enrich Keypoints to use the homography as predictor + EnrichKps(KpsTh,ArbreV,*mH,5,i); + // EnrichKps refine the homography + aNameFileHomog=mDirTest+"homographyTIRVIS-final1.xml"; + + + std::cout << "Save refined Homography TIR to VIS to " << aNameFileHomog << "\n"; + MakeFileXML(mH->Inverse().ToXmlGen(),aNameFileHomog); delete ArbreV; } } +} -/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ +void cAppliRegTIRVIS::initMSD() +{ + //===============================================================================// + // Instantiate the detector MSD // + /*===============================================================================*/ + + + msd.setPatchRadius(3); + msd.setSearchAreaRadius(5); + msd.setNMSRadius(5); + msd.setNMSScaleRadius(0); + msd.setThSaliency(0.02); + msd.setKNN(10); + msd.setScaleFactor(1.25f); + msd.setNScales(-1); + msd.setComputeOrientation(true); + msd.setCircularWindow(true); + msd.setRefinedKP(false); +} +//===============================================================================// +/* Main Function */ +//===============================================================================// + +cAppliRegTIRVIS::cAppliRegTIRVIS( int argc, char ** argv ): + mH(0) +{ + std::string Image_Pattern, Oris_VIS_dir,TestDataDIRpat, PlyFileIn, aHomogFile; + mDebug=1; + bool aTif=false; + + ElInitArgMain + ( + argc,argv, + LArgMain() + << EAMC(Image_Pattern," Images pattern : Dir+ Name: prefixes Thermal: TIR, Visual: VIS", eSAM_IsPatFile) + << EAMC(Oris_VIS_dir," Orientation files for visual images",eSAM_IsExistDirOri) + << EAMC(PlyFileIn, " Ply file that will be used as a 3d model",eSAM_IsExistFile) + ,LArgMain() << EAM(TestDataDIRpat,"ImPat4Homogr",true,"Set of images to be used to estimate the homography",eSAM_IsPatFile) + << EAM(aHomogFile,"Homog", true,"The homography, if provided, will not required to estimate it ",eSAM_IsPatFile) + << EAM(mDebug,"Debug",true,"Print more info to the terminal in order to facilitate debugging") + //<< EAM(aCalibName,"CalibTIR",true,"Calib file for thir camera, test") + ); + + if (!MMVisualMode) { + + + + + + // load of compute homography + + if (EAMIsInit(&aHomogFile)) + { + mICNM=cInterfChantierNameManipulateur::BasicAlloc("./"); + // grosse bricole + ELISE_fp tmp = ELISE_fp(aHomogFile.c_str()); + cElHomographie H=cElHomographie::read(tmp); + mH=&H; + } else { + + + if(EAMIsInit(&TestDataDIRpat) ) + { + SplitDirAndFile(mDirTest,mPatImTest,TestDataDIRpat); + std::cout<<"==============> Test data directory: "< Pattern image test: "<Terrain=>image2 Orient_Image ImV1(Oris_VIS_dir,ImCpls.at(i).N1(),aICNM); Orient_Image ImV2(Oris_VIS_dir,ImCpls.at(i).N2(),aICNM); + if (mDebug) std::cout << "image couple " << ImCpls.at(i).N1() << " and " << ImCpls.at(i).N2() << "\n"; //Apply Image1 ==> Terrain operation having into consideration a depth image @@ -1534,9 +1570,9 @@ for (uint i=0; i MSD Points "< MSD Points "<Terrain=>image2 Orient_Image ImV1(Oris_VIS_dir,ImCpls.at(i).N1(),aICNM); Orient_Image ImV2(Oris_VIS_dir,ImCpls.at(i).N2(),aICNM); - + if (mDebug) std::cout << "image couple " << ImCpls.at(i).N1() << " and " << ImCpls.at(i).N2() << "\n"; //Apply Image1 ==> Terrain operation having into consideration a depth image // from ZBufferRaster directory @@ -2440,9 +2475,9 @@ if (!ELISE_fp::IsDirectory(HomolfileInter)) // Homolfile is not created // Apply the computed homography to the thermal images Master image /****************MSD MSD MSD MSD MSD ******************/ - Kps1H=NewSetKpAfterHomog(Kps1,Hout); + Kps1H=NewSetKpAfterHomog(Kps1,*mH); /***************SIFT SIFT SIFT SIFT SIFT *************/ - Kps1HS=NewSetKpAfterHomog(Kps1S,Hout); + Kps1HS=NewSetKpAfterHomog(Kps1S,*mH); std::cout<<"===========> MSD Points "< MSD Points "<Assoc1To2(aKeyAsocHom, mIm1, mIm2,true); + if (mDebug) std::cout << "generate " << aHomolFile << "\n"; + + aPackHom.StdPutInFile(aHomolFile); + aPackHom.SelfSwap(); + aHomolFile= mICNM->Assoc1To2(aKeyAsocHom, mIm2, mIm1,true); + if (mDebug) std::cout << "generate " << aHomolFile << "\n"; + + aPackHom.StdPutInFile(aHomolFile); + std::cout << "Finished, " << count << " manual seasing of GCP have been converted in homol format\n"; + std::string aKH("NB"); + if (mExpTxt) aKH="NT"; + + std::cout << "Launch SEL for visualisation: SEL ./ " << mIm1 << " " << mIm2 << " KH=" << aKH << " SH=" << mSH << "\n"; + } else { std::cout << "I haven't found couple of 2D measure for these images pairs, sorry \n";} +} + +int GCP2Hom_main(int argc,char ** argv) +{ + cAppuis2Homol(argc,argv); + return EXIT_SUCCESS; +} diff --git a/src/TpMMPD/DIDRO/cAppuis2homol.h b/src/TpMMPD/DIDRO/cAppuis2homol.h new file mode 100644 index 0000000000..c041cec390 --- /dev/null +++ b/src/TpMMPD/DIDRO/cAppuis2homol.h @@ -0,0 +1,16 @@ +#ifndef CAPPUIS2HOMOL_H +#define CAPPUIS2HOMOL_H +#include "StdAfx.h" + +// goal: convert 2D measure of appuis (from saisie appuis tools ) to homol format +class cAppuis2Homol +{ +public: + cAppuis2Homol(int argc, char** argv); +private: + cInterfChantierNameManipulateur * mICNM; + bool mDebug,mExpTxt; + std::string mIm1,mIm2,mHomPackOut, mSH, m2DMesFileName; +}; + +#endif // CAPPUIS2HOMOL_H diff --git a/src/TpMMPD/DIDRO/cblendingorthos.cpp b/src/TpMMPD/DIDRO/cblendingorthos.cpp new file mode 100644 index 0000000000..661b749259 --- /dev/null +++ b/src/TpMMPD/DIDRO/cblendingorthos.cpp @@ -0,0 +1,11 @@ +#include "cblendingorthos.h" + +cBlendingOrthos::cBlendingOrthos() +{ + +} + +cAppli_BlendingOrthos::cAppli_cBlendingOrthos() +{ + +} diff --git a/src/TpMMPD/DIDRO/cblendingorthos.h b/src/TpMMPD/DIDRO/cblendingorthos.h new file mode 100644 index 0000000000..c862e26b55 --- /dev/null +++ b/src/TpMMPD/DIDRO/cblendingorthos.h @@ -0,0 +1,42 @@ +#ifndef CBLENDINGORTHOS_H +#define CBLENDINGORTHOS_H +#include "cero_modelonepaire.h" + +class cAppli_BlendingOrthos +{ +public: + cAppli_BlendingOrthos(int argc, char** argv); + +private: + cInterfChantierNameManipulateur * mICNM; + std::string mFullName; + std::string mDir; + std::string mPat; + std::list mLFile; + std::vector mVOs; // vecteur d'orthos + +}; + +// classe ortho contenant 1) la radiométrie 2) la zone ou on effectue le blending + +class cBlendingOrthos +{ +public: + cBlendingOrthos(cAppli_BlendingOrthos & anAppli,const std::string & aName); + Im2D_INT1 * AreaInMosaic(){return &mAreaInMosaic;} + Im2D_INT1 * blendingArea(){return &mBlendingArea;} + Im2D_REAL4 * rad(){return &mIm;} // radiometry + // receive num of label and label map and feathering distance in pixels, distance d5711 op morpho + void SetAreaInMosaicAndBlendingZ(const mLabelMap); + +private: + Im2D_INT1 mAreaInMosaic; + Im2D_INT1 mBlendingArea; + Im2D_REAL4 mIm; + std::string mName; + int mInd; + cAppli_BlendingOrthos & mAppli; + +}; + +#endif // CBLENDINGORTHOS_H diff --git a/src/TpMMPD/DIDRO/cero_appli.cpp b/src/TpMMPD/DIDRO/cero_appli.cpp index e23fcff3ed..43cdf394b7 100644 --- a/src/TpMMPD/DIDRO/cero_appli.cpp +++ b/src/TpMMPD/DIDRO/cero_appli.cpp @@ -242,14 +242,12 @@ void cERO_Appli::computeModel4EveryPairs() + mDir +cple.N2() + " Debug=" + ToString(mDebug) + " Dir=" + mDirOut - + " W1=1 W2=1" // poid des images, test avec poid egal toutes les images + + " W1=1 W2=1 WIncid=1" // poid des images, test avec poid egal toutes les images ; aLCom.push_back(aCom); if (mDebug) std::cout << aCom << "\n"; } - cEl_GPAO::DoComInParal(aLCom); - } // charger les modeles et calculer mod moyen @@ -290,7 +288,7 @@ void cERO_Appli::loadEROSmodel4OneIm(std::string aNameOrt) // LSQ matching cLSQ_2dline LSQ_mod(&aCpleRad,&aVPond); - LSQ_mod.adjustModel(); + LSQ_mod.adjustModelL2(); if (mDebug) LSQ_mod.affiche(); // sauve le modèle dans la liste des modèles diff --git a/src/TpMMPD/DIDRO/cero_modelonepaire.cpp b/src/TpMMPD/DIDRO/cero_modelonepaire.cpp index b084f12a41..e78012ea03 100644 --- a/src/TpMMPD/DIDRO/cero_modelonepaire.cpp +++ b/src/TpMMPD/DIDRO/cero_modelonepaire.cpp @@ -13,6 +13,7 @@ cERO_ModelOnePaire::cERO_ModelOnePaire(int argc, char** argv):mO1(0),mO2(0) mTmpDirEROS=""; mNbCplRadio=0; mDir="./"; + mPondIncid=0; ElInitArgMain ( @@ -24,6 +25,7 @@ cERO_ModelOnePaire::cERO_ModelOnePaire(int argc, char** argv):mO1(0),mO2(0) << EAM(mPoid2,"W2",true,"Weight we give to Ortho 2. Default 1, means we want to correct Ortho1 in order to have Ortho2 radiometry. If W1=W2, the model will try to egalize Ortho1 radiometry in order to match the mean radiometry of Ortho1 and Ortho2.") << EAM(mDebug,"Debug",true,"If true, write dozens of useless intermediates results, both in the terminal and in your disk.") << EAM(mTmpDirEROS,"Dir",true,"Directory where to save the results") + << EAM(mPondIncid,"WIncid",true,"Weighting radiometric observation with incidence value? def false") ); std::cout << "launch Radiometric Egalization of 2 Orthos \n"; @@ -56,11 +58,21 @@ cERO_ModelOnePaire::cERO_ModelOnePaire(int argc, char** argv):mO1(0),mO2(0) // Initialise les 2 orthos mO1 = new cImGeo(mDir+mNameOr1); mO2 = new cImGeo(mDir+mNameOr2); + // Determine la zone de recouvrement entre les 2 orthos mBoxOverlapTerrain=mO1->overlapBox(mO2); // clip les 2 ortho sur cette box terrain mO1clip = mO1->clipImTer(mBoxOverlapTerrain); mO2clip = mO2->clipImTer(mBoxOverlapTerrain); + + if (mDebug | mPondIncid) + { + mO1->loadIncid(); + mO2->loadIncid(); + mI1clip = mO1->clipIncidTer(mBoxOverlapTerrain); + mI2clip = mO2->clipIncidTer(mBoxOverlapTerrain); + } + /* les 2 Im2D peuvent avoir un pixels de décalage, en y dans mon exemple, pk? pas bon ça std::cout << "mO1clip size " << mO1clip.sz() << "\n"; std::cout << "mO2clip size " << mO2clip.sz() << "\n"; */ @@ -73,17 +85,21 @@ cERO_ModelOnePaire::cERO_ModelOnePaire(int argc, char** argv):mO1(0),mO2(0) ); mO2clip = tmp; + + if (mDebug){ saveIm2D(& mO1clip, mTmpDirERO+"Im1clipped.tif"); saveIm2D(& mO2clip, mTmpDirERO+"Im2clipped.tif"); + } - // on utilise ransac pour, sur chacune des tuiles, garder 2 couples radiométriques exempt de bruit et représentatif de la tuile + // on utilise ransac pour, sur chacune des tuiles, garder 2 couples radiométriques exempt de bruit et représentatif de la tuile ainsi que le produit de la moyenne de l'incidence des images ransacOnEachTiles(); + // maintenant qu'on a fait le tour de chaque tuile:dalle, on calcul les modele Global - ransacOn2Orthos(); + L1On2Orthos(); // appliquer les modèle sur les 2 orthos isolées et effectue le mosaiquage afin de pouvoir inspecter facilement le résultat if (mDebug) applyRE(); @@ -181,12 +197,14 @@ void cERO_ModelOnePaire::applyRE() // pour le moment; applique la correction rad Im2D_REAL4 im1(mO1->applyRE(mModelER1)); // je sauve le résultats saveIm2D(&im1,filename); + mO1->writeTFW(filename); // image 2 filename=mTmpDirERO+"Im2_egalized.tif"; Im2D_REAL4 im2(mO2->applyRE(mModelER2)); // je sauve le résultats saveIm2D(&im2,filename); + mO2->writeTFW(filename); // mosaique des 2 images egalizées Box2dr boxMosaic=mO1->boxEnglob(mO2); @@ -195,8 +213,7 @@ void cERO_ModelOnePaire::applyRE() // pour le moment; applique la correction rad Im2D_REAL4 mosaic=mO1->box2Im(boxMosaic); Pt2dr aCorner=Pt2dr(boxMosaic._p0.x,boxMosaic._p1.y); // xmin, ymax; Pt2di tr1=mO1->computeTrans(aCorner), tr2=mO2->computeTrans(aCorner); - mO1->loadIncid(); - mO2->loadIncid(); + //saveIm2D(mO1->Incid(),mTmpDirERO+"Incid2.tif"); // copy o1 in the mosaic @@ -251,6 +268,8 @@ void cERO_ModelOnePaire::ransacOnEachTiles() int loopX=mO1clip.sz().x/mSzTuile.x; int loopY=mO1clip.sz().y/mSzTuile.y; + int aNbObs(0); + for (int i(0) ; i * aObs, c2DLineModel aMod, int aQuartile,bool Y) { // if bool Y=true, mean we are interested in y value of the pt2d, problem is that sort method work only on first element of Pt2d. poor and dangerous Solution is to flip the observations @@ -378,6 +506,18 @@ Pt2dr cERO_ModelOnePaire::predQuantile(std::vector * aObs, c2DLineModel a } } +Pt2dr cERO_ModelOnePaire::predQuantile(std::map * aObsMap, c2DLineModel aMod, int aQuartile,bool Y) +{ + // extract the vector of Pt2dr from the map, not interrested in the key and i am not capable of manipulating map + std::vector aObs; + + for (std::map::iterator it=aObsMap->begin(); it!=aObsMap->end(); ++it) + aObs.push_back(it->second); + + return predQuantile(&aObs, aMod,aQuartile,Y); +} + +/* void cERO_ModelOnePaire::ransacOn2Orthos() { // on utilise les valeurs de pondération W1 et W2 pour savoir à quoi sert le modèle calculé @@ -409,6 +549,7 @@ void cERO_ModelOnePaire::ransacOn2Orthos() // modifie le vecteur d'information: inverse le couple radiométrique std::vector aObs; for (unsigned int i(0); i on veut corriger Im2 mais pas Im1 + // W1=W2 : on veut corriger et Im1 et Im2 afin d'obtenir une radiométrie moyenne à ces deux images. + + //plusieurs situation . 1iere, un des deux poids est 0 + if (mPoid1==0 || mPoid2==0) + { + // modèle im 1 vers im 2 + if (mPoid1==0){ + + if (mDebug) std::cout << "Radiom Im1 toward radiom Im2 \n"; + cLSQ_2dline LSQ_mod(&mObsGlob,&mPondGlob); + LSQ_mod.adjustModelL1(); + if (mDebug) LSQ_mod.affiche(); + mModelER1=LSQ_mod.getModel(); + // modèle "unité" + mModelER2=c2DLineModel(); + } + // modèle im 2 vers im 1 + if (mPoid2==0){ + if (mDebug) std::cout << "Radiom Im2 toward radiom Im1 \n"; + // modifie le vecteur d'information: inverse le couple radiométrique + for (std::map::iterator it=mObsGlob.begin(); it!=mObsGlob.end(); ++it) + { + it->second= it->second.yx(); + } + cLSQ_2dline LSQ_mod(&mObsGlob,&mPondGlob); + LSQ_mod.adjustModelL1(); + if (mDebug) LSQ_mod.affiche(); + mModelER2=LSQ_mod.getModel(); + mModelER1=c2DLineModel(); + } + } + if (mPoid1!=0 && mPoid2!=0) + { + if (mDebug) std::cout << "Correct Im1 and Im2 with weighting " << mPoid1 << " and " << mPoid2 << " respectively \n"; + // création des vecteurs d'observation pour les 2 images + std::map aObs1,aObs2; + + for (auto & obs : mObsGlob) + { + aObs1[obs.first]=Pt2dr(obs.second.x,pond(obs.second)); + aObs2[obs.first]=Pt2dr(obs.second.y,pond(obs.second)); + } + + // ajustement des modèles + + cLSQ_2dline LSQ_mod(&aObs1,&mPondGlob); + LSQ_mod.adjustModelL1(); + if (mDebug) LSQ_mod.affiche(); + mModelER1=LSQ_mod.getModel(); + + cLSQ_2dline LSQ_mod2(&aObs2,&mPondGlob); + LSQ_mod2.adjustModelL1(); + if (mDebug) LSQ_mod2.affiche(); + mModelER2=LSQ_mod2.getModel(); + } +} + + int main_ero(int argc,char ** argv) { diff --git a/src/TpMMPD/DIDRO/cero_modelonepaire.h b/src/TpMMPD/DIDRO/cero_modelonepaire.h index 26a6e66946..b056f873f7 100644 --- a/src/TpMMPD/DIDRO/cero_modelonepaire.h +++ b/src/TpMMPD/DIDRO/cero_modelonepaire.h @@ -24,13 +24,13 @@ class cERO_ModelOnePaire{ cInterfChantierNameManipulateur * mICNM; cImGeo * mO1, * mO2; // Ortho 1 et Ortho 2 Im2D_REAL4 mO1clip,mO2clip; // la zone de recouvrement - //cImGeo * mPC1,mPC2,mIncid1,mIncid2; + Im2D_REAL4 mI1clip,mI2clip; // incidence std::string mNameOr1,mNameOr2,mNameOr1AndDir,mNameOr2AndDir,mDir; bool mAddCst; // ajout de la constante a dans le modèle, default = true // les couples d'observation radiométrique sont stoquée dans un vecteur std::vector mObsDallage; // observation pour une tuille donnée - std::vector mObsGlob; // observation globale pour le modèle d'égalisation radiom à l'échelles globales - //std::vector mPond; // pour stoquer une valeur de pondération pour chacune des observations des couples radiométriques + std::map mObsGlob; // observation globale pour le modèle d'égalisation radiom à l'échelles globales + std::map mPondGlob; // pour stoquer une valeur de pondération pour chacune des observations des couples radiométriques Pt2di mSzTuile; Box2dr mBoxOverlapTerrain; Box2di mBoxOverO1, mBoxOverO2; // box en pixel de la zone de recouvrement entre les 2 orthos @@ -40,15 +40,20 @@ class cERO_ModelOnePaire{ c2DLineModel mModelER1,mModelER2; int mPoid1,mPoid2; std::string mTmpDirEROS,mTmpDirERO;// un directory global pour tout les modèles de chaques paires d'images, un local pour la paire d'ortho + bool mPondIncid; + void ransacOnEachTiles(); void ransacOn2Orthos(); + void L1On2Orthos(); void saveModels(); void boxTer2Tfw(); // sauver une geom terrain comme tfw double pond(Pt2dr aPt); - Pt2dr predQuantile(std::vector * aObs, c2DLineModel aMod, int quantile=1, bool Y=false);// on donne une liste d'observation et un modèle , renvoie le couple X(quantile n),Prediction(X) + Pt2dr predQuantile(std::map *aObsMap, c2DLineModel aMod, int quantile=1, bool Y=false);// on donne une liste d'observation et un modèle , renvoie le couple X(quantile n),Prediction(X) + Pt2dr predQuantile(std::vector *aObs, c2DLineModel aMod, int quantile=1, bool Y=false);// on donne une liste d'observation et un modèle , renvoie le couple X(quantile n),Prediction(X) + }; #endif // CERO_MODELONEPAIRE_H diff --git a/src/TpMMPD/DIDRO/cimgeo.cpp b/src/TpMMPD/DIDRO/cimgeo.cpp index 38282eff25..27d54fa486 100644 --- a/src/TpMMPD/DIDRO/cimgeo.cpp +++ b/src/TpMMPD/DIDRO/cimgeo.cpp @@ -355,6 +355,35 @@ Im2D_REAL4 cImGeo::clipImTer(Box2dr aBox) return im; } + +Im2D_REAL4 cImGeo::clipIncidTer(Box2dr aBox) +{ + + // compute box in pixel + int U(0), V(0); + U = (aBox.sz().x)/GSD(); + V = (aBox.sz().y)/GSD(); + // initialise resulting image + Im2D_REAL4 im(U,V); + + // inversion min max Y car en geom image l'axe des y pointe dans le sens contraire + Pt2dr aPMin(aBox._p0.x,aBox._p1.y), aPMax(aBox._p1.x,aBox._p0.y); + + if (containTer(aBox._p0) && containTer(aBox._p1)) + { + im=clipIncidPix(XY2UV(aPMin),XY2UV(aPMax)); + + } + else + { + std::cout << "Dans l'appel à la fonction cImGeo::clipIncidTer(Box2dr aBox); nous vous signalons que la box n'est pas comprise dans l'image \n"; + } + + return im; +} + + + Im2D_REAL4 cImGeo::clipImPix(Pt2di aMin,Pt2di aMax) { // initialise l'image retour @@ -371,6 +400,23 @@ Im2D_REAL4 cImGeo::clipImPix(Pt2di aMin,Pt2di aMax) return im; } + +Im2D_REAL4 cImGeo::clipIncidPix(Pt2di aMin,Pt2di aMax) +{ + // initialise l'image retour + Im2D_REAL4 im(aMax.x-aMin.x,aMax.y-aMin.y); + // translate l'imageGeo et clip en meme temps + Pt2di aTr(aMin.x,aMin.y); + + ELISE_COPY + ( + Incid()->all_pts(), + trans(Incid()->in(1),aTr), + im.oclip() + ); + return im; +} + bool cImGeo::containTer(Pt2dr pt) { return (pt.x<=mXmax && pt.x>=mXmin && pt.y<= mYmax && pt.y>= mYmin); diff --git a/src/TpMMPD/DIDRO/cimgeo.h b/src/TpMMPD/DIDRO/cimgeo.h index 8355404000..b0b39967f0 100644 --- a/src/TpMMPD/DIDRO/cimgeo.h +++ b/src/TpMMPD/DIDRO/cimgeo.h @@ -47,6 +47,8 @@ class cImGeo Im2D_REAL4 clipImPix(Pt2di aMin,Pt2di aMax); // clip l'image avec une box pixel Im2D_REAL4 clipImTer(Pt2dr aMin,Pt2dr aMax); // clip l'image avec une box terrain Im2D_REAL4 clipImTer(Box2dr aBox); // clip l'image avec une box terrain + Im2D_REAL4 clipIncidTer(Box2dr aBox); + Im2D_REAL4 clipIncidPix(Pt2di aMin,Pt2di aMax); Box2dr overlapBox(cImGeo * aIm2); // renvoie la box terrain du recouvrement des 2 images bool containTer(Pt2dr pt); diff --git a/src/TpMMPD/DIDRO/cimprovehomogtirvis.cpp b/src/TpMMPD/DIDRO/cimprovehomogtirvis.cpp new file mode 100644 index 0000000000..92e17b5b9b --- /dev/null +++ b/src/TpMMPD/DIDRO/cimprovehomogtirvis.cpp @@ -0,0 +1,7 @@ +#include "cimprovehomogtirvis.h" + +cImproveHomogTirVis::cImproveHomogTirVis() +{ + +} + diff --git a/src/TpMMPD/DIDRO/cimprovehomogtirvis.h b/src/TpMMPD/DIDRO/cimprovehomogtirvis.h new file mode 100644 index 0000000000..88214e55bd --- /dev/null +++ b/src/TpMMPD/DIDRO/cimprovehomogtirvis.h @@ -0,0 +1,12 @@ +#ifndef CIMPROVEHOMOGTIRVIS_H +#define CIMPROVEHOMOGTIRVIS_H + + +class cImproveHomogTirVis +{ +public: + cImproveHomogTirVis(); + +}; + +#endif // CIMPROVEHOMOGTIRVIS_H diff --git a/src/TpMMPD/DIDRO/cransac_2dline.cpp b/src/TpMMPD/DIDRO/cransac_2dline.cpp index 4b5063999c..146b7c288d 100644 --- a/src/TpMMPD/DIDRO/cransac_2dline.cpp +++ b/src/TpMMPD/DIDRO/cransac_2dline.cpp @@ -27,6 +27,17 @@ cRansac_2dline::cRansac_2dline(std::vector * aObs, ModelForm model, int a } } + +cRansac_2dline::cRansac_2dline(std::map * aObsMap, ModelForm model, int aNbInt) +{ + std::vector aObs; + + for (std::map::iterator it=aObsMap->begin(); it!=aObsMap->end(); ++it) + aObs.push_back(it->second); + cRansac_2dline(&aObs,model,aNbInt); +} + + void cRansac_2dline::affiche() { // affiche les résumés du modèle @@ -177,29 +188,41 @@ void c2DLineModel::computeCout(std::vector * aObs,std::vector * a cLSQ_2dline::cLSQ_2dline(std::vector * aObs,ModelForm model): mModelForm(model), -mObs(aObs), +mPond(1), mOk(0) { -// mettre toutes les pondérations à 1 -for (unsigned int i(0) ;i< mObs->size(); i++) -{ - mPond.push_back(1); -} + + for (unsigned int i(0) ; isize();i++) + { + pObs.push_back(std::make_pair(&mPond,&aObs->at(i))); + } } cLSQ_2dline::cLSQ_2dline(std::vector * aObs,std::vector * aPond,ModelForm model): mModelForm(model), -mObs(aObs), -mPond(*aPond), mOk(0) { // check that observation and ponderation have the same length - if(mPond.size()!=mObs->size()) + if(aPond->size()!=aObs->size()) {std::cout <<"Warning: Ponderation vector have not the same length than observations.\n";} + if((aPond->size()==0 )|(aObs->size()==0)) + {std::cout <<"Warning: Observations or ponderation vector is empty.\n";} + for (unsigned int i(0) ; isize();i++) + { + pObs.push_back(std::make_pair(&aPond->at(i),&aObs->at(i))); + } +} +cLSQ_2dline::cLSQ_2dline(std::map * aObsMap,std::map * aPondMap, ModelForm model) +{ + for (std::map::iterator it=aObsMap->begin(); it!=aObsMap->end(); ++it) + { + pObs.push_back(std::make_pair(&aPondMap->at(it->first),&it->second)); + } } -cLSQ_2dline::cLSQ_2dline():mObs(0),mOk(0) + +cLSQ_2dline::cLSQ_2dline():mPond(1),mOk(0) { } @@ -210,18 +233,17 @@ cLSQ_2dline::~cLSQ_2dline() //delete mObs; } -void cLSQ_2dline::adjustModel() +void cLSQ_2dline::adjustModelL2() { // Create L2SysSurResol to solve least square equation with 2 unknowns L2SysSurResol aSys(2); //For Each radiometric Couples, add the observations - int it(0); - for(auto & rc : *mObs){ - double aPds[2]={rc.x,1}; - double poids=mPond.at(it); - aSys.AddEquation(poids,aPds,rc.y); - it++; + for(auto & pair_pointers : pObs){ + + double aFormLin[2]={pair_pointers.second->x,1}; // b*x+a + double poids=*pair_pointers.first; + aSys.AddEquation(poids,aFormLin,pair_pointers.second->y); } Im1D_REAL8 aSol = aSys.GSSR_Solve(&mOk); @@ -232,10 +254,40 @@ void cLSQ_2dline::adjustModel() //std::cout << "solution trouvée , b =" << aData[0] << " and a " << aData[1] << " \n"; mModel=c2DLineModel(aData[1],aData[0]); } else { - std::cout << "adjustment of a 2D line by LSQ failed\n"; + std::cout << "adjustment of a 2D line by LSQ L2 failed\n"; } } + +void cLSQ_2dline::adjustModelL1() +{ + SystLinSurResolu aSys(2,pObs.size()); + + //For Each radiometric Couples, add the observations + + for(auto & pair_pointers : pObs){ + + double aFormLin[2]={pair_pointers.second->x,1}; // b*x+a + int poids=*pair_pointers.first; + //std::cout << "add equation to systlinsurREsolu : " << pair_pointers.second->x << " + a = " <y << ", weighting of " << poids << "\n"; + aSys.PushEquation(aFormLin,pair_pointers.second->y,poids); + } + + Im1D_REAL8 aSol = aSys.L1Solve(); + mOk=1; + + if (mOk) + { + double* aData = aSol.data(); + //std::cout << "solution trouvée , b =" << aData[0] << " and a " << aData[1] << " \n"; + mModel=c2DLineModel(aData[1],aData[0]); + } else { + std::cout << "adjustment of a 2D line by LSQ L1 failed\n"; + } + +} + + void cLSQ_2dline::affiche() { // affiche les résumés du modèle @@ -243,7 +295,7 @@ void cLSQ_2dline::affiche() // test si l'ajustement a déjà eu lieu if(mOk) { - std::cout << "Observations :" << mObs->size() << "\n"; + std::cout << "Observations :" << pObs.size() << "\n"; std::cout << "a + b * x = y a=" << mModel.getA() << ", b=" << mModel.getB() << ".\n"; } else std::cout << "l'ajustement as échoué ou n'as pas encore été effectué.\n"; } diff --git a/src/TpMMPD/DIDRO/cransac_2dline.h b/src/TpMMPD/DIDRO/cransac_2dline.h index 06efca6ad2..bbaf283df6 100644 --- a/src/TpMMPD/DIDRO/cransac_2dline.h +++ b/src/TpMMPD/DIDRO/cransac_2dline.h @@ -44,6 +44,7 @@ class cRansac_2dline cRansac_2dline(); // création de l'objet à partir des observations, définition de la forme du modèle, d'un vecteur de pondération cRansac_2dline(std::vector * aObs,ModelForm model,int aNbInt); + cRansac_2dline(std::map * aObsMap,ModelForm model,int aNbInt); //cRansac_2dline(std::vector aObs,ModelForm model, std::vector aPonderation); ~cRansac_2dline(); @@ -78,9 +79,11 @@ class cLSQ_2dline // création de l'objet à partir des observations, définition de la forme du modèle cLSQ_2dline(std::vector * aObs,ModelForm model=a_plus_bx); cLSQ_2dline(std::vector * aObs,std::vector * aPond,ModelForm model=a_plus_bx); + cLSQ_2dline(std::map * aObsMap, std::map * aPondMap, ModelForm model=a_plus_bx); ~cLSQ_2dline(); - void adjustModel(); + void adjustModelL2();//L2 + void adjustModelL1();//L1 // renvoie une prédiction du modèle ajusté précédement double predict(double aX) {return mModel.predict(aX);} // donne des info sur le modele @@ -90,8 +93,9 @@ class cLSQ_2dline private: c2DLineModel mModel; ModelForm mModelForm; - std::vector * mObs; - std::vector mPond; + //std::vector * mObs; + std::vector> pObs; // pointers to ponderation and observation + double mPond; bool mOk; }; diff --git a/src/TpMMPD/DIDRO/divers.cpp b/src/TpMMPD/DIDRO/divers.cpp index c100e053e3..0749f27587 100644 --- a/src/TpMMPD/DIDRO/divers.cpp +++ b/src/TpMMPD/DIDRO/divers.cpp @@ -44,6 +44,985 @@ Header-MicMac-eLiSe-25/06/2007*/ extern int RegTIRVIS_main(int , char **); +// we wish to improve coregistration between 2 orthos +class cCoreg2Ortho +{ + public: + std::string mDir; + cCoreg2Ortho(int argc,char ** argv); + private: + cImGeo * mO1; + cImGeo * mO2; + Im2D_REAL4 mO1clip,mO2clip; + std::string mNameO1, mNameO2, mNameMapOut; + Box2dr mBoxOverlapTerrain; + +}; + + + +cCoreg2Ortho::cCoreg2Ortho(int argc,char ** argv) +{ + + ElInitArgMain + ( + argc,argv, + LArgMain() << EAMC(mNameO1,"Name Ortho master", eSAM_IsExistFile) + << EAMC(mNameO2,"Name Ortho slave",eSAM_IsExistFile), + LArgMain() << EAM(mNameMapOut,"Out",true, "Name of resulting map") + ); + + if (!MMVisualMode) + { + + mDir="./"; + mNameMapOut=mNameO2 +"2"+ mNameO1; + cInterfChantierNameManipulateur * aICNM = cInterfChantierNameManipulateur::BasicAlloc(mDir); + + if (ELISE_fp::exist_file(mNameO1) & ELISE_fp::exist_file(mNameO2)) + { + // open orthos + // Initialise les 2 orthos + mO1 = new cImGeo(mDir+mNameO1); + mO2 = new cImGeo(mDir+mNameO2); + + // je teste non pas la coreg mais le blending par gaussian distance weighting + Im2D_REAL4 I1=mO1->toRAM(); + Im2D_REAL4 I2=mO2->toRAM(); + Box2dr boxMosaic=mO1->boxEnglob(mO2); + Im2D_REAL4 mosaic=mO1->box2Im(boxMosaic); + Pt2dr aCorner=Pt2dr(boxMosaic._p0.x,boxMosaic._p1.y); // xmin, ymax; + Pt2di tr1=mO1->computeTrans(aCorner), tr2=mO2->computeTrans(aCorner); + + // right now i consider that nadir position is at the center of the image + // express nadir position in pixel in geometry of mosaic + Pt2di N1((Pt2di(I1.sz()/Pt2di(2,2))+tr1)),N2((Pt2di(I2.sz()/Pt2di(2,2))+tr2)); + double constLambda(0.5); + + // test la fonction Distance de chamfer sur une image binaire bidonne + Im2D_Bits<1> mIBin(I1.sz().x,I1.sz().y,1); + std::string aName("IBin1.tif"); + ELISE_COPY( + mIBin.all_pts(), + mIBin.in(), + Tiff_Im(aName.c_str(), + mIBin.sz(), + GenIm::int1, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + ELISE_COPY( + disc(Pt2dr(128,128),100), + 0, + mIBin.out() + ); + + ELISE_COPY( + disc(Pt2dr(500,500),100), + 0, + mIBin.out() + ); + // la fonction chamfer fonctionne sur une image binaire et va calculer les distance à partir des pixels qui ont une valeur de 0. + + // la distance maximale est de 255 car la fonction chamfer + + + Im2D_U_INT1 aD(I1.sz().x,I1.sz().y); + ELISE_COPY(aD.all_pts(),mIBin.in(),aD.out()); + + Chamfer::d32.im_dist(aD); + + std::string aName3("chamferTestd32.tif"); + ELISE_COPY( + aD.all_pts(), + aD.in(), + Tiff_Im(aName3.c_str(), + aD.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + ELISE_COPY( + disc(Pt2dr(700,700),200), + 0, + aD.out() + ); + + + Chamfer::d5711.im_dist(aD); + + std::string aName4("chamferTestd5711.tif"); + ELISE_COPY( + aD.all_pts(), + aD.in(), + Tiff_Im(aName4.c_str(), + aD.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + /* + std::cout << "mosaic of size " << mosaic.sz() << ".\n"; + + for (unsigned int i(1) ; i overlapBox(mO2); + // clip les 2 ortho sur cette box terrain afin d'avoir des Im2D chargé en RAM + mO1clip = mO1->clipImTer(mBoxOverlapTerrain); + mO2clip = mO2->clipImTer(mBoxOverlapTerrain); + + std::string aOut1A("im1.tif"),aOut2A("im2.tif"); + ELISE_COPY( + mO2clip.all_pts(), + mO2clip.in(), + Tiff_Im(aOut2A.c_str(), + mO2clip.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + ELISE_COPY( + mO1clip.all_pts(), + mO1clip.in(), + Tiff_Im(aOut1A.c_str(), + mO1clip.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + + Pt2dr aCorner=Pt2dr(mBoxOverlapTerrain._p0.x,mBoxOverlapTerrain._p1.y); // xmin, ymax; + Pt2di transO1Tobox = mO1->computeTrans(aCorner); + Pt2di transO2Tobox = mO2->computeTrans(aCorner); + // Pt2di trans = mO1->computeTrans(mO2); + + Pt2di sz(25,25); + unsigned int pasX=mO1clip.sz().x/10; + unsigned int pasY=mO1clip.sz().y/10; + std::cout << "step x " << pasX << ", step Y " <sz.x) & (j*pasY>sz.y) & ((i*pasX+sz.x) mIms; // key=label + std::map mIm2Ds; + std::map mChamferDist; + std::map mBlendingArea; + std::map mMosaicArea; + std::string mNameMosaicOut, mFullDir; + int mDist; + std::list mLFile; + cInterfChantierNameManipulateur * mICNM; + Box2dr mBoxOverlapTerrain; + double mLambda; +}; + +cFeatheringAndMosaicOrtho::cFeatheringAndMosaicOrtho(int argc,char ** argv) +{ + mDist=100; // distance chamfer de 100 pour le feathering/estompage + mLambda=0.2; + ElInitArgMain + ( + argc,argv, + LArgMain() << EAMC(mFullDir,"ortho pattern", eSAM_IsPatFile) + , + LArgMain() << EAM(mNameMosaicOut,"Out",true, "Name of resulting map") + << EAM(mDist,"Dist",true, "Distance for seamline feathering blending" ) + << EAM(mLambda,"Lambda",true, "lambda value for gaussian distance weighting, def 0.2" ) + ); + + if (!MMVisualMode) + { + + mDir="./"; + mNameMosaicOut="mosaicFeatheringTest1.tif"; + mICNM = cInterfChantierNameManipulateur::BasicAlloc(mDir); + mLFile = mICNM->StdGetListOfFile(mFullDir); + + cFileOriMnt MTD = StdGetFromPCP("MTDOrtho.xml",FileOriMnt); + Box2dr boxMosaic(Pt2dr(MTD.OriginePlani().x,MTD.OriginePlani().y+MTD.ResolutionPlani().y*MTD.NombrePixels().y),Pt2dr(MTD.OriginePlani().x+MTD.ResolutionPlani().x*MTD.NombrePixels().x,MTD.OriginePlani().y)); + Im2D_REAL4 mosaic(MTD.NombrePixels().x,MTD.NombrePixels().y); + Im2D_REAL4 NbIm(MTD.NombrePixels().x,MTD.NombrePixels().y); + Im2D_INT4 SumDist(MTD.NombrePixels().x,MTD.NombrePixels().y); + Im2D_INT4 SumDist2(MTD.NombrePixels().x,MTD.NombrePixels().y); + Im2D_REAL4 PondInterne(MTD.NombrePixels().x,MTD.NombrePixels().y); + Pt2dr aCorner=Pt2dr(boxMosaic._p0.x,boxMosaic._p1.y); // xmin, ymax; + Pt2di sz(MTD.NombrePixels().x,MTD.NombrePixels().y); + + Im2D_U_INT1 label=Im2D_U_INT1::FromFileStd("Label.tif"); + /* + std::string aName("testlabel.tif"); + ELISE_COPY( + label.all_pts(), + label.in(), + Tiff_Im(aName.c_str(), + label.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + */ + + unsigned int i(0); // i is the label of the image - and the key + for (auto &im : mLFile) + { + std::cout << "Image " << im <<": load and compute feathering buffer.\n"; + // open orthos + mIms[i]= new cImGeo(mDir+im); + mIm2Ds[i]= mIms[i]->toRAM(); + mChamferDist[i]=Im2D_INT2(mIms[i]->Im().sz().x,mIms[i]->Im().sz().y,1); + Pt2di sz(mIms[i]->Im().sz()); + Pt2di tr= -mIms[i]->computeTrans(aCorner); + + // la fonction chamfer fonctionne sur une image binaire et va calculer les distance à partir des pixels qui ont une valeur de 0. + // la distance maximale est de 255 car la fonction chamfer + + //detect seam line for this image + //1) translation of label to be in ortho geometry + Im2D_U_INT1 tmp(sz.x,sz.y,1); + ELISE_COPY(select(tmp.all_pts(), trans(label.in(),tr)!=(int)i), + //trans(label.in(),tr), + 0, + tmp.out() + ); + // set value to 0 for points that are not inside area of mosaicking for this image + // compute chamfer Distance + Chamfer::d32.im_dist(tmp); + // inverse value of distance (enveloppe interne à la seamline) + ELISE_COPY(mChamferDist[i].all_pts(),-tmp.in(),mChamferDist[i].out()); + + // je réinitialise tmp + ELISE_COPY(tmp.all_pts(),1,tmp.out()); + ELISE_COPY( + select(mChamferDist[i].all_pts(),mChamferDist[i].in()==-2),// distance ==-2 c'est les pixels Sur la ligne de raccord + 0, + tmp.out()); + // ça m'ennuie qu'il retourne une distance à par rapport au bord de l'image, comment annuler ce comportement? + Chamfer::d32.im_dist(tmp); + + ELISE_COPY( + select(mChamferDist[i].all_pts(),mChamferDist[i].in()==0), + tmp.in(), + mChamferDist[i].out()); + + // apply the hidden part masq + + //std::string aNameIm = mICNM->Assoc2To1("Key-Assoc-OpIm2PC@",im,true).first; + std::string aNamePC = "PC"+ im.substr(3,im.size()-2); + + Im2D_U_INT1 masq=Im2D_U_INT1::FromFileStd(aNamePC); + // apply the mask + ELISE_COPY( + select(mChamferDist[i].all_pts(),masq.in()==255), + 255, + mChamferDist[i].out()); + + std::string aNameTmp="TestChamfer" + std::to_string(i) + ".tif"; + ELISE_COPY( + mChamferDist[i].all_pts(), + mChamferDist[i].in(), + Tiff_Im(aNameTmp.c_str(), + mChamferDist[i].sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + // warning, effet de bord du au calcul de distance chamfer!!!!!!!!!!!!! + + // comptage du nombre d'image a utiliser pour le blending (geométrie mosaic) + + ELISE_COPY(select(NbIm.all_pts(),trans(mChamferDist[i].in(mDist+1),-tr)<=mDist), + //mosaic.in()+trans(mIm2Ds[i].in(0),tr)*lut_w.in()[-trans(mChamferDist[i].in(0) ,tr)], + NbIm.in(0)+1, + NbIm.out() + ); + + // Q? pourquoi je ne peux pas renseigner juste in() sans avoir une erreur genre BITMAP : out of domain while reading (RLE mode) + + + // somme des distances de chamber dans les enveloppes externes - pour gérer les cas de blending de 3 images + ELISE_COPY(select(SumDist.all_pts(),trans(mChamferDist[i].in_proj(),-tr)<=mDist & trans(mChamferDist[i].in_proj(),-tr)>0), + SumDist.in(0)+trans(mChamferDist[i].in(0),-tr), + SumDist.out() + ); + // somme des distances de chamber dans les enveloppes inter - également pour gérer les cas de blending de 3 images + ELISE_COPY(select(SumDist2.all_pts(),trans(mChamferDist[i].in_proj(),-tr)>=-mDist & trans(mChamferDist[i].in_proj(),-tr)<0), + SumDist2.in(0)+trans(mChamferDist[i].in(0),-tr), + SumDist2.out() + ); + + i++; + } + + + // je modifie la zone du blending de 3 images afin de 1) la réduire et 2) quel aie une forme plus adéquate + // non je ne peux pas faire ça sinon la zone NbIm==3 que je supprime sera considérée comme NbIm=2 alors que c'est faux + /* ELISE_COPY(select(NbIm.all_pts(),SumDist.in()>mDist), + 2, + NbIm.out() + ); + + Im2D_U_INT1 tmp(MTD.NombrePixels().x,MTD.NombrePixels().y,0); + ELISE_COPY(select(tmp.all_pts(),NbIm.in()==3), + 1, + tmp.out() + ); + + Chamfer::d32.im_dist(tmp); + + std::string aNameTmp("TestTmp.tif"); + ELISE_COPY( + tmp.all_pts(), + tmp.in(), + Tiff_Im(aNameTmp.c_str(), + tmp.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); +*/ + std::string aNameTmp("TestCDG.tif"); + + aNameTmp="TestNbIm.tif"; + + ELISE_COPY( + NbIm.all_pts(), + NbIm.in(), + Tiff_Im(aNameTmp.c_str(), + mosaic.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + + // détection du centre de gravité des zones de recouvements 3 images = point d'intersection des 3 images + Im2D_U_INT1 Ibin(MTD.NombrePixels().x,MTD.NombrePixels().y,0); + Im2D_U_INT1 chanferZTriple(sz.x,sz.y,0); + ELISE_COPY(select(Ibin.all_pts(),NbIm.in()==3), + 1, + Ibin.out() + ); + + U_INT1 ** d = Ibin.data(); + int count(0); + Neighbourhood V8=Neighbourhood::v8(); + + std::cout << "start detection of area of 3 ortho blending\n"; + for (INT x=0; x < Ibin.sz().x; x++) + { + for (INT y=0; y < Ibin.sz().y; y++) + { + if (d[y][x] == 1) + { + count++; + Liste_Pts_INT2 cc(2); + ELISE_COPY + ( + // flux: composantes connexes du point. + conc + ( + Pt2di(x,y), + Ibin.neigh_test_and_set(V8, 1, 0, 20) ), // on change la valeur des points sélectionné comme ça à la prochaine itération on ne sélectionne plus cette zone de composante connexe + 1, // valeur bidonne, c'est juste le flux que je sauve dans cc + cc + ); + + Pt2di cdg; + ELISE_COPY + ( + cc.all_pts(), + Virgule(FX,FY), + (cdg.sigma()) // sigma: somme des position X et Y + ); + cdg=cdg/cc.card(); // centre de grav: moyenne des positions + //std::cout << " points d'une ligne de raccord intersectant 3 orthos : num " << count << " position " << cdg <<"\n"; + // fin if (d[y][x] == 1) + // ok j'ai le pt d'intersection, qu'est ce que j'en fait? + // 1distance autour de ce point. cela concerne uniquement la zone de recouvrement triple soit cc + Im2D_U_INT1 chanfer(sz.x,sz.y,1); + ELISE_COPY(cc.all_pts(), + 0, + chanfer.out() + ); + Chamfer::d32.im_dist(chanfer); + // maintenant j'utilise cette distance en combinaison avec distance Env interne + + + + + } + } + } + + aNameTmp="TestSumDist.tif"; + ELISE_COPY( + SumDist.all_pts(), + SumDist.in(), + Tiff_Im(aNameTmp.c_str(), + SumDist.sz(), + GenIm::u_int1, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + aNameTmp="TestSumDistInterne.tif"; + ELISE_COPY( + SumDist2.all_pts(), + SumDist2.in(), + Tiff_Im(aNameTmp.c_str(), + SumDist2.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + // maintenant qu'on a toute les info on applique le blending + std::cout << "Start Blending\n"; + + // look up table de weight pour enveloppe interne + Im1D_REAL4 lut_w(mDist+1); + ELISE_COPY + ( + lut_w.all_pts(), + // FX/(double)mDist, + pow(0.5,pow((FX/((double)mDist-FX+1)),2*mLambda)), + lut_w.out() + ); + + // replace by 1 distance over the threshold of mDist + ELISE_COPY(select(lut_w.all_pts(),lut_w.in()>1),1,lut_w.out()); + + /* + for (unsigned int i(0); i=-mDist), + PondInterne.in()+ (1-(1/NbIm.in())) * lut_w.in()[mDist+SumDist2.in()], + PondInterne.out() + ); + + aNameTmp="TestPondInterne.tif"; + ELISE_COPY( + PondInterne.all_pts(), + PondInterne.in(), + Tiff_Im(aNameTmp.c_str(), + PondInterne.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + + for (unsigned int i(0); i < mIms.size();i++) + { + + std::cout << "Image" << i << "\n"; + + Pt2di tr= mIms[i]->computeTrans(aCorner); + + // enveloppe interne + +/* + // toute l'enveloppe interne, partie fixe pondérée seulement par le nombre d'image + ELISE_COPY(select(mosaic.all_pts(),trans(mChamferDist[i].in(1),tr)<=0 & NbIm.in()!=0), + // ok ça ca va : mosaic.in()+ trans(mIm2Ds[i].in(0),tr)*lut_w.in()[-trans(mChamferDist[i].in(0) ,tr)], + //trans(mIm2Ds[i].in(0),tr)*lut_w.in()[-trans(mChamferDist[i].in(0) ,tr)], fonctionne pas, croppe l'image + // mosaic.in()+ trans(mIm2Ds[i].in(0),tr) * (1-(NbIm.in()-1)/NbIm.in()), + + //mosaic.in()+ trans(mIm2Ds[i].in(0),tr) *(1-((NbIm.in()-1)*2.00)/(NbIm.in()*2.00)), + // PONDERATION ONLY + mosaic.in()+ (1-((NbIm.in()-1)*2.00)/(NbIm.in()*2.00)), + mosaic.out() + ); + + + // buffer dans l'enveloppe interne , je peux utiliser SumDist2 si je le garde au finish? non car sumdist interne a des valeurs pour les 2 coté de la ligne de raccord + ELISE_COPY(select(mosaic.all_pts(),trans(mChamferDist[i].in(0),tr)<=0 & trans(mChamferDist[i].in(0),tr)>=-mDist & NbIm.in()!=0), + //mosaic.in()+ trans(mIm2Ds[i].in(0),tr) * (1-(1/NbIm.in())) * lut_w.in()[-trans(mChamferDist[i].in(0) ,tr)], + // mosaic.in()+ (1-(1/NbIm.in())) * lut_w.in()[-trans(mChamferDist[i].in(0) ,tr)], + mosaic.in()+ (1-(2.00/(2.00*NbIm.in()))) * lut_w.in()[-trans(mChamferDist[i].in(0) ,tr)], + mosaic.out() + ); + + + // enveloppe externe + + + // ligne de raccord entre 2 images + ELISE_COPY(select(mosaic.all_pts(),trans(mChamferDist[i].in(0),tr)>0 & trans(mChamferDist[i].in(0),tr)<=mDist & NbIm.in()==2), + //mosaic.in()+ trans(mIm2Ds[i].in(0),tr)* (1-lut_w.in()[trans(mChamferDist[i].in(0) ,tr)]) * (1.00-(2.00/(2*NbIm.in()))) * 2.0*trans(mChamferDist[i].in(0),tr)/(2*SumDist.in(0)),// le 2.O/2 il empeche un bug débile + //mosaic.in()+ (1-lut_w.in()[trans(mChamferDist[i].in(0) ,tr)]) * (1.00-(2.00/(2*NbIm.in()))) * 2.0*trans(mChamferDist[i].in(0),tr)/(2*SumDist.in(0)), + // ponderation only + mosaic.in()+ (1-lut_w.in()[trans(mChamferDist[i].in(0) ,tr)]) * (1.00-(2.00/(2*NbIm.in()))), + mosaic.out() + ); + + + // raccord entre 3 images + ELISE_COPY(select(mosaic.all_pts(),trans(mChamferDist[i].in(0),tr)>0 & trans(mChamferDist[i].in(0),tr)<=mDist & NbIm.in()==3 & SumDist.in()>0), + // ponderation only + // mosaic.in()+ (1-((1-(2.00/(2.00*NbIm.in()))) * lut_w.in()[-SumDist2.in()])) * 2.0*trans(mChamferDist[i].in(0),tr)/(2.0*SumDist.in(0)), + mosaic.in()+ (1-(1-(2.00/(2.00*NbIm.in())) * lut_w.in()[-SumDist2.in()])) * 2.0*trans(mChamferDist[i].in(0),tr)/(2.0*SumDist.in(0)) , + mosaic.out() + ); + */ + // enveloppe interne + ELISE_COPY(select(mosaic.all_pts(),trans(mChamferDist[i].in(0),tr)<0), + mosaic.in()+ trans(mIm2Ds[i].in(0),tr)* PondInterne.in(), + //mosaic.in()+ PondInterne.in(), + mosaic.out() + ); + + // ligne de raccord entre 2 images = partage du gateau en 2 + ELISE_COPY(select(mosaic.all_pts(),trans(mChamferDist[i].in(0),tr)>0 & trans(mChamferDist[i].in(0),tr)<=mDist & NbIm.in()==2), + mosaic.in()+ trans(mIm2Ds[i].in(0),tr)* (1-PondInterne.in()), + //mosaic.in()+(1-PondInterne.in()), + mosaic.out() + ); + + + // ligne de raccord entre 3 images = partage du gateau en 3 + for (INT x = 0 ; x< mosaic.sz().x; x++) + { + for (INT y = 0 ; y< mosaic.sz().y; y++) + { + Pt2di pos(x,y); + if (mChamferDist[i].Inside(Pt2di(pos+tr)) & NbIm.GetR(pos)==3) + { + + // enveloppe interne + /* + if (mChamferDist[i].GetR(Pt2di(pos+tr))<0 & mChamferDist[i].GetR(Pt2di(pos+tr)) >=-mDist ) + { + int lutPos = min( mDist-SumDist.GetR(pos)/2, mChamferDist[i].GetR(Pt2di(pos+tr))) ; + std::cout << " min of " << mDist-SumDist.GetR(pos)/2<< " and chamber dist " << -mChamferDist[i].GetR(Pt2di(pos+tr)) <<" give " << lutPos << "\n"; + //double val = (1-PondInterne.GetR(pos)) * lut_w.At(lutPos) + PondInterne.GetR(pos); + //PondInterne.SetR(pos,val); + } + + */ + if (mChamferDist[i].GetR(Pt2di(pos+tr))>0 & mChamferDist[i].GetR(Pt2di(pos+tr)) <=mDist & mChamferDist[i].GetR(Pt2di(pos+tr))> SumDist.GetR(pos)-mChamferDist[i].GetR(Pt2di(pos+tr)) ) + { + double val = mIm2Ds[i].GetR(Pt2di(pos+tr)) * (1-PondInterne.GetR(pos)) + mosaic.GetR(pos); + //double ratio = mChamferDist[i].GetR(Pt2di(pos+tr))/SumDist.GetR(pos); + //double val = mIm2Ds[i].GetR(Pt2di(pos+tr)) * (1-PondInterne.GetR(pos)) * ratio + mosaic.GetR(pos); + // std::cout << " Chamfer dist vaut " << mChamferDist[i].GetR(Pt2di(pos+tr)) << " et somme chanfer dist vaux " << SumDist.GetR(pos) << ", ratio de " << ratio << "\n"; + mosaic.SetR(pos,val); + } + if (mChamferDist[i].GetR(Pt2di(pos+tr))>0 & mChamferDist[i].GetR(Pt2di(pos+tr)) <=mDist & mChamferDist[i].GetR(Pt2di(pos+tr))== SumDist.GetR(pos)-mChamferDist[i].GetR(Pt2di(pos+tr)) ) + { + double val = 0.5* mIm2Ds[i].GetR(Pt2di(pos+tr)) * (1-PondInterne.GetR(pos)) + mosaic.GetR(pos); + mosaic.SetR(pos,val); + } + + } + } + } + + + } + + aNameTmp="TestPondInterne2.tif"; + ELISE_COPY( + PondInterne.all_pts(), + PondInterne.in(), + Tiff_Im(aNameTmp.c_str(), + PondInterne.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + + aNameTmp="TestMosaic.tif"; + + ELISE_COPY( + mosaic.all_pts(), + mosaic.in(), + Tiff_Im(aNameTmp.c_str(), + mosaic.sz(), + GenIm::real4, + Tiff_Im::No_Compr, + Tiff_Im::BlackIsZero).out() + ); + + } + +} + + + +// the VarioCam thermic camera record images at 16 bits, we want to convert them to 8 bits. A range of temperature is provided in order to stretch the radiometric value on this range + +class cVarioCamTo8Bits +{ + public: + std::string mDir; + cVarioCamTo8Bits(int argc,char ** argv); + private: + std::string mFullDir; + std::string mPat; + std::string mPrefix; + bool mOverwrite; + Pt2di mRangeT; + bool mCelcius; +}; + + +cVarioCamTo8Bits::cVarioCamTo8Bits(int argc,char ** argv) : + mFullDir ("img.*.tif"), + mPrefix ("8bits_"), + mOverwrite (false), + mCelcius(1) +{ + ElInitArgMain + ( + argc,argv, + LArgMain() << EAMC(mFullDir,"image pattern", eSAM_IsPatFile) + << EAMC(mRangeT,"temperature range"), + LArgMain() << EAM(mOverwrite,"F",true, "Overwrite previous output images, def false") + << EAM(mCelcius,"Celcius",true, "Is the temperature range in celcius, default true, if false, Kelvin") + << EAM(mPrefix,"Prefix",true, "Prefix for output images") + ); + + + if (!MMVisualMode) + { + + SplitDirAndFile(mDir,mPat,mFullDir); + cInterfChantierNameManipulateur * aICNM = cInterfChantierNameManipulateur::BasicAlloc(mDir); + const std::vector aSetIm = *(aICNM->Get(mPat)); + + Pt2di aRangeVario; + // convert the range to + if (mCelcius) { + aRangeVario.x=100*(273.15+mRangeT.x) ; + aRangeVario.y=100*(273.15+mRangeT.y) ; + } else {aRangeVario=mRangeT;}; + std::cout << "Range of radiometric value of variocam images : " << aRangeVario << "\n"; + + for (auto & im : aSetIm) + { + std::string NameOut(mDir+mPrefix+im); + + if (ELISE_fp::exist_file(NameOut) & !mOverwrite) + { + std::cout <<"Image " << NameOut <<" already exist, use F=1 to overwrite.\n"; + } else { + + int minRad(aRangeVario.x), rangeRad(aRangeVario.y-aRangeVario.x); + + // load input variocam images + Tiff_Im mTifIn=Tiff_Im::StdConvGen(mDir+im,1,true); + // create empty RAM image for imput image + Im2D_REAL4 imIn(mTifIn.sz().x,mTifIn.sz().y); + // create empty RAM image for output image + Im2D_U_INT1 imOut(mTifIn.sz().x,mTifIn.sz().y); + // fill it with tiff image value + ELISE_COPY( + mTifIn.all_pts(), + mTifIn.in(), + imIn.out() + ); + // change radiometry and note min and max value + int aMin(255), aMax(0),aSum(0),aNb(0); + for (int v(0); v=minRad && aVal =minRad+rangeRad) val=255.0; + } + + if (val>aMax) aMax=val; + if (val!=0){ + if (val aSetIm = *(aICNM->Get(mPat)); + + // the bad way to do this, because I cannot find how to change the optical center of a camera with micmac classes + std::string aNameTmp("Tmp-OriTrans.txt"); + std::string aCom = MMDir() + + std::string("bin/mm3d") + + + " OriExport Ori-" + + mOriIn + std::string("/Orientation-") + + mPat + std::string(".xml") + + std::string(" ") + aNameTmp; + std::cout << aCom << "\n"; + + system_call(aCom.c_str()); + + aCom= MMDir() + + std::string("bin/mm3d OriConvert '#F=N_X_Y_Z_W_P_K' ") + + aNameTmp + std::string(" ") + + std::string("OriTrans") + + std::string(" OffsetXYZ=") + ToString(mTr); + + std::cout << aCom << "\n"; + + system_call(aCom.c_str()); + + // je fais une bascule la dessus + + // je vérifie également avec un oriexport que l'offset fonctionnne + + aCom="cp Ori-"+mOriIn +"/AutoCal* Ori-"+mOriOut+"/"; + system_call(aCom.c_str()); + + /* aCom= MMDir() + + std::string("bin/mm3d GCP '#F=N_X_Y_Z_S_S_S' ") + + aNameTmp + std::string(" ") + + std::string(mOriOut) + + std::string(" OffsetXYZ=") + ToString(mTr); + + std::cout << aCom << "\n"; + + system_call(aCom.c_str()); + +*/ + + + } + +} // Applique une homographie à l'ensemble des images thermiques pour les mettres dans la géométrie des images visibles prises simultanément class cTIR2VIS_Appli; @@ -312,7 +1291,7 @@ void cTIR2VIS_Appli::changeImRadiom(std::vector aLIm) /* comparaise des orthos thermiques pour déterminer un éventuel facteur de calibration spectrale entre 2 frame successif, expliquer pouquoi tant de variabilité spectrale est présente (mosaique moche) */ -// à priori ce n'est pas ça du tout, déjà mauvaise registration TIR --> vis du coup les ortho TIR ne se superposent pas , du coup correction radiometrique ne peut pas fonctionner. +// à priori ce n'est pas ça du tout, déjà mauvaise registration TIR --> vis du coup les ortho TIR ne se superposent pas , du coup correction metrique ne peut pas fonctionner. int CmpOrthosTir_main(int argc,char ** argv) { std::string aDir, aPat="Ort_.*.tif", aPrefix="ratio"; @@ -351,7 +1330,7 @@ int CmpOrthosTir_main(int argc,char ** argv) for (auto aCurrentImGeo: mLIm) { i++; - for (int j=i ; jNON +-angle --> PAS TESTE +moins probable mais je teste quand même: +-position sur capteur --> NON +-temps écoulé depuis début du vol PAS TESTE + */ +int statRadianceVarioCam_main(int argc,char ** argv) +{ + std::string a2DMesFileName, a3DMesFileName, aOutputFile, aOri; + + ElInitArgMain + ( + argc,argv, + //mandatory arguments + LArgMain() << EAMC(a2DMesFileName, "Input mes2D file", eSAM_IsExistFile) + << EAMC(a3DMesFileName, "Input mes3D file", eSAM_IsExistFile) + << EAMC(aOri, "Orientation", eSAM_IsExistDirOri), + LArgMain() + << EAM(aOutputFile,"Out", true, "Output .txt file with radiance observation for statistic", eSAM_IsOutputFile) + + ); + cInterfChantierNameManipulateur * aICNM = cInterfChantierNameManipulateur::BasicAlloc("./"); + + // open 2D measures + cSetOfMesureAppuisFlottants aSetOfMesureAppuisFlottants=StdGetFromPCP(a2DMesFileName,SetOfMesureAppuisFlottants); + // open 3D measures + cDicoAppuisFlottant DAF= StdGetFromPCP(a3DMesFileName,DicoAppuisFlottant); + std::list & aLGCP = DAF.OneAppuisDAF(); + + // create a map of GCP and position + std::map aGCPmap; + + for (auto & GCP : aLGCP) + { + aGCPmap[GCP.NamePt()]=Pt3dr(GCP.Pt().x,GCP.Pt().y,GCP.Pt().z); + } + + + std::cout << "Image GCP U V rayon Radiance GroundDist \n" ; + + + for( std::list< cMesureAppuiFlottant1Im >::const_iterator iTmes1Im=aSetOfMesureAppuisFlottants.MesureAppuiFlottant1Im().begin(); + iTmes1Im!=aSetOfMesureAppuisFlottants.MesureAppuiFlottant1Im().end(); iTmes1Im++ ) + { + cMesureAppuiFlottant1Im anImTIR=*iTmes1Im; + + // open the image + if (ELISE_fp::exist_file(anImTIR.NameIm())) + { + + Tiff_Im mTifIn=Tiff_Im::StdConvGen(anImTIR.NameIm(),1,true); + // create empty RAM image + Im2D_REAL4 im(mTifIn.sz().x,mTifIn.sz().y); + // fill it with tiff image value + ELISE_COPY( + mTifIn.all_pts(), + mTifIn.in(), + im.out() + ); + // open the CamStenope + std::string aNameOri=aICNM->Assoc1To1("NKS-Assoc-Im2Orient@-"+aOri+"@",anImTIR.NameIm(), true); + CamStenope * aCam(CamOrientGenFromFile(aNameOri,aICNM)); + //std::cout << "Optical center Cam" << aCam->PseudoOpticalCenter() << "\n"; + // loop on the + for (auto & appuiTIR : anImTIR.OneMesureAF1I()) + { + // je garde uniquement les GCP dont le nom commence par L + if (appuiTIR.NamePt().substr(0,1)=="L") + { + Pt3dr pt = aGCPmap[appuiTIR.NamePt()]; + // std::cout << " Image " << anImTIR.NameIm() << " GCP " << appuiTIR.NamePt() << " ground position " << pt << " Image position " << appuiTIR.PtIm() << " \n"; + + double aRadiance(0); + int aNb(0); + + Pt2di UV(appuiTIR.PtIm()); + Pt2di sz(1,1); + ELISE_COPY( + rectangle(UV-sz,UV+Pt2di(2,2)*sz),// not sure how it work + Virgule(im.in(),1), + Virgule(sigma(aRadiance),sigma(aNb)) + ); + aRadiance/=aNb; + std::cout << " Radiance on windows of " << aNb << " px " << aRadiance << " \n"; + aRadiance=aRadiance-27315; + // now determine incid and distance from camera to GCP + double aDist(0); + Pt3dr vDist=aCam->PseudoOpticalCenter()-pt; + aDist=euclid(vDist); + + double aDistUV=euclid(appuiTIR.PtIm()-aCam->PP()); + + + std::cout << anImTIR.NameIm() << " " << appuiTIR.NamePt() << " " << appuiTIR.PtIm().x << " " << appuiTIR.PtIm().y << " " << aDistUV << " " << aRadiance << " " << aDist << " \n"; + + } + } + + } + // fin iter sur les mesures appuis flottant + } +} int MasqTIR_main(int argc,char ** argv) @@ -848,17 +1930,12 @@ int MasqTIR_main(int argc,char ** argv) Tiff_Im::BlackIsZero ).out() ); - - } return EXIT_SUCCESS; } - - - -int main_test(int argc,char ** argv) +int main_test2(int argc,char ** argv) { //cORT_Appli anAppli(argc,argv); //CmpOrthosTir_main(argc,argv); @@ -866,13 +1943,62 @@ int main_test(int argc,char ** argv) //RegTIRVIS_main(argc,argv); //test_main(argc,argv); //MasqTIR_main(argc,argv); - //cERO_ModelOnePaire(argc,argv); - TransfoMesureAppuisVario2TP_main(argc,argv); + //cCoreg2Ortho(argc,argv); + cFeatheringAndMosaicOrtho(argc,argv); + //cOriTran_Appli(argc,argv); + //TransfoMesureAppuisVario2TP_main(argc,argv); + //statRadianceVarioCam_main(argc,argv); return EXIT_SUCCESS; } +int main_testold(int argc,char ** argv) +{ + // manipulate the + ofstream fout("/home/lisein/data/DIDRO/lp17/GPS_RGP/GNSS_pos/test.obs"); + ifstream fin("/home/lisein/data/DIDRO/lp17/GPS_RGP/GNSS_pos/tmp.txt"); + string line; + + std::string add(" 40.000 40.000\n"); + add="\t40.000\t40.000\n"; + if (fin.is_open()) + { + unsigned int i(0); + while ( getline(fin,line) ) + { + i++; + + if (i==17){ + i=0; + fout << line << add; + + } else { + + if (i%2==1 && i>2) { + fout << line << add; + } else { + fout << line <<"\n";} + } + } + + } else { std:cout << "cannot open file in\n";} + fin.close(); + fout.close(); + + return EXIT_SUCCESS; +} + +int VarioCamTo8Bits_main(int argc,char ** argv) +{ + + cVarioCamTo8Bits(argc,argv); + + return EXIT_SUCCESS; +} + + + /*Footer-MicMac-eLiSe-25/06/2007 diff --git a/src/TpMMPD/Sources.cmake b/src/TpMMPD/Sources.cmake index 566dc20ec8..76058fe172 100644 --- a/src/TpMMPD/Sources.cmake +++ b/src/TpMMPD/Sources.cmake @@ -110,6 +110,7 @@ set(Src_TD_PPMD ${TDPPMD_DIR}/DIDRO/divers.cpp ${TDPPMD_DIR}/DIDRO/ascii2tif.cpp ${TDPPMD_DIR}/DIDRO/ctapioca_idr.cpp + ${TDPPMD_DIR}/DIDRO/cAppuis2homol.cpp ${TDPPMD_DIR}/../RegTIRVIS/Arbre.cpp ${TDPPMD_DIR}/../RegTIRVIS/Image.cpp ${TDPPMD_DIR}/../RegTIRVIS/DescriptorExtractor.cpp From 5ba122affeac1dd13c688986015d2363874f83af Mon Sep 17 00:00:00 2001 From: JoLeBelge Date: Fri, 2 Mar 2018 10:07:57 +0100 Subject: [PATCH 2/9] feathering ortho --- src/TpMMPD/DIDRO/ascii2tif.cpp | 138 ++++- src/TpMMPD/DIDRO/cfeatheringandmosaicking.cpp | 499 +++++++++++++++ src/TpMMPD/DIDRO/cfeatheringandmosaicking.h | 50 ++ src/TpMMPD/DIDRO/divers.cpp | 575 ++---------------- src/TpMMPD/Sources.cmake | 1 + 5 files changed, 721 insertions(+), 542 deletions(-) create mode 100644 src/TpMMPD/DIDRO/cfeatheringandmosaicking.cpp create mode 100644 src/TpMMPD/DIDRO/cfeatheringandmosaicking.h diff --git a/src/TpMMPD/DIDRO/ascii2tif.cpp b/src/TpMMPD/DIDRO/ascii2tif.cpp index 409c826f80..5cf59001c2 100644 --- a/src/TpMMPD/DIDRO/ascii2tif.cpp +++ b/src/TpMMPD/DIDRO/ascii2tif.cpp @@ -38,20 +38,80 @@ English : Header-MicMac-eLiSe-25/06/2007*/ #include "StdAfx.h" +template +double VarLapl(Im2D * aIm,int aSzW); +//burriness is computed as variance of Laplacian with a mean filter prior to reduce noise +//https://www.pyimagesearch.com/2015/09/07/blur-detection-with-opencv/ + + +template +double VarLapl(Im2D * aIm,int aSzW) +{ + Im2D_REAL4 aRes(aIm->sz().x,aIm->sz().y); + + Fonc_Num aF = aIm->in_proj(); + + int aNbVois = ElSquare(1+2*aSzW); + + aF = rect_som(aF,aSzW) /aNbVois; + ELISE_COPY(aIm->all_pts() + ,aF + ,aRes.out()); + + double mean; + int nb; + + ELISE_COPY(aRes.all_pts() + ,Virgule(Laplacien(aRes.in_proj()) ,1) + ,Virgule(aRes.out()|sigma(mean) ,sigma(nb))); + + std::cout << "mean of laplacian : " << mean << "\n"; + mean=mean/nb; + + std::cout << "nb of value: " << nb << "\n"; + // variance of laplacian + Fonc_Num aFVar = ElSquare((aRes.in()-mean)); + + double meanVarLap; + ELISE_COPY(aRes.all_pts() + ,aFVar + ,sigma(meanVarLap)); + // mean of variance + meanVarLap/=nb; + return meanVarLap; +} + + + int ascii2tif(int argc,char ** argv) { std::string aPat("Ree*.tif"); std::list mLFile; - bool aC2K(true); + bool aVario(true); + Pt2di aSz(1024,768); + int aHLine(0); + bool aC2K(false), aDebug(0), aOptris(0), aInt2Deg(0); + ElInitArgMain ( argc,argv, - LArgMain() << EAMC(aPat,"ascii image pattern, images generated by irbis software (variocam hd thermal camera)",eSAM_IsPatFile), + LArgMain() << EAMC(aPat,"ascii image pattern, images generated either by irbis software (variocam hd thermal camera) or with SDKdirect optris library",eSAM_IsPatFile), LArgMain() - << EAM(aC2K,"C2K", true, "conversion from cercius value to kelvin value *100, in order to got integer value instead of float.") + << EAM(aVario,"Vario", true, "Are ascii generated by irbis soft for camera Variocam? Def true.") + << EAM(aC2K,"C2K", true, "Variocam: conversion from cercius value to kelvin value *100, in order to got integer value instead of float.") + << EAM(aSz,"Sz", true, "Image size, default [1024,768] which is the size of variocam thermal imagery") + << EAM(aOptris,"Optris", true, "Are ascii generated from optris camera? Def false. if true, set proper image size.") + << EAM(aInt2Deg,"Int2Deg", true, "Optris: conversion of unsigned integer value (raw data) to degree?.") + << EAM(aHLine,"DropHL", true, "Number of Line to drop in the Header, default 10 (variocam ascii irbis soft)") + << EAM(aDebug,"Debug", true, "Display much more info in order to help debugging process") ); + if (aVario & !EAMIsInit(&aC2K)) aC2K=true; + if (aVario & !EAMIsInit(&aHLine)) aHLine=10; + + if (aOptris & !EAMIsInit(&aSz)) aSz=Pt2di(640,480); + if (aVario & !EAMIsInit(&aInt2Deg)) aInt2Deg=true; cInterfChantierNameManipulateur * aICNM = cInterfChantierNameManipulateur::BasicAlloc("./"); // create the list of images starting from the regular expression (Pattern) @@ -59,13 +119,12 @@ int ascii2tif(int argc,char ** argv) std::cout << mLFile.size() << " ascii images detected\n"; - int SzX(1024), SzY(768); - for (auto & aName : mLFile){ - std::cout << "Ascii to tif, image " <> output;} + if (aHLine>0) for (int i(0); i<(aHLine-1);i++){ ascii >> output;} for (int l(0); l < SzY; l++) { @@ -81,14 +140,19 @@ int ascii2tif(int argc,char ** argv) { Pt2di aPos(col,l); ascii >> output; - output[2]='.'; + if (aVario) output[2]='.'; FromString(aVal,output); - // std::cout << "Valeur lue: " << aVal <<" , position " << aPos << "\n"; + + if (aDebug) + { + if (col==0 | col==SzX-1) std::cout << "Valeur lue: " << aVal <<" , position " << aPos << "\n"; + } + aImRAM.SetR(aPos,aVal); } } - if (aC2K) + if (aVario & aC2K) { std::cout << "conversion celcius vers Kelvin et multiplication par 100.\n"; @@ -101,22 +165,60 @@ int ascii2tif(int argc,char ** argv) } + if (aOptris & aInt2Deg) + { + + std::cout << "conversion unsigned int Optris PI to degree value.\n"; + ELISE_COPY + ( + aImRAM.all_pts(), + (aImRAM.in()-1000)/10, + aImRAM.out() + ); + + } + + std::cout << "bluriness of this images: " << VarLapl(&aImRAM,3) << "\n"; + std::string aNameTiff=aName.substr(0, aName.size()-3) + "tif"; std::cout << "Ecriture de " <StdGetListOfFile(mFullDir); + + cFileOriMnt MTD = StdGetFromPCP("MTDOrtho.xml",FileOriMnt); + Box2dr boxMosaic(Pt2dr(MTD.OriginePlani().x,MTD.OriginePlani().y+MTD.ResolutionPlani().y*MTD.NombrePixels().y),Pt2dr(MTD.OriginePlani().x+MTD.ResolutionPlani().x*MTD.NombrePixels().x,MTD.OriginePlani().y)); + sz= Pt2di(MTD.NombrePixels().x,MTD.NombrePixels().y); + aCorner=Pt2dr(boxMosaic._p0.x,boxMosaic._p1.y); // xmin, ymax; + + mosaic=Im2D_REAL4(sz.x,sz.y); + NbIm=Im2D_REAL4(sz.x,sz.y); + label=Im2D_U_INT1::FromFileStd("Label.tif"); + + mSumDistExt=Im2D_INT4(sz.x,sz.y); + mSumDistInter= Im2D_INT4(sz.x,sz.y); + PondInterne=Im2D_REAL4(sz.x,sz.y); + mSumWeighting=Im2D_REAL4(sz.x,sz.y); + std::string aName; + + // look up table of weight + lut_w=Im1D_REAL4(mDist+1); + ELISE_COPY + ( + lut_w.all_pts(), + pow(0.5,pow((FX/((double)mDist-FX+1)),2*mLambda)), + lut_w.out() + ); + // replace by 1 distance over the threshold of mDist + ELISE_COPY(select(lut_w.all_pts(),lut_w.in()>1),1,lut_w.out()); + + if (mDebug) { + std::cout << "Look up table of weight depending on the distance from seamline\n"; + for (unsigned int i(0); i i2d) const +{ + int border(200); + Im2D_U_INT1 tmp(i2d.sz().x+2*border,i2d.sz().y+2*border,1); + ELISE_COPY(select(tmp.all_pts(),trans(i2d.in(1),-Pt2di(border,border))==0),0,tmp.out()); + Chamfer::d32.im_dist(tmp); + ELISE_COPY(i2d.all_pts(),trans(tmp.in(255),Pt2di(border,border)),i2d.oclip()); +} + + +void cFeatheringAndMosaicOrtho::ChamferDist4AllOrt() +{ + unsigned int i(0); // i is the label of the image - and the key + for (auto &im : mLFile) + { + if (mDebug) std::cout << "Image num " << i << " is " << im <<" : loading and computing feathering buffer.\n"; + // open orthos + mIms[i]= new cImGeo(mDir+im); + mIm2Ds[i]= mIms[i]->toRAM(); + mChamferDist[i]=Im2D_INT2(mIms[i]->Im().sz().x,mIms[i]->Im().sz().y,1); + Pt2di sz(mIms[i]->Im().sz()); + Pt2di tr= -mIms[i]->computeTrans(aCorner); + + // la fonction chamfer fonctionne sur une image binaire et va calculer les distance à partir des pixels qui ont une valeur de 0. + // la distance maximale est de 255 + + //detect seam line for this image + //1) translation of label to be in ortho geometry and set value to 0 for points that are not inside area of mosaicking for this image + Im2D_U_INT1 tmp(sz.x,sz.y,1); + ELISE_COPY(select(tmp.all_pts(), trans(label.in(),tr)!=(int)i), + //trans(label.in(),tr), + 0, + tmp.out() + ); + + // remove very small patch for wich we do not want to perform feathering because it is ugly otherwise + Im2D_U_INT1 Ibin(tmp.sz().x,tmp.sz().y,0); + ELISE_COPY(Ibin.all_pts(),tmp.in(),Ibin.out()); + // completely stupid but i have to ensure the border of the bin image is not ==1 otherwise I got the error out of bitmap in dilate spec Image + // so the code may suffer weakness if a small patch is located near the border--ask marc + ELISE_COPY(Ibin.border(2),0,Ibin.out()); + U_INT1 ** d = Ibin.data(); + Neighbourhood V8=Neighbourhood::v8(); + + for (INT x=0; x < Ibin.sz().x; x++) + { + for (INT y=0; y < Ibin.sz().y; y++) + { + if (d[y][x] == 1) + { + Liste_Pts_INT2 cc(2); + ELISE_COPY + ( + // flux: composantes connexes du point. + conc + ( + Pt2di(x,y), + Ibin.neigh_test_and_set(V8, 1, 0, 10) ), // on change la valeur des points sélectionné comme ça à la prochaine itération on ne sélectionne plus cette zone de composante connexe + 2, // valeur bidonne, c'est juste le flux que je sauve dans cc + cc + ); + // remove the patch + if (cc.card()