ENH: implementation of streamed label image writing done

app/otbLSGRM.cxx

@@ -17,6 +17,10 @@
 #include "lsgrmFullLambdaScheduleSegmenter.h"
 #include "lsgrmController.h"
 
+// Graph to label image (streaming version)
+#include "otbStreamingGraphToImageFilter.h"
+#include "otbStreamingImageVirtualWriter.h"
+
 // system tools
 #include <itksys/SystemTools.hxx>
 
@@ -108,6 +112,7 @@ private:
     AddChoice("tiling.none", "No tiling layout");
 
     AddParameter(ParameterType_Int, "memory", "Restrict memory use (mb)");
+
     MandatoryOff("memory");
   }
 
@@ -129,28 +134,28 @@ private:
     std::string tmpdir;
     if (HasValue("tmpdir"))
       {
-      tmpdir= GetParameterAsString("tmpdir");
+        tmpdir= GetParameterAsString("tmpdir");
       }
     else
       {
-      tmpdir = itksys::SystemTools::GetFilenamePath(outfname);
+        tmpdir = itksys::SystemTools::GetFilenamePath(outfname);
       }
 
     if (!tmpdir.empty())
       {
-      // A temporary directory is specified: we check that it ends with a POSIX separator
-      if (tmpdir[tmpdir.size()-1] != '/')
-        {
-        // If not, we add the separator
-        tmpdir.append("/");
-        }
-
-      // Check that the directory exists
-      if (!itksys::SystemTools::FileExists(tmpdir.c_str(),false))
-        {
-        otbAppLogFATAL("The directory " << tmpdir << " does not exist.");
-        }
-      otbAppLogINFO("Using temporary directory " << tmpdir);
+        // A temporary directory is specified: we check that it ends with a POSIX separator
+        if (tmpdir[tmpdir.size()-1] != '/')
+          {
+            // If not, we add the separator
+            tmpdir.append("/");
+          }
+
+        // Check that the directory exists
+        if (!itksys::SystemTools::FileExists(tmpdir.c_str(),false))
+          {
+            otbAppLogFATAL("The directory " << tmpdir << " does not exist.");
+          }
+        otbAppLogINFO("Using temporary directory " << tmpdir);
       }
 
 
@@ -163,17 +168,20 @@ private:
    * This function sets the generic parameters of a controller and runs the segmentation
    */
   template<class TSegmenter>
