ENH: JIRA 352: work on left-right disparity bijection

Code/DisparityMap/otbBijectionCoherencyFilter.h

@@ -62,9 +62,10 @@ public:
   typedef TDisparityImage  DispMapType;
   typedef TOutputImage     MaskType;
   
-  typedef TOutputImage::RegionType      OutputRegionType;
-  typedef TDisparityImage::RegionType   InputRegionType;
-  typedef TDisparityImage::SizeType     SizeType;
+  typedef typename MaskType::RegionType      OutputRegionType;
+  typedef typename DispMapType::RegionType   InputRegionType;
+  typedef typename DispMapType::SizeType     SizeType;
+  typedef typename DispMapType::IndexType    IndexType;
   
   /** Set the direct horizontal disparity map */
   void SetDirectHorizontalDisparityMapInput( const TDisparityImage * hmap );
@@ -89,6 +90,19 @@ public:
   
   /** Get the tolerance radius */
   itkGetMacro(Tolerance,double);
+  
+  /** Set/Get macro for exploration area */
+  itkSetMacro(MinHDisp,int);
+  itkGetMacro(MinHDisp,int);
+  
+  itkSetMacro(MaxHDisp,int);
+  itkGetMacro(MaxHDisp,int);
+  
+  itkSetMacro(MinVDisp,int);
+  itkGetMacro(MinVDisp,int);
+  
+  itkSetMacro(MaxVDisp,int);
+  itkGetMacro(MaxVDisp,int);
 
 protected:
   /** Constructor */
@@ -113,6 +127,17 @@ private:
   /** Tolerance radius (in pixels) */
   double m_Tolerance;
   
+  /** Minimum horizontal input disparity  */
+  int m_MinHDisp;
+  
+  /** Maximum horizontal input disparity */
+  int m_MaxHDisp;
+  
+  /** Minimum vertical input disparity */
+  int m_MinVDisp;
+  
+  /** Maximum vertical input disparity */
+  int m_MaxVDisp;
 };
 
 } // end namespace otb

Code/DisparityMap/otbBijectionCoherencyFilter.txx

@@ -34,6 +34,10 @@ BijectionCoherencyFilter<TDisparityImage,TOutputImage>
   this->SetNumberOfInputs(4);
   this->SetNumberOfRequiredInputs(1);
   this->m_Tolerance = 1.;
+  this->m_MinHDisp = -5;
+  this->m_MaxHDisp = 5;
+  this->m_MinVDisp = -5;
+  this->m_MaxVDisp = 5;
 
   // Set the outputs
   this->SetNumberOfOutputs(1);
@@ -135,6 +139,9 @@ BijectionCoherencyFilter<TDisparityImage,TOutputImage>
   const TDisparityImage * directHmap = this->GetDirectHorizontalDisparityMapInput();
   const TDisparityImage * reverseHmap = this->GetReverseHorizontalDisparityMapInput();
   
+  const TDisparityImage * directVmap = this->GetDirectVerticalDisparityMapInput();
+  const TDisparityImage * reverseVmap = this->GetReverseVerticalDisparityMapInput();
+  
   if (!directHmap)
     {
     itkExceptionMacro(<<"Direct horizontal disparity map is missing");
@@ -145,8 +152,24 @@ BijectionCoherencyFilter<TDisparityImage,TOutputImage>
     itkExceptionMacro(<<"Reverse horizontal disparity map is missing");
     }
   
-  // TODO : check sizes
+  if (directVmap && directVmap->GetLargestPossibleRegion() != directHmap->GetLargestPossibleRegion())
+    {
+    itkExceptionMacro(<<"Horizontal and vertical direct disparity maps have different sizes.");
+    }
+  
+  if (reverseVmap && reverseVmap->GetLargestPossibleRegion() != reverseHmap->GetLargestPossibleRegion())
+    {
+    itkExceptionMacro(<<"Horizontal and vertical reverse disparity maps have different sizes.");
+    }
   
+  if (this->m_MinHDisp > this->m_MaxHDisp)
+    {
+    itkExceptionMacro(<<"Wrong horizontal exploration values");
+    }
+  if (this->m_MinVDisp > this->m_MaxVDisp)
+    {
+    itkExceptionMacro(<<"Wrong horizontal exploration values");
+    }
 }
 
 template <class TDisparityImage, class TOutputImage>
@@ -157,26 +180,116 @@ BijectionCoherencyFilter<TDisparityImage,TOutputImage>
   this->Superclass::GenerateInputRequestedRegion();
   
   OutputRegionType requested = this->GetOutput()->GetRequestedRegion();
-  InputRegionType  largest = this->GetDirectHorizontalDisparityMap()->GetLargestPossibleRegion();
-  InputRegionType  inputRequested;
-  this->CallCopyOutputRegionToInputRegion(inputRequested,requested);
+  InputRegionType  directLargest = this->GetDirectHorizontalDisparityMap()->GetLargestPossibleRegion();
+  InputRegionType  directRequested;
+  InputRegionType  reverseLargest = this->GetReverseHorizontalDisparityMap()->GetLargestPossibleRegion();
+  InputRegionType  reverseRequested;
   
-  // TODO : the requested region on Reverse dispmap must be shifted with an estimate of the min/max disparities values
-  // HOW CAN I DO THAT ???
+  this->CallCopyOutputRegionToInputRegion(requested,directRequested);
   
