Fix building and 'reconstructing' mean texture from dataset

src/aam.py

@@ -5,7 +5,8 @@ import cv2
 
 # local imports
 import pca
-from utils.generate_head_texture import fill_triangle, get_colors_triangle
+from utils.generate_head_texture import fill_triangle, get_colors_triangle, \
+    get_row_colors_triangle
 import utils.triangles as tu
 
 logging.basicConfig(level=logging.INFO,
@@ -66,47 +67,40 @@ def build_shape_feature_vectors(files, get_points, flattened=False):
     points = get_points(files)
 
     if flattened:
-        points = pca.flatten_feature_vectors(points)
+        points = pca.flatten_feature_vectors(points, dim=0)
 
     return points
 
 
-def sample_from_triangles(image, points2d, triangles, n_samples=80):
-    all_triangles = []
-    h, w, c = image.shape
+def sample_from_triangles(src, points2d_src, points2d_dst, triangles):
+    # texture = np.asarray(texture, dtype=np.uint8).reshape((-1, 3))
+    triangles_pixels = []
+    pixels = 0
 
     for tri in triangles:
-        p1 = points2d[tri[0]]
-        p2 = points2d[tri[1]]
-        p3 = points2d[tri[2]]
-
-        bary_centric_range = np.linspace(0, 1, num=n_samples)
-        pixels = np.full((n_samples * n_samples, 3), fill_value=-1, dtype=np.int)
-        L = np.zeros((3, 1))
-
-        for s_i, s in enumerate(bary_centric_range):
-            for t_i, t in enumerate(bary_centric_range):
-                # make sure the coordinates are inside the triangle
-                if s + t <= 1:
-                    # build lambda's
-                    L[0] = s
-                    L[1] = t
-                    L[2] = 1 - s - t
+        src_p1, src_p2, src_p3 = points2d_src[tri]
+        dst_p1, dst_p2, dst_p3 = points2d_dst[tri]
+
+        dst = get_row_colors_triangle(
+            src,
+            src_p1[0], src_p1[1],
+            src_p2[0], src_p2[1],
+            src_p3[0], src_p3[1],
+            dst_p1[0], dst_p1[1],
+            dst_p2[0], dst_p2[1],
+            dst_p3[0], dst_p3[1]
+        )
 
-                    # cartesian x, y coordinates inside the triangle
-                    cart_x, cart_y, _ = tu.barycentric2cartesian(p1, p2, p3, L)
-                    pixels[s_i * n_samples + t_i, :] = image[cart_y, cart_x, :]
+        pixels += dst.flatten().shape[0]
 
-                    # cv2.circle(b, tuple([cart_x, cart_y]), 1, color=(0, 255, 100))
+        triangles_pixels.extend(dst.flatten())
 
-        all_triangles.append(pixels[np.where(pixels >= 0)])
+    result = np.asarray(triangles_pixels, dtype=np.uint8)
 
-    return np.asarray(all_triangles, dtype=np.uint8)
+    return result
 
 
-def build_texture_feature_vectors(files, get_image_with_landmarks, triangles,
-        flattened=True):
-    mean_texture = []
+def build_texture_feature_vectors(files, get_image_with_shape, mean_shape, triangles):
     """
     Args:
         files (list): list files
@@ -116,15 +110,21 @@ def build_texture_feature_vectors(files, get_image_with_landmarks, triangles,
     Returns:
         list: list of feature vectors
     """
+    mean_texture = []
 
-    for i, f in enumerate(files[:1]):
-        image, landmarks = get_image_with_landmarks(f)
+    mean_shape_scaled = mean_shape.reshape((58, 2))
+    mean_shape_scaled[:, 0] = mean_shape_scaled[:, 0] * 640
+    mean_shape_scaled[:, 1] = mean_shape_scaled[:, 1] * 480
+
+    for i, f in enumerate(files):
+        image, shape = get_image_with_shape(f)
         h, w, c = image.shape
-        landmarks[:, 0] = landmarks[:, 0] * w
-        landmarks[:, 1] = landmarks[:, 1] * h
+
+        shape[:, 0] = shape[:, 0] * w
+        shape[:, 1] = shape[:, 1] * h
 
         triangles_colors = sample_from_triangles(
-            image, landmarks, triangles, n_samples=80
+            image, shape, mean_shape_scaled, triangles
             )
 
         mean_texture.append(triangles_colors)
@@ -135,11 +135,8 @@ def build_texture_feature_vectors(files, get_image_with_landmarks, triangles,
 
         # if k == 27:
         #     break
-
     mean_texture = np.asarray(mean_texture)
-
-    if flattened:
-        mean_texture = pca.flatten_feature_vectors(mean_texture)
+    #mean_texture = pca.flatten_feature_vectors(mean_texture)
 
     return mean_texture
 

src/main.py

@@ -84,16 +84,17 @@ def save_pca_model_texture(args):
     assert args.model_shape_file, '--model_texture_file needs to be provided to save the pca model'
     assert args.model_texture_file, '--model_texture_file needs to be provided to save the pca model'
 
-    Vt, mean_values, triangles = pca.load(args.model_shape_file)
+    Vt, s, mean_shape, triangles = pca.load(args.model_shape_file)
 
     textures = aam.build_texture_feature_vectors(
-        args.files, imm.get_imm_image_with_landmarks, triangles, flattened=True
+        args.files, imm.get_imm_image_with_landmarks, mean_shape, triangles
     )
 
     mean_texture = aam.get_mean(textures)
-    _, _, Vt = pca.pca(textures, mean_texture)
+    _, s, Vt = pca.pca(textures, mean_texture)
+
+    pca.save(Vt, s, mean_texture, triangles, args.model_texture_file)
 
-    pca.save(Vt, mean_texture, triangles, args.model_texture_file)
     logger.info('texture pca model saved in %s', args.model_texture_file)
 
 
@@ -122,12 +123,12 @@ def save_pca_model_shape(args):
 
     mean_values = aam.get_mean(points)
 
-    _, _, Vt = pca.pca(points, mean_values)
+    _, s, Vt = pca.pca(points, mean_values)
 
     mean_xy = mean_values.reshape((-1, 2))
     triangles = aam.get_triangles(mean_xy[:, 0], mean_xy[:, 1])
 
-    pca.save(Vt, mean_values, triangles, args.model_shape_file)
+    pca.save(Vt, s, mean_values, triangles, args.model_shape_file)
     logger.info('shape pca model saved in %s', args.model_shape_file + '_shape')
 
 
@@ -164,14 +165,25 @@ def show_pca_model(args):
 
     from utils.triangles import draw_shape, draw_texture
 
-    Vt_shape, mean_values_shape, triangles = pca.load(args.model_shape_file)
-    Vt_texture, mean_values_texture, _ = pca.load(args.model_texture_file)
+    Vt_shape, s, mean_values_shape, triangles = pca.load(args.model_shape_file)
+    Vt_texture, s_texture, mean_values_texture, _ = pca.load(args.model_texture_file)
+
+    # calculate n_components which captures 90 percent of the variance
+    total = s_texture.sum()
+    subtotal = 0.0
+    i = 0
+
+    while (subtotal * 100.0) / total <= 90.0:
+        subtotal += s_texture[i]
+        i += 1
+
+    n_components = i
 
     image = np.full((480, 640, 3), fill_value=255, dtype=np.uint8)
 
-    immPoints = imm.IMMPoints(filename='data/imm_face_db/40-1m.asf')
-    input_image = immPoints.get_image()
-    input_points = immPoints.get_points()
+    imm_points = imm.IMMPoints(filename='data/imm_face_db/40-1m.asf')
+    input_image = imm_points.get_image()
+    input_points = imm_points.get_points()
     h, w, c = input_image.shape
 
     input_points[:, 0] = input_points[:, 0] * w
@@ -182,9 +194,9 @@ def show_pca_model(args):
     mean_values_shape[:, 1] = mean_values_shape[:, 1] * h
 
     while True:
-        draw_texture(input_image, image, input_points, mean_values_shape,
-                     mean_values_texture, triangles, n_samples=80)
-        draw_shape(image, mean_values_shape, triangles, multiply=False)
+        draw_texture(input_image, image, Vt_texture, input_points, mean_values_shape,
+                     mean_values_texture, triangles)
+        #draw_shape(image, mean_values_shape, triangles, multiply=False)
 
         cv2.imshow('input_image', input_image)
         cv2.imshow('image', image)

src/pca.py

@@ -36,7 +36,7 @@ def reconstruct(feature_vector, Vt, mean_values, n_components=10):
     return np.dot(Vt[:n_components].T, yk) + mean_values
 
 
-def save(Vt, mean_values, triangles, filename):
+def save(Vt, s, mean_values, triangles, filename):
     """
     Store the U, s, Vt and mean of all the asf datafiles given by the asf
     files.
@@ -52,7 +52,7 @@ def save(Vt, mean_values, triangles, filename):
         triangles = Vtm[2]
 
     """
-    saving = np.asarray([Vt, [mean_values], triangles])
+    saving = np.asarray([Vt, s, [mean_values], triangles])
     np.save(filename, saving)
 
 
@@ -76,13 +76,14 @@ def load(filename):
     Vtm = np.load(filename)
 
     Vt = Vtm[0]
-    mean_values = Vtm[1][0]
-    triangles = Vtm[2]
+    s = Vtm[1]
+    mean_values = Vtm[2][0]
+    triangles = Vtm[3]
 
-    return Vt, mean_values, triangles
+    return Vt, s, mean_values, triangles
 
 
-def flatten_feature_vectors(data):
+def flatten_feature_vectors(data, dim=0):
     """
     Flattens the feature vectors inside a ndarray
 
@@ -102,6 +103,7 @@ def flatten_feature_vectors(data):
 
     Args:
         data (numpy array): array of feature vectors
+        dim (int): dimension to flatten the data
 
     return:
         array: (numpy array): array flattened feature vectors
@@ -109,9 +111,9 @@ def flatten_feature_vectors(data):
     """
     flattened = []
 
-    rows, _, _ = data.shape
+    n = data.shape[dim]
 
-    for i in range(rows):
+    for i in range(n):
         flattened.append(np.ndarray.flatten(data[i]))
 
     return np.array(flattened)

src/utils/generate_head_texture.pyx

@@ -18,6 +18,23 @@ cdef inline float cross_product(int v1_x, int v1_y, int v2_x, int v2_y):
     return (v1_x * v2_y) - (v1_y * v2_x)
 
 
+cdef inline np.ndarray[double, ndim=1] cartesian2barycentric(
+        int r1_x, r1_y, int r2_x, int r2_y, int r3_x, int r3_y, int r_x, int r_y):
+    """
+    Given a triangle spanned by three cartesion points
+    r1, r2, r2, and point r, return the barycentric weights l1, l2, l3.
+
+    Returns:
+        ndarray (of dim 3) weights of the barycentric coordinates
+
+    """
+    a = np.array([[r1_x, r2_x, r3_x], [r1_y, r2_y, r3_y], [1, 1, 1]])
+    b = np.array([r_x, r_y, 1])
+
+    return np.linalg.solve(a, b)
+
+
+
 cdef inline np.ndarray[double, ndim=2] barycentric2cartesian(
             int x1, int x2, int x3, int y1, int y2, int y3,
             np.ndarray[long, ndim=2] matrix,
@@ -35,9 +52,10 @@ cdef inline np.ndarray[double, ndim=2] barycentric2cartesian(
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def fill_triangle(np.ndarray[unsigned char, ndim=3] src,
+def fill_triangle(np.ndarray[unsigned char, ndim=1] src,
                   np.ndarray[unsigned char, ndim=3] dst,
-                  int x1, int y1, int x2, int y2, int x3, int y3):
+                  int x1, int y1, int x2, int y2, int x3, int y3, int offset,
+                  int index):
     """
     Fill a triangle by applying the Barycentric Algorithm for deciding if a
     point lies inside or outside a triangle.
@@ -49,32 +67,42 @@ def fill_triangle(np.ndarray[unsigned char, ndim=3] src,
     cdef int y_min = min(y1, min(y2, y3))
     cdef int y_max = max(y1, max(y2, y3))
 
-    cdef int vs1_x = x2 - x1
-    cdef int vs1_y = y2 - y1
-    cdef int vs2_x = x3 - x1
-    cdef int vs2_y = y3 - y1
+    cdef np.ndarray L = np.zeros([3, 1], dtype=DTYPE_float32)
+    cdef np.ndarray matrix = np.full([3, 3], fill_value=1, dtype=DTYPE_int)
 
-    cdef float s
-    cdef float t
+    cdef np.ndarray src_loc = np.zeros([3, 1], dtype=DTYPE_float64)
+    cdef np.ndarray dst_loc = np.zeros([3, 1], dtype=DTYPE_float64)
 
-    for y in xrange(y_min, y_max):
-        for x in xrange(x_min, x_max):
-            q_x = x - x1
-            q_y = y - y1
+    cdef int w = x_max - x_min
+    cdef int h = y_max - y_min
+    cdef int new_offset
 
-            s = cross_product(q_x, q_y, vs2_x, vs2_y) / \
-                cross_product(vs1_x, vs1_y, vs2_x, vs2_y)
-            t = cross_product(vs1_x, vs1_y, q_x, q_y) / \
-                cross_product(vs1_x, vs1_y, vs2_x, vs2_y)
+    cdef np.ndarray src_reshaped = src[offset:offset + (w * h * 3)].reshape((h, w, 3))
 
-            if s >= 0 and t >= 0 and s + t <= 1:
-                dst[y, x, :] = src[y, x, :]
+    #print src_reshaped
+    #print '(', w, '*', h, '*', 3, ') * ', ' = ', offset
+
+    for j, y in enumerate(xrange(y_min, y_max)):
+        for i, x in enumerate(xrange(x_min, x_max)):
+            dst_loc = cartesian2barycentric(
+                x1, y1, x2, y2, x3, y3, x, y
+            )
+
+            s = dst_loc[0]
+            t = dst_loc[1]
+
+            # notice we have a soft margin of -0.00001, which makes sure there are no
+            # gaps due to rounding issues
+            if s >= -0.000000001 and t >= -0.000000001 and s + t <= 1.0:
+                dst[y, x, :] = src_reshaped[j, i, :]
+
+    new_offset = (w * h * 3)
+    return new_offset
 
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def get_colors_triangle(np.ndarray[unsigned char, ndim=3] src,
-                        np.ndarray[unsigned char, ndim=3] dst,
+def get_row_colors_triangle(np.ndarray[unsigned char, ndim=3] src,
                             int src_x1, int src_y1, int src_x2, int src_y2,
                             int src_x3, int src_y3,
                             int dst_x1, int dst_y1, int dst_x2, int dst_y2,
@@ -83,21 +111,32 @@ def get_colors_triangle(np.ndarray[unsigned char, ndim=3] src,
     Fill a triangle by applying the Barycentric Algorithm for deciding if a
     point lies inside or outside a triangle.
     """
-    cdef float s
-    cdef float t
+    cdef int x_min = min(dst_x1, min(dst_x2, dst_x3))
+    cdef int x_max = max(dst_x1, max(dst_x2, dst_x3))
+    cdef int y_min = min(dst_y1, min(dst_y2, dst_y3))
+    cdef int y_max = max(dst_y1, max(dst_y2, dst_y3))
 
     cdef np.ndarray L = np.zeros([3, 1], dtype=DTYPE_float32)
     cdef np.ndarray matrix = np.full([3, 3], fill_value=1, dtype=DTYPE_int)
 
     cdef np.ndarray src_loc = np.zeros([3, 1], dtype=DTYPE_float64)
     cdef np.ndarray dst_loc = np.zeros([3, 1], dtype=DTYPE_float64)
+    cdef np.ndarray dst = np.full(
+        [y_max - y_min, x_max - x_min, 3], fill_value=255, dtype=DTYPE_float64
+    )
+
+    for j, y in enumerate(xrange(y_min, y_max)):
+        for i, x in enumerate(xrange(x_min, x_max)):
+            dst_loc = cartesian2barycentric(
+                dst_x1, dst_y1, dst_x2, dst_y2, dst_x3, dst_y3, x, y
+            )
 
-    cdef np.ndarray bary_centric_range = np.linspace(0, 1, num=80)
+            s = dst_loc[0]
+            t = dst_loc[1]
 
-    # get a float value for every pixel
-    for s in bary_centric_range:
-        for t in bary_centric_range:
-            if s + t <= 1:
+            # notice we have a soft margin of -0.00001, which makes sure there are no
+            # gaps due to rounding issues
+            if s >= -0.000001 and t >= -0.000001 and s + t <= 1.0:
                 L[0] = s
                 L[1] = t
                 L[2] = 1 - s - t
@@ -109,11 +148,58 @@ def get_colors_triangle(np.ndarray[unsigned char, ndim=3] src,
                     L
                 )
 
-                dst_loc = barycentric2cartesian(
-                    dst_x1, dst_x2, dst_x3,
-                    dst_y1, dst_y2, dst_y3,
+                dst[j, i, :] = src[src_loc[1][0], src_loc[0][0], :]
+
+    return dst
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def get_colors_triangle(np.ndarray[unsigned char, ndim=3] src,
+                        np.ndarray[unsigned char, ndim=3] dst,
+                        int src_x1, int src_y1, int src_x2, int src_y2,
+                        int src_x3, int src_y3,
+                        int dst_x1, int dst_y1, int dst_x2, int dst_y2,
+                        int dst_x3, int dst_y3):
+    """
+    Fill a triangle by applying the Barycentric Algorithm for deciding if a
+    point lies inside or outside a triangle.
+    """
+    cdef int x_min = min(dst_x1, min(dst_x2, dst_x3))
+    cdef int x_max = max(dst_x1, max(dst_x2, dst_x3))
+    cdef int y_min = min(dst_y1, min(dst_y2, dst_y3))
+    cdef int y_max = max(dst_y1, max(dst_y2, dst_y3))
+
+    cdef float s
+    cdef float t
+
+    cdef np.ndarray L = np.zeros([3, 1], dtype=DTYPE_float32)
+    cdef np.ndarray matrix = np.full([3, 3], fill_value=1, dtype=DTYPE_int)
+
+    cdef np.ndarray src_loc = np.zeros([3, 1], dtype=DTYPE_float64)
+    cdef np.ndarray dst_loc = np.zeros([3, 1], dtype=DTYPE_float64)
+
+    for y in xrange(y_min, y_max):
+        for x in xrange(x_min, x_max):
+            dst_loc = cartesian2barycentric(
+                dst_x1, dst_y1, dst_x2, dst_y2, dst_x3, dst_y3, x, y
+            )
+
+            s = dst_loc[0]
+            t = dst_loc[1]
+
+            # notice we have a soft margin of -0.00001, which makes sure there are no
+            # gaps due to rounding issues
+            if s >= -0.000001 and t >= -0.000001 and s + t <= 1.0:
+                L[0] = s
+                L[1] = t
+                L[2] = 1 - s - t
+
+                src_loc = barycentric2cartesian(
+                    src_x1, src_x2, src_x3,
+                    src_y1, src_y2, src_y3,
                     matrix,
                     L
                 )
 
-                dst[dst_loc[1][0], dst_loc[0][0], :] = src[src_loc[1][0], src_loc[0][0], :]
+                dst[y, x, :] = src[src_loc[1][0], src_loc[0][0], :]

src/utils/triangles.py

@@ -67,20 +67,20 @@ def draw_shape(image, points, triangles, multiply=True):
         cv2.circle(image, tuple(p), 3, color=(0, 255, 100))
 
 
-def draw_texture(src, dest, points2d_src, points2d_dest, texture,
-                 triangles, multiply=True, n_samples=20):
-    texture = np.asarray(texture, dtype=np.uint8).reshape((-1, 3))
+def draw_texture(src, dest, Vt, points2d_src, points2d_dst, texture, triangles):
+    # texture = np.asarray(texture, dtype=np.uint8).reshape((-1, 3))
+    texture = np.asarray(texture, np.uint8)
+    offset = 0
 
     for t, tri in enumerate(triangles):
         src_p1, src_p2, src_p3 = points2d_src[tri]
-        dest_p1, dest_p2, dest_p3 = points2d_dest[tri]
-
-        get_colors_triangle(
-            src, dest,
-            src_p1[0], src_p1[1],
-            src_p2[0], src_p2[1],
-            src_p3[0], src_p3[1],
-            dest_p1[0], dest_p1[1],
-            dest_p2[0], dest_p2[1],
-            dest_p3[0], dest_p3[1]
+        dst_p1, dst_p2, dst_p3 = points2d_dst[tri]
+
+        offset += fill_triangle(
+            texture, dest,
+            dst_p1[0], dst_p1[1],
+            dst_p2[0], dst_p2[1],
+            dst_p3[0], dst_p3[1],
+            offset,
+            t
         )