-  UInt32ImageType::Pointer SetGenericParametersAndRunSegmentation(const typename TSegmenter::ParamType params){
+  void
+  SetGenericParametersAndRunSegmentation(const typename TSegmenter::ParamType params){
 
     // Instantiate the controller
     typedef typename lsgrm::Controller<TSegmenter> ControlerType;
     typename ControlerType::Pointer controller = ControlerType::New();
+    using GraphType = typename ControlerType::GraphType;
 
     // Set specific parameters
     controller->SetSpecificParameters(params);
 
     // Set input image
-    controller->SetInputImage(GetParameterFloatVectorImage("in"));
+    ImageType::Pointer inputImage = GetParameterFloatVectorImage("in");
+    controller->SetInputImage(inputImage);
 
     // Set threshold
     float thres = GetParameterFloat("threshold");
@@ -189,85 +197,92 @@ private:
     int inputTilingMode = GetParameterInt("tiling");
     if (inputTilingMode == TILING_AUTO)
       {
-      // Automatic mode
-      controller->SetTilingModeAuto();
+        // Automatic mode
+        controller->SetTilingModeAuto();
       }
     else if (inputTilingMode == TILING_USER)
       {
-      // User mode
-      controller->SetTilingModeUser();
-      controller->SetTileWidth(GetParameterInt("tiling.user.sizex"));
-      controller->SetTileHeight(GetParameterInt("tiling.user.sizey"));
-      controller->SetNumberOfFirstIterations(GetParameterInt("tiling.user.nfirstiter"));
+        // User mode
+        controller->SetTilingModeUser();
+        controller->SetTileWidth(GetParameterInt("tiling.user.sizex"));
+        controller->SetTileHeight(GetParameterInt("tiling.user.sizey"));
+        controller->SetNumberOfFirstIterations(GetParameterInt("tiling.user.nfirstiter"));
       }
     else if (inputTilingMode == TILING_NONE)
       {
-      // None mode
-      controller->SetTilingModeNone();
+        // None mode
+        controller->SetTilingModeNone();
       }
     else
       {
-      otbAppLogFATAL("Unknown tiling mode!");
+        otbAppLogFATAL("Unknown tiling mode!");
       }
 
     // Input RAM value?
     if (HasValue("memory"))
       {
-      otbAppLogINFO("Setting maximum memory to " << GetParameterInt("memory") << " MBytes");
-      controller->SetInternalMemoryAvailable(GetParameterInt("memory"));
+        otbAppLogINFO("Setting maximum memory to " << GetParameterInt("memory") << " MBytes");
+        controller->SetInternalMemoryAvailable(GetParameterInt("memory"));
       }
 
     // Run the segmentation
     controller->RunSegmentation();
 
+    // Prepare the label image source
+    typedef lsgrm::StreamingGraphToImageFilter<GraphType, UInt32ImageType> LabelImageSourceType;
+    typename LabelImageSourceType::Pointer labelImageSource = LabelImageSourceType::New();
+    labelImageSource->SetGraph(controller->GetOutputGraph());
+    labelImageSource->SetOutputSize(inputImage->GetLargestPossibleRegion().GetSize());
+    labelImageSource->SetOutputOrigin(inputImage->GetOrigin());
+    labelImageSource->SetOutputSpacing(inputImage->GetSpacing());
+    labelImageSource->SetOutputProjectionRef(inputImage->GetProjectionRef());
+    labelImageSource->GenerateOutputInformation();
+
+    m_LabelImageSource = static_cast<itk::ImageSource<UInt32ImageType>*>(labelImageSource);
+
+    SetParameterOutputImage<UInt32ImageType>("out", m_LabelImageSource->GetOutput());
+
+    // TODO: find an intelligent value of RAM
+    if (dynamic_cast<OutputImageParameter*>(GetParameterByKey("out")))
+      {
+      OutputImageParameter* paramDown = dynamic_cast<OutputImageParameter*>(GetParameterByKey("out"));
+      paramDown->SetRAMValue(256);
+      }
+
     // Get temporary files list
     m_TemporaryFilesList = controller->GetTemporaryFilesList();
 
-    // Return the label image
-    return controller->GetLabeledClusteredOutput();
   }
 
   void DoExecute()
   {
-    /*
-        TODO: the output directory in case the global graph cannot fit in memory (lsgrmController.txx)
-     */
 
-    // Input image
     ImageType::Pointer inputImage = GetParameterFloatVectorImage("in");
 
     // Switch criterion
     int inputCriterion = GetParameterInt("criterion");
     if (inputCriterion == CRITERION_BAATZ)
       {
-      grm::BaatzParam params;
-      params.m_SpectralWeight = GetParameterFloat("criterion.bs.cw");
-      params.m_ShapeWeight = GetParameterFloat("criterion.bs.sw");
-      labelImage = SetGenericParametersAndRunSegmentation<BaatzSegmenterType>(params);
+        grm::BaatzParam params;
+        params.m_SpectralWeight = GetParameterFloat("criterion.bs.cw");
+        params.m_ShapeWeight = GetParameterFloat("criterion.bs.sw");
+        SetGenericParametersAndRunSegmentation<BaatzSegmenterType>(params);
       }
     else if (inputCriterion == CRITERION_SPRING)
       {
-      grm::SpringParam params;
-      labelImage = SetGenericParametersAndRunSegmentation<SpringSegmenterType>(params);
+        grm::SpringParam params;
+        SetGenericParametersAndRunSegmentation<SpringSegmenterType>(params);
       }
     else if (inputCriterion == CRITERION_FLS)
       {
-      grm::FLSParam params;
-      labelImage = SetGenericParametersAndRunSegmentation<FLSSegmenterType>(params);
+        grm::FLSParam params;
+        SetGenericParametersAndRunSegmentation<FLSSegmenterType>(params);
       }
     else
       {
-      otbAppLogFATAL("Unknow criterion!")
+        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());
-    labelImage->SetSpacing(inputImage->GetSpacing());
-    SetParameterOutputImage<UInt32ImageType>("out", labelImage);
-
   }
 
   void AfterExecuteAndWriteOutputs()
@@ -282,16 +297,16 @@ private:
     // Delete temporary files
     for (unsigned int i = 0 ; i < m_TemporaryFilesList.size() ; i++)
       {
-      if( remove(m_TemporaryFilesList.at(i).c_str() ) != 0  )
-        {
-        otbAppLogWARNING( "Error deleting file " << m_TemporaryFilesList.at(i) );
-        }
+        if( remove(m_TemporaryFilesList.at(i).c_str() ) != 0  )
+          {
+            otbAppLogWARNING( "Error deleting file " << m_TemporaryFilesList.at(i) );
+          }
       }
   }
 
 private:
   std::vector<std::string> m_TemporaryFilesList;
