使用相邻像素叉积计算深度图像的表面法线

Question

正如标题所说，我想通过使用相邻像素的叉积来计算给定深度图像的表面法线。我想为此使用 Opencv 并避免使用 PCL，但是我并不真正了解该过程，因为我在该主题上的知识非常有限。因此，如果有人可以提供一些提示，我将不胜感激。这里提一下，除了深度图和对应的rgb图，我没有其他信息，所以没有K相机矩阵信息。

因此，假设我们有以下深度图像：

我想在对应点找到具有对应深度值的法线向量，如下图所示：

如何使用相邻像素的叉积来做到这一点？如果法线不是很准确，我不介意。

谢谢。

---

更新：

好的，我试图按照@timday 的回答将他的代码移植到 Opencv。使用以下代码：

Mat depth = <my_depth_image> of type CV_32FC1
Mat normals(depth.size(), CV_32FC3);

for(int x = 0; x < depth.rows; ++x)
{
    for(int y = 0; y < depth.cols; ++y)
    {

        float dzdx = (depth.at<float>(x+1, y) - depth.at<float>(x-1, y)) / 2.0;
        float dzdy = (depth.at<float>(x, y+1) - depth.at<float>(x, y-1)) / 2.0;

        Vec3f d(-dzdx, -dzdy, 1.0f);
        Vec3f n = normalize(d);

        normals.at<Vec3f>(x, y) = n;
    }
}

imshow("depth", depth / 255);
imshow("normals", normals);

我得到了正确的以下结果（我必须将 double 替换为 float 和 Vecd 到 Vecf，但我不知道为什么会有所不同）：

Answer 1

我认为对的代码（矩阵计算）：

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)

1. 我们应该使用名为“K”的相机内部函数。我认为 f 和 t 的值是基于相机坐标中的 3D 点。

2. 对于法线向量，(-1,-1,100) 和 (255,255,100) 在 8 位图像中是相同的颜色，但它们是完全不同的法线。所以我们应该通过normal_map=normal_map*0.5+0.5将正常值映射到(0,1)。

欢迎交流。

Answer 2

您实际上并不需要为此使用叉积，但请参见下文。

考虑你的范围图像是一个函数 z(x,y)。

表面的法线方向为 (-dz/dx,-dz/dy,1)。 （这里的 dz/dx 是指微分：z 与 x 的变化率）。然后将法线按惯例归一化为单位长度。

顺便说一句，如果您想知道 (-dz/dx,-dz/dy,1) 来自哪里...如果您将平面中的 2 个正交切向量与 x 和 y 轴平行，那么那些是 (1,0,dzdx) 和 (0,1,dzdy)。法线垂直于切线，因此应该是 (1,0,dzdx)X(0,1,dzdy) - 其中“X”是叉积 - 即 (-dzdx,-dzdy,1)。所以有你的叉积派生法线，但是当你可以直接将结果表达式用于法线时，几乎不需要在代码中如此明确地计算它。

在 (x,y) 处计算单位长度法线的伪代码类似于

dzdx=(z(x+1,y)-z(x-1,y))/2.0;
dzdy=(z(x,y+1)-z(x,y-1))/2.0;
direction=(-dzdx,-dzdy,1.0)
magnitude=sqrt(direction.x**2 + direction.y**2 + direction.z**2)
normal=direction/magnitude

根据您要执行的操作，将 NaN 值替换为一些较大的数字可能更有意义。

使用这种方法，从您的范围图像中，我可以得到：

（然后，我使用计算出的法线方向来进行一些简单的着色；请注意由于范围图像的量化导致的“阶梯式”外观；理想情况下，对于实际范围数据，您的精度应高于 8 位）。

抱歉，不是 OpenCV 或 C++ 代码，只是为了完整性：生成该图像的完整代码（嵌入在 Qt QML 文件中的 GLSL；可以使用 Qt5 的 qmlscene 运行）如下。上面的伪代码可以在片段着色器的main()函数中找到：

import QtQuick 2.2

Image {
  source: 'range.png'  // The provided image

  ShaderEffect {
    anchors.fill: parent
    blending: false

    property real dx: 1.0/parent.width
    property real dy: 1.0/parent.height
    property variant src: parent

    vertexShader: "
      uniform highp mat4 qt_Matrix;
      attribute highp vec4 qt_Vertex;
      attribute highp vec2 qt_MultiTexCoord0;
      varying highp vec2 coord;
      void main() {
        coord=qt_MultiTexCoord0;
        gl_Position=qt_Matrix*qt_Vertex;
      }"

   fragmentShader: "
     uniform highp float dx;
     uniform highp float dy;
     varying highp vec2 coord;
     uniform sampler2D src;
     void main() {
       highp float dzdx=( texture2D(src,coord+vec2(dx,0.0)).x - texture2D(src,coord+vec2(-dx,0.0)).x )/(2.0*dx);
       highp float dzdy=( texture2D(src,coord+vec2(0.0,dy)).x - texture2D(src,coord+vec2(0.0,-dy)).x )/(2.0*dy);
       highp vec3 d=vec3(-dzdx,-dzdy,1.0);
       highp vec3 n=normalize(d);
       highp vec3 lightDirection=vec3(1.0,-2.0,3.0);
       highp float shading=0.5+0.5*dot(n,normalize(lightDirection));
       gl_FragColor=vec4(shading,shading,shading,1.0);
     }"
  }
}