Merge branch 'sampling_refac_threading' into develop

Modules/Learning/Sampling/include/otbImageSampleExtractorFilter.h

@@ -69,8 +69,8 @@ public:
 
   /** Get the output samples OGR container */
   ogr::DataSource* GetOutputSamples();
-  
-  virtual void Synthetize(void);
+
+  virtual void Synthetize(void){}
 
   /** Reset method called before starting the streaming*/
   virtual void Reset(void);
@@ -95,7 +95,7 @@ protected:
   virtual void GenerateInputRequestedRegion();
 
   /** process only points */
-  virtual void ThreadedGenerateData(const RegionType&, itk::ThreadIdType threadid);
+  virtual void ThreadedGenerateVectorData(const ogr::Layer& layerForThread, itk::ThreadIdType threadid);
 
 private:
   PersistentImageSampleExtractorFilter(const Self &); //purposely not implemented

Modules/Learning/Sampling/include/otbImageSampleExtractorFilter.txx

@@ -55,15 +55,6 @@ PersistentImageSampleExtractorFilter<TInputImage>
   return static_cast<ogr::DataSource *>(this->itk::ProcessObject::GetOutput(1));
 }
 
-template<class TInputImage>
-void
-PersistentImageSampleExtractorFilter<TInputImage>
-::Synthetize(void)
-{
-  // clear temporary outputs
-  this->m_InMemoryOutputs.clear();
-}
-
 template<class TInputImage>
 void
 PersistentImageSampleExtractorFilter<TInputImage>
@@ -165,20 +156,15 @@ PersistentImageSampleExtractorFilter<TInputImage>
 template<class TInputImage>
 void
 PersistentImageSampleExtractorFilter<TInputImage>
