fix: remove face3d and fix gradio version (#26)

* feat: gradio sticthing option

* fix: face3d

* chore: add gradio
This commit is contained in:
ZhizhouZhong 2024-07-05 14:11:51 +08:00 committed by GitHub
parent 420b4a08c7
commit 2f808ca1d6
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
37 changed files with 13 additions and 24391 deletions

12
app.py
View File

@ -35,11 +35,11 @@ title_md = "assets/gradio_title.md"
example_portrait_dir = "assets/examples/source"
example_video_dir = "assets/examples/driving"
data_examples = [
[osp.join(example_portrait_dir, "s1.jpg"), osp.join(example_video_dir, "d1.mp4"), True, True, True],
[osp.join(example_portrait_dir, "s2.jpg"), osp.join(example_video_dir, "d2.mp4"), True, True, True],
[osp.join(example_portrait_dir, "s3.jpg"), osp.join(example_video_dir, "d5.mp4"), True, True, True],
[osp.join(example_portrait_dir, "s5.jpg"), osp.join(example_video_dir, "d6.mp4"), True, True, True],
[osp.join(example_portrait_dir, "s7.jpg"), osp.join(example_video_dir, "d7.mp4"), True, True, True],
[osp.join(example_portrait_dir, "s9.jpg"), osp.join(example_video_dir, "d0.mp4"), True, True, True, True],
[osp.join(example_portrait_dir, "s6.jpg"), osp.join(example_video_dir, "d0.mp4"), True, True, True, True],
[osp.join(example_portrait_dir, "s10.jpg"), osp.join(example_video_dir, "d5.mp4"), True, True, True, True],
[osp.join(example_portrait_dir, "s5.jpg"), osp.join(example_video_dir, "d6.mp4"), True, True, True, True],
[osp.join(example_portrait_dir, "s7.jpg"), osp.join(example_video_dir, "d7.mp4"), True, True, True, True],
]
#################### interface logic ####################
@ -65,8 +65,8 @@ with gr.Blocks(theme=gr.themes.Soft()) as demo:
with gr.Accordion(open=True, label="Animation Options"):
with gr.Row():
flag_relative_input = gr.Checkbox(value=True, label="relative pose")
flag_remap_input = gr.Checkbox(value=True, label="paste-back")
flag_do_crop_input = gr.Checkbox(value=True, label="do crop")
flag_remap_input = gr.Checkbox(value=True, label="paste-back")
with gr.Row():
with gr.Column():
process_button_animation = gr.Button("🚀 Animate", variant="primary")

View File

@ -105,7 +105,6 @@ python inference.py -h
We also provide a Gradio interface for a better experience. Please install `gradio` and then run `app.py`:
```bash
pip install gradio==4.37.1
python app.py
```

View File

@ -20,3 +20,4 @@ albumentations==1.4.10
matplotlib==3.9.0
imageio-ffmpeg==0.5.1
tyro==0.8.5
gradio==4.37.1

View File

@ -43,7 +43,7 @@ class GradioPipeline(LivePortraitPipeline):
input_video_path,
flag_relative_input,
flag_do_crop_input,
flag_remap_input
flag_remap_input,
):
""" for video driven potrait animation
"""
@ -53,7 +53,7 @@ class GradioPipeline(LivePortraitPipeline):
'driving_info': input_video_path,
'flag_relative': flag_relative_input,
'flag_do_crop': flag_do_crop_input,
'flag_pasteback': flag_remap_input
'flag_pasteback': flag_remap_input,
}
# update config from user input
self.args = update_args(self.args, args_user)

View File

@ -17,5 +17,4 @@ from . import model_zoo
from . import utils
from . import app
from . import data
from . import thirdparty

View File

@ -1,2 +1 @@
from .face_analysis import *
from .mask_renderer import *

View File

@ -15,9 +15,11 @@ import onnxruntime
from numpy.linalg import norm
from ..model_zoo import model_zoo
from ..utils import DEFAULT_MP_NAME, ensure_available
from ..utils import ensure_available
from .common import Face
DEFAULT_MP_NAME = 'buffalo_l'
__all__ = ['FaceAnalysis']
class FaceAnalysis:

View File