-  UInt32ImageType::Pointer labelImage;
+  itk::ImageSource<UInt32ImageType>::Pointer m_LabelImageSource;
 
 }; // app class
 } // end of namespace wrapper

include/lsgrmController.h

@@ -6,6 +6,7 @@
 #include "itkObject.h"
 #include "itkMacro.h"
 
+
 namespace lsgrm
 {
 template<class TSegmenter>
@@ -30,6 +31,7 @@ public:
   using ImageType = typename SegmenterType::ImageType;
   using LabelImageType = typename SegmenterType::LabelImageType;
   using SegmentationParameterType = typename SegmenterType::ParamType;
+  using GraphType = typename SegmenterType::GraphType;
 
   /* Default constructor and destructor. */
   Controller();
@@ -48,7 +50,6 @@ public:
   void SetTilingModeUser(){m_TilingMode = LSGRM_TILING_USER; Modified();};
   void SetTilingModeAuto(){m_TilingMode = LSGRM_TILING_AUTO; Modified();};
 
-  typename LabelImageType::Pointer GetLabeledClusteredOutput();
   std::vector<std::string> GetTemporaryFilesList();
 
   itkGetMacro(Margin, unsigned int);
@@ -77,6 +78,8 @@ public:
   itkGetMacro(TileHeight, unsigned int);
   itkSetMacro(TileHeight, unsigned int);
 
+  itkGetMacro(OutputGraph, GraphType);
+
 private:
 
   /** Enum for tiling mode */
@@ -116,7 +119,8 @@ private:
   LSGRMTilingMode                    m_TilingMode;         // tiling mode (none/user/auto)
   unsigned int                       m_Margin;             // stability margin related to m_NumberOfFirstIterations
   std::vector<ProcessingTile>        m_Tiles;              // list of tiles
-  typename LabelImageType::Pointer   m_LabelImage;         // output label image
+  GraphType                          m_OutputGraph;        // Output graph
+
 };
 } // end of namespace lsgrm
 

include/lsgrmController.txx

@@ -18,6 +18,7 @@ Controller<TSegmenter>::Controller()
   m_NbTilesY = 0;
   m_Threshold = 75;
   m_Memory = 0;
+
 }
 
 template<class TSegmenter>