-::ThreadedGenerateData(const RegionType&, itk::ThreadIdType threadid)
+::ThreadedGenerateVectorData(const ogr::Layer& layerForThread, itk::ThreadIdType threadid)
 {
   // Retrieve inputs
   TInputImage* inputImage = const_cast<TInputImage*>(this->GetInput());
   unsigned int nbBand = inputImage->GetNumberOfComponentsPerPixel();
 
-  ogr::Layer layer = this->m_InMemoryInputs[threadid]->GetLayerChecked(0);
-  if (! layer)
-    {
-    return;
-    }
-  ogr::Layer outputLayer = this->m_InMemoryOutputs[threadid][0]->GetLayerChecked(0);
+  ogr::Layer outputLayer = this->GetInMemoryOutput(threadid);
 
-  itk::ProgressReporter progress( this, threadid, layer.GetFeatureCount(true) );
+  itk::ProgressReporter progress( this, threadid, layerForThread.GetFeatureCount(true) );
 
   // Loop across the features in the layer (filtered by requested region in BeforeTGD already)
   OGRGeometry *geom;
@@ -186,8 +172,8 @@ PersistentImageSampleExtractorFilter<TInputImage>
   IndexType imgIndex;
   PixelType imgPixel;
   double imgComp;
-  ogr::Layer::const_iterator featIt = layer.begin();
-  for(; featIt!=layer.end(); ++featIt)
+  ogr::Layer::const_iterator featIt = layerForThread.begin();
+  for(; featIt!=layerForThread.end(); ++featIt)
     {
     geom = featIt->ogr().GetGeometryRef();
     switch (geom->getGeometryType())

Modules/Learning/Sampling/include/otbOGRDataToSamplePositionFilter.h

@@ -85,8 +85,7 @@ public:
   /** Runtime information support. */
   itkTypeMacro(PersistentOGRDataToSamplePositionFilter, PersistentSamplingFilterBase);
 
-  /** Synthetize the persistent filter*/
-  virtual void Synthetize(void);
+  virtual void Synthetize(void){}
 
   /** Reset method called before starting the streaming*/
   virtual void Reset(void);
@@ -135,7 +134,7 @@ protected:
   /** Method to split the input OGRDataSource
    *  according to the class partition
    */
-  virtual void DispatchInputVectors(ogr::Layer &inLayer, std::vector<ogr::Layer> &tmpLayers);
+  virtual void DispatchInputVectors(void);
 
 private:
   PersistentOGRDataToSamplePositionFilter(const Self &); //purposely not implemented

Modules/Learning/Sampling/include/otbOGRDataToSamplePositionFilter.txx

@@ -33,15 +33,6 @@ PersistentOGRDataToSamplePositionFilter<TInputImage,TMaskImage,TSampler>
   m_OriginFieldName = std::string("originfid");
 }
 
-template<class TInputImage, class TMaskImage, class TSampler>
-void
-PersistentOGRDataToSamplePositionFilter<TInputImage,TMaskImage,TSampler>
-::Synthetize(void)
-{
-  // clear temporary outputs
-  this->m_InMemoryOutputs.clear();
-}
-
 template<class TInputImage, class TMaskImage, class TSampler>
 void
 PersistentOGRDataToSamplePositionFilter<TInputImage,TMaskImage,TSampler>
@@ -185,7 +176,7 @@ PersistentOGRDataToSamplePositionFilter<TInputImage,TMaskImage,TSampler>
       ogrTmpPoint.setX(imgPoint[0]);
       ogrTmpPoint.setY(imgPoint[1]);
 
-      ogr::Layer outputLayer = this->m_InMemoryOutputs[threadid][i]->GetLayerChecked(0);
+      ogr::Layer outputLayer = this->GetInMemoryOutput(threadid,i);
       ogr::Feature feat(outputLayer.GetLayerDefn());
       feat.SetFrom(feature);
       feat[this->GetOriginFieldName()].SetValue(static_cast<int>(feature.GetFID()));
@@ -199,8 +190,43 @@ PersistentOGRDataToSamplePositionFilter<TInputImage,TMaskImage,TSampler>
 template<class TInputImage, class TMaskImage, class TSampler>
 void
 PersistentOGRDataToSamplePositionFilter<TInputImage,TMaskImage,TSampler>
-::DispatchInputVectors(ogr::Layer &inLayer, std::vector<ogr::Layer> &tmpLayers)
+::DispatchInputVectors()
 {
+  TInputImage* outputImage = this->GetOutput();
+  ogr::DataSource* vectors = const_cast<ogr::DataSource*>(this->GetOGRData());
+  ogr::Layer inLayer = vectors->GetLayer(this->GetLayerIndex());
+
+  const RegionType& requestedRegion = outputImage->GetRequestedRegion();
+  itk::ContinuousIndex<double> startIndex(requestedRegion.GetIndex());
+  itk::ContinuousIndex<double> endIndex(requestedRegion.GetUpperIndex());
+  startIndex[0] += -0.5;
+  startIndex[1] += -0.5;
+  endIndex[0] += 0.5;
+  endIndex[1] += 0.5;
+  itk::Point<double, 2> startPoint;
+  itk::Point<double, 2> endPoint;
+  outputImage->TransformContinuousIndexToPhysicalPoint(startIndex, startPoint);
+  outputImage->TransformContinuousIndexToPhysicalPoint(endIndex, endPoint);
+
+  // create geometric extent
+  OGRPolygon tmpPolygon;
+  OGRLinearRing ring;
+  ring.addPoint(startPoint[0],startPoint[1],0.0);
+  ring.addPoint(startPoint[0],endPoint[1]  ,0.0);
+  ring.addPoint(endPoint[0]  ,endPoint[1]  ,0.0);
+  ring.addPoint(endPoint[0]  ,startPoint[1],0.0);
+  ring.addPoint(startPoint[0],startPoint[1],0.0);
+  tmpPolygon.addRing(&ring);
+
+  inLayer.SetSpatialFilter(&tmpPolygon);
+
+  unsigned int numberOfThreads = this->GetNumberOfThreads();
+  std::vector<ogr::Layer> tmpLayers;
+  for (unsigned int i=0 ; i<numberOfThreads ; i++)
+    {
+    tmpLayers.push_back(this->GetInMemoryInput(i));
+    }
+
   OGRFeatureDefn &layerDefn = inLayer.GetLayerDefn();
   ogr::Layer::const_iterator featIt = inLayer.begin();
   std::string className;
@@ -212,6 +238,8 @@ PersistentOGRDataToSamplePositionFilter<TInputImage,TMaskImage,TSampler>
     className = featIt->ogr().GetFieldAsString(this->GetFieldIndex());
     tmpLayers[m_ClassPartition[className]].CreateFeature( dstFeature );
     }
+
+  inLayer.SetSpatialFilter(ITK_NULLPTR);
 }
 
 template<class TInputImage, class TMaskImage, class TSampler>