@ -1,232 +0,0 @@
import os, sys, datetime
import numpy as np
import os.path as osp
import albumentations as A
from albumentations.core.transforms_interface import ImageOnlyTransform
from .face_analysis import FaceAnalysis
from ..utils import get_model_dir
from ..thirdparty import face3d
from ..data import get_image as ins_get_image
from ..utils import DEFAULT_MP_NAME
import cv2
class MaskRenderer:
def __init__(self, name=DEFAULT_MP_NAME, root='~/.insightface', insfa=None):
#if insfa is None, enter render_only mode
self.mp_name = name
self.root = root
self.insfa = insfa
model_dir = get_model_dir(name, root)
bfm_file = osp.join(model_dir, 'BFM.mat')
assert osp.exists(bfm_file), 'should contains BFM.mat in your model directory'
self.bfm = face3d.morphable_model.MorphabelModel(bfm_file)
self.index_ind = self.bfm.kpt_ind
bfm_uv_file = osp.join(model_dir, 'BFM_UV.mat')
assert osp.exists(bfm_uv_file), 'should contains BFM_UV.mat in your model directory'
uv_coords = face3d.morphable_model.load.load_uv_coords(bfm_uv_file)
self.uv_size = (224,224)
self.mask_stxr = 0.1
self.mask_styr = 0.33
self.mask_etxr = 0.9
self.mask_etyr = 0.7
self.tex_h , self.tex_w, self.tex_c = self.uv_size[1] , self.uv_size[0],3
texcoord = np.zeros_like(uv_coords)
texcoord[:, 0] = uv_coords[:, 0] * (self.tex_h - 1)
texcoord[:, 1] = uv_coords[:, 1] * (self.tex_w - 1)
texcoord[:, 1] = self.tex_w - texcoord[:, 1] - 1
self.texcoord = np.hstack((texcoord, np.zeros((texcoord.shape[0], 1))))
self.X_ind = self.bfm.kpt_ind
self.mask_image_names = ['mask_white', 'mask_blue', 'mask_black', 'mask_green']
self.mask_aug_probs = [0.4, 0.4, 0.1, 0.1]
#self.mask_images = []
#self.mask_images_rgb = []
#for image_name in mask_image_names:
# mask_image = ins_get_image(image_name)
# self.mask_images.append(mask_image)
# mask_image_rgb = mask_image[:,:,::-1]
# self.mask_images_rgb.append(mask_image_rgb)
def prepare(self, ctx_id=0, det_thresh=0.5, det_size=(128, 128)):
self.pre_ctx_id = ctx_id
self.pre_det_thresh = det_thresh
self.pre_det_size = det_size
def transform(self, shape3D, R):
s = 1.0
shape3D[:2, :] = shape3D[:2, :]
shape3D = s * np.dot(R, shape3D)
return shape3D
def preprocess(self, vertices, w, h):
R1 = face3d.mesh.transform.angle2matrix([0, 180, 180])
t = np.array([-w // 2, -h // 2, 0])
vertices = vertices.T
vertices += t
vertices = self.transform(vertices.T, R1).T
return vertices
def project_to_2d(self,vertices,s,angles,t):
transformed_vertices = self.bfm.transform(vertices, s, angles, t)
projected_vertices = transformed_vertices.copy() # using stantard camera & orth projection
return projected_vertices[self.bfm.kpt_ind, :2]
def params_to_vertices(self,params , H , W):
fitted_sp, fitted_ep, fitted_s, fitted_angles, fitted_t = params
fitted_vertices = self.bfm.generate_vertices(fitted_sp, fitted_ep)
transformed_vertices = self.bfm.transform(fitted_vertices, fitted_s, fitted_angles,
fitted_t)
transformed_vertices = self.preprocess(transformed_vertices.T, W, H)
image_vertices = face3d.mesh.transform.to_image(transformed_vertices, H, W)
return image_vertices
def draw_lmk(self, face_image):
faces = self.insfa.get(face_image, max_num=1)
if len(faces)==0:
return face_image
return self.insfa.draw_on(face_image, faces)
def build_params(self, face_image):
#landmark = self.if3d68_handler.get(face_image)
#if landmark is None:
# return None #face not found
if self.insfa is None:
self.insfa = FaceAnalysis(name=self.mp_name, root=self.root, allowed_modules=['detection', 'landmark_3d_68'])
self.insfa.prepare(ctx_id=self.pre_ctx_id, det_thresh=self.pre_det_thresh, det_size=self.pre_det_size)
faces = self.insfa.get(face_image, max_num=1)
if len(faces)==0:
return None
landmark = faces[0].landmark_3d_68[:,:2]
fitted_sp, fitted_ep, fitted_s, fitted_angles, fitted_t = self.bfm.fit(landmark, self.X_ind, max_iter = 3)
return [fitted_sp, fitted_ep, fitted_s, fitted_angles, fitted_t]
def generate_mask_uv(self,mask, positions):
uv_size = (self.uv_size[1], self.uv_size[0], 3)
h, w, c = uv_size
uv = np.zeros(shape=(self.uv_size[1],self.uv_size[0], 3), dtype=np.uint8)
stxr, styr = positions[0], positions[1]
etxr, etyr = positions[2], positions[3]
stx, sty = int(w * stxr), int(h * styr)
etx, ety = int(w * etxr), int(h * etyr)
height = ety - sty
width = etx - stx
mask = cv2.resize(mask, (width, height))
uv[sty:ety, stx:etx] = mask
return uv
def render_mask(self,face_image, mask_image, params, input_is_rgb=False, auto_blend = True, positions=[0.1, 0.33, 0.9, 0.7]):
if isinstance(mask_image, str):
to_rgb = True if input_is_rgb else False
mask_image = ins_get_image(mask_image, to_rgb=to_rgb)
uv_mask_image = self.generate_mask_uv(mask_image, positions)
h,w,c = face_image.shape
image_vertices = self.params_to_vertices(params ,h,w)
output = (1-face3d.mesh.render.render_texture(image_vertices, self.bfm.full_triangles , uv_mask_image, self.texcoord, self.bfm.full_triangles, h , w ))*255
output = output.astype(np.uint8)
if auto_blend:
mask_bd = (output==255).astype(np.uint8)
final = face_image*mask_bd + (1-mask_bd)*output
return final
return output
#def mask_augmentation(self, face_image, label, input_is_rgb=False, p=0.1):
# if np.random.random()<p:
# assert isinstance(label, (list, np.ndarray)), 'make sure the rec dataset includes mask params'
# assert len(label)==237 or len(lable)==235, 'make sure the rec dataset includes mask params'
# if len(label)==237:
# if label[1]<0.0: #invalid label for mask aug
# return face_image
# label = label[2:]
# params = self.decode_params(label)
# mask_image_name = np.random.choice(self.mask_image_names, p=self.mask_aug_probs)
# pos = np.random.uniform(0.33, 0.5)
# face_image = self.render_mask(face_image, mask_image_name, params, input_is_rgb=input_is_rgb, positions=[0.1, pos, 0.9, 0.7])
# return face_image
@staticmethod
def encode_params(params):
p0 = list(params[0])
p1 = list(params[1])
p2 = [float(params[2])]
p3 = list(params[3])
p4 = list(params[4])
return p0+p1+p2+p3+p4
@staticmethod
def decode_params(params):
p0 = params[0:199]
p0 = np.array(p0, dtype=np.float32).reshape( (-1, 1))
p1 = params[199:228]
p1 = np.array(p1, dtype=np.float32).reshape( (-1, 1))
p2 = params[228]
p3 = tuple(params[229:232])
p4 = params[232:235]
p4 = np.array(p4, dtype=np.float32).reshape( (-1, 1))
return p0, p1, p2, p3, p4
class MaskAugmentation(ImageOnlyTransform):
def __init__(
self,
mask_names=['mask_white', 'mask_blue', 'mask_black', 'mask_green'],
mask_probs=[0.4,0.4,0.1,0.1],
h_low = 0.33,
h_high = 0.35,
always_apply=False,
p=1.0,
):
super(MaskAugmentation, self).__init__(always_apply, p)
self.renderer = MaskRenderer()
assert len(mask_names)>0
assert len(mask_names)==len(mask_probs)
self.mask_names = mask_names
self.mask_probs = mask_probs
self.h_low = h_low
self.h_high = h_high
#self.hlabel = None
def apply(self, image, hlabel, mask_name, h_pos, **params):
#print(params.keys())
#hlabel = params.get('hlabel')
assert len(hlabel)==237 or len(hlabel)==235, 'make sure the rec dataset includes mask params'
if len(hlabel)==237:
if hlabel[1]<0.0:
return image
hlabel = hlabel[2:]
#print(len(hlabel))
mask_params = self.renderer.decode_params(hlabel)
image = self.renderer.render_mask(image, mask_name, mask_params, input_is_rgb=True, positions=[0.1, h_pos, 0.9, 0.7])
return image
@property
def targets_as_params(self):
return ["image", "hlabel"]
def get_params_dependent_on_targets(self, params):
hlabel = params['hlabel']
mask_name = np.random.choice(self.mask_names, p=self.mask_probs)
h_pos = np.random.uniform(self.h_low, self.h_high)
return {'hlabel': hlabel, 'mask_name': mask_name, 'h_pos': h_pos}
def get_transform_init_args_names(self):
#return ("hlabel", 'mask_names', 'mask_probs', 'h_low', 'h_high')
return ('mask_names', 'mask_probs', 'h_low', 'h_high')
if __name__ == "__main__":
tool = MaskRenderer('antelope')
tool.prepare(det_size=(128,128))
image = cv2.imread("Tom_Hanks_54745.png")
params = tool.build_params(image)
#out = tool.draw_lmk(image)
#cv2.imwrite('output_lmk.jpg', out)
#mask_image = cv2.imread("masks/mask1.jpg")
#mask_image = cv2.imread("masks/black-mask.png")
#mask_image = cv2.imread("masks/mask2.jpg")
mask_out = tool.render_mask(image, 'mask_blue', params)# use single thread to test the time cost
cv2.imwrite('output_mask.jpg', mask_out)

View File

@ -1,13 +0,0 @@
from abc import ABC, abstractmethod
from argparse import ArgumentParser
class BaseInsightFaceCLICommand(ABC):
@staticmethod
@abstractmethod
def register_subcommand(parser: ArgumentParser):
raise NotImplementedError()
@abstractmethod
def run(self):
raise NotImplementedError()

View File

@ -1,29 +0,0 @@
#!/usr/bin/env python
from argparse import ArgumentParser
from .model_download import ModelDownloadCommand
from .rec_add_mask_param import RecAddMaskParamCommand
def main():
parser = ArgumentParser("InsightFace CLI tool", usage="insightface-cli <command> [<args>]")
commands_parser = parser.add_subparsers(help="insightface-cli command-line helpers")
# Register commands
ModelDownloadCommand.register_subcommand(commands_parser)
RecAddMaskParamCommand.register_subcommand(commands_parser)
args = parser.parse_args()
if not hasattr(args, "func"):
parser.print_help()
exit(1)
# Run
service = args.func(args)
service.run()
if __name__ == "__main__":
main()

View File

@ -1,36 +0,0 @@
from argparse import ArgumentParser
from . import BaseInsightFaceCLICommand
import os
import os.path as osp
import zipfile
import glob
from ..utils import download
def model_download_command_factory(args):
return ModelDownloadCommand(args.model, args.root, args.force)
class ModelDownloadCommand(BaseInsightFaceCLICommand):
#_url_format = '{repo_url}models/{file_name}.zip'
@staticmethod
def register_subcommand(parser: ArgumentParser):
download_parser = parser.add_parser("model.download")
download_parser.add_argument(
"--root", type=str, default='~/.insightface', help="Path to location to store the models"
)
download_parser.add_argument(
"--force", action="store_true", help="Force the model to be download even if already in root-dir"
)
download_parser.add_argument("model", type=str, help="Name of the model to download")
download_parser.set_defaults(func=model_download_command_factory)
def __init__(self, model: str, root: str, force: bool):
self._model = model
self._root = root
self._force = force
def run(self):
download('models', self._model, force=self._force, root=self._root)

View File

@ -1,94 +0,0 @@
import numbers
import os
from argparse import ArgumentParser, Namespace
import mxnet as mx
import numpy as np
from ..app import MaskRenderer
from ..data.rec_builder import RecBuilder
from . import BaseInsightFaceCLICommand
def rec_add_mask_param_command_factory(args: Namespace):
return RecAddMaskParamCommand(
args.input, args.output
)
class RecAddMaskParamCommand(BaseInsightFaceCLICommand):
@staticmethod
def register_subcommand(parser: ArgumentParser):
_parser = parser.add_parser("rec.addmaskparam")
_parser.add_argument("input", type=str, help="input rec")
_parser.add_argument("output", type=str, help="output rec, with mask param")
_parser.set_defaults(func=rec_add_mask_param_command_factory)
def __init__(
self,
input: str,
output: str,
):
self._input = input
self._output = output
def run(self):
tool = MaskRenderer()
tool.prepare(ctx_id=0, det_size=(128,128))
root_dir = self._input
path_imgrec = os.path.join(root_dir, 'train.rec')
path_imgidx = os.path.join(root_dir, 'train.idx')
imgrec = mx.recordio.MXIndexedRecordIO(path_imgidx, path_imgrec, 'r')
save_path = self._output
wrec=RecBuilder(path=save_path)
s = imgrec.read_idx(0)
header, _ = mx.recordio.unpack(s)
if header.flag > 0:
if len(header.label)==2:
imgidx = np.array(range(1, int(header.label[0])))
else:
imgidx = np.array(list(self.imgrec.keys))
else:
imgidx = np.array(list(self.imgrec.keys))
stat = [0, 0]
print('total:', len(imgidx))
for iid, idx in enumerate(imgidx):
#if iid==500000:
# break
if iid%1000==0:
print('processing:', iid)
s = imgrec.read_idx(idx)
header, img = mx.recordio.unpack(s)
label = header.label
if not isinstance(label, numbers.Number):
label = label[0]
sample = mx.image.imdecode(img).asnumpy()
bgr = sample[:,:,::-1]
params = tool.build_params(bgr)
#if iid<10:
# mask_out = tool.render_mask(bgr, 'mask_blue', params)
# cv2.imwrite('maskout_%d.jpg'%iid, mask_out)
stat[1] += 1
if params is None:
wlabel = [label] + [-1.0]*236
stat[0] += 1
else:
#print(0, params[0].shape, params[0].dtype)
#print(1, params[1].shape, params[1].dtype)
#print(2, params[2])
#print(3, len(params[3]), params[3][0].__class__)
#print(4, params[4].shape, params[4].dtype)
mask_label = tool.encode_params(params)
wlabel = [label, 0.0]+mask_label # 237 including idlabel, total mask params size is 235
if iid==0:
print('param size:', len(mask_label), len(wlabel), label)
assert len(wlabel)==237
wrec.add_image(img, wlabel)
#print(len(params))
wrec.close()
print('finished on', self._output, ', failed:', stat[0])

View File

@ -1,4 +0,0 @@
#import mesh
#import morphable_model
from . import mesh
from . import morphable_model

View File

@ -1,15 +0,0 @@
#from __future__ import absolute_import
#from cython import mesh_core_cython
#import io
#import vis
#import transform
#import light
#import render
from .cython import mesh_core_cython
from . import io
from . import vis
from . import transform
from . import light
from . import render

View File

@ -1,375 +0,0 @@
/*
functions that can not be optimazed by vertorization in python.
1. rasterization.(need process each triangle)
2. normal of each vertex.(use one-ring, need process each vertex)
3. write obj(seems that it can be verctorized? anyway, writing it in c++ is simple, so also add function here. --> however, why writting in c++ is still slow?)
Author: Yao Feng
Mail: yaofeng1995@gmail.com
*/
#include "mesh_core.h"
/* Judge whether the point is in the triangle
Method:
http://blackpawn.com/texts/pointinpoly/
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;
}
}
}
}
}
}
// ------------------------------------------------- write
// obj write
// Ref: https://github.com/patrikhuber/eos/blob/master/include/eos/core/Mesh.hpp
void _write_obj_with_colors_texture(string filename, string mtl_name,
float* vertices, int* triangles, float* colors, float* uv_coords,
int nver, int ntri, int ntexver)
{
int i;
ofstream obj_file(filename.c_str());
// first line of the obj file: the mtl name
obj_file << "mtllib " << mtl_name << endl;
// write vertices
for (i = 0; i < nver; ++i)
{
obj_file << "v " << vertices[3*i] << " " << vertices[3*i + 1] << " " << vertices[3*i + 2] << colors[3*i] << " " << colors[3*i + 1] << " " << colors[3*i + 2] << endl;
}
// write uv coordinates
for (i = 0; i < ntexver; ++i)
{
//obj_file << "vt " << uv_coords[2*i] << " " << (1 - uv_coords[2*i + 1]) << endl;
obj_file << "vt " << uv_coords[2*i] << " " << uv_coords[2*i + 1] << endl;
}
obj_file << "usemtl FaceTexture" << endl;
// write triangles
for (i = 0; i < ntri; ++i)
{
// obj_file << "f " << triangles[3*i] << "/" << triangles[3*i] << " " << triangles[3*i + 1] << "/" << triangles[3*i + 1] << " " << triangles[3*i + 2] << "/" << triangles[3*i + 2] << endl;
obj_file << "f " << triangles[3*i + 2] << "/" << triangles[3*i + 2] << " " << triangles[3*i + 1] << "/" << triangles[3*i + 1] << " " << triangles[3*i] << "/" << triangles[3*i] << endl;
}
}

View File

@ -1,83 +0,0 @@
#ifndef MESH_CORE_HPP_
#define MESH_CORE_HPP_
#include <stdio.h>
#include <cmath>
#include <algorithm>
#include <string>
#include <iostream>
#include <fstream>
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);
void _write_obj_with_colors_texture(string filename, string mtl_name,
float* vertices, int* triangles, float* colors, float* uv_coords,
int nver, int ntri, int ntexver);
#endif

View File

@ -1,109 +0,0 @@
import numpy as np
cimport numpy as np
from libcpp.string cimport string
# 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)
void _write_obj_with_colors_texture(string filename, string mtl_name,
float* vertices, int* triangles, float* colors, float* uv_coords,
int nver, int ntri, int ntexver)
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(
<float*> np.PyArray_DATA(normal), <float*> np.PyArray_DATA(tri_normal), <int*> 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(
<float*> np.PyArray_DATA(vertices), <int*> np.PyArray_DATA(triangles),
<float*> np.PyArray_DATA(depth_buffer), <int*> np.PyArray_DATA(triangle_buffer), <float*> 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(
<float*> np.PyArray_DATA(image), <float*> np.PyArray_DATA(vertices), <int*> np.PyArray_DATA(triangles),
<float*> np.PyArray_DATA(colors),
<float*> 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(
<float*> np.PyArray_DATA(image), <float*> np.PyArray_DATA(vertices), <int*> np.PyArray_DATA(triangles),
<float*> np.PyArray_DATA(texture), <float*> np.PyArray_DATA(tex_coords), <int*> np.PyArray_DATA(tex_triangles),
<float*> np.PyArray_DATA(depth_buffer),
nver, tex_nver, ntri,
h, w, c,
tex_h, tex_w, tex_c,
mapping_type)
def write_obj_with_colors_texture_core(string filename, string mtl_name,
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"] uv_coords not None,
int nver, int ntri, int ntexver
):
_write_obj_with_colors_texture(filename, mtl_name,
<float*> np.PyArray_DATA(vertices), <int*> np.PyArray_DATA(triangles), <float*> np.PyArray_DATA(colors), <float*> np.PyArray_DATA(uv_coords),
nver, ntri, ntexver)

View File

@ -1,20 +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(
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()])],
)

View File

@ -1,142 +0,0 @@
from __future__ import absolute_import
from __future__ import division
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
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])
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)
# 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)
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
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")
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)

View File

@ -1,213 +0,0 @@
'''
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.
'''
from __future__ import absolute_import
from __future__ import division
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
Args:
vertices: [nver, 3]
triangles: [ntri, 3]
Returns:
normal: [nver, 3]
'''
pt0 = vertices[triangles[:, 0], :] # [ntri, 3]
pt1 = vertices[triangles[:, 1], :] # [ntri, 3]
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])
# normalize to unit length
mag = np.sum(normal**2, 1) # [nver]
zero_ind = (mag == 0)
mag[zero_ind] = 1;
normal[zero_ind, 0] = np.ones((np.sum(zero_ind)))
normal = normal/np.sqrt(mag[:,np.newaxis])
return normal
# TODO: test
def add_light_sh(vertices, triangles, colors, sh_coeff):
'''
In 3d face, usually assume:
1. The surface of face is Lambertian(reflect only the low frequencies of lighting)
2. Lighting can be an arbitrary combination of point sources
--> can be expressed in terms of spherical harmonics(omit the lighting coefficients)
I = albedo * (sh(n) x sh_coeff)
albedo: n x 1
sh_coeff: 9 x 1
Y(n) = (1, n_x, n_y, n_z, n_xn_y, n_xn_z, n_yn_z, n_x^2 - n_y^2, 3n_z^2 - 1)': n x 9
# Y(n) = (1, n_x, n_y, n_z)': n x 4
Args:
vertices: [nver, 3]
triangles: [ntri, 3]
colors: [nver, 3] albedo
sh_coeff: [9, 1] spherical harmonics coefficients
Returns:
lit_colors: [nver, 3]
'''
assert vertices.shape[0] == colors.shape[0]
nver = vertices.shape[0]
normal = get_normal(vertices, triangles) # [nver, 3]
sh = np.array((np.ones(nver), n[:,0], n[:,1], n[:,2], n[:,0]*n[:,1], n[:,0]*n[:,2], n[:,1]*n[:,2], n[:,0]**2 - n[:,1]**2, 3*(n[:,2]**2) - 1)) # [nver, 9]
ref = sh.dot(sh_coeff) #[nver, 1]
lit_colors = colors*ref
return lit_colors
def add_light(vertices, triangles, colors, light_positions = 0, light_intensities = 0):
''' Gouraud shading. add point lights.
In 3d face, usually assume:
1. The surface of face is Lambertian(reflect only the low frequencies of lighting)
2. Lighting can be an arbitrary combination of point sources
3. No specular (unless skin is oil, 23333)
Ref: https://cs184.eecs.berkeley.edu/lecture/pipeline
Args:
vertices: [nver, 3]
triangles: [ntri, 3]
light_positions: [nlight, 3]
light_intensities: [nlight, 3]
Returns:
lit_colors: [nver, 3]
'''
nver = vertices.shape[0]
normals = get_normal(vertices, triangles) # [nver, 3]
# ambient
# La = ka*Ia
# diffuse
# Ld = kd*(I/r^2)max(0, nxl)
direction_to_lights = vertices[np.newaxis, :, :] - light_positions[:, np.newaxis, :] # [nlight, nver, 3]
direction_to_lights_n = np.sqrt(np.sum(direction_to_lights**2, axis = 2)) # [nlight, nver]
direction_to_lights = direction_to_lights/direction_to_lights_n[:, :, np.newaxis]
normals_dot_lights = normals[np.newaxis, :, :]*direction_to_lights # [nlight, nver, 3]
normals_dot_lights = np.sum(normals_dot_lights, axis = 2) # [nlight, nver]
diffuse_output = colors[np.newaxis, :, :]*normals_dot_lights[:, :, np.newaxis]*light_intensities[:, np.newaxis, :]
diffuse_output = np.sum(diffuse_output, axis = 0) # [nver, 3]
# specular
# h = (v + l)/(|v + l|) bisector
# Ls = ks*(I/r^2)max(0, nxh)^p
# increasing p narrows the reflectionlob
lit_colors = diffuse_output # only diffuse part here.
lit_colors = np.minimum(np.maximum(lit_colors, 0), 1)
return lit_colors
## TODO. estimate light(sh coeff)
## -------------------------------- estimate. can not use now.
def fit_light(image, vertices, colors, triangles, vis_ind, lamb = 10, max_iter = 3):
[h, w, c] = image.shape
# surface normal
norm = get_normal(vertices, triangles)
nver = vertices.shape[1]
# vertices --> corresponding image pixel
pt2d = vertices[:2, :]
pt2d[0,:] = np.minimum(np.maximum(pt2d[0,:], 0), w - 1)
pt2d[1,:] = np.minimum(np.maximum(pt2d[1,:], 0), h - 1)
pt2d = np.round(pt2d).astype(np.int32) # 2 x nver
image_pixel = image[pt2d[1,:], pt2d[0,:], :] # nver x 3
image_pixel = image_pixel.T # 3 x nver
# vertices --> corresponding mean texture pixel with illumination
# Spherical Harmonic Basis
harmonic_dim = 9
nx = norm[0,:];
ny = norm[1,:];
nz = norm[2,:];
harmonic = np.zeros((nver, harmonic_dim))
pi = np.pi
harmonic[:,0] = np.sqrt(1/(4*pi)) * np.ones((nver,));
harmonic[:,1] = np.sqrt(3/(4*pi)) * nx;
harmonic[:,2] = np.sqrt(3/(4*pi)) * ny;
harmonic[:,3] = np.sqrt(3/(4*pi)) * nz;
harmonic[:,4] = 1/2. * np.sqrt(3/(4*pi)) * (2*nz**2 - nx**2 - ny**2);
harmonic[:,5] = 3 * np.sqrt(5/(12*pi)) * (ny*nz);
harmonic[:,6] = 3 * np.sqrt(5/(12*pi)) * (nx*nz);
harmonic[:,7] = 3 * np.sqrt(5/(12*pi)) * (nx*ny);
harmonic[:,8] = 3/2. * np.sqrt(5/(12*pi)) * (nx*nx - ny*ny);
'''
I' = sum(albedo * lj * hj) j = 0:9 (albedo = tex)
set A = albedo*h (n x 9)
alpha = lj (9 x 1)
Y = I (n x 1)
Y' = A.dot(alpha)
opt function:
||Y - A*alpha|| + lambda*(alpha'*alpha)
result:
A'*(Y - A*alpha) + lambda*alpha = 0
==>
(A'*A*alpha - lambda)*alpha = A'*Y
left: 9 x 9
right: 9 x 1
'''
n_vis_ind = len(vis_ind)
n = n_vis_ind*c
Y = np.zeros((n, 1))
A = np.zeros((n, 9))
light = np.zeros((3, 1))
for k in range(c):
Y[k*n_vis_ind:(k+1)*n_vis_ind, :] = image_pixel[k, vis_ind][:, np.newaxis]
A[k*n_vis_ind:(k+1)*n_vis_ind, :] = texture[k, vis_ind][:, np.newaxis] * harmonic[vis_ind, :]
Ac = texture[k, vis_ind][:, np.newaxis]
Yc = image_pixel[k, vis_ind][:, np.newaxis]
light[k] = (Ac.T.dot(Yc))/(Ac.T.dot(Ac))
for i in range(max_iter):
Yc = Y.copy()
for k in range(c):
Yc[k*n_vis_ind:(k+1)*n_vis_ind, :] /= light[k]
# update alpha
equation_left = np.dot(A.T, A) + lamb*np.eye(harmonic_dim); # why + ?
equation_right = np.dot(A.T, Yc)
alpha = np.dot(np.linalg.inv(equation_left), equation_right)
# update light
for k in range(c):
Ac = A[k*n_vis_ind:(k+1)*n_vis_ind, :].dot(alpha)
Yc = Y[k*n_vis_ind:(k+1)*n_vis_ind, :]
light[k] = (Ac.T.dot(Yc))/(Ac.T.dot(Ac))
appearance = np.zeros_like(texture)
for k in range(c):
tmp = np.dot(harmonic*texture[k, :][:, np.newaxis], alpha*light[k])
appearance[k,:] = tmp.T
appearance = np.minimum(np.maximum(appearance, 0), 1)
return appearance

View File

@ -1,135 +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: 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
from .cython 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

View File

@ -1,383 +0,0 @@
'''
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
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

View File

@ -1,24 +0,0 @@
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?

View File

@ -1,10 +0,0 @@
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

View File

@ -1,170 +0,0 @@
''' 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)

View File

@ -1,215 +0,0 @@
'''
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
from __future__ import division
from __future__ import print_function
import numpy as np
def get_normal(vertices, triangles):
''' calculate normal direction in each vertex
Args:
vertices: [nver, 3]
triangles: [ntri, 3]
Returns:
normal: [nver, 3]
'''
pt0 = vertices[triangles[:, 0], :] # [ntri, 3]
pt1 = vertices[triangles[:, 1], :] # [ntri, 3]
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, :]
# normalize to unit length
mag = np.sum(normal**2, 1) # [nver]
zero_ind = (mag == 0)
mag[zero_ind] = 1;
normal[zero_ind, 0] = np.ones((np.sum(zero_ind)))
normal = normal/np.sqrt(mag[:,np.newaxis])
return normal
# TODO: test
def add_light_sh(vertices, triangles, colors, sh_coeff):
'''
In 3d face, usually assume:
1. The surface of face is Lambertian(reflect only the low frequencies of lighting)
2. Lighting can be an arbitrary combination of point sources
--> can be expressed in terms of spherical harmonics(omit the lighting coefficients)
I = albedo * (sh(n) x sh_coeff)
albedo: n x 1
sh_coeff: 9 x 1
Y(n) = (1, n_x, n_y, n_z, n_xn_y, n_xn_z, n_yn_z, n_x^2 - n_y^2, 3n_z^2 - 1)': n x 9
# Y(n) = (1, n_x, n_y, n_z)': n x 4
Args:
vertices: [nver, 3]
triangles: [ntri, 3]
colors: [nver, 3] albedo
sh_coeff: [9, 1] spherical harmonics coefficients
Returns:
lit_colors: [nver, 3]
'''
assert vertices.shape[0] == colors.shape[0]
nver = vertices.shape[0]
normal = get_normal(vertices, triangles) # [nver, 3]
sh = np.array((np.ones(nver), n[:,0], n[:,1], n[:,2], n[:,0]*n[:,1], n[:,0]*n[:,2], n[:,1]*n[:,2], n[:,0]**2 - n[:,1]**2, 3*(n[:,2]**2) - 1)) # [nver, 9]
ref = sh.dot(sh_coeff) #[nver, 1]
lit_colors = colors*ref
return lit_colors
def add_light(vertices, triangles, colors, light_positions = 0, light_intensities = 0):
''' Gouraud shading. add point lights.
In 3d face, usually assume:
1. The surface of face is Lambertian(reflect only the low frequencies of lighting)
2. Lighting can be an arbitrary combination of point sources
3. No specular (unless skin is oil, 23333)
Ref: https://cs184.eecs.berkeley.edu/lecture/pipeline
Args:
vertices: [nver, 3]
triangles: [ntri, 3]
light_positions: [nlight, 3]
light_intensities: [nlight, 3]
Returns:
lit_colors: [nver, 3]
'''
nver = vertices.shape[0]
normals = get_normal(vertices, triangles) # [nver, 3]
# ambient
# La = ka*Ia
# diffuse
# Ld = kd*(I/r^2)max(0, nxl)
direction_to_lights = vertices[np.newaxis, :, :] - light_positions[:, np.newaxis, :] # [nlight, nver, 3]
direction_to_lights_n = np.sqrt(np.sum(direction_to_lights**2, axis = 2)) # [nlight, nver]
direction_to_lights = direction_to_lights/direction_to_lights_n[:, :, np.newaxis]
normals_dot_lights = normals[np.newaxis, :, :]*direction_to_lights # [nlight, nver, 3]
normals_dot_lights = np.sum(normals_dot_lights, axis = 2) # [nlight, nver]
diffuse_output = colors[np.newaxis, :, :]*normals_dot_lights[:, :, np.newaxis]*light_intensities[:, np.newaxis, :]
diffuse_output = np.sum(diffuse_output, axis = 0) # [nver, 3]
# specular
# h = (v + l)/(|v + l|) bisector
# Ls = ks*(I/r^2)max(0, nxh)^p
# increasing p narrows the reflectionlob
lit_colors = diffuse_output # only diffuse part here.
lit_colors = np.minimum(np.maximum(lit_colors, 0), 1)
return lit_colors
## TODO. estimate light(sh coeff)
## -------------------------------- estimate. can not use now.
def fit_light(image, vertices, colors, triangles, vis_ind, lamb = 10, max_iter = 3):
[h, w, c] = image.shape
# surface normal
norm = get_normal(vertices, triangles)
nver = vertices.shape[1]
# vertices --> corresponding image pixel
pt2d = vertices[:2, :]
pt2d[0,:] = np.minimum(np.maximum(pt2d[0,:], 0), w - 1)
pt2d[1,:] = np.minimum(np.maximum(pt2d[1,:], 0), h - 1)
pt2d = np.round(pt2d).astype(np.int32) # 2 x nver
image_pixel = image[pt2d[1,:], pt2d[0,:], :] # nver x 3
image_pixel = image_pixel.T # 3 x nver
# vertices --> corresponding mean texture pixel with illumination
# Spherical Harmonic Basis
harmonic_dim = 9
nx = norm[0,:];
ny = norm[1,:];
nz = norm[2,:];
harmonic = np.zeros((nver, harmonic_dim))
pi = np.pi
harmonic[:,0] = np.sqrt(1/(4*pi)) * np.ones((nver,));
harmonic[:,1] = np.sqrt(3/(4*pi)) * nx;
harmonic[:,2] = np.sqrt(3/(4*pi)) * ny;
harmonic[:,3] = np.sqrt(3/(4*pi)) * nz;
harmonic[:,4] = 1/2. * np.sqrt(3/(4*pi)) * (2*nz**2 - nx**2 - ny**2);
harmonic[:,5] = 3 * np.sqrt(5/(12*pi)) * (ny*nz);
harmonic[:,6] = 3 * np.sqrt(5/(12*pi)) * (nx*nz);
harmonic[:,7] = 3 * np.sqrt(5/(12*pi)) * (nx*ny);
harmonic[:,8] = 3/2. * np.sqrt(5/(12*pi)) * (nx*nx - ny*ny);
'''
I' = sum(albedo * lj * hj) j = 0:9 (albedo = tex)
set A = albedo*h (n x 9)
alpha = lj (9 x 1)
Y = I (n x 1)
Y' = A.dot(alpha)
opt function:
||Y - A*alpha|| + lambda*(alpha'*alpha)
result:
A'*(Y - A*alpha) + lambda*alpha = 0
==>
(A'*A*alpha - lambda)*alpha = A'*Y
left: 9 x 9
right: 9 x 1
'''
n_vis_ind = len(vis_ind)
n = n_vis_ind*c
Y = np.zeros((n, 1))
A = np.zeros((n, 9))
light = np.zeros((3, 1))
for k in range(c):
Y[k*n_vis_ind:(k+1)*n_vis_ind, :] = image_pixel[k, vis_ind][:, np.newaxis]
A[k*n_vis_ind:(k+1)*n_vis_ind, :] = texture[k, vis_ind][:, np.newaxis] * harmonic[vis_ind, :]
Ac = texture[k, vis_ind][:, np.newaxis]
Yc = image_pixel[k, vis_ind][:, np.newaxis]
light[k] = (Ac.T.dot(Yc))/(Ac.T.dot(Ac))
for i in range(max_iter):
Yc = Y.copy()
for k in range(c):
Yc[k*n_vis_ind:(k+1)*n_vis_ind, :] /= light[k]
# update alpha
equation_left = np.dot(A.T, A) + lamb*np.eye(harmonic_dim); # why + ?
equation_right = np.dot(A.T, Yc)
alpha = np.dot(np.linalg.inv(equation_left), equation_right)
# update light
for k in range(c):
Ac = A[k*n_vis_ind:(k+1)*n_vis_ind, :].dot(alpha)
Yc = Y[k*n_vis_ind:(k+1)*n_vis_ind, :]
light[k] = (Ac.T.dot(Yc))/(Ac.T.dot(Ac))
appearance = np.zeros_like(texture)
for k in range(c):
tmp = np.dot(harmonic*texture[k, :][:, np.newaxis], alpha*light[k])
appearance[k,:] = tmp.T
appearance = np.minimum(np.maximum(appearance, 0), 1)
return appearance

View File

@ -1,287 +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)
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<umin or vmax<vmin:
continue
for u in range(umin, umax+1):
for v in range(vmin, vmax+1):
if not isPointInTri([u,v], vertices[tri, :2]):
continue
w0, w1, w2 = get_point_weight([u, v], vertices[tri, :2]) # barycentric weight
point_depth = w0*vertices[tri[0], 2] + w1*vertices[tri[1], 2] + w2*vertices[tri[2], 2]
if point_depth > 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<umin or vmax<vmin:
continue
for u in range(umin, umax+1):
for v in range(vmin, vmax+1):
if not isPointInTri([u,v], vertices[tri, :2]):
continue
w0, w1, w2 = get_point_weight([u, v], vertices[tri, :2])
point_depth = w0*vertices[tri[0], 2] + w1*vertices[tri[1], 2] + w2*vertices[tri[2], 2]
if point_depth > 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<umin or vmax<vmin:
continue
for u in range(umin, umax+1):
for v in range(vmin, vmax+1):
if not isPointInTri([u,v], vertices[tri, :2]):
continue
w0, w1, w2 = get_point_weight([u, v], vertices[tri, :2])
point_depth = w0*vertices[tri[0], 2] + w1*vertices[tri[1], 2] + w2*vertices[tri[2], 2]
if point_depth > 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

View File

@ -1,385 +0,0 @@
'''
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

View File

@ -1,24 +0,0 @@
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?

View File

@ -1,7 +0,0 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from .. import mesh
from .morphabel_model import MorphabelModel
from . import load

View File

@ -1,272 +0,0 @@
'''
Estimating parameters about vertices: shape para, exp para, pose para(s, R, t)
'''
import numpy as np
from .. import mesh
''' TODO: a clear document.
Given: image_points, 3D Model, Camera Matrix(s, R, t2d)
Estimate: shape parameters, expression parameters
Inference:
projected_vertices = s*P*R(mu + shape + exp) + t2d --> image_points
s*P*R*shape + s*P*R(mu + exp) + t2d --> image_poitns
# Define:
X = vertices
x_hat = projected_vertices
x = image_points
A = s*P*R
b = s*P*R(mu + exp) + t2d
==>
x_hat = A*shape + b (2 x n)
A*shape (2 x n)
shape = reshape(shapePC * sp) (3 x n)
shapePC*sp : (3n x 1)
* flatten:
x_hat_flatten = A*shape + b_flatten (2n x 1)
A*shape (2n x 1)
--> A*shapePC (2n x 199) sp: 199 x 1
# Define:
pc_2d = A* reshape(shapePC)
pc_2d_flatten = flatten(pc_2d) (2n x 199)
=====>
x_hat_flatten = pc_2d_flatten * sp + b_flatten ---> x_flatten (2n x 1)
Goals:
(ignore flatten, pc_2d-->pc)
min E = || x_hat - x || + lambda*sum(sp/sigma)^2
= || pc * sp + b - x || + lambda*sum(sp/sigma)^2
Solve:
d(E)/d(sp) = 0
2 * pc' * (pc * sp + b - x) + 2 * lambda * sp / (sigma' * sigma) = 0
Get:
(pc' * pc + lambda / (sigma'* sigma)) * sp = pc' * (x - b)
'''
def estimate_shape(x, shapeMU, shapePC, shapeEV, expression, s, R, t2d, lamb = 3000):
'''
Args:
x: (2, n). image points (to be fitted)
shapeMU: (3n, 1)
shapePC: (3n, n_sp)
shapeEV: (n_sp, 1)
expression: (3, n)
s: scale
R: (3, 3). rotation matrix
t2d: (2,). 2d translation
lambda: regulation coefficient
Returns:
shape_para: (n_sp, 1) shape parameters(coefficients)
'''
x = x.copy()
assert(shapeMU.shape[0] == shapePC.shape[0])
assert(shapeMU.shape[0] == x.shape[1]*3)
dof = shapePC.shape[1]
n = x.shape[1]
sigma = shapeEV
t2d = np.array(t2d)
P = np.array([[1, 0, 0], [0, 1, 0]], dtype = np.float32)
A = s*P.dot(R)
# --- calc pc
pc_3d = np.resize(shapePC.T, [dof, n, 3]) # 199 x n x 3
pc_3d = np.reshape(pc_3d, [dof*n, 3])
pc_2d = pc_3d.dot(A.T.copy()) # 199 x n x 2
pc = np.reshape(pc_2d, [dof, -1]).T # 2n x 199
# --- calc b
# shapeMU
mu_3d = np.resize(shapeMU, [n, 3]).T # 3 x n
# expression
exp_3d = expression
#
b = A.dot(mu_3d + exp_3d) + np.tile(t2d[:, np.newaxis], [1, n]) # 2 x n
b = np.reshape(b.T, [-1, 1]) # 2n x 1
# --- solve
equation_left = np.dot(pc.T, pc) + lamb * np.diagflat(1/sigma**2)
x = np.reshape(x.T, [-1, 1])
equation_right = np.dot(pc.T, x - b)
shape_para = np.dot(np.linalg.inv(equation_left), equation_right)
return shape_para
def estimate_expression(x, shapeMU, expPC, expEV, shape, s, R, t2d, lamb = 2000):
'''
Args:
x: (2, n). image points (to be fitted)
shapeMU: (3n, 1)
expPC: (3n, n_ep)
expEV: (n_ep, 1)
shape: (3, n)
s: scale
R: (3, 3). rotation matrix
t2d: (2,). 2d translation
lambda: regulation coefficient
Returns:
exp_para: (n_ep, 1) shape parameters(coefficients)
'''
x = x.copy()
assert(shapeMU.shape[0] == expPC.shape[0])
assert(shapeMU.shape[0] == x.shape[1]*3)
dof = expPC.shape[1]
n = x.shape[1]
sigma = expEV
t2d = np.array(t2d)
P = np.array([[1, 0, 0], [0, 1, 0]], dtype = np.float32)
A = s*P.dot(R)
# --- calc pc
pc_3d = np.resize(expPC.T, [dof, n, 3])
pc_3d = np.reshape(pc_3d, [dof*n, 3])
pc_2d = pc_3d.dot(A.T)
pc = np.reshape(pc_2d, [dof, -1]).T # 2n x 29
# --- calc b
# shapeMU
mu_3d = np.resize(shapeMU, [n, 3]).T # 3 x n
# expression
shape_3d = shape
#
b = A.dot(mu_3d + shape_3d) + np.tile(t2d[:, np.newaxis], [1, n]) # 2 x n
b = np.reshape(b.T, [-1, 1]) # 2n x 1
# --- solve
equation_left = np.dot(pc.T, pc) + lamb * np.diagflat(1/sigma**2)
x = np.reshape(x.T, [-1, 1])
equation_right = np.dot(pc.T, x - b)
exp_para = np.dot(np.linalg.inv(equation_left), equation_right)
return exp_para
# ---------------- fit
def fit_points(x, X_ind, model, n_sp, n_ep, max_iter = 4):
'''
Args:
x: (n, 2) image points
X_ind: (n,) corresponding Model vertex indices
model: 3DMM
max_iter: iteration
Returns:
sp: (n_sp, 1). shape parameters
ep: (n_ep, 1). exp parameters
s, R, t
'''
x = x.copy().T
#-- init
sp = np.zeros((n_sp, 1), dtype = np.float32)
ep = np.zeros((n_ep, 1), dtype = np.float32)
#-------------------- estimate
X_ind_all = np.tile(X_ind[np.newaxis, :], [3, 1])*3
X_ind_all[1, :] += 1
X_ind_all[2, :] += 2
valid_ind = X_ind_all.flatten('F')
shapeMU = model['shapeMU'][valid_ind, :]
shapePC = model['shapePC'][valid_ind, :n_sp]
expPC = model['expPC'][valid_ind, :n_ep]
for i in range(max_iter):
X = shapeMU + shapePC.dot(sp) + expPC.dot(ep)
X = np.reshape(X, [int(len(X)/3), 3]).T
#----- estimate pose
P = mesh.transform.estimate_affine_matrix_3d22d(X.T, x.T)
s, R, t = mesh.transform.P2sRt(P)
rx, ry, rz = mesh.transform.matrix2angle(R)
#print('Iter:{}; estimated pose: s {}, rx {}, ry {}, rz {}, t1 {}, t2 {}'.format(i, s, rx, ry, rz, t[0], t[1]))
#----- estimate shape
# expression
shape = shapePC.dot(sp)
shape = np.reshape(shape, [int(len(shape)/3), 3]).T
ep = estimate_expression(x, shapeMU, expPC, model['expEV'][:n_ep,:], shape, s, R, t[:2], lamb = 20)
# shape
expression = expPC.dot(ep)
expression = np.reshape(expression, [int(len(expression)/3), 3]).T
if i == 0 :
sp = estimate_shape(x, shapeMU, shapePC, model['shapeEV'][:n_sp,:], expression, s, R, t[:2], lamb = 40)
return sp, ep, s, R, t
# ---------------- fitting process
def fit_points_for_show(x, X_ind, model, n_sp, n_ep, max_iter = 4):
'''
Args:
x: (n, 2) image points
X_ind: (n,) corresponding Model vertex indices
model: 3DMM
max_iter: iteration
Returns:
sp: (n_sp, 1). shape parameters
ep: (n_ep, 1). exp parameters
s, R, t
'''
x = x.copy().T
#-- init
sp = np.zeros((n_sp, 1), dtype = np.float32)
ep = np.zeros((n_ep, 1), dtype = np.float32)
#-------------------- estimate
X_ind_all = np.tile(X_ind[np.newaxis, :], [3, 1])*3
X_ind_all[1, :] += 1
X_ind_all[2, :] += 2
valid_ind = X_ind_all.flatten('F')
shapeMU = model['shapeMU'][valid_ind, :]
shapePC = model['shapePC'][valid_ind, :n_sp]
expPC = model['expPC'][valid_ind, :n_ep]
s = 4e-04
R = mesh.transform.angle2matrix([0, 0, 0])
t = [0, 0, 0]
lsp = []; lep = []; ls = []; lR = []; lt = []
for i in range(max_iter):
X = shapeMU + shapePC.dot(sp) + expPC.dot(ep)
X = np.reshape(X, [int(len(X)/3), 3]).T
lsp.append(sp); lep.append(ep); ls.append(s), lR.append(R), lt.append(t)
#----- estimate pose
P = mesh.transform.estimate_affine_matrix_3d22d(X.T, x.T)
s, R, t = mesh.transform.P2sRt(P)
lsp.append(sp); lep.append(ep); ls.append(s), lR.append(R), lt.append(t)
#----- estimate shape
# expression
shape = shapePC.dot(sp)
shape = np.reshape(shape, [int(len(shape)/3), 3]).T
ep = estimate_expression(x, shapeMU, expPC, model['expEV'][:n_ep,:], shape, s, R, t[:2], lamb = 20)
lsp.append(sp); lep.append(ep); ls.append(s), lR.append(R), lt.append(t)
# shape
expression = expPC.dot(ep)
expression = np.reshape(expression, [int(len(expression)/3), 3]).T
sp = estimate_shape(x, shapeMU, shapePC, model['shapeEV'][:n_sp,:], expression, s, R, t[:2], lamb = 40)
# print('ls', ls)
# print('lR', lR)
return np.array(lsp), np.array(lep), np.array(ls), np.array(lR), np.array(lt)

View File

@ -1,110 +0,0 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import numpy as np
import scipy.io as sio
### --------------------------------- load BFM data
def load_BFM(model_path):
''' load BFM 3DMM model
Args:
model_path: path to BFM model.
Returns:
model: (nver = 53215, ntri = 105840). nver: number of vertices. ntri: number of triangles.
'shapeMU': [3*nver, 1]
'shapePC': [3*nver, 199]
'shapeEV': [199, 1]
'expMU': [3*nver, 1]
'expPC': [3*nver, 29]
'expEV': [29, 1]
'texMU': [3*nver, 1]
'texPC': [3*nver, 199]
'texEV': [199, 1]
'tri': [ntri, 3] (start from 1, should sub 1 in python and c++)
'tri_lip': [114, 3] (start from 1, as a supplement to lip triangles)
'kpt_ind': [68,] (start from 1)
PS:
You can change codes according to your own saved data.
Just make sure the model has corresponding attributes.
'''
C = sio.loadmat(model_path)
model = C['model']
model = model[0,0]
# change dtype from double(np.float64) to np.float32,
# since big matrix process(espetially matrix dot) is too slow in python.
model['shapeMU'] = (model['shapeMU'] + model['expMU']).astype(np.float32)
model['shapePC'] = model['shapePC'].astype(np.float32)
model['shapeEV'] = model['shapeEV'].astype(np.float32)
model['expEV'] = model['expEV'].astype(np.float32)
model['expPC'] = model['expPC'].astype(np.float32)
# matlab start with 1. change to 0 in python.
model['tri'] = model['tri'].T.copy(order = 'C').astype(np.int32) - 1
model['tri_lip'] = model['tri_lip'].T.copy(order = 'C').astype(np.int32) - 1
# kpt ind
model['kpt_ind'] = (np.squeeze(model['kpt_ind']) - 1).astype(np.int32)
return model
def load_BFM_info(path = 'BFM_info.mat'):
''' load 3DMM model extra information
Args:
path: path to BFM info.
Returns:
model_info:
'symlist': 2 x 26720
'symlist_tri': 2 x 52937
'segbin': 4 x n (0: nose, 1: eye, 2: lip, 3: cheek)
'segbin_tri': 4 x ntri
'face_contour': 1 x 28
'face_contour_line': 1 x 512
'face_contour_front': 1 x 28
'face_contour_front_line': 1 x 512
'nose_hole': 1 x 142
'nose_hole_right': 1 x 71
'nose_hole_left': 1 x 71
'parallel': 17 x 1 cell
'parallel_face_contour': 28 x 1 cell
'uv_coords': n x 2
'''
C = sio.loadmat(path)
model_info = C['model_info']
model_info = model_info[0,0]
return model_info
def load_uv_coords(path = 'BFM_UV.mat'):
''' load uv coords of BFM
Args:
path: path to data.
Returns:
uv_coords: [nver, 2]. range: 0-1
'''
C = sio.loadmat(path)
uv_coords = C['UV'].copy(order = 'C')
return uv_coords
def load_pncc_code(path = 'pncc_code.mat'):
''' load pncc code of BFM
PNCC code: Defined in 'Face Alignment Across Large Poses: A 3D Solution Xiangyu'
download at http://www.cbsr.ia.ac.cn/users/xiangyuzhu/projects/3DDFA/main.htm.
Args:
path: path to data.
Returns:
pncc_code: [nver, 3]
'''
C = sio.loadmat(path)
pncc_code = C['vertex_code'].T
return pncc_code
##
def get_organ_ind(model_info):
''' get nose, eye, lip index
'''
valid_bin = model_info['segbin'].astype(bool)
organ_ind = np.nonzero(valid_bin[0,:])[0]
for i in range(1, valid_bin.shape[0] - 1):
organ_ind = np.union1d(organ_ind, np.nonzero(valid_bin[i,:])[0])
return organ_ind.astype(np.int32)

View File

@ -1,143 +0,0 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import numpy as np
import scipy.io as sio
from .. import mesh
from . import fit
from . import load
class MorphabelModel(object):
"""docstring for MorphabelModel
model: nver: number of vertices. ntri: number of triangles. *: must have. ~: can generate ones array for place holder.
'shapeMU': [3*nver, 1]. *
'shapePC': [3*nver, n_shape_para]. *
'shapeEV': [n_shape_para, 1]. ~
'expMU': [3*nver, 1]. ~
'expPC': [3*nver, n_exp_para]. ~
'expEV': [n_exp_para, 1]. ~
'texMU': [3*nver, 1]. ~
'texPC': [3*nver, n_tex_para]. ~
'texEV': [n_tex_para, 1]. ~
'tri': [ntri, 3] (start from 1, should sub 1 in python and c++). *
'tri_lip': [114, 3] (start from 1, as a supplement to lip triangles). ~
'kpt_ind': [68,] (start from 1). ~
"""
def __init__(self, model_path, model_type = 'BFM'):
super( MorphabelModel, self).__init__()
if model_type=='BFM':
self.model = load.load_BFM(model_path)
else:
print('sorry, not support other 3DMM model now')
exit()
# fixed attributes
self.nver = self.model['shapePC'].shape[0]/3
self.ntri = self.model['tri'].shape[0]
self.n_shape_para = self.model['shapePC'].shape[1]
self.n_exp_para = self.model['expPC'].shape[1]
self.n_tex_para = self.model['texMU'].shape[1]
self.kpt_ind = self.model['kpt_ind']
self.triangles = self.model['tri']
self.full_triangles = np.vstack((self.model['tri'], self.model['tri_lip']))
# ------------------------------------- shape: represented with mesh(vertices & triangles(fixed))
def get_shape_para(self, type = 'random'):
if type == 'zero':
sp = np.random.zeros((self.n_shape_para, 1))
elif type == 'random':
sp = np.random.rand(self.n_shape_para, 1)*1e04
return sp
def get_exp_para(self, type = 'random'):
if type == 'zero':
ep = np.zeros((self.n_exp_para, 1))
elif type == 'random':
ep = -1.5 + 3*np.random.random([self.n_exp_para, 1])
ep[6:, 0] = 0
return ep
def generate_vertices(self, shape_para, exp_para):
'''
Args:
shape_para: (n_shape_para, 1)
exp_para: (n_exp_para, 1)
Returns:
vertices: (nver, 3)
'''
vertices = self.model['shapeMU'] + self.model['shapePC'].dot(shape_para) + self.model['expPC'].dot(exp_para)
vertices = np.reshape(vertices, [int(3), int(len(vertices)/3)], 'F').T
return vertices
# -------------------------------------- texture: here represented with rgb value(colors) in vertices.
def get_tex_para(self, type = 'random'):
if type == 'zero':
tp = np.zeros((self.n_tex_para, 1))
elif type == 'random':
tp = np.random.rand(self.n_tex_para, 1)
return tp
def generate_colors(self, tex_para):
'''
Args:
tex_para: (n_tex_para, 1)
Returns:
colors: (nver, 3)
'''
colors = self.model['texMU'] + self.model['texPC'].dot(tex_para*self.model['texEV'])
colors = np.reshape(colors, [int(3), int(len(colors)/3)], 'F').T/255.
return colors
# ------------------------------------------- transformation
# ------------- transform
def rotate(self, vertices, angles):
''' rotate face
Args:
vertices: [nver, 3]
angles: [3] x, y, z rotation angle(degree)
x: pitch. positive for looking down
y: yaw. positive for looking left
z: roll. positive for tilting head right
Returns:
vertices: rotated vertices
'''
return mesh.transform.rotate(vertices, angles)
def transform(self, vertices, s, angles, t3d):
R = mesh.transform.angle2matrix(angles)
return mesh.transform.similarity_transform(vertices, s, R, t3d)
def transform_3ddfa(self, vertices, s, angles, t3d): # only used for processing 300W_LP data
R = mesh.transform.angle2matrix_3ddfa(angles)
return mesh.transform.similarity_transform(vertices, s, R, t3d)
# --------------------------------------------------- fitting
def fit(self, x, X_ind, max_iter = 4, isShow = False):
''' fit 3dmm & pose parameters
Args:
x: (n, 2) image points
X_ind: (n,) corresponding Model vertex indices
max_iter: iteration
isShow: whether to reserve middle results for show
Returns:
fitted_sp: (n_sp, 1). shape parameters
fitted_ep: (n_ep, 1). exp parameters
s, angles, t
'''
if isShow:
fitted_sp, fitted_ep, s, R, t = fit.fit_points_for_show(x, X_ind, self.model, n_sp = self.n_shape_para, n_ep = self.n_exp_para, max_iter = max_iter)
angles = np.zeros((R.shape[0], 3))
for i in range(R.shape[0]):
angles[i] = mesh.transform.matrix2angle(R[i])
else:
fitted_sp, fitted_ep, s, R, t = fit.fit_points(x, X_ind, self.model, n_sp = self.n_shape_para, n_ep = self.n_exp_para, max_iter = max_iter)
angles = mesh.transform.matrix2angle(R)
return fitted_sp, fitted_ep, s, angles, t