(61d00a474) v0.9.7.1

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/MotionEST/Anandan.py

@@ -0,0 +1,185 @@
+#!/usr/bin/env python
+# coding: utf-8
+import numpy as np
+import numpy.linalg as LA
+from scipy.ndimage.filters import gaussian_filter
+from scipy.sparse import csc_matrix
+from scipy.sparse.linalg import inv
+from MotionEST import MotionEST
+"""Anandan Model"""
+
+
+class Anandan(MotionEST):
+  """
+    constructor:
+        cur_f: current frame
+        ref_f: reference frame
+        blk_sz: block size
+        beta: smooth constrain weight
+        k1,k2,k3: confidence coefficients
+        max_iter: maximum number of iterations
+    """
+
+  def __init__(self, cur_f, ref_f, blk_sz, beta, k1, k2, k3, max_iter=100):
+    super(Anandan, self).__init__(cur_f, ref_f, blk_sz)
+    self.levels = int(np.log2(blk_sz))
+    self.intensity_hierarchy()
+    self.c_maxs = []
+    self.c_mins = []
+    self.e_maxs = []
+    self.e_mins = []
+    for l in xrange(self.levels + 1):
+      c_max, c_min, e_max, e_min = self.get_curvature(self.cur_Is[l])
+      self.c_maxs.append(c_max)
+      self.c_mins.append(c_min)
+      self.e_maxs.append(e_max)
+      self.e_mins.append(e_min)
+    self.beta = beta
+    self.k1, self.k2, self.k3 = k1, k2, k3
+    self.max_iter = max_iter
+
+  """
+    build intensity hierarchy
+    """
+
+  def intensity_hierarchy(self):
+    level = 0
+    self.cur_Is = []
+    self.ref_Is = []
+    #build each level itensity by using gaussian filters
+    while level <= self.levels:
+      cur_I = gaussian_filter(self.cur_yuv[:, :, 0], sigma=(2**level) * 0.56)
+      ref_I = gaussian_filter(self.ref_yuv[:, :, 0], sigma=(2**level) * 0.56)
+      self.ref_Is.append(ref_I)
+      self.cur_Is.append(cur_I)
+      level += 1
+
+  """
+    get curvature of each block
+    """
+
+  def get_curvature(self, I):
+    c_max = np.zeros((self.num_row, self.num_col))
+    c_min = np.zeros((self.num_row, self.num_col))
+    e_max = np.zeros((self.num_row, self.num_col, 2))
+    e_min = np.zeros((self.num_row, self.num_col, 2))
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        h11, h12, h21, h22 = 0, 0, 0, 0
+        for i in xrange(r * self.blk_sz, r * self.blk_sz + self.blk_sz):
+          for j in xrange(c * self.blk_sz, c * self.blk_sz + self.blk_sz):
+            if 0 <= i < self.height - 1 and 0 <= j < self.width - 1:
+              Ix = I[i][j + 1] - I[i][j]
+              Iy = I[i + 1][j] - I[i][j]
+              h11 += Iy * Iy
+              h12 += Ix * Iy
+              h21 += Ix * Iy
+              h22 += Ix * Ix
+        U, S, _ = LA.svd(np.array([[h11, h12], [h21, h22]]))
+        c_max[r, c], c_min[r, c] = S[0], S[1]
+        e_max[r, c] = U[:, 0]
+        e_min[r, c] = U[:, 1]
+    return c_max, c_min, e_max, e_min
+
+  """
+    get ssd of motion vector:
+      cur_I: current intensity
+      ref_I: reference intensity
+      center: current position
+      mv: motion vector
+    """
+
+  def get_ssd(self, cur_I, ref_I, center, mv):
+    ssd = 0
+    for r in xrange(int(center[0]), int(center[0]) + self.blk_sz):
+      for c in xrange(int(center[1]), int(center[1]) + self.blk_sz):
+        if 0 <= r < self.height and 0 <= c < self.width:
+          tr, tc = r + int(mv[0]), c + int(mv[1])
+          if 0 <= tr < self.height and 0 <= tc < self.width:
+            ssd += (ref_I[tr, tc] - cur_I[r, c])**2
+          else:
+            ssd += cur_I[r, c]**2
+    return ssd
+
+  """
+    get region match of level l
+      l: current level
+      last_mvs: matchine results of last level
+      radius: movenment radius
+    """
+
+  def region_match(self, l, last_mvs, radius):
+    mvs = np.zeros((self.num_row, self.num_col, 2))
+    min_ssds = np.zeros((self.num_row, self.num_col))
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        center = np.array([r * self.blk_sz, c * self.blk_sz])
+        #use overlap hierarchy policy
+        init_mvs = []
+        if last_mvs is None:
+          init_mvs = [np.array([0, 0])]
+        else:
+          for i, j in {(r, c), (r, c + 1), (r + 1, c), (r + 1, c + 1)}:
+            if 0 <= i < last_mvs.shape[0] and 0 <= j < last_mvs.shape[1]:
+              init_mvs.append(last_mvs[i, j])
+        #use last matching results as the start postion as current level
+        min_ssd = None
+        min_mv = None
+        for init_mv in init_mvs:
+          for i in xrange(-2, 3):
+            for j in xrange(-2, 3):
+              mv = init_mv + np.array([i, j]) * radius
+              ssd = self.get_ssd(self.cur_Is[l], self.ref_Is[l], center, mv)
+              if min_ssd is None or ssd < min_ssd:
+                min_ssd = ssd
+                min_mv = mv
+        min_ssds[r, c] = min_ssd
+        mvs[r, c] = min_mv
+    return mvs, min_ssds
+
+  """
+    smooth motion field based on neighbor constraint
+      uvs: current estimation
+      mvs: matching results
+      min_ssds: minimum ssd of matching results
+      l: current level
+    """
+
+  def smooth(self, uvs, mvs, min_ssds, l):
+    sm_uvs = np.zeros((self.num_row, self.num_col, 2))
+    c_max = self.c_maxs[l]
+    c_min = self.c_mins[l]
+    e_max = self.e_maxs[l]
+    e_min = self.e_mins[l]
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        w_max = c_max[r, c] / (
+            self.k1 + self.k2 * min_ssds[r, c] + self.k3 * c_max[r, c])
+        w_min = c_min[r, c] / (
+            self.k1 + self.k2 * min_ssds[r, c] + self.k3 * c_min[r, c])
+        w = w_max * w_min / (w_max + w_min + 1e-6)
+        if w < 0:
+          w = 0
+        avg_uv = np.array([0.0, 0.0])
+        for i, j in {(r - 1, c), (r + 1, c), (r, c - 1), (r, c + 1)}:
+          if 0 <= i < self.num_row and 0 <= j < self.num_col:
+            avg_uv += 0.25 * uvs[i, j]
+        sm_uvs[r, c] = (w * w * mvs[r, c] + self.beta * avg_uv) / (
+            self.beta + w * w)
+    return sm_uvs
+
+  """
+    motion field estimation
+    """
+
+  def motion_field_estimation(self):
+    last_mvs = None
+    for l in xrange(self.levels, -1, -1):
+      mvs, min_ssds = self.region_match(l, last_mvs, 2**l)
+      uvs = np.zeros(mvs.shape)
+      for _ in xrange(self.max_iter):
+        uvs = self.smooth(uvs, mvs, min_ssds, l)
+      last_mvs = uvs
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        self.mf[r, c] = uvs[r, c]

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/MotionEST/Exhaust.py

