intensity correction for Siemens .dat files (of reconstructed epi brain images)

[mfaizan-ai]

I am working with siemens raw data, and I read the data using twixtools package, there are multiple .dat files and I read each with the twixtools. So twixtools first measurement (measurement index 0 ) is the prescan with head coil and body coil each with the (32, 32, 38, 128) and (32, 32, 2, 128) respectively. And my reconstructed data is in this format (Time, Slices, Phase, Coils, Readout), for instance (10, 36, 64, 38, 128), how can I customise your code to this? I tried but it didn't work for me. Any recommendations?

[LajerQAQ]

From the prescan data from the first measurement:

Head coil: (32, 32, 38, 128)

Body coil: (32, 32, 2, 128)

These are sufficient for creating the low-resolution correction maps.

From the geometry information in the second measurement, we can determine where to interpolate the plane for the correction map and at what resolution.

Did you run the demo code with the demo data? If so, did it work?

When you ran the code with your data, what error did you see? Can you paste the error message?

This code should work for your data. I assume you just need to add an index (e.g., an extra index for the time dimension) to make it run.

Once you are able to obtain the correction maps, applying them to your data may only require a transpose to align the dimensions.

My suggestions are:

First, check the dimension information of the demo data by printing:

dim_info_data, dim_info_ref

after reading the data with the rawdata_reader() function.

Second, use this information as guidance to modify the code.

Please let me know whether the code works or if you have any further questions.

[mfaizan-ai]

@LajerQAQ Thank you so much for you response. :)

I tried to run your noteobook that is provided in the repo for intensity inhomoengity correction. So I had to customised the code to load my twix data as it does not load the prescan data which include the head kspace and the body kspace. so I write a custom function to extract it into a numpy arrays. Therefore, each prescan array has the following shapes:

kspace_body_shape: (32, 32, 2, 128)
kspace_head_shape: (32, 32, 38, 128)

And the I load the EPI.

And load single band EPI data which has the shape : (T, 36, 64, 38, 128), the last dimension is oversampled non removed. After removing the oversampling teh shape becomes the 64 for readout. Then I run your code and the I got the images that I will attach here. this is the code is I running after some customisation in the helpers functions and in the notebooks.

auto_rotation = 'Dicom' # 'Dicom' or 'LGE', for using LGE rotation you need to have "SAX, LAX or 2CH" in the file names
apply_correction_to_sensitivity_maps = False # Select mode for reconstruction
channel_keep = 20 # Number of channels to keep after coil compression if you don't want to do coil compression then set this to None
lamb = 5e-2      # Regularization parameter for generating the correction maps, the larger the value the smoother the maps
oversampling_phase_factor = 1 # Set this to 1 if you don't want to do oversampling phase correction other wise set it to 3. Other values larger than 3 may or may not work
remove_readout_oversampling = True # For removing the readout oversampling for low resolution data, set this to True can speed up the code, but if the object is large (example: part of the object is outside of the image), you may want to set this to False
par_index_method = 1 # par_index_method: Method to determine the method for fetching quaternions and offsets for different slices, try set it to 0 if the output results do not match your data. (default is 1).

