Skip to content

Commit

Permalink
Merge pull request #832 from alicevision/dev/ImproveDelaunayGraphCut
Browse files Browse the repository at this point in the history 
[fuseCut] improve delaunayGraphCut
  • Loading branch information
@fabiencastan
fabiencastan authored Jul 29, 2020
2 parents c995e01 + de2a440 commit aad5edc
Show file tree
Hide file tree
Showing 12 changed files with 566 additions and 1,001 deletions.
1,059 changes: 445 additions & 614 deletions src/aliceVision/fuseCut/DelaunayGraphCut.cpp
View file

Large diffs are not rendered by default.

140 changes: 45 additions & 95 deletions src/aliceVision/fuseCut/DelaunayGraphCut.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,17 +97,37 @@ class DelaunayGraphCut
DelaunayGraphCut(mvsUtils::MultiViewParams* _mp);
virtual ~DelaunayGraphCut();

/// Get absolute opposite vertex index
/**
* @brief Retrieve the global vertex index of the localVertexIndex of the facet.
*
* @param f the facet
* @return the global vertex index
*/
inline VertexIndex getOppositeVertexIndex(const Facet& f) const
{
return _tetrahedralization->cell_vertex(f.cellIndex, f.localVertexIndex);
}
/// Get absolute vertex index

/**
* @brief Retrieve the global vertex index of a vertex from a facet and an relative index
* compared to the localVertexIndex of the facet.
*
* @param f the facet
* @param i the relative index (relative to the localVertexIndex of the facet)
* @return the global vertex index
*/
inline VertexIndex getVertexIndex(const Facet& f, int i) const
{
return _tetrahedralization->cell_vertex(f.cellIndex, ((f.localVertexIndex + i + 1) % 4));
}

inline const std::array<const Point3d*, 3> getFacetsPoints(const Facet& f) const
{
return {&(_verticesCoords[getVertexIndex(f, 0)]),
&(_verticesCoords[getVertexIndex(f, 1)]),
&(_verticesCoords[getVertexIndex(f, 2)])};
}

inline std::size_t getNbVertices() const
{
return _verticesAttr.size();
Expand Down Expand Up @@ -234,6 +254,17 @@ class DelaunayGraphCut
return localCells[lvi];
}

/**
* @brief Retrieves the global indexes of neighboring cells using the global index of a vertex.
*
* @param vi the global vertexIndex
* @return a vector of neighboring cell indexes
*/
inline const std::vector<CellIndex>& getNeighboringCellsByVertexIndex(VertexIndex vi) const
{
return _neighboringCellsPerVertex.at(vi);
}

void initVertices();
void computeDelaunay();
void initCells();
Expand All @@ -251,32 +282,23 @@ class DelaunayGraphCut

void addPointsFromSfM(const Point3d hexah[8], const StaticVector<int>& cams, const sfmData::SfMData& sfmData);
void addPointsFromCameraCenters(const StaticVector<int>& cams, float minDist);
void addPointsToPreventSingularities(const Point3d Voxel[], float minDist);
void addPointsToPreventSingularities(const Point3d Voxel[8], float minDist);

/**
* @brief Add volume points to prevent singularities
*/
void addHelperPoints(int nGridHelperVolumePointsDim, const Point3d Voxel[], float minDist);
void addHelperPoints(int nGridHelperVolumePointsDim, const Point3d Voxel[8], float minDist);

void fuseFromDepthMaps(const StaticVector<int>& cams, const Point3d voxel[8], const FuseParams& params);
void loadPrecomputedDensePoints(const StaticVector<int>* voxelsIds, const Point3d voxel[8], VoxelsGrid* ls);

void createTetrahedralizationFromDepthMapsCamsVoxel(const StaticVector<int>& allCams,
StaticVector<int>* voxelsIds, Point3d Voxel[8], VoxelsGrid* ls);

void computeVerticesSegSize(bool allPoints, float alpha = 0.0f);
void removeSmallSegs(int minSegSize);

bool rayCellIntersection(const Point3d &camC, const Point3d &p, int c, Facet& out_facet, bool nearestFarest,
Point3d& out_nlpi) const;
inline bool nearestNeighCellToTheCamOnTheRay(const Point3d &camC, Point3d& out_p, int tetrahedron, Facet &out_f1, Facet &f2,
Point3d& out_lpi) const;
inline bool farestNeighCellToTheCamOnTheRay(Point3d& camC, Point3d& p, int tetrahedron, Facet &f1, Facet &f2,
Point3d& lpi) const;
inline Facet getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex, int cam) const;
Facet getFacetInFrontVertexOnTheRayToThePoint3d(VertexIndex vi, Point3d& ptt) const;
Facet getFacetBehindVertexOnTheRayToTheCam(VertexIndex vi, int cam) const;
int getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& lpi) const;
bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, const Facet& inFacet, Facet& outFacet,
bool nearestFarest, Point3d& outIntersectPt) const;

