diff --git a/get_proj.py b/get_proj.py
new file mode 100644
index 0000000..be6a2d1
--- /dev/null
+++ b/get_proj.py
@@ -0,0 +1,114 @@
+# https://blender.stackexchange.com/questions/38009/3x4-camera-matrix-from-blender-camera
+
+import bpy
+from mathutils import Matrix
+from mathutils import Vector
+
+#---------------------------------------------------------------
+# 3x4 P matrix from Blender camera
+#---------------------------------------------------------------
+
+# Build intrinsic camera parameters from Blender camera data
+#
+# See notes on this in 
+# blender.stackexchange.com/questions/15102/what-is-blenders-camera-projection-matrix-model
+def get_calibration_matrix_K_from_blender(camd):
+    f_in_mm = camd.lens
+    scene = bpy.context.scene
+    resolution_x_in_px = scene.render.resolution_x
+    resolution_y_in_px = scene.render.resolution_y
+    scale = scene.render.resolution_percentage / 100
+    sensor_width_in_mm = camd.sensor_width
+    sensor_height_in_mm = camd.sensor_height
+    pixel_aspect_ratio = scene.render.pixel_aspect_x / scene.render.pixel_aspect_y
+    if (camd.sensor_fit == 'VERTICAL'):
+        # the sensor height is fixed (sensor fit is horizontal), 
+        # the sensor width is effectively changed with the pixel aspect ratio
+        s_u = resolution_x_in_px * scale / sensor_width_in_mm / pixel_aspect_ratio 
+        s_v = resolution_y_in_px * scale / sensor_height_in_mm
+    else: # 'HORIZONTAL' and 'AUTO'
+        # the sensor width is fixed (sensor fit is horizontal), 
+        # the sensor height is effectively changed with the pixel aspect ratio
+        pixel_aspect_ratio = scene.render.pixel_aspect_x / scene.render.pixel_aspect_y
+        s_u = resolution_x_in_px * scale / sensor_width_in_mm
+        s_v = resolution_y_in_px * scale * pixel_aspect_ratio / sensor_height_in_mm
+    
+
+    # Parameters of intrinsic calibration matrix K
+    alpha_u = f_in_mm * s_u
+    alpha_v = f_in_mm * s_v
+    u_0 = resolution_x_in_px * scale / 2
+    v_0 = resolution_y_in_px * scale / 2
+    skew = 0 # only use rectangular pixels
+
+    K = Matrix(
+        ((alpha_u, skew,    u_0),
+        (    0  , alpha_v, v_0),
+        (    0  , 0,        1 )))
+    return K
+
+# Returns camera rotation and translation matrices from Blender.
+# 
+# There are 3 coordinate systems involved:
+#    1. The World coordinates: "world"
+#       - right-handed
+#    2. The Blender camera coordinates: "bcam"
+#       - x is horizontal
+#       - y is up
+#       - right-handed: negative z look-at direction
+#    3. The desired computer vision camera coordinates: "cv"
+#       - x is horizontal
+#       - y is down (to align to the actual pixel coordinates 
+#         used in digital images)
+#       - right-handed: positive z look-at direction
+def get_3x4_RT_matrix_from_blender(cam):
+    # bcam stands for blender camera
+    R_bcam2cv = Matrix(
+        ((1, 0,  0),
+         (0, -1, 0),
+         (0, 0, -1)))
+
+    # Transpose since the rotation is object rotation, 
+    # and we want coordinate rotation
+    # R_world2bcam = cam.rotation_euler.to_matrix().transposed()
+    # T_world2bcam = -1*R_world2bcam * location
+    #
+    # Use matrix_world instead to account for all constraints
+    location, rotation = cam.matrix_world.decompose()[0:2]
+    R_world2bcam = rotation.to_matrix().transposed()
+
+    # Convert camera location to translation vector used in coordinate changes
+    # T_world2bcam = -1*R_world2bcam*cam.location
+    # Use location from matrix_world to account for constraints:     
+    T_world2bcam = -1*R_world2bcam @ location
+
+    # Build the coordinate transform matrix from world to computer vision camera
+    # NOTE: Use * instead of @ here for older versions of Blender
+    # TODO: detect Blender version
+    R_world2cv = R_bcam2cv@R_world2bcam
+    T_world2cv = R_bcam2cv@T_world2bcam
+
+    # put into 3x4 matrix
+    RT = Matrix((
+        R_world2cv[0][:] + (T_world2cv[0],),
+        R_world2cv[1][:] + (T_world2cv[1],),
+        R_world2cv[2][:] + (T_world2cv[2],)
+         ))
+    return RT
+
+def get_3x4_P_matrix_from_blender(cam):
+    K = get_calibration_matrix_K_from_blender(cam.data)
+    RT = get_3x4_RT_matrix_from_blender(cam)
+    return K@RT, K, RT
+
+
+
+def run_script(scene):
+    # projection_matrix = scene.camera.matrix_world
+    projection_matrix, _, _ = get_3x4_P_matrix_from_blender(scene.camera)
+    with open('/home/damien/Documents/3A/projet-be/imgs/torus/matrices.txt', 'a') as f:
+        f.write(projection_matrix.__repr__() + '\n\n')
+
+f = open('/home/damien/Documents/3A/projet-be/imgs/torus/matrices.txt', 'w')
+f.close()
+bpy.app.handlers.frame_change_post.append(run_script)
diff --git a/main.py b/main.py
new file mode 100644
index 0000000..d81b802
--- /dev/null
+++ b/main.py
@@ -0,0 +1,30 @@
+import cv2
+import numpy as np
+
+from matrices_reader import *
+
+VOXEL_SIZE = 1e-3
+X_MIN, X_MAX = -2.0, 2.0
+Y_MIN, Y_MAX = -2.0, 2.0
+Z_MIN, Z_MAX = -2.0, 2.0
+
+# grid = [[[
+#     1 for z in np.arange(Z_MIN, Z_MAX, VOXEL_SIZE)
+#     ] for y in np.arange(Y_MIN, Y_MAX, VOXEL_SIZE)
+#     ] for x in np.arange(X_MIN, X_MAX, VOXEL_SIZE)
+# ]
+
+projection_matrices = matrices_reader('data/torus/matrices.txt')
+nb_frame = len(projection_matrices)
+point = np.array([1.0, 0.0, 0.0, 1.0])
+
+for k in range(nb_frame):
+    
+    proj_mat = projection_matrices[k]
+    cam_point = proj_mat @ point
+    cam_point /= cam_point[2]
+
+    frame = cv2.imread(f'data/torus/torus{k+1:04}.png')
+    cv2.circle(frame, (int(cam_point[0]), int(cam_point[1])), 2, (0, 0, 255))
+    cv2.imshow('Frame', frame)
+    cv2.waitKey(0)
\ No newline at end of file
diff --git a/matrices_reader.py b/matrices_reader.py
new file mode 100644
index 0000000..527cc08
--- /dev/null
+++ b/matrices_reader.py
@@ -0,0 +1,36 @@
+import re
+import numpy as np
+
+
+def matrices_reader(path: str) -> list[np.ndarray]:
+    """Read projection matrices.
+
+    Args:
+        path (str): path to matrices.txt
+
+    Returns:
+        list[np.ndarray]: list of projection matrix
+    """
+
+    with open(path, 'r') as f:
+        lines = f.readlines()
+
+    k = 0
+    world_matrices = []
+    while k+3 < len(lines):
+        # Match matrices one by one
+        mat_str = ""
+        for line in lines[k:k+4]:
+            mat_str += line
+        float_reg = r"(-|\d|\.|e)+"
+        res = re.search(
+            f"Matrix\(\(\(({float_reg}), ({float_reg}), ({float_reg}), ({float_reg})\),\n +\(({float_reg}), ({float_reg}), ({float_reg}), ({float_reg})\),\n\ +\(({float_reg}), ({float_reg}), ({float_reg}), ({float_reg})\)\)\)", mat_str)
+        
+        # Convert string to np.ndarray
+        values = [float(res.group(i)) for i in range(1,len(res.groups()) + 1, 2)]
+        world_mat = np.array([[values[4*i + j] for j in range(4)] for i in range(3)])
+        world_matrices.append(world_mat)
+        
+        k += 4
+    
+    return world_matrices[1:]
diff --git a/torus.blend b/torus.blend
new file mode 100644
index 0000000..0d1bfda
Binary files /dev/null and b/torus.blend differ