@@ -0,0 +1,251 @@
+#!/usr/bin/env python
+# coding: utf-8
+import numpy as np
+import numpy.linalg as LA
+from Util import MSE
+from MotionEST import MotionEST
+"""Exhaust Search:"""
+
+
+class Exhaust(MotionEST):
+  """
+    Constructor:
+        cur_f: current frame
+        ref_f: reference frame
+        blk_sz: block size
+        wnd_size: search window size
+        metric: metric to compare the blocks distrotion
+    """
+
+  def __init__(self, cur_f, ref_f, blk_size, wnd_size, metric=MSE):
+    self.name = 'exhaust'
+    self.wnd_sz = wnd_size
+    self.metric = metric
+    super(Exhaust, self).__init__(cur_f, ref_f, blk_size)
+
+  """
+    search method:
+        cur_r: start row
+        cur_c: start column
+    """
+
+  def search(self, cur_r, cur_c):
+    min_loss = self.block_dist(cur_r, cur_c, [0, 0], self.metric)
+    cur_x = cur_c * self.blk_sz
+    cur_y = cur_r * self.blk_sz
+    ref_x = cur_x
+    ref_y = cur_y
+    #search all validate positions and select the one with minimum distortion
+    for y in xrange(cur_y - self.wnd_sz, cur_y + self.wnd_sz):
+      for x in xrange(cur_x - self.wnd_sz, cur_x + self.wnd_sz):
+        if 0 <= x < self.width - self.blk_sz and 0 <= y < self.height - self.blk_sz:
+          loss = self.block_dist(cur_r, cur_c, [y - cur_y, x - cur_x],
+                                 self.metric)
+          if loss < min_loss:
+            min_loss = loss
+            ref_x = x
+            ref_y = y
+    return ref_x, ref_y
+
+  def motion_field_estimation(self):
+    for i in xrange(self.num_row):
+      for j in xrange(self.num_col):
+        ref_x, ref_y = self.search(i, j)
+        self.mf[i, j] = np.array(
+            [ref_y - i * self.blk_sz, ref_x - j * self.blk_sz])
+
+
+"""Exhaust with Neighbor Constraint"""
+
+
+class ExhaustNeighbor(MotionEST):
+  """
+    Constructor:
+        cur_f: current frame
+        ref_f: reference frame
+        blk_sz: block size
+        wnd_size: search window size
+        beta: neigbor loss weight
+        metric: metric to compare the blocks distrotion
+    """
+
+  def __init__(self, cur_f, ref_f, blk_size, wnd_size, beta, metric=MSE):
+    self.name = 'exhaust + neighbor'
+    self.wnd_sz = wnd_size
+    self.beta = beta
+    self.metric = metric
+    super(ExhaustNeighbor, self).__init__(cur_f, ref_f, blk_size)
+    self.assign = np.zeros((self.num_row, self.num_col), dtype=np.bool)
+
+  """
+    estimate neighbor loss:
+        cur_r: current row
+        cur_c: current column
+        mv: current motion vector
+    """
+
+  def neighborLoss(self, cur_r, cur_c, mv):
+    loss = 0
+    #accumulate difference between current block's motion vector with neighbors'
+    for i, j in {(-1, 0), (1, 0), (0, 1), (0, -1)}:
+      nb_r = cur_r + i
+      nb_c = cur_c + j
+      if 0 <= nb_r < self.num_row and 0 <= nb_c < self.num_col and self.assign[
+          nb_r, nb_c]:
+        loss += LA.norm(mv - self.mf[nb_r, nb_c])
+    return loss
+
+  """
+    search method:
+        cur_r: start row
+        cur_c: start column
+    """
+
+  def search(self, cur_r, cur_c):
+    dist_loss = self.block_dist(cur_r, cur_c, [0, 0], self.metric)
+    nb_loss = self.neighborLoss(cur_r, cur_c, np.array([0, 0]))
+    min_loss = dist_loss + self.beta * nb_loss
+    cur_x = cur_c * self.blk_sz
+    cur_y = cur_r * self.blk_sz
+    ref_x = cur_x
+    ref_y = cur_y
+    #search all validate positions and select the one with minimum distortion
+    # as well as weighted neighbor loss
+    for y in xrange(cur_y - self.wnd_sz, cur_y + self.wnd_sz):
+      for x in xrange(cur_x - self.wnd_sz, cur_x + self.wnd_sz):
+        if 0 <= x < self.width - self.blk_sz and 0 <= y < self.height - self.blk_sz:
+          dist_loss = self.block_dist(cur_r, cur_c, [y - cur_y, x - cur_x],
+                                      self.metric)
+          nb_loss = self.neighborLoss(cur_r, cur_c, [y - cur_y, x - cur_x])
+          loss = dist_loss + self.beta * nb_loss
+          if loss < min_loss:
+            min_loss = loss
+            ref_x = x
+            ref_y = y
+    return ref_x, ref_y
+
+  def motion_field_estimation(self):
+    for i in xrange(self.num_row):
+      for j in xrange(self.num_col):
+        ref_x, ref_y = self.search(i, j)
+        self.mf[i, j] = np.array(
+            [ref_y - i * self.blk_sz, ref_x - j * self.blk_sz])
+        self.assign[i, j] = True
+
+
+"""Exhaust with Neighbor Constraint and Feature Score"""
+
+
+class ExhaustNeighborFeatureScore(MotionEST):
+  """
+    Constructor:
+        cur_f: current frame
+        ref_f: reference frame
+        blk_sz: block size
+        wnd_size: search window size
+        beta: neigbor loss weight
+        max_iter: maximum number of iterations
+        metric: metric to compare the blocks distrotion
+    """
+
+  def __init__(self,
+               cur_f,
+               ref_f,
+               blk_size,
+               wnd_size,
+               beta=1,
+               max_iter=100,
+               metric=MSE):
+    self.name = 'exhaust + neighbor+feature score'
+    self.wnd_sz = wnd_size
+    self.beta = beta
+    self.metric = metric
+    self.max_iter = max_iter
+    super(ExhaustNeighborFeatureScore, self).__init__(cur_f, ref_f, blk_size)
+    self.fs = self.getFeatureScore()
+
+  """
+    get feature score of each block
+    """
+
+  def getFeatureScore(self):
+    fs = np.zeros((self.num_row, self.num_col))
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        IxIx = 0
+        IyIy = 0
+        IxIy = 0
+        #get ssd surface
+        for x in xrange(self.blk_sz - 1):
+          for y in xrange(self.blk_sz - 1):
+            ox = c * self.blk_sz + x
+            oy = r * self.blk_sz + y
+            Ix = self.cur_yuv[oy, ox + 1, 0] - self.cur_yuv[oy, ox, 0]
+            Iy = self.cur_yuv[oy + 1, ox, 0] - self.cur_yuv[oy, ox, 0]
+            IxIx += Ix * Ix
+            IyIy += Iy * Iy
+            IxIy += Ix * Iy
+        #get maximum and minimum eigenvalues
+        lambda_max = 0.5 * ((IxIx + IyIy) + np.sqrt(4 * IxIy * IxIy +
+                                                    (IxIx - IyIy)**2))
+        lambda_min = 0.5 * ((IxIx + IyIy) - np.sqrt(4 * IxIy * IxIy +
+                                                    (IxIx - IyIy)**2))
+        fs[r, c] = lambda_max * lambda_min / (1e-6 + lambda_max + lambda_min)
+        if fs[r, c] < 0:
+          fs[r, c] = 0
+    return fs
+
+  """
+    do exhaust search
+    """
+
+  def search(self, cur_r, cur_c):
+    min_loss = self.block_dist(cur_r, cur_c, [0, 0], self.metric)
+    cur_x = cur_c * self.blk_sz
+    cur_y = cur_r * self.blk_sz
+    ref_x = cur_x
+    ref_y = cur_y
+    #search all validate positions and select the one with minimum distortion
+    for y in xrange(cur_y - self.wnd_sz, cur_y + self.wnd_sz):
+      for x in xrange(cur_x - self.wnd_sz, cur_x + self.wnd_sz):
+        if 0 <= x < self.width - self.blk_sz and 0 <= y < self.height - self.blk_sz:
+          loss = self.block_dist(cur_r, cur_c, [y - cur_y, x - cur_x],
+                                 self.metric)
+          if loss < min_loss:
+            min_loss = loss
+            ref_x = x
+            ref_y = y
+    return ref_x, ref_y
+
+  """
+    add smooth constraint
+    """
+
+  def smooth(self, uvs, mvs):
+    sm_uvs = np.zeros(uvs.shape)
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        avg_uv = np.array([0.0, 0.0])
+        for i, j in {(r - 1, c), (r + 1, c), (r, c - 1), (r, c + 1)}:
+          if 0 <= i < self.num_row and 0 <= j < self.num_col:
+            avg_uv += uvs[i, j] / 6.0
+        for i, j in {(r - 1, c - 1), (r - 1, c + 1), (r + 1, c - 1),
+                     (r + 1, c + 1)}:
+          if 0 <= i < self.num_row and 0 <= j < self.num_col:
+            avg_uv += uvs[i, j] / 12.0
+        sm_uvs[r, c] = (self.fs[r, c] * mvs[r, c] + self.beta * avg_uv) / (
+            self.beta + self.fs[r, c])
+    return sm_uvs
+
+  def motion_field_estimation(self):
+    #get matching results
+    mvs = np.zeros(self.mf.shape)
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        ref_x, ref_y = self.search(r, c)
+        mvs[r, c] = np.array([ref_y - r * self.blk_sz, ref_x - c * self.blk_sz])
+    #add smoothness constraint
+    uvs = np.zeros(self.mf.shape)
+    for _ in xrange(self.max_iter):
+      uvs = self.smooth(uvs, mvs)
+    self.mf = uvs

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/MotionEST/GroundTruth.py

@@ -0,0 +1,40 @@
+#!/ usr / bin / env python
+#coding : utf - 8
+import numpy as np
+import numpy.linalg as LA
+from MotionEST import MotionEST
+"""Ground Truth:
+
+  Load in ground truth motion field and mask
+"""
+
+
+class GroundTruth(MotionEST):
+  """constructor:
+
+    cur_f:current
+    frame ref_f:reference
+    frame blk_sz:block size
+    gt_path:ground truth motion field file path
+    """
+
+  def __init__(self, cur_f, ref_f, blk_sz, gt_path, mf=None, mask=None):
+    self.name = 'ground truth'
+    super(GroundTruth, self).__init__(cur_f, ref_f, blk_sz)
+    self.mask = np.zeros((self.num_row, self.num_col), dtype=np.bool)
+    if gt_path:
+      with open(gt_path) as gt_file:
+        lines = gt_file.readlines()
+        for i in xrange(len(lines)):
+          info = lines[i].split(';')
+          for j in xrange(len(info)):
+            x, y = info[j].split(',')
+            #-, - stands for nothing
+            if x == '-' or y == '-':
+              self.mask[i, -j - 1] = True
+              continue
+            #the order of original file is flipped on the x axis
+            self.mf[i, -j - 1] = np.array([float(y), -float(x)], dtype=np.int)
+    else:
+      self.mf = mf
+      self.mask = mask

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/MotionEST/HornSchunck.py

