[depthMap] Add similarity volumes debug export

src/aliceVision/depthMap/SemiGlobalMatchingRc.cpp

@@ -475,16 +475,31 @@ bool SemiGlobalMatchingRc::sgmrc(bool checkIfExists)
         delete simVolume;
     }
 
+    //if(_sp->exportIntermediateResults)
+    //  svol->exportVolume(*_depths, _rc, _scale, _step, _sp->mp->getDepthMapsFolder() + std::to_string(_sp->mp->getViewId(_rc)) + "_vol_beforeReduction.abc");
+
     // reduction of 'volume' (X, Y, Z) into 'volumeStepZ' (X, Y, Z/step)
     svol->cloneVolumeSecondStepZ();
 
+    if(_sp->exportIntermediateResults)
+    {
+      svol->exportVolumeStep(*_depths, _rc, _scale, _step, _sp->mp->getDepthMapsFolder() + std::to_string(_sp->mp->getViewId(_rc)) + "_vol_afterReduction.abc");
+      svol->export9PCSV(*_depths, _rc, _scale, _step, "afterReduction", _sp->mp->getDepthMapsFolder() + std::to_string(_sp->mp->getViewId(_rc)) + "_9p.csv");
+    }
+
     // filter on the 3D volume to weight voxels based on their neighborhood strongness.
     // so it downweights local minimums that are not supported by their neighborhood.
     // this is here for experimental reason ... to show how SGGC work on non
     // optimized depthmaps ... it must equals to true in normal case
     if(_sp->doSGMoptimizeVolume)
         svol->SGMoptimizeVolumeStepZ(_rc, _step, 0, 0, _scale);
 
+    if(_sp->exportIntermediateResults)
+    {
+      svol->exportVolumeStep(*_depths, _rc, _scale, _step, _sp->mp->getDepthMapsFolder() + std::to_string(_sp->mp->getViewId(_rc)) + "_vol_afterFiltering.abc");
+      svol->export9PCSV(*_depths, _rc, _scale, _step, "afterFiltering", _sp->mp->getDepthMapsFolder() + std::to_string(_sp->mp->getViewId(_rc)) + "_9p.csv");
+    }
+
     // for each pixel: choose the voxel with the minimal similarity value
     const int zborder = 2;
     _volumeBestIdVal = svol->getOrigVolumeBestIdValFromVolumeStepZ(zborder);

src/aliceVision/depthMap/SemiGlobalMatchingVolume.cpp

@@ -7,8 +7,12 @@
 #include "SemiGlobalMatchingVolume.hpp"
 #include <aliceVision/system/Logger.hpp>
 #include <aliceVision/mvsData/Point3d.hpp>
+#include <aliceVision/mvsData/jetColorMap.hpp>
 #include <aliceVision/mvsUtils/common.hpp>
 