-  SizeType radius;
-  radius.Fill(static_cast<unsigned int>(vcl_ceil(this->m_Tolerance)));
+  reverseRequested.SetIndex(0,requested.GetIndex(0) + this->m_MinHDisp);
+  reverseRequested.SetIndex(1,requested.GetIndex(1) + this->m_MinVDisp);
+  reverseRequested.SetSize(0,requested.GetSize(0) + this->m_MaxHDisp - this->m_MinHDisp);
+  reverseRequested.SetSize(1,requested.GetSize(1) + this->m_MaxVDisp - this->m_MinVDisp);
   
-  inputRequested.PadByRadius(radius);
-  
-  inputRequested.Crop(largest);
+  reverseRequested.Crop(reverseLargest);
   
   TDisparityImage * directHmap = const_cast<TDisparityImage *>(this->GetDirectHorizontalDisparityMapInput());
   TDisparityImage * directVmap = const_cast<TDisparityImage *>(this->GetDirectVerticalDisparityMapInput());
+  TDisparityImage * reverseHmap = const_cast<TDisparityImage *>(this->GetDirectHorizontalDisparityMapInput());
+  TDisparityImage * reverseVmap = const_cast<TDisparityImage *>(this->GetDirectVerticalDisparityMapInput());
   
+  directHmap->SetRequestedRegion(directRequested);
+  if (directVmap) directVmap->SetRequestedRegion(directRequested);
   
+  reverseHmap->SetRequestedRegion(reverseRequested);
+  if (reverseVmap) reverseVmap->SetRequestedRegion(reverseRequested);
 }
 
+template <class TDisparityImage, class TOutputImage>
+void
+BijectionCoherencyFilter<TDisparityImage,TOutputImage>
+::ThreadedGenerateData(const RegionType & outputRegionForThread, int threadId)
+{
+  const TDisparityImage * directHmap = this->GetDirectHorizontalDisparityMapInput();
+  const TDisparityImage * directVmap = this->GetDirectVerticalDisparityMapInput();
+  const TDisparityImage * reverseHmap = this->GetDirectHorizontalDisparityMapInput();
+  const TDisparityImage * reverseVmap = this->GetDirectVerticalDisparityMapInput();
+  
+  TOutputImage * output = this->GetOutput();
+  
+  InputRegionType buffered = reverseHmap->GetBufferedRegion();
+  
+  typedef itk::ImageRegionIterator<TOutputImage> MaskIteratorType;
+  MaskIteratorType outIter = MaskIteratorType(output,outputRegionForThread);
+  
+  typedef itk::ImageRegionConstIteratorWithIndex<TDisparityImage> DispIteratorType;
+  DispIteratorType directHorizIter = DispIteratorType(directHmap,outputRegionForThread);
+  
+  DispIteratorType directVertiIter;
+  if (directVmap)
+    {
+    directVertiIter = DispIteratorType(directVmap,outputRegionForThread);
+    directVertiIter.GoToBegin();
+    }
+  
+  outIter.GoToBegin();
+  directHorizIter.GoToBegin();
+  
+  while (!outIter.IsAtEnd())
+  {
+    IndexType startIndex = directHorizIter.GetIndex();
+    itk::ContinuousIndex<double,2> tmpIndex(startIndex);
+    
+    tmpIndex[0] += directHorizIter.Get();
+    if (directVmap) tmpIndex[1] += directVertiIter.Get();
+    
+    // Interpolate in reverse disparity map
+    IndexType ul,ur,ll,lr;
+    ul[0] = static_cast<long>(vcl_floor(tmpIndex[0]));
+    ul[1] = static_cast<long>(vcl_floor(tmpIndex[1]));
+    if (ul[0]<buffered.GetIndex()[0]) ul[0]=buffered.GetIndex()[0];
+    if (ul[1]<buffered.GetIndex()[1]) ul[1]=buffered.GetIndex()[1];
+    if (ul[0]>(buffered.GetIndex()[0]+buffered.GetSize()[0]-1)) ul[0]=(buffered.GetIndex()[0]+buffered.GetSize()[0]-1);
+    if (ul[1]>(buffered.GetIndex()[1]+buffered.GetSize()[1]-1)) ul[1]=(buffered.GetIndex()[1]+buffered.GetSize()[1]-1);
+  
+    ur = ul;
+    ur[0] += 1;
+    ll = ul;
+    ll[1] += 1;
+    lr = ul;
+    lr[0] += 1;
+    lr[1] += 1;
+    
+    double rx = tmpIndex[0] - static_cast<double>(ul[0]);
+    double ry = tmpIndex[1] - static_cast<double>(ul[1]);
+    
+    itk::ContinuousIndex<double,2> backIndex(tmpIndex);
+    backIndex[0] += (1. - ry) * ((1. - rx) * reverseHmap->GetPixel(ul) + rx * reverseHmap->GetPixel(ur)) +
+                            ry * ((1. - rx) * reverseHmap->GetPixel(ll) + rx * reverseHmap->GetPixel(lr));
+    if (reverseVmap)
+      {
+      backIndex[1] += (1. - ry) * ((1. - rx) * reverseVmap->GetPixel(ul) + rx * reverseVmap->GetPixel(ur)) +
+                          ry * ((1. - rx) * reverseVmap->GetPixel(ll) + rx * reverseVmap->GetPixel(lr));
+      }
+    
+    if (vcl_abs(backIndex[0] - static_cast<double>(startIndex[0]))< this->m_Tolerance &&
+        vcl_abs(backIndex[1] - static_cast<double>(startIndex[1]))< this->m_Tolerance)
+      {
+      outIter.Set(255);
+      }
+    else
+      {
+      outIter.Set(0);
+      }
+    
+    ++outIter;
+    ++directHorizIter;
+    if (directVmap) ++directVertiIter;
+  }
+  
+  
+}
 
 }