@@ -0,0 +1,204 @@
+#!/usr/bin/env python
+# coding: utf-8
+import numpy as np
+import numpy.linalg as LA
+from scipy.ndimage.filters import gaussian_filter
+from scipy.sparse import csc_matrix
+from scipy.sparse.linalg import inv
+from MotionEST import MotionEST
+"""Horn & Schunck Model"""
+
+
+class HornSchunck(MotionEST):
+  """
+    constructor:
+        cur_f: current frame
+        ref_f: reference frame
+        blk_sz: block size
+        alpha: smooth constrain weight
+        sigma: gaussian blur parameter
+    """
+
+  def __init__(self, cur_f, ref_f, blk_sz, alpha, sigma, max_iter=100):
+    super(HornSchunck, self).__init__(cur_f, ref_f, blk_sz)
+    self.cur_I, self.ref_I = self.getIntensity()
+    #perform gaussian blur to smooth the intensity
+    self.cur_I = gaussian_filter(self.cur_I, sigma=sigma)
+    self.ref_I = gaussian_filter(self.ref_I, sigma=sigma)
+    self.alpha = alpha
+    self.max_iter = max_iter
+    self.Ix, self.Iy, self.It = self.intensityDiff()
+
+  """
+    Build Frame Intensity
+    """
+
+  def getIntensity(self):
+    cur_I = np.zeros((self.num_row, self.num_col))
+    ref_I = np.zeros((self.num_row, self.num_col))
+    #use average intensity as block's intensity
+    for i in xrange(self.num_row):
+      for j in xrange(self.num_col):
+        r = i * self.blk_sz
+        c = j * self.blk_sz
+        cur_I[i, j] = np.mean(self.cur_yuv[r:r + self.blk_sz, c:c + self.blk_sz,
+                                           0])
+        ref_I[i, j] = np.mean(self.ref_yuv[r:r + self.blk_sz, c:c + self.blk_sz,
+                                           0])
+    return cur_I, ref_I
+
+  """
+    Get First Order Derivative
+    """
+
+  def intensityDiff(self):
+    Ix = np.zeros((self.num_row, self.num_col))
+    Iy = np.zeros((self.num_row, self.num_col))
+    It = np.zeros((self.num_row, self.num_col))
+    sz = self.blk_sz
+    for i in xrange(self.num_row - 1):
+      for j in xrange(self.num_col - 1):
+        """
+                Ix:
+                (i  ,j) <--- (i  ,j+1)
+                (i+1,j) <--- (i+1,j+1)
+                """
+        count = 0
+        for r, c in {(i, j + 1), (i + 1, j + 1)}:
+          if 0 <= r < self.num_row and 0 < c < self.num_col:
+            Ix[i, j] += (
+                self.cur_I[r, c] - self.cur_I[r, c - 1] + self.ref_I[r, c] -
+                self.ref_I[r, c - 1])
+            count += 2
+        Ix[i, j] /= count
+        """
+                Iy:
+                (i  ,j)      (i  ,j+1)
+                   ^             ^
+                   |             |
+                (i+1,j)      (i+1,j+1)
+                """
+        count = 0
+        for r, c in {(i + 1, j), (i + 1, j + 1)}:
+          if 0 < r < self.num_row and 0 <= c < self.num_col:
+            Iy[i, j] += (
+                self.cur_I[r, c] - self.cur_I[r - 1, c] + self.ref_I[r, c] -
+                self.ref_I[r - 1, c])
+            count += 2
+        Iy[i, j] /= count
+        count = 0
+        #It:
+        for r in xrange(i, i + 2):
+          for c in xrange(j, j + 2):
+            if 0 <= r < self.num_row and 0 <= c < self.num_col:
+              It[i, j] += (self.ref_I[r, c] - self.cur_I[r, c])
+              count += 1
+        It[i, j] /= count
+    return Ix, Iy, It
+
+  """
+    Get weighted average of neighbor motion vectors
+    for evaluation of laplacian
+    """
+
+  def averageMV(self):
+    avg = np.zeros((self.num_row, self.num_col, 2))
+    """
+        1/12 ---  1/6 --- 1/12
+         |         |       |
+        1/6  --- -1/8 --- 1/6
+         |         |       |
+        1/12 ---  1/6 --- 1/12
+        """
+    for i in xrange(self.num_row):
+      for j in xrange(self.num_col):
+        for r, c in {(-1, 0), (1, 0), (0, -1), (0, 1)}:
+          if 0 <= i + r < self.num_row and 0 <= j + c < self.num_col:
+            avg[i, j] += self.mf[i + r, j + c] / 6.0
+        for r, c in {(-1, -1), (-1, 1), (1, -1), (1, 1)}:
+          if 0 <= i + r < self.num_row and 0 <= j + c < self.num_col:
+            avg[i, j] += self.mf[i + r, j + c] / 12.0
+    return avg
+
+  def motion_field_estimation(self):
+    count = 0
+    """
+        u_{n+1} = ~u_n - Ix(Ix.~u_n+Iy.~v+It)/(IxIx+IyIy+alpha^2)
+        v_{n+1} = ~v_n - Iy(Ix.~u_n+Iy.~v+It)/(IxIx+IyIy+alpha^2)
+        """
+    denom = self.alpha**2 + np.power(self.Ix, 2) + np.power(self.Iy, 2)
+    while count < self.max_iter:
+      avg = self.averageMV()
+      self.mf[:, :, 1] = avg[:, :, 1] - self.Ix * (
+          self.Ix * avg[:, :, 1] + self.Iy * avg[:, :, 0] + self.It) / denom
+      self.mf[:, :, 0] = avg[:, :, 0] - self.Iy * (
+          self.Ix * avg[:, :, 1] + self.Iy * avg[:, :, 0] + self.It) / denom
+      count += 1
+    self.mf *= self.blk_sz
+
+  def motion_field_estimation_mat(self):
+    row_idx = []
+    col_idx = []
+    data = []
+
+    N = 2 * self.num_row * self.num_col
+    b = np.zeros((N, 1))
+    for i in xrange(self.num_row):
+      for j in xrange(self.num_col):
+        """(IxIx+alpha^2)u+IxIy.v-alpha^2~u IxIy.u+(IyIy+alpha^2)v-alpha^2~v"""
+        u_idx = i * 2 * self.num_col + 2 * j
+        v_idx = u_idx + 1
+        b[u_idx, 0] = -self.Ix[i, j] * self.It[i, j]
+        b[v_idx, 0] = -self.Iy[i, j] * self.It[i, j]
+        #u: (IxIx+alpha^2)u
+        row_idx.append(u_idx)
+        col_idx.append(u_idx)
+        data.append(self.Ix[i, j] * self.Ix[i, j] + self.alpha**2)
+        #IxIy.v
+        row_idx.append(u_idx)
+        col_idx.append(v_idx)
+        data.append(self.Ix[i, j] * self.Iy[i, j])
+
+        #v: IxIy.u
+        row_idx.append(v_idx)
+        col_idx.append(u_idx)
+        data.append(self.Ix[i, j] * self.Iy[i, j])
+        #(IyIy+alpha^2)v
+        row_idx.append(v_idx)
+        col_idx.append(v_idx)
+        data.append(self.Iy[i, j] * self.Iy[i, j] + self.alpha**2)
+
+        #-alpha^2~u
+        #-alpha^2~v
+        for r, c in {(-1, 0), (1, 0), (0, -1), (0, 1)}:
+          if 0 <= i + r < self.num_row and 0 <= j + c < self.num_col:
+            u_nb = (i + r) * 2 * self.num_col + 2 * (j + c)
+            v_nb = u_nb + 1
+
+            row_idx.append(u_idx)
+            col_idx.append(u_nb)
+            data.append(-1 * self.alpha**2 / 6.0)
+
+            row_idx.append(v_idx)
+            col_idx.append(v_nb)
+            data.append(-1 * self.alpha**2 / 6.0)
+        for r, c in {(-1, -1), (-1, 1), (1, -1), (1, 1)}:
+          if 0 <= i + r < self.num_row and 0 <= j + c < self.num_col:
+            u_nb = (i + r) * 2 * self.num_col + 2 * (j + c)
+            v_nb = u_nb + 1
+
+            row_idx.append(u_idx)
+            col_idx.append(u_nb)
+            data.append(-1 * self.alpha**2 / 12.0)
+
+            row_idx.append(v_idx)
+            col_idx.append(v_nb)
+            data.append(-1 * self.alpha**2 / 12.0)
+    M = csc_matrix((data, (row_idx, col_idx)), shape=(N, N))
+    M_inv = inv(M)
+    uv = M_inv.dot(b)
+
+    for i in xrange(self.num_row):
+      for j in xrange(self.num_col):
+        self.mf[i, j, 0] = uv[i * 2 * self.num_col + 2 * j + 1, 0] * self.blk_sz
+        self.mf[i, j, 1] = uv[i * 2 * self.num_col + 2 * j, 0] * self.blk_sz

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/MotionEST/MotionEST.py

