Point cloud visibility filter issue

[originlake]

Hi,

I noticed there could be a problem in the point cloud visibility filter. The problem is described OpenDroneMap/ODM#1785. I believe the current algorithm cannot handle well for the case like the red shot in the below image, where closer points will occlude further points and leads to unpredictable results in the algorithm. I think introducing normal into consideration could help resolving this problem, I tried to ignore the points that have large angle compared to the shot view direction (smaller than 100° for example) which works for solving the issue in my dataset.

[cdcseacave]

Thx for posting this issue. Pls clarify which is the filtering algorithm you are referring to? Point me to the function or command in OpenMVS

[originlake]

The command is shown below(example), it's executed right after the dense point cloud reconstruction(example)

DensifyPointCloud --filter-point-cloud -20 -i <mvs file>.mvs

Function

openMVS/libs/MVS/SceneDensify.cpp

Lines 2225 to 2359 in ce03889

	// filter point-cloud based on camera-point visibility intersections
	void Scene::PointCloudFilter(int thRemove)
	{
	TD_TIMER_STARTD();
	typedef TOctree<PointCloud::PointArr,PointCloud::Point::Type,3,uint32_t> Octree;
	struct Collector {
	typedef Octree::IDX_TYPE IDX;
	typedef PointCloud::Point::Type Real;
	typedef TCone<Real,3> Cone;
	typedef TSphere<Real,3> Sphere;
	typedef TConeIntersect<Real,3> ConeIntersect;
	Cone cone;
	const ConeIntersect coneIntersect;
	const PointCloud& pointcloud;
	IntArr& visibility;
	PointCloud::Index idxPoint;
	Real distance;
	int weight;
	#ifdef DENSE_USE_OPENMP
	uint8_t pcs[sizeof(CriticalSection)];
	#endif
	Collector(const Cone::RAY& ray, Real angle, const PointCloud& _pointcloud, IntArr& _visibility)
	: cone(ray, angle), coneIntersect(cone), pointcloud(_pointcloud), visibility(_visibility)
	#ifdef DENSE_USE_OPENMP
	{ new(pcs) CriticalSection; }
	~Collector() { reinterpret_cast<CriticalSection*>(pcs)->~CriticalSection(); }
	inline CriticalSection& GetCS() { return *reinterpret_cast<CriticalSection*>(pcs); }
	#else
	{}
	#endif
	inline void Init(PointCloud::Index _idxPoint, const PointCloud::Point& X, int _weight) {
	const Real thMaxDepth(1.02f);
	idxPoint =_idxPoint;
	const PointCloud::Point::EVec D((PointCloud::Point::EVec&)X-cone.ray.m_pOrig);
	distance = D.norm();
	cone.ray.m_vDir = D/distance;
	cone.maxHeight = MaxDepthDifference(distance, thMaxDepth);
	weight = _weight;
	}
	inline bool Intersects(const Octree::POINT_TYPE& center, Octree::Type radius) const {
	return coneIntersect(Sphere(center, radius*Real(SQRT_3)));
	}
	inline void operator() (const IDX* idices, IDX size) {
	const Real thSimilar(0.01f);
	Real dist;
	FOREACHRAWPTR(pIdx, idices, size) {
	const PointCloud::Index idx(*pIdx);
	if (coneIntersect.Classify(pointcloud.points[idx], dist) == VISIBLE && !IsDepthSimilar(distance, dist, thSimilar)) {
	if (dist > distance)
	visibility[idx] += pointcloud.pointViews[idx].size();
	else
	visibility[idx] -= weight;
	}
	}
	}
	};
	typedef CLISTDEF2(Collector) Collectors;
	// create octree to speed-up search
	Octree octree(pointcloud.points, [](Octree::IDX_TYPE size, Octree::Type /*radius*/) {
	return size > 128;
	});
	IntArr visibility(pointcloud.GetSize()); visibility.Memset(0);
	Collectors collectors; collectors.reserve(images.size());
	FOREACH(idxView, images) {
	const Image& image = images[idxView];
	const Ray3f ray(Cast<float>(image.camera.C), Cast<float>(image.camera.Direction()));
	const float angle(float(image.ComputeFOV(0)/image.width));
	collectors.emplace_back(ray, angle, pointcloud, visibility);
	}
	// run all camera-point visibility intersections
	Util::Progress progress(_T("Point visibility checks"), pointcloud.GetSize());
	#ifdef DENSE_USE_OPENMP
	#pragma omp parallel for //schedule(dynamic)
	for (int64_t i=0; i<(int64_t)pointcloud.GetSize(); ++i) {
	const PointCloud::Index idxPoint((PointCloud::Index)i);
	#else
	FOREACH(idxPoint, pointcloud.points) {
	#endif
	const PointCloud::Point& X = pointcloud.points[idxPoint];
	const PointCloud::ViewArr& views = pointcloud.pointViews[idxPoint];
	for (PointCloud::View idxView: views) {
	Collector& collector = collectors[idxView];
	#ifdef DENSE_USE_OPENMP
	Lock l(collector.GetCS());
	#endif
	collector.Init(idxPoint, X, (int)views.size());
	octree.Collect(collector, collector);
	}
	++progress;
	}
	progress.close();
	#if TD_VERBOSE != TD_VERBOSE_OFF
	if (g_nVerbosityLevel > 2) {
	// print visibility stats
	UnsignedArr counts(0, 64);
	for (int views: visibility) {
	if (views > 0)
	continue;
	while (counts.size() <= IDX(-views))
	counts.push_back(0);
	++counts[-views];
	}
	String msg;
	msg.reserve(64*counts.size());
	FOREACH(c, counts)
	if (counts[c])
	msg += String::FormatString("\n\t% 3u - % 9u", c, counts[c]);
	VERBOSE("Visibility lengths (%u points):%s", pointcloud.GetSize(), msg.c_str());
	// save outlier points
	PointCloud pc;
	RFOREACH(idxPoint, pointcloud.points) {
	if (visibility[idxPoint] <= thRemove) {
	pc.points.push_back(pointcloud.points[idxPoint]);
	pc.colors.push_back(pointcloud.colors[idxPoint]);
	}
	}
	pc.Save(MAKE_PATH("scene_dense_outliers.ply"));
	}
	#endif
	// filter points
	const size_t numInitPoints(pointcloud.GetSize());
	RFOREACH(idxPoint, pointcloud.points) {
	if (visibility[idxPoint] <= thRemove)
	pointcloud.RemovePoint(idxPoint);
	}
	DEBUG_EXTRA("Point-cloud filtered: %u/%u points (%d%%%%) (%s)", pointcloud.points.size(), numInitPoints, ROUND2INT((100.f*pointcloud.points.GetSize())/numInitPoints), TD_TIMER_GET_FMT().c_str());
	} // PointCloudFilter