diff --git a/README.md b/README.md
index 74ccf97..4447770 100644
--- a/README.md
+++ b/README.md
@@ -631,8 +631,23 @@ You should check the quality of the reconstructed block at this stage, and make
 
 ### Inter-block registration
 
-The registration strategy is described in the data paper, but the details of registration can differ for different specimens. As a result, no registration scripts are distributed at the moment.
+The inter-block registration process is described in detail in the data descriptor. It comprises several steps that can be applied in a different order if necessary. Thus, the code associated with this step consists of a “toolbox” rather than a complete pipeline.
 
+The tools are integrated into the `interblock_registration.py` code and can be called using the following command:
+```shell
+# Overview of the tools available
+phcp-interblock-tools -h
+
+# Selction of the tools
+phcp-interblock-tools name_tool -h
+```
+The following tools can be called:
+1. `extract_hull`: Create the shell of a segmentation by applying a combination of dilation and erosion and subtracting the results. This tool is useful to help for extracting the raw cut surface.
+2. `run_svd`: Returns the normal vector and centroid of the estimated mean plane calculated on the raw cut surface. This allows the plane equation to be reconstructed in order to project the raw surface.
+3. `project`: Project voxels from the raw surface orthogonally in the given equation plane. **inputs**: NifTI file of the manual segmentation of the cut surface, output filename, centroïd coordinates, normal vector.  **output**: NifTI file of the raw cut surface projected on the plane.
+4. `extract_intensity': Create a 2D NifTI image of the voxels within the estimated plane. **inputs**: NifTI file of the reference 3D space, output filename, centroïd coordinates, normal vector, in which direction to project the voxels (--usexz, --usexy, --useyz).
+5. `backproj`: Create a 3D NifTI image from a 2D image by projecting it within the given equation plane. **inputs**: NifTI file of the reference 3D space, NifTI file of the 2D image, output filename, centroïd coordinates, normal vector, in which direction to project he voxels (--usexz, --usexy, --useyz). **output**: 3D NifTI file of the 2D image projected in the plane according to the given equation.
+6. `TPS_interp`: Create a deformation field from two points cloud of landmarks by applying TPS interpolation. **inputs**: NifTI image of the landmarks in the reference image (7T), NifTI image of the landmarks in the reconstructed block (11.7T), NifTI image of the space of reference, output filename. **output**: Deformation field (NifTI).
 
 ### Concatenation of transformations (optional)
 
diff --git a/phcp/scripts/interblock_registration.py b/phcp/scripts/interblock_registration.py
new file mode 100755
index 0000000..32c0c93
--- /dev/null
+++ b/phcp/scripts/interblock_registration.py
@@ -0,0 +1,490 @@
+#!/usr/bin/env python3
+
+import logging
+import sys
+from collections.abc import Sequence
+from pathlib import Path
+
+import ants
+import nibabel as ni
+import numpy as np
+import scipy.ndimage as nd
+from scipy.ndimage import center_of_mass, map_coordinates
+
+logger = logging.getLogger(__name__)
+
+
+def spherical_structure(radius):
+    size = 2 * radius + 1
+    center = radius
+    X, Y, Z = np.ogrid[:size, :size, :size]
+    sphere = (X - center) ** 2 + (Y - center) ** 2 + (Z - center) ** 2
+    return sphere.astype(np.uint8)
+
+
+def extract_hull(
+    segmentation_filename: str | Path, output_filename: str | Path
+) -> None:
+    meta = ni.load(segmentation_filename)
+    arr = meta.get_fdata()
+    structure = spherical_structure(3)
+    erode = nd.binary_erosion(arr, structure=structure, iterations=2)
+    dilate = nd.binary_dilation(arr, structure=structure, iterations=2)
+    newarr = dilate ^ erode
+    ni_im = ni.Nifti1Image(newarr, meta.affine, meta.header)
+    ni.save(ni_im, output_filename)
+
+
+def plane_extractor_SVD(filename: str | Path) -> None:
+    meta = ni.load(filename)
+    arr = meta.get_fdata()
+    points = np.argwhere(arr > 0)
+
+    logger.info("Downsampled part")
+    sampling_ratio = 0.005
+    indices = np.arange(len(points))
+    np.random.shuffle(indices)
+    sample_size = int(len(points) * sampling_ratio)
+    sampled_indices = indices[:sample_size]
+    downsampled_points = points[sampled_indices]
+
+    centroid = np.mean(downsampled_points, axis=0)
+    centered_points = downsampled_points - centroid
+
+    logger.info("SVD part")
+    _, _, vh = np.linalg.svd(centered_points)
+    normal_vector = vh[-1, :]
+
+    logger.info("Centroid =" + str(centroid))
+    logger.info("Normal vector =" + str(normal_vector))
+
+
+def get_2dplane(
+    block_filename: str | Path,
+    output_filename: str | Path,
+    centroid: Sequence[float],
+    normal_vector: Sequence[float],
+    usexz: bool,
+    usexy: bool,
+    useyz: bool,
+) -> None:
+    meta = ni.load(block_filename)
+    arr = meta.get_fdata()
+
+    centroid = np.asarray(centroid, dtype=float)
+    normal_vector = np.asarray(normal_vector, dtype=float)
+
+    d = -centroid.dot(normal_vector)
+
+    x_dim, y_dim, z_dim = arr.shape
+
+    if usexz:
+        x = np.linspace(0, x_dim - 1, x_dim)
+        z = np.linspace(0, z_dim - 1, z_dim)
+
+        xx, zz = np.meshgrid(x, z)
+
+        yy = (-normal_vector[0] * xx - normal_vector[2] * zz - d) / normal_vector[1]
+
+    if usexy:
+        x = np.linspace(0, x_dim - 1, x_dim)
+        y = np.linspace(0, y_dim - 1, y_dim)
+
+        xx, yy = np.meshgrid(x, y)
+
+        zz = (-normal_vector[0] * xx - normal_vector[1] * yy - d) / normal_vector[2]
+
+    if useyz:
+        z = np.linspace(0, z_dim - 1, z_dim)
+        y = np.linspace(0, y_dim - 1, y_dim)
+
+        yy, zz = np.meshgrid(y, z)
+
+        xx = (-normal_vector[1] - normal_vector[2] * zz * yy - d) / normal_vector[0]
+
+    # Stack the coordinates for interpolation
+    coords = np.stack([xx, yy, zz], axis=-1)
+
+    # Interpolate the values at the grid points
+    plane_2d_image = map_coordinates(
+        arr, [coords[..., 0], coords[..., 1], coords[..., 2]], order=3, mode="nearest"
+    )
+
+    affine = np.eye(4) * meta.header["pixdim"][1]
+    affine[3][3] = 1
+
+    ni.save(ni.Nifti1Image(plane_2d_image, affine), output_filename)
+
+
+def back_proj(
+    block_filename: str | Path,
+    slice_to_project_filename: str | Path,
+    output_filename: str | Path,
+    centroid: Sequence[float],
+    normal_vector: Sequence[float],
+    usexz: bool,
+    usexy: bool,
+    useyz: bool,
+) -> None:
+    logger.info("Load data")
+    meta = ni.load(block_filename)
+    arr = meta.get_fdata()
+
+    meta_seg = ni.load(slice_to_project_filename)
+    modified_plane_2d_image = meta_seg.get_fdata()
+
+    centroid = np.asarray(centroid, dtype=float)
+    normal_vector = np.asarray(normal_vector, dtype=float)
+
+    d = -centroid.dot(normal_vector)
+
+    # Define the bounding box for the grid points
+    x_dim, y_dim, z_dim = arr.shape
+
+    # Map the modified 2D plane data back to 3D space
+    modified_3d_plane = np.zeros_like(arr)
+
+    logger.info("Fill 3D space")
+    if usexz:
+        x = np.linspace(0, x_dim - 1, x_dim)
+        z = np.linspace(0, z_dim - 1, z_dim)
+
+        xx, zz = np.meshgrid(x, z)
+
+        yy = (-normal_vector[0] * xx - normal_vector[2] * zz - d) / normal_vector[1]
+
+        # Iterate over the grid and update the 3D image
+        for i in range(z_dim - 1):
+            for j in range(x_dim - 1):
+                x, y, z = int(xx[i, j]), int(yy[i, j]), int(zz[i, j])
+                if 0 <= x < x_dim and 0 <= y < y_dim and 0 <= z < z_dim:
+                    modified_3d_plane[x, y, z] = modified_plane_2d_image[i, j]
+    if usexy:
+        x = np.linspace(0, x_dim - 1, x_dim)
+        y = np.linspace(0, y_dim - 1, y_dim)
+
+        xx, yy = np.meshgrid(x, y)
+
+        zz = (-normal_vector[0] * xx - normal_vector[1] * yy - d) / normal_vector[2]
+
+        for i in range(y_dim - 1):
+            for j in range(x_dim - 1):
+                x, y, z = int(xx[i, j]), int(yy[i, j]), int(zz[i, j])
+                if 0 <= x < x_dim and 0 <= y < y_dim and 0 <= z < z_dim:
+                    modified_3d_plane[x, y, z] = modified_plane_2d_image[i, j]
+
+    if useyz:
+        z = np.linspace(0, z_dim - 1, z_dim)
+        y = np.linspace(0, y_dim - 1, y_dim)
+
+        yy, zz = np.meshgrid(y, z)
+
+        xx = (-normal_vector[1] - normal_vector[2] * zz * yy - d) / normal_vector[0]
+
+        for i in range(z_dim - 1):
+            for j in range(y_dim - 1):
+                x, y, z = int(xx[i, j]), int(yy[i, j]), int(zz[i, j])
+                if 0 <= x < x_dim and 0 <= y < y_dim and 0 <= z < z_dim:
+                    modified_3d_plane[x, y, z] = modified_plane_2d_image[i, j]
+    logger.info("Save 3D file")
+    ni.save(
+        ni.Nifti1Image(modified_3d_plane, meta.affine, meta.header), output_filename
+    )
+
+
+def distance_projection(normal_vector, centroid, xa, ya, za):
+    d = -centroid.dot(normal_vector)
+    a, b, c = normal_vector[0], normal_vector[1], normal_vector[2]
+    return np.abs(a * xa + b * ya + c * za + d) / (np.linalg.norm(normal_vector))
+
+
+def lambda_proj(normal_vector, centroid, xa, ya, za):
+    d = -centroid.dot(normal_vector)
+    a, b, c = normal_vector[0], normal_vector[1], normal_vector[2]
+    return (a * xa + b * ya + c * za + d) / (np.linalg.norm(normal_vector) ** 2)
+
+
+def projection(
+    cutting_surface_filename: str | Path,
+    output_filename: str | Path,
+    centroid: Sequence[float],
+    normal_vector: Sequence[float],
+) -> None:
+    meta = ni.load(cutting_surface_filename)
+    arr = meta.get_fdata()
+
+    centroid = np.asarray(centroid, dtype=float)
+    normal_vector = np.asarray(normal_vector, dtype=float)
+    size = arr.shape
+    res = np.zeros(size)
+    a, b, c = normal_vector[0], normal_vector[1], normal_vector[2]
+    for xa in range(size[0]):
+        for ya in range(size[1]):
+            for za in range(size[2]):
+                if (
+                    arr[xa, ya, za] == 1
+                    and distance_projection(normal_vector, centroid, xa, ya, za) < 30
+                ):  # np.abs(value)<5 and
+                    const = lambda_proj(normal_vector, centroid, xa, ya, za)
+                    xh, yh, zh = (
+                        int(xa - const * a),
+                        int(ya - const * b),
+                        int(za - const * c),
+                    )
+
+                    if xh > 0 and yh > 0 and zh > 0:
+                        res[xh][yh][zh] = 1
+
+    ni.save(ni.Nifti1Image(res, meta.affine, meta.header), output_filename)
+
+
+def extract_centroids(seg_img_path):
+    meta = ni.load(seg_img_path)
+    arr = meta.get_fdata()
+    labels = np.unique(arr)[1:]
+    logger.info(labels)
+    centroid_voxel = [center_of_mass(arr == el) for el in labels]
+    return np.array(centroid_voxel), meta
+
+
+def run_tps_interpolation(
+    ref_7T_filename: str | Path,
+    hr_filename: str | Path,
+    reconstructed_block_filename: str | Path,
+    output_filename: str | Path,
+) -> None:
+    fixed_ants = ants.image_read(ref_7T_filename)
+    moving_ants = ants.image_read(hr_filename)
+
+    logger.info("Extraction of the centroid for landmarks in fixed space")
+    fixed_centroid, meta_fixed = extract_centroids(ref_7T_filename)
+
+    liste_centroid = [list(el) for el in fixed_centroid]
+    liste_centroid_phys_ants = np.asarray(
+        [
+            ants.transform_index_to_physical_point(fixed_ants, el)
+            for el in np.int16(liste_centroid)
+        ]
+    )
+
+    landmark_7T_npy_filename = (
+        Path(ref_7T_filename).parent / "Segmentation7T_Landmarks.npy"
+    )
+    np.save(landmark_7T_npy_filename, liste_centroid_phys_ants)
+
+    logger.info("Extraction of the centroid for landmarks in moving space")
+    moving_centroid, meta_moving = extract_centroids(hr_filename)
+
+    liste_centroid = [list(el) for el in moving_centroid]
+    liste_centroid_phys_ants = np.asarray(
+        [
+            ants.transform_index_to_physical_point(moving_ants, el)
+            for el in np.int16(liste_centroid)
+        ]
+    )
+
+    landmark_HR_npy_filename = Path(hr_filename).parent / "SegmentationHR_Landmarks.npy"
+    np.save(landmark_HR_npy_filename, liste_centroid_phys_ants)
+
+    logger.info("TPS interpolation")
+    moving_points_phySpace = np.load(landmark_HR_npy_filename)
+    fixed_points_phySpace = np.load(landmark_7T_npy_filename)
+    tx = ants.landmark_transforms.fit_transform_to_paired_points(
+        moving_points=moving_points_phySpace,
+        fixed_points=fixed_points_phySpace,
+        domain_image=fixed_ants,
+        transform_type="tps",
+        verbose=True,
+    )
+
+    logger.info("Convert and save the warp field in ANTS fomat")
+    tx_field = ants.transform_to_displacement_field(tx, fixed_ants)
+    tx_field.to_file(Path(hr_filename).parent / "Landmark68_TPS0Warp.nii.gz")
+
+    logger.info("Apply the deformation field on the reconstructed block given")
+    img = ants.image_read(reconstructed_block_filename)
+    img2 = tx.apply_to_image(img)
+    ants.image_write(img2, output_filename)
+
+
+def parse_command_line(argv):
+    """Parse the script's command line."""
+    import argparse
+
+    parser = argparse.ArgumentParser(
+        description="Inter-bloc registration toolbox",
+    )
+    subparsers = parser.add_subparsers(dest="command", required=True)
+
+    # sub parser skull extraction
+    p1 = subparsers.add_parser(
+        "extract_hull", help="Extract hull of a binary segmentation"
+    )
+    p1.add_argument(
+        "binary_segmentation_filename", help="Filename of the binary segmentation"
+    )
+    p1.add_argument("output_filename", help="Ouput filename of the hull")
+    p1.set_defaults(
+        func=lambda args: extract_hull(
+            args.binary_segmentation_filename, args.output_filename
+        )
+    )
+
+    # sub parser plane_extractor_SVD
+    p2 = subparsers.add_parser(
+        "run_svd",
+        help="Extraction of the centroid and the normal vector to the plane from the SVD",
+    )
+    p2.add_argument("filename", help="Name of the file on which to perform the SVD")
+    p2.set_defaults(func=lambda args: plane_extractor_SVD(args.filename))
+
+    # sub parser intensity extraction
+    p3 = subparsers.add_parser(
+        "extract_intensity", help="Extract intensity within the equation plan"
+    )
+    p3.add_argument("input_filename", help="Filename of the reconstructed block")
+    p3.add_argument("output_filename", help="Filename of the 2D slice extracted")
+    p3.add_argument(
+        "--centroid",
+        nargs=3,
+        type=float,
+        required=True,
+        help="Centroid coordinates (x y z)",
+    )
+    p3.add_argument(
+        "--normal",
+        nargs=3,
+        type=float,
+        required=True,
+        help='Normal vector (nx ny nz). WARNING: if you have negative values, you must enclose the value in quotation marks and leave a space at the beginning (e.g. " -2.8456e-09")',
+    )
+    p3.add_argument(
+        "--usexz", action="store_true", help="Use the algorithm to extract xz plane"
+    )
+    p3.add_argument(
+        "--usexy", action="store_true", help="Use the algorithm to extract xy plane"
+    )
+    p3.add_argument(
+        "--useyz", action="store_true", help="Use the algorithm to extract yz plane"
+    )
+    p3.set_defaults(
+        func=lambda args: get_2dplane(
+            args.input_filename,
+            args.output_filename,
+            args.centroid,
+            args.normal,
+            args.usexz,
+            args.usexy,
+            args.useyz,
+        )
+    )
+
+    # sub parser back projection
+    p4 = subparsers.add_parser(
+        "backproj", help="Project back the intensity in the 3D space"
+    )
+    p4.add_argument("input_filename", help="Filename of the reconstructed block")
+    p4.add_argument("seg_filename", help="Filename of the manual segmentation")
+    p4.add_argument("output_filename", help="Filename of the 2D slice extracted")
+    p4.add_argument(
+        "--centroid",
+        nargs=3,
+        type=float,
+        required=True,
+        help="Centroid coordinates (x y z)",
+    )
+    p4.add_argument(
+        "--normal",
+        nargs=3,
+        type=float,
+        required=True,
+        help='Normal vector (nx ny nz). WARNING: if you have negative values, you must enclose the value in quotation marks and leave a space at the beginning (e.g. " -2.8456e-09")',
+    )
+    p4.add_argument(
+        "--usexz", action="store_true", help="Use the algorithm to extract xz plane"
+    )
+    p4.add_argument(
+        "--usexy", action="store_true", help="Use the algorithm to extract xy plane"
+    )
+    p4.add_argument(
+        "--useyz", action="store_true", help="Use the algorithm to extract yz plane"
+    )
+    p4.set_defaults(
+        func=lambda args: back_proj(
+            args.input_filename,
+            args.seg_filename,
+            args.output_filename,
+            args.centroid,
+            args.normal,
+            args.usexz,
+            args.usexy,
+            args.useyz,
+        )
+    )
+
+    # sub parser intensity extraction
+    p5 = subparsers.add_parser(
+        "project",
+        help="Project in the space of filename the binary cutting surface segmentation in the estimated plane",
+    )
+    p5.add_argument(
+        "input_filename",
+        help="Filename of the manual segmentation of the cutting surface",
+    )
+    p5.add_argument(
+        "output_filename",
+        help="Filename of the cutting surface projected in the estimated plane",
+    )
+    p5.add_argument(
+        "--centroid",
+        nargs=3,
+        type=float,
+        required=True,
+        help="Centroid coordinates (x y z)",
+    )
+    p5.add_argument(
+        "--normal",
+        nargs=3,
+        type=float,
+        required=True,
+        help='Normal vector (nx ny nz). WARNING: if you have negative values, you must enclose the value in quotation marks and leave a space at the beginning (e.g. " -2.8456e-09")',
+    )
+    p5.set_defaults(
+        func=lambda args: projection(
+            args.input_filename, args.output_filename, args.centroid, args.normal
+        )
+    )
+
+    # sub parser run TPS interpolation
+    p6 = subparsers.add_parser(
+        "TPS_interp",
+        help="run TPS interpolation and apply on the reconstructed block given",
+    )
+    p6.add_argument(
+        "landmark_ref",
+        help="Filename of the file containing the landmark in the fixed space (7T)",
+    )
+    p6.add_argument(
+        "landmark_HR",
+        help="Filename of the file containing the landmark in the moving space (T2w)",
+    )
+    p6.add_argument("HR_block", help="Filename of the reconstructed block to correct")
+    p6.add_argument("output_filename", help="Filename of HR_block corrected")
+    p6.set_defaults(
+        func=lambda args: run_tps_interpolation(
+            args.landmark_ref, args.landmark_HR, args.HR_block, args.output_filename
+        )
+    )
+    args = parser.parse_args()
+    return args
+
+
+def main(argv=sys.argv):
+    """The script's entry point."""
+    logging.basicConfig(level=logging.INFO)
+    args = parse_command_line(argv)
+    return args.func(args) or 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/phcp/scripts/interfov_fusion.py b/phcp/scripts/interfov_fusion.py
index 637489e..bc5c35d 100755
--- a/phcp/scripts/interfov_fusion.py
+++ b/phcp/scripts/interfov_fusion.py
@@ -41,6 +41,48 @@ def insert_transformation_file(
             )
 
 
+def add_padding_to_fov(
+    input_filename: str | os.PathLike, padding_1: int, padding_0: int
+) -> None:
+    FOV_filename = input_filename.replace("_NLMF_N4_rec-unwarp.nii.gz", "_FOV.nii.gz")
+
+    meta = ni.load(input_filename)
+    ni.save(ni.Nifti1Image(np.ones(meta.shape), meta.affine, meta.header), FOV_filename)
+
+    FOV_weight_filename = FOV_filename.replace(".nii.gz", "_weight.nii.gz")
+
+    fov_mapping = np.ones(meta.shape)
+
+    pad_fov_mapping = np.pad(
+        fov_mapping,
+        pad_width=((padding_1, padding_1), (padding_1, padding_1), (0, 0)),
+        mode="constant",
+        constant_values=1,
+    )
+    pad_fov_mapping = np.pad(
+        pad_fov_mapping,
+        pad_width=((padding_0, padding_0), (padding_0, padding_0), (5, 5)),
+        mode="constant",
+        constant_values=0,
+    )
+
+    distance_mapping = nd.distance_transform_edt(
+        pad_fov_mapping, meta.header["pixdim"][1]
+    )
+    weight_mapping = sigmoid(distance_mapping, 2, 4)
+
+    # Change the header to integrate padding
+    affine = meta.affine
+    affine[0][3] -= meta.header["pixdim"][1] * (padding_1 + padding_0)
+    affine[1][3] -= meta.header["pixdim"][3] * 5
+    affine[2][3] -= meta.header["pixdim"][2] * (padding_1 + padding_0)
+
+    ni.save(
+        ni.Nifti1Image(weight_mapping, affine, meta.header),
+        FOV_weight_filename,
+    )
+
+
 def merge_FOVs(File_directory: str | os.PathLike, JSON_filename: str) -> None:
     preparation_directory = os.path.join(File_directory, "02-Preparation")
     transformfile_directory = os.path.join(File_directory, "03-TransformFiles")
@@ -85,18 +127,36 @@ def merge_FOVs(File_directory: str | os.PathLike, JSON_filename: str) -> None:
 
         logger.info(f"Send {filename} to the intermediate space : Done")
 
-        logger.info("========= Creation of FOV mapping =========")
-        meta = ni.load(input_filename)
+        logger.info(
+            "========= Creation of the geometric penalty in the initial space ========="
+        )
+
         FOV_filename = input_filename.replace(
             "_NLMF_N4_rec-unwarp.nii.gz", "_FOV.nii.gz"
         )
+
+        FOV_weight_filename = FOV_filename.replace(".nii.gz", "_weight.nii.gz")
+
         output_FOV_filename = os.path.join(
             distribution_directory,
             FOV_filename.split("/")[-1].replace(".nii.gz", "_IS.nii.gz"),
         )
-        ni.save(
-            ni.Nifti1Image(np.ones(meta.shape), meta.affine, meta.header), FOV_filename
+
+        output_FOV_weight_filename = os.path.join(
+            distribution_directory,
+            FOV_weight_filename.split("/")[-1].replace(".nii.gz", "_IS.nii.gz"),
         )
+
+        add_padding_to_fov(input_filename, padding_1=30, padding_0=5)
+
+        run_ants_apply_registration(
+            intermediate_space_filename,
+            FOV_weight_filename,
+            output_FOV_weight_filename,
+            transforms_list,
+            interpolation="NearestNeighbor",
+        )
+
         run_ants_apply_registration(
             intermediate_space_filename,
             FOV_filename,
@@ -104,10 +164,6 @@ def merge_FOVs(File_directory: str | os.PathLike, JSON_filename: str) -> None:
             transforms_list,
             interpolation="NearestNeighbor",
         )
-        logger.info("Creation of FOV mapping : Done")
-        logger.info("========= Creation of FOV weight =========")
-        calcul_weight_FOVs(output_FOV_filename)
-        logger.info("Creation of FOV weight : Done")
 
         if not (i > 0):
             reconstructed_block_nifti_im = ni.load(output_filename)
@@ -156,7 +212,7 @@ def merge_consecutive_FOVs(
     )
     act_fov_mask_filename = "".join([actual_subject_filename, "_FOV_IS.nii.gz"])
     act_fov_weight_filename = "".join(
-        [actual_subject_filename, "_FOV_IS_weight.nii.gz"]
+        [actual_subject_filename, "_FOV_weight_IS.nii.gz"]
     )
 
     meta_mask_prec = ni.load(prec_fov_mask_filename)
@@ -206,17 +262,6 @@ def find_coefficient_to_correct_intensity(
     return gmm2.means_.flatten().max() / gmm.means_.flatten().max()
 
 
-def calcul_weight_FOVs(fov_mapping_filename: str | os.PathLike) -> None:
-    meta = ni.load(fov_mapping_filename)
-    arr = meta.get_fdata()
-    distance_mapping = nd.distance_transform_edt(arr, meta.header["pixdim"][1])
-    weight_mapping = sigmoid(distance_mapping, 2, 4)
-    ni.save(
-        ni.Nifti1Image(weight_mapping, meta.affine, meta.header),
-        fov_mapping_filename.replace(".nii.gz", "_weight.nii.gz"),
-    )
-
-
 def parse_command_line(argv):
     """Parse the script's command line."""
     import argparse
diff --git a/pyproject.toml b/pyproject.toml
index e742887..7085b18 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -28,6 +28,7 @@ phcp-diffusion-pipeline = "phcp.fov_scripts.diffusion_pipeline:main"
 phcp-intrafov-registration = "phcp.scripts.intrafov_registration:main"
 phcp-interfov-registration = "phcp.scripts.interfov_registration:main"
 phcp-interfov-fusion = "phcp.scripts.interfov_fusion:main"
+phcp-interblock-tools = "phcp.scripts.interblock_registration:main"
 phcp-concat-transforms = "phcp.scripts.concat_transforms:main"
 phcp-transform-and-fuse = "phcp.scripts.transform_and_fuse:main"
 phcp-prepare-published-dataset = "phcp.scripts.prepare_published_dataset:main"