diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..a726827b905b85d2243032438c5200a7c9d29b06
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,8 @@
+car_models/*
+output/*
+mh_models/*
+backgrounds/
+.idea/
+
+*.pyc
+*.blend1
diff --git a/base.blend b/base.blend
new file mode 100644
index 0000000000000000000000000000000000000000..f75511705b22c7bbc7a69877a0e7ca6b358ee7db
Binary files /dev/null and b/base.blend differ
diff --git a/global_script.py b/global_script.py
new file mode 100644
index 0000000000000000000000000000000000000000..9bc51ca8c5e0fb4a86080725565ec7d6e586fc98
--- /dev/null
+++ b/global_script.py
@@ -0,0 +1,327 @@
+import random
+import glob
+import json
+import os
+import sys
+import shutil
+from collections import OrderedDict
+
+import bpy
+from mathutils import Vector
+
+print('#' * 30)
+
+WORKING_DIR = r"D:\Mingming\synthe_dripe"
+os.chdir(WORKING_DIR)
+sys.path.append(WORKING_DIR)
+
+from scripts import random_pose
+from scripts import utils
+from scripts import camera_proj
+from scripts import human
+
+import importlib
+
+importlib.reload(random_pose)
+importlib.reload(utils)
+importlib.reload(camera_proj)
+importlib.reload(human)
+
+from scripts.human import Human, HumanLoader, set_shrinkwraps
+
+
+def abs_path(rel_path):
+    return os.path.join(os.path.dirname(os.path.dirname(__file__)), rel_path)
+
+
+random.seed()
+
+C = bpy.context
+D = bpy.data
+
+# Clean scene
+try:
+    bpy.ops.object.mode_set(mode='OBJECT')
+    utils.unselect_all()
+except RuntimeError:
+    pass
+
+for col in D.collections:
+    D.collections.remove(col)
+
+for bpy_data_iter in (
+        D.objects,
+        D.meshes,
+        D.lights,
+        D.cameras,
+        D.armatures,
+        D.images,
+):
+    for id_data in bpy_data_iter:
+        bpy_data_iter.remove(id_data)
+
+for ob in D.objects:
+    C.scene.collection.objects.unlink(ob)
+
+# Import car
+car_collection = D.collections.new("Cars")
+C.scene.collection.children.link(car_collection)
+
+cars = []
+
+for src_path, name in {
+    r"car_models\suv_car\car.blend": 'SUV',
+    r"car_models\red_car\red.blend": 'Red',
+    r"car_models\pickup_car\pickup.blend": 'PickUp',
+    r"car_models\family_car\family_car.blend": 'Family',
+    r"car_models\coupe_car\coupe_car.blend": 'Coupe',
+    r"car_models\truck\truck_open.blend": 'Truck',
+
+}.items():
+    with D.libraries.load(abs_path(src_path)) as (data_from, data_to):
+        data_to.objects = data_from.objects
+
+    for obj in data_to.objects:
+        car_collection.objects.link(obj)
+
+    D.objects['Car'].name = name
+    cars.append(D.objects[name])
+car_picker = utils.Randomizer(cars)
+
+# import humans
+human_loader = HumanLoader('mh_models/exports')
+
+# Creation scene
+# add camera
+C.scene.render.engine = 'BLENDER_EEVEE'
+
+camera_data = D.cameras.new(name='Camera')
+camera_data.type = 'PERSP'
+camera_data.lens_unit = 'FOV'
+# camera_data.angle = utils.r(68)
+camera_data.angle = utils.r(68)
+
+C.scene.render.resolution_x = 640
+C.scene.render.resolution_y = 480
+
+camera_object = D.objects.new('Camera', camera_data)
+C.scene.collection.objects.link(camera_object)
+
+# camera_object.location = [-1.1, -0.21, 0.54]
+# camera_object.rotation_euler = utils.r([76, 12, -76])
+camera_object.location = [-.97, -0.11, 0.68]
+camera_object.rotation_euler = utils.r([72, 8, -82])
+
+# set background
+back_folder = abs_path("backgrounds")
+back_imgs = {}
+for key in ['night', 'day']:
+    list_imgs = glob.glob(os.path.join(back_folder, f'{key}_*'))
+    back_imgs[key] = []
+    for img in list_imgs:
+        img_name = os.path.basename(img)
+        bpy.ops.image.open(filepath=img, directory=back_folder,
+                           files=[{"name": img_name}], relative_path=False, show_multiview=False)
+        back_imgs[key].append(img_name)
+
+# Create image holder
+bpy.ops.import_image.to_plane(
+    directory=back_folder,
+    files=[{"name": "default_green.png"}],
+    shader='SHADELESS',
+    use_transparency=False,
+    offset=False,
+    height=round(10 / 2.25, 1),
+    align_axis="X-"
+)
+
+image_holder = C.active_object
+image_holder.name = 'Image_holder'
+image_holder.location = (4, 1.5, 1.3)
+image_holder.rotation_euler.z = utils.r(-100)
+image_holder.active_material.shadow_method = 'NONE'
+
+# C.scene.node_tree.nodes["Rotate"].inputs[1].default_value = camera_object.rotation_euler[1] - utils.r(5)
+
+# add light
+sun_collection = D.collections.new("Sun")
+C.scene.collection.children.link(sun_collection)
+
+light_params = {
+    'day': {
+        'energy_bounds': (10, 35),
+        'back_color_bounds': (0.4, 0.65),
+        'sun_color': (1, 1, 1)
+    },
+    'night': {
+        'energy_bounds': (5, 15),
+        'back_color_bounds': (0.05, 0.15),
+        'sun_color': (1, 0.5, 0.2)
+    }
+
+}
+
+light_data = D.lights.new(name="Sun", type='SUN')
+light_data.energy = 20
+light_object = D.objects.new(name="Sun", object_data=light_data)
+
+sun_collection.objects.link(light_object)
+light_object.location = (5, 0, 3)
+
+# Sun position
+C.scene.sun_pos_properties.usage_mode = 'NORMAL'
+C.scene.sun_pos_properties.sun_object = light_object
+C.scene.sun_pos_properties.object_collection = sun_collection
+C.scene.sun_pos_properties.object_collection_type = 'DIURNAL'
+C.scene.sun_pos_properties.co_parser = "41°22′14″N 2°09′00″E"
+C.scene.sun_pos_properties.sun_distance = 3
+C.scene.sun_pos_properties.use_day_of_year = True
+C.scene.sun_pos_properties.year = 2022
+C.scene.sun_pos_properties.day_of_year = 182
+
+fp = abs_path(r"output")
+if os.path.isdir(fp):
+    shutil.rmtree(fp)
+os.mkdir(fp)
+
+fp_img = os.path.join(fp, 'images')
+fp_ann_2D = os.path.join(fp, 'annots_2D')
+fp_ann_3D = os.path.join(fp, 'annots_3D')
+os.mkdir(fp_img)
+os.mkdir(fp_ann_2D)
+os.mkdir(fp_ann_3D)
+
+frame_rate = 25
+nb_scene = 2
+nb_pose = 1
+human_loader.max_len = min(human_loader.max_len, nb_scene)
+ratio_conf_man = int(nb_scene / len(human_loader.human_paths))
+
+scenes_ids = OrderedDict()
+
+for sc in range(nb_scene):
+    ## Random car
+    car = car_picker()
+    ## Random personne
+    if ratio_conf_man < 1:
+        if not sc % 10:
+            human_loader.load_next()
+        man = human_loader(car=car)
+    else:
+        if not sc % ratio_conf_man:
+            man = human_loader.next(car=car)
+        else:
+            human.set_bounds(man, car)
+
+    human_path = human_loader.paths[man]
+    scenes_ids.setdefault(human_path, -1)
+    scenes_ids[human_path] += 1
+
+    ## Random time
+    C.scene.sun_pos_properties.north_offset = utils.r(random.randint(-179, 180))
+    time_day = random.randint(0, 23)
+    if 6 <= time_day <= 21:
+        day_night = 'day'
+        C.scene.sun_pos_properties.time = time_day
+    else:
+        day_night = 'night'
+        C.scene.sun_pos_properties.time = (time_day + 12) % 24
+
+    light_param = light_params[day_night]
+    light_data.energy = random.randint(*light_param['energy_bounds'])
+    light_data.color = light_param['sun_color']
+    back_val = random.uniform(*light_param['back_color_bounds'])
+    bpy.data.worlds["World"].node_tree.nodes["Background.001"].inputs[0].default_value = \
+        (back_val, back_val, back_val, 1)
+
+    ## Random background
+    back_img = random.choice(back_imgs[day_night])
+    # C.scene.node_tree.nodes["Image"].image = D.images[back_img]
+    image_holder.active_material.node_tree.nodes['Image Texture'].image = D.images[back_img]
+
+    # Camera movement
+    camera_object.rotation_euler = utils.r([72 + random.randint(-2, 2), 8, -82])
+
+    C.scene.render.filepath = fp
+    C.scene.render.image_settings.file_format = 'PNG'  # set output format to .png
+    C.scene.camera = camera_object
+
+    P, K, RT = camera_proj.get_3x4_P_matrix_from_blender(camera_object)
+
+    file_root_name = f'{list(scenes_ids).index(human_path)}_{scenes_ids[human_path]}'
+
+    with open(os.path.join(fp, 'cameras.json'), 'a+') as f_cam:
+        previous_cameras = f_cam.read()
+        if previous_cameras:
+            previous_cameras = json.loads(previous_cameras)
+        else:
+            previous_cameras = {}
+        previous_cameras[file_root_name] = {
+            'P': utils.mat_to_list(P),
+            'K': utils.mat_to_list(K),
+            'RT': utils.mat_to_list(RT),
+        }
+        json.dump(previous_cameras, f_cam, indent=4)
+
+    man.animation_data_clear()
+
+    ## Pour chaque paire, on veut NbTotalFrames / NbPersonnes / NbCar poses
+    ## Exemple: 150k / 200 / 2 = 1500 pose
+    with open(os.path.join(fp_ann_2D, f'annotations_{file_root_name}.csv'), 'w') as annot_file_2D, \
+            open(os.path.join(fp_ann_3D, f'annotations_{file_root_name}.csv'), 'w') as annot_file_3D:
+        bone_lbls = list(man.pose.bones.keys())
+        annot_file_2D.write(';'.join([lbl for bone in bone_lbls for lbl in [bone + k for k in ['_x', '_y']]]) + '\n')
+        annot_file_3D.write(
+            ';'.join([lbl for bone in bone_lbls for lbl in [bone + k for k in ['_X', '_Y', '_Z']]]) + '\n')
+
+        for po in range(nb_pose):
+            C.scene.frame_set(po * frame_rate)
+            ### On random une pose à chaque fois
+            random_pose.random_pose_ik(man)
+
+            ### On crée une frame
+            bpy.ops.object.mode_set(mode='OBJECT')
+
+            man.keyframe_insert(data_path="location", index=-1)
+            man.keyframe_insert(data_path="rotation_euler", index=-1)
+
+            bpy.ops.object.mode_set(mode='POSE')
+
+            for bone in man.pose.bones:
+                bone.keyframe_insert(data_path="rotation_euler", index=-1)
+                if bone.name[-3:] == '_IK':
+                    bone.keyframe_insert(data_path="location", index=-1)
+
+            bpy.ops.object.mode_set(mode='OBJECT')
+
+            ### Random caméra et lumières  (Blender proc)
+            ### On génère les images au bon  FPS
+            ### On sauvegarde les poses 3D + paramètres caméra (extra et intra)
+
+            # set output path so render won't get overwritten
+            C.scene.render.filepath = os.path.join(fp_img, f"{file_root_name}_{po}")
+            bpy.ops.render.render(write_still=True)  # render still
+
+            annotations_2D = []
+            annotations_3D = []
+            for lbl in bone_lbls:
+                bone_3d = utils.get_head_pose(lbl, man)
+                annotations_3D.append(f"{bone_3d[0]:.3f};{bone_3d[1]:.3f};{bone_3d[2]:.3f}")
+
+                bone_2d = P @ bone_3d
+                bone_2d /= bone_2d[-1]
+                annotations_2D.append(f"{bone_2d[0]:.2f};{bone_2d[1]:.2f}")
+            annot_file_2D.write(';'.join(annotations_2D) + '\n')
+            annot_file_3D.write(';'.join(annotations_3D) + '\n')
+
+C.scene.frame_end = int(frame_rate * (nb_pose - 0.5))
+utils.select_only(man)
+bpy.ops.object.mode_set(mode='POSE')
+
+with open(os.path.join(fp, 'infos.json'), 'w') as f:
+    json.dump({
+        'models': list(scenes_ids),
+        'id_max_scenes': scenes_ids
+    }, f, indent=4)
+
+print('Done', '#' * 25)
diff --git a/scripts/__init__.py b/scripts/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/scripts/camera_proj.py b/scripts/camera_proj.py
new file mode 100644
index 0000000000000000000000000000000000000000..76329f8b826f7a225f48e66ac2decb8eed5b7c44
--- /dev/null
+++ b/scripts/camera_proj.py
@@ -0,0 +1,118 @@
+import bpy
+import bpy_extras
+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]
+    location, rotation = cam.matrix_basis.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
+    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
+
+
+# ----------------------------------------------------------
+# Alternate 3D coordinates to 2D pixel coordinate projection code
+# adapted from https://blender.stackexchange.com/questions/882/how-to-find-image-coordinates-of-the-rendered-vertex?lq=1
+# to have the y axes pointing up and origin at the top-left corner
+def project_by_object_utils(cam, point):
+    scene = bpy.context.scene
+    co_2d = bpy_extras.object_utils.world_to_camera_view(scene, cam, point)
+    render_scale = scene.render.resolution_percentage / 100
+    render_size = (
+        int(scene.render.resolution_x * render_scale),
+        int(scene.render.resolution_y * render_scale),
+    )
+    return Vector((co_2d.x * render_size[0], render_size[1] - co_2d.y * render_size[1]))
diff --git a/scripts/human.py b/scripts/human.py
new file mode 100644
index 0000000000000000000000000000000000000000..718df66c591d70f569b1c3fd06d0f3a44226e8b2
--- /dev/null
+++ b/scripts/human.py
@@ -0,0 +1,297 @@
+import sys
+from math import cos, sin, pi, acos
+import random
+import json
+import glob
+import os
+
+import numpy as np
+from mathutils import Vector
+
+from scripts import utils
+
+import importlib
+
+importlib.reload(utils)
+from scripts.utils import *
+
+D = bpy.data
+
+
+# Load clothes list
+with open(r"mh_models\mh_mass_produce.json") as f:
+    data = json.load(f)
+    hairs = data['hair']['allNames']
+    tops = []
+    bottoms = []
+    suits = []
+    for name, cloth in data['clothes']['clothesInfos'].items():
+        if cloth["maleUpper"] or cloth["femaleUpper"] or cloth["mixedUpper"]:
+            tops.append(name.replace(' ', '_'))
+        if cloth["maleLower"] or cloth["femaleLower"] or cloth["mixedLower"]:
+            bottoms.append(name.replace(' ', '_'))
+        if cloth["maleFull"] or cloth["femaleFull"] or cloth["mixedFull"]:
+            bottoms.append(name.replace(' ', '_'))
+
+
+class Human:
+    def __init__(self, model):
+        self.model = model
+        self.name = self.model.name
+
+        self.init_textures()
+        self.init_bones()
+        self.top = self.init_color_ramp(tops + suits)
+        self.bot = self.init_color_ramp(bottoms + suits)
+
+    def init_color_ramp(self, clothes):
+        mat = None
+        for obj in self.model.children:
+            for clo in clothes:
+                if clo.lower() in obj.name.lower():
+                    if not len(obj.material_slots):
+                        break
+                    utils.select_only(obj)
+                    mat = obj.material_slots[0].material
+                    mat.use_nodes = True
+                    obj.active_material = mat
+                    nodes = mat.node_tree.nodes
+                    node_ramp = nodes.new("ShaderNodeValToRGB")
+                    node_ramp.color_ramp.elements.remove(node_ramp.color_ramp.elements[0])
+                    node_ramp.color_ramp.elements[0].position = 1.
+
+                    input_links = nodes['Principled BSDF'].inputs[0].links
+                    if len(input_links):
+                        img_node = nodes['Principled BSDF'].inputs[0].links[0].from_node
+
+                        node_mix = nodes.new("ShaderNodeMixRGB")
+                        node_mix.blend_type = 'MIX'
+                        node_mix.inputs[0].default_value = 0.5
+
+                        mat.node_tree.links.new(img_node.outputs[0], node_mix.inputs[1])
+                        mat.node_tree.links.new(node_ramp.outputs[0], node_mix.inputs[2])
+                        mat.node_tree.links.new(node_mix.outputs[0], nodes['Principled BSDF'].inputs[0])
+                    else:
+                        mat.node_tree.links.new(node_ramp.outputs[0], nodes['Principled BSDF'].inputs[0])
+
+                    node_ramp.color_ramp.elements[0].color = [random.random() for _ in range(3)] + [1]
+                    node_ramp.color_ramp.elements[0].color = [random.random() for _ in range(3)] + [1]
+                    break
+
+        return mat
+
+    def init_bones(self):
+        armature = bpy.data.armatures[self.model.name]
+        bpy.ops.object.mode_set(mode='EDIT')
+        for s in ('l', 'r'):
+            # Disconnect clavicle
+            armature.edit_bones[f'clavicle_{s}'].use_connect = False
+
+            # Create IK
+            hand = armature.edit_bones[f'hand_{s}']
+            hand_ik = armature.edit_bones.new(f'hand_{s}_IK')
+            utils.select_only_edit_bone(armature, hand_ik)
+            hand_ik.length = hand.length
+            bpy.ops.transform.resize(value=Vector((1.5, 1.5, 1.5)))
+
+        bpy.ops.object.mode_set(mode='POSE')
+        ik_constraints = {
+            'upperarm': [(-45, 120), (-90, 90), (-80, 80)],
+            'lowerarm': [(-30, 120), None, None]
+        }
+        for s in ('l', 'r'):
+            lowerarm = self.model.pose.bones[f'lowerarm_{s}']
+            lowerarm.constraints.new('IK')
+            lowerarm.constraints['IK'].target = self.model
+            lowerarm.constraints['IK'].subtarget = f'hand_{s}_IK'
+            lowerarm.constraints['IK'].chain_count = 2
+
+            for name, consts in ik_constraints.items():
+                bone = self.model.pose.bones[f'{name}_{s}']
+                for axe, const in zip(('x', 'y', 'z'), consts):
+                    if const is None:
+                        setattr(bone, f'lock_ik_{axe}', True)
+                    else:
+                        setattr(bone, f'lock_ik_{axe}', False)
+                        setattr(bone, f'use_ik_limit_{axe}', True)
+                        setattr(bone, f'ik_min_{axe}', utils.r(const[0]))
+                        setattr(bone, f'ik_max_{axe}', utils.r(const[1]))
+
+            # self.model.pose.bones[f'hand_{s}_IK'].location = Vector((-15, 10, 0))
+            self.model.pose.bones[f'hand_{s}_IK'].location = Vector((0, 5, -10))
+
+        bpy.ops.object.mode_set(mode='OBJECT')
+
+    def init_textures(self):
+        for o in self.model.children:
+            if "eye" in o.name or "high" in o.name:
+                continue
+            try:
+                for hair in hairs:
+                    if hair.lower().replace(' ', '_') in o.name.lower():
+                        o.active_material.blend_method = 'HASHED'
+                        break
+                else:
+                    o.active_material.blend_method = 'OPAQUE'
+
+            except AttributeError:
+                continue
+
+    def __call__(self, car=None):
+        if self.top is not None:
+            node = self.top.node_tree.nodes["ColorRamp"]
+            node.color_ramp.elements[0].color = [random.random() for _ in range(3)] + [1, ]
+        else:
+            print('Top ', self.model.name)
+
+        if self.bot is not None:
+            node = self.bot.node_tree.nodes["ColorRamp"]
+            node.color_ramp.elements[0].color = [random.random() for _ in range(3)] + [1, ]
+        else:
+            print('Bot ', self.model.name)
+
+        set_bounds(self.model, car)
+
+        return self.model
+
+
+def set_bounds(model, car=None):
+    # set_floors(model, car)
+    set_shrinkwraps(model, car)
+
+
+def set_floors(model, car=None):
+    original_mode = bpy.context.mode
+    utils.select_only(model)
+    bpy.ops.object.mode_set(mode='POSE')
+
+    if car is not None:
+        planes = [ch for ch in car.children_recursive if ch.name[:5] == 'Plane']
+
+    for s in 'lr':
+        bone = model.pose.bones[f'hand_{s}_IK']
+        [bone.constraints.remove(constr) for constr in bone.constraints if 'Floor' in constr.name]
+        if car is None:
+            continue
+        for obj in planes:
+            constr = bone.constraints.new('FLOOR')
+            constr.target = obj
+            constr.floor_location = 'FLOOR_Z'
+            constr.use_rotation = True
+            constr.offset = 1
+
+    bpy.ops.object.mode_set(mode=original_mode)
+
+
+def set_shrinkwraps(model, car=None):
+    original_mode = bpy.context.mode
+    utils.select_only(model)
+    bpy.ops.object.mode_set(mode='POSE')
+
+    if car is not None:
+        out_objects = [ch for ch in car.children_recursive if ch.name[:4] == 'OUT_']
+
+    for s in 'lr':
+        bone = model.pose.bones[f'hand_{s}_IK']
+        [bone.constraints.remove(constr) for constr in bone.constraints if 'Shrinkwrap' in constr.name]
+        if car is None:
+            continue
+        for obj in out_objects:
+            constr = bone.constraints.new('SHRINKWRAP')
+            constr.target = obj
+            constr.wrap_mode = 'OUTSIDE'
+            if 'Back' in obj.name:
+                constr.distance = model.pose.bones['lowerarm_l'].length * model.scale.x / obj.scale.y * 1
+            if 'Side' in obj.name:
+                constr.distance = model.pose.bones[
+                                      'lowerarm_l'].length * model.scale.x / obj.scale.z * 0.2
+            else:
+                constr.distance = model.pose.bones['hand_l'].length * model.scale.x / obj.scale.z * 1.5
+
+    bpy.ops.object.mode_set(mode=original_mode)
+
+
+class HumanLoader:
+    def __init__(self, dir_path, max_len=10):
+        self.human_paths = glob.glob(os.path.join(dir_path, '*.fbx'))
+        # random.shuffle(self.human_paths)
+        self.paths = {}
+        self.humans = {}
+        self.max_len = max_len
+        self.start_idx = 0
+        self.current_idx = -1
+
+        self.collection = D.collections.new("Humans")
+        C.scene.collection.children.link(self.collection)
+
+        self.picker = None
+        # self.load_next()
+
+    def load_next(self):
+        if self.max_len >= len(self.human_paths):
+            # If asked too much
+            self.start_idx = 0
+            if len(self.paths) == len(self.human_paths):
+                # If  everything already loaded
+                return
+            else:
+                next_idx = len(self.human_paths)
+        else:
+            next_idx = self.start_idx + self.max_len
+
+        human_paths = (self.human_paths * 2)[self.start_idx:next_idx]
+        human_paths = list({k: None for k in human_paths})
+
+        already_loaded = [hp for hp in human_paths if hp in self.paths.values()]
+        self.paths = {k: v for k, v in self.paths.items() if v in already_loaded}
+        self.humans = {k: v for k, v in self.humans.items() if k in already_loaded}
+        clear_humans(exceptions=list(self.paths))
+        for human_path in human_paths:
+            if human_path in already_loaded:
+                continue
+            bpy.ops.import_scene.fbx(filepath=human_path)
+
+            model = C.active_object
+            self.paths[model] = human_path
+            self.move_human(model)
+            self.humans[human_path] = Human(model)
+
+        self.picker = utils.Randomizer(list(self.humans.values()))
+        self.start_idx = next_idx % len(self.human_paths)
+
+    def next(self, *args, **kwargs):
+        self.current_idx += 1
+        if not self.current_idx % self.max_len:
+            self.start_idx = self.current_idx
+            self.load_next()
+            self.current_idx = 0
+        return self.picker.get(self.current_idx, *args, **kwargs)
+
+    def move_human(self, obj):
+        for ch in obj.children:
+            self.move_human(ch)
+
+        if obj.name in C.collection.objects:
+            self.collection.objects.link(obj)
+            C.collection.objects.unlink(obj)
+
+    def __call__(self, *args, **kwargs):
+        return self.picker(*args, **kwargs)
+
+
+def clear_humans(exceptions=[]):
+    collection = D.collections["Humans"]
+    exceptions = exceptions + [obj.children_recursive for obj in exceptions]
+    for obj in collection.objects:
+        if obj in exceptions:
+            continue
+        for ch in obj.children_recursive + [obj]:
+            obj_data = ch.data
+            if isinstance(obj_data, bpy.types.Armature):
+                D.armatures.remove(obj_data)
+            elif isinstance(obj_data, bpy.types.Mesh):
+                D.meshes.remove(obj_data)
+            try:
+                D.objects.remove(obj_data)
+            except ReferenceError:
+                pass
diff --git a/scripts/random_pose.py b/scripts/random_pose.py
new file mode 100644
index 0000000000000000000000000000000000000000..227f336c64493dc016852aafcf0ea595d11f92b6
--- /dev/null
+++ b/scripts/random_pose.py
@@ -0,0 +1,324 @@
+import sys
+from math import cos, sin, pi, acos
+import random
+
+from mathutils import Vector
+import numpy as np
+
+from scripts import utils
+
+import importlib
+
+importlib.reload(utils)
+from scripts.utils import *
+
+
+def default_rots():
+    return {
+        'thigh_r': (r(-70), r(-5 - random.random() * 15), 0),
+        'thigh_l': (r(-70), r(5 + random.random() * 15), 0),
+        'calf_r': (r(20 + random.random() * 50), r(30), 0),
+        'calf_l': (r(20 + random.random() * 50), r(-30), 0)
+    }
+
+
+def bounds():
+    values = {
+        'spine_03': [[0, 40], 45, None],
+        'upperarm_l': [None, 10, 10],
+        'upperarm_r': [None, 25, 20],
+        'neck_01': [[0, 40], 0, 10],
+        'lowerarm_X': [[-25, 50], [0, 35], None]
+    }
+
+    temp = values.copy()
+    for k, v in temp.items():
+        if '_X' in k:
+            values.update({k.replace('X', u): v for u in ['l', 'r']})
+            del (values[k])
+
+    for k, v in values.items():
+        values[k] = [[-a, a] if isinstance(a, int) else a for a in v]
+
+    return values
+
+
+def get_angles(angs):
+    new_angs = []
+    for ang in angs:
+        if isinstance(ang, list):
+            new_angs.append(random.randint(*ang))
+        elif isinstance(ang, int):
+            new_angs.append(random.randint(-ang, ang))
+        elif ang is None:
+            new_angs.append(0)
+
+    return Vector(r(new_angs))
+
+
+def reset_subject_old(subject):
+    bpy.ops.object.mode_set(mode='OBJECT')
+    subject.scale = (0.01, 0.01, 0.01)
+    subject.location = [-0.15, 0.5, 0.08]
+    subject.rotation_euler = r([-18, 0, 0])
+
+    bpy.ops.object.mode_set(mode='POSE')
+    for bone in subject.pose.bones:
+        bone.location = (0, 0, 0)
+        bone.rotation_euler = (0, 0, 0)
+
+
+def reset_subject(subject):
+    bpy.ops.object.mode_set(mode='POSE')
+    for bone in subject.pose.bones:
+        bone.rotation_mode = 'XYZ'
+        if bone.name in ["Root", "pelvis"]:
+            continue
+        bone.location = (0, 0, 0)
+        bone.rotation_euler = (0, 0, 0)
+
+    bpy.ops.object.mode_set(mode='OBJECT')
+    scale = round(subject.scale[0] * (1.6 + 0.3 * random.random()) / (subject.dimensions.y), 3)
+    # subject.scale = [scale] * 3
+
+    subject.scale = [0.1] * 3
+    height = subject.dimensions.y
+    # subject.location = [0, -0.16, -0.5 * 1.66]
+    subject.location = [0, -0.16, -0.5 * height]
+    subject.rotation_euler = r([65, 0, 0])
+
+
+def hand_pose(pose):
+    for s in ['l', 'r']:
+        hand_ratio = random.uniform(0.1, 0.8)
+        for finger in ['thumb', 'index', 'middle', 'ring', 'pinky']:
+            angles = r([40, 40, 40]) if finger == 'thumb' else r([70, 90, 40])
+            solo_ratio = random.uniform(0.7, 1) * hand_ratio
+            for i in range(3):
+                pose.bones[f'{finger}_{i + 1:02}_{s}'].rotation_euler.x = angles[i] * solo_ratio
+
+
+def random_pose_ik(subject, auto_ik=False):
+    """
+        1- reset and fix legs
+        2- randomize back
+        2b- bend back when too much twisted
+        3- randomize neck (backward proportional to back bending)
+        4- move arms with IK
+        :param subject: subject Object
+        :return:
+        """
+    # 1
+    bpy.ops.object.mode_set(mode='OBJECT')
+
+    pose = subject.pose
+    select_only(subject)
+
+    def rota(bone):
+        return Vector(pose.bones[bone].rotation_euler)
+
+    reset_subject(subject)
+    arm_length = dist(get_head_pose('upperarm_l', subject), get_head_pose('hand_l', subject))
+
+    base_rots = default_rots()
+    matrix_world = C.active_object.matrix_world.copy()
+    bpy.ops.object.mode_set(mode='POSE')
+
+    bounds_vals = bounds()
+    for bone, angles in base_rots.items():
+        pose.bones[bone].rotation_euler = angles
+
+    # 2
+    pose.bones['spine_03'].rotation_euler = get_angles(bounds_vals['spine_03'])
+
+    # 2b Compensate for shoulder in seat by bending back to front
+    pose.bones['spine_01'].rotation_euler = r(
+        Vector((random.randint(0, 10) + max(d(abs(rota('spine_03').y)) - 20, 0) / 2, 0,
+                0))) + rota('spine_01')
+
+    # 3
+    pose.bones['neck_01'].rotation_euler = (
+            get_angles(bounds_vals['neck_01']) +
+            get_angles(
+                [[d((rota('spine_01').x + rota('spine_03').x) * -0.5), 0], None, None]) +
+            Vector((0, random.uniform(0, 0.5) * rota('spine_03').y, 0))
+    )
+
+    hand_pose(pose)
+
+    pose.use_auto_ik = auto_ik
+
+    targets = {
+        'l': Vector((0.3, -0.1, 0.4)),
+        'r': Vector((-0.4, -0.1, 0.2))
+    }
+
+    back_rota_fact = sin(rota('spine_03').y) / sin(r(30))
+    for s in ['l', 'r']:
+        # Disconnect clavicle
+        armature = bpy.data.armatures[subject.name]
+
+        if auto_ik:
+            pose_bone = pose.bones[f'hand_{s}']
+        else:
+            pose_bone = pose.bones[f'hand_{s}_IK']
+
+        bone = pose_bone.bone
+        select_only_bone(armature, bone)
+
+        target = Vector()
+        min_arm_factor = 0.2
+        if False:
+            # target.x = random.random() * (0.4 if s == 'l' else -0.7)
+            target.x = random.uniform(max(-0.5, back_rota_fact * 0.4 + 0.4 if s == 'l' else 0),
+                                      min(0.4, back_rota_fact * 0.4 + 0 if s == 'l' else -0.4))
+            target.z = random.uniform(0.05, 0.7)
+            if abs(target.x) < 0.2:
+                target.y = random.uniform(-0.04, 0.02)
+            else:
+                target.y = random.uniform(-0.24, 0.02)
+            # sloap
+            target.y -= random.random() * (target.z * 0.1)
+        elif False:
+            phi = random.uniform(r(-180 + 45), 0) - (r(45) if s == 'r' else 0) + rota('spine_03').y
+            costheta = random.uniform(-1, 1)
+            u = random.uniform(0.1, 1)
+
+            theta = acos(costheta)
+            # u = arm_length * (u ** 1 / 3)
+
+            target.x = u * sin(theta) * cos(phi)
+            target.y = u * sin(theta) * sin(phi)
+            target.z = u * cos(theta)
+
+            shoulder_pose = get_head_pose(f'upperarm_{s}', subject)
+            target *= arm_length
+            target += shoulder_pose
+
+            # target.x = max(min(target.x, 0.32), -0.55)
+            # target.y = max(min(target.y, 0.02), -0.04 if abs(target.x) < 0.2 else -0.24)
+            # target.z = max(min(target.z, 0.7), 0.15)
+
+            bounded_width = Vector(((0.32 if s == 'l' else 0.55) + 0.35,
+                                    2 * arm_length,
+                                    0.7 - 0.15))
+            target = (
+                    (target - shoulder_pose) *
+                    (bounded_width / (2 * arm_length)) +
+                    shoulder_pose + 0.5 * (bounded_width - Vector([2 * arm_length, ] * 3))
+            )
+
+            # y done afterward because of the wheel
+            bounded_width_y = 0.02 - (- 0.04 if abs(target.x) < 0.2 else -0.24)
+            target.y = (
+                    (target.y - shoulder_pose.y) *
+                    (bounded_width_y / (0.9 * arm_length)) +
+                    shoulder_pose.y + 0.5 * (bounded_width_y - 0.9 * arm_length)
+            )
+        elif False:
+            phi = random.uniform(
+                r(-180) + rota('spine_03').y if s == 'r' else r(-120) + rota('spine_03').y,
+                r(-60) + rota('spine_03').y if s == 'r' else r(-60) + rota('spine_03').y
+            )
+            costheta = random.uniform(max(-0.8, 0.2 - cos(rota('spine_03').x + rota('spine_01').x)),
+                                      min(0.8, -0.2 + cos(rota('spine_03').x - rota('spine_01').x)))
+
+            theta = acos(costheta)
+            bounds_u = []
+            shoulder_pose = get_head_pose(f'upperarm_{s}', subject)
+            if not sin(theta) * cos(phi):
+                min_u_x = 0
+                max_u_x = arm_length
+            else:
+                bounds_u_x = (
+                    (-0.55 - shoulder_pose.x) / (sin(theta) * cos(phi)),
+                    (0.32 - shoulder_pose.x) / (sin(theta) * cos(phi))
+                )
+                min_u_x = (-0.55 - shoulder_pose.x) / (sin(theta) * cos(phi))
+                max_u_x = (0.32 - shoulder_pose.x) / (sin(theta) * cos(phi))
+                bounds_u += [min_u_x, max_u_x]
+
+            if not cos(theta):
+                min_u_z = 0
+                max_u_z = arm_length
+            else:
+                bounds_u_z = (
+                    (0.15 - shoulder_pose.z) / cos(theta),
+                    (0.7 - shoulder_pose.z) / cos(theta)
+                )
+                min_u_z = (0.15 - shoulder_pose.z) / cos(theta)
+                max_u_z = (0.7 - shoulder_pose.z) / cos(theta)
+
+                bounds_u += [min_u_z, max_u_z]
+
+            max_u_xz = min([m for m in bounds_u + [arm_length] if m > min_arm_factor * arm_length])
+            max_x = max_u_xz * sin(theta) * cos(phi) + shoulder_pose.x
+
+            if not sin(theta) * sin(phi):
+                min_u_y = 0
+                max_u_y = arm_length
+            else:
+                bounds_u_y = (
+                    ((-0.04 if abs(max_x) < 0.2 else -0.24) - shoulder_pose.y) / (sin(theta) * sin(phi)),
+                    (0.02 - shoulder_pose.y) / (sin(theta) * sin(phi))
+                )
+                min_u_y = ((-0.04 if abs(max_x) < 0.2 else -0.24) - shoulder_pose.y) / (sin(theta) * sin(phi))
+                max_u_y = (0.02 - shoulder_pose.y) / (sin(theta) * sin(phi))
+
+                bounds_u += [min_u_y, max_u_y]
+
+            max_u = min([m for m in bounds_u + [arm_length] if m > min_arm_factor * arm_length])
+
+            # print(d(phi), d(theta))
+            # print(min_u_x, min_u_y, min_u_z)
+            # print(max_u_x, max_u_y, max_u_z, max_u / arm_length)
+
+            # u = random.uniform((min_arm_factor) ** 1/2, (max_u / arm_length) ** 1/2) **  2 * arm_length
+            u = random.uniform(min_arm_factor * arm_length, max_u)
+            # print(u / arm_length)
+
+            target.x = u * sin(theta) * cos(phi)
+            target.y = u * sin(theta) * sin(phi)
+            target.z = u * cos(theta)
+
+            target += shoulder_pose
+
+        else:
+            back_forward_angle = rota('spine_03').x + rota('spine_01').x - r(30)  # 0 = straight
+
+            phi = random.uniform(
+                max((r(-160) if s == 'r' else r(-100)) + rota('spine_03').y, r(-160)),
+                min((r(-80) if s == 'r' else r(-40)) + rota('spine_03').y, r(-20))
+            )
+            theta_bound = 0.8
+            costheta = random.uniform(max(-0.8, -cos(theta_bound - back_forward_angle - rota('neck_01').x)),
+                                      min(0.8, cos(theta_bound + back_forward_angle + rota('neck_01').x)))
+
+            theta = acos(costheta)
+            bounds_u = []
+            shoulder_pose = get_head_pose(f'upperarm_{s}', subject)
+
+            # u = random.uniform((min_arm_factor) ** 1/2, (max_u / arm_length) ** 1/2) **  2 * arm_length
+            min_arm_factor = 0.2 + max(sin(back_forward_angle), 0)
+            # print(min_arm_factor, d(phi), d(theta))
+            u = random.uniform(min_arm_factor, 1) * arm_length
+
+            target.x = u * sin(theta) * cos(phi)
+            target.y = u * sin(theta) * sin(phi)
+            target.z = u * cos(theta)
+
+            target += shoulder_pose
+
+        temp_rota = rota(f'upperarm_{s}') + rota(f'lowerarm_{s}')
+        # target = targets[s]
+
+        location = get_head_pose(bone.name, subject)
+        # print(location)
+        # print(target)
+        bpy.ops.transform.translate(value=target - location)
+
+    bpy.ops.object.mode_set(mode='OBJECT')
+
+
+if __name__ == '__main__':
+    subject = get_object('Subject')
diff --git a/scripts/utils.py b/scripts/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..4320460998b9cfb575c4d60df5f7fc09e2da6fc2
--- /dev/null
+++ b/scripts/utils.py
@@ -0,0 +1,147 @@
+import random
+from math import pi
+
+import bpy
+import numpy as np
+from mathutils import Matrix, Vector, Euler
+
+
+# Maths utils
+def rd(a, d):
+    if type(a) in [float, int]:
+        if d:
+            return int(180 * a / pi)
+        else:
+            return pi * a / 180.
+    else:
+        try:
+            iterator = iter(a)
+        except TypeError as te:
+            raise ValueError('Cant convert to radians ', a)
+        return type(a)([r(k) for k in a])
+
+
+def r(a):
+    return rd(a, d=False)
+
+
+def d(a):
+    return rd(a, d=True)
+
+
+def dist(a, b):
+    return sum([k ** 2 for k in (a - b)]) ** 0.5
+
+
+# Blender utils
+C = bpy.context
+D = bpy.data
+
+
+def hide_object(obj, hide=True):
+    for o in obj.children:
+        hide_object(o, hide=hide)
+    obj.hide_set(hide)
+    obj.hide_select = hide
+    if 'Plane' in obj.name or 'OUT_' in obj.name:
+        obj.hide_render = True
+        obj.hide_set(True)
+    else:
+        obj.hide_render = hide
+
+
+def get_object(name):
+    return bpy.data.objects[name]
+
+
+def unselect_all():
+    select_only(None)
+
+
+def select_only(obj=None):
+    for ob in bpy.data.objects:
+        ob.select_set(False)
+
+    if obj is None:
+        return
+
+    obj.select_set(True)
+    bpy.context.view_layer.objects.active = obj
+
+
+def select_only_bone(armature, bone=None):
+    for bo in armature.bones:
+        bo.select = False
+
+    if bone is None:
+        return
+
+    bone.select = True
+
+
+def select_only_edit_bone(armature, bone=None):
+    for bo in armature.edit_bones:
+        bo.select = False
+        bo.select_head = False
+        bo.select_tail = False
+
+    if bone is None:
+        return
+
+    bone.select = True
+    bone.select_head = True
+    bone.select_tail = True
+    armature.edit_bones.active = bone
+
+
+def get_head_pose(bone, struct='Subject'):
+    if isinstance(struct, str):
+        struct = get_object(struct)
+
+    current_mode = bpy.context.mode
+    bpy.ops.object.mode_set(mode='OBJECT')
+    matrix_world = struct.matrix_world
+    bpy.ops.object.mode_set(mode=current_mode)
+    return (matrix_world @ struct.pose.bones[bone].head.to_4d()).to_3d()
+
+
+# Other utils
+def mat_to_list(mat):
+    return np.array(mat).tolist()
+
+
+class Randomizer:
+    def __init__(self, objects):
+        self.objects = objects
+        for obj in self:
+            hide_object(obj)
+
+    def to_list(self):
+        return [obj if isinstance(obj, bpy.types.Object) else obj.model for obj in self.objects]
+
+    def __iter__(self):
+        return iter(self.to_list())
+
+    def __getitem__(self, item):
+        return self.to_list()[item]
+
+    def __len__(self):
+        return len(self.objects)
+
+    def get(self, pick_idx, *args, **kwargs):
+        pick = self.objects[pick_idx]
+        pick = pick if isinstance(pick, bpy.types.Object) else pick(*args, **kwargs)
+        self.swap_object(self[pick_idx])
+        return pick
+
+    def __call__(self, *args, **kwargs):
+        pick_idx = random.randint(0, len(self) - 1)
+        return self.get(pick_idx)
+
+    def swap_object(self, obj=None):
+        for o in self:
+            if not o.hide_get():
+                hide_object(o)
+
+        if obj is not None:
+            hide_object(obj, hide=False)