Surface normal calculation from depth map in python

Question

I try to implement following c++ code in python:

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 code:

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)

input image:

c++ code output:

my python code output:

i can not find reason of these differences. How i can get c++ code output with python?

Actually when I test it, my Python output looks reasonable; what versions of the used modules do you have installed? (pip list --local) — user8408080
– user8408080, Commented Nov 17, 2018 at 11:19
okay, my versions are: cv2: 3.4.3, numpy: 1.13.1. So I don't think, that this is the problem. I have to admit, that my image still doesn't look as smooth as your C++ image. I have one question: In Python you construct your vectors with d_im[j-1,i,0]; in C++ you write depth.at<double>(y-1, x) why don't you need three indices in C++?` — user8408080
– user8408080, Commented Nov 17, 2018 at 11:35
Okay after tinkering around with the code and eventually supersampling the image, I had the idea to look what depth the original image has. The image I downloaded from you is an 8-bit image! So what you are seeing there with your Python script is a combination of very low bit depth and horrible jpeg artefacts. You should try an get a better testing image — user8408080
– user8408080, Commented Nov 17, 2018 at 12:12
This also explains why I did get a better result: images on SO are png; so there were no jpeg artefacts for me — user8408080
– user8408080, Commented Nov 17, 2018 at 12:16

keineahnung2345 · Accepted Answer · 2021-03-05 08:26:34Z

6

As user8408080 said you output seems to have artifacts caused by the jpeg format. Also keep in mind that importing an 8-bit image as a depth map will not give the same results than using directly the depth map matrix.

Regarding your python code, my advice would be to use vectorized functions and avoid loops as much as you can (it's very slow).

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])

edited Mar 5, 2021 at 8:26

keineahnung2345

2,6914 gold badges16 silver badges30 bronze badges

answered Nov 17, 2018 at 12:45

sowlosc

4803 silver badges9 bronze badges

Sign up to request clarification or add additional context in comments.

6 Comments

user8408080 Over a year ago

Very good answer, as it also points out the problems that come up when converting C/C++ to Python. Also a little sidenote: If you try to only imshow the image, just comment out the normal *= 255

Bedrick Kiq Over a year ago

Thanks for vectorized version. However, in this post , they also used depth-map image and they got better results.

user8408080 Over a year ago

As stated on the link you posted: The original image they had was 32-bit. They only uploaded 8-bit, so we can never reproduce what they did there

Bedrick Kiq Over a year ago

oh, sorry, i didn't consider this. Thank for your all efforts.

Mohammad Over a year ago

why is the third dimension is one? normal = np.dstack((-zx, -zy, np.ones_like(d_im)))

|

Yang · Accepted Answer · 2023-01-29 12:46:13Z

@sowlosc's answer don't take into account the camera intrinsics
@Baichuan's answer has obvious bugs

I implemented a brief, correct, fast, easy-to-use function get_surface_normal_by_depthnd compared with the revised @Baichuan's scheme, the correctness of the results was verified.

import cv2
import numpy as np

def get_surface_normal_by_depth(depth, K=None):
    """
    depth: (h, w) of float, the unit of depth is meter
    K: (3, 3) of float, the depth camere's intrinsic
    """
    K = [[1, 0], [0, 1]] if K is None else K
    fx, fy = K[0][0], K[1][1]

    dz_dv, dz_du = np.gradient(depth)  # u, v mean the pixel coordinate in the image
    # u*depth = fx*x + cx --> du/dx = fx / depth
    du_dx = fx / depth  # x is xyz of camera coordinate
    dv_dy = fy / depth

    dz_dx = dz_du * du_dx
    dz_dy = dz_dv * dv_dy
    # cross-product (1,0,dz_dx)X(0,1,dz_dy) = (-dz_dx, -dz_dy, 1)
    normal_cross = np.dstack((-dz_dx, -dz_dy, np.ones_like(depth)))
    # normalize to unit vector
    normal_unit = normal_cross / np.linalg.norm(normal_cross, axis=2, keepdims=True)
    # set default normal to [0, 0, 1]
    normal_unit[~np.isfinite(normal_unit).all(2)] = [0, 0, 1]
    return normal_unit

def get_normal_map_by_point_cloud(depth, K):
    height, width = depth.shape

    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 * 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, t, axisa=0, axisb=0)
    normal_map = normalization(normal_map)
    return normal_map


vis_normal = lambda normal: np.uint8((normal + 1) / 2 * 255)[..., ::-1]

normal1 = get_surface_normal_by_depth(depth, K)    #  spend time: 60ms
normal2 = get_normal_map_by_point_cloud(depth, K)  #  spend time: 90ms

cv2.imwrite("normal1.png", vis_normal(normal1))
cv2.imwrite("normal2.png", vis_normal(normal2))

The result: normal1.png and normal2.png

Baichuan · Accepted Answer · 2020-08-23 11:43:23Z

The code (matrix calculation) should be:

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)

We should use the camera intrinsics named 'K'. I think the value f and t is based on 3D points in camera coordinate.
For the normal vector, the (-1,-1,100) and (255,255,100) are the same color in 8 bites images but they are totally different normal. So we should map the normal values to (0,1) by normal_map=normal_map*0.5+0.5.

Collectives™ on Stack Overflow

Surface normal calculation from depth map in python

3 Answers 3

6 Comments

2 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

6 Comments

2 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related