Merge pull request #2566 from akashsharma02:master

[GSoC] Add new Hashtable based TSDF volume to RGBD module

* - Add HashTSDF class
    - Implement Integrate function (untested)

* Integration seems to be working, raycasting does not

* Update integration code

* Integration and Raycasting fixes, (both work now)

* - Format code
- Clean up comments and few fixes

* Add Kinect Fusion backup file

* - Add interpolation for vertices and normals (slow and unreliable!)
- Format code
- Delete kinfu_back.cpp

* Bug fix for integration and noisy odometry

* - Create volume abstract class
- Address Review comments

* - Add getPoints and getNormals function
- Fix formatting according to comments
- Move volume abstract class to include/opencv2/rgbd/
- Write factory method for creating TSDFVolumes
- Small bug fixes
- Minor fixes according to comments

* - Add tests for hashTSDF
- Fix raycasting bug causing to loop forever
- Suppress warnings by explicit conversion
- Disable hashTsdf test until we figure out memory leak
- style changes
- Add missing license in a few files, correct precomp.hpp usage
pull/2595/head
Akash Sharma 5 years ago committed by GitHub
parent 52200a82e7
commit 468345511f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
  1. 81
      modules/rgbd/include/opencv2/rgbd/intrinsics.hpp
  2. 8
      modules/rgbd/include/opencv2/rgbd/kinfu.hpp
  3. 58
      modules/rgbd/include/opencv2/rgbd/volume.hpp
  4. 9
      modules/rgbd/samples/kinfu_demo.cpp
  5. 578
      modules/rgbd/src/hash_tsdf.cpp
  6. 112
      modules/rgbd/src/hash_tsdf.hpp
  7. 96
      modules/rgbd/src/kinfu.cpp
  8. 14
      modules/rgbd/src/kinfu_frame.cpp
  9. 1
      modules/rgbd/src/kinfu_frame.hpp
  10. 467
      modules/rgbd/src/tsdf.cpp
  11. 111
      modules/rgbd/src/tsdf.hpp
  12. 2
      modules/rgbd/src/utils.cpp
  13. 62
      modules/rgbd/src/utils.hpp
  14. 34
      modules/rgbd/src/volume.cpp
  15. 11
      modules/rgbd/test/test_kinfu.cpp

@ -0,0 +1,81 @@
// This file is part of OpenCV project.
// It is subject to the license terms in the LICENSE file found in the top-level directory
// of this distribution and at http://opencv.org/license.html
#ifndef __OPENCV_RGBD_INTRINSICS_HPP__
#define __OPENCV_RGBD_INTRINSICS_HPP__
#include "opencv2/core/matx.hpp"
namespace cv
{
namespace kinfu
{
struct Intr
{
/** @brief Camera intrinsics */
/** Reprojects screen point to camera space given z coord. */
struct Reprojector
{
Reprojector() {}
inline Reprojector(Intr intr)
{
fxinv = 1.f/intr.fx, fyinv = 1.f/intr.fy;
cx = intr.cx, cy = intr.cy;
}
template<typename T>
inline cv::Point3_<T> operator()(cv::Point3_<T> p) const
{
T x = p.z * (p.x - cx) * fxinv;
T y = p.z * (p.y - cy) * fyinv;
return cv::Point3_<T>(x, y, p.z);
}
float fxinv, fyinv, cx, cy;
};
/** Projects camera space vector onto screen */
struct Projector
{
inline Projector(Intr intr) : fx(intr.fx), fy(intr.fy), cx(intr.cx), cy(intr.cy) { }
template<typename T>
inline cv::Point_<T> operator()(cv::Point3_<T> p) const
{
T invz = T(1)/p.z;
T x = fx*(p.x*invz) + cx;
T y = fy*(p.y*invz) + cy;
return cv::Point_<T>(x, y);
}
template<typename T>
inline cv::Point_<T> operator()(cv::Point3_<T> p, cv::Point3_<T>& pixVec) const
{
T invz = T(1)/p.z;
pixVec = cv::Point3_<T>(p.x*invz, p.y*invz, 1);
T x = fx*pixVec.x + cx;
T y = fy*pixVec.y + cy;
return cv::Point_<T>(x, y);
}
float fx, fy, cx, cy;
};
Intr() : fx(), fy(), cx(), cy() { }
Intr(float _fx, float _fy, float _cx, float _cy) : fx(_fx), fy(_fy), cx(_cx), cy(_cy) { }
Intr(cv::Matx33f m) : fx(m(0, 0)), fy(m(1, 1)), cx(m(0, 2)), cy(m(1, 2)) { }
// scale intrinsics to pyramid level
inline Intr scale(int pyr) const
{
float factor = (1.f /(1 << pyr));
return Intr(fx*factor, fy*factor, cx*factor, cy*factor);
}
inline Reprojector makeReprojector() const { return Reprojector(*this); }
inline Projector makeProjector() const { return Projector(*this); }
inline cv::Matx33f getMat() const { return Matx33f(fx, 0, cx, 0, fy, cy, 0, 0, 1); }
float fx, fy, cx, cy;
};
} // namespace rgbd
} // namespace cv
#endif

