Create normals from triangulated surface

Question

Objective

I have a 3D face mesh triangulated. I have computed the midpoint of each triangle in the mesh, and I have the normal vector (n1, n2,n3) for each triangle. The objective is to create a 2D grid with [width x Height] resembling an image (with z assumed to be a constant). This is placed on top of the face mesh.

Refer to this image:

Figure 1

From these normal vectors, I then derive an RGB triplet for the resulting image.

My code is functionally correct, but it is very very slow. I cannot afford such huge time. I need help in optimizing my code in such a way that this would be computed faster and more efficient.

Things to replicate:

1. Download all the mesh related 3 files -Download_Files
2. Link those 3 files in the location correctly
3. Run the code

The slow block is highlighted with the comment Need Optimization HERE - VERY SLOW.

import os
import math
import numpy as np
import cv2
import matplotlib.pyplot as plt
def data_cleaning(filename):
 # Using readlines()
 data_list = []
 file1 = open(filename, 'r')
 Lines = file1.readlines()
 Data_dic = dict()
 for line in Lines:
 line = line.split("(")[-1].split(")")[0]
 line = line.split(",")
 no = int(line[0])
 x = float(line[1])
 y = float(line[2])
 z = float(line[3])
 Data_dic[no] = [x,y,z]
 # df = df.append({'No': no, 'X': x, 'Y': y, 'Z': z}, ignore_index=False)
 # print(no, x, y, z)
 data_list.append([x,y,z])
 return Data_dic, data_list
def distance_2points(P1, P2):
 p1 = np.array(P1)
 p2 = np.array(P2)
 squared_dist = np.sum((p1 - p2) ** 2, axis=0)
 dist = np.sqrt(squared_dist)
 return dist
def Pixel_Grid(Width, Height):
 # Grid Generation Position
 TOP_LEFT_X = -0.9296
 TOP_RIGHT_X = 1.053
 TOP_LEFT_Y = -0.8783
 BOTTOM_LEFT_Y = 1.311
 x = np.linspace(TOP_LEFT_X, TOP_RIGHT_X, Width)
 y = np.linspace(TOP_LEFT_Y, BOTTOM_LEFT_Y, Height)
 XX, YY = np.meshgrid(x, y, sparse=True)
 return XX, YY
def midpoint(x1, y1, x2, y2):
 x = (x1 + x2) / 2
 y = (y1 + y2) / 2
 return x, y
# Provide Meshgrid Sparse Data
def mid_points(xx, yy, w, h):
 Width = w
 Height = h
 xx = xx.reshape(Width, 1)
 l_x = len(xx)
 l_y = len(yy)
 all_midpoints = []
 for i in range(l_y):
 for j in range(l_x):
 if (j + 1 < Width):
 current_x = xx[j]
 next_x = xx[j + 1]
 if (i + 1 < Height):
 current_y = yy[i]
 next_y = yy[i + 1]
 x_value, y_value = midpoint(current_x, current_y, next_x, next_y)
 x_value = x_value.tolist()[0]
 y_value = y_value.tolist()[0]
 mid_point_xy = [x_value, y_value]
 all_midpoints.append(mid_point_xy)
 return all_midpoints
def Normal_map_fn(Mapping_dic, Mid_Points, indexing, Width, Height, Channels):
 Normal_map = np.zeros((Height, Width, Channels), np.float32)
 # ---------------------------------Need Optimization HERE - VERY SLOW ----------------------
 for e, i in enumerate(Mid_Points):
 p_y, p_x = i
 dis = []
 vala = []
 all = []
 for k, v in Mapping_dic.items():
 Mid_point_triangle = v[0]
 normal = v[1]
 M_x, M_y, M_z = Mid_point_triangle
 P1 = [p_x, p_y]
 P2 = [M_x, M_y]
 dist = distance_2points(P1, P2)
 dis.append(dist)
 vala.append(k)
 all.append([k, dist, normal])
 Min_value = min(dis)
 Index = np.argmin(dis)
 Key_Value = vala[Index]
 k, dist, normal = all[Index]
 val1, val2 = indexing[e]
 Normal_map[val1, val2] = (normal[0], normal[1], normal[2])
 dis = []
 vala = []
 print("done !", e)
 # ----------------------------------------------------------------------------------
 fig, ax = plt.subplots()
 cax = plt.imshow(Normal_map, cmap='gray')
 # cbar = fig.colorbar(cax)
 plt.show()
 name0 = os.path.join("Normal_map.png")
 cv2.imwrite(name0, Normal_map)
 return Normal_map
