From ed409742086f7e1264f33bf12543482a1f26d88b Mon Sep 17 00:00:00 2001 From: YadiraF Date: Tue, 9 Oct 2018 19:05:29 +0800 Subject: [PATCH] change structure --- examples/1_pipeline.py | 16 +- examples/2_3dmm.py | 7 +- examples/3_transform.py | 4 +- examples/4_light.py | 6 +- examples/5_render.py | 65 ++-- examples/6_generate_image_map.py | 9 +- examples/7_generate_uv_map.py | 7 +- examples/8_generate_posmap_300WLP.py | 3 +- examples/Data/BFM/readme.md | 2 +- face3d/__init__.py | 1 - face3d/mesh/__init__.py | 5 +- face3d/mesh/io.py | 38 +- face3d/mesh/light.py | 14 +- face3d/mesh/render.py | 291 ++++----------- face3d/mesh/transform.py | 9 + face3d/mesh_cython/__init__.py | 4 - face3d/mesh_cython/mesh_core.cpp | 328 ----------------- face3d/mesh_cython/mesh_core.h | 76 ---- face3d/mesh_cython/mesh_core_cython.pyx | 93 ----- face3d/mesh_cython/render.py | 131 ------- face3d/mesh_cython/setup.py | 30 -- face3d/mesh_numpy/__init__.py | 10 + face3d/mesh_numpy/io.py | 170 +++++++++ face3d/{mesh_cython => mesh_numpy}/light.py | 18 +- face3d/mesh_numpy/render.py | 287 +++++++++++++++ face3d/mesh_numpy/transform.py | 385 ++++++++++++++++++++ face3d/mesh_numpy/vis.py | 24 ++ 27 files changed, 1054 insertions(+), 979 deletions(-) delete mode 100644 face3d/mesh_cython/__init__.py delete mode 100644 face3d/mesh_cython/mesh_core.cpp delete mode 100644 face3d/mesh_cython/mesh_core.h delete mode 100644 face3d/mesh_cython/mesh_core_cython.pyx delete mode 100644 face3d/mesh_cython/render.py delete mode 100644 face3d/mesh_cython/setup.py create mode 100644 face3d/mesh_numpy/__init__.py create mode 100644 face3d/mesh_numpy/io.py rename face3d/{mesh_cython => mesh_numpy}/light.py (92%) create mode 100644 face3d/mesh_numpy/render.py create mode 100644 face3d/mesh_numpy/transform.py create mode 100644 face3d/mesh_numpy/vis.py diff --git a/examples/1_pipeline.py b/examples/1_pipeline.py index a0bc6e7..30c1b11 100644 --- a/examples/1_pipeline.py +++ b/examples/1_pipeline.py @@ -11,7 +11,6 @@ sys.path.append('..') import face3d from face3d import mesh -from face3d import mesh_cython # ------------------------------ 1. load mesh data # -- mesh data consists of: vertices, triangles, color(optinal), texture(optional) @@ -50,11 +49,16 @@ # change to image coords for rendering image_vertices = mesh.transform.to_image(projected_vertices, h, w) # render -rendering = face3d.mesh_cython.render.render_colors(image_vertices, triangles, lit_colors, h, w) +rendering = mesh.render.render_colors(image_vertices, triangles, lit_colors, h, w) # ---- show rendering -plt.imshow(rendering) -plt.show() +# plt.imshow(rendering) +# plt.show() +save_folder = 'results/pipeline' +if not os.path.exists(save_folder): + os.mkdir(save_folder) +io.imsave('{}/rendering.jpg'.format(save_folder), rendering) + # ---- show mesh -mesh.vis.plot_mesh(camera_vertices, triangles) -plt.show() +# mesh.vis.plot_mesh(camera_vertices, triangles) +# plt.show() diff --git a/examples/2_3dmm.py b/examples/2_3dmm.py index ffe1e95..bcf750d 100644 --- a/examples/2_3dmm.py +++ b/examples/2_3dmm.py @@ -13,7 +13,6 @@ sys.path.append('..') import face3d from face3d import mesh -from face3d import mesh_cython from face3d.morphable_model import MorphabelModel # --------------------- Forward: parameters(shape, expression, pose) --> 3D obj --> 2D image --------------- @@ -41,7 +40,7 @@ # set prop of rendering h = w = 256; c = 3 image_vertices = mesh.transform.to_image(projected_vertices, h, w) -image = mesh_cython.render.render_colors(image_vertices, bfm.triangles, colors, h, w) +image = mesh.render.render_colors(image_vertices, bfm.triangles, colors, h, w) # -------------------- Back: 2D image points and corresponding 3D vertex indices--> parameters(pose, shape, expression) ------ ## only use 68 key points to fit @@ -56,7 +55,7 @@ transformed_vertices = bfm.transform(fitted_vertices, fitted_s, fitted_angles, fitted_t) image_vertices = mesh.transform.to_image(transformed_vertices, h, w) -fitted_image = mesh_cython.render.render_colors(image_vertices, bfm.triangles, colors, h, w) +fitted_image = mesh.render.render_colors(image_vertices, bfm.triangles, colors, h, w) # ------------- print & show @@ -81,7 +80,7 @@ transformed_vertices = bfm.transform(fitted_vertices, fitted_s[i], fitted_angles[i], fitted_t[i]) image_vertices = mesh.transform.to_image(transformed_vertices, h, w) - fitted_image = mesh_cython.render.render_colors(image_vertices, bfm.triangles, colors, h, w) + fitted_image = mesh.render.render_colors(image_vertices, bfm.triangles, colors, h, w) io.imsave('{}/show_{:0>2d}.jpg'.format(save_folder, i), fitted_image) options = '-delay 20 -loop 0 -layers optimize' # gif. need ImageMagick. diff --git a/examples/3_transform.py b/examples/3_transform.py index fdb8f38..79d1b86 100644 --- a/examples/3_transform.py +++ b/examples/3_transform.py @@ -11,8 +11,6 @@ sys.path.append('..') import face3d from face3d import mesh -from face3d import mesh_cython - def transform_test(vertices, obj, camera, h = 256, w = 256): ''' @@ -35,7 +33,7 @@ def transform_test(vertices, obj, camera, h = 256, w = 256): ## to image coords(position in image) image_vertices = mesh.transform.to_image(projected_vertices, h, w, True) - rendering = mesh_cython.render.render_colors(image_vertices, triangles, colors, h, w) + rendering = mesh.render.render_colors(image_vertices, triangles, colors, h, w) rendering = np.minimum((np.maximum(rendering, 0)), 1) return rendering diff --git a/examples/4_light.py b/examples/4_light.py index 297a5b4..ae7535c 100644 --- a/examples/4_light.py +++ b/examples/4_light.py @@ -11,13 +11,11 @@ sys.path.append('..') import face3d from face3d import mesh -from face3d import mesh_cython - def light_test(vertices, light_positions, light_intensities, h = 256, w = 256): - lit_colors = mesh_cython.light.add_light(vertices, triangles, colors, light_positions, light_intensities) + lit_colors = mesh.light.add_light(vertices, triangles, colors, light_positions, light_intensities) image_vertices = mesh.transform.to_image(vertices, h, w) - rendering = mesh_cython.render.render_colors(image_vertices, triangles, lit_colors, h, w) + rendering = mesh.render.render_colors(image_vertices, triangles, lit_colors, h, w) rendering = np.minimum((np.maximum(rendering, 0)), 1) return rendering diff --git a/examples/5_render.py b/examples/5_render.py index aad35fe..9b4757d 100644 --- a/examples/5_render.py +++ b/examples/5_render.py @@ -12,7 +12,7 @@ sys.path.append('..') import face3d from face3d import mesh -from face3d import mesh_cython +from face3d import mesh_numpy from face3d.morphable_model import MorphabelModel # load data @@ -30,19 +30,19 @@ image_vertices = mesh.transform.to_image(vertices, h, w) # -----------------------------------------render colors -# render colors python +# render colors python numpy st = time() -rendering_cp = face3d.mesh.render.render_colors(image_vertices, triangles, colors, h, w) +rendering_cp = mesh_numpy.render.render_colors(image_vertices, triangles, colors, h, w) print('----------colors python: ', time() - st) -# render colors python ras +# render colors python ras numpy st = time() -rendering_cpr = face3d.mesh.render.render_colors_ras(image_vertices, triangles, colors, h, w) +rendering_cpr = mesh_numpy.render.render_colors_ras(image_vertices, triangles, colors, h, w) print('----------colors python ras: ', time() - st) # render colors python c++ st = time() -rendering_cc = face3d.mesh_cython.render.render_colors(image_vertices, triangles, colors, h, w) +rendering_cc = mesh.render.render_colors(image_vertices, triangles, colors, h, w) print('----------colors c++: ', time() - st) # ----------------------------------------- render texture(texture mapping) @@ -51,30 +51,51 @@ texture = io.imread('Data/uv_texture_map.jpg')/255. tex_h, tex_w, _ = texture.shape # load texcoord(uv coords) -uv_coords = face3d.morphable_model.load.load_uv_coords('Data/BFM/Out/BFM_UV.mat') # -uv_coords[:,0] = uv_coords[:,0]*(tex_h - 1) -uv_coords[:,1] = uv_coords[:,1]*(tex_w - 1) -uv_coords[:,1] = tex_w - uv_coords[:,1] - 1 -texcoord = np.hstack((uv_coords, np.zeros((uv_coords.shape[0], 1)))) # add z# render texture python +uv_coords = face3d.morphable_model.load.load_uv_coords('Data/BFM/Out/BFM_UV.mat') # general UV coords: range [0-1] +# to texture size +texcoord = np.zeros_like(uv_coords) +texcoord[:,0] = uv_coords[:,0]*(tex_h - 1) +texcoord[:,1] = uv_coords[:,1]*(tex_w - 1) +texcoord[:,1] = tex_w - texcoord[:,1] - 1 +texcoord = np.hstack((texcoord, np.zeros((texcoord.shape[0], 1)))) # add z# render texture python # tex_triangles tex_triangles = triangles # render texture python st = time() -rendering_tp = face3d.mesh.render.render_texture(image_vertices, triangles, texture, texcoord, tex_triangles, h, w) +rendering_tp = face3d.mesh_numpy.render.render_texture(image_vertices, triangles, texture, texcoord, tex_triangles, h, w) print('----------texture python: ', time() - st) # render texture c++ st = time() -rendering_tc = face3d.mesh_cython.render.render_texture(image_vertices, triangles, texture, texcoord, tex_triangles, h, w) +rendering_tc = face3d.mesh.render.render_texture(image_vertices, triangles, texture, texcoord, tex_triangles, h, w) print('----------texture c++: ', time() - st) -plt.subplot(2,2,1) -plt.imshow(rendering_cp) -plt.subplot(2,2,2) -plt.imshow(rendering_cc) -plt.subplot(2,2,3) -plt.imshow(rendering_tp) -plt.subplot(2,2,4) -plt.imshow(rendering_tc) -plt.show() +# plot +# plt.subplot(2,2,1) +# plt.imshow(rendering_cp) +# plt.subplot(2,2,2) +# plt.imshow(rendering_cc) +# plt.subplot(2,2,3) +# plt.imshow(rendering_tp) +# plt.subplot(2,2,4) +# plt.imshow(rendering_tc) +# plt.show() + + +## --------------------------- write obj +# mesh.vis.plot_mesh(vertices, triangles) +# plt.show() +save_folder = os.path.join('results', 'io') +if not os.path.exists(save_folder): + os.mkdir(save_folder) + +vertices, colors, uv_coords = image_vertices.astype(np.float32).copy(), colors.astype(np.float32).copy(), uv_coords.astype(np.float32).copy() + +st = time() +mesh_numpy.io.write_obj_with_colors_texture(os.path.join(save_folder, 'numpy.obj'), vertices, triangles, colors, texture, uv_coords) +print('----------write obj numpy: ', time() - st) + +st = time() +mesh.io.write_obj_with_colors_texture(os.path.join(save_folder, 'cython.obj'), vertices, triangles, colors, texture, uv_coords) +print('----------write obj cython: ', time() - st) diff --git a/examples/6_generate_image_map.py b/examples/6_generate_image_map.py index 0ef2bcc..552e8f6 100644 --- a/examples/6_generate_image_map.py +++ b/examples/6_generate_image_map.py @@ -12,7 +12,6 @@ sys.path.append('..') import face3d from face3d import mesh -from face3d import mesh_cython # ------------------------------ load mesh data C = sio.loadmat('Data/example1.mat') @@ -41,7 +40,7 @@ ## 0. color map attribute = colors -color_image = mesh_cython.render.render_colors(image_vertices, triangles, attribute, h, w, c=3) +color_image = mesh.render.render_colors(image_vertices, triangles, attribute, h, w, c=3) io.imsave('{}/color.jpg'.format(save_folder), np.squeeze(color_image)) ## 1. depth map @@ -49,19 +48,19 @@ z = z - np.min(z) z = z/np.max(z) attribute = z -depth_image = mesh_cython.render.render_colors(image_vertices, triangles, attribute, h, w, c=1) +depth_image = mesh.render.render_colors(image_vertices, triangles, attribute, h, w, c=1) io.imsave('{}/depth.jpg'.format(save_folder), np.squeeze(depth_image)) ## 2. pncc in 'Face Alignment Across Large Poses: A 3D Solution'. for dense correspondences pncc = face3d.morphable_model.load.load_pncc_code('Data/BFM/Out/pncc_code.mat') attribute = pncc -pncc_image = mesh_cython.render.render_colors(image_vertices, triangles, attribute, h, w, c=3) +pncc_image = mesh.render.render_colors(image_vertices, triangles, attribute, h, w, c=3) io.imsave('{}/pncc.jpg'.format(save_folder), np.squeeze(pncc_image)) ## 3. uv coordinates in 'DenseReg: Fully convolutional dense shape regression in-the-wild'. for dense correspondences uv_coords = face3d.morphable_model.load.load_uv_coords('Data/BFM/Out/BFM_UV.mat') # attribute = uv_coords # note that: original paper used quantized coords, here not -uv_coords_image = mesh_cython.render.render_colors(image_vertices, triangles, attribute, h, w, c=2) # two channels: u, v +uv_coords_image = mesh.render.render_colors(image_vertices, triangles, attribute, h, w, c=2) # two channels: u, v # add one channel for show uv_coords_image = np.concatenate((np.zeros((h, w, 1)), uv_coords_image), 2) io.imsave('{}/uv_coords.jpg'.format(save_folder), np.squeeze(uv_coords_image)) diff --git a/examples/7_generate_uv_map.py b/examples/7_generate_uv_map.py index f8aff04..3fcf6a9 100644 --- a/examples/7_generate_uv_map.py +++ b/examples/7_generate_uv_map.py @@ -13,7 +13,6 @@ sys.path.append('..') import face3d from face3d import mesh -from face3d import mesh_cython from face3d.morphable_model import MorphabelModel def process_uv(uv_coords, uv_h = 256, uv_w = 256): @@ -47,7 +46,7 @@ def process_uv(uv_coords, uv_h = 256, uv_w = 256): #-- 1. uv texture map attribute = colors -uv_texture_map = mesh_cython.render.render_colors(uv_coords, triangles, attribute, uv_h, uv_w, c=3) +uv_texture_map = mesh.render.render_colors(uv_coords, triangles, attribute, uv_h, uv_w, c=3) io.imsave('{}/uv_texture_map.jpg'.format(save_folder), np.squeeze(uv_texture_map)) #-- 2. uv position map in 'Joint 3D Face Reconstruction and Dense Alignment with Position Map Regression Network' @@ -61,8 +60,8 @@ def process_uv(uv_coords, uv_h = 256, uv_w = 256): position[:,2] = position[:,2] - np.min(position[:,2]) # translate z attribute = position # corresponding 2d facial image -image = mesh_cython.render.render_colors(image_vertices, triangles, colors, image_h, image_w, c=3) -uv_position_map = mesh_cython.render.render_colors(uv_coords, triangles, attribute, uv_h, uv_w, c=3) +image = mesh.render.render_colors(image_vertices, triangles, colors, image_h, image_w, c=3) +uv_position_map = mesh.render.render_colors(uv_coords, triangles, attribute, uv_h, uv_w, c=3) io.imsave('{}/image.jpg'.format(save_folder), np.squeeze(image)) np.save('{}/uv_position_map.npy'.format(save_folder), uv_position_map) io.imsave('{}/uv_position_map.jpg'.format(save_folder), (uv_position_map)/max(image_h, image_w)) # only for show diff --git a/examples/8_generate_posmap_300WLP.py b/examples/8_generate_posmap_300WLP.py index 0bc4055..3d48579 100644 --- a/examples/8_generate_posmap_300WLP.py +++ b/examples/8_generate_posmap_300WLP.py @@ -12,7 +12,6 @@ sys.path.append('..') import face3d from face3d import mesh -from face3d import mesh_cython from face3d.morphable_model import MorphabelModel def process_uv(uv_coords, uv_h = 256, uv_w = 256): @@ -76,7 +75,7 @@ def run_posmap_300W_LP(bfm, image_path, mat_path, save_folder, uv_h = 256, uv_w position[:, 2] = position[:, 2] - np.min(position[:, 2]) # translate z # 4. uv position map: render position in uv space - uv_position_map = mesh_cython.render.render_colors(uv_coords, bfm.full_triangles, position, uv_h, uv_w, c = 3) + uv_position_map = mesh.render.render_colors(uv_coords, bfm.full_triangles, position, uv_h, uv_w, c = 3) # 5. save files io.imsave('{}/{}'.format(save_folder, image_name), np.squeeze(cropped_image)) diff --git a/examples/Data/BFM/readme.md b/examples/Data/BFM/readme.md index 1eeded1..4485bce 100644 --- a/examples/Data/BFM/readme.md +++ b/examples/Data/BFM/readme.md @@ -2,7 +2,7 @@ 1. Download raw BFM model - website: https://faces.cs.unibas.ch/bfm/index.php?nav=1-0&id=basel_face_model. + website: https://faces.dmi.unibas.ch/bfm/main.php?nav=1-2&id=downloads copy 01_MorphabelModel.mat to raw/ diff --git a/face3d/__init__.py b/face3d/__init__.py index 4d49a30..6db1992 100644 --- a/face3d/__init__.py +++ b/face3d/__init__.py @@ -1,3 +1,2 @@ from . import mesh -from . import mesh_cython from . import morphable_model diff --git a/face3d/mesh/__init__.py b/face3d/mesh/__init__.py index 65d5039..7f8b3f4 100644 --- a/face3d/mesh/__init__.py +++ b/face3d/mesh/__init__.py @@ -1,7 +1,4 @@ -from __future__ import absolute_import -from __future__ import division -from __future__ import print_function - +from .cython import mesh_core_cython from . import io from . import vis from . import transform diff --git a/face3d/mesh/io.py b/face3d/mesh/io.py index 58a1b3e..4bb1460 100644 --- a/face3d/mesh/io.py +++ b/face3d/mesh/io.py @@ -3,7 +3,11 @@ from __future__ import print_function import numpy as np +import os +from skimage import io +from time import time +from .cython import mesh_core_cython ## TODO ## TODO: c++ version @@ -104,15 +108,15 @@ def write_obj_with_texture(obj_name, vertices, triangles, texture, uv_coords): # write texture as png imsave(texture_name, texture) - -def write_obj_with_colors_texture(obj_name, vertices, colors, triangles, texture, uv_coords): +# c++ version +def write_obj_with_colors_texture(obj_name, vertices, triangles, colors, texture, uv_coords): ''' Save 3D face model with texture. Ref: https://github.com/patrikhuber/eos/blob/bd00155ebae4b1a13b08bf5a991694d682abbada/include/eos/core/Mesh.hpp Args: obj_name: str vertices: shape = (nver, 3) - colors: shape = (nver, 3) triangles: shape = (ntri, 3) + colors: shape = (nver, 3) texture: shape = (256,256,3) uv_coords: shape = (nver, 3) max value<=1 ''' @@ -125,29 +129,9 @@ def write_obj_with_colors_texture(obj_name, vertices, colors, triangles, texture triangles += 1 # mesh lab start with 1 # write obj - with open(obj_name, 'w') as f: - # first line: write mtlib(material library) - s = "mtllib {}\n".format(os.path.abspath(mtl_name)) - f.write(s) - - # write vertices - for i in range(vertices.shape[0]): - s = 'v {} {} {} {} {} {}\n'.format(vertices[i, 0], vertices[i, 1], vertices[i, 2], colors[i, 0], colors[i, 1], colors[i, 2]) - f.write(s) - - # write uv coords - for i in range(uv_coords.shape[0]): - s = 'vt {} {}\n'.format(uv_coords[i,0], 1 - uv_coords[i,1]) - f.write(s) - - f.write("usemtl FaceTexture\n") - - # write f: ver ind/ uv ind - for i in range(triangles.shape[0]): - # s = 'f {}/{} {}/{} {}/{}\n'.format(triangles[i,0], triangles[i,0], triangles[i,1], triangles[i,1], triangles[i,2], triangles[i,2]) - s = 'f {}/{} {}/{} {}/{}\n'.format(triangles[i,2], triangles[i,2], triangles[i,1], triangles[i,1], triangles[i,0], triangles[i,0]) - f.write(s) - + vertices, colors, uv_coords = vertices.astype(np.float32).copy(), colors.astype(np.float32).copy(), uv_coords.astype(np.float32).copy() + mesh_core_cython.write_obj_with_colors_texture_core(str.encode(obj_name), str.encode(os.path.abspath(mtl_name)), vertices, triangles, colors, uv_coords, vertices.shape[0], triangles.shape[0], uv_coords.shape[0]) + # write mtl with open(mtl_name, 'w') as f: f.write("newmtl FaceTexture\n") @@ -155,4 +139,4 @@ def write_obj_with_colors_texture(obj_name, vertices, colors, triangles, texture f.write(s) # write texture as png - imsave(texture_name, texture) \ No newline at end of file + io.imsave(texture_name, texture) \ No newline at end of file diff --git a/face3d/mesh/light.py b/face3d/mesh/light.py index 5463fc5..5f0da63 100644 --- a/face3d/mesh/light.py +++ b/face3d/mesh/light.py @@ -9,6 +9,7 @@ from __future__ import print_function import numpy as np +from .cython import mesh_core_cython def get_normal(vertices, triangles): ''' calculate normal direction in each vertex @@ -23,12 +24,13 @@ def get_normal(vertices, triangles): pt2 = vertices[triangles[:, 2], :] # [ntri, 3] tri_normal = np.cross(pt0 - pt1, pt0 - pt2) # [ntri, 3]. normal of each triangle - normal = np.zeros_like(vertices) # [nver, 3] - for i in range(triangles.shape[0]): - normal[triangles[i, 0], :] = normal[triangles[i, 0], :] + tri_normal[i, :] - normal[triangles[i, 1], :] = normal[triangles[i, 1], :] + tri_normal[i, :] - normal[triangles[i, 2], :] = normal[triangles[i, 2], :] + tri_normal[i, :] - + normal = np.zeros_like(vertices, dtype = np.float32).copy() # [nver, 3] + # for i in range(triangles.shape[0]): + # normal[triangles[i, 0], :] = normal[triangles[i, 0], :] + tri_normal[i, :] + # normal[triangles[i, 1], :] = normal[triangles[i, 1], :] + tri_normal[i, :] + # normal[triangles[i, 2], :] = normal[triangles[i, 2], :] + tri_normal[i, :] + mesh_core_cython.get_normal_core(normal, tri_normal.astype(np.float32).copy(), triangles.copy(), triangles.shape[0]) + # normalize to unit length mag = np.sum(normal**2, 1) # [nver] zero_ind = (mag == 0) diff --git a/face3d/mesh/render.py b/face3d/mesh/render.py index efae442..1995722 100644 --- a/face3d/mesh/render.py +++ b/face3d/mesh/render.py @@ -5,10 +5,10 @@ 1. Generally, render func includes camera, light, raterize. Here no camera and light(I write these in other files) 2. Generally, the input vertices are normalized to [-1,1] and cetered on [0, 0]. (in world space) Here, the vertices are using image coords, which centers on [w/2, h/2] with the y-axis pointing to oppisite direction. -Means: render here only conducts interpolation.(I just want to make the input flexible) + Means: render here only conducts interpolation.(I just want to make the input flexible) -Author: YadiraF -Mail: fengyao@sjtu.edu.cn +Author: Yao Feng +Mail: yaofeng1995@gmail.com ''' from __future__ import absolute_import from __future__ import division @@ -17,82 +17,7 @@ import numpy as np from time import time -def isPointInTri(point, tri_points): - ''' Judge whether the point is in the triangle - Method: - http://blackpawn.com/texts/pointinpoly/ - Args: - point: (2,). [u, v] or [x, y] - tri_points: (3 vertices, 2 coords). three vertices(2d points) of a triangle. - Returns: - bool: true for in triangle - ''' - tp = tri_points - - # vectors - v0 = tp[2,:] - tp[0,:] - v1 = tp[1,:] - tp[0,:] - v2 = point - tp[0,:] - - # dot products - dot00 = np.dot(v0.T, v0) - dot01 = np.dot(v0.T, v1) - dot02 = np.dot(v0.T, v2) - dot11 = np.dot(v1.T, v1) - dot12 = np.dot(v1.T, v2) - - # barycentric coordinates - if dot00*dot11 - dot01*dot01 == 0: - inverDeno = 0 - else: - inverDeno = 1/(dot00*dot11 - dot01*dot01) - - u = (dot11*dot02 - dot01*dot12)*inverDeno - v = (dot00*dot12 - dot01*dot02)*inverDeno - - # check if point in triangle - return (u >= 0) & (v >= 0) & (u + v < 1) - -def get_point_weight(point, tri_points): - ''' Get the weights of the position - Methods: https://gamedev.stackexchange.com/questions/23743/whats-the-most-efficient-way-to-find-barycentric-coordinates - -m1.compute the area of the triangles formed by embedding the point P inside the triangle - -m2.Christer Ericson's book "Real-Time Collision Detection". faster.(used) - Args: - point: (2,). [u, v] or [x, y] - tri_points: (3 vertices, 2 coords). three vertices(2d points) of a triangle. - Returns: - w0: weight of v0 - w1: weight of v1 - w2: weight of v3 - ''' - tp = tri_points - # vectors - v0 = tp[2,:] - tp[0,:] - v1 = tp[1,:] - tp[0,:] - v2 = point - tp[0,:] - - # dot products - dot00 = np.dot(v0.T, v0) - dot01 = np.dot(v0.T, v1) - dot02 = np.dot(v0.T, v2) - dot11 = np.dot(v1.T, v1) - dot12 = np.dot(v1.T, v2) - - # barycentric coordinates - if dot00*dot11 - dot01*dot01 == 0: - inverDeno = 0 - else: - inverDeno = 1/(dot00*dot11 - dot01*dot01) - - u = (dot11*dot02 - dot01*dot12)*inverDeno - v = (dot00*dot12 - dot01*dot02)*inverDeno - - w0 = 1 - u - v - w1 = v - w2 = u - - return w0, w1, w2 +from .cython import mesh_core_cython def rasterize_triangles(vertices, triangles, h, w): ''' @@ -109,116 +34,63 @@ def rasterize_triangles(vertices, triangles, h, w): # Each triangle has 3 vertices & Each vertex has 3 coordinates x, y, z. # h, w is the size of rendering ''' + # initial - depth_buffer = np.zeros([h, w]) - 999999. #+ np.min(vertices[2,:]) - 999999. # set the initial z to the farest position + depth_buffer = np.zeros([h, w]) - 999999. #set the initial z to the farest position triangle_buffer = np.zeros([h, w], dtype = np.int32) - 1 # if tri id = -1, the pixel has no triangle correspondance barycentric_weight = np.zeros([h, w, 3], dtype = np.float32) # - for i in range(triangles.shape[0]): - tri = triangles[i, :] # 3 vertex indices - - # the inner bounding box - umin = max(int(np.ceil(np.min(vertices[tri, 0]))), 0) - umax = min(int(np.floor(np.max(vertices[tri, 0]))), w-1) - - vmin = max(int(np.ceil(np.min(vertices[tri, 1]))), 0) - vmax = min(int(np.floor(np.max(vertices[tri, 1]))), h-1) - - if umax depth_buffer[v, u]: - depth_buffer[v, u] = point_depth - triangle_buffer[v, u] = i - barycentric_weight[v, u, :] = np.array([w0, w1, w2]) - - return depth_buffer, triangle_buffer, barycentric_weight + vertices = vertices.astype(np.float32).copy() + triangles = triangles.astype(np.int32).copy() + mesh_core_cython.rasterize_triangles_core( + vertices, triangles, + depth_buffer, triangle_buffer, barycentric_weight, + vertices.shape[0], triangles.shape[0], + h, w) -def render_colors_ras(vertices, triangles, colors, h, w, c = 3): - ''' render mesh with colors(rasterize triangle first) +def render_colors(vertices, triangles, colors, h, w, c = 3, BG = None): + ''' render mesh with colors Args: vertices: [nver, 3] triangles: [ntri, 3] colors: [nver, 3] h: height - w: width + w: width c: channel + BG: background image Returns: - image: [h, w, c]. rendering. + image: [h, w, c]. rendered image./rendering. ''' - assert vertices.shape[0] == colors.shape[0] - - depth_buffer, triangle_buffer, barycentric_weight = rasterize_triangles(vertices, triangles, h, w) - - triangle_buffer_flat = np.reshape(triangle_buffer, [-1]) # [h*w] - barycentric_weight_flat = np.reshape(barycentric_weight, [-1, c]) #[h*w, c] - weight = barycentric_weight_flat[:, :, np.newaxis] # [h*w, 3(ver in tri), 1] - - colors_flat = colors[triangles[triangle_buffer_flat, :], :] # [h*w(tri id in pixel), 3(ver in tri), c(color in ver)] - colors_flat = weight*colors_flat # [h*w, 3, 3] - colors_flat = np.sum(colors_flat, 1) #[h*w, 3]. add tri. - - image = np.reshape(colors_flat, [h, w, c]) - # mask = (triangle_buffer[:,:] > -1).astype(np.float32) - # image = image*mask[:,:,np.newaxis] - return image - -def render_colors(vertices, triangles, colors, h, w, c = 3): - ''' render mesh with colors - Args: - vertices: [nver, 3] - triangles: [ntri, 3] - colors: [nver, 3] - h: height - w: width - Returns: - image: [h, w, c]. - ''' - assert vertices.shape[0] == colors.shape[0] - # initial - image = np.zeros((h, w, c)) - depth_buffer = np.zeros([h, w]) - 999999. - - for i in range(triangles.shape[0]): - tri = triangles[i, :] # 3 vertex indices - - # the inner bounding box - umin = max(int(np.ceil(np.min(vertices[tri, 0]))), 0) - umax = min(int(np.floor(np.max(vertices[tri, 0]))), w-1) - - vmin = max(int(np.ceil(np.min(vertices[tri, 1]))), 0) - vmax = min(int(np.floor(np.max(vertices[tri, 1]))), h-1) - - if umax depth_buffer[v, u]: - depth_buffer[v, u] = point_depth - image[v, u, :] = w0*colors[tri[0], :] + w1*colors[tri[1], :] + w2*colors[tri[2], :] + if BG is None: + image = np.zeros((h, w, c), dtype = np.float32) + else: + assert BG.shape[0] == h and BG.shape[1] == w and BG.shape[2] == c + image = BG + depth_buffer = np.zeros([h, w], dtype = np.float32, order = 'C') - 999999. + + # change orders. --> C-contiguous order(column major) + vertices = vertices.astype(np.float32).copy() + triangles = triangles.astype(np.int32).copy() + colors = colors.astype(np.float32).copy() + ### + st = time() + mesh_core_cython.render_colors_core( + image, vertices, triangles, + colors, + depth_buffer, + vertices.shape[0], triangles.shape[0], + h, w, c) return image -def render_texture(vertices, triangles, texture, tex_coords, tex_triangles, h, w, c = 3, mapping_type = 'nearest'): +def render_texture(vertices, triangles, texture, tex_coords, tex_triangles, h, w, c = 3, mapping_type = 'nearest', BG = None): ''' render mesh with texture map Args: - vertices: [nver], 3 - triangles: [ntri, 3] + vertices: [3, nver] + triangles: [3, ntri] texture: [tex_h, tex_w, 3] tex_coords: [ntexcoords, 3] tex_triangles: [ntri, 3] @@ -227,58 +99,37 @@ def render_texture(vertices, triangles, texture, tex_coords, tex_triangles, h, w c: channel mapping_type: 'bilinear' or 'nearest' ''' - assert triangles.shape[0] == tex_triangles.shape[0] - tex_h, tex_w, _ = texture.shape - # initial - image = np.zeros((h, w, c)) - depth_buffer = np.zeros([h, w]) - 999999. - - for i in range(triangles.shape[0]): - tri = triangles[i, :] # 3 vertex indices - tex_tri = tex_triangles[i, :] # 3 tex indice - - # the inner bounding box - umin = max(int(np.ceil(np.min(vertices[tri, 0]))), 0) - umax = min(int(np.floor(np.max(vertices[tri, 0]))), w-1) - - vmin = max(int(np.ceil(np.min(vertices[tri, 1]))), 0) - vmax = min(int(np.floor(np.max(vertices[tri, 1]))), h-1) - - if umax depth_buffer[v, u]: - # update depth - depth_buffer[v, u] = point_depth - - # tex coord - tex_xy = w0*tex_coords[tex_tri[0], :] + w1*tex_coords[tex_tri[1], :] + w2*tex_coords[tex_tri[2], :] - tex_xy[0] = max(min(tex_xy[0], float(tex_w - 1)), 0.0); - tex_xy[1] = max(min(tex_xy[1], float(tex_h - 1)), 0.0); - - # nearest - if mapping_type == 'nearest': - tex_xy = np.round(tex_xy).astype(np.int32) - tex_value = texture[tex_xy[1], tex_xy[0], :] - - # bilinear - elif mapping_type == 'bilinear': - # next 4 pixels - ul = texture[int(np.floor(tex_xy[1])), int(np.floor(tex_xy[0])), :] - ur = texture[int(np.floor(tex_xy[1])), int(np.ceil(tex_xy[0])), :] - dl = texture[int(np.ceil(tex_xy[1])), int(np.floor(tex_xy[0])), :] - dr = texture[int(np.ceil(tex_xy[1])), int(np.ceil(tex_xy[0])), :] + if BG is None: + image = np.zeros((h, w, c), dtype = np.float32) + else: + assert BG.shape[0] == h and BG.shape[1] == w and BG.shape[2] == c + image = BG - yd = tex_xy[1] - np.floor(tex_xy[1]) - xd = tex_xy[0] - np.floor(tex_xy[0]) - tex_value = ul*(1-xd)*(1-yd) + ur*xd*(1-yd) + dl*(1-xd)*yd + dr*xd*yd + depth_buffer = np.zeros([h, w], dtype = np.float32, order = 'C') - 999999. + + tex_h, tex_w, tex_c = texture.shape + if mapping_type == 'nearest': + mt = int(0) + elif mapping_type == 'bilinear': + mt = int(1) + else: + mt = int(0) + + # -> C order + vertices = vertices.astype(np.float32).copy() + triangles = triangles.astype(np.int32).copy() + texture = texture.astype(np.float32).copy() + tex_coords = tex_coords.astype(np.float32).copy() + tex_triangles = tex_triangles.astype(np.int32).copy() + + mesh_core_cython.render_texture_core( + image, vertices, triangles, + texture, tex_coords, tex_triangles, + depth_buffer, + vertices.shape[0], tex_coords.shape[0], triangles.shape[0], + h, w, c, + tex_h, tex_w, tex_c, + mt) + return image - image[v, u, :] = tex_value - return image \ No newline at end of file diff --git a/face3d/mesh/transform.py b/face3d/mesh/transform.py index 5eda05b..d91b09f 100644 --- a/face3d/mesh/transform.py +++ b/face3d/mesh/transform.py @@ -2,6 +2,9 @@ Functions about transforming mesh(changing the position: modify vertices). 1. forward: transform(transform, camera, project). 2. backward: estimate transform matrix from correspondences. + +Author: Yao Feng +Mail: yaofeng1995@gmail.com ''' from __future__ import absolute_import @@ -200,6 +203,12 @@ def to_image(vertices, h, w, is_perspective = False): ''' change vertices to image coord system 3d system: XYZ, center(0, 0, 0) 2d image: x(u), y(v). center(w/2, h/2), flip y-axis. + Args: + vertices: [nver, 3] + h: height of the rendering + w : width of the rendering + Returns: + projected_vertices: [nver, 3] ''' image_vertices = vertices.copy() if is_perspective: diff --git a/face3d/mesh_cython/__init__.py b/face3d/mesh_cython/__init__.py deleted file mode 100644 index 349121c..0000000 --- a/face3d/mesh_cython/__init__.py +++ /dev/null @@ -1,4 +0,0 @@ -from . import mesh_core_cython -from . import light -from . import render - diff --git a/face3d/mesh_cython/mesh_core.cpp b/face3d/mesh_cython/mesh_core.cpp deleted file mode 100644 index 7ea88af..0000000 --- a/face3d/mesh_cython/mesh_core.cpp +++ /dev/null @@ -1,328 +0,0 @@ -#include "mesh_core.h" - - -/* Judge whether the point is in the triangle -Method: - http://blackpawn.com/texts/pointinpoly/ - https://blogs.msdn.microsoft.com/rezanour/2011/08/07/barycentric-coordinates-and-point-in-triangle-tests/ -Args: - point: [x, y] - tri_points: three vertices(2d points) of a triangle. 2 coords x 3 vertices -Returns: - bool: true for in triangle -*/ -bool isPointInTri(point p, point p0, point p1, point p2) -{ - // vectors - point v0, v1, v2; - v0 = p2 - p0; - v1 = p1 - p0; - v2 = p - p0; - - // dot products - float dot00 = v0.dot(v0); //v0.x * v0.x + v0.y * v0.y //np.dot(v0.T, v0) - float dot01 = v0.dot(v1); //v0.x * v1.x + v0.y * v1.y //np.dot(v0.T, v1) - float dot02 = v0.dot(v2); //v0.x * v2.x + v0.y * v2.y //np.dot(v0.T, v2) - float dot11 = v1.dot(v1); //v1.x * v1.x + v1.y * v1.y //np.dot(v1.T, v1) - float dot12 = v1.dot(v2); //v1.x * v2.x + v1.y * v2.y//np.dot(v1.T, v2) - - // barycentric coordinates - float inverDeno; - if(dot00*dot11 - dot01*dot01 == 0) - inverDeno = 0; - else - inverDeno = 1/(dot00*dot11 - dot01*dot01); - - float u = (dot11*dot02 - dot01*dot12)*inverDeno; - float v = (dot00*dot12 - dot01*dot02)*inverDeno; - - // check if point in triangle - return (u >= 0) && (v >= 0) && (u + v < 1); -} - - -void get_point_weight(float* weight, point p, point p0, point p1, point p2) -{ - // vectors - point v0, v1, v2; - v0 = p2 - p0; - v1 = p1 - p0; - v2 = p - p0; - - // dot products - float dot00 = v0.dot(v0); //v0.x * v0.x + v0.y * v0.y //np.dot(v0.T, v0) - float dot01 = v0.dot(v1); //v0.x * v1.x + v0.y * v1.y //np.dot(v0.T, v1) - float dot02 = v0.dot(v2); //v0.x * v2.x + v0.y * v2.y //np.dot(v0.T, v2) - float dot11 = v1.dot(v1); //v1.x * v1.x + v1.y * v1.y //np.dot(v1.T, v1) - float dot12 = v1.dot(v2); //v1.x * v2.x + v1.y * v2.y//np.dot(v1.T, v2) - - // barycentric coordinates - float inverDeno; - if(dot00*dot11 - dot01*dot01 == 0) - inverDeno = 0; - else - inverDeno = 1/(dot00*dot11 - dot01*dot01); - - float u = (dot11*dot02 - dot01*dot12)*inverDeno; - float v = (dot00*dot12 - dot01*dot02)*inverDeno; - - // weight - weight[0] = 1 - u - v; - weight[1] = v; - weight[2] = u; -} - - -void _get_normal_core( - float* normal, float* tri_normal, int* triangles, - int ntri) -{ - int i, j; - int tri_p0_ind, tri_p1_ind, tri_p2_ind; - - for(i = 0; i < ntri; i++) - { - tri_p0_ind = triangles[3*i]; - tri_p1_ind = triangles[3*i + 1]; - tri_p2_ind = triangles[3*i + 2]; - - for(j = 0; j < 3; j++) - { - normal[3*tri_p0_ind + j] = normal[3*tri_p0_ind + j] + tri_normal[3*i + j]; - normal[3*tri_p1_ind + j] = normal[3*tri_p1_ind + j] + tri_normal[3*i + j]; - normal[3*tri_p2_ind + j] = normal[3*tri_p2_ind + j] + tri_normal[3*i + j]; - } - } -} - - -void _rasterize_triangles_core( - float* vertices, int* triangles, - float* depth_buffer, int* triangle_buffer, float* barycentric_weight, - int nver, int ntri, - int h, int w) -{ - int i; - int x, y, k; - int tri_p0_ind, tri_p1_ind, tri_p2_ind; - point p0, p1, p2, p; - int x_min, x_max, y_min, y_max; - float p_depth, p0_depth, p1_depth, p2_depth; - float weight[3]; - - for(i = 0; i < ntri; i++) - { - tri_p0_ind = triangles[3*i]; - tri_p1_ind = triangles[3*i + 1]; - tri_p2_ind = triangles[3*i + 2]; - - p0.x = vertices[3*tri_p0_ind]; p0.y = vertices[3*tri_p0_ind + 1]; p0_depth = vertices[3*tri_p0_ind + 2]; - p1.x = vertices[3*tri_p1_ind]; p1.y = vertices[3*tri_p1_ind + 1]; p1_depth = vertices[3*tri_p1_ind + 2]; - p2.x = vertices[3*tri_p2_ind]; p2.y = vertices[3*tri_p2_ind + 1]; p2_depth = vertices[3*tri_p2_ind + 2]; - - x_min = max((int)ceil(min(p0.x, min(p1.x, p2.x))), 0); - x_max = min((int)floor(max(p0.x, max(p1.x, p2.x))), w - 1); - - y_min = max((int)ceil(min(p0.y, min(p1.y, p2.y))), 0); - y_max = min((int)floor(max(p0.y, max(p1.y, p2.y))), h - 1); - - if(x_max < x_min || y_max < y_min) - { - continue; - } - - for(y = y_min; y <= y_max; y++) //h - { - for(x = x_min; x <= x_max; x++) //w - { - p.x = x; p.y = y; - if(p.x < 2 || p.x > w - 3 || p.y < 2 || p.y > h - 3 || isPointInTri(p, p0, p1, p2)) - { - get_point_weight(weight, p, p0, p1, p2); - p_depth = weight[0]*p0_depth + weight[1]*p1_depth + weight[2]*p2_depth; - - if((p_depth > depth_buffer[y*w + x])) - { - depth_buffer[y*w + x] = p_depth; - triangle_buffer[y*w + x] = i; - for(k = 0; k < 3; k++) - { - barycentric_weight[y*w*3 + x*3 + k] = weight[k]; - } - } - } - } - } - } -} - - -void _render_colors_core( - float* image, float* vertices, int* triangles, - float* colors, - float* depth_buffer, - int nver, int ntri, - int h, int w, int c) -{ - int i; - int x, y, k; - int tri_p0_ind, tri_p1_ind, tri_p2_ind; - point p0, p1, p2, p; - int x_min, x_max, y_min, y_max; - float p_depth, p0_depth, p1_depth, p2_depth; - float p_color, p0_color, p1_color, p2_color; - float weight[3]; - - for(i = 0; i < ntri; i++) - { - tri_p0_ind = triangles[3*i]; - tri_p1_ind = triangles[3*i + 1]; - tri_p2_ind = triangles[3*i + 2]; - - p0.x = vertices[3*tri_p0_ind]; p0.y = vertices[3*tri_p0_ind + 1]; p0_depth = vertices[3*tri_p0_ind + 2]; - p1.x = vertices[3*tri_p1_ind]; p1.y = vertices[3*tri_p1_ind + 1]; p1_depth = vertices[3*tri_p1_ind + 2]; - p2.x = vertices[3*tri_p2_ind]; p2.y = vertices[3*tri_p2_ind + 1]; p2_depth = vertices[3*tri_p2_ind + 2]; - - x_min = max((int)ceil(min(p0.x, min(p1.x, p2.x))), 0); - x_max = min((int)floor(max(p0.x, max(p1.x, p2.x))), w - 1); - - y_min = max((int)ceil(min(p0.y, min(p1.y, p2.y))), 0); - y_max = min((int)floor(max(p0.y, max(p1.y, p2.y))), h - 1); - - if(x_max < x_min || y_max < y_min) - { - continue; - } - - for(y = y_min; y <= y_max; y++) //h - { - for(x = x_min; x <= x_max; x++) //w - { - p.x = x; p.y = y; - if(p.x < 2 || p.x > w - 3 || p.y < 2 || p.y > h - 3 || isPointInTri(p, p0, p1, p2)) - { - get_point_weight(weight, p, p0, p1, p2); - p_depth = weight[0]*p0_depth + weight[1]*p1_depth + weight[2]*p2_depth; - - if((p_depth > depth_buffer[y*w + x])) - { - for(k = 0; k < c; k++) // c - { - p0_color = colors[c*tri_p0_ind + k]; - p1_color = colors[c*tri_p1_ind + k]; - p2_color = colors[c*tri_p2_ind + k]; - - p_color = weight[0]*p0_color + weight[1]*p1_color + weight[2]*p2_color; - image[y*w*c + x*c + k] = p_color; - } - - depth_buffer[y*w + x] = p_depth; - } - } - } - } - } -} - - -void _render_texture_core( - float* image, float* vertices, int* triangles, - float* texture, float* tex_coords, int* tex_triangles, - float* depth_buffer, - int nver, int tex_nver, int ntri, - int h, int w, int c, - int tex_h, int tex_w, int tex_c, - int mapping_type) -{ - int i; - int x, y, k; - int tri_p0_ind, tri_p1_ind, tri_p2_ind; - int tex_tri_p0_ind, tex_tri_p1_ind, tex_tri_p2_ind; - point p0, p1, p2, p; - point tex_p0, tex_p1, tex_p2, tex_p; - int x_min, x_max, y_min, y_max; - float weight[3]; - float p_depth, p0_depth, p1_depth, p2_depth; - float xd, yd; - float ul, ur, dl, dr; - for(i = 0; i < ntri; i++) - { - // mesh - tri_p0_ind = triangles[3*i]; - tri_p1_ind = triangles[3*i + 1]; - tri_p2_ind = triangles[3*i + 2]; - - p0.x = vertices[3*tri_p0_ind]; p0.y = vertices[3*tri_p0_ind + 1]; p0_depth = vertices[3*tri_p0_ind + 2]; - p1.x = vertices[3*tri_p1_ind]; p1.y = vertices[3*tri_p1_ind + 1]; p1_depth = vertices[3*tri_p1_ind + 2]; - p2.x = vertices[3*tri_p2_ind]; p2.y = vertices[3*tri_p2_ind + 1]; p2_depth = vertices[3*tri_p2_ind + 2]; - - // texture - tex_tri_p0_ind = tex_triangles[3*i]; - tex_tri_p1_ind = tex_triangles[3*i + 1]; - tex_tri_p2_ind = tex_triangles[3*i + 2]; - - tex_p0.x = tex_coords[3*tex_tri_p0_ind]; tex_p0.y = tex_coords[3*tri_p0_ind + 1]; - tex_p1.x = tex_coords[3*tex_tri_p1_ind]; tex_p1.y = tex_coords[3*tri_p1_ind + 1]; - tex_p2.x = tex_coords[3*tex_tri_p2_ind]; tex_p2.y = tex_coords[3*tri_p2_ind + 1]; - - - x_min = max((int)ceil(min(p0.x, min(p1.x, p2.x))), 0); - x_max = min((int)floor(max(p0.x, max(p1.x, p2.x))), w - 1); - - y_min = max((int)ceil(min(p0.y, min(p1.y, p2.y))), 0); - y_max = min((int)floor(max(p0.y, max(p1.y, p2.y))), h - 1); - - - if(x_max < x_min || y_max < y_min) - { - continue; - } - - for(y = y_min; y <= y_max; y++) //h - { - for(x = x_min; x <= x_max; x++) //w - { - p.x = x; p.y = y; - if(p.x < 2 || p.x > w - 3 || p.y < 2 || p.y > h - 3 || isPointInTri(p, p0, p1, p2)) - { - get_point_weight(weight, p, p0, p1, p2); - p_depth = weight[0]*p0_depth + weight[1]*p1_depth + weight[2]*p2_depth; - - if((p_depth > depth_buffer[y*w + x])) - { - // -- color from texture - // cal weight in mesh tri - get_point_weight(weight, p, p0, p1, p2); - // cal coord in texture - tex_p = tex_p0*weight[0] + tex_p1*weight[1] + tex_p2*weight[2]; - tex_p.x = max(min(tex_p.x, float(tex_w - 1)), float(0)); - tex_p.y = max(min(tex_p.y, float(tex_h - 1)), float(0)); - - yd = tex_p.y - floor(tex_p.y); - xd = tex_p.x - floor(tex_p.x); - for(k = 0; k < c; k++) - { - if(mapping_type==0)// nearest - { - image[y*w*c + x*c + k] = texture[int(round(tex_p.y))*tex_w*tex_c + int(round(tex_p.x))*tex_c + k]; - } - else//bilinear interp - { - ul = texture[(int)floor(tex_p.y)*tex_w*tex_c + (int)floor(tex_p.x)*tex_c + k]; - ur = texture[(int)floor(tex_p.y)*tex_w*tex_c + (int)ceil(tex_p.x)*tex_c + k]; - dl = texture[(int)ceil(tex_p.y)*tex_w*tex_c + (int)floor(tex_p.x)*tex_c + k]; - dr = texture[(int)ceil(tex_p.y)*tex_w*tex_c + (int)ceil(tex_p.x)*tex_c + k]; - - image[y*w*c + x*c + k] = ul*(1-xd)*(1-yd) + ur*xd*(1-yd) + dl*(1-xd)*yd + dr*xd*yd; - } - - } - - depth_buffer[y*w + x] = p_depth; - } - } - } - } - } -} - diff --git a/face3d/mesh_cython/mesh_core.h b/face3d/mesh_cython/mesh_core.h deleted file mode 100644 index 24ec46d..0000000 --- a/face3d/mesh_cython/mesh_core.h +++ /dev/null @@ -1,76 +0,0 @@ -#ifndef MESH_CORE_HPP_ -#define MESH_CORE_HPP_ - -#include -#include -#include - -using namespace std; - -class point -{ - public: - float x; - float y; - - float dot(point p) - { - return this->x * p.x + this->y * p.y; - } - - point operator-(const point& p) - { - point np; - np.x = this->x - p.x; - np.y = this->y - p.y; - return np; - } - - point operator+(const point& p) - { - point np; - np.x = this->x + p.x; - np.y = this->y + p.y; - return np; - } - - point operator*(float s) - { - point np; - np.x = s * this->x; - np.y = s * this->y; - return np; - } -}; - - -bool isPointInTri(point p, point p0, point p1, point p2, int h, int w); -void get_point_weight(float* weight, point p, point p0, point p1, point p2); - -void _get_normal_core( - float* normal, float* tri_normal, int* triangles, - int ntri); - -void _rasterize_triangles_core( - float* vertices, int* triangles, - float* depth_buffer, int* triangle_buffer, float* barycentric_weight, - int nver, int ntri, - int h, int w); - -void _render_colors_core( - float* image, float* vertices, int* triangles, - float* colors, - float* depth_buffer, - int nver, int ntri, - int h, int w, int c); - -void _render_texture_core( - float* image, float* vertices, int* triangles, - float* texture, float* tex_coords, int* tex_triangles, - float* depth_buffer, - int nver, int tex_nver, int ntri, - int h, int w, int c, - int tex_h, int tex_w, int tex_c, - int mapping_type); - -#endif \ No newline at end of file diff --git a/face3d/mesh_cython/mesh_core_cython.pyx b/face3d/mesh_cython/mesh_core_cython.pyx deleted file mode 100644 index d3f90ca..0000000 --- a/face3d/mesh_cython/mesh_core_cython.pyx +++ /dev/null @@ -1,93 +0,0 @@ -import numpy as np -cimport numpy as np - -# use the Numpy-C-API from Cython -np.import_array() - -# cdefine the signature of our c function -cdef extern from "mesh_core.h": - void _rasterize_triangles_core( - float* vertices, int* triangles, - float* depth_buffer, int* triangle_buffer, float* barycentric_weight, - int nver, int ntri, - int h, int w) - - void _render_colors_core( - float* image, float* vertices, int* triangles, - float* colors, - float* depth_buffer, - int nver, int ntri, - int h, int w, int c) - - void _render_texture_core( - float* image, float* vertices, int* triangles, - float* texture, float* tex_coords, int* tex_triangles, - float* depth_buffer, - int nver, int tex_nver, int ntri, - int h, int w, int c, - int tex_h, int tex_w, int tex_c, - int mapping_type) - - void _get_normal_core( - float* normal, float* tri_normal, int* triangles, - int ntri) - -def get_normal_core(np.ndarray[float, ndim=2, mode = "c"] normal not None, - np.ndarray[float, ndim=2, mode = "c"] tri_normal not None, - np.ndarray[int, ndim=2, mode="c"] triangles not None, - int ntri - ): - _get_normal_core( - np.PyArray_DATA(normal), np.PyArray_DATA(tri_normal), np.PyArray_DATA(triangles), - ntri) - -def rasterize_triangles_core( - np.ndarray[float, ndim=2, mode = "c"] vertices not None, - np.ndarray[int, ndim=2, mode="c"] triangles not None, - np.ndarray[float, ndim=2, mode = "c"] depth_buffer not None, - np.ndarray[int, ndim=2, mode = "c"] triangle_buffer not None, - np.ndarray[float, ndim=2, mode = "c"] barycentric_weight not None, - int nver, int ntri, - int h, int w - ): - _rasterize_triangles_core( - np.PyArray_DATA(vertices), np.PyArray_DATA(triangles), - np.PyArray_DATA(depth_buffer), np.PyArray_DATA(triangle_buffer), np.PyArray_DATA(barycentric_weight), - nver, ntri, - h, w) - -def render_colors_core(np.ndarray[float, ndim=3, mode = "c"] image not None, - np.ndarray[float, ndim=2, mode = "c"] vertices not None, - np.ndarray[int, ndim=2, mode="c"] triangles not None, - np.ndarray[float, ndim=2, mode = "c"] colors not None, - np.ndarray[float, ndim=2, mode = "c"] depth_buffer not None, - int nver, int ntri, - int h, int w, int c - ): - _render_colors_core( - np.PyArray_DATA(image), np.PyArray_DATA(vertices), np.PyArray_DATA(triangles), - np.PyArray_DATA(colors), - np.PyArray_DATA(depth_buffer), - nver, ntri, - h, w, c) - -def render_texture_core(np.ndarray[float, ndim=3, mode = "c"] image not None, - np.ndarray[float, ndim=2, mode = "c"] vertices not None, - np.ndarray[int, ndim=2, mode="c"] triangles not None, - np.ndarray[float, ndim=3, mode = "c"] texture not None, - np.ndarray[float, ndim=2, mode = "c"] tex_coords not None, - np.ndarray[int, ndim=2, mode="c"] tex_triangles not None, - np.ndarray[float, ndim=2, mode = "c"] depth_buffer not None, - int nver, int tex_nver, int ntri, - int h, int w, int c, - int tex_h, int tex_w, int tex_c, - int mapping_type - ): - _render_texture_core( - np.PyArray_DATA(image), np.PyArray_DATA(vertices), np.PyArray_DATA(triangles), - np.PyArray_DATA(texture), np.PyArray_DATA(tex_coords), np.PyArray_DATA(tex_triangles), - np.PyArray_DATA(depth_buffer), - nver, tex_nver, ntri, - h, w, c, - tex_h, tex_w, tex_c, - mapping_type) diff --git a/face3d/mesh_cython/render.py b/face3d/mesh_cython/render.py deleted file mode 100644 index 31fe0e4..0000000 --- a/face3d/mesh_cython/render.py +++ /dev/null @@ -1,131 +0,0 @@ -''' -functions about rendering mesh(from 3d obj to 2d image). -only use rasterization render here. -Note that: -1. Generally, render func includes camera, light, raterize. Here no camera and light(I write these in other files) -2. Generally, the input vertices are normalized to [-1,1] and cetered on [0, 0]. (in world space) - Here, the vertices are using image coords, which centers on [w/2, h/2] with the y-axis pointing to oppisite direction. -Means: render here only conducts interpolation.(I just want to make the input flexible) - -Author: YadiraF -Mail: fengyao@sjtu.edu.cn -''' -import numpy as np -from time import time - -import mesh_core_cython - -def rasterize_triangles(vertices, triangles, h, w): - ''' - Args: - vertices: [nver, 3] - triangles: [ntri, 3] - h: height - w: width - Returns: - depth_buffer: [h, w] saves the depth, here, the bigger the z, the fronter the point. - triangle_buffer: [h, w] saves the tri id(-1 for no triangle). - barycentric_weight: [h, w, 3] saves corresponding barycentric weight. - - # Each triangle has 3 vertices & Each vertex has 3 coordinates x, y, z. - # h, w is the size of rendering - ''' - - # initial - depth_buffer = np.zeros([h, w]) - 999999. #set the initial z to the farest position - triangle_buffer = np.zeros([h, w], dtype = np.int32) - 1 # if tri id = -1, the pixel has no triangle correspondance - barycentric_weight = np.zeros([h, w, 3], dtype = np.float32) # - - vertices = vertices.astype(np.float32).copy() - triangles = triangles.astype(np.int32).copy() - - mesh_core_cython.rasterize_triangles_core( - vertices, triangles, - depth_buffer, triangle_buffer, barycentric_weight, - vertices.shape[0], triangles.shape[0], - h, w) - -def render_colors(vertices, triangles, colors, h, w, c = 3, BG = None): - ''' render mesh with colors - Args: - vertices: [nver, 3] - triangles: [ntri, 3] - colors: [nver, 3] - h: height - w: width - c: channel - BG: background image - Returns: - image: [h, w, c]. rendered image./rendering. - ''' - - # initial - if BG is None: - image = np.zeros((h, w, c), dtype = np.float32) - else: - assert BG.shape[0] == h and BG.shape[1] == w and BG.shape[2] == c - image = BG - depth_buffer = np.zeros([h, w], dtype = np.float32, order = 'C') - 999999. - - # change orders. --> C-contiguous order(column major) - vertices = vertices.astype(np.float32).copy() - triangles = triangles.astype(np.int32).copy() - colors = colors.astype(np.float32).copy() - ### - st = time() - mesh_core_cython.render_colors_core( - image, vertices, triangles, - colors, - depth_buffer, - vertices.shape[0], triangles.shape[0], - h, w, c) - return image - - -def render_texture(vertices, triangles, texture, tex_coords, tex_triangles, h, w, c = 3, mapping_type = 'nearest', BG = None): - ''' render mesh with texture map - Args: - vertices: [3, nver] - triangles: [3, ntri] - texture: [tex_h, tex_w, 3] - tex_coords: [ntexcoords, 3] - tex_triangles: [ntri, 3] - h: height of rendering - w: width of rendering - c: channel - mapping_type: 'bilinear' or 'nearest' - ''' - # initial - if BG is None: - image = np.zeros((h, w, c), dtype = np.float32) - else: - assert BG.shape[0] == h and BG.shape[1] == w and BG.shape[2] == c - image = BG - - depth_buffer = np.zeros([h, w], dtype = np.float32, order = 'C') - 999999. - - tex_h, tex_w, tex_c = texture.shape - if mapping_type == 'nearest': - mt = int(0) - elif mapping_type == 'bilinear': - mt = int(1) - else: - mt = int(0) - - # -> C order - vertices = vertices.astype(np.float32).copy() - triangles = triangles.astype(np.int32).copy() - texture = texture.astype(np.float32).copy() - tex_coords = tex_coords.astype(np.float32).copy() - tex_triangles = tex_triangles.astype(np.int32).copy() - - mesh_core_cython.render_texture_core( - image, vertices, triangles, - texture, tex_coords, tex_triangles, - depth_buffer, - vertices.shape[0], tex_coords.shape[0], triangles.shape[0], - h, w, c, - tex_h, tex_w, tex_c, - mt) - return image - diff --git a/face3d/mesh_cython/setup.py b/face3d/mesh_cython/setup.py deleted file mode 100644 index c20dab9..0000000 --- a/face3d/mesh_cython/setup.py +++ /dev/null @@ -1,30 +0,0 @@ -''' -python setup.py build_ext -i -to compile -''' - -# setup.py -from distutils.core import setup, Extension -from Cython.Build import cythonize -from Cython.Distutils import build_ext -import numpy -# # setup(ext_modules = cythonize(Extension( -# # 'render_texture', -# # sources=['render_texture.pyx'], -# # language='c++', -# # include_dirs=[numpy.get_include()], -# # library_dirs=[], -# # libraries=[], -# # extra_compile_args=[], -# # extra_link_args=[] -# # ))) - -setup( - name = 'mesh_core_cython', - cmdclass={'build_ext': build_ext}, - ext_modules=[Extension("mesh_core_cython", - sources=["mesh_core_cython.pyx", "mesh_core.cpp"], - language='c++', - include_dirs=[numpy.get_include()])], -) - diff --git a/face3d/mesh_numpy/__init__.py b/face3d/mesh_numpy/__init__.py new file mode 100644 index 0000000..65d5039 --- /dev/null +++ b/face3d/mesh_numpy/__init__.py @@ -0,0 +1,10 @@ +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function + +from . import io +from . import vis +from . import transform +from . import light +from . import render + diff --git a/face3d/mesh_numpy/io.py b/face3d/mesh_numpy/io.py new file mode 100644 index 0000000..0e48a76 --- /dev/null +++ b/face3d/mesh_numpy/io.py @@ -0,0 +1,170 @@ +''' io: read&write mesh +1. read obj as array(TODO) +2. write arrays to obj + +Preparation knowledge: +representations of 3d face: mesh, point cloud... +storage format: obj, ply, bin, asc, mat... +''' + +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function + +import numpy as np +import os +from skimage import io + +## TODO +## TODO: c++ version +def read_obj(obj_name): + ''' read mesh + ''' + return 0 + +# ------------------------- write +def write_asc(path, vertices): + ''' + Args: + vertices: shape = (nver, 3) + ''' + if path.split('.')[-1] == 'asc': + np.savetxt(path, vertices) + else: + np.savetxt(path + '.asc', vertices) + +def write_obj_with_colors(obj_name, vertices, triangles, colors): + ''' Save 3D face model with texture represented by colors. + Args: + obj_name: str + vertices: shape = (nver, 3) + triangles: shape = (ntri, 3) + colors: shape = (nver, 3) + ''' + triangles = triangles.copy() + triangles += 1 # meshlab start with 1 + + if obj_name.split('.')[-1] != 'obj': + obj_name = obj_name + '.obj' + + # write obj + with open(obj_name, 'w') as f: + + # write vertices & colors + for i in range(vertices.shape[0]): + # s = 'v {} {} {} \n'.format(vertices[0,i], vertices[1,i], vertices[2,i]) + s = 'v {} {} {} {} {} {}\n'.format(vertices[i, 0], vertices[i, 1], vertices[i, 2], colors[i, 0], colors[i, 1], colors[i, 2]) + f.write(s) + + # write f: ver ind/ uv ind + [k, ntri] = triangles.shape + for i in range(triangles.shape[0]): + # s = 'f {} {} {}\n'.format(triangles[i, 0], triangles[i, 1], triangles[i, 2]) + s = 'f {} {} {}\n'.format(triangles[i, 2], triangles[i, 1], triangles[i, 0]) + f.write(s) + +## TODO: c++ version +def write_obj_with_texture(obj_name, vertices, triangles, texture, uv_coords): + ''' Save 3D face model with texture represented by texture map. + Ref: https://github.com/patrikhuber/eos/blob/bd00155ebae4b1a13b08bf5a991694d682abbada/include/eos/core/Mesh.hpp + Args: + obj_name: str + vertices: shape = (nver, 3) + triangles: shape = (ntri, 3) + texture: shape = (256,256,3) + uv_coords: shape = (nver, 3) max value<=1 + ''' + if obj_name.split('.')[-1] != 'obj': + obj_name = obj_name + '.obj' + mtl_name = obj_name.replace('.obj', '.mtl') + texture_name = obj_name.replace('.obj', '_texture.png') + + triangles = triangles.copy() + triangles += 1 # mesh lab start with 1 + + # write obj + with open(obj_name, 'w') as f: + # first line: write mtlib(material library) + s = "mtllib {}\n".format(os.path.abspath(mtl_name)) + f.write(s) + + # write vertices + for i in range(vertices.shape[0]): + s = 'v {} {} {}\n'.format(vertices[i, 0], vertices[i, 1], vertices[i, 2]) + f.write(s) + + # write uv coords + for i in range(uv_coords.shape[0]): + # s = 'vt {} {}\n'.format(uv_coords[i,0], 1 - uv_coords[i,1]) + s = 'vt {} {}\n'.format(uv_coords[i,0], uv_coords[i,1]) + f.write(s) + + f.write("usemtl FaceTexture\n") + + # write f: ver ind/ uv ind + for i in range(triangles.shape[0]): + s = 'f {}/{} {}/{} {}/{}\n'.format(triangles[i,2], triangles[i,2], triangles[i,1], triangles[i,1], triangles[i,0], triangles[i,0]) + f.write(s) + + # write mtl + with open(mtl_name, 'w') as f: + f.write("newmtl FaceTexture\n") + s = 'map_Kd {}\n'.format(os.path.abspath(texture_name)) # map to image + f.write(s) + + # write texture as png + imsave(texture_name, texture) + + +def write_obj_with_colors_texture(obj_name, vertices, triangles, colors, texture, uv_coords): + ''' Save 3D face model with texture. + Ref: https://github.com/patrikhuber/eos/blob/bd00155ebae4b1a13b08bf5a991694d682abbada/include/eos/core/Mesh.hpp + Args: + obj_name: str + vertices: shape = (nver, 3) + triangles: shape = (ntri, 3) + colors: shape = (nver, 3) + texture: shape = (256,256,3) + uv_coords: shape = (nver, 3) max value<=1 + ''' + if obj_name.split('.')[-1] != 'obj': + obj_name = obj_name + '.obj' + mtl_name = obj_name.replace('.obj', '.mtl') + texture_name = obj_name.replace('.obj', '_texture.png') + + triangles = triangles.copy() + triangles += 1 # mesh lab start with 1 + + # write obj + with open(obj_name, 'w') as f: + # first line: write mtlib(material library) + s = "mtllib {}\n".format(os.path.abspath(mtl_name)) + f.write(s) + + # write vertices + for i in range(vertices.shape[0]): + s = 'v {} {} {} {} {} {}\n'.format(vertices[i, 0], vertices[i, 1], vertices[i, 2], colors[i, 0], colors[i, 1], colors[i, 2]) + f.write(s) + + # write uv coords + for i in range(uv_coords.shape[0]): + # s = 'vt {} {}\n'.format(uv_coords[i,0], 1 - uv_coords[i,1]) + s = 'vt {} {}\n'.format(uv_coords[i,0], uv_coords[i,1]) + f.write(s) + + f.write("usemtl FaceTexture\n") + + # write f: ver ind/ uv ind + for i in range(triangles.shape[0]): + # s = 'f {}/{} {}/{} {}/{}\n'.format(triangles[i,0], triangles[i,0], triangles[i,1], triangles[i,1], triangles[i,2], triangles[i,2]) + s = 'f {}/{} {}/{} {}/{}\n'.format(triangles[i,2], triangles[i,2], triangles[i,1], triangles[i,1], triangles[i,0], triangles[i,0]) + f.write(s) + + # write mtl + with open(mtl_name, 'w') as f: + f.write("newmtl FaceTexture\n") + s = 'map_Kd {}\n'.format(os.path.abspath(texture_name)) # map to image + f.write(s) + + # write texture as png + io.imsave(texture_name, texture) \ No newline at end of file diff --git a/face3d/mesh_cython/light.py b/face3d/mesh_numpy/light.py similarity index 92% rename from face3d/mesh_cython/light.py rename to face3d/mesh_numpy/light.py index 661c2be..bde5711 100644 --- a/face3d/mesh_cython/light.py +++ b/face3d/mesh_numpy/light.py @@ -2,6 +2,10 @@ Functions about lighting mesh(changing colors/texture of mesh). 1. add light to colors/texture (shade each vertex) 2. fit light according to colors/texture & image. + +Preparation knowledge: +lighting: https://cs184.eecs.berkeley.edu/lecture/pipeline +spherical harmonics in human face: '3D Face Reconstruction from a Single Image Using a Single Reference Face Shape' ''' from __future__ import absolute_import @@ -9,7 +13,6 @@ from __future__ import print_function import numpy as np -from . import mesh_core_cython def get_normal(vertices, triangles): ''' calculate normal direction in each vertex @@ -24,13 +27,12 @@ def get_normal(vertices, triangles): pt2 = vertices[triangles[:, 2], :] # [ntri, 3] tri_normal = np.cross(pt0 - pt1, pt0 - pt2) # [ntri, 3]. normal of each triangle - normal = np.zeros_like(vertices, dtype = np.float32).copy() # [nver, 3] - # for i in range(triangles.shape[0]): - # normal[triangles[i, 0], :] = normal[triangles[i, 0], :] + tri_normal[i, :] - # normal[triangles[i, 1], :] = normal[triangles[i, 1], :] + tri_normal[i, :] - # normal[triangles[i, 2], :] = normal[triangles[i, 2], :] + tri_normal[i, :] - mesh_core_cython.get_normal_core(normal, tri_normal.astype(np.float32).copy(), triangles.copy(), triangles.shape[0]) - + normal = np.zeros_like(vertices) # [nver, 3] + for i in range(triangles.shape[0]): + normal[triangles[i, 0], :] = normal[triangles[i, 0], :] + tri_normal[i, :] + normal[triangles[i, 1], :] = normal[triangles[i, 1], :] + tri_normal[i, :] + normal[triangles[i, 2], :] = normal[triangles[i, 2], :] + tri_normal[i, :] + # normalize to unit length mag = np.sum(normal**2, 1) # [nver] zero_ind = (mag == 0) diff --git a/face3d/mesh_numpy/render.py b/face3d/mesh_numpy/render.py new file mode 100644 index 0000000..eb0d92d --- /dev/null +++ b/face3d/mesh_numpy/render.py @@ -0,0 +1,287 @@ +''' +functions about rendering mesh(from 3d obj to 2d image). +only use rasterization render here. +Note that: +1. Generally, render func includes camera, light, raterize. Here no camera and light(I write these in other files) +2. Generally, the input vertices are normalized to [-1,1] and cetered on [0, 0]. (in world space) + Here, the vertices are using image coords, which centers on [w/2, h/2] with the y-axis pointing to oppisite direction. +Means: render here only conducts interpolation.(I just want to make the input flexible) + +Preparation knowledge: +z-buffer: https://cs184.eecs.berkeley.edu/lecture/pipeline + +Author: Yao Feng +Mail: yaofeng1995@gmail.com +''' +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function + +import numpy as np +from time import time + +def isPointInTri(point, tri_points): + ''' Judge whether the point is in the triangle + Method: + http://blackpawn.com/texts/pointinpoly/ + Args: + point: (2,). [u, v] or [x, y] + tri_points: (3 vertices, 2 coords). three vertices(2d points) of a triangle. + Returns: + bool: true for in triangle + ''' + tp = tri_points + + # vectors + v0 = tp[2,:] - tp[0,:] + v1 = tp[1,:] - tp[0,:] + v2 = point - tp[0,:] + + # dot products + dot00 = np.dot(v0.T, v0) + dot01 = np.dot(v0.T, v1) + dot02 = np.dot(v0.T, v2) + dot11 = np.dot(v1.T, v1) + dot12 = np.dot(v1.T, v2) + + # barycentric coordinates + if dot00*dot11 - dot01*dot01 == 0: + inverDeno = 0 + else: + inverDeno = 1/(dot00*dot11 - dot01*dot01) + + u = (dot11*dot02 - dot01*dot12)*inverDeno + v = (dot00*dot12 - dot01*dot02)*inverDeno + + # check if point in triangle + return (u >= 0) & (v >= 0) & (u + v < 1) + +def get_point_weight(point, tri_points): + ''' Get the weights of the position + Methods: https://gamedev.stackexchange.com/questions/23743/whats-the-most-efficient-way-to-find-barycentric-coordinates + -m1.compute the area of the triangles formed by embedding the point P inside the triangle + -m2.Christer Ericson's book "Real-Time Collision Detection". faster.(used) + Args: + point: (2,). [u, v] or [x, y] + tri_points: (3 vertices, 2 coords). three vertices(2d points) of a triangle. + Returns: + w0: weight of v0 + w1: weight of v1 + w2: weight of v3 + ''' + tp = tri_points + # vectors + v0 = tp[2,:] - tp[0,:] + v1 = tp[1,:] - tp[0,:] + v2 = point - tp[0,:] + + # dot products + dot00 = np.dot(v0.T, v0) + dot01 = np.dot(v0.T, v1) + dot02 = np.dot(v0.T, v2) + dot11 = np.dot(v1.T, v1) + dot12 = np.dot(v1.T, v2) + + # barycentric coordinates + if dot00*dot11 - dot01*dot01 == 0: + inverDeno = 0 + else: + inverDeno = 1/(dot00*dot11 - dot01*dot01) + + u = (dot11*dot02 - dot01*dot12)*inverDeno + v = (dot00*dot12 - dot01*dot02)*inverDeno + + w0 = 1 - u - v + w1 = v + w2 = u + + return w0, w1, w2 + +def rasterize_triangles(vertices, triangles, h, w): + ''' + Args: + vertices: [nver, 3] + triangles: [ntri, 3] + h: height + w: width + Returns: + depth_buffer: [h, w] saves the depth, here, the bigger the z, the fronter the point. + triangle_buffer: [h, w] saves the tri id(-1 for no triangle). + barycentric_weight: [h, w, 3] saves corresponding barycentric weight. + + # Each triangle has 3 vertices & Each vertex has 3 coordinates x, y, z. + # h, w is the size of rendering + ''' + # initial + depth_buffer = np.zeros([h, w]) - 999999. #+ np.min(vertices[2,:]) - 999999. # set the initial z to the farest position + triangle_buffer = np.zeros([h, w], dtype = np.int32) - 1 # if tri id = -1, the pixel has no triangle correspondance + barycentric_weight = np.zeros([h, w, 3], dtype = np.float32) # + + for i in range(triangles.shape[0]): + tri = triangles[i, :] # 3 vertex indices + + # the inner bounding box + umin = max(int(np.ceil(np.min(vertices[tri, 0]))), 0) + umax = min(int(np.floor(np.max(vertices[tri, 0]))), w-1) + + vmin = max(int(np.ceil(np.min(vertices[tri, 1]))), 0) + vmax = min(int(np.floor(np.max(vertices[tri, 1]))), h-1) + + if umax depth_buffer[v, u]: + depth_buffer[v, u] = point_depth + triangle_buffer[v, u] = i + barycentric_weight[v, u, :] = np.array([w0, w1, w2]) + + return depth_buffer, triangle_buffer, barycentric_weight + + +def render_colors_ras(vertices, triangles, colors, h, w, c = 3): + ''' render mesh with colors(rasterize triangle first) + Args: + vertices: [nver, 3] + triangles: [ntri, 3] + colors: [nver, 3] + h: height + w: width + c: channel + Returns: + image: [h, w, c]. rendering. + ''' + assert vertices.shape[0] == colors.shape[0] + + depth_buffer, triangle_buffer, barycentric_weight = rasterize_triangles(vertices, triangles, h, w) + + triangle_buffer_flat = np.reshape(triangle_buffer, [-1]) # [h*w] + barycentric_weight_flat = np.reshape(barycentric_weight, [-1, c]) #[h*w, c] + weight = barycentric_weight_flat[:, :, np.newaxis] # [h*w, 3(ver in tri), 1] + + colors_flat = colors[triangles[triangle_buffer_flat, :], :] # [h*w(tri id in pixel), 3(ver in tri), c(color in ver)] + colors_flat = weight*colors_flat # [h*w, 3, 3] + colors_flat = np.sum(colors_flat, 1) #[h*w, 3]. add tri. + + image = np.reshape(colors_flat, [h, w, c]) + # mask = (triangle_buffer[:,:] > -1).astype(np.float32) + # image = image*mask[:,:,np.newaxis] + return image + + +def render_colors(vertices, triangles, colors, h, w, c = 3): + ''' render mesh with colors + Args: + vertices: [nver, 3] + triangles: [ntri, 3] + colors: [nver, 3] + h: height + w: width + Returns: + image: [h, w, c]. + ''' + assert vertices.shape[0] == colors.shape[0] + + # initial + image = np.zeros((h, w, c)) + depth_buffer = np.zeros([h, w]) - 999999. + + for i in range(triangles.shape[0]): + tri = triangles[i, :] # 3 vertex indices + + # the inner bounding box + umin = max(int(np.ceil(np.min(vertices[tri, 0]))), 0) + umax = min(int(np.floor(np.max(vertices[tri, 0]))), w-1) + + vmin = max(int(np.ceil(np.min(vertices[tri, 1]))), 0) + vmax = min(int(np.floor(np.max(vertices[tri, 1]))), h-1) + + if umax depth_buffer[v, u]: + depth_buffer[v, u] = point_depth + image[v, u, :] = w0*colors[tri[0], :] + w1*colors[tri[1], :] + w2*colors[tri[2], :] + return image + + +def render_texture(vertices, triangles, texture, tex_coords, tex_triangles, h, w, c = 3, mapping_type = 'nearest'): + ''' render mesh with texture map + Args: + vertices: [nver], 3 + triangles: [ntri, 3] + texture: [tex_h, tex_w, 3] + tex_coords: [ntexcoords, 3] + tex_triangles: [ntri, 3] + h: height of rendering + w: width of rendering + c: channel + mapping_type: 'bilinear' or 'nearest' + ''' + assert triangles.shape[0] == tex_triangles.shape[0] + tex_h, tex_w, _ = texture.shape + + # initial + image = np.zeros((h, w, c)) + depth_buffer = np.zeros([h, w]) - 999999. + + for i in range(triangles.shape[0]): + tri = triangles[i, :] # 3 vertex indices + tex_tri = tex_triangles[i, :] # 3 tex indice + + # the inner bounding box + umin = max(int(np.ceil(np.min(vertices[tri, 0]))), 0) + umax = min(int(np.floor(np.max(vertices[tri, 0]))), w-1) + + vmin = max(int(np.ceil(np.min(vertices[tri, 1]))), 0) + vmax = min(int(np.floor(np.max(vertices[tri, 1]))), h-1) + + if umax depth_buffer[v, u]: + # update depth + depth_buffer[v, u] = point_depth + + # tex coord + tex_xy = w0*tex_coords[tex_tri[0], :] + w1*tex_coords[tex_tri[1], :] + w2*tex_coords[tex_tri[2], :] + tex_xy[0] = max(min(tex_xy[0], float(tex_w - 1)), 0.0); + tex_xy[1] = max(min(tex_xy[1], float(tex_h - 1)), 0.0); + + # nearest + if mapping_type == 'nearest': + tex_xy = np.round(tex_xy).astype(np.int32) + tex_value = texture[tex_xy[1], tex_xy[0], :] + + # bilinear + elif mapping_type == 'bilinear': + # next 4 pixels + ul = texture[int(np.floor(tex_xy[1])), int(np.floor(tex_xy[0])), :] + ur = texture[int(np.floor(tex_xy[1])), int(np.ceil(tex_xy[0])), :] + dl = texture[int(np.ceil(tex_xy[1])), int(np.floor(tex_xy[0])), :] + dr = texture[int(np.ceil(tex_xy[1])), int(np.ceil(tex_xy[0])), :] + + yd = tex_xy[1] - np.floor(tex_xy[1]) + xd = tex_xy[0] - np.floor(tex_xy[0]) + tex_value = ul*(1-xd)*(1-yd) + ur*xd*(1-yd) + dl*(1-xd)*yd + dr*xd*yd + + image[v, u, :] = tex_value + return image \ No newline at end of file diff --git a/face3d/mesh_numpy/transform.py b/face3d/mesh_numpy/transform.py new file mode 100644 index 0000000..ab56e37 --- /dev/null +++ b/face3d/mesh_numpy/transform.py @@ -0,0 +1,385 @@ +''' +Functions about transforming mesh(changing the position: modify vertices). +1. forward: transform(transform, camera, project). +2. backward: estimate transform matrix from correspondences. + +Preparation knowledge: +transform&camera model: +https://cs184.eecs.berkeley.edu/lecture/transforms-2 +Part I: camera geometry and single view geometry in MVGCV +''' + +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function + +import numpy as np +import math +from math import cos, sin + +def angle2matrix(angles): + ''' get rotation matrix from three rotation angles(degree). right-handed. + Args: + angles: [3,]. x, y, z angles + x: pitch. positive for looking down. + y: yaw. positive for looking left. + z: roll. positive for tilting head right. + Returns: + R: [3, 3]. rotation matrix. + ''' + x, y, z = np.deg2rad(angles[0]), np.deg2rad(angles[1]), np.deg2rad(angles[2]) + # x + Rx=np.array([[1, 0, 0], + [0, cos(x), -sin(x)], + [0, sin(x), cos(x)]]) + # y + Ry=np.array([[ cos(y), 0, sin(y)], + [ 0, 1, 0], + [-sin(y), 0, cos(y)]]) + # z + Rz=np.array([[cos(z), -sin(z), 0], + [sin(z), cos(z), 0], + [ 0, 0, 1]]) + + R=Rz.dot(Ry.dot(Rx)) + return R.astype(np.float32) + +def angle2matrix_3ddfa(angles): + ''' get rotation matrix from three rotation angles(radian). The same as in 3DDFA. + Args: + angles: [3,]. x, y, z angles + x: pitch. + y: yaw. + z: roll. + Returns: + R: 3x3. rotation matrix. + ''' + # x, y, z = np.deg2rad(angles[0]), np.deg2rad(angles[1]), np.deg2rad(angles[2]) + x, y, z = angles[0], angles[1], angles[2] + + # x + Rx=np.array([[1, 0, 0], + [0, cos(x), sin(x)], + [0, -sin(x), cos(x)]]) + # y + Ry=np.array([[ cos(y), 0, -sin(y)], + [ 0, 1, 0], + [sin(y), 0, cos(y)]]) + # z + Rz=np.array([[cos(z), sin(z), 0], + [-sin(z), cos(z), 0], + [ 0, 0, 1]]) + R = Rx.dot(Ry).dot(Rz) + return R.astype(np.float32) + + +## ------------------------------------------ 1. transform(transform, project, camera). +## ---------- 3d-3d transform. Transform obj in world space +def rotate(vertices, angles): + ''' rotate vertices. + X_new = R.dot(X). X: 3 x 1 + Args: + vertices: [nver, 3]. + rx, ry, rz: degree angles + rx: pitch. positive for looking down + ry: yaw. positive for looking left + rz: roll. positive for tilting head right + Returns: + rotated vertices: [nver, 3] + ''' + R = angle2matrix(angles) + rotated_vertices = vertices.dot(R.T) + + return rotated_vertices + +def similarity_transform(vertices, s, R, t3d): + ''' similarity transform. dof = 7. + 3D: s*R.dot(X) + t + Homo: M = [[sR, t],[0^T, 1]]. M.dot(X) + Args:(float32) + vertices: [nver, 3]. + s: [1,]. scale factor. + R: [3,3]. rotation matrix. + t3d: [3,]. 3d translation vector. + Returns: + transformed vertices: [nver, 3] + ''' + t3d = np.squeeze(np.array(t3d, dtype = np.float32)) + transformed_vertices = s * vertices.dot(R.T) + t3d[np.newaxis, :] + + return transformed_vertices + + +## -------------- Camera. from world space to camera space +# Ref: https://cs184.eecs.berkeley.edu/lecture/transforms-2 +def normalize(x): + epsilon = 1e-12 + norm = np.sqrt(np.sum(x**2, axis = 0)) + norm = np.maximum(norm, epsilon) + return x/norm + +def lookat_camera(vertices, eye, at = None, up = None): + """ 'look at' transformation: from world space to camera space + standard camera space: + camera located at the origin. + looking down negative z-axis. + vertical vector is y-axis. + Xcam = R(X - C) + Homo: [[R, -RC], [0, 1]] + Args: + vertices: [nver, 3] + eye: [3,] the XYZ world space position of the camera. + at: [3,] a position along the center of the camera's gaze. + up: [3,] up direction + Returns: + transformed_vertices: [nver, 3] + """ + if at is None: + at = np.array([0, 0, 0], np.float32) + if up is None: + up = np.array([0, 1, 0], np.float32) + + eye = np.array(eye).astype(np.float32) + at = np.array(at).astype(np.float32) + z_aixs = -normalize(at - eye) # look forward + x_aixs = normalize(np.cross(up, z_aixs)) # look right + y_axis = np.cross(z_aixs, x_aixs) # look up + + R = np.stack((x_aixs, y_axis, z_aixs))#, axis = 0) # 3 x 3 + transformed_vertices = vertices - eye # translation + transformed_vertices = transformed_vertices.dot(R.T) # rotation + return transformed_vertices + +## --------- 3d-2d project. from camera space to image plane +# generally, image plane only keeps x,y channels, here reserve z channel for calculating z-buffer. +def orthographic_project(vertices): + ''' scaled orthographic projection(just delete z) + assumes: variations in depth over the object is small relative to the mean distance from camera to object + x -> x*f/z, y -> x*f/z, z -> f. + for point i,j. zi~=zj. so just delete z + ** often used in face + Homo: P = [[1,0,0,0], [0,1,0,0], [0,0,1,0]] + Args: + vertices: [nver, 3] + Returns: + projected_vertices: [nver, 3] if isKeepZ=True. [nver, 2] if isKeepZ=False. + ''' + return vertices.copy() + +def perspective_project(vertices, fovy, aspect_ratio = 1., near = 0.1, far = 1000.): + ''' perspective projection. + Args: + vertices: [nver, 3] + fovy: vertical angular field of view. degree. + aspect_ratio : width / height of field of view + near : depth of near clipping plane + far : depth of far clipping plane + Returns: + projected_vertices: [nver, 3] + ''' + fovy = np.deg2rad(fovy) + top = near*np.tan(fovy) + bottom = -top + right = top*aspect_ratio + left = -right + + #-- homo + P = np.array([[near/right, 0, 0, 0], + [0, near/top, 0, 0], + [0, 0, -(far+near)/(far-near), -2*far*near/(far-near)], + [0, 0, -1, 0]]) + vertices_homo = np.hstack((vertices, np.ones((vertices.shape[0], 1)))) # [nver, 4] + projected_vertices = vertices_homo.dot(P.T) + projected_vertices = projected_vertices/projected_vertices[:,3:] + projected_vertices = projected_vertices[:,:3] + projected_vertices[:,2] = -projected_vertices[:,2] + + #-- non homo. only fovy + # projected_vertices = vertices.copy() + # projected_vertices[:,0] = -(near/right)*vertices[:,0]/vertices[:,2] + # projected_vertices[:,1] = -(near/top)*vertices[:,1]/vertices[:,2] + return projected_vertices + + +def to_image(vertices, h, w, is_perspective = False): + ''' change vertices to image coord system + 3d system: XYZ, center(0, 0, 0) + 2d image: x(u), y(v). center(w/2, h/2), flip y-axis. + Args: + vertices: [nver, 3] + h: height of the rendering + w : width of the rendering + Returns: + projected_vertices: [nver, 3] + ''' + image_vertices = vertices.copy() + if is_perspective: + # if perspective, the projected vertices are normalized to [-1, 1]. so change it to image size first. + image_vertices[:,0] = image_vertices[:,0]*w/2 + image_vertices[:,1] = image_vertices[:,1]*h/2 + # move to center of image + image_vertices[:,0] = image_vertices[:,0] + w/2 + image_vertices[:,1] = image_vertices[:,1] + h/2 + # flip vertices along y-axis. + image_vertices[:,1] = h - image_vertices[:,1] - 1 + return image_vertices + + +#### -------------------------------------------2. estimate transform matrix from correspondences. +def estimate_affine_matrix_3d23d(X, Y): + ''' Using least-squares solution + Args: + X: [n, 3]. 3d points(fixed) + Y: [n, 3]. corresponding 3d points(moving). Y = PX + Returns: + P_Affine: (3, 4). Affine camera matrix (the third row is [0, 0, 0, 1]). + ''' + X_homo = np.hstack((X, np.ones([X.shape[1],1]))) #n x 4 + P = np.linalg.lstsq(X_homo, Y)[0].T # Affine matrix. 3 x 4 + return P + +def estimate_affine_matrix_3d22d(X, x): + ''' Using Golden Standard Algorithm for estimating an affine camera + matrix P from world to image correspondences. + See Alg.7.2. in MVGCV + Code Ref: https://github.com/patrikhuber/eos/blob/master/include/eos/fitting/affine_camera_estimation.hpp + x_homo = X_homo.dot(P_Affine) + Args: + X: [n, 3]. corresponding 3d points(fixed) + x: [n, 2]. n>=4. 2d points(moving). x = PX + Returns: + P_Affine: [3, 4]. Affine camera matrix + ''' + X = X.T; x = x.T + assert(x.shape[1] == X.shape[1]) + n = x.shape[1] + assert(n >= 4) + + #--- 1. normalization + # 2d points + mean = np.mean(x, 1) # (2,) + x = x - np.tile(mean[:, np.newaxis], [1, n]) + average_norm = np.mean(np.sqrt(np.sum(x**2, 0))) + scale = np.sqrt(2) / average_norm + x = scale * x + + T = np.zeros((3,3), dtype = np.float32) + T[0, 0] = T[1, 1] = scale + T[:2, 2] = -mean*scale + T[2, 2] = 1 + + # 3d points + X_homo = np.vstack((X, np.ones((1, n)))) + mean = np.mean(X, 1) # (3,) + X = X - np.tile(mean[:, np.newaxis], [1, n]) + m = X_homo[:3,:] - X + average_norm = np.mean(np.sqrt(np.sum(X**2, 0))) + scale = np.sqrt(3) / average_norm + X = scale * X + + U = np.zeros((4,4), dtype = np.float32) + U[0, 0] = U[1, 1] = U[2, 2] = scale + U[:3, 3] = -mean*scale + U[3, 3] = 1 + + # --- 2. equations + A = np.zeros((n*2, 8), dtype = np.float32); + X_homo = np.vstack((X, np.ones((1, n)))).T + A[:n, :4] = X_homo + A[n:, 4:] = X_homo + b = np.reshape(x, [-1, 1]) + + # --- 3. solution + p_8 = np.linalg.pinv(A).dot(b) + P = np.zeros((3, 4), dtype = np.float32) + P[0, :] = p_8[:4, 0] + P[1, :] = p_8[4:, 0] + P[-1, -1] = 1 + + # --- 4. denormalization + P_Affine = np.linalg.inv(T).dot(P.dot(U)) + return P_Affine + +def P2sRt(P): + ''' decompositing camera matrix P + Args: + P: (3, 4). Affine Camera Matrix. + Returns: + s: scale factor. + R: (3, 3). rotation matrix. + t: (3,). translation. + ''' + t = P[:, 3] + R1 = P[0:1, :3] + R2 = P[1:2, :3] + s = (np.linalg.norm(R1) + np.linalg.norm(R2))/2.0 + r1 = R1/np.linalg.norm(R1) + r2 = R2/np.linalg.norm(R2) + r3 = np.cross(r1, r2) + + R = np.concatenate((r1, r2, r3), 0) + return s, R, t + +#Ref: https://www.learnopencv.com/rotation-matrix-to-euler-angles/ +def isRotationMatrix(R): + ''' checks if a matrix is a valid rotation matrix(whether orthogonal or not) + ''' + Rt = np.transpose(R) + shouldBeIdentity = np.dot(Rt, R) + I = np.identity(3, dtype = R.dtype) + n = np.linalg.norm(I - shouldBeIdentity) + return n < 1e-6 + +def matrix2angle(R): + ''' get three Euler angles from Rotation Matrix + Args: + R: (3,3). rotation matrix + Returns: + x: pitch + y: yaw + z: roll + ''' + assert(isRotationMatrix) + sy = math.sqrt(R[0,0] * R[0,0] + R[1,0] * R[1,0]) + + singular = sy < 1e-6 + + if not singular : + x = math.atan2(R[2,1] , R[2,2]) + y = math.atan2(-R[2,0], sy) + z = math.atan2(R[1,0], R[0,0]) + else : + x = math.atan2(-R[1,2], R[1,1]) + y = math.atan2(-R[2,0], sy) + z = 0 + + # rx, ry, rz = np.rad2deg(x), np.rad2deg(y), np.rad2deg(z) + rx, ry, rz = x*180/np.pi, y*180/np.pi, z*180/np.pi + return rx, ry, rz + +# def matrix2angle(R): +# ''' compute three Euler angles from a Rotation Matrix. Ref: http://www.gregslabaugh.net/publications/euler.pdf +# Args: +# R: (3,3). rotation matrix +# Returns: +# x: yaw +# y: pitch +# z: roll +# ''' +# # assert(isRotationMatrix(R)) + +# if R[2,0] !=1 or R[2,0] != -1: +# x = math.asin(R[2,0]) +# y = math.atan2(R[2,1]/cos(x), R[2,2]/cos(x)) +# z = math.atan2(R[1,0]/cos(x), R[0,0]/cos(x)) + +# else:# Gimbal lock +# z = 0 #can be anything +# if R[2,0] == -1: +# x = np.pi/2 +# y = z + math.atan2(R[0,1], R[0,2]) +# else: +# x = -np.pi/2 +# y = -z + math.atan2(-R[0,1], -R[0,2]) + +# return x, y, z \ No newline at end of file diff --git a/face3d/mesh_numpy/vis.py b/face3d/mesh_numpy/vis.py new file mode 100644 index 0000000..9db972f --- /dev/null +++ b/face3d/mesh_numpy/vis.py @@ -0,0 +1,24 @@ +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function + +import numpy as np +import matplotlib.pyplot as plt +from skimage import measure +from mpl_toolkits.mplot3d import Axes3D + +def plot_mesh(vertices, triangles, subplot = [1,1,1], title = 'mesh', el = 90, az = -90, lwdt=.1, dist = 6, color = "grey"): + ''' + plot the mesh + Args: + vertices: [nver, 3] + triangles: [ntri, 3] + ''' + ax = plt.subplot(subplot[0], subplot[1], subplot[2], projection = '3d') + ax.plot_trisurf(vertices[:, 0], vertices[:, 1], vertices[:, 2], triangles = triangles, lw = lwdt, color = color, alpha = 1) + ax.axis("off") + ax.view_init(elev = el, azim = az) + ax.dist = dist + plt.title(title) + +### -------------- Todo: use vtk to visualize mesh? or visvis? or VisPy?