python中深度图的表面法线计算答案

【问题标题】：Surface normal calculation from depth map in pythonpython中深度图的表面法线计算
【发布时间】：2018-11-17 10:36:47
【问题描述】：

我尝试在 python 中实现以下 c++ 代码：

depth.convertTo(depth, CV_64FC1); // I do not know why it is needed to be 
transformed to 64bit image my input is 32bit

Mat nor(depth.size(), CV_64FC3);

for(int x = 1; x < depth.cols - 1; ++x)
{
   for(int y = 1; y < depth.rows - 1; ++y)
   {
      Vec3d t(x,y-1,depth.at<double>(y-1, x)/*depth(y-1,x)*/);
      Vec3d l(x-1,y,depth.at<double>(y, x-1)/*depth(y,x-1)*/);
      Vec3d c(x,y,depth.at<double>(y, x)/*depth(y,x)*/);
      Vec3d d = (l-c).cross(t-c);
      Vec3d n = normalize(d);
      nor.at<Vec3d>(y,x) = n;
   }
}

imshow("normals", nor);

python 代码：

d_im = cv2.imread("depth.jpg")
d_im = d_im.astype("float64")

normals = np.array(d_im, dtype="float32")
h,w,d = d_im.shape
for i in range(1,w-1):
  for j in range(1,h-1):
    t = np.array([i,j-1,d_im[j-1,i,0]],dtype="float64")
    f = np.array([i-1,j,d_im[j,i-1,0]],dtype="float64")
    c = np.array([i,j,d_im[j,i,0]] , dtype = "float64")
    d = np.cross(f-c,t-c)
    n = d / np.sqrt((np.sum(d**2)))
    normals[j,i,:] = n

cv2.imwrite("normal.jpg",normals*255)

输入图像：

c++代码输出：

我的python代码输出：

我找不到这些差异的原因。如何使用 python 获取 c++ 代码输出？

【问题讨论】：

实际上当我测试它时，我的 Python 输出看起来很合理；您安装了哪些版本的使用模块？ (pip list --local)
opencv-python 3.4.3.18，numpy 1.15.2
好的，我的版本是：cv2：3.4.3，numpy：1.13.1。所以我不认为这是问题所在。我不得不承认，my image 看起来仍然没有你的 C++ 图像那么流畅。我有一个问题：在 Python 中，你用d_im[j-1,i,0] 构造你的向量；在 C++ 中你写 depth.at<double>(y-1, x) 为什么在 C++ 中不需要三个索引？`
好吧，在修改了代码并最终对图像进行了超级采样之后，我想到了查看原始图像的深度。我从您那里下载的图像是 8 位图像！因此，您在 Python 脚本中看到的是非常低的位深度和可怕的 jpeg 伪影的组合。您应该尝试获得更好的测试图像
这也解释了为什么我确实得到了更好的结果：SO 上的图像是 png；所以我没有 jpeg 人工制品

标签： python opencv normals

【解决方案1】：

正如 user8408080 所说，您的输出似乎有由 jpeg 格式引起的伪影。另请记住，将 8 位图像作为深度图导入与直接使用深度图矩阵的结果不同。

关于您的 Python 代码，我的建议是使用矢量化函数并尽可能避免循环（这很慢）。

zy, zx = np.gradient(d_im)  
# You may also consider using Sobel to get a joint Gaussian smoothing and differentation
# to reduce noise
#zx = cv2.Sobel(d_im, cv2.CV_64F, 1, 0, ksize=5)     
#zy = cv2.Sobel(d_im, cv2.CV_64F, 0, 1, ksize=5)

normal = np.dstack((-zx, -zy, np.ones_like(d_im)))
n = np.linalg.norm(normal, axis=2)
normal[:, :, 0] /= n
normal[:, :, 1] /= n
normal[:, :, 2] /= n

# offset and rescale values to be in 0-255
normal += 1
normal /= 2
normal *= 255

cv2.imwrite("normal.png", normal[:, :, ::-1])

【讨论】：

非常好的答案，因为它还指出了将 C/C++ 转换为 Python 时出现的问题。还有一点旁注：如果您尝试只使用imshow 图像，只需注释掉normal *= 255
感谢矢量化版本。但是，在this post 中，他们也使用了深度图图像，并且得到了更好的结果。
如您发布的链接所述：他们拥有的原始图像是 32 位的。他们只上传了 8 位，所以我们永远无法重现他们在那里所做的事情
哦，对不起，我没有考虑到这一点。感谢您的所有努力。
为什么第三维是一？ normal = np.dstack((-zx, -zy, np.ones_like(d_im)))

【解决方案2】：

代码（矩阵计算）应该是：

def normalization(data):
   mo_chang =np.sqrt(np.multiply(data[:,:,0],data[:,:,0])+np.multiply(data[:,:,1],data[:,:,1])+np.multiply(data[:,:,2],data[:,:,2]))
   mo_chang = np.dstack((mo_chang,mo_chang,mo_chang))
   return data/mo_chang

x,y=np.meshgrid(np.arange(0,width),np.arange(0,height))
x=x.reshape([-1])
y=y.reshape([-1])
xyz=np.vstack((x,y,np.ones_like(x)))
pts_3d=np.dot(np.linalg.inv(K),xyz*img1_depth.reshape([-1]))
pts_3d_world=pts_3d.reshape((3,height,width))
f= pts_3d_world[:,1:height-1,2:width]-pts_3d_world[:,1:height-1,1:width-1]
t= pts_3d_world[:,2:height,1:width-1]-pts_3d_world[:,1:height-1,1:width-1]
normal_map=np.cross(f,l,axisa=0,axisb=0)
normal_map=normalization(normal_map)
normal_map=normal_map*0.5+0.5
alpha = np.full((height-2,width-2,1), (1.), dtype="float32")
normal_map=np.concatenate((normal_map,alpha),axis=2)

我们应该使用名为“K”的相机内部函数。我认为 f 和 t 的值是基于相机坐标中的 3D 点。
对于法线向量，(-1,-1,100) 和 (255,255,100) 在 8 位图像中是相同的颜色，但它们是完全不同的法线。所以我们应该通过normal_map=normal_map*0.5+0.5将法线值映射到(0,1)。

【讨论】：