# Task 0 - Cleaning the Mesh, Normal and Triangle data
Vertex_Data = '../3D_Info/Vertex_Data.txt'
Triangle_Index = '../3D_Info/Mesh_Data.txt'
Normal_Vector = '../3D_Info/Face_Normal.txt'
Vertex, v_list = data_cleaning(Vertex_Data)
Triangle, t_list = data_cleaning(Triangle_Index)
Normal, n_list = data_cleaning(Normal_Vector)
# Task 1 - Removing Negative Normals and Triangles
Camera_Vector = [0, 0, 1]
l = len(Normal)
Negative_dic = dict()
angle = []
for i in range(l):
 # print(Normal[i])
 uv = np.dot(Normal[i], Camera_Vector)
 u_mag = np.sqrt(pow(Normal[i][0], 2) + pow(Normal[i][1], 2) + pow(Normal[i][2], 2))
 v_mag = np.sqrt(pow(Camera_Vector[0], 2) + pow(Camera_Vector[1], 2) + pow(Camera_Vector[2], 2))
 res = math.acos(uv/u_mag*v_mag) * (180.0 / math.pi)
 if (uv <= 0):
 Negative_dic[i] = Normal[i]
 angle.append(res)
nl = len(Negative_dic.items())
# Task 2 - Delete all the negative normals
for k, v in Negative_dic.items():
 del Triangle[k]
 del Normal[k]
# Task 3 - Mapping
Mapping_dic = dict()
for k, v in Triangle.items():
 one, two, three = v
 point_m = Vertex[one]
 point_n = Vertex[two]
 point_o = Vertex[three]
 mid_x = (point_m[0] + point_n[0] + point_o[0]) / 3
 mid_y = (point_m[1] + point_n[1] + point_o[1]) / 3
 mid_z = (point_m[2] + point_n[2] + point_o[2]) / 3
 Mid_point_triangle = [mid_x, mid_y, mid_z]
 # point_list.append(Mid_point_triangle)
 Mapping_dic[k] = [Mid_point_triangle, Normal[k]]
# Task 4- Baking Normal Map
Width = 600
Height = 800
Channels = 3
indexing = dict()
count = 0
for j in range(Height):
 for i in range(Width):
 indexing[count] = [j, i]
 count = count + 1
XX, YY = Pixel_Grid(Width, Height)
Mid_Points = mid_points(XX, YY, Width, Height)
Normal = Normal_map_fn(Mapping_dic, Mid_Points, indexing, Width, Height, Channels)

• Reinderien: Working on it. You should post excerpts (perhaps 20 lines each) of your three files inline here.

Answer 1

It's not only the "slow block" you indicated that needs refactoring; it's basically everything. As a beginner, assume that you may never use for loops when writing Numpy code otherwise you're going to land in deep trouble.

Add PEP484 type hints.

Your data_cleaning is poorly named because it doesn't clean data, per se - it loads data. Also, none of that code needs to exist if you make one call to np.loadtxt.

distance_2points needs to go away. If you actually needed to calculate distance yourself, call np.linalg.norm, but as you'll see, you don't.

If there's one method here that's almost basically kind of OK, it's Pixel_Grid (which should be named pixel_grid).

midpoint needs to go away and be vectorised.

mid_points has a logic problem. Your boundary condition doesn't make any sense. From a 600x800 image, if you're doing a mean filter, you can only meaningfully produce a 599x799 image.

Don't incrementally append to lists when writing Numpy code.

Broadly speaking, the algorithm in Normal_map_fn is a true disaster. Even if it were properly vectorised - which it very much isn't - it's about as brute force as conceivably possible. This is where Google helps you - use one of the spatial index features of Scipy. You're effectively making a Voronoi diagram, but the actual Voronoi interface is inconvenient to use so just use a KDTree instead.