@@ -150,7 +151,7 @@ void Controller<TSegmenter>::RunSegmentation()
     if(accumulatedMemory <= m_Memory)
       {
       // Merge all the graphs
-      m_LabelImage = MergeAllGraphsAndAchieveSegmentation<TSegmenter>(
+        m_OutputGraph = MergeAllGraphsAndAchieveSegmentation<TSegmenter>(
           m_SpecificParameters,
           m_Threshold,
           m_Tiles,
@@ -195,14 +196,15 @@ void Controller<TSegmenter>::RunSegmentation()
     segmenter.SetInput(m_InputImage);
     segmenter.Update();
 
-    // Get label image
-    m_LabelImage = segmenter.GetLabeledClusteredOutput();
+    m_OutputGraph = segmenter.m_Graph;
     }
   else
     {
     itkExceptionMacro(<<"Unknow tiling mode!");
     }
 
+  // TODO: [MPI] broadcast the graph to other nodes
+
 }
 
 /*
@@ -433,16 +435,16 @@ void Controller<TSegmenter>::SetSpecificParameters(const SegmentationParameterTy
   m_SpecificParameters = params;
 }
 
-template<class TSegmenter>
-typename Controller<TSegmenter>::LabelImageType::Pointer
-Controller<TSegmenter>::GetLabeledClusteredOutput()
-{
-#ifdef OTB_USE_MPI
-  // Get the label image from the master process (the one which achieves segmentation)
-  BroadcastImage<typename TSegmenter::LabelImageType>(m_LabelImage);
-#endif
-  return m_LabelImage;
-}
+//template<class TSegmenter>
+//typename Controller<TSegmenter>::LabelImageType::Pointer
+//Controller<TSegmenter>::GetLabeledClusteredOutput()
+//{
+//#ifdef OTB_USE_MPI
+//  // Get the label image from the master process (the one which achieves segmentation)
+//  BroadcastImage<typename TSegmenter::LabelImageType>(m_LabelImage);
+//#endif
+//  return m_LabelImage;
+//}
 
 template <class TSegmenter>
 std::vector<std::string> Controller<TSegmenter>::GetTemporaryFilesList()

include/lsgrmGraphOperations.h

@@ -5,9 +5,6 @@
 #include "grmGraphOperations.h"
 #include "otbVectorImage.h"
 #include "otbMultiChannelExtractROI.h"
-#include "itkGrayscaleFillholeImageFilter.h"
-#include "otbStreamingGraphToImageFilter.h"
-#include "otbStreamingImageVirtualWriter.h"
 
 namespace lsgrm
 {
@@ -34,7 +31,7 @@ typename TSegmenter::ImageType::Pointer ReadImageRegion(
     typename TSegmenter::ImageType::RegionType region);
 
 template<class TSegmenter>
-typename TSegmenter::LabelImageType::Pointer
+typename TSegmenter::GraphType
 MergeAllGraphsAndAchieveSegmentation(
     const typename TSegmenter::ParamType& params,
     const float& threshold,

include/lsgrmGraphOperations.txx

@@ -25,7 +25,7 @@ typename TSegmenter::ImageType::Pointer ReadImageRegion(
 }
 
 template<class TSegmenter>
-typename TSegmenter::LabelImageType::Pointer
+typename TSegmenter::GraphType
 MergeAllGraphsAndAchieveSegmentation(
     const typename TSegmenter::ParamType& params,
     const float& threshold,
@@ -91,86 +91,30 @@ MergeAllGraphsAndAchieveSegmentation(
   segmenter.SetThreshold(threshold);
   segmenter.SetDoFastSegmentation(false); // was true
   segmenter.SetNumberOfIterations(numberOfIterations);
-
   grm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
 
-  //  // Write output graph to the output graph directory
-  //  WriteGraph<TSegmenter>(segmenter.m_Graph, tmpDir, 0, 0);
-  
-  // TODO: use std::sort to sort the graph (m_Id must be increasing) before the polygon filling
-  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;
-  index.Fill(0);
-  typename TSegmenter::LabelImageType::SizeType size;
-  size[0] = imageWidth;
-  size[1] = imageHeight;
-  typename TSegmenter::LabelImageType::RegionType region(index, size);
-  const typename TSegmenter::LabelImageType::InternalPixelType noDataLabel = 0;
-  typename TSegmenter::LabelImageType::Pointer labelImage = TSegmenter::LabelImageType::New();
-  labelImage->SetRegions(region);
-  labelImage->Allocate();
-  labelImage->FillBuffer(noDataLabel);
-
-  using LabelImageIterator = itk::ImageRegionIterator<typename TSegmenter::LabelImageType>;
-  LabelImageIterator it(labelImage, labelImage->GetLargestPossibleRegion());
-
-  // Start at 1 (value 0 can be used for invalid pixels)
-  unsigned int label = 1;
-  for(auto& node : segmenter.m_Graph.m_Nodes)
-    {
-    lp::CellLists borderPixels;
-    lp::ContourOperations::GenerateBorderCells(borderPixels, node->m_Contour, node->m_Id, imageWidth);
-
-    for (auto& pix: borderPixels)
-      {
-      index[0] = pix % imageWidth;
-      index[1] = pix / imageWidth;
-      labelImage->SetPixel(index, label);
-      }
-    ++label;
-    }
+  // Sort the nodes top-->down and left-->right
+  std::sort(segmenter.m_Graph.m_Nodes.begin(), segmenter.m_Graph.m_Nodes.end(),
+      [imageWidth](const auto & a, const auto & b) -> bool
+  {
+    lp::CellLists borderPixelsA, borderPixelsB;
+    lp::ContourOperations::GenerateBorderCells(borderPixelsA, a->m_Contour, a->m_Id, imageWidth);
+    lp::ContourOperations::GenerateBorderCells(borderPixelsB, b->m_Contour, b->m_Id, imageWidth);
 
-  // Re-order nodes labels (left->right to top->bottom)
-  vnl_vector<typename TSegmenter::LabelImageType::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++;
-      }
-    }
+    unsigned int offA = 0;
+    unsigned int offB = 0;
+    for (auto& pix: borderPixelsA)
+      if (pix>offA)
+        offA=pix;
 
-  // Apply LUT
-  for(it.GoToBegin();!it.IsAtEnd(); ++it)
-    it.Set(lut[it.Get()]);
+    for (auto& pix: borderPixelsB)
+      if (pix>offB)
+        offB=pix;
 
-  // Fill holes
-  typedef itk::GrayscaleFillholeImageFilter<typename TSegmenter::LabelImageType,
-      typename TSegmenter::LabelImageType> FillholeFilterType;
-  typename FillholeFilterType::Pointer fillFilter = FillholeFilterType::New();
-  fillFilter->SetInput(labelImage);
-  fillFilter->Update();
+      return offA > offB;
+  });
 
-  return fillFilter->GetOutput();
+  return segmenter.m_Graph;
 }
 
 template<class TSegmenter>

include/otbStreamingGraphToImageFilter.h

@@ -54,6 +54,8 @@ public:
   typedef typename TLabelImage::RegionType              RegionType;
   typedef typename TLabelImage::IndexType               IndexType;
   typedef typename TLabelImage::SizeType                SizeType;
+  typedef typename TLabelImage::SpacingType             SpacingType;
+  typedef typename TLabelImage::PointType               PointType;
   typedef typename TGraph::NodePointerType              NodePointerType;
   typedef itk::GrayscaleFillholeImageFilter<TLabelImage,TLabelImage> FillholeFilterType;
   typedef bg::model::point<float, 2, bg::cs::cartesian> point;
@@ -70,9 +72,12 @@ public:
   virtual void GenerateData();
 
   itkSetMacro(OutputSize, SizeType);
+  itkSetMacro(OutputOrigin, PointType);
+  itkSetMacro(OutputSpacing, SpacingType);
+  itkSetMacro(OutputProjectionRef, std::string);
 
 protected:
-  StreamingGraphToImageFilter();
+  StreamingGraphToImageFilter(){};
   virtual ~StreamingGraphToImageFilter(){};
 
 private:
@@ -80,9 +85,12 @@ private:
   StreamingGraphToImageFilter(const Self &); //purposely not implemented
   void operator =(const Self&); //purposely not implemented
 
-  TGraph m_Graph;
+  TGraph       m_Graph;
   bgi::rtree< value, bgi::quadratic<16> > rtree;
-  SizeType m_OutputSize;
+  SizeType     m_OutputSize;
+  SpacingType  m_OutputSpacing;
+  PointType    m_OutputOrigin;
+  std::string  m_OutputProjectionRef;
 };
 
 } /* namespace lsgrm */

include/otbStreamingGraphToImageFilter.txx

@@ -13,13 +13,6 @@
 namespace lsgrm
 {
 
-template <class TGraph, class TLabelImage>
-StreamingGraphToImageFilter<TGraph, TLabelImage>
-::StreamingGraphToImageFilter()
- {
-
- }
-
 template <class TGraph, class TLabelImage>
 void
 StreamingGraphToImageFilter<TGraph, TLabelImage>
@@ -28,23 +21,27 @@ StreamingGraphToImageFilter<TGraph, TLabelImage>
   m_Graph = graph;
 
   // Build a R-Tree
-  std::cout << "Building R-Tree" << std::endl;
+  itkDebugMacro( << "Building R-Tree" );
+
   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));