Modules/Learning/Sampling/include/otbPersistentSamplingFilterBase.h

@@ -104,15 +104,14 @@ protected:
    *  the requested region for the mask */
   virtual void GenerateInputRequestedRegion();
 
-  /** Prepare temporary input and output OGR data sources */
-  virtual void BeforeThreadedGenerateData(void);
+  /** Generate data should thread over */
+  virtual void GenerateData(void);
 
-  /** Gather data from multiple threads and
-   *  write to output OGRDataSource (if any) */
-  virtual void AfterThreadedGenerateData(void);
+  /** Allocate in-memory layers for input and outputs */
+  virtual void AllocateOutputs(void);
 
   /** Start of main processing loop */
-  virtual void ThreadedGenerateData(const RegionType&, itk::ThreadIdType threadid);
+  virtual void ThreadedGenerateVectorData(const ogr::Layer& layerForThread, itk::ThreadIdType threadid);
 
   /** Process a geometry, recursive method when the geometry is a collection */
   void ExploreGeometry(const ogr::Feature& feature,
@@ -155,28 +154,17 @@ protected:
   /** Get the region bounding a set of features */
   RegionType FeatureBoundingRegion(const TInputImage* image, otb::ogr::Layer::const_iterator& featIt) const;
 
-  /** Prepares in-memory layers for each thread,
-   *  then calls DispatchInputVectors() for the actual dispatch */
-  virtual void PrepareInputVectors();
-
   /** Method to split the input OGRDataSource between several containers
    *  for each thread. Default is to put the same number of features for
    *  each thread.*/
-  virtual void DispatchInputVectors(ogr::Layer &inLayer, std::vector<ogr::Layer> &tmpLayers);
+  virtual void DispatchInputVectors(void);
 
-  /** Prepare output feature containers for each thread and each output
-   *  OGRDataSource*/
-  virtual void PrepareOutputVectors();
+  /** Gather the content of in-memory output layer into the filter outputs */
+  virtual void GatherOutputVectors(void);
 
   /** Utility method to add new fields on an output layer */
   virtual void InitializeOutputDataSource(ogr::DataSource* inputDS, ogr::DataSource* outputDS);
 
-  /** In-memory containers storing input geometries for each thread*/
-  std::vector<OGRDataPointer> m_InMemoryInputs;
-
-  /** In-memory containers storing position during iteration loop*/
-  std::vector<std::vector<OGRDataPointer> > m_InMemoryOutputs;
-
   typedef struct {
     std::string Name;
     OGRFieldType Type;
@@ -196,6 +184,21 @@ protected:
   /** Get a reference over the additional fields */
   const std::vector<SimpleFieldDefn>& GetAdditionalFields();
 
+  /** Callback function to launch VectorThreadedGenerateData in each thread */
+  static ITK_THREAD_RETURN_TYPE VectorThreaderCallback(void *arg);
+
+  /** basically the same struct as itk::ImageSource::ThreadStruct */
+  struct VectorThreadStruct
+    {
+      Pointer Filter;
+    };
+
+  /** Give access to in-memory input layers */
+  ogr::Layer GetInMemoryInput(unsigned int threadId);
+
+  /** Give access to in-memory output layers */
+  ogr::Layer GetInMemoryOutput(unsigned int threadId, unsigned int index=0);
+
 private:
   PersistentSamplingFilterBase(const Self &); //purposely not implemented
   void operator =(const Self&); //purposely not implemented
@@ -218,6 +221,12 @@ private:
   /** Additional field definitions to add in output data sources */
   std::vector<SimpleFieldDefn> m_AdditionalFields;
 
+  /** In-memory containers storing input geometries for each thread*/
+  std::vector<OGRDataPointer> m_InMemoryInputs;
+
+  /** In-memory containers storing position during iteration loop*/
+  std::vector<std::vector<OGRDataPointer> > m_InMemoryOutputs;
+
 };
 } // End namespace otb
 