Each of your "tasks" would be well-suited to a separate function. For ease of comprehension I have not changed their names, but you should.

There's a lot of dead code that can be deleted, such as your acos() calls.

This code does not call for a dictionary. Dictionaries and lists in a Numpy context are usually code smells.

Do not call into OpenCV just to create an image; use PIL instead. This will require that you normalise the array.

Your left, right, top and bottom boundary names are not well-chosen, because they describe corners when they should describe edges.

Suggested

I am not patient enough to run your code to completion for comparison, but this stands a reasonable chance of producing something matching your intent. It runs in about one second.

from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import KDTree
from PIL import Image
def data_cleaning(stem: str, dtype: object = np.float64) -> np.ndarray:
 parent = Path('../3D_Info')
 filename = (parent / stem).with_suffix('.txt')
 return np.loadtxt(filename, usecols=(1, 2, 3), delimiter=', ', comments=')', dtype=dtype)
def pixel_grid(width: int, height: int) -> tuple[
 np.ndarray, # xx
 np.ndarray, # yy
]:
 # Grid Generation Position
 left = -0.9296
 right = 1.053
 top = -0.8783
 bottom = 1.311
 x = np.linspace(left, right, width)
 y = np.linspace(top, bottom, height)
 return np.meshgrid(x, y, sparse=True)
def mid_points(xx: np.ndarray, yy: np.ndarray) -> np.ndarray:
 """Provide Meshgrid Sparse Data"""
 xmid = (xx[0, :-1] + xx[0, 1:])/2
 ymid = (yy[:-1] + yy[1:])[:, 0]/2
 return np.stack(
 np.meshgrid(ymid, xmid)
 ).T
def map_normals(
 triangle_midpoints: np.ndarray, # 7149 x 2
 normals: np.ndarray, # 7149 x 3
 midpoints: np.ndarray, # 799 x 599 x 2
):
 tree = KDTree(triangle_midpoints)
 _, idx = tree.query(x=midpoints, workers=-1) # 799 x 599
 return normals[idx, :] # 799 x 599 x 3
def task_0() -> tuple[np.ndarray, np.ndarray, np.ndarray]:
 """Cleaning the Mesh, Normal and Triangle data"""
 vertices = data_cleaning('Vertex_Data')
 triangles = data_cleaning('Mesh_Data', dtype=int)
 normals = data_cleaning('Face_Normal')
 return vertices, triangles, normals
def task_1(normals: np.ndarray) -> np.ndarray:
 """Removing Negative Normals and Triangles"""
 camera_vector = [0, 0, 1]
 return np.dot(normals, camera_vector) > 0
def task_2(triangles: np.ndarray, normals: np.ndarray, positive: np.ndarray) -> tuple[
 np.ndarray, # triangles
 np.ndarray, # normals
]:
 """Delete all the negative normals"""
 new_indices = np.cumsum(positive) - 1
 new_triangles = new_indices[triangles[positive, :]]
 new_normals = normals[positive, :]
 return new_triangles, new_normals
def task_3(triangles: np.ndarray, vertices: np.ndarray) -> np.ndarray:
 """Mapping"""
 yx = vertices[triangles, 1::-1]
 return yx.mean(axis=1)
def task_4(triangle_midpoints: np.ndarray, normals: np.ndarray):
 """Baking normals Map"""
 xx, yy = pixel_grid(width=600, height=800)
 midpoints = mid_points(xx, yy)
 return map_normals(triangle_midpoints, normals, midpoints)
def main() -> None:
 vertices, triangles, normals = task_0()
 positive = task_1(normals)
 triangles, normals = task_2(triangles, normals, positive)
 triangle_midpoints = task_3(triangles, vertices)
 normal_map = task_4(triangle_midpoints, normals)
 normalised = ((normal_map + 1) * (255/2)).astype(np.uint8)
 image = Image.fromarray(normalised, mode='RGB')
 image.save("normal_map.png")
 fig, ax = plt.subplots()
 ax.imshow(normal_map, cmap='gray')
 plt.show()
if __name__ == '__main__':
 main()

voronoi