@@ -0,0 +1,109 @@
+#!/ usr / bin / env python
+#coding : utf - 8
+import numpy as np
+import numpy.linalg as LA
+import matplotlib.pyplot as plt
+from Util import drawMF, MSE
+"""The Base Class of Estimators"""
+
+
+class MotionEST(object):
+  """
+    constructor:
+        cur_f: current frame
+        ref_f: reference frame
+        blk_sz: block size
+    """
+
+  def __init__(self, cur_f, ref_f, blk_sz):
+    self.cur_f = cur_f
+    self.ref_f = ref_f
+    self.blk_sz = blk_sz
+    #convert RGB to YUV
+    self.cur_yuv = np.array(self.cur_f.convert('YCbCr'), dtype=np.int)
+    self.ref_yuv = np.array(self.ref_f.convert('YCbCr'), dtype=np.int)
+    #frame size
+    self.width = self.cur_f.size[0]
+    self.height = self.cur_f.size[1]
+    #motion field size
+    self.num_row = self.height // self.blk_sz
+    self.num_col = self.width // self.blk_sz
+    #initialize motion field
+    self.mf = np.zeros((self.num_row, self.num_col, 2))
+
+  """estimation function Override by child classes"""
+
+  def motion_field_estimation(self):
+    pass
+
+  """
+    distortion of a block:
+      cur_r: current row
+      cur_c: current column
+      mv: motion vector
+      metric: distortion metric
+  """
+
+  def block_dist(self, cur_r, cur_c, mv, metric=MSE):
+    cur_x = cur_c * self.blk_sz
+    cur_y = cur_r * self.blk_sz
+    h = min(self.blk_sz, self.height - cur_y)
+    w = min(self.blk_sz, self.width - cur_x)
+    cur_blk = self.cur_yuv[cur_y:cur_y + h, cur_x:cur_x + w, :]
+    ref_x = int(cur_x + mv[1])
+    ref_y = int(cur_y + mv[0])
+    if 0 <= ref_x < self.width - w and 0 <= ref_y < self.height - h:
+      ref_blk = self.ref_yuv[ref_y:ref_y + h, ref_x:ref_x + w, :]
+    else:
+      ref_blk = np.zeros((h, w, 3))
+    return metric(cur_blk, ref_blk)
+
+  """
+    distortion of motion field
+  """
+
+  def distortion(self, mask=None, metric=MSE):
+    loss = 0
+    count = 0
+    for i in xrange(self.num_row):
+      for j in xrange(self.num_col):
+        if mask is not None and mask[i, j]:
+          continue
+        loss += self.block_dist(i, j, self.mf[i, j], metric)
+        count += 1
+    return loss / count
+
+  """evaluation compare the difference with ground truth"""
+
+  def motion_field_evaluation(self, ground_truth):
+    loss = 0
+    count = 0
+    gt = ground_truth.mf
+    mask = ground_truth.mask
+    for i in xrange(self.num_row):
+      for j in xrange(self.num_col):
+        if mask is not None and mask[i][j]:
+          continue
+        loss += LA.norm(gt[i, j] - self.mf[i, j])
+        count += 1
+    return loss / count
+
+  """render the motion field"""
+
+  def show(self, ground_truth=None, size=10):
+    cur_mf = drawMF(self.cur_f, self.blk_sz, self.mf)
+    if ground_truth is None:
+      n_row = 1
+    else:
+      gt_mf = drawMF(self.cur_f, self.blk_sz, ground_truth)
+      n_row = 2
+    plt.figure(figsize=(n_row * size, size * self.height / self.width))
+    plt.subplot(1, n_row, 1)
+    plt.imshow(cur_mf)
+    plt.title('Estimated Motion Field')
+    if ground_truth is not None:
+      plt.subplot(1, n_row, 2)
+      plt.imshow(gt_mf)
+      plt.title('Ground Truth')
+    plt.tight_layout()
+    plt.show()

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/MotionEST/SearchSmooth.py

@@ -0,0 +1,213 @@
+#!/usr/bin/env python
+# coding: utf-8
+import numpy as np
+import numpy.linalg as LA
+from Util import MSE
+from MotionEST import MotionEST
+"""Search & Smooth Model with Adapt Weights"""
+
+
+class SearchSmoothAdapt(MotionEST):
+  """
+    Constructor:
+        cur_f: current frame
+        ref_f: reference frame
+        blk_sz: block size
+        wnd_size: search window size
+        beta: neigbor loss weight
+        max_iter: maximum number of iterations
+        metric: metric to compare the blocks distrotion
+    """
+
+  def __init__(self, cur_f, ref_f, blk_size, search, max_iter=100):
+    self.search = search
+    self.max_iter = max_iter
+    super(SearchSmoothAdapt, self).__init__(cur_f, ref_f, blk_size)
+
+  """
+    get local diffiencial of refernce
+    """
+
+  def getRefLocalDiff(self, mvs):
+    m, n = self.num_row, self.num_col
+    localDiff = [[] for _ in xrange(m)]
+    blk_sz = self.blk_sz
+    for r in xrange(m):
+      for c in xrange(n):
+        I_row = 0
+        I_col = 0
+        #get ssd surface
+        count = 0
+        center = self.cur_yuv[r * blk_sz:(r + 1) * blk_sz,
+                              c * blk_sz:(c + 1) * blk_sz, 0]
+        ty = np.clip(r * blk_sz + int(mvs[r, c, 0]), 0, self.height - blk_sz)
+        tx = np.clip(c * blk_sz + int(mvs[r, c, 1]), 0, self.width - blk_sz)
+        target = self.ref_yuv[ty:ty + blk_sz, tx:tx + blk_sz, 0]
+        for y, x in {(ty - blk_sz, tx), (ty + blk_sz, tx)}:
+          if 0 <= y < self.height - blk_sz and 0 <= x < self.width - blk_sz:
+            nb = self.ref_yuv[y:y + blk_sz, x:x + blk_sz, 0]
+            I_row += np.sum(np.abs(nb - center)) - np.sum(
+                np.abs(target - center))
+            count += 1
+        I_row //= (count * blk_sz * blk_sz)
+        count = 0
+        for y, x in {(ty, tx - blk_sz), (ty, tx + blk_sz)}:
+          if 0 <= y < self.height - blk_sz and 0 <= x < self.width - blk_sz:
+            nb = self.ref_yuv[y:y + blk_sz, x:x + blk_sz, 0]
+            I_col += np.sum(np.abs(nb - center)) - np.sum(
+                np.abs(target - center))
+            count += 1
+        I_col //= (count * blk_sz * blk_sz)
+        localDiff[r].append(
+            np.array([[I_row * I_row, I_row * I_col],
+                      [I_col * I_row, I_col * I_col]]))
+    return localDiff
+
+  """
+    add smooth constraint
+    """
+
+  def smooth(self, uvs, mvs):
+    sm_uvs = np.zeros(uvs.shape)
+    blk_sz = self.blk_sz
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        nb_uv = np.array([0.0, 0.0])
+        for i, j in {(r - 1, c), (r + 1, c), (r, c - 1), (r, c + 1)}:
+          if 0 <= i < self.num_row and 0 <= j < self.num_col:
+            nb_uv += uvs[i, j] / 6.0
+          else:
+            nb_uv += uvs[r, c] / 6.0
+        for i, j in {(r - 1, c - 1), (r - 1, c + 1), (r + 1, c - 1),
+                     (r + 1, c + 1)}:
+          if 0 <= i < self.num_row and 0 <= j < self.num_col:
+            nb_uv += uvs[i, j] / 12.0
+          else:
+            nb_uv += uvs[r, c] / 12.0
+        ssd_nb = self.block_dist(r, c, self.blk_sz * nb_uv)
+        mv = mvs[r, c]
+        ssd_mv = self.block_dist(r, c, mv)
+        alpha = (ssd_nb - ssd_mv) / (ssd_mv + 1e-6)
+        M = alpha * self.localDiff[r][c]
+        P = M + np.identity(2)
+        inv_P = LA.inv(P)
+        sm_uvs[r, c] = np.dot(inv_P, nb_uv) + np.dot(
+            np.matmul(inv_P, M), mv / blk_sz)
+    return sm_uvs
+
+  def block_matching(self):
+    self.search.motion_field_estimation()
+
+  def motion_field_estimation(self):
+    self.localDiff = self.getRefLocalDiff(self.search.mf)
+    #get matching results
+    mvs = self.search.mf
+    #add smoothness constraint
+    uvs = mvs / self.blk_sz
+    for _ in xrange(self.max_iter):
+      uvs = self.smooth(uvs, mvs)
+    self.mf = uvs * self.blk_sz
+
+
+"""Search & Smooth Model with Fixed Weights"""
+
+
+class SearchSmoothFix(MotionEST):
+  """
+    Constructor:
+        cur_f: current frame
+        ref_f: reference frame
+        blk_sz: block size
+        wnd_size: search window size
+        beta: neigbor loss weight
+        max_iter: maximum number of iterations
+        metric: metric to compare the blocks distrotion
+    """
+
+  def __init__(self, cur_f, ref_f, blk_size, search, beta, max_iter=100):
+    self.search = search
+    self.max_iter = max_iter
+    self.beta = beta
+    super(SearchSmoothFix, self).__init__(cur_f, ref_f, blk_size)
+
+  """
+    get local diffiencial of refernce
+    """
+
+  def getRefLocalDiff(self, mvs):
+    m, n = self.num_row, self.num_col
+    localDiff = [[] for _ in xrange(m)]
+    blk_sz = self.blk_sz
+    for r in xrange(m):
+      for c in xrange(n):
+        I_row = 0
+        I_col = 0
+        #get ssd surface
+        count = 0
+        center = self.cur_yuv[r * blk_sz:(r + 1) * blk_sz,
+                              c * blk_sz:(c + 1) * blk_sz, 0]
+        ty = np.clip(r * blk_sz + int(mvs[r, c, 0]), 0, self.height - blk_sz)
+        tx = np.clip(c * blk_sz + int(mvs[r, c, 1]), 0, self.width - blk_sz)
+        target = self.ref_yuv[ty:ty + blk_sz, tx:tx + blk_sz, 0]
+        for y, x in {(ty - blk_sz, tx), (ty + blk_sz, tx)}:
+          if 0 <= y < self.height - blk_sz and 0 <= x < self.width - blk_sz:
+            nb = self.ref_yuv[y:y + blk_sz, x:x + blk_sz, 0]
+            I_row += np.sum(np.abs(nb - center)) - np.sum(
+                np.abs(target - center))
+            count += 1
+        I_row //= (count * blk_sz * blk_sz)
+        count = 0
+        for y, x in {(ty, tx - blk_sz), (ty, tx + blk_sz)}:
+          if 0 <= y < self.height - blk_sz and 0 <= x < self.width - blk_sz:
+            nb = self.ref_yuv[y:y + blk_sz, x:x + blk_sz, 0]
+            I_col += np.sum(np.abs(nb - center)) - np.sum(
+                np.abs(target - center))
+            count += 1
+        I_col //= (count * blk_sz * blk_sz)
+        localDiff[r].append(
+            np.array([[I_row * I_row, I_row * I_col],
+                      [I_col * I_row, I_col * I_col]]))
+    return localDiff
+
+  """
+    add smooth constraint
+    """
+
+  def smooth(self, uvs, mvs):
+    sm_uvs = np.zeros(uvs.shape)
+    blk_sz = self.blk_sz
+    for r in xrange(self.num_row):
+      for c in xrange(self.num_col):
+        nb_uv = np.array([0.0, 0.0])
+        for i, j in {(r - 1, c), (r + 1, c), (r, c - 1), (r, c + 1)}:
+          if 0 <= i < self.num_row and 0 <= j < self.num_col:
+            nb_uv += uvs[i, j] / 6.0
+          else:
+            nb_uv += uvs[r, c] / 6.0
+        for i, j in {(r - 1, c - 1), (r - 1, c + 1), (r + 1, c - 1),
+                     (r + 1, c + 1)}:
+          if 0 <= i < self.num_row and 0 <= j < self.num_col:
+            nb_uv += uvs[i, j] / 12.0
+          else:
+            nb_uv += uvs[r, c] / 12.0
+        mv = mvs[r, c] / blk_sz
+        M = self.localDiff[r][c]
+        P = M + self.beta * np.identity(2)
+        inv_P = LA.inv(P)
+        sm_uvs[r, c] = np.dot(inv_P, self.beta * nb_uv) + np.dot(
+            np.matmul(inv_P, M), mv)
+    return sm_uvs
+
+  def block_matching(self):
+    self.search.motion_field_estimation()
+
+  def motion_field_estimation(self):
+    #get local structure
+    self.localDiff = self.getRefLocalDiff(self.search.mf)
+    #get matching results
+    mvs = self.search.mf
+    #add smoothness constraint
+    uvs = mvs / self.blk_sz
+    for _ in xrange(self.max_iter):
+      uvs = self.smooth(uvs, mvs)
+    self.mf = uvs * self.blk_sz

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/MotionEST/Util.py