Modules/Learning/Sampling/include/otbPersistentSamplingFilterBase.txx

@@ -158,34 +158,119 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
     }
 }
 
-
 template <class TInputImage, class TMaskImage>
 void
 PersistentSamplingFilterBase<TInputImage,TMaskImage>
-::BeforeThreadedGenerateData(void)
+::GenerateData(void)
 {
-  this->PrepareInputVectors();
+  this->AllocateOutputs();
+  this->BeforeThreadedGenerateData();
+
+  // Split the data into in-memory layers
+  this->DispatchInputVectors();
+
+  // struct to store filter pointer
+  VectorThreadStruct str;
+  str.Filter = this;
+
+  // Get the output pointer
+  //const InputImageType *outputPtr = this->GetOutput();
+
+  this->GetMultiThreader()->SetNumberOfThreads( this->GetNumberOfThreads() );
+  this->GetMultiThreader()->SetSingleMethod(this->VectorThreaderCallback, &str);
+
+  // multithread the execution
+  this->GetMultiThreader()->SingleMethodExecute();
 
-  this->PrepareOutputVectors();
+  // gather the data from in-memory output layers
+  this->GatherOutputVectors();
+
+  this->AfterThreadedGenerateData();
 }
 
 template <class TInputImage, class TMaskImage>
 void
 PersistentSamplingFilterBase<TInputImage,TMaskImage>