+#include <aliceVision/sfmData/SfMData.hpp>
+#include <aliceVision/sfmDataIO/sfmDataIO.hpp>
+
 namespace aliceVision {
 namespace depthMap {
 
@@ -137,6 +141,122 @@ void SemiGlobalMatchingVolume::cloneVolumeSecondStepZ()
         mvsUtils::printfElapsedTime(tall, "SemiGlobalMatchingVolume::cloneVolumeSecondStepZ ");
 }
 
+void SemiGlobalMatchingVolume::exportVolume(StaticVector<float>& depths, int camIndex,  int scale, int step, const std::string& filepath) const
+{
+    sfmData::SfMData pointCloud;
+    const unsigned char* _volumeSecondBestPtr = _volumeSecondBest->getData().data();
+    const mvsUtils::MultiViewParams* mp = sp->mp;
+    int eStep = 10;
+
+    IndexT landmarkId;
+    for(int z = 0; z < volDimZ; ++z)
+    {
+        for(int y = 0; y < volDimY; y+=eStep)
+        {
+            for(int x = 0; x < volDimX; x+=eStep)
+            {
+                const double planeDepth = depths[z];
+                const Point3d planen = (mp->iRArr[camIndex] * Point3d(0.0f, 0.0f, 1.0f)).normalize();
+                const Point3d planep = mp->CArr[camIndex] + planen * planeDepth;
+                const Point3d v = (mp->iCamArr[camIndex] * Point2d(x * scale * step, y * scale * step)).normalize();
+                const Point3d p = linePlaneIntersect(mp->CArr[camIndex], v, planep, planen);
+
+                const int index = z * volDimX * volDimY + y * volDimX + x;
+                const int maxValue = 80;
+                if(_volumeSecondBestPtr[index] > maxValue)
+                  continue;
+                const rgb c = getRGBFromJetColorMap(static_cast<double>(_volumeSecondBestPtr[index]) / double(maxValue));
+                pointCloud.getLandmarks()[landmarkId] = sfmData::Landmark(Vec3(p.x, p.y, p.z), feature::EImageDescriberType::UNKNOWN, sfmData::Observations(), image::RGBColor(c.r, c.g, c.b));
+
+                ++landmarkId;
+            }
+        }
+    }
+
+    sfmDataIO::Save(pointCloud, filepath, sfmDataIO::ESfMData::STRUCTURE);
+}
+
+void SemiGlobalMatchingVolume::exportVolumeStep(StaticVector<float>& depths, int camIndex,  int scale, int step, const std::string& filepath) const
+{
+    sfmData::SfMData pointCloud;
+    const unsigned char* volumePtr = _volumeStepZ->getData().data();
+    const mvsUtils::MultiViewParams* mp = sp->mp;
+    int eStep = 10;
+
+    IndexT landmarkId;
+    for(int stepZ = 0; stepZ < (volDimZ / volStepZ); ++stepZ)
+    {
+        for(int y = 0; y < volDimY; y+=eStep)
+        {
+            for(int x = 0; x < volDimX; x+=eStep)
+            {
+                const int z = (*_volumeBestZ)[stepZ * volDimX * volDimY + y * volDimX + x];
+                const double planeDepth = depths[z];
+                const Point3d planen = (mp->iRArr[camIndex] * Point3d(0.0f, 0.0f, 1.0f)).normalize();
+                const Point3d planep = mp->CArr[camIndex] + planen * planeDepth;
+                const Point3d v = (mp->iCamArr[camIndex] * Point2d(x * scale * step, y * scale * step)).normalize();
+                const Point3d p = linePlaneIntersect(mp->CArr[camIndex], v, planep, planen);
+
+                const int index = stepZ * volDimX * volDimY + y * volDimX + x;
+                const int maxValue = 80;
+                if(volumePtr[index] > maxValue)
+                  continue;
+                const rgb c = getRGBFromJetColorMap(static_cast<double>(volumePtr[index]) / double(maxValue));
+                pointCloud.getLandmarks()[landmarkId] = sfmData::Landmark(Vec3(p.x, p.y, p.z), feature::EImageDescriberType::UNKNOWN, sfmData::Observations(), image::RGBColor(c.r, c.g, c.b));
+
+                ++landmarkId;
+            }
+        }
+    }
+
+    sfmDataIO::Save(pointCloud, filepath, sfmDataIO::ESfMData::STRUCTURE);
+}
+
+void SemiGlobalMatchingVolume::export9PCSV(StaticVector<float>& depths, int camIndex,  int scale, int step, const std::string& name, const std::string& filepath) const
+{
+    const unsigned char* volumePtr = _volumeStepZ->getData().data();
+
+    const int xOffset = std::floor(volDimX / 4.0f);
+    const int yOffset = std::floor(volDimY / 4.0f);
+
+    std::array<std::vector<double>, 9> ptsDepths;
+
+    for(int iy = 0; iy < 3; ++iy)
+    {
+        for(int ix = 0; ix < 3; ++ix)
+        {
+            const int x = (ix + 1) * xOffset;
+            const int y = (iy + 1) * yOffset;
+
+            std::vector<double>& pDepths = ptsDepths.at(iy * 3 + ix);
+
+            for(int stepZ = 0; stepZ < (volDimZ / volStepZ); ++stepZ)
+            {
+                pDepths.push_back(volumePtr[stepZ * volDimX * volDimY + y * volDimX + x]);
+            }
+        }
+    }
+
+    std::stringstream ss;
+    {
+      ss << name << "\n";
+      int ptId = 1;
+      for(const std::vector<double>& pDepths : ptsDepths)
+      {
+        ss << "p" << ptId << ";";
+        for(const double& depth : pDepths)
+          ss << depth << ";";
+        ss << "\n";
+        ++ptId;
+      }
+    }
+
+    std::ofstream file;
+    file.open(filepath, std::ios_base::app);
+    if(file.is_open())
+      file << ss.str();
+}
+
 /**
  * @param[in] volStepXY step in the image space
  */

src/aliceVision/depthMap/SemiGlobalMatchingVolume.hpp

@@ -30,6 +30,10 @@ class SemiGlobalMatchingVolume
     void SGMoptimizeVolumeStepZ(int rc, int volStepXY, int volLUX, int volLUY, int scale);
     StaticVector<IdValue>* getOrigVolumeBestIdValFromVolumeStepZ(int zborder);
 
+    void exportVolume(StaticVector<float>& depths, int camIndex, int scale, int step, const std::string& filepath) const;
+    void exportVolumeStep(StaticVector<float>& depths, int camIndex, int scale, int step, const std::string& filepath) const;
+    void export9PCSV(StaticVector<float>& depths, int camIndex,  int scale, int step, const std::string& name, const std::string& filepath) const;
+
 private:
     SemiGlobalMatchingParams* sp;
 

src/aliceVision/fuseCut/DelaunayGraphCut.cpp

@@ -325,7 +325,7 @@ void createVerticesWithVisibilities(const StaticVector<int>& cams, std::vector<P
                 if(depth <= 0.0f)
                     continue;
 
-                const Point3d p = mp->CArr[c] + (mp->iCamArr[c] * Point2d((float)x, (float)y)).normalize() * depth;
+                const Point3d p = mp->backproject(c, Point2d(x, y), depth);
                 const double pixSize = mp->getCamPixelSize(p, c);
 #ifdef USE_GEOGRAM_KDTREE
                 const std::size_t nearestVertexIndex = kdTree.get_nearest_neighbor(p.m);

src/aliceVision/mvsUtils/MultiViewParams.hpp

@@ -109,6 +109,11 @@ class MultiViewParams
 
     ~MultiViewParams();
 
+    inline Point3d backproject(const int camIndex, const Point2d& pix, double depth) const
+    {
+        const Point3d p = CArr[camIndex] + (iCamArr[camIndex] * pix).normalize() * depth;
+        return p;
+    }
     inline const std::string& getImagePath(int index) const
     {
         return _imagesParams.at(index).path;