RAMAC Example on Custom Volumetric images

RAMAC repository is tailored for implementation of the algorithm on Monaleesa Computed Tomography (CT) data set as outlined in the below paper : Image registration based automated lesion correspondence pipeline for longitudinal CT data

Hence users can use this tool for their own volumetric image dataset. The following function provides a convenient format for users to establish correspondence within their own volumetric image datasets. By simply providing the directory path containing the images and the corresponding Region of Interest (ROI) coordinate information, users can effortlessly establish the necessary correspondence.

import sys
sys.path.append("../ramac")
from registration import *
from correspondence_csv_input import *
from transform_coordinates import *
from merge_dataframe import *
from phantominator import shepp_logan
from plots import *
from input_transform import *
from preprocessing import mask_air
from utils import *
#import warnings
import scipy
import numpy as np
import SimpleITK as sitk
import pandas as pd
import os

The below function is to perform a test using custom data. This function shows the process to load DICOM series for longitudinal image data at two timepoints, performs rigid 3D registration,perform adaaptive hungarian and saves result files, and plots triplanar views of the images.

def test_custom_data():

    ## Directory 1 = '/projects01/didsr-novartis/reference/Novartis_Monaleesa/Images/ML3_deid/1003003/Screening/CT - Chest - Abdomen - Pelvis - Venous - Spiral'
    ## Directory 2 = '/projects01/didsr-novartis/reference/Novartis_Monaleesa/Images/ML3_deid/1003003/Week 8/CT - Chest - Abdomen - Pelvis - Venous - Spiral'

    Screening_image = load_dicom_series('Directory 1') ## User needs to give the directory path
    Week_8_image = load_dicom_series('Directory 2')

    ## Preprocess the images if needed
    fixed_image = mask_air(Screening_image)
    moving_image = mask_air(Week_8_image)
    # # Apply the 3D rigid registration
    moving_resampled, [init_transform, final_transform] = registration_3d_rigid_series(fixed_image, moving_image)

    result_folder = "monaleesa"
    os.makedirs(result_folder, exist_ok=True)  # Create the result folder if it doesn't exist

    # File paths for saving CSV files which contains the ROI coordinates as inputs in CSV format (to be placed by user)

    filename_fixed_R1 = os.path.join(result_folder, 'df_lesions_Screening_1.csv')
    filename_fixed_R2 = os.path.join(result_folder, 'df_lesions_Screening_2.csv')
    filename_moving_R1 = os.path.join(result_folder, 'df_lesions_Week_8_1.csv')
    filename_moving_R2 = os.path.join(result_folder, 'df_lesions_Week_8_2.csv')

    ## Apply Adaptive Hungarian
    threshold = 30
    correspondences_screening, unmatched_names_df1_screening, unmatched_names_df2_screening = process_lesion_timepoints(filename_fixed_R1, filename_fixed_R2, threshold)
    plot_lesion_correspondences_timepoints(filename_fixed_R1, filename_fixed_R2)
    final_df_screening = create_final_dataframe_timepoints(correspondences_screening, unmatched_names_df1_screening, unmatched_names_df2_screening, filename_fixed_R1, filename_fixed_R2)

    correspondences_week_8, unmatched_names_df1_week_8, unmatched_names_df2_week_8 = process_lesion_timepoints(filename_moving_R1, filename_moving_R2, threshold)
    plot_lesion_correspondences_timepoints(filename_moving_R1, filename_moving_R2)
    final_df_week_8 = create_final_dataframe_timepoints(correspondences_week_8, unmatched_names_df1_week_8, unmatched_names_df2_week_8, filename_moving_R1, filename_moving_R2)

    final_df_S = create_trim_dataframe(final_df_screening)
    final_df_8 = create_trim_dataframe(final_df_week_8)

    ## Save the dataframes of the ROI coordiantes in CSV format
    filename_fixed = os.path.join(result_folder, 'fixed_coordinates.csv')
    filename_moving = os.path.join(result_folder, 'moving_coordinates.csv')
    filename_registered = os.path.join(result_folder, 'registered_coordinates.csv')
    filename_correspondence = os.path.join(result_folder, 'correspondence_indices.csv')

    save_transformed_dataframe(final_df_S, filename_fixed)
    save_transformed_dataframe(final_df_8, filename_moving)
    ## Apply the transformation matrix on ROI cooridnates of moving image (Image B) to obtain registered ROI cooridnates
    final_transform_1 = final_transform.GetInverse()

    registered_coordinates = create_transformed_dataframe(filename_moving, final_transform_1)
    save_transformed_dataframe(registered_coordinates, filename_registered)

    ## Apply Adaptive Hungarian on the fixed ROI cooridnates and registered coordinates
    threshold = 30
    correspondences, unmatched_names_df1, unmatched_names_df2 = process_lesion_timepoints(filename_fixed, filename_registered, threshold)
    plot_lesion_correspondences_timepoints(filename_fixed, filename_registered)
    final_df = create_final_dataframe_timepoints(correspondences, unmatched_names_df1, unmatched_names_df2, filename_fixed, filename_registered)

    ## Saving the final correspondence thereby achieved in a CSV format
    df_transformed = merge_indices(final_df)
    df_transformed.to_csv(filename_correspondence, index=False)

    ## Transform the ROI coordinates from physical to voxel unit for plotting

    fixed_voxel = transform_physical_to_index(fixed_image, filename_fixed)
    fixed_voxel_tuples = voxel_tuples_from_dataframe(fixed_voxel)
    moving_voxel = transform_physical_to_index(moving_image, filename_moving)
    moving_voxel_tuples = voxel_tuples_from_dataframe(moving_voxel)
    registered_voxel = transform_physical_to_index(moving_resampled, filename_registered)
    registered_voxel_tuples = voxel_tuples_from_dataframe(registered_voxel)

    ## Triplanar visualization after registration

    window_level = 100
    window_width = 2000

    images = [fixed_image, moving_image, moving_resampled]
    voxel_tuples_list = [fixed_voxel_tuples, moving_voxel_tuples, registered_voxel_tuples]
    image_names = ['Fixed', 'Moving', 'Registered']

    for image, voxel_tuples, image_name in zip(images, voxel_tuples_list, image_names):
        plot_triplanar(image, voxel_tuples, window_level, window_width, image_name)


test_custom_data()