diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..88213e6
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+# C extensions
+*.so
+*.o
+*.pyc
+
+utils/fast_render/build
diff --git a/README.md b/README.md
index f7bfcd3..832beca 100644
--- a/README.md
+++ b/README.md
@@ -173,6 +173,12 @@ What you can custom:
 
    you can refer to  [c++ version](https://github.com/YadiraF/face3d/blob/master/face3d/mesh/render.py). 
 
+   Another alternative is to use the C++ already in this repo, integrated with ctypes. You must install both CMake and a C++ compiler, then run the following command: 
+   ```
+   sh build.sh
+   ```
+   This code only implements the hotpath of rendering the face depth map. The speedup from a 256x256 image was from 2681ms to 1.8ms, quite a difference.
+
    c. other parts like detecting face, writing obj
 
    the best way is to rewrite them in c++.
diff --git a/build.sh b/build.sh
new file mode 100644
index 0000000..8c624ac
--- /dev/null
+++ b/build.sh
@@ -0,0 +1,5 @@
+cd utils/fast_render
+mkdir build
+cd build
+cmake ..
+make
diff --git a/utils/fast_render/CMakeLists.txt b/utils/fast_render/CMakeLists.txt
new file mode 100644
index 0000000..7ee5f2e
--- /dev/null
+++ b/utils/fast_render/CMakeLists.txt
@@ -0,0 +1,17 @@
+cmake_minimum_required(VERSION 3.5)
+
+project(fast_render VERSION 0.0.1)
+
+set_property(GLOBAL PROPERTY USE_FOLDERS ON)
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wfatal-errors -O2")
+
+include_directories(
+    "src/"
+)
+
+file(GLOB all_fast_render_src
+    "src/*.cpp"
+)
+
+add_library(fast_render SHARED ${all_fast_render_src})
diff --git a/utils/fast_render/__init__.py b/utils/fast_render/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/utils/fast_render/fast_render.py b/utils/fast_render/fast_render.py
new file mode 100644
index 0000000..a39de4d
--- /dev/null
+++ b/utils/fast_render/fast_render.py
@@ -0,0 +1,42 @@
+import ctypes
+import numpy as np
+from numpy.ctypeslib import ndpointer
+
+lib = ctypes.CDLL('./utils/fast_render/build/libfast_render.so')
+
+c_render_texture_loop = lib.render_texture_loop
+c_render_texture_loop.argtypes = [
+    ndpointer(ctypes.c_double, flags='C_CONTIGUOUS'), # tri_depth
+    ctypes.c_int,
+    ndpointer(ctypes.c_double, flags='C_CONTIGUOUS'), # tri_tex
+    ctypes.c_int,
+    ctypes.c_int,
+    ndpointer(ctypes.c_int, flags='C_CONTIGUOUS'), # triangles
+    ctypes.c_int,
+    ctypes.c_int,
+    ndpointer(ctypes.c_double, flags='C_CONTIGUOUS'), # vertices
+    ctypes.c_int,
+    ctypes.c_int,
+    ndpointer(ctypes.c_double, flags='C_CONTIGUOUS'), # depth_buffer
+    ctypes.c_int,
+    ctypes.c_int,
+    ndpointer(ctypes.c_double, flags='C_CONTIGUOUS'), # image
+    ctypes.c_int,
+    ctypes.c_int,
+]
+
+
+def render_texture_loop(tri_depth, tri_tex, triangles, vertices, depth_buffer, image):
+    if len(image.shape) == 2:
+        image_channels = 1
+    else:
+        image_channels = image.shape[2]
+
+    c_render_texture_loop(
+            tri_depth, tri_depth.shape[0],
+            tri_tex, tri_tex.shape[1], tri_tex.shape[0],
+            triangles, triangles.shape[1], triangles.shape[0],
+            vertices, vertices.shape[1], vertices.shape[0],
+            depth_buffer, depth_buffer.shape[1], depth_buffer.shape[0],
+            image, image.shape[1], image.shape[0], image_channels
+        )
diff --git a/utils/fast_render/src/render.cpp b/utils/fast_render/src/render.cpp
new file mode 100644
index 0000000..ecdbe08
--- /dev/null
+++ b/utils/fast_render/src/render.cpp
@@ -0,0 +1,64 @@
+#include "render.h"
+
+
+bool PointInTriangle(double p_x, double p_y, double p0_x, double p1_x, double p2_x, double p0_y, double p1_y,  double p2_y) {
+    // Original code from here: https://github.com/SebLague/Gamedev-Maths/blob/master/PointInTriangle.cs
+    double s1 = p2_y - p0_y;
+    double s2 = p2_x - p0_x;
+    double s3 = p1_y - p0_y;
+    double s4 = p_y - p0_y;
+
+    double w1 = (p0_x * s1 + s4 * s2 - p_x * s1) / (s3 * s2 - (p1_x-p0_x) * s1);
+    double w2 = (s4- w1 * s3) / s1;
+    return w1 >= 0 && w2 >= 0 && (w1 + w2) <= 1;
+}
+
+void render_texture_loop(const double *tri_depth, int tri_depth_height,
+        const double *tri_tex, int tri_tex_width, int tri_tex_height,
+        const int *triangles, int triangles_width, int triangles_height,
+        const double *vertices, int vertices_width, int vertices_height,
+        double *depth_buffer, int depth_buffer_width, int depth_buffer_height,
+        double *image, int image_width, int image_height, int image_channels)
+{
+    for (int i=0; i<triangles_width; i++) {
+        // 3 vertex indices
+        int tri1 = triangles[0*triangles_width+i];
+        int tri2 = triangles[1*triangles_width+i];
+        int tri3 = triangles[2*triangles_width+i];
+
+        double v1_u = vertices[0*vertices_width+tri1];
+        double v2_u = vertices[0*vertices_width+tri2];
+        double v3_u = vertices[0*vertices_width+tri3];
+
+        double v1_v = vertices[1*vertices_width+tri1];
+        double v2_v = vertices[1*vertices_width+tri2];
+        double v3_v = vertices[1*vertices_width+tri3];
+
+
+        // the inner bounding box
+        int umin = std::max(ceil(std::min(std::min(v1_u,v2_u),v3_u)), 0.0);
+        int umax = std::min(floor(std::max(std::max(v1_u,v2_u),v3_u)), double(image_width-1));
+
+        int vmin = std::max(ceil(std::min(std::min(v1_v,v2_v),v3_v)), 0.0);
+        int vmax = std::min(floor(std::max(std::max(v1_v,v2_v),v3_v)), double(image_height-1));
+
+        if (umax<umin || vmax<vmin)
+            continue;
+
+        for (int u = umin; u<umax+1; u++) {
+            for (int v = vmin; v<vmax+1; v++) {
+                if (tri_depth[i] > depth_buffer[v*depth_buffer_width + u]) {
+                    if (PointInTriangle((double)u, (double)v, v1_u, v2_u, v3_u, v1_v, v2_v, v3_v)) {
+                        depth_buffer[v*depth_buffer_width + u] = tri_depth[i];
+                        for (int aux=0; aux<tri_tex_height; aux++) {
+                            for (int c=0; c<image_channels; c++) {
+                                image[v*image_width*image_channels + u*image_channels + c] = tri_tex[aux*tri_tex_width + i];
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
diff --git a/utils/fast_render/src/render.h b/utils/fast_render/src/render.h
new file mode 100644
index 0000000..1cafee4
--- /dev/null
+++ b/utils/fast_render/src/render.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include <string.h>
+#include <stdio.h>
+#include <math.h> // ceil
+#include <algorithm> // min max
+
+extern "C" {
+    void render_texture_loop(
+            const double *tri_depth, int tri_depth_height,
+            const double *tri_tex, int tri_tex_width, int tri_tex_height,
+            const int *triangles, int triangles_width, int triangles_height,
+            const double *vertices, int vertices_width, int vertices_height,
+            double *depth_buffer, int depth_buffer_width, int depth_buffer_height,
+            double *image, int image_width, int image_height, int image_channels
+        );
+}
diff --git a/utils/render.py b/utils/render.py
index 348ae62..14fccd0 100644
--- a/utils/render.py
+++ b/utils/render.py
@@ -3,6 +3,13 @@
 Mail: fengyao@sjtu.edu.cn
 '''
 import numpy as np
+try:
+    from utils.fast_render import fast_render
+    using_fast_render = True
+except BaseException as e:
+    print("INFO - You can run 'sh build.sh' to build a faster "
+            "rendering function for the depthmap images")
+    using_fast_render = False
 
 def isPointInTri(point, tri_points):
     ''' Judge whether the point is in the triangle
@@ -99,24 +106,31 @@ def render_texture(vertices, colors, triangles, h, w, c = 3):
     tri_depth = (vertices[2, triangles[0,:]] + vertices[2,triangles[1,:]] + vertices[2, triangles[2,:]])/3. 
     tri_tex = (colors[:, triangles[0,:]] + colors[:,triangles[1,:]] + colors[:, triangles[2,:]])/3.
 
-    for i in range(triangles.shape[1]):
-        tri = triangles[:, i] # 3 vertex indices
+    triangles = np.ascontiguousarray(triangles)
+    vertices = np.ascontiguousarray(vertices)
 
-        # the inner bounding box
-        umin = max(int(np.ceil(np.min(vertices[0,tri]))), 0)
-        umax = min(int(np.floor(np.max(vertices[0,tri]))), w-1)
+    if using_fast_render:
+        fast_render.render_texture_loop(tri_depth, tri_tex, triangles, vertices, depth_buffer, image)
+    else:
+        for i in range(triangles.shape[1]):
+            tri = triangles[:, i] # 3 vertex indices
 
-        vmin = max(int(np.ceil(np.min(vertices[1,tri]))), 0)
-        vmax = min(int(np.floor(np.max(vertices[1,tri]))), h-1)
+            # the inner bounding box
+            umin = max(int(np.ceil(np.min(vertices[0,tri]))), 0)
+            umax = min(int(np.floor(np.max(vertices[0,tri]))), w-1)
 
-        if umax<umin or vmax<vmin:
-            continue
+            vmin = max(int(np.ceil(np.min(vertices[1,tri]))), 0)
+            vmax = min(int(np.floor(np.max(vertices[1,tri]))), h-1)
+
+            if umax<umin or vmax<vmin:
+                continue
+
+            for u in range(umin, umax+1):
+                for v in range(vmin, vmax+1):
+                    if tri_depth[i] > depth_buffer[v, u] and isPointInTri([u,v], vertices[:2, tri]): 
+                        depth_buffer[v, u] = tri_depth[i]
+                        image[v, u, :] = tri_tex[:, i]
 
-        for u in range(umin, umax+1):
-            for v in range(vmin, vmax+1):
-                if tri_depth[i] > depth_buffer[v, u] and isPointInTri([u,v], vertices[:2, tri]): 
-                    depth_buffer[v, u] = tri_depth[i]
-                    image[v, u, :] = tri_tex[:, i]
     return image