diff --git a/python/_project.pyx b/python/_project.pyx
index fe4f3a9..61ca0dc 100644
--- a/python/_project.pyx
+++ b/python/_project.pyx
@@ -9,6 +9,11 @@ DTYPE=np.float64
 
 __all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE","interp3d","interp2d"]
 
+cdef extern from "project_tools.hpp" namespace "":
+
+  DTYPE_t compute_projection(DTYPE_t *vertex_value, DTYPE_t *u, DTYPE_t *u0, DTYPE_t rho)
+
+
 @cython.boundscheck(False)
 @cython.cdivision(True)
 cdef DTYPE_t interp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y, 
@@ -520,49 +525,46 @@ cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]):
 
     return alpha_max
 
-#cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]):
-#    cdef DTYPE_t alpha_max
-#    cdef DTYPE_t tmp_a
-#    cdef DTYPE_t v[3], term[4]
-#    cdef int i, j, q
+@cython.boundscheck(False)
+@cython.cdivision(True)
+cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil:
+  cdef int i
+  cdef DTYPE_t s
+  
+  s = 0
+  
+  for i in xrange(q):
+    s += v[i]
+  return s
 
-#    alpha_max = 10.0 # A big number
+@cython.boundscheck(False)
+@cython.cdivision(True)
+cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]):
+    cdef DTYPE_t alpha_max
+    cdef DTYPE_t tmp_a
+    cdef DTYPE_t v[3], term[4]
+    cdef int i, j, q
 
-#    j = 0
-#    for i in range(3):
-#        if u[i] == 0.:
-#            continue
+    alpha_max = 10.0 # A big number
 
-#        if u[i] < 0:
-#            tmp_a = -u0[i]/u[i]
-#        else:
-#            tmp_a = (1-u0[i])/u[i]
+    j = 0
+    for i in range(3):
+        if u[i] == 0.:
+            continue
 
-#        if tmp_a < alpha_max:
-#            alpha_max = tmp_a
-#            j = i
+        if u[i] < 0:
+            tmp_a = -u0[i]/u[i]
+        else:
+            tmp_a = (1-u0[i])/u[i]
+
+        if tmp_a < alpha_max:
+            alpha_max = tmp_a
+            j = i
 
     # alpha_max is the integration length
-    # now we compute the integration of a trilinearly interpolated field
-    # There are four terms.
+    # we integrate between 0 and alpha_max (curvilinear coordinates)
     
-    
-    # First term
-#    term[0]= (u0[0]*u0[1]*u0[2])*sum(vertex_value)
-    
-    # Second term
-#    term[1] = 0
-#
-#    for q in range(3):
-#      for r in range(8):
-#        pass
-
-#    for i in range(3):
-#        u0[i] += u[i]*alpha_max
-
-#    r[0] = j
-
-#    return 0#alpha_max
+    return compute_projection(vertex_value, u, u0, alpha_max)
     
 @cython.boundscheck(False)
 def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
diff --git a/python/project_tool.hpp b/python/project_tool.hpp
new file mode 100644
index 0000000..9322b5f
--- /dev/null
+++ b/python/project_tool.hpp
@@ -0,0 +1,113 @@
+
+// Only in 3d
+
+template<typename T, typename ProdType>
+static T project_tool(T *vertex_value, T *u, T *u0)
+{
+  T ret0 = 0;
+  for (int i = 0; i < 8; i++)
+    {
+      int c[3] = { i & 1, (i>>1)&1, (i>>2)&1 };
+      int epsilon[3];
+      T ret = 0;
+      
+      for (int q = 0; q < 3; q++)
+        epsilon[q] = 2*c[q] - 1;
+
+      for (int q = 0; q < ProdType::numProducts; q++)
+        ret += ProdType::product(u, u0, epsilon, q);
+      ret *= vertex_value[i];
+      ret0 += ret;
+    }
+    
+  return ret0;
+}
+
+
+template<typename T>
+struct ProductTerm0
+{
+  static const int numProducts = 1;
+  
+  static T product(T *u, T *u0, int *epsilon, int q)
+  {
+    T a = 1;
+    
+    for (int r = 0; r < 3; r++)
+      a *= (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
+    return a;
+  }
+};
+
+template<typename T>
+struct ProductTerm1
+{
+  static const int numProducts = 3;
+  
+  static T product(T *u, T *u0, int *epsilon, int q)
+  {
+    T a = 1;
+    double G[3];
+    
+    for (int r = 0; r < 3; r++)
+      {
+        G[r] = (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
+      }
+
+    double F[3] = { G[0]*u[1]*u[2], u[0]*G[1]*u[2], u[0]*u[1]*G[2] };    
+
+    return F[q] * epsilon[q];
+  }
+};
+
+template<typename T>
+struct ProductTerm2
+{
+  static const int numProducts = 3;
+  
+  static T product(T *u, T *u0, int *epsilon, int q)
+  {
+    T a = 1;
+    double G[3];
+    
+    for (int r = 0; r < 3; r++)
+      {
+        G[r] = (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
+      }
+
+    double F[3] = { u[0]*G[1]*G[2], G[0]*u[1]*G[2], G[0]*G[1]*u[2] };    
+
+    return F[q] * epsilon[q];
+  }
+};
+
+
+template<typename T>
+struct ProductTerm3
+{
+  static const int numProducts = 1;
+  
+  static T product(T *u, T *u0, int *epsilon, int q)
+  {
+    T a = 1;
+    
+    return epsilon[0]*epsilon[1]*epsilon[2];
+  }
+};
+
+
+template<typename T>
+T compute_projection(T *vertex_value, T *u, T *u0, T rho)
+{
+  T ret;
+  
+  ret = project_tool<T, ProductTerm0<T> >(vertex_value, u, u0) * rho;
+  ret += project_tool<T, ProductTerm1<T> >(vertex_value, u, u0) * rho * rho / 2;
+  ret += project_tool<T, ProductTerm2<T> >(vertex_value, u, u0) * rho * rho * rho / 3;
+  ret += project_tool<T, ProductTerm3<T> >(vertex_value, u, u0) * rho * rho * rho * rho / 4;
+  return ret;
+}
+
+template
+double compute_projection(double *vertex_value, double *u, double *u0, double rho);
+