for i, data_path_filename in enumerate(path_input):
    
    twix, image_3D_body_coils, image_3D_surface_coils, data, dim_info_data, data_ref, dim_info_ref, num_sli = rawdata_reader(data_path_filename)

    #coil compression with PCA, you can use the u parameter to preform the same coil compression on extra data
    try:
        data, data_ref, image_3D_surface_coils, u = compress_data_with_pca(data = data, channel_keep = channel_keep, axis_CH= dim_info_data.index('Cha'), 
                                extra_data=data_ref, extra_axis_CH=dim_info_ref.index('Cha'), 
                                extra_data2=image_3D_surface_coils, extra_axis_CH2=-1)
    except:
        #when there is no reference data
        data, image_3D_surface_coils, u = compress_data_with_pca(data = data, channel_keep = channel_keep, axis_CH= dim_info_data.index('Cha'),
                            extra_data=image_3D_surface_coils, extra_axis_CH=-2)

    #coil compression end
    img_correction_map, sens_correction_map, low_resolution_surface_coil_imgs, img_quats = correction_map_generator(twix, 
                                                                                                                    image_3D_body_coils, image_3D_surface_coils, 
                                                                                                                    data, dim_info_data, num_sli , 
                                                                                                                    auto_rotation = auto_rotation, 
                                                                                                                    lamb = lamb, tol = 1e-4, maxiter=500,
                                                                                                                    apply_correction_to_sensitivity_maps = apply_correction_to_sensitivity_maps,
                                                                                                                    oversampling_phase_factor = oversampling_phase_factor,
                                                                                                                    remove_readout_oversampling = remove_readout_oversampling,
                                                                                                                    par_index_method = par_index_method)

    # Initialize SENSE reconstruction results array
    try:
        sense_recon_results = np.zeros((num_sli,data.shape[dim_info_data.index('Lin')],data.shape[dim_info_data.index('Col')]//2,data.shape[dim_info_data.index('Par')]),dtype=np.complex64)
    except:
        sense_recon_results = np.zeros((num_sli,data.shape[dim_info_data.index('Lin')],data.shape[dim_info_data.index('Col')]//2),dtype=np.complex64)
    
    # Only do recon one image for each different slices (this is for demo purpose)
    data,dim_info_data = data_reduction(data,data_dimensions = dim_info_data, dims_to_keep = ['Sli', 'Lin', 'Cha', 'Col', 'Par'])

    try:
        for n in range(num_sli):
            ############################ Custom Reconstruction Process Start ############################
            #preprocessing for the sense reconstruction
            ksp,ref_padded = pad_ref(data,data_ref,n,dim_info_ref = dim_info_ref,dim_info_org=dim_info_data)
            ksp = remove_RO_oversamling(ksp,axis_RO= [item for item in dim_info_data if item != 'Sli'].index('Col') ) # we are indexing the Sli dimension here, so we need to input the position of the Col dimension without the Sli dimension
            ref_padded = remove_RO_oversamling(ref_padded,axis_RO= [item for item in dim_info_ref if item != 'Sli'].index('Col') )

            #rotate the ksp and ref_padded
            ksp = sense_input_rotation(auto_rotation, img_quats[n], ksp)
            ref_padded = sense_input_rotation(auto_rotation, img_quats[n], ref_padded)


            if apply_correction_to_sensitivity_maps:
                #if we want to apply the correction during the sense reconstruction
                sense_reconstructed_img = sense_reconstruction(ksp,ref_padded,sens_correction_map[n,...])
            else:
                sense_reconstructed_img = sense_reconstruction(ksp,ref_padded)

            try:
                sense_recon_results[n,...] = sense_reconstructed_img
            except:
                sense_recon_results = np.zeros((num_sli,sense_reconstructed_img.shape[0],sense_reconstructed_img.shape[1]),dtype=np.complex64)
                sense_recon_results[n,...] = sense_reconstructed_img
            ############################ Custom Reconstruction Process End ############################
    except:
        print("\033[91m" + "If you want to apply the correction to the sensitivity maps, please make sure you have modified Custom Reconstruction part for your datasets.")
        print("Using low resolution surface coil images instead.\033[0m")

        try:
            for n in range(num_sli):
                sense_recon_results[n,...] = low_resolution_surface_coil_imgs
        except:
            sense_recon_results = low_resolution_surface_coil_imgs
        
    #we are rotating the images here (optional, may not apply to other datasets)
    sense_recon_results, img_correction_map, sens_correction_map = auto_image_rotation(sense_recon_results, img_quats, 
                                                                            auto_rotation = auto_rotation, 
                                                                            img_correction_map = img_correction_map, 
                                                                            sens_correction_map = sens_correction_map,
                                                                            filename = data_path_filename)
    
    #save the results
    save_sense_recon_results(path_output[i], sense_recon_results, img_correction_map, sens_correction_map, img_quats, apply_correction_to_sensitivity_maps)

and I am not sure what does that auto_rotation means here. and other parameters and not sure if I am using them correctly ore not. I will attach the images that I get form the running this code, I am not sure if its a correct maps of the images or not. these are images of sense reconstruction slices for 36 slices and the image correction maps for the 36 slices, do they make any sense?

your guidance on that? do I need to run the second cell that you provided? or is it sufficient to normalise?

[LajerQAQ]

"auto_rotation" is a function that automatically aligns cardiac image orientations. You can leave it at the default setting for now.

Did you obtain these SENSE reconstruction results by running the cell you pasted above? In the cell you shared, the "Custom Reconstruction Process" section can only process the demo code because SENSE reconstruction varies across datasets. For custom datasets, it will output low-resolution surface-coil images. The second figure looks like low-resolution surface-coil images—please correct me if I’m wrong.

A quick way to check whether these correction maps are correct, and whether you need to change their orientation, is to compare the low-resolution surface-coil images with the actual reconstruction of your data. The low-resolution surface-coil images should look very similar to your reconstruction, but at lower resolution.

Once you confirm the correction maps are what you need, you can simply perform pointwise multiplication between img_correction_map and your reconstructions.

One more thing that may help is to run the following 4D flow demo dataset. If I remember correctly, I didn’t customize the SENSE reconstruction part for this dataset. You should see low-resolution outputs (if you see high-resolution outputs, then I may have modified the SENSE code at some point).

print("Demo dataset(s) are being downloaded, please wait......")
#download_file_from_figshare(saving_path,url_2CH) # 461.95 MB May take a couple of minutes to download <- commet this line . You may want to move or clean the downloaded dataset in the data folder manually
#download_file_from_figshare(saving_path,url_SAX) # 3.63 GB May take tens of minutes to download
#download_file_from_figshare(saving_path,url_LAX) # 457.44 MB May take a couple of minutes to download
download_file_from_figshare(saving_path,url_4DFlow) # 636.57 MB May take a couple of minutes to download <- uncommet this line in the demo ipynb

[mfaizan-ai]

@LajerQAQ these files are not downloadable for me?

[LajerQAQ]

You should be able to download all datasets listed in the demo ipynb file. Please let me know if you encounter any error.