ADD: add a otbStreamingGraphToImageFilter

app/otbLSGRM.cxx

@@ -65,6 +65,7 @@ private:
 
   void DoInit()
   {
+
     SetName("GenericRegionMerging");
     SetDescription("This application allows to use the Large Scale Generic Region Merging library "
         "(LSGRM) and provides currently 3 homogeneity criteria: Euclidean Distance, "
@@ -235,9 +236,6 @@ private:
     // Input image
     ImageType::Pointer inputImage = GetParameterFloatVectorImage("in");
 
-    // Output image
-    UInt32ImageType::Pointer labelImage = UInt32ImageType::New();
-
     // Switch criterion
     int inputCriterion = GetParameterInt("criterion");
     if (inputCriterion == CRITERION_BAATZ)
@@ -262,6 +260,8 @@ private:
       otbAppLogFATAL("Unknow criterion!")
       }
 
+    std::cout << "Label image size: " << labelImage->GetLargestPossibleRegion() << std::endl;
+
     // Set output image projection, origin and spacing for labelImage
     labelImage->SetProjectionRef(inputImage->GetProjectionRef());
     labelImage->SetOrigin(inputImage->GetOrigin());
@@ -291,6 +291,7 @@ private:
 
 private:
   std::vector<std::string> m_TemporaryFilesList;
+  UInt32ImageType::Pointer labelImage;
 
 }; // app class
 } // end of namespace wrapper

include/lsgrmGraphOperations.h

@@ -6,6 +6,8 @@
 #include "otbVectorImage.h"
 #include "otbMultiChannelExtractROI.h"
 #include "itkGrayscaleFillholeImageFilter.h"