-::AfterThreadedGenerateData(void)
+::AllocateOutputs(void)
 {
-  // clean temporary inputs
-  this->m_InMemoryInputs.clear();
+  Superclass::AllocateOutputs();
+
+  ogr::DataSource* vectors = const_cast<ogr::DataSource*>(this->GetOGRData());
+  ogr::Layer inLayer = vectors->GetLayer(m_LayerIndex);
 
   unsigned int numberOfThreads = this->GetNumberOfThreads();
-  
-  unsigned int actualNumberOfThreads = numberOfThreads;
-  
-  if(numberOfThreads > this->GetOutput()->GetRequestedRegion().GetSize()[1])
+
+  // Prepare temporary input
+  this->m_InMemoryInputs.clear();
+  this->m_InMemoryInputs.reserve(numberOfThreads);
+  std::string tmpLayerName("thread");
+  OGRSpatialReference * oSRS = ITK_NULLPTR;
+  if (inLayer.GetSpatialRef())
     {
-    actualNumberOfThreads = this->GetOutput()->GetRequestedRegion().GetSize()[1];
+    oSRS = inLayer.GetSpatialRef()->Clone();
+    }
+  OGRFeatureDefn &layerDefn = inLayer.GetLayerDefn();
+  //std::vector<ogr::Layer> tmpLayers;
+  for (unsigned int i=0 ; i < numberOfThreads ; i++)
+    {
+    ogr::DataSource::Pointer tmpOgrDS = ogr::DataSource::New();
+    ogr::Layer tmpLayer = tmpOgrDS->CreateLayer(
+      tmpLayerName,
+      oSRS,
+      inLayer.GetGeomType());
+    // add field definitions
+    for (int k=0 ; k < layerDefn.GetFieldCount() ; k++)
+      {
+      OGRFieldDefn originDefn(layerDefn.GetFieldDefn(k));
+      ogr::FieldDefn fieldDefn(originDefn);
+      tmpLayer.CreateField(fieldDefn);
+      }
+    this->m_InMemoryInputs.push_back(tmpOgrDS);
     }
 
+  // Prepare in-memory outputs
+  this->m_InMemoryOutputs.clear();
+  this->m_InMemoryOutputs.reserve(numberOfThreads);
+  tmpLayerName = std::string("threadOut");
+  for (unsigned int i=0 ; i < numberOfThreads ; i++)
+    {
+    std::vector<OGRDataPointer> tmpContainer;
+    // iterate over outputs, only process ogr::DataSource
+    for (unsigned int k=0 ; k < this->GetNumberOfOutputs() ; k++)
+      {
+      ogr::DataSource* realOutput = dynamic_cast<ogr::DataSource *>(
+        this->itk::ProcessObject::GetOutput(k));
+      if (realOutput)
+        {
+        ogr::Layer realLayer = realOutput->GetLayersCount() == 1
+                               ? realOutput->GetLayer(0)
+                               : realOutput->GetLayer(m_OutLayerName);
+        OGRFeatureDefn &outLayerDefn = realLayer.GetLayerDefn();
+        ogr::DataSource::Pointer tmpOutput = ogr::DataSource::New();
+        ogr::Layer tmpLayer = tmpOutput->CreateLayer(
+          tmpLayerName, oSRS,  realLayer.GetGeomType());
+        // add field definitions
+        for (int f=0 ; f < outLayerDefn.GetFieldCount() ; f++)
+          {
+          OGRFieldDefn originDefn(outLayerDefn.GetFieldDefn(f));
+          tmpLayer.CreateField(originDefn);
+          }
+        tmpContainer.push_back(tmpOutput);
+        }
+      }
+    this->m_InMemoryOutputs.push_back(tmpContainer);
+    }
+}
+
+template <class TInputImage, class TMaskImage>
+void
+PersistentSamplingFilterBase<TInputImage,TMaskImage>
+::GatherOutputVectors(void)
+{
+  // clean temporary inputs
+  this->m_InMemoryInputs.clear();
+
+  unsigned int numberOfThreads = this->GetNumberOfThreads();
   
   // gather temporary outputs and write to output
   const otb::ogr::DataSource* vectors = this->GetOGRData();
@@ -208,7 +293,7 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
         itkExceptionMacro(<< "Unable to start transaction for OGR layer " << outLayer.ogr().GetName() << ".");
         }
   
-      for (unsigned int thread=0 ; thread < actualNumberOfThreads ; thread++)
+      for (unsigned int thread=0 ; thread < numberOfThreads ; thread++)
         {
         ogr::Layer inLayer = this->m_InMemoryOutputs[thread][count]->GetLayerChecked(0);
         if (!inLayer)
@@ -248,29 +333,24 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
     }
   chrono.Stop();
   otbMsgDebugMacro(<< "write ogr points took " << chrono.GetTotal() << " sec");
+  this->m_InMemoryOutputs.clear();
 }
 
 template <class TInputImage, class TMaskImage>
 void
 PersistentSamplingFilterBase<TInputImage,TMaskImage>
