本文主要是对opendr部分源码阅读。
chumpy代码地址:https://github.com/mattloper/chumpy
chumpy使用说明:http://files.is.tue.mpg.de/black/papers/chumpy_tutorial.pdf
阅读SMPL等三维重建代码时,经常能看到引用opendr,chumpy第三方库(SMPL,OpenDR,和chumpy出自同个作者,opendr发表于2014,smpl发表于2015)关于opendr的源码阅读,网上没有相关信息,所以在这尽最大所能对部分代码进行说明。不得不说,作者对opendr中变量命名真是蛋疼啊。。。本文首先介绍chumpy,随后介绍opendr里的最最最最常用的三个文件camer.py。
1、chumpy
作者说为了方便计算微分,自己造了一个轮子,命名为chumpy。emmm可惜使用的人不是很多,github上star为52,估计都是做三维重建“被迫”研究的,网上关于chumpy的资料也很少,不过好在作者提供了一份说明文档(链接见文首),不用从头看源代码。
我简单的对文档里几个重要的函数或类进行说明。
1.1 chumpy objects
我们可以自定义自己的chumpy类,加入sin函数不存在的话,可以创建自己的sin类。
from chumpy import Ch import scipy.sparse as sp class Sin(Ch): #继承了Ch类 dterms = (’x’,) #定义所有可微参数变量 def compute_r(self): #r应该就是result的缩写。这个是继承Ch后要重写的函数。 return np.sin(self.x.r) #self.x.r 中的x是dterm里的参数。 def compute_dr_wrt(self, wrt): #计算r关于wrt参数的导数。这个是继承Ch后也要重写的函数 if wrt is self.x: #判断输入参数是不是上述列的可微参数集合。 values = np.cos(self.x.r) return sp.diags([values.ravel()], [0]) # 使用这个类的实例 x1 = Ch(10) result = Sin(x1) # or "result = Sin(x=x1)" print result print result.dr_wrt(x1)
作者还提到我们可能会有一些函数,它的参数不可微,此时可以用term参数。例如,The following example allows the specification of values and indices to select, and returns the unraveled, selected input:(关于csc的用法,网上有很多教程,可以参看 https://blog.csdn.net/sinat_33741547/article/details/79878547 文章)
from chumpy import Ch import scipy.sparse as sp class Select(Ch): dterms = (’x’,) #可微的参数 terms = (’indices’,) #不可微的参数 def compute_r(self): return self.x.r.ravel()[self.indices] # 平铺一行矩阵,输出选定的值 def compute_dr_wrt(self, wrt): # np.arrang 输入x长度的递增序列,比如[1,1.5,2,2.5,3,...,10] # lself.indices = [1 , 2 , 5] # np.ones = [1 , 1 , 1] return sp.csc_matrix(((np.arange(self.x.size), self.indices), np.ones(self.indices.size))) # 返回的结果应该是 # 使用类 x1 = ch.arange(20)/2. selected = Select(x=x1, indices=[1,2,5]) print selected # prints .5, 1, 2.5 selected.indices = [10,20] #改变索引 print selected # prints 5, 10 print selected.dr_wrt(selected.indices) # prints None print selected.dr_wrt(selected.x) # returns sparse matrix
1.2 caching
chumpy 可以使用“depend_on”装饰器,它公开一个属性,该属性仅在其依赖项的一个或多个发生更改时根据需要重新计算。例如上例中,如果频繁运行np.ones,则会导致重复的内存申请,此时我们可以这样更改:
from chumpy import Ch import scipy.sparse as sp class Select(Ch): dterms = (’x’,) terms = (’indices’,) def compute_r(self): return self.x.r.ravel()[self.indices] @depends_on(’indices’) def ones(self): return np.ones(self.indices.size) def compute_dr_wrt(self, wrt): return sp.csc_matrix(((np.arange(self.x.size), self.indices), self.ones))
另一个是使用“on_changed”方法,它提供了不同的风格。 假设我们有一个隐藏值,只想在更改外部术语时重新计算,并且只在调用compute_r或compute_dr_wrt之前。 方法如下:
from chumpy import Ch import scipy.sparse as sp class Select(Ch): dterms = (’x’,) terms = (’indices’,) def compute_r(self): return self.x.r.ravel()[self.indices] def on_changed(self, which): if ’indices’ in which: # which is always a nonempty list of strings self.ones = np.ones(self.indices.size) def compute_dr_wrt(self, wrt): return sp.csc_matrix(((np.arange(self.x.size), self.indices), self.ones))
只当在compute_r or compute_dr_wrt被调用时,on_changed函数才开始执行。
2、opendr代码理解
opendr就比较惨了,没文档,注释少,网上的信息也少,需要我们阅读源码。其中最常用的三个文件:camera.py lighting.py 和 render.py。以下逐一介绍和注释。
2.1、camera.py
需要的储备知识:
刚体变换:刚体变换用于刚性物体的变换,只改变物体的方向和位置,不改变形状。可以将刚体矩阵X写成一个平移矩阵和一个旋转矩阵的级联。
相机内外参数
/usr/bin/env python # encoding: utf-8 """ Author(s): Matthew Loper See LICENCE.txt for licensing and contact information. """ __all__ = ['ProjectPoints3D', 'ProjectPoints', 'RigidTransform'] import chumpy as ch from chumpy import depends_on, Ch from cvwrap import cv2 import numpy as np import scipy.sparse as sp from chumpy.utils import row, col from geometry import Rodrigues def RigidTransformSlow(**kwargs): # 刚体变换函数,速度慢 # Returns a Ch object with dterms 'v', 'rt', and 't' #定义可微变量 v是刚体,rt是旋转向量,t是平移量 result = Ch(lambda v, rt, t : v.dot(Rodrigues(rt=rt)) + t) # 刚体变换公式,同时夹杂Rodrigues公式 if len(kwargs) > 0: result.set(**kwargs) return result class RigidTransform(Ch): # 刚体变换类 dterms = 'v', 'rt', 't' #定义可微变量 v是刚体,rt是旋转向量,t是平移量 def compute_r(self): return (cv2.Rodrigues(self.rt.r)[0].dot(self.v.r.T) + col(self.t.r)).T.copy() #self.rt.r就是rt的结果,看过chumpy才知道,之前一直在这晕着 def compute_dr_wrt(self, wrt): if wrt not in (self.v, self.rt, self.t): #如果输入的参数不是这些参数 return if wrt is self.t: #如果是对平移量t求微 if not hasattr(self, '_drt') or self._drt.shape[0] != self.v.r.size: #如果没有dir属性或者shape不相等 IS = np.arange(self.v.r.size) # 行数,比如[1,2,3,4,5] JS = IS % 3 # 筛选的列数 [1,2,0,1,2] 模三 指定三列 data = np.ones(len(IS)) self._drt = sp.csc_matrix((data, (IS, JS))) #生成_drt return self._drt if wrt is self.rt: #如果是对旋转向量求微 rot, rot_dr = cv2.Rodrigues(self.rt.r) #第一项为旋转矩阵或向量,第二项为雅克比矩阵,输出的偏导 rot_dr = rot_dr.reshape((3,3,3)) dr = np.einsum('abc, zc ->zba', rot_dr, self.v.r).reshape((-1,3)) return dr if wrt is self.v: rot = cv2.Rodrigues(self.rt.r)[0] # 旋转矩阵 IS = np.repeat(np.arange(self.v.r.size), 3) JS = np.repeat(np.arange(self.v.r.size).reshape((-1,3)), 3, axis=0) data = np.vstack([rot for i in range(self.v.r.size/3)]) result = sp.csc_matrix((data.ravel(), (IS.ravel(), JS.ravel()))) return result class ProjectPoints(Ch): # 给定内参数和外参数计算射影点 dterms = 'v', 'rt', 't', 'f', 'c', 'k' # v 刚体坐标 rt 旋转矩阵 t 平移向量 这些是相机外参数 # f 也不算焦距 c 也不算主点坐标 k 畸变参数 这些是相机内参数 def is_valid(self): if any([len(v.r.shape) > 1 for v in [self.rt, self.t, self.f, self.c, self.k]]): return False, 'rt, t, f, c, and k must be 1D' #必须是1维向量 if any([v.r.size != 3 for v in [self.rt, self.t]]): # 对应着三维坐标 return False, 'rt and t must have size=3' if any([v.r.size != 2 for v in [self.f, self.c]]): # fx fy cx cy return False, 'f and c must have size=2' return True, '' def compute_r(self): return self.r_and_derivatives[0].squeeze() # 返回二维平面点 #return self.get_r_and_derivatives(self.v.r, self.rt.r, self.t.r, self.f.r, self.c.r, self.k.r)[0].squeeze() def compute_dr_wrt(self, wrt): if wrt not in [self.v, self.rt, self.t, self.f, self.c, self.k]: return None j = self.r_and_derivatives[1] # 返回雅克比矩阵 if wrt is self.rt: return j[:, :3] elif wrt is self.t: return j[:, 3:6] elif wrt is self.f: return j[:, 6:8] elif wrt is self.c: return j[:, 8:10] elif wrt is self.k: return j[:, 10:10+self.k.size] elif wrt is self.v: rot = cv2.Rodrigues(self.rt.r)[0] data = np.asarray(j[:, 3:6].dot(rot), order='C').ravel() # 雅克比矩阵点乘旋转矩阵 IS = np.repeat(np.arange(self.v.r.size*2/3), 3) JS = np.asarray(np.repeat(np.arange(self.v.r.size).reshape((-1,3)), 2, axis=0), order='C').ravel() result = sp.csc_matrix((data, (IS, JS))) #暂时还不理解这行,等以后填坑 return result def unproject_points(self, uvd, camera_space=False): cam = ProjectPoints3D(**{k: getattr(self, k) for k in self.dterms if hasattr(self, k)}) try: xy_undistorted_camspace = cv2.undistortPoints(np.asarray(uvd[:,:2].reshape((1,-1,2)).copy()), np.asarray(cam.camera_mtx), cam.k.r) xyz_camera_space = np.hstack((xy_undistorted_camspace.squeeze(), col(uvd[:,2]))) xyz_camera_space[:,:2] *= col(xyz_camera_space[:,2]) # scale x,y by z if camera_space: return xyz_camera_space other_answer = xyz_camera_space - row(cam.view_mtx[:,3]) # translate result = other_answer.dot(cam.view_mtx[:,:3]) # rotate except: # slow way, probably not so good. But doesn't require cv2.undistortPoints. cam.v = np.ones_like(uvd) ch.minimize(cam - uvd, x0=[cam.v], method='dogleg', options={'disp': 0}) result = cam.v.r return result def unproject_depth_image(self, depth_image, camera_space=False): us = np.arange(depth_image.size) % depth_image.shape[1] vs = np.arange(depth_image.size) // depth_image.shape[1] ds = depth_image.ravel() uvd = ch.array(np.vstack((us.ravel(), vs.ravel(), ds.ravel())).T) xyz = self.unproject_points(uvd, camera_space=camera_space) return xyz.reshape((depth_image.shape[0], depth_image.shape[1], -1)) @depends_on('f','c') def camera_mtx(self): # 按照一定顺序构建相机内参数 return np.array([[self.f.r[0], 0, self.c.r[0]],[0., self.f.r[1], self.c.r[1]],[0.,0.,1.]], dtype=np.float64) @depends_on('t', 'rt') def view_mtx(self): R = cv2.Rodrigues(self.rt.r)[0] return np.hstack((R,col(self.t.r))) @depends_on('v', 'rt', 't', 'f', 'c', 'k') def r_and_derivatives(self): v = self.v.r.reshape((-1,3)).copy() # 调用opencv的projectpoint函数,将3d点变换到平面,输入参数以此为 3d坐标,旋转矩阵,平移向量,相机内参数,畸形参数,返回值为imagepoint和雅克比矩阵 return cv2.projectPoints(v, self.rt.r, self.t.r, self.camera_mtx, self.k.r) @property def view_matrix(self): R = cv2.Rodrigues(self.rt.r)[0] return np.hstack((R, col(self.t.r))) class ProjectPoints3D(ProjectPoints): dterms = 'v', 'rt', 't', 'f', 'c', 'k' def compute_r(self): result = ProjectPoints.compute_r(self) return np.hstack((result, col(self.z_coords.r))) @property def z_coords(self): assert(self.v.r.shape[1]==3) return RigidTransform(v=self.v, rt=self.rt, t=self.t)[:,2] def compute_dr_wrt(self, wrt): result = ProjectPoints.compute_dr_wrt(self, wrt) if result is None: return None if sp.issparse(result): drz = self.z_coords.dr_wrt(wrt).tocoo() result = result.tocoo() result.row = result.row*3/2 IS = np.concatenate((result.row, drz.row*3+2)) JS = np.concatenate((result.col, drz.col)) data = np.concatenate((result.data, drz.data)) result = sp.csc_matrix((data, (IS, JS)), shape=(self.v.r.size, wrt.r.size)) else: bigger = np.zeros((result.shape[0]/2, 3, result.shape[1])) bigger[:, :2, :] = result.reshape((-1, 2, result.shape[-1])) drz = self.z_coords.dr_wrt(wrt) if drz is not None: if sp.issparse(drz): drz = drz.todense() bigger[:,2,:] = drz.reshape(bigger[:,2,:].shape) result = bigger.reshape((-1, bigger.shape[-1])) return result def main(): import unittest from test_camera import TestCamera suite = unittest.TestLoader().loadTestsFromTestCase(TestCamera) unittest.TextTestRunner(verbosity=2).run(suite) if __name__ == '__main__': main()
2.2 render.py
主要是渲染器,这里介绍使用的color渲染器。emmm更难懂,尽力
</pre><pre>pixel_center_offset = 0.5 class BaseRenderer(Ch): # 基础渲染器 </pre><pre> terms = ['f', 'frustum','overdraw'] dterms = ['camera', 'v'] @depends_on('f') # not v: specifically, it depends only on the number of vertices, not on the values in v def primitives_per_edge(self): v=self.v.r.reshape((-1,3)) f=self.f vpe = get_vertices_per_edge(v, f) # 得到相连接的顶点坐标。比如(11,13)意味着11和13个顶点相连接 fpe = get_faces_per_edge(v, f, vpe) return fpe, vpe @depends_on('f', 'frustum', 'camera', 'overdraw') def barycentric_image(self): self._call_on_changed() return draw_barycentric_image(self.glf, self.v.r, self.f, self.boundarybool_image if self.overdraw else None) @depends_on(terms+dterms) def boundaryid_image(self): self._call_on_changed() return draw_boundaryid_image(self.glb, self.v.r, self.f, self.vpe, self.fpe, self.camera) @depends_on('f', 'frustum', 'camera', 'overdraw') def visibility_image(self): # 可视化图像 self._call_on_changed() return draw_visibility_image(self.glb, self.v.r, self.f, self.boundarybool_image if self.overdraw else None) @depends_on(terms+dterms) def boundarybool_image(self): self._call_on_changed() boundaryid_image = self.boundaryid_image return np.asarray(boundaryid_image != 4294967295, np.uint32).reshape(boundaryid_image.shape) @property def shape(self): raise NotImplementedError('Should be implemented in inherited class.') @property def v(self): return self.camera.v @v.setter def v(self, newval): self.camera.v = newval @property def vpe(self): return self.primitives_per_edge[1] @property def fpe(self): return self.primitives_per_edge[0] class ColoredRenderer(BaseRenderer): terms = 'f', 'frustum', 'background_image', 'overdraw', 'num_channels' dterms = 'vc', 'camera', 'bgcolor' @property def shape(self): if not hasattr(self, 'num_channels'): self.num_channels = 3 if self.num_channels > 1: return (self.frustum['height'], self.frustum['width'], self.num_channels) else: return (self.frustum['height'], self.frustum['width']) def compute_r(self): tmp = self.camera.r return self.color_image # .reshape((self.frustum['height'], self.frustum['width'], -1)).squeeze() def compute_dr_wrt(self, wrt): if wrt is not self.camera and wrt is not self.vc and wrt is not self.bgcolor: return None visibility = self.visibility_image shape = visibility.shape color = self.color_image visible = np.nonzero(visibility.ravel() != 4294967295)[0] num_visible = len(visible) barycentric = self.barycentric_image if wrt is self.camera: if self.overdraw: return common.dImage_wrt_2dVerts_bnd(color, visible, visibility, barycentric, self.frustum['width'], self.frustum['height'], self.v.r.size/3, self.f, self.boundaryid_image != 4294967295) else: return common.dImage_wrt_2dVerts(color, visible, visibility, barycentric, self.frustum['width'], self.frustum['height'], self.v.r.size/3, self.f) elif wrt is self.vc: return common.dr_wrt_vc(visible, visibility, self.f, barycentric, self.frustum, self.vc.size, num_channels=self.num_channels) elif wrt is self.bgcolor: return common.dr_wrt_bgcolor(visibility, self.frustum, num_channels=self.num_channels) def on_changed(self, which): if 'frustum' in which: w = self.frustum['width'] h = self.frustum['height'] self.glf = OsContext(w, h, typ=GL_FLOAT) self.glf.Viewport(0, 0, w, h) self.glb = OsContext(w, h, typ=GL_UNSIGNED_BYTE) self.glb.Viewport(0, 0, w, h) if 'frustum' in which or 'camera' in which: setup_camera(self.glb, self.camera, self.frustum) setup_camera(self.glf, self.camera, self.frustum) if not hasattr(self, 'num_channels'): self.num_channels = 3 if not hasattr(self, 'bgcolor'): self.bgcolor = Ch(np.array([.5]*self.num_channels)) which.add('bgcolor') if not hasattr(self, 'overdraw'): self.overdraw = True if 'bgcolor' in which or ('frustum' in which and hasattr(self, 'bgcolor')): self.glf.ClearColor(self.bgcolor.r[0], self.bgcolor.r[1%self.num_channels], self.bgcolor.r[2%self.num_channels], 1.) def flow_to(self, v_next, cam_next=None): return common.flow_to(self, v_next, cam_next) def filter_for_triangles(self, which_triangles): cim = self.color_image vim = self.visibility_image+1 arr = np.zeros(len(self.f)+1) arr[which_triangles+1] = 1 relevant_pixels = arr[vim.ravel()] cim2 = cim.copy() * np.atleast_3d(relevant_pixels.reshape(vim.shape)) relevant_pixels = np.nonzero(arr[vim.ravel()])[0] xs = relevant_pixels % vim.shape[1] ys = relevant_pixels / vim.shape[1] return cim2[np.min(ys):np.max(ys), np.min(xs):np.max(xs), :] @depends_on('f', 'camera', 'vc') def boundarycolor_image(self): return self.draw_boundarycolor_image(with_vertex_colors=True) def draw_color_image(self, gl): self._call_on_changed() gl.Clear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); # use face colors if given # FIXME: this won't work for 2 channels draw_colored_verts(gl, self.v.r, self.f, self.vc.r) result = np.asarray(deepcopy(gl.getImage()[:,:,:self.num_channels].squeeze()), np.float64) if hasattr(self, 'background_image'): bg_px = np.tile(np.atleast_3d(self.visibility_image) == 4294967295, (1,1,self.num_channels)).squeeze() fg_px = 1 - bg_px result = bg_px * self.background_image + fg_px * result return result @depends_on(dterms+terms) def color_image(self): gl = self.glf gl.PolygonMode(GL_FRONT_AND_BACK, GL_FILL) no_overdraw = self.draw_color_image(gl) if not self.overdraw: return no_overdraw gl.PolygonMode(GL_FRONT_AND_BACK, GL_LINE) overdraw = self.draw_color_image(gl) gl.PolygonMode(GL_FRONT_AND_BACK, GL_FILL) boundarybool_image = self.boundarybool_image if self.num_channels > 1: boundarybool_image = np.atleast_3d(boundarybool_image) return np.asarray((overdraw*boundarybool_image + no_overdraw*(1-boundarybool_image)), order='C') @depends_on('f', 'frustum', 'camera') def boundary_images(self): self._call_on_changed() return draw_boundary_images(self.glb, self.v.r, self.f, self.vpe, self.fpe, self.camera) @depends_on(terms+dterms) def boundarycolor_image(self): self._call_on_changed() gl = self.glf colors = self.vc.r.reshape((-1,3))[self.vpe.ravel()] gl.Clear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); draw_colored_primitives(gl, self.v.r.reshape((-1,3)), self.vpe, colors) return np.asarray(deepcopy(gl.getImage()), np.float64)
2.3 lighting.py
感觉没啥可注释的。。。唯一需要注释的就是参数的含义,后续再补坑
</pre> <pre>class LambertianPointLight(Ch) : # 郎博提点光 terms = 'f', 'num_verts', 'light_color' # dterms = 'light_pos', 'v', 'vc', 'vn' # def on_changed(self, which): # which就是上述7个参数和下面的_lpl构成的元组 if not hasattr(self, '_lpl'): self.add_dterm('_lpl', maximum(multiply(a=multiply()), 0.0)) if not hasattr(self, 'ldn'): self.ldn = LightDotNormal(self.v.r.size/3) if not hasattr(self, 'vn'): logger.info('LambertianPointLight using auto-normals. This will be slow for derivative-free computations.') self.vn = VertNormals(f=self.f, v=self.v) self.vn.needs_autoupdate = True if 'v' in which and hasattr(self.vn, 'needs_autoupdate') and self.vn.needs_autoupdate: self.vn.v = self.v ldn_args = {k: getattr(self, k) for k in which if k in ('light_pos', 'v', 'vn')} if len(ldn_args) > 0: self.ldn.set(**ldn_args) self._lpl.a.a.a = self.ldn.reshape((-1,1)) if 'num_verts' in which or 'light_color' in which: # nc = self.num_channels # IS = np.arange(self.num_verts*nc) # JS = np.repeat(np.arange(self.num_verts), 3) # data = (row(self.light_color)*np.ones((self.num_verts, 3))).ravel() # mtx = sp.csc_matrix((data, (IS,JS)), shape=(self.num_verts*3, self.num_verts)) self._lpl.a.a.b = self.light_color.reshape((1,self.num_channels)) if 'vc' in which: self._lpl.a.b = self.vc.reshape((-1,self.num_channels)) @property def num_channels(self): return self.light_color.size def compute_r(self): return self._lpl.r def compute_dr_wrt(self, wrt): if wrt is self._lpl: return 1 # def compute_light_repeat(num_verts): # IS = np.arange(num_verts*3) # JS = IS % 3 # data = np.ones_like(IS, dtype=np.float64) # ij = np.vstack((row(IS), row(JS))) # return sp.csc_matrix((data, ij), shape=(num_verts*3, 3))</pre> <pre>
23
没有评论
你先离开吧:)