@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+# coding: utf-8
+import numpy as np
+import numpy.linalg as LA
+import matplotlib.pyplot as plt
+from scipy.ndimage import filters
+from PIL import Image, ImageDraw
+
+
+def MSE(blk1, blk2):
+  return np.mean(
+      LA.norm(
+          np.array(blk1, dtype=np.int) - np.array(blk2, dtype=np.int), axis=2))
+
+
+def drawMF(img, blk_sz, mf):
+  img_rgba = img.convert('RGBA')
+  mf_layer = Image.new(mode='RGBA', size=img_rgba.size, color=(0, 0, 0, 0))
+  draw = ImageDraw.Draw(mf_layer)
+  width = img_rgba.size[0]
+  height = img_rgba.size[1]
+  num_row = height // blk_sz
+  num_col = width // blk_sz
+  for i in xrange(num_row):
+    left = (0, i * blk_sz)
+    right = (width, i * blk_sz)
+    draw.line([left, right], fill=(0, 0, 255, 255))
+  for j in xrange(num_col):
+    up = (j * blk_sz, 0)
+    down = (j * blk_sz, height)
+    draw.line([up, down], fill=(0, 0, 255, 255))
+  for i in xrange(num_row):
+    for j in xrange(num_col):
+      center = (j * blk_sz + 0.5 * blk_sz, i * blk_sz + 0.5 * blk_sz)
+      """mf[i,j][0] is the row shift and mf[i,j][1] is the column shift In PIL coordinates, head[0] is x (column shift) and head[1] is y (row shift)."""
+      head = (center[0] + mf[i, j][1], center[1] + mf[i, j][0])
+      draw.line([center, head], fill=(255, 0, 0, 255))
+  return Image.alpha_composite(img_rgba, mf_layer)

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/genY4M/genY4M.py

