Ajout de geolocalisation dans le résultat

MakefileNoOTB

@@ -34,7 +34,8 @@ TTH_OUT = $(patsubst %, $(OUT_DIR)/%, $(TTH_PRG))
 
 ## FLAGS ###############################
 # fast and no control
-DFLAGS = -O2 -DNDEBUG -DNO_OTB -DINTEL_TBB_THREAD -DBOOST_DISABLE_ASSERTS
+#DFLAGS = -O2 -DNDEBUG -DNO_OTB -DINTEL_TBB_THREAD -DBOOST_DISABLE_ASSERTS
+DFLAGS = -O2 -DNDEBUG -DNO_OTB -DBOOST_DISABLE_ASSERTS
 # fast but control
 #DFLAGS = -O2 -DNDEBUG -DNO_OTB -DINTEL_TBB_THREAD
 # debug multi-threaded

include/IImage.hpp

@@ -19,6 +19,7 @@ namespace triskele {
   private:
     Size		size;
     vector<PixelT>	pixels;
+
   public:
     inline void		setSize (const Size &size);
     inline const Size	&getSize () const;
@@ -39,6 +40,8 @@ namespace triskele {
   class IImage {
   public:
     void			setFileName (string fileName);
+    inline void			getGeo (string& projectionRef, vector<double> &geoTransform) const;
+    inline void			setGeo (const string& projectionRef, vector<double> geoTransform, const Point &topLeft);
     inline const string		&getFileName () const;
     inline const Size		&getSize () const;
     inline const DimChanel	&getBandCount () const;
@@ -66,6 +69,8 @@ namespace triskele {
     
   private:
     static size_t	gdalCount;
+    string		projectionRef;
+    vector<double>	geoTransform;
     string		fileName;
     Size		size;
     DimChanel		bandCount;

include/IImage.tpp

@@ -86,6 +86,37 @@ IImage::setFileName (string fileName) {
     close ();
 }
 
+inline void
+IImage::getGeo (string& projectionRef, vector<double> &geoTransform) const {
+  projectionRef = this->projectionRef;
+  geoTransform = this->geoTransform;
+}
+
+inline void
+IImage::setGeo (const string& projectionRef, vector<double> geoTransform, const Point &topLeft) {
+  DEF_LOG ("IImage::setGeo", "projectionRef:" << projectionRef);
+  // move TopLeft cf WordFile
+  double
+    x (geoTransform[0]),
+    xScale (geoTransform[1]),
+    xSkew (geoTransform[2]),
+    y (geoTransform[3]),
+    ySkew (geoTransform[4]),
+    yScale (geoTransform[5]);
+  double
+    xTrans = x + topLeft.x*xScale + topLeft.y*xSkew,
+    yTrans = y + topLeft.x*ySkew + topLeft.y*yScale;
+  geoTransform[0] = xTrans;
+  geoTransform[3] = yTrans;
+  CPLErr err = gdalOutputDataset->SetGeoTransform (&geoTransform[0]);
+  if (err != CE_None)
+    cerr << "IImage::setGeo: can't set geoTransform for " << fileName << endl;
+  err = gdalOutputDataset->SetProjection (projectionRef.c_str ());
+  if (err != CE_None)
+    cerr << "IImage::setGeo: can't set projection for " << fileName << endl;
+}
+
+
 inline const string &
 IImage::getFileName () const {
   return fileName;

include/TreeStats.hpp

@@ -11,6 +11,7 @@
 #include <boost/accumulators/statistics/max.hpp>
 #include <boost/accumulators/statistics/mean.hpp>
 
+#include "triskeleDebug.hpp"
 #include "ArrayTree/triskeleArrayTreeBase.hpp"
 
 namespace otb {

include/triskeleDebug.hpp

@@ -6,6 +6,8 @@
 #include <string>
 #include <sstream>
 #include <boost/date_time/posix_time/posix_time.hpp>
+#include <boost/chrono.hpp>
+
 
 #ifdef ENABLE_SMART_LOG
 
@@ -101,7 +103,11 @@ namespace triskele {
     }
   };
 
+  inline string ns2string (double delta);
+
   // ========================================
+#include "triskeleDebug.tpp"
+
 }//namespace triskele
 
 #endif //_TRISKELE_DEBUG_HPP

src/IImage.cpp

@@ -11,7 +11,8 @@ size_t
 IImage::gdalCount = 0;
 
 IImage::IImage (const string &imageFileName)
-  : fileName (imageFileName),
+  : geoTransform (6, 0),
+    fileName (imageFileName),
     gdalInputDataset (nullptr),
     gdalOutputDataset (nullptr),
     read (false)
@@ -36,11 +37,18 @@ IImage::readImage () {
     cerr << "GDALError: can't define dataset" << endl;
     return;
   }
+
   size = Size (gdalInputDataset->GetRasterXSize (),
 	       gdalInputDataset->GetRasterYSize ());
   bandCount = gdalInputDataset->GetRasterCount ();
   LOG ("size: " << size << " x " << bandCount);
   read = true;
+  projectionRef = gdalInputDataset->GetProjectionRef ();
+  CPLErr err = gdalInputDataset->GetGeoTransform (&geoTransform[0]);
+  if (err != CE_None)
+    cerr << "IImage::readImage: can't read geoTransform from " << fileName << endl;
+  else
+    LOG ("geoTransform: " << geoTransform [0] << " " << geoTransform [1] << " " <<  geoTransform [2] << " " <<  geoTransform [3] << " " <<  geoTransform [4] << " " <<  geoTransform [5]);
 
   for (DimChanel band = 0; band < bandCount; ++band) {
     GDALRasterBand &poBand = *gdalInputDataset->GetRasterBand (band+1);
@@ -81,6 +89,7 @@ IImage::createImage (const Size &size, const GDALDataType &dataType, const DimCh
   GDALDriver *driverTiff = GetGDALDriverManager ()->GetDriverByName ("GTiff");
   remove (fileName.c_str ());
   gdalOutputDataset = driverTiff->Create (fileName.c_str (), size.width, size.height, nbBands, dataType, NULL);
+
   LOG("gdalCount: " << gdalCount);
 }
 
@@ -88,7 +97,6 @@ void
 IImage::close () {
   DEF_LOG ("IImage::close", "fileName: " << fileName);
   if (gdalOutputDataset) {
-    // XXX pour écriture gdalOutputDataset->SetProjection ("WGS84");
     GDALClose (gdalOutputDataset);
     gdalOutputDataset = nullptr;
     if (!--gdalCount)

src/TreeStats.cpp

@@ -1,7 +1,5 @@
 #include <iomanip>
 
-#include <boost/chrono.hpp>
-
 #include "TreeStats.hpp"
 
 using namespace otb::triskele;
@@ -47,26 +45,6 @@ TreeStats::reset () {
 }
 
 // ========================================
-using namespace boost::chrono;
-inline string
-ns2string (double delta) {
-  ostringstream oss;
-  duration<double> ns (delta);
-  oss.fill ('0');
-  // typedef duration<int, ratio<86400> > days;
-  // auto d = duration_cast<days>(ns);
-  // ns -= d;
-  auto h = duration_cast<hours> (ns);
-  ns -= h;
-  auto m = duration_cast<minutes> (ns);
-  ns -= m;
-  oss << setw (2) << h.count () << ":"
-      << setw (2) << m.count () << ":"
-      << setw (9) << fixed << setprecision (6) << ns.count ();
-  return oss.str ();
-}
-
-// ----------------------------------------
 TreeStats::CPrintDim::CPrintDim (const TreeStats &treeStats)
   : treeStats (treeStats) {
 }

src/apGenerator.cpp

@@ -45,6 +45,12 @@ writeBand (Option &option, PixelT *pixels, DimChanel band) {
   fileNameStream << outputBaseName << "-" << std::setfill ('0') << std::setw (3) << (band) << outputExtension;
   IImage outputImage (fileNameStream.str ());
   outputImage.createImage (option.size, option.inputImage.getDataType (), 1);
+  // XXX en double
+  string projectionRef;
+  vector<double> geoTransform;
+  option.inputImage.getGeo (projectionRef, geoTransform);
+  outputImage.setGeo (projectionRef, geoTransform, option.topLeft);
+
   outputImage.writeBand (pixels, 0);
 }
 
@@ -73,8 +79,14 @@ void apGenerator (Option &option) {
   vector <vector <PixelT> > allBands (maxThresholds, vector<PixelT> (leafCount, 0));
 
   DimChanel outputBandsCard = option.selectedBand.getSet ().size ()*(1+treeTypesCard*(option.areaThresholds.size ()+option.levelThresholds.size ()+option.sdThresholds.size ()+option.moiThresholds.size ()));
-  if (!option.oneBand)
+  if (!option.oneBand) {
     option.outputImage.createImage (option.size, option.inputImage.getDataType (), outputBandsCard);
+    // XXX en double
+    string projectionRef;
+    vector<double> geoTransform;
+    option.inputImage.getGeo (projectionRef, geoTransform);
+    option.outputImage.setGeo (projectionRef, geoTransform, option.topLeft);
+  }
 
   Raster<PixelT> raster;
   if (option.border) {

test/TestDeal.cpp

@@ -0,0 +1,91 @@
+// faire test 10000 + lamba + 1...coreCount
+#include <boost/chrono.hpp>
+#include <boost/thread.hpp>
+
+#include "triskeleDealThreads.hpp"
+#include "TestThread.hpp"
+
+using namespace std;
+using namespace boost::chrono;
+using namespace triskele;
+using namespace otb;
+using namespace otb::triskele;
+
+inline void
+fDirect (const DimImg &nbItem, const vector<DimImg> &global) {
+  DimImg sum = 0;
+  for (DimImg x = 0; x < nbItem; ++x)
+    sum += global[x];
+}
+
+template<typename FunctId>
+inline void
+fLambda (DimImg &nbItem, const FunctId &functId/* functId (id) */) {
+  for (DimImg x = 0; x < nbItem; ++x)
+    functId (x);
+}
+
+template<typename FunctId>
+inline void
+fThread (double &inDuration, const DimImg &nbItem, const FunctId &functId/* functId (id) */) {
+#if INTEL_TBB_THREAD
+  using namespace tbb;
+#pragma warning(disable: 588)
+  parallel_for (size_t (0), size_t (1), [&nbItem, &functId] (size_t idCopyValInThread) {
+      auto start = high_resolution_clock::now ();
+      for (DimImg x = 0; x < nbItem; ++x)
+	functId (x);
+      auto end = high_resolution_clock::now ();
+      inDuration = duration_cast<duration<double> > (end-start).count ();
+    });
+#else /* BOOST thread */
+  std::vector<boost::thread> tasks;
+  for (unsigned int idCopyValInThread = 0; idCopyValInThread < 1; ++idCopyValInThread) {
+    tasks.push_back (boost::thread ([/*no ref!!!*/idCopyValInThread, &nbItem, &functId, &inDuration] () {
+	  auto start = high_resolution_clock::now ();
+	  for (DimImg x = 0; x < nbItem; ++x)
+	    functId (x);
+	  auto end = high_resolution_clock::now ();
+	  inDuration += duration_cast<duration<double> > (end-start).count ();
+	}));
+  }
+  for (unsigned int i = 0; i < 1; ++i)
+    tasks[i].join ();
+#endif
+}
+
+void
+TestThread::testDeal () {
+
+  DimImg nbItem = global.size ();
+  double inDuration = 0;
+  // XXX direct
+  auto startDirect = high_resolution_clock::now ();
+  fDirect (nbItem, global);
+
+  // XXX lambda
+  auto startLamba = high_resolution_clock::now ();
+  DimImg sum = 0;
+  fLambda (nbItem, [this, &sum] (const DimImg &x) {
+      sum += global [x];
+    });
+
+  // XXX thread
+  auto startThread = high_resolution_clock::now ();
+  sum = 0;
+  fThread (inDuration, nbItem, [this, &sum] (const DimImg &x) {
+      sum += global [x];
+    });
+
+  // XXX n direct
+
+  // XXX n lambda
+
+  // XXX n thread
+
+  auto end = high_resolution_clock::now ();
+  addTime (directDealStats, duration_cast<duration<double> > (startLamba-startDirect).count ());
+  addTime (lambdaDealStats, duration_cast<duration<double> > (startThread-startLamba).count ());
+  addTime (threadDealStats, duration_cast<duration<double> > (end-startThread).count ());
+  addTime (inThreadDealStats, inDuration);  
+}

test/TestThread.cpp

@@ -12,6 +12,12 @@ using namespace otb;
 using namespace otb::triskele;
 
 static string timeTypeLabels [TimeTypeCard] = {
+  "directDeal",
+  "lambdaDeal",
+  "threadDeal",
+  "inThreadDeal",
+
+
   "initStats",
   "seqReadStats",
   "parReadStats",
@@ -21,6 +27,9 @@ static string timeTypeLabels [TimeTypeCard] = {
   "parRWStats"
 };
 
+const unsigned int
+TestThread::maxCoreCount = boost::thread::hardware_concurrency ();
+
 // ========================================
 inline string
 TestThread::ns2string (double delta) {
@@ -81,8 +90,9 @@ nodealThread (const DimImg &maxId, unsigned int coreCount, const FunctThreadMinM
 }
 
 // ========================================
-TestThread::TestThread ()
+TestThread::TestThread (const unsigned int &coreCount)
-  : global (nbItem*coreCount, 0)
+  : coreCount (coreCount),
+    global (nbItem*coreCount, 0)
 {
 }
 
@@ -120,23 +130,22 @@ TestThread::fillVector (vector<T> &vect) {
 
 void
 TestThread::multiTest () {
-  for (unsigned int nbCore = 1; nbCore <= coreCount; ++nbCore) {
   for (int i = 0; i < 100; ++i) {
-      vector<long> sumSeq (nbCore, 0);
+    vector<long> sumSeq (coreCount, 0);
-      vector<long> sumPar (nbCore, 0);
+    vector<long> sumPar (coreCount, 0);
 
     auto start = high_resolution_clock::now ();
     fillVector (global);
 
     // lecture seq => faire somme
     auto startSeqRead = high_resolution_clock::now ();
-      nodealThreadThreadRange (nbItem, nbCore, [this, &sumSeq] (const unsigned int &threadId, const DimImg &item) {
+    nodealThreadThreadRange (nbItem, coreCount, [this, &sumSeq] (const unsigned int &threadId, const DimImg &item) {
 	sumSeq[threadId] += global[item];
       });
 
     // lecture // => faire somme
     auto startParRead = high_resolution_clock::now ();
-      dealThreadThreadRange (nbItem, nbCore, [this, &sumPar] (const unsigned int &threadId, const DimImg &item) {
+    dealThreadThreadRange (nbItem, coreCount, [this, &sumPar] (const unsigned int &threadId, const DimImg &item) {
 	sumPar[threadId] += global[item];
       });
     // XXX vérifier égalité de sumSeq sumPar
@@ -154,15 +163,19 @@ TestThread::multiTest () {
   }
   print (cout, timeStats);
 }
-}
 
 // ========================================
 int
 main (int argc, char** argv) {
   cout << "start test" << endl;
   srand (time (NULL));
+
   TestThread tt;
-  tt.multiTest ();
+  //tt.multiTest ();
+  for (int i = 0; i < 100; ++i)
+    tt.testDeal ();
+  tt.print (cout, tt.timeStats);
+
   return 0;
 }