-::ThreadedGenerateData(const RegionType&, itk::ThreadIdType threadid)
+::ThreadedGenerateVectorData(const ogr::Layer& layerForThread, itk::ThreadIdType threadid)
 {
   // Retrieve inputs
   TInputImage* inputImage = const_cast<TInputImage*>(this->GetInput());
   TInputImage* outputImage = this->GetOutput();
   RegionType requestedRegion = outputImage->GetRequestedRegion();
 
-  ogr::Layer layer = this->m_InMemoryInputs.at(threadid)->GetLayerChecked(0);
-  if (! layer)
-    {
-    return;
-    }
-
-  itk::ProgressReporter progress( this, threadid, layer.GetFeatureCount(true) );
+  itk::ProgressReporter progress( this, threadid, layerForThread.GetFeatureCount(true) );
 
   // Loop across the features in the layer (filtered by requested region in BeforeTGD already)
-  ogr::Layer::const_iterator featIt = layer.begin();
-  for(; featIt!=layer.end(); ++featIt)
+  ogr::Layer::const_iterator featIt = layerForThread.begin();
+  for(; featIt!=layerForThread.end(); ++featIt)
     {
     // Compute the intersection of thread region and polygon bounding region, called "considered region"
     // This need not be done in ThreadedGenerateData and could be pre-processed and cached before filter execution if needed
@@ -603,7 +683,7 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
 template<class TInputImage, class TMaskImage>
 void
 PersistentSamplingFilterBase<TInputImage,TMaskImage>
-::PrepareInputVectors()
+::DispatchInputVectors()
 {
   TInputImage* outputImage = this->GetOutput();
   ogr::DataSource* vectors = const_cast<ogr::DataSource*>(this->GetOGRData());
@@ -634,52 +714,13 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
   inLayer.SetSpatialFilter(&tmpPolygon);
 
   unsigned int numberOfThreads = this->GetNumberOfThreads();
-
-  unsigned int actualNumberOfThreads = numberOfThreads;
-
-  if(numberOfThreads > this->GetOutput()->GetRequestedRegion().GetSize()[1])
-    {
-    actualNumberOfThreads = this->GetOutput()->GetRequestedRegion().GetSize()[1];
-    }
-  
-  // prepare temporary input : split input features between available threads
-  this->m_InMemoryInputs.clear();
-  std::string tmpLayerName("thread");
-  OGRSpatialReference * oSRS = ITK_NULLPTR;
-  if (inLayer.GetSpatialRef())
-    {
-    oSRS = inLayer.GetSpatialRef()->Clone();
-    }
-  OGRFeatureDefn &layerDefn = inLayer.GetLayerDefn();
   std::vector<ogr::Layer> tmpLayers;
-  for (unsigned int i=0 ; i < actualNumberOfThreads ; i++)
+  tmpLayers.reserve(numberOfThreads);
+  for (unsigned int i=0 ; i<numberOfThreads ; i++)
     {
-    ogr::DataSource::Pointer tmpOgrDS = ogr::DataSource::New();
-    ogr::Layer tmpLayer = tmpOgrDS->CreateLayer(
-      tmpLayerName,
-      oSRS,
-      inLayer.GetGeomType());
-    // add field definitions
-    for (int k=0 ; k < layerDefn.GetFieldCount() ; k++)
-      {
-      OGRFieldDefn originDefn(layerDefn.GetFieldDefn(k));
-      ogr::FieldDefn fieldDefn(originDefn);
-      tmpLayer.CreateField(fieldDefn);
-      }
-    this->m_InMemoryInputs.push_back(tmpOgrDS);
-    tmpLayers.push_back(tmpLayer);
+    tmpLayers.push_back(this->GetInMemoryInput(i));
     }
-
-  this->DispatchInputVectors(inLayer,tmpLayers);
-
-  inLayer.SetSpatialFilter(ITK_NULLPTR);
-}
-
-template<class TInputImage, class TMaskImage>
-void
-PersistentSamplingFilterBase<TInputImage,TMaskImage>
-::DispatchInputVectors(ogr::Layer &inLayer, std::vector<ogr::Layer> &tmpLayers)
-{
+  
   OGRFeatureDefn &layerDefn = inLayer.GetLayerDefn();
   ogr::Layer::const_iterator featIt = inLayer.begin();
   unsigned int counter=0;
@@ -693,58 +734,8 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
     if (counter >= tmpLayers.size())
       counter = 0;
     }
-}
 
-template<class TInputImage, class TMaskImage>
-void
-PersistentSamplingFilterBase<TInputImage,TMaskImage>
-::PrepareOutputVectors()
-{
-  // Prepare in-memory outputs
-  unsigned int numberOfThreads = this->GetNumberOfThreads();
-
-  unsigned int actualNumberOfThreads = numberOfThreads;
-
-  if(numberOfThreads > this->GetOutput()->GetRequestedRegion().GetSize()[1])
-    {
-    actualNumberOfThreads = this->GetOutput()->GetRequestedRegion().GetSize()[1];
-    }
-  
-  this->m_InMemoryOutputs.clear();
-  std::string tmpLayerName("threadOut");
-  for (unsigned int i=0 ; i < actualNumberOfThreads ; i++)
-    {
-    std::vector<OGRDataPointer> tmpContainer;
-    // iterate over outputs, only process ogr::DataSource
-    for (unsigned int k=0 ; k < this->GetNumberOfOutputs() ; k++)
-      {
-      ogr::DataSource* realOutput = dynamic_cast<ogr::DataSource *>(
-        this->itk::ProcessObject::GetOutput(k));
-      if (realOutput)
-        {
-        ogr::Layer realLayer = realOutput->GetLayersCount() == 1
-                               ? realOutput->GetLayer(0)
-                               : realOutput->GetLayer(m_OutLayerName);
-        OGRSpatialReference * oSRS = ITK_NULLPTR;
-        if (realLayer.GetSpatialRef())
-          {
-          oSRS = realLayer.GetSpatialRef()->Clone();
-          }
-        OGRFeatureDefn &layerDefn = realLayer.GetLayerDefn();
-        ogr::DataSource::Pointer tmpOutput = ogr::DataSource::New();
-        ogr::Layer tmpLayer = tmpOutput->CreateLayer(
-          tmpLayerName, oSRS,  realLayer.GetGeomType());
-        // add field definitions
-        for (int f=0 ; f < layerDefn.GetFieldCount() ; f++)
-          {
-          OGRFieldDefn originDefn(layerDefn.GetFieldDefn(f));
-          tmpLayer.CreateField(originDefn);
-          }
-        tmpContainer.push_back(tmpOutput);
-        }
-      }
-    this->m_InMemoryOutputs.push_back(tmpContainer);
-    }
+  inLayer.SetSpatialFilter(ITK_NULLPTR);
 }
 
 template<class TInputImage, class TMaskImage>
@@ -857,6 +848,53 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
   return this->m_AdditionalFields;
 }
 