@@ -0,0 +1,76 @@
+import argparse
+from os import listdir, path
+from PIL import Image
+import sys
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--frame_path", default="../data/frame/", type=str)
+parser.add_argument("--frame_rate", default="25:1", type=str)
+parser.add_argument("--interlacing", default="Ip", type=str)
+parser.add_argument("--pix_ratio", default="0:0", type=str)
+parser.add_argument("--color_space", default="4:2:0", type=str)
+parser.add_argument("--output", default="output.y4m", type=str)
+
+
+def generate(args, frames):
+  if len(frames) == 0:
+    return
+  #sort the frames based on the frame index
+  frames = sorted(frames, key=lambda x: x[0])
+  #convert the frames to YUV form
+  frames = [f.convert("YCbCr") for _, f in frames]
+  #write the header
+  header = "YUV4MPEG2 W%d H%d F%s %s A%s" % (frames[0].width, frames[0].height,
+                                             args.frame_rate, args.interlacing,
+                                             args.pix_ratio)
+  cs = args.color_space.split(":")
+  header += " C%s%s%s\n" % (cs[0], cs[1], cs[2])
+  #estimate the sample step based on subsample value
+  subsamples = [int(c) for c in cs]
+  r_step = [1, int(subsamples[2] == 0) + 1, int(subsamples[2] == 0) + 1]
+  c_step = [1, 4 // subsamples[1], 4 // subsamples[1]]
+  #write in frames
+  with open(args.output, "wb") as y4m:
+    y4m.write(header)
+    for f in frames:
+      y4m.write("FRAME\n")
+      px = f.load()
+      for k in xrange(3):
+        for i in xrange(0, f.height, r_step[k]):
+          for j in xrange(0, f.width, c_step[k]):
+            yuv = px[j, i]
+            y4m.write(chr(yuv[k]))
+
+
+if __name__ == "__main__":
+  args = parser.parse_args()
+  frames = []
+  frames_mv = []
+  for filename in listdir(args.frame_path):
+    name, ext = filename.split(".")
+    if ext == "png":
+      name_parse = name.split("_")
+      idx = int(name_parse[-1])
+      img = Image.open(path.join(args.frame_path, filename))
+      if name_parse[-2] == "mv":
+        frames_mv.append((idx, img))
+      else:
+        frames.append((idx, img))
+  if len(frames) == 0:
+    print("No frames in directory: " + args.frame_path)
+    sys.exit()
+  print("----------------------Y4M Info----------------------")
+  print("width:  %d" % frames[0][1].width)
+  print("height: %d" % frames[0][1].height)
+  print("#frame: %d" % len(frames))
+  print("frame rate: %s" % args.frame_rate)
+  print("interlacing: %s" % args.interlacing)
+  print("pixel ratio: %s" % args.pix_ratio)
+  print("color space: %s" % args.color_space)
+  print("----------------------------------------------------")
+
+  print("Generating ...")
+  generate(args, frames)
+  if len(frames_mv) != 0:
+    args.output = args.output.replace(".y4m", "_mv.y4m")
+    generate(args, frames_mv)

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/BVH.pde

@@ -0,0 +1,163 @@
+/*
+ *AABB bounding box
+ *Bouding Volume Hierarchy
+ */
+class BoundingBox {
+  float min_x, min_y, min_z, max_x, max_y, max_z;
+  PVector center;
+  BoundingBox() {
+    min_x = Float.POSITIVE_INFINITY;
+    min_y = Float.POSITIVE_INFINITY;
+    min_z = Float.POSITIVE_INFINITY;
+    max_x = Float.NEGATIVE_INFINITY;
+    max_y = Float.NEGATIVE_INFINITY;
+    max_z = Float.NEGATIVE_INFINITY;
+    center = new PVector();
+  }
+  // build a bounding box for a triangle
+  void create(Triangle t) {
+    min_x = min(t.p1.x, min(t.p2.x, t.p3.x));
+    max_x = max(t.p1.x, max(t.p2.x, t.p3.x));
+
+    min_y = min(t.p1.y, min(t.p2.y, t.p3.y));
+    max_y = max(t.p1.y, max(t.p2.y, t.p3.y));
+
+    min_z = min(t.p1.z, min(t.p2.z, t.p3.z));
+    max_z = max(t.p1.z, max(t.p2.z, t.p3.z));
+    center.x = (max_x + min_x) / 2;
+    center.y = (max_y + min_y) / 2;
+    center.z = (max_z + min_z) / 2;
+  }
+  // merge two bounding boxes
+  void add(BoundingBox bbx) {
+    min_x = min(min_x, bbx.min_x);
+    min_y = min(min_y, bbx.min_y);
+    min_z = min(min_z, bbx.min_z);
+
+    max_x = max(max_x, bbx.max_x);
+    max_y = max(max_y, bbx.max_y);
+    max_z = max(max_z, bbx.max_z);
+    center.x = (max_x + min_x) / 2;
+    center.y = (max_y + min_y) / 2;
+    center.z = (max_z + min_z) / 2;
+  }
+  // get bounding box center axis value
+  float getCenterAxisValue(int axis) {
+    if (axis == 1) {
+      return center.x;
+    } else if (axis == 2) {
+      return center.y;
+    }
+    // when axis == 3
+    return center.z;
+  }
+  // check if a ray is intersected with the bounding box
+  boolean intersect(Ray r) {
+    float tmin, tmax;
+    if (r.dir.x >= 0) {
+      tmin = (min_x - r.ori.x) * (1.0f / r.dir.x);
+      tmax = (max_x - r.ori.x) * (1.0f / r.dir.x);
+    } else {
+      tmin = (max_x - r.ori.x) * (1.0f / r.dir.x);
+      tmax = (min_x - r.ori.x) * (1.0f / r.dir.x);
+    }
+
+    float tymin, tymax;
+    if (r.dir.y >= 0) {
+      tymin = (min_y - r.ori.y) * (1.0f / r.dir.y);
+      tymax = (max_y - r.ori.y) * (1.0f / r.dir.y);
+    } else {
+      tymin = (max_y - r.ori.y) * (1.0f / r.dir.y);
+      tymax = (min_y - r.ori.y) * (1.0f / r.dir.y);
+    }
+
+    if (tmax < tymin || tymax < tmin) {
+      return false;
+    }
+
+    tmin = tmin < tymin ? tymin : tmin;
+    tmax = tmax > tymax ? tymax : tmax;
+
+    float tzmin, tzmax;
+    if (r.dir.z >= 0) {
+      tzmin = (min_z - r.ori.z) * (1.0f / r.dir.z);
+      tzmax = (max_z - r.ori.z) * (1.0f / r.dir.z);
+    } else {
+      tzmin = (max_z - r.ori.z) * (1.0f / r.dir.z);
+      tzmax = (min_z - r.ori.z) * (1.0f / r.dir.z);
+    }
+    if (tmax < tzmin || tmin > tzmax) {
+      return false;
+    }
+    return true;
+  }
+}
+// Bounding Volume Hierarchy
+class BVH {
+  // Binary Tree
+  BVH left, right;
+  BoundingBox overall_bbx;
+  ArrayList<Triangle> mesh;
+  BVH(ArrayList<Triangle> mesh) {
+    this.mesh = mesh;
+    overall_bbx = new BoundingBox();
+    left = null;
+    right = null;
+    int mesh_size = this.mesh.size();
+    if (mesh_size <= 1) {
+      return;
+    }
+    // random select an axis
+    int axis = int(random(100)) % 3 + 1;
+    // build bounding box and save the selected center component
+    float[] axis_values = new float[mesh_size];
+    for (int i = 0; i < mesh_size; i++) {
+      Triangle t = this.mesh.get(i);
+      overall_bbx.add(t.bbx);
+      axis_values[i] = t.bbx.getCenterAxisValue(axis);
+    }
+    // find the median value of selected center component as pivot
+    axis_values = sort(axis_values);
+    float pivot;
+    if (mesh_size % 2 == 1) {
+      pivot = axis_values[mesh_size / 2];
+    } else {
+      pivot =
+          0.5f * (axis_values[mesh_size / 2 - 1] + axis_values[mesh_size / 2]);
+    }
+    // Build left node and right node by partitioning the mesh based on triangle
+    // bounding box center component value
+    ArrayList<Triangle> left_mesh = new ArrayList<Triangle>();
+    ArrayList<Triangle> right_mesh = new ArrayList<Triangle>();
+    for (int i = 0; i < mesh_size; i++) {
+      Triangle t = this.mesh.get(i);
+      if (t.bbx.getCenterAxisValue(axis) < pivot) {
+        left_mesh.add(t);
+      } else if (t.bbx.getCenterAxisValue(axis) > pivot) {
+        right_mesh.add(t);
+      } else if (left_mesh.size() < right_mesh.size()) {
+        left_mesh.add(t);
+      } else {
+        right_mesh.add(t);
+      }
+    }
+    left = new BVH(left_mesh);
+    right = new BVH(right_mesh);
+  }
+  // check if a ray intersect with current volume
+  boolean intersect(Ray r, float[] param) {
+    if (mesh.size() == 0) {
+      return false;
+    }
+    if (mesh.size() == 1) {
+      Triangle t = mesh.get(0);
+      return t.intersect(r, param);
+    }
+    if (!overall_bbx.intersect(r)) {
+      return false;
+    }
+    boolean left_res = left.intersect(r, param);
+    boolean right_res = right.intersect(r, param);
+    return left_res || right_res;
+  }
+}

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/Camera.pde

@@ -0,0 +1,138 @@
+class Camera {
+  // camera's field of view
+  float fov;
+  // camera's position, look at point and axis
+  PVector pos, center, axis;
+  PVector init_pos, init_center, init_axis;
+  float move_speed;
+  float rot_speed;
+  Camera(float fov, PVector pos, PVector center, PVector axis) {
+    this.fov = fov;
+    this.pos = pos;
+    this.center = center;
+    this.axis = axis;
+    this.axis.normalize();
+    move_speed = 0.001;
+    rot_speed = 0.01 * PI;
+    init_pos = pos.copy();
+    init_center = center.copy();
+    init_axis = axis.copy();
+  }
+
+  Camera copy() {
+    Camera cam = new Camera(fov, pos.copy(), center.copy(), axis.copy());
+    return cam;
+  }
+
+  PVector project(PVector pos) {
+    PVector proj = MatxVec3(getCameraMat(), PVector.sub(pos, this.pos));
+    proj.x = (float)height / 2.0 * proj.x / proj.z / tan(fov / 2.0f);
+    proj.y = (float)height / 2.0 * proj.y / proj.z / tan(fov / 2.0f);
+    proj.z = proj.z;
+    return proj;
+  }
+
+  float[] getCameraMat() {
+    float[] mat = new float[9];
+    PVector dir = PVector.sub(center, pos);
+    dir.normalize();
+    PVector left = dir.cross(axis);
+    left.normalize();
+    // processing camera system does not follow right hand rule
+    mat[0] = -left.x;
+    mat[1] = -left.y;
+    mat[2] = -left.z;
+    mat[3] = axis.x;
+    mat[4] = axis.y;
+    mat[5] = axis.z;
+    mat[6] = dir.x;
+    mat[7] = dir.y;
+    mat[8] = dir.z;
+
+    return mat;
+  }
+
+  void run() {
+    PVector dir, left;
+    if (mousePressed) {
+      float angleX = (float)mouseX / width * PI - PI / 2;
+      float angleY = (float)mouseY / height * PI - PI;
+      PVector diff = PVector.sub(center, pos);
+      float radius = diff.mag();
+      pos.x = radius * sin(angleY) * sin(angleX) + center.x;
+      pos.y = radius * cos(angleY) + center.y;
+      pos.z = radius * sin(angleY) * cos(angleX) + center.z;
+      dir = PVector.sub(center, pos);
+      dir.normalize();
+      PVector up = new PVector(0, 1, 0);
+      left = up.cross(dir);
+      left.normalize();
+      axis = dir.cross(left);
+      axis.normalize();
+    }
+
+    if (keyPressed) {
+      switch (key) {
+        case 'w':
+          dir = PVector.sub(center, pos);
+          dir.normalize();
+          pos = PVector.add(pos, PVector.mult(dir, move_speed));
+          center = PVector.add(center, PVector.mult(dir, move_speed));
+          break;
+        case 's':
+          dir = PVector.sub(center, pos);
+          dir.normalize();
+          pos = PVector.sub(pos, PVector.mult(dir, move_speed));
+          center = PVector.sub(center, PVector.mult(dir, move_speed));
+          break;
+        case 'a':
+          dir = PVector.sub(center, pos);
+          dir.normalize();
+          left = axis.cross(dir);
+          left.normalize();
+          pos = PVector.add(pos, PVector.mult(left, move_speed));
+          center = PVector.add(center, PVector.mult(left, move_speed));
+          break;
+        case 'd':
+          dir = PVector.sub(center, pos);
+          dir.normalize();
+          left = axis.cross(dir);
+          left.normalize();
+          pos = PVector.sub(pos, PVector.mult(left, move_speed));
+          center = PVector.sub(center, PVector.mult(left, move_speed));
+          break;
+        case 'r':
+          dir = PVector.sub(center, pos);
+          dir.normalize();
+          float[] mat = getRotationMat3x3(rot_speed, dir.x, dir.y, dir.z);
+          axis = MatxVec3(mat, axis);
+          axis.normalize();
+          break;
+        case 'b':
+          pos = init_pos.copy();
+          center = init_center.copy();
+          axis = init_axis.copy();
+          break;
+        case '+': move_speed *= 2.0f; break;
+        case '-': move_speed /= 2.0; break;
+        case CODED:
+          if (keyCode == UP) {
+            pos = PVector.add(pos, PVector.mult(axis, move_speed));
+            center = PVector.add(center, PVector.mult(axis, move_speed));
+          } else if (keyCode == DOWN) {
+            pos = PVector.sub(pos, PVector.mult(axis, move_speed));
+            center = PVector.sub(center, PVector.mult(axis, move_speed));
+          }
+      }
+    }
+  }
+  void open() {
+    perspective(fov, float(width) / height, 1e-6, 1e5);
+    camera(pos.x, pos.y, pos.z, center.x, center.y, center.z, axis.x, axis.y,
+           axis.z);
+  }
+  void close() {
+    ortho(-width, 0, -height, 0);
+    camera(0, 0, 0, 0, 0, 1, 0, 1, 0);
+  }
+}

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/MotionField.pde

@@ -0,0 +1,102 @@
+class MotionField {
+  int block_size;
+  ArrayList<PVector> motion_field;
+  MotionField(int block_size) {
+    this.block_size = block_size;
+    motion_field = new ArrayList<PVector>();
+  }
+
+  void update(Camera last_cam, Camera current_cam, PointCloud point_cloud,
+              BVH bvh) {
+    // clear motion field
+    motion_field = new ArrayList<PVector>();
+    int r_num = height / block_size, c_num = width / block_size;
+    for (int i = 0; i < r_num * c_num; i++)
+      motion_field.add(new PVector(0, 0, 0));
+    // estimate motion vector of each point in point cloud
+    for (int i = 0; i < point_cloud.size(); i++) {
+      PVector p = point_cloud.getPosition(i);
+      PVector p0 = current_cam.project(p);
+      PVector p1 = last_cam.project(p);
+      int row = int((p0.y + height / 2.0f) / block_size);
+      int col = int((p0.x + width / 2.0f) / block_size);
+      if (row >= 0 && row < r_num && col >= 0 && col < c_num) {
+        PVector accu = motion_field.get(row * c_num + col);
+        accu.x += p1.x - p0.x;
+        accu.y += p1.y - p0.y;
+        accu.z += 1;
+      }
+    }
+    // if some blocks do not have point, then use ray tracing to see if they are
+    // in triangles
+    for (int i = 0; i < r_num; i++)
+      for (int j = 0; j < c_num; j++) {
+        PVector accu = motion_field.get(i * c_num + j);
+        if (accu.z > 0) {
+          continue;
+        }
+        // use the center of the block to generate view ray
+        float cx = j * block_size + block_size / 2.0f - width / 2.0f;
+        float cy = i * block_size + block_size / 2.0f - height / 2.0f;
+        float cz = 0.5f * height / tan(current_cam.fov / 2.0f);
+        PVector dir = new PVector(cx, cy, cz);
+        float[] camMat = current_cam.getCameraMat();
+        dir = MatxVec3(transpose3x3(camMat), dir);
+        dir.normalize();
+        Ray r = new Ray(current_cam.pos, dir);
+        // ray tracing
+        float[] param = new float[4];
+        param[0] = Float.POSITIVE_INFINITY;
+        if (bvh.intersect(r, param)) {
+          PVector p = new PVector(param[1], param[2], param[3]);
+          PVector p0 = current_cam.project(p);
+          PVector p1 = last_cam.project(p);
+          accu.x += p1.x - p0.x;
+          accu.y += p1.y - p0.y;
+          accu.z += 1;
+        }
+      }
+    // estimate the motion vector of each block
+    for (int i = 0; i < r_num * c_num; i++) {
+      PVector mv = motion_field.get(i);
+      if (mv.z > 0) {
+        motion_field.set(i, new PVector(mv.x / mv.z, mv.y / mv.z, 0));
+      } else  // there is nothing in the block, use -1 to mark it.
+      {
+        motion_field.set(i, new PVector(0.0, 0.0, -1));
+      }
+    }
+  }
+
+  void render() {
+    int r_num = height / block_size, c_num = width / block_size;
+    for (int i = 0; i < r_num; i++)
+      for (int j = 0; j < c_num; j++) {
+        PVector mv = motion_field.get(i * c_num + j);
+        float ox = j * block_size + 0.5f * block_size;
+        float oy = i * block_size + 0.5f * block_size;
+        stroke(255, 0, 0);
+        line(ox, oy, ox + mv.x, oy + mv.y);
+      }
+  }
+
+  void save(String path) {
+    int r_num = height / block_size;
+    int c_num = width / block_size;
+    String[] mvs = new String[r_num];
+    for (int i = 0; i < r_num; i++) {
+      mvs[i] = "";
+      for (int j = 0; j < c_num; j++) {
+        PVector mv = motion_field.get(i * c_num + j);
+        if (mv.z != -1) {
+          mvs[i] += str(mv.x) + "," + str(mv.y);
+        } else  // there is nothing
+        {
+          mvs[i] += "-,-";
+        }
+        if (j != c_num - 1) mvs[i] += ";";
+      }
+    }
+    saveStrings(path, mvs);
+  }
+}

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/PointCloud.pde

@@ -0,0 +1,138 @@
+class PointCloud {
+  ArrayList<PVector> points;  // array to save points
+  IntList point_colors;       // array to save points color
+  PVector cloud_mass;
+  float[] depth;
+  boolean[] real;
+  PointCloud() {
+    // initialize
+    points = new ArrayList<PVector>();
+    point_colors = new IntList();
+    cloud_mass = new PVector(0, 0, 0);
+    depth = new float[width * height];
+    real = new boolean[width * height];
+  }
+
+  void generate(PImage rgb_img, PImage depth_img, Transform trans) {
+    if (depth_img.width != width || depth_img.height != height ||
+        rgb_img.width != width || rgb_img.height != height) {
+      println("rgb and depth file dimension should be same with window size");
+      exit();
+    }
+    // clear depth and real
+    for (int i = 0; i < width * height; i++) {
+      depth[i] = 0;
+      real[i] = false;
+    }
+    for (int v = 0; v < height; v++)
+      for (int u = 0; u < width; u++) {
+        // get depth value (red channel)
+        color depth_px = depth_img.get(u, v);
+        depth[v * width + u] = depth_px & 0x0000FFFF;
+        if (int(depth[v * width + u]) != 0) {
+          real[v * width + u] = true;
+        }
+        point_colors.append(rgb_img.get(u, v));
+      }
+    for (int v = 0; v < height; v++)
+      for (int u = 0; u < width; u++) {
+        if (int(depth[v * width + u]) == 0) {
+          interpolateDepth(v, u);
+        }
+        // add transformed pixel as well as pixel color to the list
+        PVector pos = trans.transform(u, v, int(depth[v * width + u]));
+        points.add(pos);
+        // accumulate z value
+        cloud_mass = PVector.add(cloud_mass, pos);
+      }
+  }
+  void fillInDepthAlongPath(float d, Node node) {
+    node = node.parent;
+    while (node != null) {
+      int i = node.row;
+      int j = node.col;
+      if (depth[i * width + j] == 0) {
+        depth[i * width + j] = d;
+      }
+      node = node.parent;
+    }
+  }
+  // interpolate
+  void interpolateDepth(int row, int col) {
+    if (row < 0 || row >= height || col < 0 || col >= width ||
+        int(depth[row * width + col]) != 0) {
+      return;
+    }
+    ArrayList<Node> queue = new ArrayList<Node>();
+    queue.add(new Node(row, col, null));
+    boolean[] visited = new boolean[width * height];
+    for (int i = 0; i < width * height; i++) visited[i] = false;
+    visited[row * width + col] = true;
+    // Using BFS to Find the Nearest Neighbor
+    while (queue.size() > 0) {
+      // pop
+      Node node = queue.get(0);
+      queue.remove(0);
+      int i = node.row;
+      int j = node.col;
+      // if current position have a real depth
+      if (depth[i * width + j] != 0 && real[i * width + j]) {
+        fillInDepthAlongPath(depth[i * width + j], node);
+        break;
+      } else {
+        // search unvisited 8 neighbors
+        for (int r = max(0, i - 1); r < min(height, i + 2); r++) {
+          for (int c = max(0, j - 1); c < min(width, j + 2); c++) {
+            if (!visited[r * width + c]) {
+              visited[r * width + c] = true;
+              queue.add(new Node(r, c, node));
+            }
+          }
+        }
+      }
+    }
+  }
+  // get point cloud size
+  int size() { return points.size(); }
+  // get ith position
+  PVector getPosition(int i) {
+    if (i >= points.size()) {
+      println("point position: index " + str(i) + " exceeds");
+      exit();
+    }
+    return points.get(i);
+  }
+  // get ith color
+  color getColor(int i) {
+    if (i >= point_colors.size()) {
+      println("point color: index " + str(i) + " exceeds");
+      exit();
+    }
+    return point_colors.get(i);
+  }
+  // get cloud center
+  PVector getCloudCenter() {
+    if (points.size() > 0) {
+      return PVector.div(cloud_mass, points.size());
+    }
+    return new PVector(0, 0, 0);
+  }
+  // merge two clouds
+  void merge(PointCloud point_cloud) {
+    for (int i = 0; i < point_cloud.size(); i++) {
+      points.add(point_cloud.getPosition(i));
+      point_colors.append(point_cloud.getColor(i));
+    }
+    cloud_mass = PVector.add(cloud_mass, point_cloud.cloud_mass);
+  }
+}
+
+class Node {
+  int row, col;
+  Node parent;
+  Node(int row, int col, Node parent) {
+    this.row = row;
+    this.col = col;
+    this.parent = parent;
+  }
+}

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/Ray_Tracing.pde

@@ -0,0 +1,61 @@
+// Triangle
+class Triangle {
+  // position
+  PVector p1, p2, p3;
+  // color
+  color c1, c2, c3;
+  BoundingBox bbx;
+  Triangle(PVector p1, PVector p2, PVector p3, color c1, color c2, color c3) {
+    this.p1 = p1;
+    this.p2 = p2;
+    this.p3 = p3;
+    this.c1 = c1;
+    this.c2 = c2;
+    this.c3 = c3;
+    bbx = new BoundingBox();
+    bbx.create(this);
+  }
+  // check to see if a ray intersects with the triangle
+  boolean intersect(Ray r, float[] param) {
+    PVector p21 = PVector.sub(p2, p1);
+    PVector p31 = PVector.sub(p3, p1);
+    PVector po1 = PVector.sub(r.ori, p1);
+
+    PVector dxp31 = r.dir.cross(p31);
+    PVector po1xp21 = po1.cross(p21);
+    float denom = p21.dot(dxp31);
+    float t = p31.dot(po1xp21) / denom;
+    float alpha = po1.dot(dxp31) / denom;
+    float beta = r.dir.dot(po1xp21) / denom;
+
+    boolean res = t > 0 && alpha > 0 && alpha < 1 && beta > 0 && beta < 1 &&
+                  alpha + beta < 1;
+    // depth test
+    if (res && t < param[0]) {
+      param[0] = t;
+      param[1] = alpha * p1.x + beta * p2.x + (1 - alpha - beta) * p3.x;
+      param[2] = alpha * p1.y + beta * p2.y + (1 - alpha - beta) * p3.y;
+      param[3] = alpha * p1.z + beta * p2.z + (1 - alpha - beta) * p3.z;
+    }
+    return res;
+  }
+  void render() {
+    beginShape(TRIANGLES);
+    fill(c1);
+    vertex(p1.x, p1.y, p1.z);
+    fill(c2);
+    vertex(p2.x, p2.y, p2.z);
+    fill(c3);
+    vertex(p3.x, p3.y, p3.z);
+    endShape();
+  }
+}
+// Ray
+class Ray {
+  // origin and direction
+  PVector ori, dir;
+  Ray(PVector ori, PVector dir) {
+    this.ori = ori;
+    this.dir = dir;
+  }
+}

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/Scene.pde

@@ -0,0 +1,59 @@
+class Scene {
+  PointCloud point_cloud;
+  ArrayList<Triangle> mesh;
+  BVH bvh;
+  MotionField motion_field;
+  Camera last_cam;
+  Camera current_cam;
+  int frame_count;
+
+  Scene(Camera camera, PointCloud point_cloud, MotionField motion_field) {
+    this.point_cloud = point_cloud;
+    this.motion_field = motion_field;
+    mesh = new ArrayList<Triangle>();
+    for (int v = 0; v < height - 1; v++)
+      for (int u = 0; u < width - 1; u++) {
+        PVector p1 = point_cloud.getPosition(v * width + u);
+        PVector p2 = point_cloud.getPosition(v * width + u + 1);
+        PVector p3 = point_cloud.getPosition((v + 1) * width + u + 1);
+        PVector p4 = point_cloud.getPosition((v + 1) * width + u);
+        color c1 = point_cloud.getColor(v * width + u);
+        color c2 = point_cloud.getColor(v * width + u + 1);
+        color c3 = point_cloud.getColor((v + 1) * width + u + 1);
+        color c4 = point_cloud.getColor((v + 1) * width + u);
+        mesh.add(new Triangle(p1, p2, p3, c1, c2, c3));
+        mesh.add(new Triangle(p3, p4, p1, c3, c4, c1));
+      }
+    bvh = new BVH(mesh);
+    last_cam = camera.copy();
+    current_cam = camera;
+    frame_count = 0;
+  }
+
+  void run() {
+    last_cam = current_cam.copy();
+    current_cam.run();
+    motion_field.update(last_cam, current_cam, point_cloud, bvh);
+    frame_count += 1;
+  }
+
+  void render(boolean show_motion_field) {
+    // build mesh
+    current_cam.open();
+    noStroke();
+    for (int i = 0; i < mesh.size(); i++) {
+      Triangle t = mesh.get(i);
+      t.render();
+    }
+    if (show_motion_field) {
+      current_cam.close();
+      motion_field.render();
+    }
+  }
+
+  void save(String path) { saveFrame(path + "_" + str(frame_count) + ".png"); }
+
+  void saveMotionField(String path) {
+    motion_field.save(path + "_" + str(frame_count) + ".txt");
+  }
+}

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/Transform.pde

@@ -0,0 +1,82 @@
+class Transform {
+  float[] inv_rot;  // inverse of rotation matrix
+  PVector inv_mov;  // inverse of movement vector
+  float focal;      // the focal distacne of real camera
+  int w, h;         // the width and height of the frame
+  float normalier;  // nomalization factor of depth
+  Transform(float tx, float ty, float tz, float qx, float qy, float qz,
+            float qw, float fov, int w, int h, float normalier) {
+    // currently, we did not use the info of real camera's position and
+    // quaternion maybe we will use it in the future when combine all frames
+    float[] rot = quaternion2Mat3x3(qx, qy, qz, qw);
+    inv_rot = transpose3x3(rot);
+    inv_mov = new PVector(-tx, -ty, -tz);
+    this.focal = 0.5f * h / tan(fov / 2.0);
+    this.w = w;
+    this.h = h;
+    this.normalier = normalier;
+  }
+
+  PVector transform(int i, int j, float d) {
+    // transfer from camera view to world view
+    float z = d / normalier;
+    float x = (i - w / 2.0f) * z / focal;
+    float y = (j - h / 2.0f) * z / focal;
+    return new PVector(x, y, z);
+  }
+}
+
+// get rotation matrix by using rotation axis and angle
+float[] getRotationMat3x3(float angle, float ax, float ay, float az) {
+  float[] mat = new float[9];
+  float c = cos(angle);
+  float s = sin(angle);
+  mat[0] = c + ax * ax * (1 - c);
+  mat[1] = ax * ay * (1 - c) - az * s;
+  mat[2] = ax * az * (1 - c) + ay * s;
+  mat[3] = ay * ax * (1 - c) + az * s;
+  mat[4] = c + ay * ay * (1 - c);
+  mat[5] = ay * az * (1 - c) - ax * s;
+  mat[6] = az * ax * (1 - c) - ay * s;
+  mat[7] = az * ay * (1 - c) + ax * s;
+  mat[8] = c + az * az * (1 - c);
+  return mat;
+}
+
+// get rotation matrix by using quaternion
+float[] quaternion2Mat3x3(float qx, float qy, float qz, float qw) {
+  float[] mat = new float[9];
+  mat[0] = 1 - 2 * qy * qy - 2 * qz * qz;
+  mat[1] = 2 * qx * qy - 2 * qz * qw;
+  mat[2] = 2 * qx * qz + 2 * qy * qw;
+  mat[3] = 2 * qx * qy + 2 * qz * qw;
+  mat[4] = 1 - 2 * qx * qx - 2 * qz * qz;
+  mat[5] = 2 * qy * qz - 2 * qx * qw;
+  mat[6] = 2 * qx * qz - 2 * qy * qw;
+  mat[7] = 2 * qy * qz + 2 * qx * qw;
+  mat[8] = 1 - 2 * qx * qx - 2 * qy * qy;
+  return mat;
+}
+
+// tranpose a 3x3 matrix
+float[] transpose3x3(float[] mat) {
+  float[] Tmat = new float[9];
+  for (int i = 0; i < 3; i++)
+    for (int j = 0; j < 3; j++) {
+      Tmat[i * 3 + j] = mat[j * 3 + i];
+    }
+  return Tmat;
+}
+
+// multiply a matrix with vector
+PVector MatxVec3(float[] mat, PVector v) {
+  float[] vec = v.array();
+  float[] res = new float[3];
+  for (int i = 0; i < 3; i++) {
+    res[i] = 0.0f;
+    for (int j = 0; j < 3; j++) {
+      res[i] += mat[i * 3 + j] * vec[j];
+    }
+  }
+  return new PVector(res[0], res[1], res[2]);
+}

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/Util.pde

@@ -0,0 +1,28 @@
+// show grids
+void showGrids(int block_size) {
+  ortho(-width, 0, -height, 0);
+  camera(0, 0, 0, 0, 0, 1, 0, 1, 0);
+  stroke(0, 0, 255);
+  for (int i = 0; i < height; i += block_size) {
+    line(0, i, width, i);
+  }
+  for (int i = 0; i < width; i += block_size) {
+    line(i, 0, i, height);
+  }
+}
+
+// save the point clould information
+void savePointCloud(PointCloud point_cloud, String file_name) {
+  String[] positions = new String[point_cloud.points.size()];
+  String[] colors = new String[point_cloud.points.size()];
+  for (int i = 0; i < point_cloud.points.size(); i++) {
+    PVector point = point_cloud.getPosition(i);
+    color point_color = point_cloud.getColor(i);
+    positions[i] = str(point.x) + ' ' + str(point.y) + ' ' + str(point.z);
+    colors[i] = str(((point_color >> 16) & 0xFF) / 255.0) + ' ' +
+                str(((point_color >> 8) & 0xFF) / 255.0) + ' ' +
+                str((point_color & 0xFF) / 255.0);
+  }
+  saveStrings(file_name + "_pos.txt", positions);
+  saveStrings(file_name + "_color.txt", colors);
+}

Libraries/webm_mem_playback/libvpx/tools/3D-Reconstruction/sketch_3D_reconstruction/sketch_3D_reconstruction.pde

@@ -0,0 +1,74 @@
+/*The dataset is from
+ *Computer Vision Group
+ *TUM Department of Informatics Technical
+ *University of Munich
+ *https://vision.in.tum.de/data/datasets/rgbd-dataset/download#freiburg1_xyz
+ */
+Scene scene;
+void setup() {
+  size(640, 480, P3D);
+  // default settings
+  int frame_no = 0;            // frame number
+  float fov = PI / 3;          // field of view
+  int block_size = 8;          // block size
+  float normalizer = 5000.0f;  // normalizer
+  // initialize
+  PointCloud point_cloud = new PointCloud();
+  // synchronized rgb, depth and ground truth
+  String head = "../data/";
+  String[] rgb_depth_gt = loadStrings(head + "rgb_depth_groundtruth.txt");
+  // read in rgb and depth image file paths as well as corresponding camera
+  // posiiton and quaternion
+  String[] info = split(rgb_depth_gt[frame_no], ' ');
+  String rgb_path = head + info[1];
+  String depth_path = head + info[3];
+  float tx = float(info[7]), ty = float(info[8]),
+        tz = float(info[9]);  // real camera position
+  float qx = float(info[10]), qy = float(info[11]), qz = float(info[12]),
+        qw = float(info[13]);  // quaternion
+
+  // build transformer
+  Transform trans =
+      new Transform(tx, ty, tz, qx, qy, qz, qw, fov, width, height, normalizer);
+  PImage rgb = loadImage(rgb_path);
+  PImage depth = loadImage(depth_path);
+  // generate point cloud
+  point_cloud.generate(rgb, depth, trans);
+  // initialize camera
+  Camera camera = new Camera(fov, new PVector(0, 0, 0), new PVector(0, 0, 1),
+                             new PVector(0, 1, 0));
+  // initialize motion field
+  MotionField motion_field = new MotionField(block_size);
+  // initialize scene
+  scene = new Scene(camera, point_cloud, motion_field);
+}
+boolean inter = false;
+void draw() {
+  background(0);
+  // run camera dragged mouse to rotate camera
+  // w: go forward
+  // s: go backward
+  // a: go left
+  // d: go right
+  // up arrow: go up
+  // down arrow: go down
+  //+ increase move speed
+  //- decrease move speed
+  // r: rotate the camera
+  // b: reset to initial position
+  scene.run();  // true: make interpolation; false: do not make
+                // interpolation
+  if (keyPressed && key == 'o') {
+    inter = true;
+  }
+  scene.render(
+      false);  // true: turn on motion field; false: turn off motion field
+  // save frame with no motion field
+  scene.save("../data/frame/raw");
+  background(0);
+  scene.render(true);
+  showGrids(scene.motion_field.block_size);
+  // save frame with motion field
+  scene.save("../data/frame/raw_mv");
+  scene.saveMotionField("../data/frame/mv");
+}