@ -9,6 +9,7 @@
#include "opencv2/core.hpp"
#include "opencv2/core/affine.hpp"
#include <opencv2/rgbd/volume.hpp>
namespace cv {
namespace kinfu {
@ -68,9 +69,16 @@ struct CV_EXPORTS_W Params
*/
CV_WRAP static Ptr<Params> coarseParams();
/** @brief HashTSDF parameters
A set of parameters suitable for use with HashTSDFVolume
*/
CV_WRAP static Ptr<Params> hashTSDFParams(bool isCoarse);
/** @brief frame size in pixels */
CV_PROP_RW Size frameSize;
CV_PROP_RW cv::kinfu::VolumeType volumeType;
/** @brief camera intrinsics */
CV_PROP_RW Matx33f intr;

@ -0,0 +1,58 @@
// This file is part of OpenCV project.
// It is subject to the license terms in the LICENSE file found in the top-level directory
// of this distribution and at http://opencv.org/license.html
// This code is also subject to the license terms in the LICENSE_KinectFusion.md file found in this
// module's directory
#ifndef __OPENCV_RGBD_VOLUME_H__
#define __OPENCV_RGBD_VOLUME_H__
#include "intrinsics.hpp"
#include "opencv2/core/affine.hpp"
namespace cv
{
namespace kinfu
{
class Volume
{
public:
Volume(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor)
: voxelSize(_voxelSize),
voxelSizeInv(1.0f / voxelSize),
pose(_pose),
raycastStepFactor(_raycastStepFactor)
{
}
virtual ~Volume(){};
virtual void integrate(InputArray _depth, float depthFactor, const cv::Affine3f& cameraPose,
const cv::kinfu::Intr& intrinsics) = 0;
virtual void raycast(const cv::Affine3f& cameraPose, const cv::kinfu::Intr& intrinsics,
cv::Size frameSize, cv::OutputArray points,
cv::OutputArray normals) const = 0;
virtual void fetchNormals(cv::InputArray points, cv::OutputArray _normals) const = 0;
virtual void fetchPointsNormals(cv::OutputArray points, cv::OutputArray normals) const = 0;
virtual void reset() = 0;
public:
const float voxelSize;
const float voxelSizeInv;
const cv::Affine3f pose;
const float raycastStepFactor;
};
enum class VolumeType
{
TSDF = 0,
HASHTSDF = 1
};
cv::Ptr<Volume> makeVolume(VolumeType _volumeType, float _voxelSize, cv::Affine3f _pose,
float _raycastStepFactor, float _truncDist, int _maxWeight,
float _truncateThreshold, Point3i _resolution);
} // namespace kinfu
} // namespace cv
#endif

@ -312,6 +312,7 @@ static const char* keys =
"{camera |0| Index of depth camera to be used as a depth source }"
"{coarse | | Run on coarse settings (fast but ugly) or on default (slow but looks better),"
" in coarse mode points and normals are displayed }"
"{useHashTSDF | | Use the newer hashtable based TSDFVolume (relatively fast) and for larger reconstructions}"
"{idle | | Do not run KinFu, just display depth frames }"
"{record | | Write depth frames to specified file list"
" (the same format as for the 'depth' key) }"
@ -327,6 +328,7 @@ int main(int argc, char **argv)
{
bool coarse = false;
bool idle = false;
bool useHashTSDF = false;
string recordPath;
CommandLineParser parser(argc, argv, keys);
@ -352,6 +354,10 @@ int main(int argc, char **argv)
{
recordPath = parser.get<String>("record");
}
if(parser.has("useHashTSDF"))
{
useHashTSDF = true;
}
if(parser.has("idle"))
{
idle = true;
@ -382,6 +388,9 @@ int main(int argc, char **argv)
else
params = Params::defaultParams();
if(useHashTSDF)
params = Params::hashTSDFParams(coarse);
// These params can be different for each depth sensor
ds->updateParams(*params);

@ -0,0 +1,578 @@
// This file is part of OpenCV project.
// It is subject to the license terms in the LICENSE file found in the top-level directory
// of this distribution and at http://opencv.org/license.html
#include "precomp.hpp"
#include "hash_tsdf.hpp"
#include <atomic>
#include <functional>
#include <iostream>
#include <limits>
#include <vector>
#include "kinfu_frame.hpp"
#include "opencv2/core/cvstd.hpp"
#include "opencv2/core/utility.hpp"
#include "opencv2/core/utils/trace.hpp"
#include "utils.hpp"
namespace cv
{
namespace kinfu
{
HashTSDFVolume::HashTSDFVolume(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor,
float _truncDist, int _maxWeight, float _truncateThreshold,
int _volumeUnitRes, bool _zFirstMemOrder)
: Volume(_voxelSize, _pose, _raycastStepFactor),
maxWeight(_maxWeight),
truncateThreshold(_truncateThreshold),
volumeUnitResolution(_volumeUnitRes),
volumeUnitSize(voxelSize * volumeUnitResolution),
zFirstMemOrder(_zFirstMemOrder)
{
truncDist = std::max(_truncDist, 4.0f * voxelSize);
}
HashTSDFVolumeCPU::HashTSDFVolumeCPU(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor,
float _truncDist, int _maxWeight, float _truncateThreshold,
int _volumeUnitRes, bool _zFirstMemOrder)
: HashTSDFVolume(_voxelSize, _pose, _raycastStepFactor, _truncDist, _maxWeight,
_truncateThreshold, _volumeUnitRes, _zFirstMemOrder)
{
}
// zero volume, leave rest params the same
void HashTSDFVolumeCPU::reset()
{
CV_TRACE_FUNCTION();
volumeUnits.clear();
}
struct AllocateVolumeUnitsInvoker : ParallelLoopBody
{
AllocateVolumeUnitsInvoker(HashTSDFVolumeCPU& _volume, const Depth& _depth, Intr intrinsics,
cv::Affine3f cameraPose, float _depthFactor, int _depthStride = 4)
: ParallelLoopBody(),
volume(_volume),
depth(_depth),
reproj(intrinsics.makeReprojector()),
cam2vol(_volume.pose.inv() * cameraPose),
depthFactor(1.0f / _depthFactor),
depthStride(_depthStride)
{
}
virtual void operator()(const Range& range) const override
{
for (int y = range.start; y < range.end; y += depthStride)
{
const depthType* depthRow = depth[y];
for (int x = 0; x < depth.cols; x += depthStride)
{
depthType z = depthRow[x] * depthFactor;
if (z <= 0)
continue;
Point3f camPoint = reproj(Point3f((float)x, (float)y, z));
Point3f volPoint = cam2vol * camPoint;
//! Find accessed TSDF volume unit for valid 3D vertex
cv::Vec3i lower_bound = volume.volumeToVolumeUnitIdx(
volPoint - cv::Point3f(volume.truncDist, volume.truncDist, volume.truncDist));
cv::Vec3i upper_bound = volume.volumeToVolumeUnitIdx(
volPoint + cv::Point3f(volume.truncDist, volume.truncDist, volume.truncDist));
VolumeUnitIndexSet localAccessVolUnits;
for (int i = lower_bound[0]; i <= upper_bound[0]; i++)
for (int j = lower_bound[1]; j <= upper_bound[1]; j++)
for (int k = lower_bound[2]; k <= lower_bound[2]; k++)
{
const cv::Vec3i tsdf_idx = cv::Vec3i(i, j, k);
if (!localAccessVolUnits.count(tsdf_idx))
{
localAccessVolUnits.emplace(tsdf_idx);
}
}
AutoLock al(mutex);
for (const auto& tsdf_idx : localAccessVolUnits)
{
//! If the insert into the global set passes
if (!volume.volumeUnits.count(tsdf_idx))
{
VolumeUnit volumeUnit;
cv::Point3i volumeDims(volume.volumeUnitResolution,
volume.volumeUnitResolution,
volume.volumeUnitResolution);
cv::Affine3f subvolumePose =
volume.pose.translate(volume.volumeUnitIdxToVolume(tsdf_idx));
volumeUnit.pVolume = cv::makePtr<TSDFVolumeCPU>(
volume.voxelSize, subvolumePose, volume.raycastStepFactor,
volume.truncDist, volume.maxWeight, volumeDims);
//! This volume unit will definitely be required for current integration
volumeUnit.index = tsdf_idx;
volumeUnit.isActive = true;
volume.volumeUnits[tsdf_idx] = volumeUnit;
}
}
}
}
}
HashTSDFVolumeCPU& volume;
const Depth& depth;
const Intr::Reprojector reproj;
const cv::Affine3f cam2vol;
const float depthFactor;
const int depthStride;
mutable Mutex mutex;
};
void HashTSDFVolumeCPU::integrate(InputArray _depth, float depthFactor,
const cv::Affine3f& cameraPose, const Intr& intrinsics)
{
CV_TRACE_FUNCTION();
CV_Assert(_depth.type() == DEPTH_TYPE);
Depth depth = _depth.getMat();
//! Compute volumes to be allocated
AllocateVolumeUnitsInvoker allocate_i(*this, depth, intrinsics, cameraPose, depthFactor);
Range allocateRange(0, depth.rows);
parallel_for_(allocateRange, allocate_i);
//! Get volumes that are in the current camera frame
std::vector<Vec3i> totalVolUnits;
for (const auto& keyvalue : volumeUnits)
{
totalVolUnits.push_back(keyvalue.first);
}
//! Mark volumes in the camera frustum as active
Range inFrustumRange(0, (int)volumeUnits.size());
parallel_for_(inFrustumRange, [&](const Range& range) {
const Affine3f vol2cam(cameraPose.inv() * pose);
const Intr::Projector proj(intrinsics.makeProjector());
for (int i = range.start; i < range.end; ++i)
{
cv::Vec3i tsdf_idx = totalVolUnits[i];
VolumeUnitMap::iterator it = volumeUnits.find(tsdf_idx);
if (it == volumeUnits.end())
return;
Point3f volumeUnitPos = volumeUnitIdxToVolume(it->first);
Point3f volUnitInCamSpace = vol2cam * volumeUnitPos;
if (volUnitInCamSpace.z < 0 || volUnitInCamSpace.z > truncateThreshold)
{
it->second.isActive = false;
return;
}
Point2f cameraPoint = proj(volUnitInCamSpace);
if (cameraPoint.x >= 0 && cameraPoint.y >= 0 && cameraPoint.x < depth.cols &&
cameraPoint.y < depth.rows)
{
assert(it != volumeUnits.end());
it->second.isActive = true;
}
}
});
//! Integrate the correct volumeUnits
parallel_for_(Range(0, (int)totalVolUnits.size()), [&](const Range& range) {
for (int i = range.start; i < range.end; i++)
{
cv::Vec3i tsdf_idx = totalVolUnits[i];
VolumeUnitMap::iterator it = volumeUnits.find(tsdf_idx);
if (it == volumeUnits.end())
return;
VolumeUnit& volumeUnit = it->second;
if (volumeUnit.isActive)
{
//! The volume unit should already be added into the Volume from the allocator
volumeUnit.pVolume->integrate(depth, depthFactor, cameraPose, intrinsics);
//! Ensure all active volumeUnits are set to inactive for next integration
volumeUnit.isActive = false;
}
}
});
}
cv::Vec3i HashTSDFVolumeCPU::volumeToVolumeUnitIdx(cv::Point3f p) const
{
return cv::Vec3i(cvFloor(p.x / volumeUnitSize), cvFloor(p.y / volumeUnitSize),
cvFloor(p.z / volumeUnitSize));
}
cv::Point3f HashTSDFVolumeCPU::volumeUnitIdxToVolume(cv::Vec3i volumeUnitIdx) const
{
return cv::Point3f(volumeUnitIdx[0] * volumeUnitSize, volumeUnitIdx[1] * volumeUnitSize,
volumeUnitIdx[2] * volumeUnitSize);
}
cv::Point3f HashTSDFVolumeCPU::voxelCoordToVolume(cv::Vec3i voxelIdx) const
{
return cv::Point3f(voxelIdx[0] * voxelSize, voxelIdx[1] * voxelSize, voxelIdx[2] * voxelSize);
}
cv::Vec3i HashTSDFVolumeCPU::volumeToVoxelCoord(cv::Point3f point) const
{
return cv::Vec3i(cvFloor(point.x * voxelSizeInv), cvFloor(point.y * voxelSizeInv),
cvFloor(point.z * voxelSizeInv));
}
inline TsdfVoxel HashTSDFVolumeCPU::at(const cv::Vec3i& volumeIdx) const
{
cv::Vec3i volumeUnitIdx = cv::Vec3i(cvFloor(volumeIdx[0] / volumeUnitResolution),
cvFloor(volumeIdx[1] / volumeUnitResolution),
cvFloor(volumeIdx[2] / volumeUnitResolution));
VolumeUnitMap::const_iterator it = volumeUnits.find(volumeUnitIdx);
if (it == volumeUnits.end())
{
TsdfVoxel dummy;
dummy.tsdf = 1.f;
dummy.weight = 0;
return dummy;
}
cv::Ptr<TSDFVolumeCPU> volumeUnit =
std::dynamic_pointer_cast<TSDFVolumeCPU>(it->second.pVolume);
cv::Vec3i volUnitLocalIdx = volumeIdx - cv::Vec3i(volumeUnitIdx[0] * volumeUnitResolution,
volumeUnitIdx[1] * volumeUnitResolution,
volumeUnitIdx[2] * volumeUnitResolution);
volUnitLocalIdx =
cv::Vec3i(abs(volUnitLocalIdx[0]), abs(volUnitLocalIdx[1]), abs(volUnitLocalIdx[2]));
return volumeUnit->at(volUnitLocalIdx);
}
inline TsdfVoxel HashTSDFVolumeCPU::at(const cv::Point3f& point) const
{
cv::Vec3i volumeUnitIdx = volumeToVolumeUnitIdx(point);
VolumeUnitMap::const_iterator it = volumeUnits.find(volumeUnitIdx);
if (it == volumeUnits.end())
{
TsdfVoxel dummy;
dummy.tsdf = 1.f;
dummy.weight = 0;
return dummy;
}
cv::Ptr<TSDFVolumeCPU> volumeUnit =
std::dynamic_pointer_cast<TSDFVolumeCPU>(it->second.pVolume);
cv::Point3f volumeUnitPos = volumeUnitIdxToVolume(volumeUnitIdx);
cv::Vec3i volUnitLocalIdx = volumeToVoxelCoord(point - volumeUnitPos);
volUnitLocalIdx =
cv::Vec3i(abs(volUnitLocalIdx[0]), abs(volUnitLocalIdx[1]), abs(volUnitLocalIdx[2]));
return volumeUnit->at(volUnitLocalIdx);
}
inline Point3f HashTSDFVolumeCPU::getNormalVoxel(Point3f point) const
{
Vec3f pointVec(point);
Vec3f normal = Vec3f(0, 0, 0);
Vec3f pointPrev = point;
Vec3f pointNext = point;
for (int c = 0; c < 3; c++)
{
pointPrev[c] -= voxelSize * 0.5f;
pointNext[c] += voxelSize * 0.5f;
normal[c] = at(Point3f(pointNext)).tsdf - at(Point3f(pointPrev)).tsdf;
normal[c] *= 0.5f;
pointPrev[c] = pointVec[c];
pointNext[c] = pointVec[c];
}
return normalize(normal);
}
struct HashRaycastInvoker : ParallelLoopBody
{
HashRaycastInvoker(Points& _points, Normals& _normals, const Affine3f& cameraPose,
const Intr& intrinsics, const HashTSDFVolumeCPU& _volume)
: ParallelLoopBody(),
points(_points),
normals(_normals),
volume(_volume),
tstep(_volume.truncDist * _volume.raycastStepFactor),
cam2vol(volume.pose.inv() * cameraPose),
vol2cam(cameraPose.inv() * volume.pose),
reproj(intrinsics.makeReprojector())
{
}
virtual void operator()(const Range& range) const override
{
const Point3f cam2volTrans = cam2vol.translation();
const Matx33f cam2volRot = cam2vol.rotation();
const Matx33f vol2camRot = vol2cam.rotation();
const float blockSize = volume.volumeUnitSize;
for (int y = range.start; y < range.end; y++)
{
ptype* ptsRow = points[y];
ptype* nrmRow = normals[y];
for (int x = 0; x < points.cols; x++)
{
//! Initialize default value
Point3f point = nan3, normal = nan3;
//! Ray origin and direction in the volume coordinate frame
Point3f orig = cam2volTrans;
Point3f rayDirV =
normalize(Vec3f(cam2volRot * reproj(Point3f(float(x), float(y), 1.f))));
float tmin = 0;
float tmax = volume.truncateThreshold;
float tcurr = tmin;
cv::Vec3i prevVolumeUnitIdx =
cv::Vec3i(std::numeric_limits<int>::min(), std::numeric_limits<int>::min(),
std::numeric_limits<int>::min());
float tprev = tcurr;
TsdfType prevTsdf = volume.truncDist;
cv::Ptr<TSDFVolumeCPU> currVolumeUnit;
while (tcurr < tmax)
{
Point3f currRayPos = orig + tcurr * rayDirV;
cv::Vec3i currVolumeUnitIdx = volume.volumeToVolumeUnitIdx(currRayPos);
VolumeUnitMap::const_iterator it = volume.volumeUnits.find(currVolumeUnitIdx);
TsdfType currTsdf = prevTsdf;
int currWeight = 0;
float stepSize = 0.5f * blockSize;
cv::Vec3i volUnitLocalIdx;
//! Does the subvolume exist in hashtable
if (it != volume.volumeUnits.end())
{
currVolumeUnit =
std::dynamic_pointer_cast<TSDFVolumeCPU>(it->second.pVolume);
cv::Point3f currVolUnitPos =
volume.volumeUnitIdxToVolume(currVolumeUnitIdx);
volUnitLocalIdx = volume.volumeToVoxelCoord(currRayPos - currVolUnitPos);
//! TODO: Figure out voxel interpolation
TsdfVoxel currVoxel = currVolumeUnit->at(volUnitLocalIdx);
currTsdf = currVoxel.tsdf;
currWeight = currVoxel.weight;
stepSize = tstep;
}
//! Surface crossing
if (prevTsdf > 0.f && currTsdf <= 0.f && currWeight > 0)
{
float tInterp =
(tcurr * prevTsdf - tprev * currTsdf) / (prevTsdf - currTsdf);
if (!cvIsNaN(tInterp) && !cvIsInf(tInterp))
{
Point3f pv = orig + tInterp * rayDirV;
Point3f nv = volume.getNormalVoxel(pv);
if (!isNaN(nv))
{
normal = vol2camRot * nv;
point = vol2cam * pv;
}
}
break;
}
prevVolumeUnitIdx = currVolumeUnitIdx;
prevTsdf = currTsdf;
tprev = tcurr;
tcurr += stepSize;
}
ptsRow[x] = toPtype(point);
nrmRow[x] = toPtype(normal);
}
}
}
Points& points;
Normals& normals;
const HashTSDFVolumeCPU& volume;
const float tstep;
const Affine3f cam2vol;
const Affine3f vol2cam;
const Intr::Reprojector reproj;
};
void HashTSDFVolumeCPU::raycast(const cv::Affine3f& cameraPose, const cv::kinfu::Intr& intrinsics,
cv::Size frameSize, cv::OutputArray _points,
cv::OutputArray _normals) const
{
CV_TRACE_FUNCTION();
CV_Assert(frameSize.area() > 0);
_points.create(frameSize, POINT_TYPE);
_normals.create(frameSize, POINT_TYPE);
Points points = _points.getMat();
Normals normals = _normals.getMat();
HashRaycastInvoker ri(points, normals, cameraPose, intrinsics, *this);
const int nstripes = -1;
parallel_for_(Range(0, points.rows), ri, nstripes);
}
struct FetchPointsNormalsInvoker : ParallelLoopBody
{
FetchPointsNormalsInvoker(const HashTSDFVolumeCPU& _volume,
const std::vector<Vec3i>& _totalVolUnits,
std::vector<std::vector<ptype>>& _pVecs,
std::vector<std::vector<ptype>>& _nVecs, bool _needNormals)
: ParallelLoopBody(),
volume(_volume),
totalVolUnits(_totalVolUnits),
pVecs(_pVecs),
nVecs(_nVecs),
needNormals(_needNormals)
{
}
virtual void operator()(const Range& range) const override
{
std::vector<ptype> points, normals;
for (int i = range.start; i < range.end; i++)
{
cv::Vec3i tsdf_idx = totalVolUnits[i];
VolumeUnitMap::const_iterator it = volume.volumeUnits.find(tsdf_idx);
Point3f base_point = volume.volumeUnitIdxToVolume(tsdf_idx);
if (it != volume.volumeUnits.end())
{
cv::Ptr<TSDFVolumeCPU> volumeUnit =
std::dynamic_pointer_cast<TSDFVolumeCPU>(it->second.pVolume);
std::vector<ptype> localPoints;
std::vector<ptype> localNormals;
for (int x = 0; x < volume.volumeUnitResolution; x++)
for (int y = 0; y < volume.volumeUnitResolution; y++)
for (int z = 0; z < volume.volumeUnitResolution; z++)
{
cv::Vec3i voxelIdx(x, y, z);
TsdfVoxel voxel = volumeUnit->at(voxelIdx);
if (voxel.tsdf != 1.f && voxel.weight != 0)
{
Point3f point = base_point + volume.voxelCoordToVolume(voxelIdx);
localPoints.push_back(toPtype(point));
if (needNormals)
{
Point3f normal = volume.getNormalVoxel(point);
localNormals.push_back(toPtype(normal));
}
}
}
AutoLock al(mutex);
pVecs.push_back(localPoints);
nVecs.push_back(localNormals);
}
}
}
const HashTSDFVolumeCPU& volume;
std::vector<cv::Vec3i> totalVolUnits;
std::vector<std::vector<ptype>>& pVecs;
std::vector<std::vector<ptype>>& nVecs;
const TsdfVoxel* volDataStart;
bool needNormals;
mutable Mutex mutex;
};
void HashTSDFVolumeCPU::fetchPointsNormals(OutputArray _points, OutputArray _normals) const
{
CV_TRACE_FUNCTION();
if (_points.needed())
{
std::vector<std::vector<ptype>> pVecs, nVecs;
std::vector<Vec3i> totalVolUnits;
for (const auto& keyvalue : volumeUnits)
{
totalVolUnits.push_back(keyvalue.first);
}
FetchPointsNormalsInvoker fi(*this, totalVolUnits, pVecs, nVecs, _normals.needed());
Range range(0, (int)totalVolUnits.size());
const int nstripes = -1;
parallel_for_(range, fi, nstripes);
std::vector<ptype> points, normals;
for (size_t i = 0; i < pVecs.size(); i++)
{
points.insert(points.end(), pVecs[i].begin(), pVecs[i].end());
normals.insert(normals.end(), nVecs[i].begin(), nVecs[i].end());
}
_points.create((int)points.size(), 1, POINT_TYPE);
if (!points.empty())
Mat((int)points.size(), 1, POINT_TYPE, &points[0]).copyTo(_points.getMat());
if (_normals.needed())
{
_normals.create((int)normals.size(), 1, POINT_TYPE);
if (!normals.empty())
Mat((int)normals.size(), 1, POINT_TYPE, &normals[0]).copyTo(_normals.getMat());
}
}
}
struct PushNormals
{
PushNormals(const HashTSDFVolumeCPU& _volume, Normals& _normals)
: volume(_volume), normals(_normals), invPose(volume.pose.inv())
{
}
void operator()(const ptype& point, const int* position) const
{
Point3f p = fromPtype(point);
Point3f n = nan3;
if (!isNaN(p))
{
Point3f voxelPoint = invPose * p;
n = volume.pose.rotation() * volume.getNormalVoxel(voxelPoint);
}
normals(position[0], position[1]) = toPtype(n);
}
const HashTSDFVolumeCPU& volume;
Normals& normals;
Affine3f invPose;
};
void HashTSDFVolumeCPU::fetchNormals(cv::InputArray _points, cv::OutputArray _normals) const
{
CV_TRACE_FUNCTION();
if (_normals.needed())
{
Points points = _points.getMat();
CV_Assert(points.type() == POINT_TYPE);
_normals.createSameSize(_points, _points.type());
Normals normals = _normals.getMat();
points.forEach(PushNormals(*this, normals));
}
}
cv::Ptr<HashTSDFVolume> makeHashTSDFVolume(float _voxelSize, cv::Affine3f _pose,
float _raycastStepFactor, float _truncDist,
int _maxWeight, float _truncateThreshold,
int _volumeUnitResolution)
{
return cv::makePtr<HashTSDFVolumeCPU>(_voxelSize, _pose, _raycastStepFactor, _truncDist,
_maxWeight, _truncateThreshold, _volumeUnitResolution);
}
} // namespace kinfu
} // namespace cv

@ -0,0 +1,112 @@
// This file is part of OpenCV project.
// It is subject to the license terms in the LICENSE file found in the top-level directory
// of this distribution and at http://opencv.org/license.html
#ifndef __OPENCV_HASH_TSDF_H__
#define __OPENCV_HASH_TSDF_H__
#include <opencv2/rgbd/volume.hpp>
#include <unordered_map>
#include <unordered_set>
#include "tsdf.hpp"
namespace cv
{
namespace kinfu
{
class HashTSDFVolume : public Volume
{
public:
// dimension in voxels, size in meters
//! Use fixed volume cuboid
HashTSDFVolume(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor, float _truncDist,
int _maxWeight, float _truncateThreshold, int _volumeUnitRes,
bool zFirstMemOrder = true);
virtual ~HashTSDFVolume() = default;
public:
int maxWeight;
float truncDist;
float truncateThreshold;
int volumeUnitResolution;
float volumeUnitSize;
bool zFirstMemOrder;
};
struct VolumeUnit
{
VolumeUnit() : pVolume(nullptr){};
~VolumeUnit() = default;
cv::Ptr<TSDFVolume> pVolume;
cv::Vec3i index;
bool isActive;
};
//! Spatial hashing
struct tsdf_hash
{
size_t operator()(const cv::Vec3i& x) const noexcept
{
size_t seed = 0;
constexpr uint32_t GOLDEN_RATIO = 0x9e3779b9;
for (uint16_t i = 0; i < 3; i++)
{
seed ^= std::hash<int>()(x[i]) + GOLDEN_RATIO + (seed << 6) + (seed >> 2);
}
return seed;
}
};
typedef std::unordered_set<cv::Vec3i, tsdf_hash> VolumeUnitIndexSet;
typedef std::unordered_map<cv::Vec3i, VolumeUnit, tsdf_hash> VolumeUnitMap;
class HashTSDFVolumeCPU : public HashTSDFVolume
{
public:
// dimension in voxels, size in meters
HashTSDFVolumeCPU(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor,
float _truncDist, int _maxWeight, float _truncateThreshold,
int _volumeUnitRes, bool zFirstMemOrder = true);
virtual void integrate(InputArray _depth, float depthFactor, const cv::Affine3f& cameraPose,
const cv::kinfu::Intr& intrinsics) override;
virtual void raycast(const cv::Affine3f& cameraPose, const cv::kinfu::Intr& intrinsics,
cv::Size frameSize, cv::OutputArray points,
cv::OutputArray normals) const override;
virtual void fetchNormals(cv::InputArray points, cv::OutputArray _normals) const override;
virtual void fetchPointsNormals(cv::OutputArray points, cv::OutputArray normals) const override;
virtual void reset() override;
//! Return the voxel given the voxel index in the universal volume (1 unit = 1 voxel_length)
virtual TsdfVoxel at(const cv::Vec3i& volumeIdx) const;
//! Return the voxel given the point in volume coordinate system i.e., (metric scale 1 unit =
//! 1m)
virtual TsdfVoxel at(const cv::Point3f& point) const;
inline TsdfType interpolateVoxel(const cv::Point3f& point) const;
Point3f getNormalVoxel(cv::Point3f p) const;
//! Utility functions for coordinate transformations
cv::Vec3i volumeToVolumeUnitIdx(cv::Point3f point) const;
cv::Point3f volumeUnitIdxToVolume(cv::Vec3i volumeUnitIdx) const;
cv::Point3f voxelCoordToVolume(cv::Vec3i voxelIdx) const;
cv::Vec3i volumeToVoxelCoord(cv::Point3f point) const;
public:
//! Hashtable of individual smaller volume units
VolumeUnitMap volumeUnits;
};
cv::Ptr<HashTSDFVolume> makeHashTSDFVolume(float _voxelSize, cv::Affine3f _pose,
float _raycastStepFactor, float _truncDist,
int _maxWeight, float truncateThreshold,
int volumeUnitResolution = 16);
} // namespace kinfu
} // namespace cv
#endif

@ -7,6 +7,7 @@
#include "precomp.hpp"
#include "fast_icp.hpp"
#include "tsdf.hpp"
#include "hash_tsdf.hpp"
#include "kinfu_frame.hpp"
namespace cv {
@ -28,6 +29,8 @@ Ptr<Params> Params::defaultParams()
p.frameSize = Size(640, 480);
p.volumeType = VolumeType::TSDF;
float fx, fy, cx, cy;
fx = fy = 525.f;
cx = p.frameSize.width/2 - 0.5f;
@ -60,7 +63,7 @@ Ptr<Params> Params::defaultParams()
// default pose of volume cube
p.volumePose = Affine3f().translate(Vec3f(-volSize/2.f, -volSize/2.f, 0.5f));
p.tsdf_trunc_dist = 0.04f; //meters;
p.tsdf_trunc_dist = 7 * p.voxelSize; // about 0.04f in meters
p.tsdf_max_weight = 64; //frames
p.raycast_step_factor = 0.25f; //in voxel sizes
@ -86,14 +89,26 @@ Ptr<Params> Params::coarseParams()
float volSize = 3.f;
p->volumeDims = Vec3i::all(128); //number of voxels
p->voxelSize = volSize/128.f;
p->tsdf_trunc_dist = 2 * p->voxelSize; // 0.04f in meters
p->raycast_step_factor = 0.75f; //in voxel sizes
return p;
}
Ptr<Params> Params::hashTSDFParams(bool isCoarse)
{
Ptr<Params> p;
if(isCoarse)
p = coarseParams();
else
p = defaultParams();
p->volumeType = VolumeType::HASHTSDF;
p->truncateThreshold = rgbd::Odometry::DEFAULT_MAX_DEPTH();
return p;
}
// T should be Mat or UMat
template< typename T >
// MatType should be Mat or UMat
template< typename MatType>
class KinFuImpl : public KinFu
{
public:
@ -104,7 +119,8 @@ public:
void render(OutputArray image, const Matx44f& cameraPose) const CV_OVERRIDE;
void getCloud(OutputArray points, OutputArray normals) const CV_OVERRIDE;
//! TODO(Akash): Add back later
virtual void getCloud(OutputArray points, OutputArray normals) const CV_OVERRIDE;
void getPoints(OutputArray points) const CV_OVERRIDE;
void getNormals(InputArray points, OutputArray normals) const CV_OVERRIDE;
@ -114,53 +130,52 @@ public:
bool update(InputArray depth) CV_OVERRIDE;
bool updateT(const T& depth);
bool updateT(const MatType& depth);
private:
Params params;
cv::Ptr<ICP> icp;
cv::Ptr<TSDFVolume> volume;
cv::Ptr<Volume> volume;
int frameCounter;
Affine3f pose;
std::vector<T> pyrPoints;
std::vector<T> pyrNormals;
std::vector<MatType> pyrPoints;
std::vector<MatType> pyrNormals;
};
template< typename T >
KinFuImpl<T>::KinFuImpl(const Params &_params) :
template< typename MatType >
KinFuImpl<MatType>::KinFuImpl(const Params &_params) :
params(_params),
icp(makeICP(params.intr, params.icpIterations, params.icpAngleThresh, params.icpDistThresh)),
volume(makeTSDFVolume(params.volumeDims, params.voxelSize, params.volumePose,
params.tsdf_trunc_dist, params.tsdf_max_weight,
params.raycast_step_factor)),
pyrPoints(), pyrNormals()
{
volume = makeVolume(params.volumeType, params.voxelSize, params.volumePose, params.raycast_step_factor,
params.tsdf_trunc_dist, params.tsdf_max_weight, params.truncateThreshold, params.volumeDims);
reset();
}
template< typename T >
void KinFuImpl<T>::reset()
template< typename MatType >
void KinFuImpl<MatType >::reset()
{
frameCounter = 0;
pose = Affine3f::Identity();
volume->reset();
}
template< typename T >
KinFuImpl<T>::~KinFuImpl()
template< typename MatType >
KinFuImpl<MatType>::~KinFuImpl()
{ }
template< typename T >
const Params& KinFuImpl<T>::getParams() const
template< typename MatType >
const Params& KinFuImpl<MatType>::getParams() const
{
return params;
}
template< typename T >
const Affine3f KinFuImpl<T>::getPose() const
template< typename MatType >
const Affine3f KinFuImpl<MatType>::getPose() const
{
return pose;
}
@ -202,18 +217,19 @@ bool KinFuImpl<UMat>::update(InputArray _depth)
}
template< typename T >
bool KinFuImpl<T>::updateT(const T& _depth)
template< typename MatType >
bool KinFuImpl<MatType>::updateT(const MatType& _depth)
{
CV_TRACE_FUNCTION();
T depth;
MatType depth;
if(_depth.type() != DEPTH_TYPE)
_depth.convertTo(depth, DEPTH_TYPE);
else
depth = _depth;
std::vector<T> newPoints, newNormals;
std::vector<MatType> newPoints, newNormals;
makeFrameFromDepth(depth, newPoints, newNormals, params.intr,
params.pyramidLevels,
params.depthFactor,
@ -221,12 +237,10 @@ bool KinFuImpl<T>::updateT(const T& _depth)
params.bilateral_sigma_spatial,
params.bilateral_kernel_size,
params.truncateThreshold);
if(frameCounter == 0)
{
// use depth instead of distance
volume->integrate(depth, params.depthFactor, pose, params.intr);
pyrPoints = newPoints;
pyrNormals = newNormals;
}
@ -247,11 +261,9 @@ bool KinFuImpl<T>::updateT(const T& _depth)
// use depth instead of distance
volume->integrate(depth, params.depthFactor, pose, params.intr);
}
T& points = pyrPoints [0];
T& normals = pyrNormals[0];
MatType& points = pyrPoints [0];
MatType& normals = pyrNormals[0];
volume->raycast(pose, params.intr, params.frameSize, points, normals);
// build a pyramid of points and normals
buildPyramidPointsNormals(points, normals, pyrPoints, pyrNormals,
params.pyramidLevels);
}
@ -261,8 +273,8 @@ bool KinFuImpl<T>::updateT(const T& _depth)
}
template< typename T >
void KinFuImpl<T>::render(OutputArray image, const Matx44f& _cameraPose) const
template< typename MatType >
void KinFuImpl<MatType>::render(OutputArray image, const Matx44f& _cameraPose) const
{
CV_TRACE_FUNCTION();
@ -276,29 +288,29 @@ void KinFuImpl<T>::render(OutputArray image, const Matx44f& _cameraPose) const
}
else
{
T points, normals;
MatType points, normals;
volume->raycast(cameraPose, params.intr, params.frameSize, points, normals);
renderPointsNormals(points, normals, image, params.lightPose);
}
}
template< typename T >
void KinFuImpl<T>::getCloud(OutputArray p, OutputArray n) const
template< typename MatType >
void KinFuImpl<MatType>::getCloud(OutputArray p, OutputArray n) const
{
volume->fetchPointsNormals(p, n);
}
template< typename T >
void KinFuImpl<T>::getPoints(OutputArray points) const
template< typename MatType >
void KinFuImpl<MatType>::getPoints(OutputArray points) const
{
volume->fetchPointsNormals(points, noArray());
}
template< typename T >
void KinFuImpl<T>::getNormals(InputArray points, OutputArray normals) const
template< typename MatType >
void KinFuImpl<MatType>::getNormals(InputArray points, OutputArray normals) const
{
volume->fetchNormals(points, normals);
}
@ -315,11 +327,11 @@ Ptr<KinFu> KinFu::create(const Ptr<Params>& params)
if(cv::ocl::useOpenCL())
return makePtr< KinFuImpl<UMat> >(*params);
#endif
return makePtr< KinFuImpl<Mat> >(*params);
return makePtr< KinFuImpl<Mat> >(*params);
}
#else
Ptr<KinFu> KinFu::create(const Ptr<Params>& /*params*/)
Ptr<KinFu> KinFu::create(const Ptr<Params>& /* params */)
{
CV_Error(Error::StsNotImplemented,
"This algorithm is patented and is excluded in this configuration; "

@ -505,10 +505,13 @@ static bool ocl_makeFrameFromDepth(const UMat depth, OutputArrayOfArrays points,
return false;
// depth truncation can be used in some scenes
UMat depthThreshold;
if(truncateThreshold > 0.f)
threshold(depth, depth, truncateThreshold*depthFactor, 0.0, THRESH_TOZERO_INV);
threshold(smooth, depthThreshold, truncateThreshold * depthFactor, 0.0, THRESH_TOZERO_INV);
else
depthThreshold = smooth;
UMat scaled = smooth;
UMat scaled = depthThreshold;
Size sz = smooth.size();
points.create(levels, 1, POINT_TYPE);
normals.create(levels, 1, POINT_TYPE);
@ -624,13 +627,16 @@ void makeFrameFromDepth(InputArray _depth,
bilateralFilter(depth, smooth, kernelSize, sigmaDepth*depthFactor, sigmaSpatial);
// depth truncation can be used in some scenes
Depth depthThreshold;
if(truncateThreshold > 0.f)
threshold(depth, depth, truncateThreshold, 0.0, THRESH_TOZERO_INV);
threshold(smooth, depthThreshold, truncateThreshold * depthFactor, 0.0, THRESH_TOZERO_INV);
else
depthThreshold = smooth;
// we don't need depth pyramid outside this method
// if we do, the code is to be refactored
Depth scaled = smooth;
Depth scaled = depthThreshold;
Size sz = smooth.size();
pyrPoints.create(levels, 1, POINT_TYPE);
pyrNormals.create(levels, 1, POINT_TYPE);

@ -7,6 +7,7 @@
#ifndef __OPENCV_KINFU_FRAME_H__
#define __OPENCV_KINFU_FRAME_H__
#include <opencv2/core/affine.hpp>
#include "utils.hpp"
namespace cv {

@ -12,63 +12,15 @@ namespace cv {
namespace kinfu {
// TODO: Optimization possible:
// * volumeType can be FP16
// * weight can be int16
typedef float volumeType;
struct Voxel
{
volumeType v;
int weight;
};
typedef Vec<uchar, sizeof(Voxel)> VecT;
class TSDFVolumeCPU : public TSDFVolume
{
public:
// dimension in voxels, size in meters
TSDFVolumeCPU(Point3i _res, float _voxelSize, cv::Affine3f _pose, float _truncDist, int _maxWeight,
float _raycastStepFactor, bool zFirstMemOrder = true);
virtual void integrate(InputArray _depth, float depthFactor, cv::Affine3f cameraPose, cv::kinfu::Intr intrinsics) override;
virtual void raycast(cv::Affine3f cameraPose, cv::kinfu::Intr intrinsics, cv::Size frameSize,
cv::OutputArray points, cv::OutputArray normals) const override;
virtual void fetchNormals(cv::InputArray points, cv::OutputArray _normals) const override;
virtual void fetchPointsNormals(cv::OutputArray points, cv::OutputArray normals) const override;
virtual void reset() override;
volumeType interpolateVoxel(cv::Point3f p) const;
Point3f getNormalVoxel(cv::Point3f p) const;
#if USE_INTRINSICS
volumeType interpolateVoxel(const v_float32x4& p) const;
v_float32x4 getNormalVoxel(const v_float32x4& p) const;
#endif
// See zFirstMemOrder arg of parent class constructor
// for the array layout info
// Consist of Voxel elements
Mat volume;
};
TSDFVolume::TSDFVolume(Point3i _res, float _voxelSize, Affine3f _pose, float _truncDist, int _maxWeight,
float _raycastStepFactor, bool zFirstMemOrder) :
voxelSize(_voxelSize),
voxelSizeInv(1.f/_voxelSize),
volResolution(_res),
maxWeight(_maxWeight),
pose(_pose),
raycastStepFactor(_raycastStepFactor)
TSDFVolume::TSDFVolume(float _voxelSize, Affine3f _pose, float _raycastStepFactor, float _truncDist,
int _maxWeight, Point3i _resolution, bool zFirstMemOrder)
: Volume(_voxelSize, _pose, _raycastStepFactor),
volResolution(_resolution),
maxWeight(_maxWeight)
{
// Unlike original code, this should work with any volume size
// Not only when (x,y,z % 32) == 0
volSize = Point3f(volResolution) * voxelSize;
volSize = Point3f(volResolution) * voxelSize;
truncDist = std::max(_truncDist, 2.1f * voxelSize);
// (xRes*yRes*zRes) array
@ -103,11 +55,13 @@ TSDFVolume::TSDFVolume(Point3i _res, float _voxelSize, Affine3f _pose, float _tr
}
// dimension in voxels, size in meters
TSDFVolumeCPU::TSDFVolumeCPU(Point3i _res, float _voxelSize, cv::Affine3f _pose, float _truncDist, int _maxWeight,
float _raycastStepFactor, bool zFirstMemOrder) :
TSDFVolume(_res, _voxelSize, _pose, _truncDist, _maxWeight, _raycastStepFactor, zFirstMemOrder)
TSDFVolumeCPU::TSDFVolumeCPU(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor,
float _truncDist, int _maxWeight, Point3i _resolution,
bool zFirstMemOrder)
: TSDFVolume(_voxelSize, _pose, _raycastStepFactor, _truncDist, _maxWeight, _resolution,
zFirstMemOrder)
{
volume = Mat(1, volResolution.x * volResolution.y * volResolution.z, rawType<Voxel>());
volume = Mat(1, volResolution.x * volResolution.y * volResolution.z, rawType<TsdfVoxel>());
reset();
}
@ -117,13 +71,32 @@ void TSDFVolumeCPU::reset()
{
CV_TRACE_FUNCTION();
volume.forEach<VecT>([](VecT& vv, const int* /* position */)
volume.forEach<VecTsdfVoxel>([](VecTsdfVoxel& vv, const int* /* position */)
{
Voxel& v = reinterpret_cast<Voxel&>(vv);
v.v = 0; v.weight = 0;
TsdfVoxel& v = reinterpret_cast<TsdfVoxel&>(vv);
v.tsdf = 0.0f; v.weight = 0;
});
}
TsdfVoxel TSDFVolumeCPU::at(const cv::Vec3i& volumeIdx) const
{
//! Out of bounds
if ((volumeIdx[0] >= volResolution.x || volumeIdx[0] < 0) ||
(volumeIdx[1] >= volResolution.y || volumeIdx[1] < 0) ||
(volumeIdx[2] >= volResolution.z || volumeIdx[2] < 0))
{
TsdfVoxel dummy;
dummy.tsdf = 1.0f;
dummy.weight = 0;
return dummy;
}
const TsdfVoxel* volData = volume.ptr<TsdfVoxel>();
int coordBase =
volumeIdx[0] * volDims[0] + volumeIdx[1] * volDims[1] + volumeIdx[2] * volDims[2];
return volData[coordBase];
}
// SIMD version of that code is manually inlined
#if !USE_INTRINSICS
static const bool fixMissingData = false;
@ -204,7 +177,7 @@ static inline depthType bilinearDepth(const Depth& m, cv::Point2f pt)
struct IntegrateInvoker : ParallelLoopBody
{
IntegrateInvoker(TSDFVolumeCPU& _volume, const Depth& _depth, Intr intrinsics, cv::Affine3f cameraPose,
IntegrateInvoker(TSDFVolumeCPU& _volume, const Depth& _depth, const Intr& intrinsics, const cv::Affine3f& cameraPose,
float depthFactor) :
ParallelLoopBody(),
volume(_volume),
@ -214,7 +187,7 @@ struct IntegrateInvoker : ParallelLoopBody
truncDistInv(1.f/_volume.truncDist),
dfac(1.f/depthFactor)
{
volDataStart = volume.volume.ptr<Voxel>();
volDataStart = volume.volume.ptr<TsdfVoxel>();
}
#if USE_INTRINSICS
@ -231,10 +204,10 @@ struct IntegrateInvoker : ParallelLoopBody
for(int x = range.start; x < range.end; x++)
{
Voxel* volDataX = volDataStart + x*volume.volDims[0];
TsdfVoxel* volDataX = volDataStart + x*volume.volDims[0];
for(int y = 0; y < volume.volResolution.y; y++)
{
Voxel* volDataY = volDataX + y*volume.volDims[1];
TsdfVoxel* volDataY = volDataX + y*volume.volDims[1];
// optimization of camSpace transformation (vector addition instead of matmul at each z)
Point3f basePt = vol2cam*(Point3f((float)x, (float)y, 0)*volume.voxelSize);
v_float32x4 camSpacePt(basePt.x, basePt.y, basePt.z, 0);
@ -246,19 +219,20 @@ struct IntegrateInvoker : ParallelLoopBody
if(zStepPt.z > 0)
{
startZ = baseZ;
endZ = volume.volResolution.z;
endZ = volume.volResolution.z;
}
else
{
startZ = 0;
endZ = baseZ;
endZ = baseZ;
}
}
else
{
if(basePt.z > 0)
if (basePt.z > 0)
{
startZ = 0; endZ = volume.volResolution.z;
startZ = 0;
endZ = volume.volResolution.z;
}
else
{
@ -267,7 +241,7 @@ struct IntegrateInvoker : ParallelLoopBody
}
}
startZ = max(0, startZ);
endZ = min(volume.volResolution.z, endZ);
endZ = min(volume.volResolution.z, endZ);
for(int z = startZ; z < endZ; z++)
{
// optimization of the following:
@ -276,7 +250,6 @@ struct IntegrateInvoker : ParallelLoopBody
camSpacePt += zStep;
float zCamSpace = v_reinterpret_as_f32(v_rotate_right<2>(v_reinterpret_as_u32(camSpacePt))).get0();
if(zCamSpace <= 0.f)
continue;
@ -298,11 +271,11 @@ struct IntegrateInvoker : ParallelLoopBody
continue;
// xi, yi = floor(pt)
v_int32x4 ip = v_floor(pt);
v_int32x4 ip = v_floor(pt);
v_int32x4 ipshift = ip;
int xi = ipshift.get0();
ipshift = v_rotate_right<1>(ipshift);
int yi = ipshift.get0();
int xi = ipshift.get0();
ipshift = v_rotate_right<1>(ipshift);
int yi = ipshift.get0();
const depthType* row0 = depth[yi+0];
const depthType* row1 = depth[yi+1];
@ -319,15 +292,15 @@ struct IntegrateInvoker : ParallelLoopBody
if(v_check_all(vall > v_setzero_f32()))
{
v_float32x4 t = pt - v_cvt_f32(ip);
float tx = t.get0();
float tx = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
v_float32x4 ty = v_setall_f32(t.get0());
// vx is y-interpolated between rows 0 and 1
v_float32x4 vx = v001 + ty*(v101 - v001);
float v0 = vx.get0();
float v0 = vx.get0();
vx = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(vx)));
float v1 = vx.get0();
v = v0 + tx*(v1 - v0);
v = v0 + tx*(v1 - v0);
}
else
continue;
@ -336,20 +309,19 @@ struct IntegrateInvoker : ParallelLoopBody
// norm(camPixVec) produces double which is too slow
float pixNorm = sqrt(v_reduce_sum(camPixVec*camPixVec));
// difference between distances of point and of surface to camera
volumeType sdf = pixNorm*(v*dfac - zCamSpace);
TsdfType sdf = pixNorm*(v*dfac - zCamSpace);
// possible alternative is:
// kftype sdf = norm(camSpacePt)*(v*dfac/camSpacePt.z - 1);
if(sdf >= -volume.truncDist)
{
volumeType tsdf = fmin(1.f, sdf * truncDistInv);
TsdfType tsdf = fmin(1.f, sdf * truncDistInv);
Voxel& voxel = volDataY[z*volume.volDims[2]];
int& weight = voxel.weight;
volumeType& value = voxel.v;
TsdfVoxel& voxel = volDataY[z*volume.volDims[2]];
int& weight = voxel.weight;
TsdfType& value = voxel.tsdf;
// update TSDF
value = (value*weight+tsdf) / (weight + 1);
value = (value*weight+tsdf) / (weight + 1);
weight = min(weight + 1, volume.maxWeight);
}
}
@ -361,18 +333,18 @@ struct IntegrateInvoker : ParallelLoopBody
{
for(int x = range.start; x < range.end; x++)
{
Voxel* volDataX = volDataStart + x*volume.volDims[0];
TsdfVoxel* volDataX = volDataStart + x*volume.volDims[0];
for(int y = 0; y < volume.volResolution.y; y++)
{
Voxel* volDataY = volDataX+y*volume.volDims[1];
TsdfVoxel* volDataY = volDataX+y*volume.volDims[1];
// optimization of camSpace transformation (vector addition instead of matmul at each z)
Point3f basePt = vol2cam*(Point3f(x, y, 0)*volume.voxelSize);
Point3f camSpacePt = basePt;
// zStep == vol2cam*(Point3f(x, y, 1)*voxelSize) - basePt;
// zStep == vol2cam*[Point3f(x, y, 1) - Point3f(x, y, 0)]*voxelSize
Point3f zStep = Point3f(vol2cam.matrix(0, 2),
vol2cam.matrix(1, 2),
vol2cam.matrix(2, 2))*volume.voxelSize;
int startZ, endZ;
if(abs(zStep.z) > 1e-5)
{
@ -380,19 +352,20 @@ struct IntegrateInvoker : ParallelLoopBody
if(zStep.z > 0)
{
startZ = baseZ;
endZ = volume.volResolution.z;
endZ = volume.volResolution.z;
}
else
{
startZ = 0;
endZ = baseZ;
endZ = baseZ;
}
}
else
{
if(basePt.z > 0)
{
startZ = 0; endZ = volume.volResolution.z;
startZ = 0;
endZ = volume.volResolution.z;
}
else
{
@ -401,14 +374,13 @@ struct IntegrateInvoker : ParallelLoopBody
}
}
startZ = max(0, startZ);
endZ = min(volume.volResolution.z, endZ);
endZ = min(volume.volResolution.z, endZ);
for(int z = startZ; z < endZ; z++)
{
// optimization of the following:
//Point3f volPt = Point3f(x, y, z)*volume.voxelSize;
//Point3f camSpacePt = vol2cam * volPt;
camSpacePt += zStep;
if(camSpacePt.z <= 0)
continue;
@ -422,20 +394,20 @@ struct IntegrateInvoker : ParallelLoopBody
// norm(camPixVec) produces double which is too slow
float pixNorm = sqrt(camPixVec.dot(camPixVec));
// difference between distances of point and of surface to camera
volumeType sdf = pixNorm*(v*dfac - camSpacePt.z);
TsdfType sdf = pixNorm*(v*dfac - camSpacePt.z);
// possible alternative is:
// kftype sdf = norm(camSpacePt)*(v*dfac/camSpacePt.z - 1);
if(sdf >= -volume.truncDist)
{
volumeType tsdf = fmin(1.f, sdf * truncDistInv);
TsdfType tsdf = fmin(1.f, sdf * truncDistInv);
Voxel& voxel = volDataY[z*volume.volDims[2]];
int& weight = voxel.weight;
volumeType& value = voxel.v;
TsdfVoxel& voxel = volDataY[z*volume.volDims[2]];
int& weight = voxel.weight;
TsdfType& value = voxel.tsdf;
// update TSDF
value = (value*weight+tsdf) / (weight + 1);
value = (value*weight+tsdf) / (weight + 1);
weight = min(weight + 1, volume.maxWeight);
}
}
@ -450,17 +422,17 @@ struct IntegrateInvoker : ParallelLoopBody
const cv::Affine3f vol2cam;
const float truncDistInv;
const float dfac;
Voxel* volDataStart;
TsdfVoxel* volDataStart;
};
// use depth instead of distance (optimization)
void TSDFVolumeCPU::integrate(InputArray _depth, float depthFactor, cv::Affine3f cameraPose, Intr intrinsics)
void TSDFVolumeCPU::integrate(InputArray _depth, float depthFactor, const cv::Affine3f& cameraPose,
const Intr& intrinsics)
{
CV_TRACE_FUNCTION();
CV_Assert(_depth.type() == DEPTH_TYPE);
Depth depth = _depth.getMat();
IntegrateInvoker ii(*this, depth, intrinsics, cameraPose, depthFactor);
Range range(0, volResolution.x);
parallel_for_(range, ii);
@ -468,35 +440,37 @@ void TSDFVolumeCPU::integrate(InputArray _depth, float depthFactor, cv::Affine3f
#if USE_INTRINSICS
// all coordinate checks should be done in inclosing cycle
inline volumeType TSDFVolumeCPU::interpolateVoxel(Point3f _p) const
inline TsdfType TSDFVolumeCPU::interpolateVoxel(Point3f _p) const
{
v_float32x4 p(_p.x, _p.y, _p.z, 0);
return interpolateVoxel(p);
}
inline volumeType TSDFVolumeCPU::interpolateVoxel(const v_float32x4& p) const
inline TsdfType TSDFVolumeCPU::interpolateVoxel(const v_float32x4& p) const
{
// tx, ty, tz = floor(p)
v_int32x4 ip = v_floor(p);
v_int32x4 ip = v_floor(p);
v_float32x4 t = p - v_cvt_f32(ip);
float tx = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
float ty = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
float tz = t.get0();
float tx = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
float ty = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
float tz = t.get0();
int xdim = volDims[0], ydim = volDims[1], zdim = volDims[2];
const Voxel* volData = volume.ptr<Voxel>();
const TsdfVoxel* volData = volume.ptr<TsdfVoxel>();
int ix = ip.get0(); ip = v_rotate_right<1>(ip);
int iy = ip.get0(); ip = v_rotate_right<1>(ip);
int ix = ip.get0();
ip = v_rotate_right<1>(ip);
int iy = ip.get0();
ip = v_rotate_right<1>(ip);
int iz = ip.get0();
int coordBase = ix*xdim + iy*ydim + iz*zdim;
volumeType vx[8];
TsdfType vx[8];
for(int i = 0; i < 8; i++)
vx[i] = volData[neighbourCoords[i] + coordBase].v;
vx[i] = volData[neighbourCoords[i] + coordBase].tsdf;
v_float32x4 v0246(vx[0], vx[2], vx[4], vx[6]), v1357(vx[1], vx[3], vx[5], vx[7]);
v_float32x4 vxx = v0246 + v_setall_f32(tz)*(v1357 - v0246);
@ -505,14 +479,14 @@ inline volumeType TSDFVolumeCPU::interpolateVoxel(const v_float32x4& p) const
v_float32x4 v01_11 = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(vxx)));
v_float32x4 v0_1 = v00_10 + v_setall_f32(ty)*(v01_11 - v00_10);
float v0 = v0_1.get0();
v0_1 = v_reinterpret_as_f32(v_rotate_right<2>(v_reinterpret_as_u32(v0_1)));
float v1 = v0_1.get0();
float v0 = v0_1.get0();
v0_1 = v_reinterpret_as_f32(v_rotate_right<2>(v_reinterpret_as_u32(v0_1)));
float v1 = v0_1.get0();
return v0 + tx*(v1 - v0);
}
#else
inline volumeType TSDFVolumeCPU::interpolateVoxel(Point3f p) const
inline TsdfType TSDFVolumeCPU::interpolateVoxel(Point3f p) const
{
int xdim = volDims[0], ydim = volDims[1], zdim = volDims[2];
@ -525,19 +499,19 @@ inline volumeType TSDFVolumeCPU::interpolateVoxel(Point3f p) const
float tz = p.z - iz;
int coordBase = ix*xdim + iy*ydim + iz*zdim;
const Voxel* volData = volume.ptr<Voxel>();
const TsdfVoxel* volData = volume.ptr<TsdfVoxel>();
volumeType vx[8];
TsdfType vx[8];
for(int i = 0; i < 8; i++)
vx[i] = volData[neighbourCoords[i] + coordBase].v;
vx[i] = volData[neighbourCoords[i] + coordBase].tsdf;
volumeType v00 = vx[0] + tz*(vx[1] - vx[0]);
volumeType v01 = vx[2] + tz*(vx[3] - vx[2]);
volumeType v10 = vx[4] + tz*(vx[5] - vx[4]);
volumeType v11 = vx[6] + tz*(vx[7] - vx[6]);
TsdfType v00 = vx[0] + tz*(vx[1] - vx[0]);
TsdfType v01 = vx[2] + tz*(vx[3] - vx[2]);
TsdfType v10 = vx[4] + tz*(vx[5] - vx[4]);
TsdfType v11 = vx[6] + tz*(vx[7] - vx[6]);
volumeType v0 = v00 + ty*(v01 - v00);
volumeType v1 = v10 + ty*(v11 - v10);
TsdfType v0 = v00 + ty*(v01 - v00);
TsdfType v1 = v10 + ty*(v11 - v10);
return v0 + tx*(v1 - v0);
}
@ -563,16 +537,16 @@ inline v_float32x4 TSDFVolumeCPU::getNormalVoxel(const v_float32x4& p) const
))
return nanv;
v_int32x4 ip = v_floor(p);
v_int32x4 ip = v_floor(p);
v_float32x4 t = p - v_cvt_f32(ip);
float tx = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
float ty = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
float tz = t.get0();
float tx = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
float ty = t.get0();
t = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(t)));
float tz = t.get0();
const int xdim = volDims[0], ydim = volDims[1], zdim = volDims[2];
const Voxel* volData = volume.ptr<Voxel>();
const TsdfVoxel* volData = volume.ptr<TsdfVoxel>();
int ix = ip.get0(); ip = v_rotate_right<1>(ip);
int iy = ip.get0(); ip = v_rotate_right<1>(ip);
@ -585,12 +559,12 @@ inline v_float32x4 TSDFVolumeCPU::getNormalVoxel(const v_float32x4& p) const
for(int c = 0; c < 3; c++)
{
const int dim = volDims[c];
float& nv = an[c];
float& nv = an[c];
volumeType vx[8];
TsdfType vx[8];
for(int i = 0; i < 8; i++)
vx[i] = volData[neighbourCoords[i] + coordBase + 1*dim].v -
volData[neighbourCoords[i] + coordBase - 1*dim].v;
vx[i] = volData[neighbourCoords[i] + coordBase + 1*dim].tsdf -
volData[neighbourCoords[i] + coordBase - 1*dim].tsdf;
v_float32x4 v0246(vx[0], vx[2], vx[4], vx[6]), v1357(vx[1], vx[3], vx[5], vx[7]);
v_float32x4 vxx = v0246 + v_setall_f32(tz)*(v1357 - v0246);
@ -599,14 +573,14 @@ inline v_float32x4 TSDFVolumeCPU::getNormalVoxel(const v_float32x4& p) const
v_float32x4 v01_11 = v_reinterpret_as_f32(v_rotate_right<1>(v_reinterpret_as_u32(vxx)));
v_float32x4 v0_1 = v00_10 + v_setall_f32(ty)*(v01_11 - v00_10);
float v0 = v0_1.get0();
v0_1 = v_reinterpret_as_f32(v_rotate_right<2>(v_reinterpret_as_u32(v0_1)));
float v1 = v0_1.get0();
float v0 = v0_1.get0();
v0_1 = v_reinterpret_as_f32(v_rotate_right<2>(v_reinterpret_as_u32(v0_1)));
float v1 = v0_1.get0();
nv = v0 + tx*(v1 - v0);
}
v_float32x4 n = v_load_aligned(an);
v_float32x4 n = v_load_aligned(an);
v_float32x4 invNorm = v_invsqrt(v_setall_f32(v_reduce_sum(n*n)));
return n*invNorm;
}
@ -614,7 +588,7 @@ inline v_float32x4 TSDFVolumeCPU::getNormalVoxel(const v_float32x4& p) const
inline Point3f TSDFVolumeCPU::getNormalVoxel(Point3f p) const
{
const int xdim = volDims[0], ydim = volDims[1], zdim = volDims[2];
const Voxel* volData = volume.ptr<Voxel>();
const TsdfVoxel* volData = volume.ptr<TsdfVoxel>();
if(p.x < 1 || p.x >= volResolution.x - 2 ||
p.y < 1 || p.y >= volResolution.y - 2 ||
@ -637,18 +611,18 @@ inline Point3f TSDFVolumeCPU::getNormalVoxel(Point3f p) const
const int dim = volDims[c];
float& nv = an[c];
volumeType vx[8];
TsdfType vx[8];
for(int i = 0; i < 8; i++)
vx[i] = volData[neighbourCoords[i] + coordBase + 1*dim].v -
volData[neighbourCoords[i] + coordBase - 1*dim].v;
vx[i] = volData[neighbourCoords[i] + coordBase + 1*dim].tsdf -
volData[neighbourCoords[i] + coordBase - 1*dim].tsdf;
volumeType v00 = vx[0] + tz*(vx[1] - vx[0]);
volumeType v01 = vx[2] + tz*(vx[3] - vx[2]);
volumeType v10 = vx[4] + tz*(vx[5] - vx[4]);
volumeType v11 = vx[6] + tz*(vx[7] - vx[6]);
TsdfType v00 = vx[0] + tz*(vx[1] - vx[0]);
TsdfType v01 = vx[2] + tz*(vx[3] - vx[2]);
TsdfType v10 = vx[4] + tz*(vx[5] - vx[4]);
TsdfType v11 = vx[6] + tz*(vx[7] - vx[6]);
volumeType v0 = v00 + ty*(v01 - v00);
volumeType v1 = v10 + ty*(v11 - v10);
TsdfType v0 = v00 + ty*(v01 - v00);
TsdfType v1 = v10 + ty*(v11 - v10);
nv = v0 + tx*(v1 - v0);
}
@ -657,11 +631,10 @@ inline Point3f TSDFVolumeCPU::getNormalVoxel(Point3f p) const
}
#endif
struct RaycastInvoker : ParallelLoopBody
{
RaycastInvoker(Points& _points, Normals& _normals, Affine3f cameraPose,
Intr intrinsics, const TSDFVolumeCPU& _volume) :
RaycastInvoker(Points& _points, Normals& _normals, const Affine3f& cameraPose,
const Intr& intrinsics, const TSDFVolumeCPU& _volume) :
ParallelLoopBody(),
points(_points),
normals(_normals),
@ -686,8 +659,8 @@ struct RaycastInvoker : ParallelLoopBody
const v_float32x4 vcxy(reproj.cx, reproj.cy, 0, 0);
const float (&cm)[16] = cam2vol.matrix.val;
const v_float32x4 camRot0(cm[0], cm[4], cm[ 8], 0);
const v_float32x4 camRot1(cm[1], cm[5], cm[ 9], 0);
const v_float32x4 camRot0(cm[0], cm[4], cm[8], 0);
const v_float32x4 camRot1(cm[1], cm[5], cm[9], 0);
const v_float32x4 camRot2(cm[2], cm[6], cm[10], 0);
const v_float32x4 camTrans(cm[3], cm[7], cm[11], 0);
@ -699,8 +672,8 @@ struct RaycastInvoker : ParallelLoopBody
volume.voxelSizeInv, 1.f);
const float (&vm)[16] = vol2cam.matrix.val;
const v_float32x4 volRot0(vm[0], vm[4], vm[ 8], 0);
const v_float32x4 volRot1(vm[1], vm[5], vm[ 9], 0);
const v_float32x4 volRot0(vm[0], vm[4], vm[8], 0);
const v_float32x4 volRot1(vm[1], vm[5], vm[9], 0);
const v_float32x4 volRot2(vm[2], vm[6], vm[10], 0);
const v_float32x4 volTrans(vm[3], vm[7], vm[11], 0);
@ -719,14 +692,14 @@ struct RaycastInvoker : ParallelLoopBody
// 1. reproject (x, y) on projecting plane where z = 1.f
v_float32x4 planed = (v_float32x4((float)x, (float)y, 0.f, 0.f) - vcxy)*vfxy;
planed = v_combine_low(planed, v_float32x4(1.f, 0.f, 0.f, 0.f));
planed = v_combine_low(planed, v_float32x4(1.f, 0.f, 0.f, 0.f));
// 2. rotate to volume space
planed = v_matmuladd(planed, camRot0, camRot1, camRot2, v_setzero_f32());
// 3. normalize
v_float32x4 invNorm = v_invsqrt(v_setall_f32(v_reduce_sum(planed*planed)));
v_float32x4 dir = planed*invNorm;
v_float32x4 dir = planed*invNorm;
// compute intersection of ray with all six bbox planes
v_float32x4 rayinv = v_setall_f32(1.f)/dir;
@ -740,8 +713,8 @@ struct RaycastInvoker : ParallelLoopBody
// near clipping plane
const float clip = 0.f;
float tmin = max(v_reduce_max(minAx), clip);
float tmax = v_reduce_min(maxAx);
float tmin = max(v_reduce_max(minAx), clip);
float tmax = v_reduce_min(maxAx);
// precautions against getting coordinates out of bounds
tmin = tmin + tstep;
@ -753,12 +726,12 @@ struct RaycastInvoker : ParallelLoopBody
orig *= invVoxelSize;
dir *= invVoxelSize;
int xdim = volume.volDims[0];
int ydim = volume.volDims[1];
int zdim = volume.volDims[2];
int xdim = volume.volDims[0];
int ydim = volume.volDims[1];
int zdim = volume.volDims[2];
v_float32x4 rayStep = dir * v_setall_f32(tstep);
v_float32x4 next = (orig + dir * v_setall_f32(tmin));
volumeType f = volume.interpolateVoxel(next), fnext = f;
v_float32x4 next = (orig + dir * v_setall_f32(tmin));
TsdfType f = volume.interpolateVoxel(next), fnext = f;
//raymarch
int steps = 0;
@ -772,7 +745,7 @@ struct RaycastInvoker : ParallelLoopBody
int iz = ip.get0();
int coord = ix*xdim + iy*ydim + iz*zdim;
fnext = volume.volume.at<Voxel>(coord).v;
fnext = volume.volume.at<TsdfVoxel>(coord).tsdf;
if(fnext != f)
{
fnext = volume.interpolateVoxel(next);
@ -790,10 +763,8 @@ struct RaycastInvoker : ParallelLoopBody
if(f > 0.f && fnext < 0.f)
{
v_float32x4 tp = next - rayStep;
volumeType ft = volume.interpolateVoxel(tp);
volumeType ftdt = volume.interpolateVoxel(next);
// float t = tmin + steps*tstep;
// float ts = t - tstep*ft/(ftdt - ft);
TsdfType ft = volume.interpolateVoxel(tp);
TsdfType ftdt = volume.interpolateVoxel(next);
float ts = tmin + tstep*(steps - ft/(ftdt - ft));
// avoid division by zero
@ -867,7 +838,7 @@ struct RaycastInvoker : ParallelLoopBody
Point3f rayStep = dir * tstep;
Point3f next = (orig + dir * tmin);
volumeType f = volume.interpolateVoxel(next), fnext = f;
TsdfType f = volume.interpolateVoxel(next), fnext = f;
//raymarch
int steps = 0;
@ -881,7 +852,7 @@ struct RaycastInvoker : ParallelLoopBody
int ix = cvRound(next.x);
int iy = cvRound(next.y);
int iz = cvRound(next.z);
fnext = volume.volume.at<Voxel>(ix*xdim + iy*ydim + iz*zdim).v;
fnext = volume.volume.at<TsdfVoxel>(ix*xdim + iy*ydim + iz*zdim).tsdf;
if(fnext != f)
{
fnext = volume.interpolateVoxel(next);
@ -898,9 +869,9 @@ struct RaycastInvoker : ParallelLoopBody
// linearly interpolate t between two f values
if(f > 0.f && fnext < 0.f)
{
Point3f tp = next - rayStep;
volumeType ft = volume.interpolateVoxel(tp);
volumeType ftdt = volume.interpolateVoxel(next);
Point3f tp = next - rayStep;
TsdfType ft = volume.interpolateVoxel(tp);
TsdfType ftdt = volume.interpolateVoxel(next);
// float t = tmin + steps*tstep;
// float ts = t - tstep*ft/(ftdt - ft);
float ts = tmin + tstep*(steps - ft/(ftdt - ft));
@ -944,7 +915,7 @@ struct RaycastInvoker : ParallelLoopBody
};
void TSDFVolumeCPU::raycast(cv::Affine3f cameraPose, Intr intrinsics, Size frameSize,
void TSDFVolumeCPU::raycast(const cv::Affine3f& cameraPose, const Intr& intrinsics, Size frameSize,
cv::OutputArray _points, cv::OutputArray _normals) const
{
CV_TRACE_FUNCTION();
@ -967,8 +938,8 @@ void TSDFVolumeCPU::raycast(cv::Affine3f cameraPose, Intr intrinsics, Size frame
struct FetchPointsNormalsInvoker : ParallelLoopBody
{
FetchPointsNormalsInvoker(const TSDFVolumeCPU& _volume,
std::vector< std::vector<ptype> >& _pVecs,
std::vector< std::vector<ptype> >& _nVecs,
std::vector<std::vector<ptype>>& _pVecs,
std::vector<std::vector<ptype>>& _nVecs,
bool _needNormals) :
ParallelLoopBody(),
vol(_volume),
@ -976,7 +947,7 @@ struct FetchPointsNormalsInvoker : ParallelLoopBody
nVecs(_nVecs),
needNormals(_needNormals)
{
volDataStart = vol.volume.ptr<Voxel>();
volDataStart = vol.volume.ptr<TsdfVoxel>();
}
inline void coord(std::vector<ptype>& points, std::vector<ptype>& normals,
@ -988,37 +959,37 @@ struct FetchPointsNormalsInvoker : ParallelLoopBody
float Vc = 0.f;
if(axis == 0)
{
shift = Point3i(1, 0, 0);
shift = Point3i(1, 0, 0);
limits = (x + 1 < vol.volResolution.x);
Vc = V.x;
Vc = V.x;
}
if(axis == 1)
{
shift = Point3i(0, 1, 0);
shift = Point3i(0, 1, 0);
limits = (y + 1 < vol.volResolution.y);
Vc = V.y;
Vc = V.y;
}
if(axis == 2)
{
shift = Point3i(0, 0, 1);
shift = Point3i(0, 0, 1);
limits = (z + 1 < vol.volResolution.z);
Vc = V.z;
Vc = V.z;
}
if(limits)
{
const Voxel& voxeld = volDataStart[(x+shift.x)*vol.volDims[0] +
(y+shift.y)*vol.volDims[1] +
(z+shift.z)*vol.volDims[2]];
volumeType vd = voxeld.v;
const TsdfVoxel& voxeld = volDataStart[(x+shift.x)*vol.volDims[0] +
(y+shift.y)*vol.volDims[1] +
(z+shift.z)*vol.volDims[2]];
TsdfType vd = voxeld.tsdf;
if(voxeld.weight != 0 && vd != 1.f)
{
if((v0 > 0 && vd < 0) || (v0 < 0 && vd > 0))
{
//linearly interpolate coordinate
float Vn = Vc + vol.voxelSize;
float dinv = 1.f/(abs(v0)+abs(vd));
float Vn = Vc + vol.voxelSize;
float dinv = 1.f/(abs(v0)+abs(vd));
float inter = (Vc*abs(vd) + Vn*abs(v0))*dinv;
Point3f p(shift.x ? inter : V.x,
@ -1040,14 +1011,14 @@ struct FetchPointsNormalsInvoker : ParallelLoopBody
std::vector<ptype> points, normals;
for(int x = range.start; x < range.end; x++)
{
const Voxel* volDataX = volDataStart + x*vol.volDims[0];
const TsdfVoxel* volDataX = volDataStart + x*vol.volDims[0];
for(int y = 0; y < vol.volResolution.y; y++)
{
const Voxel* volDataY = volDataX + y*vol.volDims[1];
const TsdfVoxel* volDataY = volDataX + y*vol.volDims[1];
for(int z = 0; z < vol.volResolution.z; z++)
{
const Voxel& voxel0 = volDataY[z*vol.volDims[2]];
volumeType v0 = voxel0.v;
const TsdfVoxel& voxel0 = volDataY[z*vol.volDims[2]];
TsdfType v0 = voxel0.tsdf;
if(voxel0.weight != 0 && v0 != 1.f)
{
Point3f V(Point3f((float)x + 0.5f, (float)y + 0.5f, (float)z + 0.5f)*vol.voxelSize);
@ -1067,9 +1038,9 @@ struct FetchPointsNormalsInvoker : ParallelLoopBody
}
const TSDFVolumeCPU& vol;
std::vector< std::vector<ptype> >& pVecs;
std::vector< std::vector<ptype> >& nVecs;
const Voxel* volDataStart;
std::vector<std::vector<ptype>>& pVecs;
std::vector<std::vector<ptype>>& nVecs;
const TsdfVoxel* volDataStart;
bool needNormals;
mutable Mutex mutex;
};
@ -1080,7 +1051,7 @@ void TSDFVolumeCPU::fetchPointsNormals(OutputArray _points, OutputArray _normals
if(_points.needed())
{
std::vector< std::vector<ptype> > pVecs, nVecs;
std::vector<std::vector<ptype>> pVecs, nVecs;
FetchPointsNormalsInvoker fi(*this, pVecs, nVecs, _normals.needed());
Range range(0, volResolution.x);
const int nstripes = -1;
@ -1105,10 +1076,9 @@ void TSDFVolumeCPU::fetchPointsNormals(OutputArray _points, OutputArray _normals
}
}
struct PushNormals
{
PushNormals(const TSDFVolumeCPU& _vol, Mat_<ptype>& _nrm) :
PushNormals(const TSDFVolumeCPU& _vol, Normals& _nrm) :
vol(_vol), normals(_nrm), invPose(vol.pose.inv())
{ }
void operator ()(const ptype &pp, const int * position) const
@ -1118,13 +1088,13 @@ struct PushNormals
if(!isNaN(p))
{
Point3f voxPt = (invPose * p);
voxPt = voxPt * vol.voxelSizeInv;
n = vol.pose.rotation() * vol.getNormalVoxel(voxPt);
voxPt = voxPt * vol.voxelSizeInv;
n = vol.pose.rotation() * vol.getNormalVoxel(voxPt);
}
normals(position[0], position[1]) = toPtype(n);
}
const TSDFVolumeCPU& vol;
Mat_<ptype>& normals;
Normals& normals;
Affine3f invPose;
};
@ -1140,7 +1110,7 @@ void TSDFVolumeCPU::fetchNormals(InputArray _points, OutputArray _normals) const
CV_Assert(points.type() == POINT_TYPE);
_normals.createSameSize(_points, _points.type());
Mat_<ptype> normals = _normals.getMat();
Normals normals = _normals.getMat();
points.forEach(PushNormals(*this, normals));
}
@ -1149,34 +1119,9 @@ void TSDFVolumeCPU::fetchNormals(InputArray _points, OutputArray _normals) const
///////// GPU implementation /////////
#ifdef HAVE_OPENCL
class TSDFVolumeGPU : public TSDFVolume
{
public:
// dimension in voxels, size in meters
TSDFVolumeGPU(Point3i _res, float _voxelSize, cv::Affine3f _pose, float _truncDist, int _maxWeight,
float _raycastStepFactor);
virtual void integrate(InputArray _depth, float depthFactor, cv::Affine3f cameraPose, cv::kinfu::Intr intrinsics) override;
virtual void raycast(cv::Affine3f cameraPose, cv::kinfu::Intr intrinsics, cv::Size frameSize,
cv::OutputArray _points, cv::OutputArray _normals) const override;
virtual void fetchPointsNormals(cv::OutputArray points, cv::OutputArray normals) const override;
virtual void fetchNormals(cv::InputArray points, cv::OutputArray normals) const override;
virtual void reset() override;
// See zFirstMemOrder arg of parent class constructor
// for the array layout info
// Array elem is CV_32FC2, read as (float, int)
// TODO: optimization possible to (fp16, int16), see Voxel definition
UMat volume;
};
TSDFVolumeGPU::TSDFVolumeGPU(Point3i _res, float _voxelSize, cv::Affine3f _pose, float _truncDist, int _maxWeight,
float _raycastStepFactor) :
TSDFVolume(_res, _voxelSize, _pose, _truncDist, _maxWeight, _raycastStepFactor, false)
TSDFVolumeGPU::TSDFVolumeGPU(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor, float _truncDist, int _maxWeight,
Point3i _resolution) :
TSDFVolume(_voxelSize, _pose, _raycastStepFactor, _truncDist, _maxWeight, _resolution, false)
{
volume = UMat(1, volResolution.x * volResolution.y * volResolution.z, CV_32FC2);
@ -1195,16 +1140,16 @@ void TSDFVolumeGPU::reset()
// use depth instead of distance (optimization)
void TSDFVolumeGPU::integrate(InputArray _depth, float depthFactor,
cv::Affine3f cameraPose, Intr intrinsics)
const cv::Affine3f& cameraPose, const Intr& intrinsics)
{
CV_TRACE_FUNCTION();
UMat depth = _depth.getUMat();
cv::String errorStr;
cv::String name = "integrate";
cv::String name = "integrate";
ocl::ProgramSource source = ocl::rgbd::tsdf_oclsrc;
cv::String options = "-cl-mad-enable";
cv::String options = "-cl-mad-enable";
ocl::Kernel k;
k.create(name.c_str(), source, options, &errorStr);
@ -1240,7 +1185,7 @@ void TSDFVolumeGPU::integrate(InputArray _depth, float depthFactor,
}
void TSDFVolumeGPU::raycast(cv::Affine3f cameraPose, Intr intrinsics, Size frameSize,
void TSDFVolumeGPU::raycast(const cv::Affine3f& cameraPose, const Intr& intrinsics, Size frameSize,
cv::OutputArray _points, cv::OutputArray _normals) const
{
CV_TRACE_FUNCTION();
@ -1248,9 +1193,9 @@ void TSDFVolumeGPU::raycast(cv::Affine3f cameraPose, Intr intrinsics, Size frame
CV_Assert(frameSize.area() > 0);
cv::String errorStr;
cv::String name = "raycast";
cv::String name = "raycast";
ocl::ProgramSource source = ocl::rgbd::tsdf_oclsrc;
cv::String options = "-cl-mad-enable";
cv::String options = "-cl-mad-enable";
ocl::Kernel k;
k.create(name.c_str(), source, options, &errorStr);
@ -1316,9 +1261,9 @@ void TSDFVolumeGPU::fetchNormals(InputArray _points, OutputArray _normals) const
UMat normals = _normals.getUMat();
cv::String errorStr;
cv::String name = "getNormals";
cv::String name = "getNormals";
ocl::ProgramSource source = ocl::rgbd::tsdf_oclsrc;
cv::String options = "-cl-mad-enable";
cv::String options = "-cl-mad-enable";
ocl::Kernel k;
k.create(name.c_str(), source, options, &errorStr);
@ -1365,7 +1310,7 @@ void TSDFVolumeGPU::fetchPointsNormals(OutputArray points, OutputArray normals)
cv::String errorStr;
ocl::ProgramSource source = ocl::rgbd::tsdf_oclsrc;
cv::String options = "-cl-mad-enable";
cv::String options = "-cl-mad-enable";
kscan.create("scanSize", source, options, &errorStr);
@ -1378,15 +1323,15 @@ void TSDFVolumeGPU::fetchPointsNormals(OutputArray points, OutputArray normals)
globalSize[2] = (size_t)volResolution.z;
const ocl::Device& device = ocl::Device::getDefault();
size_t wgsLimit = device.maxWorkGroupSize();
size_t memSize = device.localMemSize();
size_t wgsLimit = device.maxWorkGroupSize();
size_t memSize = device.localMemSize();
// local mem should keep a point (and a normal) for each thread in a group
// use 4 float per each point and normal
size_t elemSize = (sizeof(float)*4)*(needNormals ? 2 : 1);
const size_t lcols = 8;
const size_t lrows = 8;
size_t lplanes = min(memSize/elemSize, wgsLimit)/lcols/lrows;
lplanes = roundDownPow2(lplanes);
size_t elemSize = (sizeof(float)*4)*(needNormals ? 2 : 1);
const size_t lcols = 8;
const size_t lrows = 8;
size_t lplanes = min(memSize/elemSize, wgsLimit)/lcols/lrows;
lplanes = roundDownPow2(lplanes);
size_t localSize[3] = {lcols, lrows, lplanes};
Vec3i ngroups((int)divUp(globalSize[0], (unsigned int)localSize[0]),
(int)divUp(globalSize[1], (unsigned int)localSize[1]),
@ -1416,7 +1361,7 @@ void TSDFVolumeGPU::fetchPointsNormals(OutputArray points, OutputArray normals)
throw std::runtime_error("Failed to run kernel");
Mat groupedSumCpu = groupedSum.getMat(ACCESS_READ);
int gpuSum = (int)cv::sum(groupedSumCpu)[0];
int gpuSum = (int)cv::sum(groupedSumCpu)[0];
// should be no CPU copies when new kernel is executing
groupedSumCpu.release();
@ -1472,14 +1417,16 @@ void TSDFVolumeGPU::fetchPointsNormals(OutputArray points, OutputArray normals)
#endif
cv::Ptr<TSDFVolume> makeTSDFVolume(Point3i _res, float _voxelSize, cv::Affine3f _pose, float _truncDist, int _maxWeight,
float _raycastStepFactor)
cv::Ptr<TSDFVolume> makeTSDFVolume(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor,
float _truncDist, int _maxWeight, Point3i _resolution)
{
#ifdef HAVE_OPENCL
if(cv::ocl::useOpenCL())
return cv::makePtr<TSDFVolumeGPU>(_res, _voxelSize, _pose, _truncDist, _maxWeight, _raycastStepFactor);
if (cv::ocl::useOpenCL())
return cv::makePtr<TSDFVolumeGPU>(_voxelSize, _pose, _raycastStepFactor, _truncDist, _maxWeight,
_resolution);
#endif
return cv::makePtr<TSDFVolumeCPU>(_res, _voxelSize, _pose, _truncDist, _maxWeight, _raycastStepFactor);
return cv::makePtr<TSDFVolumeCPU>(_voxelSize, _pose, _raycastStepFactor, _truncDist, _maxWeight,
_resolution);
}
} // namespace kinfu

@ -7,36 +7,37 @@
#ifndef __OPENCV_KINFU_TSDF_H__
#define __OPENCV_KINFU_TSDF_H__
#include "kinfu_frame.hpp"
#include <opencv2/rgbd/volume.hpp>
namespace cv {
namespace kinfu {
#include "kinfu_frame.hpp"
#include "utils.hpp"
namespace cv
{
namespace kinfu
{
// TODO: Optimization possible:
// * TsdfType can be FP16
// * weight can be uint16
typedef float TsdfType;
struct TsdfVoxel
{
TsdfType tsdf;
int weight;
};
typedef Vec<uchar, sizeof(TsdfVoxel)> VecTsdfVoxel;
class TSDFVolume
class TSDFVolume : public Volume
{
public:
public:
// dimension in voxels, size in meters
TSDFVolume(Point3i _res, float _voxelSize, cv::Affine3f _pose, float _truncDist, int _maxWeight,
float _raycastStepFactor, bool zFirstMemOrder = true);
virtual void integrate(InputArray _depth, float depthFactor, cv::Affine3f cameraPose, cv::kinfu::Intr intrinsics) = 0;
virtual void raycast(cv::Affine3f cameraPose, cv::kinfu::Intr intrinsics, cv::Size frameSize,
cv::OutputArray points, cv::OutputArray normals) const = 0;
TSDFVolume(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor, float _truncDist,
int _maxWeight, Point3i _resolution, bool zFirstMemOrder = true);
virtual ~TSDFVolume() = default;
virtual void fetchPointsNormals(cv::OutputArray points, cv::OutputArray normals) const = 0;
virtual void fetchNormals(cv::InputArray points, cv::OutputArray _normals) const = 0;
virtual void reset() = 0;
virtual ~TSDFVolume() { }
float voxelSize;
float voxelSizeInv;
public:
Point3i volResolution;
int maxWeight;
cv::Affine3f pose;
float raycastStepFactor;
Point3f volSize;
float truncDist;
@ -44,9 +45,67 @@ public:
Vec8i neighbourCoords;
};
cv::Ptr<TSDFVolume> makeTSDFVolume(Point3i _res, float _voxelSize, cv::Affine3f _pose, float _truncDist, int _maxWeight,
float _raycastStepFactor);
class TSDFVolumeCPU : public TSDFVolume
{
public:
// dimension in voxels, size in meters
TSDFVolumeCPU(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor, float _truncDist,
int _maxWeight, Point3i _resolution, bool zFirstMemOrder = true);
virtual void integrate(InputArray _depth, float depthFactor, const cv::Affine3f& cameraPose,
const cv::kinfu::Intr& intrinsics) override;
virtual void raycast(const cv::Affine3f& cameraPose, const cv::kinfu::Intr& intrinsics,
cv::Size frameSize, cv::OutputArray points,
cv::OutputArray normals) const override;
virtual void fetchNormals(cv::InputArray points, cv::OutputArray _normals) const override;
virtual void fetchPointsNormals(cv::OutputArray points, cv::OutputArray normals) const override;
virtual void reset() override;
virtual TsdfVoxel at(const cv::Vec3i& volumeIdx) const;
TsdfType interpolateVoxel(cv::Point3f p) const;
Point3f getNormalVoxel(cv::Point3f p) const;
#if USE_INTRINSICS
TsdfType interpolateVoxel(const v_float32x4& p) const;
v_float32x4 getNormalVoxel(const v_float32x4& p) const;
#endif
// See zFirstMemOrder arg of parent class constructor
// for the array layout info
// Consist of Voxel elements
Mat volume;
};
#ifdef HAVE_OPENCL
class TSDFVolumeGPU : public TSDFVolume
{
public:
// dimension in voxels, size in meters
TSDFVolumeGPU(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor, float _truncDist,
int _maxWeight, Point3i _resolution);
virtual void integrate(InputArray _depth, float depthFactor, const cv::Affine3f& cameraPose,
const cv::kinfu::Intr& intrinsics) override;
virtual void raycast(const cv::Affine3f& cameraPose, const cv::kinfu::Intr& intrinsics,
cv::Size frameSize, cv::OutputArray _points,
cv::OutputArray _normals) const override;
virtual void fetchPointsNormals(cv::OutputArray points, cv::OutputArray normals) const override;
virtual void fetchNormals(cv::InputArray points, cv::OutputArray normals) const override;
} // namespace kinfu
} // namespace cv
virtual void reset() override;
// See zFirstMemOrder arg of parent class constructor
// for the array layout info
// Array elem is CV_32FC2, read as (float, int)
// TODO: optimization possible to (fp16, int16), see Voxel definition
UMat volume;
};
#endif
cv::Ptr<TSDFVolume> makeTSDFVolume(float _voxelSize, cv::Affine3f _pose, float _raycastStepFactor,
float _truncDist, int _maxWeight, Point3i _resolution);
} // namespace kinfu
} // namespace cv
#endif

@ -12,7 +12,7 @@
namespace cv
{
namespace rgbd
{
{
/** If the input image is of type CV_16UC1 (like the Kinect one), the image is converted to floats, divided
* by 1000 to get a depth in meters, and the values 0 are converted to std::numeric_limits<float>::quiet_NaN()
* Otherwise, the image is simply converted to floats

@ -22,6 +22,8 @@ namespace rgbd
* @param the desired output depth (floats or double)
* @param out The rescaled float depth image
*/
/* void rescaleDepth(InputArray in_in, int depth, OutputArray out_out); */
template<typename T>
void
rescaleDepthTemplated(const Mat& in, Mat& out);
@ -68,66 +70,6 @@ static inline bool isNaN(const cv::v_float32x4& p)
}
#endif
/** @brief Camera intrinsics */
struct Intr
{
/** Reprojects screen point to camera space given z coord. */
struct Reprojector
{
Reprojector() {}
inline Reprojector(Intr intr)
{
fxinv = 1.f/intr.fx, fyinv = 1.f/intr.fy;
cx = intr.cx, cy = intr.cy;
}
template<typename T>
inline cv::Point3_<T> operator()(cv::Point3_<T> p) const
{
T x = p.z * (p.x - cx) * fxinv;
T y = p.z * (p.y - cy) * fyinv;
return cv::Point3_<T>(x, y, p.z);
}
float fxinv, fyinv, cx, cy;
};
/** Projects camera space vector onto screen */
struct Projector
{
inline Projector(Intr intr) : fx(intr.fx), fy(intr.fy), cx(intr.cx), cy(intr.cy) { }
template<typename T>
inline cv::Point_<T> operator()(cv::Point3_<T> p) const
{
T invz = T(1)/p.z;
T x = fx*(p.x*invz) + cx;
T y = fy*(p.y*invz) + cy;
return cv::Point_<T>(x, y);
}
template<typename T>
inline cv::Point_<T> operator()(cv::Point3_<T> p, cv::Point3_<T>& pixVec) const
{
T invz = T(1)/p.z;
pixVec = cv::Point3_<T>(p.x*invz, p.y*invz, 1);
T x = fx*pixVec.x + cx;
T y = fy*pixVec.y + cy;
return cv::Point_<T>(x, y);
}
float fx, fy, cx, cy;
};
Intr() : fx(), fy(), cx(), cy() { }
Intr(float _fx, float _fy, float _cx, float _cy) : fx(_fx), fy(_fy), cx(_cx), cy(_cy) { }
Intr(cv::Matx33f m) : fx(m(0, 0)), fy(m(1, 1)), cx(m(0, 2)), cy(m(1, 2)) { }
// scale intrinsics to pyramid level
inline Intr scale(int pyr) const
{
float factor = (1.f /(1 << pyr));
return Intr(fx*factor, fy*factor, cx*factor, cy*factor);
}
inline Reprojector makeReprojector() const { return Reprojector(*this); }
inline Projector makeProjector() const { return Projector(*this); }
float fx, fy, cx, cy;
};
inline size_t roundDownPow2(size_t x)
{
size_t shift = 0;

@ -0,0 +1,34 @@
// This file is part of OpenCV project.
// It is subject to the license terms in the LICENSE file found in the top-level directory
// of this distribution and at http://opencv.org/license.html
#include "precomp.hpp"
#include <opencv2/rgbd/volume.hpp>
#include "tsdf.hpp"
#include "hash_tsdf.hpp"
namespace cv
{
namespace kinfu
{
cv::Ptr<Volume> makeVolume(VolumeType _volumeType, float _voxelSize, cv::Affine3f _pose,
float _raycastStepFactor, float _truncDist, int _maxWeight,
float _truncateThreshold, Point3i _resolution)
{
if (_volumeType == VolumeType::TSDF)
{
return makeTSDFVolume(_voxelSize, _pose, _raycastStepFactor, _truncDist, _maxWeight,
_resolution);
}
else if (_volumeType == VolumeType::HASHTSDF)
{
return makeHashTSDFVolume(_voxelSize, _pose, _raycastStepFactor, _truncDist, _maxWeight,
_truncateThreshold);
}
else
return nullptr;
}
} // namespace kinfu
} // namespace cv

@ -276,7 +276,7 @@ Ptr<Scene> Scene::create(int nScene, Size sz, Matx33f _intr, float _depthFactor)
static const bool display = false;
void flyTest(bool hiDense, bool inequal)
void flyTest(bool hiDense, bool inequal, bool hashTsdf = false)
{
Ptr<kinfu::Params> params;
if(hiDense)
@ -284,6 +284,9 @@ void flyTest(bool hiDense, bool inequal)
else
params = kinfu::Params::coarseParams();
if(hashTsdf)
params = kinfu::Params::hashTSDFParams(!hiDense);
if(inequal)
{
params->volumeDims[0] += 32;
@ -369,4 +372,10 @@ TEST(KinectFusion, DISABLED_OCL)
}
#endif
TEST( KinectFusion, DISABLED_hashTsdf )
{
flyTest(false, false, true);
//! hashTSDF does not support non-equal volumeDims
flyTest(true, false, true);
}
}} // namespace

Loading…
Cancel
Save