Facet getFacetFromVertexOnTheRayToTheCam(VertexIndex globalVertexIndex, int cam, bool nearestFarest) const;
Facet getFirstFacetOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const;

float distFcn(float maxDist, float dist, float distFcnHeight) const;

Expand All @@ -285,24 +307,17 @@ class DelaunayGraphCut
double maxEdgeLength() const;
Point3d cellCircumScribedSphereCentre(CellIndex ci) const;
double getFaceWeight(const Facet &f1) const;
float weightFromSim(float sim);

float weightFcn(float nrc, bool labatutWeights, int ncams);

virtual void fillGraph(bool fixesSigma, float nPixelSizeBehind, bool allPoints, bool behind, bool labatutWeights,
virtual void fillGraph(bool fixesSigma, float nPixelSizeBehind, bool labatutWeights,
bool fillOut, float distFcnHeight = 0.0f);
void fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBehind, int vertexIndex, int cam, float weight,
bool fixesSigma, float nPixelSizeBehind, bool allPoints, bool behind, bool fillOut,
bool fixesSigma, float nPixelSizeBehind, bool fillOut,
float distFcnHeight);

void forceTedgesByGradientCVPR11(bool fixesSigma, float nPixelSizeBehind);
void forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSizeBehind);

void updateGraphFromTmpPtsCamsHexah(const StaticVector<int>& incams, Point3d hexah[8], std::string tmpCamsPtsFolderName,
bool labatutWeights, float distFcnHeight = 0.0f);
void updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8], std::string tmpCamsPtsFolderName,
bool labatutWeights, float distFcnHeight);

int setIsOnSurface();

void addToInfiniteSw(float sW);
Expand All @@ -313,16 +328,11 @@ class DelaunayGraphCut

void maxflow();

void reconstructExpetiments(const StaticVector<int>& cams, const std::string& folderName,
bool update, Point3d hexahInflated[8], const std::string& tmpCamsPtsFolderName,
const Point3d& spaceSteps);

void voteFullEmptyScore(const StaticVector<int>& cams, const std::string& folderName);

void createDensePointCloud(Point3d hexah[8], const StaticVector<int>& cams, const sfmData::SfMData* sfmData, const FuseParams* depthMapsFuseParams);
void createDensePointCloudFromPrecomputedDensePoints(Point3d hexah[8], const StaticVector<int>& cams, StaticVector<int>* voxelsIds, VoxelsGrid* ls);

void createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, VoxelsGrid* ls, const std::string& folderName, const std::string& tmpCamsPtsFolderName,
bool removeSmallSegments, const Point3d& spaceSteps);
void createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments);