+      box b(
+          point(
+              node->m_Bbox.m_UX,
+              node->m_Bbox.m_UY),
+          point(
+              node->m_Bbox.m_UX + 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
+  itkDebugMacro( << "Building R-Tree finished" );
 
  }
 
@@ -53,14 +50,19 @@ void
 StreamingGraphToImageFilter<TGraph, TLabelImage>
 ::GenerateOutputInformation()
  {
-  TLabelImage * outputPtr = this->GetOutput();
-  // TODO:
-  //  outputPtr->SetOrigin ( m_OutputOrigin );
-  //  outputPtr->SetSpacing ( m_OutputSpacing );
+  itkDebugMacro( << "Entering GenerateOutputInformation()" );
+
+  // Output Largest Possible Region
   IndexType index;
   index.Fill(0);
   RegionType outputRegion(index, m_OutputSize);
+
+  // Set output informations
+  TLabelImage * outputPtr = this->GetOutput();
+  outputPtr->SetOrigin ( m_OutputOrigin );
+  outputPtr->SetSpacing ( m_OutputSpacing );
   outputPtr->SetLargestPossibleRegion( outputRegion );
+  outputPtr->SetProjectionRef(m_OutputProjectionRef);
 
  }
 
@@ -70,6 +72,8 @@ void
 StreamingGraphToImageFilter<TGraph, TLabelImage>
 ::GenerateData()
  {
+  itkDebugMacro( << "Entering GenerateData()" );
+
   // Allocate the output buffer
   TLabelImage * outputPtr = this->GetOutput();
   RegionType outReqRegion = outputPtr->GetRequestedRegion();
@@ -77,12 +81,17 @@ StreamingGraphToImageFilter<TGraph, TLabelImage>
   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)));
+  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;
+  itkDebugMacro( << "R-Tree query on output requested region " << outReqRegion );
   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;
+  itkDebugMacro( << "R-Tree query done. Number of nodes: " << result_s.size() );
 
   // Retrieve the bounding box of the intersecting nodes (a kind of "Input requested region")
   box realBBox(query_box);
@@ -97,23 +106,22 @@ StreamingGraphToImageFilter<TGraph, TLabelImage>
   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;
+  itkDebugMacro( << "Input Requested region: " << inputRequestedRegion );
 
   // Generate the label image
-  std::cout << "Allocate " << std::endl;
+  itkDebugMacro( << "Allocate buffered region ");
   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;
+  itkDebugMacro( << "Allocate buffered region ok" );
 
   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 ;
+  // Burn boundaries
+  itkDebugMacro( << "Burning boundaries " );
   for(auto& res : result_s)
     {
       NodePointerType node = m_Graph.m_Nodes[res.second];
@@ -125,43 +133,20 @@ StreamingGraphToImageFilter<TGraph, TLabelImage>