+#include "otbStreamingGraphToImageFilter.h"
+#include "otbStreamingImageVirtualWriter.h"
 
 namespace lsgrm
 {

include/lsgrmGraphOperations.txx

@@ -96,6 +96,21 @@ MergeAllGraphsAndAchieveSegmentation(
 
   //  // Write output graph to the output graph directory
   //  WriteGraph<TSegmenter>(segmenter.m_Graph, tmpDir, 0, 0);
+  typename TSegmenter::LabelImageType::SizeType outSize;
+  outSize[0] = imageWidth;
+  outSize[1] = imageHeight;
+  typedef StreamingGraphToImageFilter<typename TSegmenter::GraphType,
+      typename TSegmenter::LabelImageType> GraphToImageFilterType;
+  typename GraphToImageFilterType::Pointer graphToImageFilter = GraphToImageFilterType::New();
+  graphToImageFilter->SetGraph(segmenter.m_Graph);
+  graphToImageFilter->SetOutputSize(outSize);
+  typedef otb::ImageFileWriter<typename TSegmenter::LabelImageType> VWriterType;
+
+  typename VWriterType::Pointer wrt = VWriterType::New();
+  wrt->SetInput(graphToImageFilter->GetOutput());
+  wrt->SetNumberOfDivisionsStrippedStreaming(5);
+  wrt->SetFileName("/homeL/cresson/data/tmp/spot6_surrech_segm2.tif");
+  wrt->Update();
 
   // Generate the label image
   typename TSegmenter::LabelImageType::IndexType index;
@@ -153,7 +168,8 @@ MergeAllGraphsAndAchieveSegmentation(
   fillFilter->SetInput(labelImage);
   fillFilter->Update();
 
-  return fillFilter->GetOutput();}
+  return fillFilter->GetOutput();
+}
 
 template<class TSegmenter>
 long long unsigned int RunPartialSegmentation(const typename TSegmenter::ParamType& params,

include/otbStreamingGraphToImageFilter.h

+/*
+ * otbStreamingGraphToImageFilter.h
+ *
+ *  Created on: 6 nov. 2017
+ *      Author: cresson
+ */
+
+#ifndef MODULES_REMOTE_LSGRM_INCLUDE_OTBSTREAMINGGRAPHTOIMAGEFILTER_H_
+#define MODULES_REMOTE_LSGRM_INCLUDE_OTBSTREAMINGGRAPHTOIMAGEFILTER_H_
+
+#include "itkImageSource.h"
+#include "itkExceptionObject.h"
+#include "itkImageRegion.h"
+#include "itkGrayscaleFillholeImageFilter.h"
+
+#include "grmGraphOperations.h"
+
+// Boost R-Tree
+#include <boost/geometry.hpp>
+#include <boost/geometry/geometries/point.hpp>
+#include <boost/geometry/geometries/box.hpp>
+
+#include <boost/geometry/index/rtree.hpp>
+
+// to store queries results
+#include <vector>
+
+// just for output
+#include <iostream>
+#include <boost/foreach.hpp>
+
+namespace bg = boost::geometry;
+namespace bgi = boost::geometry::index;
+
+namespace lsgrm
+{
+
+template <class TGraph, class TLabelImage>
+class ITK_EXPORT StreamingGraphToImageFilter : public itk::ImageSource<TLabelImage>
+{
+public:
+  /** Standard class typedefs. */
+  typedef StreamingGraphToImageFilter                   Self;
+  typedef itk::ImageSource<TLabelImage>                 Superclass;
+  typedef itk::SmartPointer<Self>                       Pointer;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(StreamingGraphToImageFilter, ImageSource);
+
+  /** Typedefs for processing */
+  typedef typename TLabelImage::RegionType              RegionType;
+  typedef typename TLabelImage::IndexType               IndexType;
+  typedef typename TLabelImage::SizeType                SizeType;
+  typedef typename TGraph::NodePointerType              NodePointerType;
+  typedef itk::GrayscaleFillholeImageFilter<TLabelImage,TLabelImage> FillholeFilterType;
+  typedef bg::model::point<float, 2, bg::cs::cartesian> point;
+  typedef bg::model::box<point> box;
+  typedef std::pair<box, unsigned> value;
+
+  /** Set the input graph and build the R-Tree */
+  void SetGraph(const TGraph graph);
+
+  /** Prepare image allocation at the first call of the pipeline processing */
+  virtual void GenerateOutputInformation(void);
+
+  /** Does the real work. */
+  virtual void GenerateData();
+
+  itkSetMacro(OutputSize, SizeType);
+
+protected:
+  StreamingGraphToImageFilter();
+  virtual ~StreamingGraphToImageFilter(){};
+
+private:
+
+  StreamingGraphToImageFilter(const Self &); //purposely not implemented
+  void operator =(const Self&); //purposely not implemented
+
+  TGraph m_Graph;
+  bgi::rtree< value, bgi::quadratic<16> > rtree;
+  SizeType m_OutputSize;
+};
+
+} /* namespace lsgrm */
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include <otbStreamingGraphToImageFilter.txx>
+#endif
+
+#endif /* MODULES_REMOTE_LSGRM_INCLUDE_OTBSTREAMINGGRAPHTOIMAGEFILTER_H_ */

include/otbStreamingGraphToImageFilter.txx

+/*
+ * otbStreamingGraphToImageFilter.txx
+ *
+ *  Created on: 6 nov. 2017
+ *      Author: cresson
+ */
+
+#ifndef MODULES_REMOTE_LSGRM_INCLUDE_OTBSTREAMINGGRAPHTOIMAGEFILTER_TXX_
+#define MODULES_REMOTE_LSGRM_INCLUDE_OTBSTREAMINGGRAPHTOIMAGEFILTER_TXX_
+
+#include <otbStreamingGraphToImageFilter.h>
+
+namespace lsgrm
+{
+
+template <class TGraph, class TLabelImage>
+StreamingGraphToImageFilter<TGraph, TLabelImage>
+::StreamingGraphToImageFilter()
+ {
+
+ }
+
+template <class TGraph, class TLabelImage>
+void
+StreamingGraphToImageFilter<TGraph, TLabelImage>
+::SetGraph(const TGraph graph)
+ {
+  m_Graph = graph;
+
+  // Build a R-Tree
+  std::cout << "Building R-Tree" << std::endl;
+  rtree.clear();
+  unsigned int count = 0;
+  for(auto& node : m_Graph.m_Nodes)
+    {
+      // create a box
+      box b(point(node->m_Bbox.m_UX, node->m_Bbox.m_UY),
+          point(node->m_Bbox.m_W, node->m_Bbox.m_UY+node->m_Bbox.m_H));
+
+      // insert new value
+      rtree.insert(std::make_pair(b, count));
+      count++;
+    }
+
+  std::cout << "Building R-Tree finished" << std::endl;
+
+  // TODO: create a LUT for reordering labels from top-->down and left-->right
+
+ }
+
+template <class TGraph, class TLabelImage>
+void
+StreamingGraphToImageFilter<TGraph, TLabelImage>
+::GenerateOutputInformation()
+ {
+  TLabelImage * outputPtr = this->GetOutput();
+  // TODO:
+  //  outputPtr->SetOrigin ( m_OutputOrigin );
+  //  outputPtr->SetSpacing ( m_OutputSpacing );
+  IndexType index;
+  index.Fill(0);
+  RegionType outputRegion(index, m_OutputSize);
+  outputPtr->SetLargestPossibleRegion( outputRegion );
+
+ }
+
+
+template <class TGraph, class TLabelImage>
+void
+StreamingGraphToImageFilter<TGraph, TLabelImage>
+::GenerateData()
+ {
+  // Allocate the output buffer
+  TLabelImage * outputPtr = this->GetOutput();
+  RegionType outReqRegion = outputPtr->GetRequestedRegion();
+  outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
+  outputPtr->Allocate();
+
+  // Find nodes intersecting find the output requested region
+  box query_box(point(outReqRegion.GetIndex(0), outReqRegion.GetIndex(1)),
+      point(outReqRegion.GetIndex(0)+outReqRegion.GetSize(0), outReqRegion.GetIndex(1)+outReqRegion.GetSize(1)));
+  std::vector<value> result_s;
+  std::cout << "R-Tree query on output requested region " << outReqRegion << std::endl;
+  rtree.query(bgi::intersects(query_box), std::back_inserter(result_s));
+  std::cout << "R-Tree query done. Number of nodes: " << result_s.size() << std::endl;
+
+  // Retrieve the bounding box of the intersecting nodes (a kind of "Input requested region")
+  box realBBox(query_box);
+  for(auto& res : result_s)
+    {
+      boost::geometry::expand(realBBox, res.first);
+    }
+  IndexType index;
+  index[0] = realBBox.min_corner().get<0>();
+  index[1] = realBBox.min_corner().get<1>();
+  SizeType size;
+  size[0] = realBBox.max_corner().get<0>() - realBBox.min_corner().get<0>();
+  size[1] = realBBox.max_corner().get<1>() - realBBox.min_corner().get<1>();
+  RegionType inputRequestedRegion(index, size);
+  std::cout << "Input Requested region: " << inputRequestedRegion << std::endl;
+
+  // Generate the label image
+  std::cout << "Allocate " << std::endl;
+  const typename TLabelImage::InternalPixelType noDataLabel = 0;
+  typename TLabelImage::Pointer labelImage = TLabelImage::New();
+  labelImage->SetRegions(inputRequestedRegion);
+  labelImage->Allocate();
+  labelImage->FillBuffer(noDataLabel);
+  std::cout << "Allocate ok" << std::endl;
+
+  using LabelImageIterator = itk::ImageRegionIterator<TLabelImage>;
+  LabelImageIterator it(labelImage, inputRequestedRegion);
+
+  // Start at 1 (value 0 can be used for invalid pixels)
+  std::cout << "Iterate " << std::endl;
+  unsigned int label = 1 ;
+  for(auto& res : result_s)
+    {
+      NodePointerType node = m_Graph.m_Nodes[res.second];
+
+      lp::CellLists borderPixels;
+      lp::ContourOperations::GenerateBorderCells(borderPixels, node->m_Contour, node->m_Id, m_OutputSize[0]);
+
+      for (auto& pix: borderPixels)
+        {
+          index[0] = pix % m_OutputSize[0];
+          index[1] = pix / m_OutputSize[0];
+          labelImage->SetPixel(index, label);
+        }
+      ++label;
+    }
+  std::cout << "Iterate ok " << std::endl;
+
+  // Re-order nodes labels (left->right to top->bottom)
+  std::cout << "re-order " << std::endl;
+  vnl_vector<typename TLabelImage::InternalPixelType> lut(label,noDataLabel);
+  label = 1;
+  for(it.GoToBegin();!it.IsAtEnd(); ++it)
+    {
+      unsigned int inputLabel = it.Get();
+      if (lut[ inputLabel ] == noDataLabel && inputLabel != noDataLabel)
+        {
+          lut[ inputLabel ] = label;
+          label++;
+        }
+    }
+  std::cout << "re-order ok" << std::endl;
+
+  // Apply LUT
+  std::cout << "apply LUT" << std::endl;
+  for(it.GoToBegin();!it.IsAtEnd(); ++it)
+    it.Set(lut[it.Get()]);
+  std::cout << "apply LUT ok" << std::endl;
+
+  // Fill holes
+  std::cout << "fill Hole filter" << std::endl;
+  typename FillholeFilterType::Pointer fillFilter = FillholeFilterType::New();
+  fillFilter->SetInput(labelImage);
+  fillFilter->Update();
+  std::cout << "Fill Hole filter OutputLargestPossibleRegion:" << fillFilter->GetOutput()->GetLargestPossibleRegion() << std::endl;
+  std::cout << "fill Hole filter ok" << std::endl;
+
+  // Extract the stable region
+
+  LabelImageIterator outIt(outputPtr, outReqRegion);
+  LabelImageIterator inIt (fillFilter->GetOutput(), outReqRegion);
+  for (inIt.GoToBegin(), outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt, ++inIt)
+    outIt.Set(inIt.Get());
+ }
+
+
+
+} /* namespace lsgrm */
+
+#endif