/**
* @brief Invert full/empty status of cells if they represent a too small group after labelling.
Expand All @@ -331,74 +341,14 @@ class DelaunayGraphCut

void graphCutPostProcessing();

void clearAllPointsInFreeSpace();
void clearAllPointsNotOnSurface();
void addNewPointsToOccupiedSpace();
void clearOutAddIn();
StaticVector<int>* getNearestTrisFromMeshTris(mesh::Mesh* otherMesh);

void segmentFullOrFree(bool full, StaticVector<int>** inColors, int& nsegments);
int removeBubbles();
int removeDust(int minSegSize);
void leaveLargestFullSegmentOnly();

mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true);
mesh::Mesh* createTetrahedralMesh() const;
};


inline bool DelaunayGraphCut::nearestNeighCellToTheCamOnTheRay(const Point3d& camC, Point3d& out_p, int tetrahedron, Facet& out_f1,
Facet& out_f2, Point3d& out_lpi) const
{
out_f1.cellIndex = GEO::NO_CELL;
out_f1.localVertexIndex = GEO::NO_VERTEX;
out_f2.cellIndex = GEO::NO_CELL;
out_f2.localVertexIndex = GEO::NO_VERTEX;
out_lpi = out_p;

if(rayCellIntersection(camC, out_p, tetrahedron, out_f1, true, out_lpi) == true)
{
out_f2 = mirrorFacet(out_f1);
out_p = out_lpi;
return true;
}
return false;
}

inline bool DelaunayGraphCut::farestNeighCellToTheCamOnTheRay(Point3d& camC, Point3d& p, int tetrahedron, Facet& f1,
Facet& f2, Point3d& lpi) const
{
f1.cellIndex = GEO::NO_CELL;
f1.localVertexIndex = GEO::NO_VERTEX;
f2.cellIndex = GEO::NO_CELL;
f2.localVertexIndex = GEO::NO_VERTEX;
lpi = p;

if(rayCellIntersection(camC, p, tetrahedron, f1, false, lpi) == true)
{
f2 = mirrorFacet(f1);
p = lpi;
return true;
}

return false;
}

inline DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex,
int cam) const
{
// if (btest) {
// printf("cam %i pt %f %f %f\n",cam,p.x,p.y,p.z);
// printf("ptid %i\n",vp->info().id);
// printf("campos %f %f
//%f\n",mp->CArr[cam].x,mp->CArr[cam].y,mp->CArr[cam].z);
//};
if((cam < 0) || (cam >= mp->ncams))
{
ALICEVISION_LOG_WARNING("Bad camId, cam: " << cam << ", ptid: " << vertexIndex);
}

return getFacetInFrontVertexOnTheRayToThePoint3d(vertexIndex, mp->CArr[cam]);
}

} // namespace fuseCut
} // namespace aliceVision
29 changes: 20 additions & 9 deletions src/aliceVision/fuseCut/Fuser.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -629,12 +629,24 @@ void Fuser::divideSpaceFromSfM(const sfmData::SfMData& sfmData, Point3d* hexah,
accMaxZ(z);
}

const double xMin = quantile(accMinX, quantile_probability = 1.0 - percentile);
const double yMin = quantile(accMinY, quantile_probability = 1.0 - percentile);
const double zMin = quantile(accMinZ, quantile_probability = 1.0 - percentile);
const double xMax = quantile(accMaxX, quantile_probability = percentile);
const double yMax = quantile(accMaxY, quantile_probability = percentile);
const double zMax = quantile(accMaxZ, quantile_probability = percentile);
// Remove a percentile of the observations (to remove unstable points)
double xMin = quantile(accMinX, quantile_probability = 1.0 - percentile);
double yMin = quantile(accMinY, quantile_probability = 1.0 - percentile);
double zMin = quantile(accMinZ, quantile_probability = 1.0 - percentile);
double xMax = quantile(accMaxX, quantile_probability = percentile);
double yMax = quantile(accMaxY, quantile_probability = percentile);
double zMax = quantile(accMaxZ, quantile_probability = percentile);

// Add a margin on the result
const double xMargin = (xMax - xMin) * 0.05;
const double yMargin = (yMax - yMin) * 0.05;
const double zMargin = (zMax - zMin) * 0.05;
xMin -= xMargin;
yMin -= yMargin;
zMin -= zMargin;
xMax += xMargin;
yMax += yMargin;
zMax += zMargin;

hexah[0] = Point3d(xMax, yMax, zMax);
hexah[1] = Point3d(xMin, yMax, zMax);
Expand Down Expand Up @@ -766,7 +778,7 @@ std::string generateTempPtsSimsFiles(std::string tmpDir, mvsUtils::MultiViewPara
}

int nnoisePts = ((percNoisePts / 100.0f) * (float)(idsAlive->size()));
StaticVector<int>* randIdsAlive = mvsUtils::createRandomArrayOfIntegers(idsAlive->size());
const std::vector<int> randIdsAlive = mvsUtils::createRandomArrayOfIntegers(idsAlive->size());

srand(time(nullptr));

Expand All @@ -775,7 +787,7 @@ std::string generateTempPtsSimsFiles(std::string tmpDir, mvsUtils::MultiViewPara
for(int x = 0; x < w; ++x)
{
int id = y * w + x;
int i = (*idsAlive)[(*randIdsAlive)[id]];
int i = (*idsAlive)[randIdsAlive[id]];
double depth = depthMap[i];

double sim = simMap[i];
Expand Down Expand Up @@ -823,7 +835,6 @@ std::string generateTempPtsSimsFiles(std::string tmpDir, mvsUtils::MultiViewPara
mvsUtils::printfElapsedTime(t1);

delete idsAlive;
delete randIdsAlive;
}
else
{
Expand Down
71 changes: 0 additions & 71 deletions src/aliceVision/fuseCut/ReconstructionPlan.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,73 +248,6 @@ void ReconstructionPlan::getHexahedronForID(float dist, int id, Point3d* out)
mvsUtils::inflateHexahedron(&(*voxels)[id * 8], out, dist);
}

void reconstructSpaceAccordingToVoxelsArray(const std::string& voxelsArrayFileName, LargeScale* ls)
{
StaticVector<Point3d>* voxelsArray = loadArrayFromFile<Point3d>(voxelsArrayFileName);

ReconstructionPlan* rp =
new ReconstructionPlan(ls->dimensions, &ls->space[0], ls->mp, ls->spaceVoxelsFolderName);

StaticVector<Point3d>* hexahsToExcludeFromResultingMesh = new StaticVector<Point3d>();
hexahsToExcludeFromResultingMesh->reserve(voxelsArray->size());

for(int i = 0; i < voxelsArray->size() / 8; i++)
{
ALICEVISION_LOG_INFO("Reconstructing " << (voxelsArray->size() / 8) << "-th Voxel of " << i << ".");

const std::string folderName = ls->getReconstructionVoxelFolder(i);
bfs::create_directory(folderName);

const std::string meshBinFilepath = folderName + "mesh.bin";
if(!mvsUtils::FileExists(meshBinFilepath))
{
StaticVector<int>* voxelsIds = rp->voxelsIdsIntersectingHexah(&(*voxelsArray)[i * 8]);
DelaunayGraphCut delaunayGC(ls->mp);
Point3d* hexah = &(*voxelsArray)[i * 8];

StaticVector<int> cams = ls->mp->findCamsWhichIntersectsHexahedron(hexah);

if(cams.empty())
throw std::logic_error("No camera to make the reconstruction");

delaunayGC.createDensePointCloudFromPrecomputedDensePoints(hexah, cams, voxelsIds, (VoxelsGrid*)rp);
delaunayGC.createGraphCut(hexah, cams,(VoxelsGrid*)rp, folderName, ls->getSpaceCamsTracksDir(), false, ls->getSpaceSteps());
delete voxelsIds;

// Save mesh as .bin and .obj
mesh::Mesh* mesh = delaunayGC.createMesh();
StaticVector<StaticVector<int>> ptsCams;
delaunayGC.createPtsCams(ptsCams);

mesh::meshPostProcessing(mesh, ptsCams, *ls->mp, folderName, hexahsToExcludeFromResultingMesh, hexah);
mesh->saveToBin(folderName + "mesh.bin");
mesh->saveToObj(folderName + "mesh.obj");

saveArrayOfArraysToFile<int>(folderName + "meshPtsCamsFromDGC.bin", ptsCams);

delete mesh;
}

/*
if(doComputeColoredMeshes)
{
std::string resultFolderName = folderName + "/";
computeColoredMesh(resultFolderName, ls);
}
*/
Point3d hexahThin[8];
mvsUtils::inflateHexahedron(&(*voxelsArray)[i * 8], hexahThin, 0.9);
for(int k = 0; k < 8; k++)
{
hexahsToExcludeFromResultingMesh->push_back(hexahThin[k]);
}
delete hexahsToExcludeFromResultingMesh;
}
delete rp;
delete voxelsArray;
}