+template<class TInputImage, class TMaskImage>
+ITK_THREAD_RETURN_TYPE
+PersistentSamplingFilterBase<TInputImage,TMaskImage>
+::VectorThreaderCallback(void *arg)
+{
+  VectorThreadStruct *str = (VectorThreadStruct*)(((itk::MultiThreader::ThreadInfoStruct *)(arg))->UserData);
+
+  int threadId = ((itk::MultiThreader::ThreadInfoStruct *)(arg))->ThreadID;
+  int threadCount = ((itk::MultiThreader::ThreadInfoStruct *)(arg))->NumberOfThreads;
+
+  ogr::Layer layer = str->Filter->GetInMemoryInput(threadId);
+
+  if (threadId < threadCount)
+    {
+    str->Filter->ThreadedGenerateVectorData(layer,threadId);
+    }
+
+  return ITK_THREAD_RETURN_VALUE;
+}
+
+template<class TInputImage, class TMaskImage>
+ogr::Layer
+PersistentSamplingFilterBase<TInputImage,TMaskImage>
+::GetInMemoryInput(unsigned int threadId)
+{
+  if (threadId >= m_InMemoryInputs.size())
+    {
+    itkExceptionMacro(<< "Requested in-memory input layer not available " << threadId << " (total size : "<< m_InMemoryInputs.size() <<").");
+    }
+  return m_InMemoryInputs[threadId]->GetLayerChecked(0);
+}
+
+template<class TInputImage, class TMaskImage>
+ogr::Layer
+PersistentSamplingFilterBase<TInputImage,TMaskImage>
+::GetInMemoryOutput(unsigned int threadId, unsigned int index)
+{
+  if (threadId >= m_InMemoryOutputs.size())
+    {
+    itkExceptionMacro(<< "Requested in-memory output layer not available " << threadId << " (total size : "<< m_InMemoryOutputs.size() <<").");
+    }
+  if (index >= m_InMemoryOutputs[threadId].size())
+    {
+    itkExceptionMacro(<< "Requested output dataset not available " << index << " (available : "<< m_InMemoryOutputs[threadId].size() <<").");
+    }
+  return m_InMemoryOutputs[threadId][index]->GetLayerChecked(0);
+}
 
 } // end namespace otb