本文主要是对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