StaticVector<StaticVector<int>*>* loadLargeScalePtsCams(const std::vector<std::string>& recsDirs)
{
StaticVector<StaticVector<int>*>* ptsCamsFromDct = new StaticVector<StaticVector<int>*>();
Expand Down Expand Up @@ -499,10 +432,6 @@ mesh::Mesh* joinMeshes(int gl, LargeScale* ls)
{
subFolderName = "LabatutCFG09";
}
if(ls->mp->userParams.get<bool>("global.JancosekCVPR11", true))
{
subFolderName = "JancosekCVPR11";
}
}
subFolderName = subFolderName + "/";

Expand Down
1 change: 0 additions & 1 deletion src/aliceVision/fuseCut/ReconstructionPlan.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,6 @@ class ReconstructionPlan : public VoxelsGrid
};

void reconstructAccordingToOptimalReconstructionPlan(int gl, LargeScale* ls);
void reconstructSpaceAccordingToVoxelsArray(const std::string& voxelsArrayFileName, LargeScale* ls);
mesh::Mesh* joinMeshes(const std::vector<std::string>& recsDirs, StaticVector<Point3d>* voxelsArray, LargeScale* ls);
mesh::Mesh* joinMeshes(int gl, LargeScale* ls);
mesh::Mesh* joinMeshes(const std::string& voxelsArrayFileName, LargeScale* ls);
Expand Down
Loading

0 comments on commit aad5edc

Please sign in to comment.