diff --git a/.gitignore b/.gitignore index 541c02e..aacf428 100644 --- a/.gitignore +++ b/.gitignore @@ -28,6 +28,9 @@ test/integration/*.xprot test/integration/*.xml .atom-build.json #.vscode/ +test/nhlbi_integration_tests/test +test/nhlbi_integration_tests/baselines +*.seq dep-build/ stats.csv junit.xml diff --git a/Dockerfile b/Dockerfile index 749d977..46df286 100644 --- a/Dockerfile +++ b/Dockerfile @@ -3,7 +3,7 @@ ARG USERNAME="vscode" ARG USER_UID=1000 ARG USER_GID=1000 -FROM ubuntu:22.04 AS gadgetron_baseimage +FROM ubuntu:24.04 AS gadgetron_baseimage LABEL org.opencontainers.image.source=https://github.com/gadgetron/gadgetron ARG USERNAME @@ -12,9 +12,13 @@ ARG USER_GID ARG HOME=/home/$USERNAME RUN apt-get update \ - && DEBIAN_FRONTEND=noninteractive apt-get install -y sudo wget git-core git-lfs rsync curl net-tools libxml2 \ + && DEBIAN_FRONTEND=noninteractive apt-get install -y sudo wget git-core git-lfs rsync curl net-tools libxml2 emacs-nox \ && apt-get clean +# Remove default ubuntu user (Ubuntu 24.04 ships with UID/GID 1000 taken) +RUN userdel -r ubuntu 2>/dev/null || true \ + && groupdel ubuntu 2>/dev/null || true + # Create the user RUN groupadd --gid $USER_GID $USERNAME \ && useradd --uid $USER_UID --gid $USER_GID -m $USERNAME -s /bin/bash \ @@ -41,20 +45,7 @@ RUN wget --no-hsts --quiet https://github.com/conda-forge/miniforge/releases/dow && chmod -R g+w /opt/conda \ && find /opt -type d | xargs -n 1 chmod g+s -#RUN wget --no-hsts --quiet https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh -O /tmp/miniforge.sh \ -# && /bin/bash /tmp/miniforge.sh -b -p /opt/conda \ -# && rm /tmp/miniforge.sh \ -# && /opt/conda/bin/mamba clean --tarballs --index-cache --packages --yes \ -# && find /opt/conda -follow -type f -name '*.a' -delete \ -# && find /opt/conda -follow -type f -name '*.pyc' -delete \ -# && /opt/conda/bin/mamba clean --force-pkgs-dirs --all --yes \ -# && groupadd -r conda --gid ${CONDA_GID} \ -# && usermod -aG conda ${USERNAME} \ -# && chown -R :conda /opt/conda \ -# && chmod -R g+w /opt/conda \ -# && find /opt -type d | xargs -n 1 chmod g+s - -# Copy environment, which will be filtered for later staged +# Copy environment COPY --chown=$USER_UID:conda environment.yml /tmp/build/ # Create mount points for tests @@ -76,8 +67,8 @@ ADD https://github.com/krallin/tini/releases/download/${TINI_VERSION}/tini /tini RUN chmod +x /tini USER root -RUN mkdir -p /opt/GIRF/ COPY --chown=$USER_UID:conda toolboxes/nhlbi_gt_toolbox/GIRF/ /opt/GIRF/ +#COPY --chown=$USER_UID:conda toolboxes/nhlbi_gt_toolbox/models_dl/ /opt/models/ FROM gadgetron_baseimage AS gadgetron_dev_cuda ARG USER_UID @@ -85,8 +76,12 @@ ARG HOME USER ${USER_UID} RUN mkdir -p ${HOME}/.cache/conda/notices && sudo chown -R ${USER_UID}:conda ${HOME}/.cache/conda/notices -RUN grep -v "#.*\" /tmp/build/environment.yml > /tmp/build/filtered_environment.yml -RUN umask 0002 && /opt/conda/bin/mamba env remove -n gadgetron || true && /opt/conda/bin/mamba env create -f /tmp/build/filtered_environment.yml && /opt/conda/bin/mamba clean -afy && sudo chown -R :conda /opt/conda +RUN umask 0002 && /opt/conda/bin/mamba env remove -n gadgetron || true \ + && /opt/conda/bin/mamba env create -f /tmp/build/environment.yml \ + && /opt/conda/bin/mamba clean -afy \ + && sudo chown -R :conda /opt/conda + + RUN echo 'export PATH="/opt/conda/envs/gadgetron/nvvm/bin/:$PATH"' >> ${HOME}/.bashrc @@ -98,8 +93,7 @@ ARG USER_UID ARG HOME USER ${USER_UID} RUN mkdir -p ${HOME}/.cache/conda/notices && sudo chown -R ${USER_UID}:conda ${HOME}/.cache/conda/notices -RUN grep -v "#.*\" /tmp/build/environment.yml > /tmp/build/filtered_environment.yml -RUN umask 0002 && /opt/conda/bin/mamba env create -f /tmp/build/filtered_environment.yml && /opt/conda/bin/mamba clean -afy && sudo chown -R :conda /opt/conda +RUN umask 0002 && /opt/conda/bin/mamba env create -f /tmp/build/environment.yml && /opt/conda/bin/mamba clean -afy && sudo chown -R :conda /opt/conda USER root FROM gadgetron_dev_cuda AS gadgetron_cudabuild @@ -108,28 +102,16 @@ USER ${USER_UID} WORKDIR /opt ARG USER_GID=1000 -RUN echo "UID=$USER_UID GID=$USER_GID" - -RUN sudo chown $USER_UID:$USER_GID /opt && mkdir -p /opt/code/gadgetron && mkdir -p /opt/package +RUN sudo chown $USER_UID:$USER_GID /opt && mkdir -p /opt/code/gadgetron && mkdir -p /opt/package COPY --chown=$USER_UID:conda . /opt/code/gadgetron/ - -# Ensure NHLBI-GT-Non-Cartesian is inside the build context before building -#COPY --chown=$USER_UID:conda NHLBI-GT-Non-Cartesian/ /opt/code/NHLBI-GT-Non-Cartesian/ - - ENV CONDA_ENV_PATH=/opt/conda/envs/gadgetron RUN echo "export NUMBA_CUDA_ENABLE_PYNVJITLINK=1" >> $CONDA_ENV_PATH/etc/conda/activate.d/env_vars.sh -#COPY --chown=$USER_UID:conda OpticalFlow3d/ /opt/code/OpticalFlow3d/ RUN . /opt/conda/etc/profile.d/conda.sh && umask 0002 && conda activate gadgetron && \ sh -x && \ pip install git+https://github.com/ahsanjav/OpticalFlow3d.git - # cd /opt/code/OpticalFlow3d/ && \ - # pip install -e . - # # && \ - # pip install numpy==1.23 - + SHELL ["/bin/bash", "-c"] RUN . /opt/conda/etc/profile.d/conda.sh && umask 0002 && conda activate gadgetron && \ export PATH="/opt/conda/envs/gadgetron/nvvm/bin/:$PATH" && \ @@ -165,28 +147,25 @@ USER ${USER_UID} ARG USER_GID=1000 RUN mkdir -p ${HOME}/.cache/conda/notices && sudo chown -R ${USER_UID}:conda ${HOME}/.cache/conda/notices -RUN grep -v "#.*\" /tmp/build/environment.yml > /tmp/build/filtered_environment.yml -RUN umask 0002 && /opt/conda/bin/mamba env create -f /tmp/build/filtered_environment.yml && /opt/conda/bin/mamba clean -afy && sudo chown -R :conda /opt/conda +# Copy the already-solved conda env from the build stage (eliminates redundant ~20 min conda solve) +COPY --from=gadgetron_cudabuild --chown=$USER_UID:conda /opt/conda/envs/gadgetron /opt/conda/envs/gadgetron + +# Copy build artifacts and entrypoint COPY --from=gadgetron_cudabuild --chown=$USER_UID:conda /opt/package /opt/conda/envs/gadgetron/ COPY --from=gadgetron_cudabuild --chown=$USER_UID:conda /opt/code/gadgetron/docker/entrypoint.sh /opt/ -#RUN mkdir -p /opt/models/ #COPY --from=gadgetron_cudabuild --chown=$USER_UID:conda /opt/models/ /opt/models/ ENV CONDA_ENV_PATH=/opt/conda/envs/gadgetron RUN echo "export NUMBA_CUDA_ENABLE_PYNVJITLINK=1" >> $CONDA_ENV_PATH/etc/conda/activate.d/env_vars.sh -#COPY --chown=$USER_UID:conda OpticalFlow3d/ /opt/code/OpticalFlow3d/ -RUN . /opt/conda/etc/profile.d/conda.sh && umask 0002 && conda activate gadgetron && \ - sh -x && \ - pip install git+https://github.com/ahsanjav/OpticalFlow3d.git - # cd /opt/code/OpticalFlow3d/ && \ - # pip install . -RUN sudo rm -rf /opt/code/ -RUN chmod +x /opt/entrypoint.sh -RUN sudo mkdir -p /opt/integration-test && sudo chown ${USER_GID}:${USER_UID} /opt/integration-test +RUN sudo rm -rf /opt/code/ \ + && chmod +x /opt/entrypoint.sh \ + && sudo mkdir -p /opt/integration-test \ + && sudo chown ${USER_GID}:${USER_UID} /opt/integration-test COPY --from=gadgetron_cudabuild --chown=$USER_UID:conda /opt/code/gadgetron/test/integration /opt/integration-test/ +COPY --from=gadgetron_cudabuild --chown=$USER_UID:conda /opt/code/gadgetron/test/nhlbi_integration_tests /opt/nhlbi-integration-test/ ENTRYPOINT [ "/tini", "--", "/opt/entrypoint.sh" ] # FROM gadgetron_baseimage AS gadgetron_rt_nocuda diff --git a/README.md b/README.md index 037237e..4179387 100644 --- a/README.md +++ b/README.md @@ -1,15 +1,12 @@ -# “Two-for-one”: 3D cardiac and pulmonary MR imaging from a single acquisition using bSTAR - -This project is a Gadgetron extension that provides a dual inline image reconstruction pipeline for simultaneous cardiac-resolved whole heart and respiratory-resolved lung imaging. +# Volumetric real-time MRI for the guidance of cardiac catheterization at 0.55 T +This project is a Gadgetron extension that provides a inline image reconstruction pipeline for volumetric, stack-of-spirals real-time acquisitions to guide invasive procedures such as right-heart catheterization. ## Publication -This open-source is described in : - -- Daudé P.,Taylor J., Plummer J., Ramasawmy R., Tsakirellis A., Moss J., Bandettini W.P., Chen M., Campbell-Washburn A.E.,Javed A. “Two-for-one”: 3D cardiac and pulmonary MR imaging from a single acquisition using bSTAR -. In preparation +This project is described in : +- Kumar P., Ramasawmy R., Javed A., Le D., Daudé P., O'Brien K., Jaimes A.E., Franson D., Chow K., Dong F., Maier F., Lederman R.J., Nayak K.S., Campbell-Washburn A.E., “Volumetric real-time MRI for the guidance of cardiac catheterization at 0.55 T". In preparation ## Documentation diff --git a/apps/gadgetron/connection/HeaderConnection.cpp b/apps/gadgetron/connection/HeaderConnection.cpp index 2dc1b89..673535a 100644 --- a/apps/gadgetron/connection/HeaderConnection.cpp +++ b/apps/gadgetron/connection/HeaderConnection.cpp @@ -48,11 +48,15 @@ namespace { /* tSequenceVariant in the header give information about the hash of the pulseq seq file Loop through all seq files in a directory and find the correct seq file (hash is in the last line ) */ + + // PK commenting this out ecause the presence of tSequenceVariant is not guaranteed for only bSTAR... + /* if (header.userParameters) { ISMRMRD::UserParameters user_params = header.userParameters.get(); std::vector strings = user_params.userParameterString; std::vector::iterator it; + for (it = strings.begin(); it != strings.end(); ++it) { if (it->name == "tSequenceVariant") { GDEBUG("Hash bstar found %s found\n", it->value.c_str()); @@ -159,9 +163,11 @@ namespace { } } } + */ //PD header_callback(header); } + private: std::function header_callback; diff --git a/environment.yml b/environment.yml index a8d009c..92ddac3 100644 --- a/environment.yml +++ b/environment.yml @@ -15,7 +15,8 @@ dependencies: - numpy=1.24.4 - scipy=1.10.1 - pyfftw=0.15.0 - - h5py>=3.6,<3.8 # Must match conda HDF5 1.10.x + - hdf5>=1.14.3 + - h5py=3.12.1 # Requires hdf5 >=1.14.3 (provided by ismrmrd=1.15.0) - pyyaml=6.0 - jq=1.6 - yq=2.13 @@ -52,27 +53,25 @@ dependencies: # Math and Signal Processing - fftw=3.3.9 - - armadillo=9.900.5 + - armadillo>=9.900 - mkl>=2024.2.2,<2025.0a0 - mkl-include>=2024.2.2,<2025.0a0 - - libxml2=2.9 - - libcurl=7.79.1 + - libxml2>=2.9.14,<2.14 + - libcurl>=7.79 - plplot=5.15.0 - - pugixml=1.12.1 - - onnxruntime=1.10.0 + - pugixml=1.15 + - onnxruntime>=1.10 - xsdata>=22.2 - # Gadgetron + ISMRMRD - - ismrmrd::ismrmrd=1.13.7 - - ismrmrd::ismrmrd-python>=1.9.8 - - ismrmrd::siemens_to_ismrmrd=1.2.11 + # Gadgetron + ISMRMRD (C libraries via conda; Python bindings via pip below) + - ismrmrd::ismrmrd=1.15.0 + - ismrmrd::siemens_to_ismrmrd=1.2.13 - ismrmrd::mrd-storage-server=0.0.12 - - gadgetron::gadgetron-python=1.4.1 - gadgetron::range-v3=0.11.0 # Boost - boost=1.80.0 - - icu=70.1 + - icu>=70.1 # Docs - sphinx=4.4.0 @@ -94,6 +93,10 @@ dependencies: - pynvjitlink - pip: + # ismrmrd and gadgetron Python bindings installed via pip + # (conda ismrmrd-python channel only has up to 1.14.2) + - ismrmrd==1.14.2 + - gadgetron==1.4.1 - nvidia-cuda-runtime-cu12==12.8.57 # Required for CuPy JIT compilation (must match PyTorch cu128 requirement) - sigpy - scikit-image @@ -108,3 +111,8 @@ dependencies: - torchvision==0.22.0 - torchaudio==2.7.0 - onnxruntime-gpu + # NHLBI integration tests (Azure Blob Storage for private test data) + - azure-storage-blob + - azure-identity + - matplotlib + - pyigtl diff --git a/gadgets/python/legacy/gadgets/accumulate_and_recon.py b/gadgets/python/legacy/gadgets/accumulate_and_recon.py index 7e5c552..68078c3 100644 --- a/gadgets/python/legacy/gadgets/accumulate_and_recon.py +++ b/gadgets/python/legacy/gadgets/accumulate_and_recon.py @@ -38,7 +38,7 @@ def process(self, acq, data, *args): if (acq.flags & (1<<7)): #Is this the last scan in slice image = cifftn(self.myBuffer, axes=(0, 1, 2)) - image = image * np.product(image.shape)*100 #Scaling for the scanner + image = image * np.prod(image.shape)*100 #Scaling for the scanner #Create a new image header and transfer value img_head = ismrmrd.ImageHeader() img_head.version = 1 diff --git a/test/nhlbi_integration_tests/cases/3drt.cfg b/test/nhlbi_integration_tests/cases/3drt.cfg new file mode 100644 index 0000000..506c0d4 --- /dev/null +++ b/test/nhlbi_integration_tests/cases/3drt.cfg @@ -0,0 +1,52 @@ +[dependency.siemens] +data_file = 3drt/noise_data.h5 +measurement = 0 +additional_arguments = skip_converstion + +[dependency.client] +configuration = default_measurement_dependencies.xml + +[dependency.siemens.1] +data_file = 3drt/prescan-reference_noise.h5 +measurement = 0 +additional_arguments = skip_converstion + +[dependency.client.1] +configuration = default_measurement_dependencies.xml + +[dependency.siemens.2] +data_file = 3drt/prescan-reference.h5 +measurement = 0 +additional_arguments = skip_converstion + +[dependency.client.2] +configuration = spiral_3DRT_CSM_test.xml + +[reconstruction.siemens] +data_file = 3drt/recon_data.h5 +measurement = 0 +additional_arguments = skip_converstion + +[reconstruction.client] +configuration = spiral_3DRT_test.xml + +[reconstruction.test] +reference_file = 3drt/baseline_output.h5 +reference_images = spiral_3DRT_test.xml/image_0 +output_images = spiral_3DRT_test.xml/image_0 +value_comparison_threshold = 0.1 +scale_comparison_threshold = 0.01 + +[requirements] +system_memory = 8192 +gpu_support = 1 +gpu_memory = 8192 + +[tags] +tags = nhlbi,3drt + +[nhlbi] +description = spiral 3drt +noise_file = 3drt/noise_data.h5 +baseline_recon_time = 42.3 + diff --git a/test/nhlbi_integration_tests/cases/csm_prescan.cfg b/test/nhlbi_integration_tests/cases/csm_prescan.cfg new file mode 100644 index 0000000..12d8613 --- /dev/null +++ b/test/nhlbi_integration_tests/cases/csm_prescan.cfg @@ -0,0 +1,36 @@ +[dependency.siemens] +data_file = csm_prescan/noise_data.h5 +measurement = 0 +additional_arguments = skip_converstion + +[dependency.client] +configuration = default_measurement_dependencies.xml + +[reconstruction.siemens] +data_file = csm_prescan/recon_data.h5 +measurement = 0 +additional_arguments = skip_converstion + +[reconstruction.client] +configuration = spiral_3DRT_CSM_test.xml + +[reconstruction.test] +reference_file = csm_prescan/baseline_output.h5 +reference_images = spiral_3DRT_CSM_test.xml/image_0 +output_images = spiral_3DRT_CSM_test.xml/image_0 +value_comparison_threshold = 0.01 +scale_comparison_threshold = 0.01 + +[requirements] +system_memory = 8192 +gpu_support = 1 +gpu_memory = 8192 + +[tags] +tags = nhlbi,3drt + +[nhlbi] +description = Prescan CSM generation +noise_file = csm_prescan/noise_data.h5 +baseline_recon_time = 56.2 + diff --git a/test/nhlbi_integration_tests/cases/imoco_vds.cfg b/test/nhlbi_integration_tests/cases/imoco_vds.cfg new file mode 100644 index 0000000..e846bea --- /dev/null +++ b/test/nhlbi_integration_tests/cases/imoco_vds.cfg @@ -0,0 +1,36 @@ +[dependency.siemens] +data_file = imoco_vds/noise_data.h5 +measurement = 0 +additional_arguments = skip_converstion + +[dependency.client] +configuration = default_measurement_dependencies.xml + +[reconstruction.siemens] +data_file = imoco_vds/recon_data.h5 +measurement = 0 +additional_arguments = skip_converstion + +[reconstruction.client] +configuration = imoco_recon_vds.xml + +[reconstruction.test] +reference_file = imoco_vds/baseline_output.h5 +reference_images = imoco_recon_vds.xml/image_0 +output_images = imoco_recon_vds.xml/image_0 +value_comparison_threshold = 0.01 +scale_comparison_threshold = 0.01 + +[requirements] +system_memory = 8192 +gpu_support = 1 +gpu_memory = 8192 + +[tags] +tags = nhlbi,imoco + +[nhlbi] +description = iMOCO VDS 3D lung reconstruction +noise_file = imoco_vds/noise_data.h5 +baseline_recon_time = 386.0 + diff --git a/test/nhlbi_integration_tests/delete_test.py b/test/nhlbi_integration_tests/delete_test.py new file mode 100644 index 0000000..d011ad6 --- /dev/null +++ b/test/nhlbi_integration_tests/delete_test.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python3 +"""Remove an NHLBI integration test case. + +Deletes the .cfg file, removes manifest entries, and optionally deletes +the associated Azure Blob Storage data. + +Usage: + python delete_test.py imoco_vds + python delete_test.py imoco_vds --keep-data +""" + +import argparse +import sys +from pathlib import Path + +from get_nhlbi_data import get_container_client, load_manifest, save_manifest + +CASES_DIR = Path(__file__).parent / "cases" +BASELINES_DIR = Path(__file__).parent / "baselines" + + +def main(): + parser = argparse.ArgumentParser( + description="Remove an NHLBI integration test", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument('name', help="Test case name to delete") + parser.add_argument('--keep-data', action='store_true', + help="Keep Azure blobs; only remove local .cfg and manifest entries") + parser.add_argument('--yes', '-y', action='store_true', + help="Skip confirmation prompt") + + args = parser.parse_args() + name = args.name + + cfg_path = CASES_DIR / f"{name}.cfg" + manifest = load_manifest() + test_entries = [e for e in manifest if e.get('test') == name] + + if not cfg_path.exists() and not test_entries: + print(f"Error: Test '{name}' not found") + sys.exit(1) + + # Show what will be deleted + print(f"Test: {name}") + if cfg_path.exists(): + print(f" Config: {cfg_path}") + if test_entries: + print(f" Manifest entries: {len(test_entries)}") + for entry in test_entries: + print(f" - {entry['file']} ({entry.get('type', 'unknown')})") + if not args.keep_data and test_entries: + print(f" Azure blobs: {len(test_entries)} will be deleted") + else: + print(f" Azure blobs: kept") + + baseline_dir = BASELINES_DIR / name + if baseline_dir.exists(): + print(f" Local baselines: {baseline_dir}") + + # Confirm + if not args.yes: + action = "and Azure blobs" if not args.keep_data else "(keeping Azure data)" + response = input(f"\nDelete test '{name}' {action}? [y/N] ").strip().lower() + if response not in ('y', 'yes'): + print("Cancelled.") + sys.exit(0) + + # Delete Azure blobs + if not args.keep_data and test_entries: + try: + container_client = get_container_client() + for entry in test_entries: + blob_name = entry['file'] + try: + blob_client = container_client.get_blob_client(blob_name) + blob_client.delete_blob() + print(f" Deleted blob: {blob_name}") + except Exception as e: + print(f" Warning: Could not delete blob {blob_name}: {e}") + except Exception as e: + print(f" Warning: Could not connect to Azure: {e}") + print(" Local files will still be removed.") + + # Remove from manifest + remaining = [e for e in manifest if e.get('test') != name] + save_manifest(remaining) + print(f" Removed {len(test_entries)} manifest entries") + + # Delete .cfg + if cfg_path.exists(): + cfg_path.unlink() + print(f" Deleted {cfg_path}") + + # Delete local baseline directory + if baseline_dir.exists(): + import shutil + shutil.rmtree(baseline_dir) + print(f" Deleted {baseline_dir}") + + print(f"\nTest '{name}' removed successfully.") + + +if __name__ == '__main__': + main() diff --git a/test/nhlbi_integration_tests/generate_baseline.py b/test/nhlbi_integration_tests/generate_baseline.py new file mode 100644 index 0000000..af4b252 --- /dev/null +++ b/test/nhlbi_integration_tests/generate_baseline.py @@ -0,0 +1,554 @@ +#!/usr/bin/env python3 +"""Generate and validate baselines for NHLBI integration tests. + +Runs a reconstruction, generates preview artifacts, and prompts the user +to accept or reject the output as a baseline. Accepted baselines are +uploaded to Azure Blob Storage and registered in the manifest. + +Usage: + python generate_baseline.py --test imoco_vds [--port 9002] +""" + +import argparse +import configparser +import json +import os +import subprocess +import sys +import tempfile +import time +from datetime import datetime +from pathlib import Path +from get_nhlbi_data import download_data +import h5py +import numpy as np +from test_utils import get_gadgetron_bin_path +from get_nhlbi_data import ( + calc_sha256, + get_container_client, + load_manifest, + save_manifest, + upload_blob, +) +from test_utils import read_h5 +from collections import OrderedDict +import os.path as op + +CASES_DIR = Path(__file__).parent / "cases" +BASELINES_DIR = Path(__file__).parent / "baselines" + +# Ensure gadgetron binaries are on PATH +_gadgetron_bin = get_gadgetron_bin_path() +if _gadgetron_bin not in os.environ.get("PATH", ""): + os.environ["PATH"] = _gadgetron_bin + ":" + os.environ.get("PATH", "") + +# Import from existing integration test framework +sys.path.insert(0, str(Path(__file__).parent.parent / 'integration')) +from run_gadgetron_test import ( + send_data_to_gadgetron, + start_gadgetron_instance, + start_storage_server, +) + + +def get_data_dir(): + return os.environ.get('NHLBI_DATA_CACHE', str(Path(__file__).parent / 'data')) + + +def echo_handler(cmd): + print(' '.join(cmd)) + + +class GadgetronInstance: + def __init__(self, host, port): + self.host = host + self.port = port + + +def generate_preview(output_file, preview_dir): + """Generate text summary and optional PNG montage of reconstruction output.""" + os.makedirs(preview_dir, exist_ok=True) + summary_lines = [] + + try: + with h5py.File(output_file, 'r') as f: + summary_lines.append(f"File: {output_file}") + summary_lines.append(f"Groups: {list(f.keys())}") + + def visit_datasets(name, obj): + if isinstance(obj, h5py.Dataset): + summary_lines.append(f" Dataset: {name}") + summary_lines.append(f" Shape: {obj.shape}") + summary_lines.append(f" Dtype: {obj.dtype}") + if np.issubdtype(obj.dtype, np.number) and obj.size > 0: + data = obj[...] + if np.iscomplexobj(data): + data = np.abs(data) + summary_lines.append(f" Min: {np.min(data):.6e}") + summary_lines.append(f" Max: {np.max(data):.6e}") + summary_lines.append(f" Mean: {np.mean(data):.6e}") + summary_lines.append(f" Std: {np.std(data):.6e}") + + f.visititems(visit_datasets) + except Exception as e: + summary_lines.append(f"Error reading output: {e}") + + summary_text = '\n'.join(summary_lines) + summary_path = os.path.join(preview_dir, 'summary.txt') + with open(summary_path, 'w') as f: + f.write(summary_text) + + print("\n=== Baseline Preview ===") + print(summary_text) + print("========================\n") + + # Attempt PNG montage of central slices + try: + _generate_montage(output_file, preview_dir) + except Exception as e: + print(f"Note: Could not generate PNG montage: {e}") + print("Install matplotlib for visual previews: pip install matplotlib") + + return summary_path + + +def _generate_montage(output_file, preview_dir): + """Generate a PNG montage showing central slices from each dimension.""" + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + img_list,header_list=read_h5(output_file) + + if len(img_list) == 0: + print("No image datasets found for montage.") + return + for k in range(len(img_list)): + image_data = img_list[k].squeeze() + print(image_data.shape) + header= header_list[k] + print(header[0]) + img_serie_index=header[0]['image_series_index'] + if image_data.ndim == 2: + fig, ax = plt.subplots(1, 1, figsize=(6, 6)) + ax.imshow(image_data, cmap='gray') + ax.set_title('2D Output') + ax.axis('off') + elif image_data.ndim == 3: + nslices = image_data.shape[0] + # Show up to 9 evenly spaced slices + n_show = min(9, nslices) + indices = np.linspace(0, nslices - 1, n_show, dtype=int) + cols = min(3, n_show) + rows = (n_show + cols - 1) // cols + fig, axes = plt.subplots(rows, cols, figsize=(4 * cols, 4 * rows)) + axes = np.atleast_2d(axes) + for i, idx in enumerate(indices): + r, c = divmod(i, cols) + axes[r, c].imshow(image_data[idx], cmap='gray') + axes[r, c].set_title(f'Slice {idx}') + axes[r, c].axis('off') + for i in range(n_show, rows * cols): + r, c = divmod(i, cols) + axes[r, c].axis('off') + elif image_data.ndim >= 4: + # Show central slice of last two dims across first dimension + shape = image_data.shape + # Flatten to 3D: combine all leading dims + idx_0=image_data.shape[0] + flat = image_data.reshape(-1, shape[-2], shape[-1],order="C") + nslices = flat.shape[0] + central_slices=(nslices/idx_0) // 2 + indices = np.arange(central_slices,nslices,idx_0).astype(np.int32) + n_show = min(9, len(indices)) + if len(indices) > n_show: + indices = indices[:n_show] + cols = min(3, n_show) + rows = (n_show + cols - 1) // cols + fig, axes = plt.subplots(rows, cols, figsize=(4 * cols, 4 * rows)) + axes = np.atleast_2d(axes) + for i, idx in enumerate(indices): + r, c = divmod(i, cols) + axes[r, c].imshow(flat[idx], cmap='gray') + axes[r, c].set_title(f'Frame {idx}') + axes[r, c].axis('off') + for i in range(n_show, rows * cols): + r, c = divmod(i, cols) + axes[r, c].axis('off') + else: + print("Data is 1D or scalar, skipping montage") + continue + + fig.suptitle(f'Baseline Preview image {img_serie_index}', fontsize=14) + fig.tight_layout() + preview_path = os.path.join(preview_dir, f'preview_{img_serie_index}.png') + fig.savefig(preview_path, dpi=150, bbox_inches='tight') + plt.close(fig) + print(f"Preview montage saved: {preview_path}") + + +def run_reconstruction(test_name, port, storage_port): + """Run noise dependency + reconstruction and return the output file path.""" + cfg_path = CASES_DIR / f"{test_name}.cfg" + if not cfg_path.exists(): + print(f"Error: Test case not found: {cfg_path}") + sys.exit(1) + + config = configparser.ConfigParser() + config.read_dict({ + "DEFAULT": { + 'parameter_xml': 'IsmrmrdParameterMap_Siemens.xml', + 'parameter_xsl': 'IsmrmrdParameterMap_Siemens.xsl', + 'value_comparison_threshold': '0.01', + 'scale_comparison_threshold': '0.01', + } + }) + config.read(cfg_path) + + data_dir = get_data_dir() + + test_dir = str(BASELINES_DIR / test_name) + os.makedirs(test_dir, exist_ok=True) + + # Resolve data file paths + noise_file = os.path.join(data_dir, config['dependency.siemens']['data_file']) + recon_file = os.path.join(data_dir, config['reconstruction.siemens']['data_file']) + + if not os.path.isfile(noise_file): + print(f"Error: Noise file not found: {noise_file}") + print("Run 'python get_nhlbi_data.py download --test {}' first.".format(test_name)) + sys.exit(1) + if not os.path.isfile(recon_file): + print(f"Error: Recon data file not found: {recon_file}") + print("Run 'python get_nhlbi_data.py download --test {}' first.".format(test_name)) + sys.exit(1) + + # additional dependency files + + add_dep_sections_number=[dep_sec.split(".")[-1] for dep_sec in config.sections() if "dependency.siemens." in dep_sec] + + for num_dep in add_dep_sections_number: + dep_file=os.path.join(data_dir,config[f"dependency.siemens.{num_dep}"]['data_file']) + if not os.path.isfile(dep_file): + print(f"Error: Recon data file not found: {dep_file}") + print("Run 'python get_nhlbi_data.py download --test {}' first.".format(test_name)) + sys.exit(1) + + output_file = os.path.join(test_dir, 'output.h5') + + noise_config = config['dependency.client']['configuration'] + recon_config = config['reconstruction.client']['configuration'] + + gadgetron_instance = GadgetronInstance("localhost", str(port)) + + with tempfile.TemporaryDirectory() as storage_folder: + storage_log = open(os.path.join(test_dir, 'storage.log'), 'w') + try: + storage_proc = start_storage_server( + log=storage_log, + port=str(storage_port), + storage_folder=storage_folder, + ) + except Exception as e: + storage_log.close() + print(f"Error starting storage server: {e}") + sys.exit(1) + + try: + gt_log_out = open(os.path.join(test_dir, 'gadgetron.log.out'), 'w') + gt_log_err = open(os.path.join(test_dir, 'gadgetron.log.err'), 'w') + storage_address = f"http://localhost:{storage_port}" + + gt_proc = start_gadgetron_instance( + log_stdout=gt_log_out, + log_stderr=gt_log_err, + port=str(port), + storage_address=storage_address, + ) + + try: + # Wait briefly for gadgetron to start + time.sleep(2) + + if len(add_dep_sections_number)>0: + print("\n--- Sending additional dependencies ---") + print(add_dep_sections_number) + for num_dep in add_dep_sections_number: + dep_siem=f"dependency.siemens.{num_dep}" + input_dep_file=os.path.join(data_dir,config[f"dependency.siemens.{num_dep}"]['data_file']) + input_config_file=config[f"dependency.client.{num_dep}"]['configuration'] + dep_name=f"dep_{num_dep}" + dep_log = open(os.path.join(test_dir, f'{dep_name}.log'), 'w') + print("\n--- Sending dependency {input_dep_file} with config {input_config_file} ---") + send_data_to_gadgetron( + echo_handler, gadgetron_instance, + input=input_dep_file, + output=os.path.join(test_dir, f'{dep_name}_output.h5'), + configuration=['-c', input_config_file], + group=noise_config, + log=dep_log, + additional_arguments=config[f"dependency.siemens.{num_dep}"].get('additional_arguments'), + ) + dep_log.close() + # Send noise data + print(f"\n--- Sending noise data ({noise_config}) ---") + noise_log = open(os.path.join(test_dir, 'noise_client.log'), 'w') + send_data_to_gadgetron( + echo_handler, gadgetron_instance, + input=noise_file, + output=os.path.join(test_dir, 'noise_output.h5'), + configuration=['-c', noise_config], + group=noise_config, + log=noise_log, + additional_arguments=config['dependency.siemens'].get('additional_arguments'), + ) + noise_log.close() + + # Send reconstruction data + print(f"\n--- Sending reconstruction data ({recon_config}) ---") + recon_log = open(os.path.join(test_dir, 'recon_client.log'), 'w') + start_time = time.time() + send_data_to_gadgetron( + echo_handler, gadgetron_instance, + input=recon_file, + output=output_file, + configuration=['-c', recon_config], + group=recon_config, + log=recon_log, + additional_arguments=config['reconstruction.siemens'].get('additional_arguments'), + ) + recon_log.close() + elapsed = time.time() - start_time + print(f"Reconstruction completed in {elapsed:.1f}s") + + finally: + gt_proc.kill() + gt_log_out.close() + gt_log_err.close() + finally: + storage_proc.kill() + storage_log.close() + + if not os.path.isfile(output_file): + print("Error: No output file was produced.") + print(f"Check logs in {test_dir}/") + sys.exit(1) + + return output_file, elapsed + + +def update_cfg(cfg_path, baseline_file, recon_time): + """Update the .cfg file with the new baseline reference and recon time.""" + test_name=op.basename(cfg_path).replace('.cfg','') + config = configparser.ConfigParser() + config.read(cfg_path) + images,headers=read_h5(baseline_file) + print(config.sections()) + reconstruction_tests=[section_name for section_name in config.sections() if section_name.startswith('reconstruction.test')] + value_comparison_threshold_initial=config[reconstruction_tests[0]]['value_comparison_threshold'] + scale_comparison_threshold_initial=config[reconstruction_tests[0]]['scale_comparison_threshold'] + for section_name in reconstruction_tests: + config.remove_section(section_name) + + key_names=[f'reconstruction.test.{i+1}' for i in range(len(images))] + if len(key_names)==1: + key_names=['reconstruction.test'] + for key_name,image,header in zip(key_names,images,headers): + image_serie_index=header[0]['image_series_index'] + config[key_name] = { + 'reference_file': f'{test_name}/baseline_output.h5', + 'reference_images': f"{config['reconstruction.client']['configuration']}/image_{image_serie_index}", + 'output_images': f"{config['reconstruction.client']['configuration']}/image_{image_serie_index}", + 'value_comparison_threshold': value_comparison_threshold_initial, + 'scale_comparison_threshold': scale_comparison_threshold_initial, + } + config['nhlbi']['baseline_recon_time'] = str(round(recon_time, 1)) + # Sorting keys to ensure deterministic order in the .cfg file + desired_order = [section_name for section_name in config.sections() if section_name.startswith('dependency')] + desired_order.extend(['reconstruction.siemens','reconstruction.client']) + desired_order.extend([section_name for section_name in config.sections() if section_name.startswith('reconstruction.test')]) + desired_order.extend(['requirements','tags','nhlbi']) + + ordered_config = configparser.ConfigParser() + + section_order = config.sections() + for section in desired_order: + if config.has_section(section): + ordered_config.add_section(section) + for key, value in config.items(section): + ordered_config.set(section, key, value) + + # Optionally, add any remaining sections not in the desired order + for section in config.sections(): + if section not in desired_order: + print(f"Warning: Section '{section}' not in desired order list, adding at the end.") + ordered_config.add_section(section) + for key, value in config.items(section): + ordered_config.set(section, key, value) + + with open(cfg_path, 'w') as f: + ordered_config.write(f) + +def main(): + parser = argparse.ArgumentParser( + description="Generate and validate NHLBI test baselines", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument('--test', required=True, help="Test case name") + parser.add_argument('--port', type=int, default=9003, help="Gadgetron port") + parser.add_argument('--storage-port', type=int, default=9113, help="Storage server port") + parser.add_argument('--auto-accept', action='store_true', + help="Accept baseline without interactive prompt (for CI)") + parser.add_argument('--skip-upload', action='store_true', + help="Skip uploading baseline to Azure") + parser.add_argument('--accept-existing', action='store_true', + help="Accept an existing output in baselines// without re-running reconstruction") + parser.add_argument('--recon-time', type=float, default=None, + help="Override reconstruction time in seconds (use with --accept-existing)") + + args = parser.parse_args() + + test_name = args.test + preview_dir = str(BASELINES_DIR / test_name) + + if args.accept_existing: + # Accept an already-generated baseline without re-running + output_file = str(BASELINES_DIR / test_name / 'output.h5') + if not os.path.isfile(output_file): + print(f"Error: No existing output at {output_file}") + print("Run without --accept-existing to generate it first.") + sys.exit(1) + + # Try to extract recon time from gadgetron log + recon_time = args.recon_time + if recon_time is None: + recon_time = _extract_recon_time_from_log(test_name) + if recon_time is None: + recon_time = 0.0 + print("Warning: Could not determine reconstruction time. Use --recon-time to set it.") + + print(f"Using existing output: {output_file}") + print(f"Reconstruction time: {recon_time:.1f}s") + generate_preview(output_file, preview_dir) + else: + # Download test data if needed + print(f"Ensuring test data is available for '{test_name}'...") + dl_args = argparse.Namespace( + destination=get_data_dir(), + test=test_name, + list=str(Path(__file__).parent / 'nhlbi_data.json'), + ) + try: + download_data(dl_args) + except Exception as e: + print(f"Warning: Could not download data: {e}") + print("Continuing with locally available data...") + + # Run reconstruction + print(f"\nRunning reconstruction for '{test_name}'...") + output_file, recon_time = run_reconstruction(test_name, args.port, args.storage_port) + + # Generate preview + generate_preview(output_file, preview_dir) + + # Interactive validation + if args.auto_accept: + accept = True + else: + print(f"\nOutput file: {output_file}") + print(f"Preview dir: {preview_dir}/") + print(f"Reconstruction time: {recon_time:.1f}s") + response = input("Accept this output as baseline? [y/N] ").strip().lower() + accept = response in ('y', 'yes') + + if not accept: + print("\nBaseline rejected. Output kept for inspection at:") + print(f" {preview_dir}/") + print(f"\nTo re-run: python generate_baseline.py --test {test_name}") + sys.exit(0) + + # Upload baseline + baseline_sha256 = calc_sha256(output_file) + + if not args.skip_upload: + print("\nUploading baseline to Azure...") + container_client = get_container_client() + upload_blob(container_client, output_file, f"{test_name}/baseline_output.h5") + else: + print("Skipping Azure upload (--skip-upload)") + + # Update manifest + manifest = load_manifest() + # Remove any existing baseline entry for this test + manifest = [e for e in manifest if not (e.get('test') == test_name and e.get('type') == 'baseline')] + manifest.append({ + 'file': f'{test_name}/baseline_output.h5', + 'sha256': baseline_sha256, + 'type': 'baseline', + 'test': test_name, + 'validated_by': os.environ.get('USER', 'unknown'), + 'validated_date': datetime.now().strftime('%Y-%m-%d'), + 'git_sha': _get_git_sha(), + 'recon_time_seconds': round(recon_time, 1), + }) + save_manifest(manifest) + + # Update .cfg with reference_file and baseline timing + + cfg_path = CASES_DIR / f"{test_name}.cfg" + update_cfg(cfg_path, output_file, recon_time) + + """config = configparser.ConfigParser() + config.read(cfg_path) + + # Set the reference file to point to the baseline output in Azure Blob Storage + + + config['reconstruction.test']['reference_file'] = f'{test_name}/baseline_output.h5' + config['nhlbi']['baseline_recon_time'] = str(round(recon_time, 1)) + with open(cfg_path, 'w') as f: + config.write(f) + """ + print(f"\nBaseline accepted and registered for '{test_name}'") + print(f"Baseline reconstruction time: {recon_time:.1f}s") + print(f"Run 'python run_nhlbi_tests.py cases/{test_name}.cfg' to verify") + + +def _extract_recon_time_from_log(test_name): + """Try to extract reconstruction time from gadgetron server log timestamps.""" + log_path = BASELINES_DIR / test_name / 'gadgetron.log.err' + if not log_path.exists(): + return None + try: + import re + timestamps = [] + with open(log_path, 'r') as f: + for line in f: + m = re.match(r'^(\d{2}-\d{2} \d{2}:\d{2}:\d{2}\.\d+)', line) + if m: + timestamps.append(m.group(1)) + if len(timestamps) >= 2: + from datetime import datetime + fmt = '%m-%d %H:%M:%S.%f' + start = datetime.strptime(timestamps[0], fmt) + end = datetime.strptime(timestamps[-1], fmt) + elapsed = (end - start).total_seconds() + if elapsed > 0: + return elapsed + except Exception: + pass + return None + + +def _get_git_sha(): + try: + result = subprocess.run( + ['git', 'rev-parse', '--short', 'HEAD'], + capture_output=True, text=True, cwd=Path(__file__).parent, + ) + return result.stdout.strip() if result.returncode == 0 else 'unknown' + except Exception: + return 'unknown' + + +if __name__ == '__main__': + main() diff --git a/test/nhlbi_integration_tests/get_nhlbi_data.py b/test/nhlbi_integration_tests/get_nhlbi_data.py new file mode 100644 index 0000000..299960a --- /dev/null +++ b/test/nhlbi_integration_tests/get_nhlbi_data.py @@ -0,0 +1,170 @@ +#!/usr/bin/env python3 +"""Download and upload NHLBI test data from/to private Azure Blob Storage. + +Authentication uses DefaultAzureCredential (picks up `az login` for developers, +managed identity or service principal for CI). Fallback: set NHLBI_AZURE_SAS_TOKEN +environment variable for SAS-token-based access. + +Dependencies: pip install azure-storage-blob azure-identity +""" + +import argparse +import hashlib +import json +import os +import sys +from concurrent.futures import ThreadPoolExecutor +from pathlib import Path + +STORAGE_ACCOUNT = "gadgetrondata" +CONTAINER_NAME = "nhlbitestdata" +ACCOUNT_URL = f"https://{STORAGE_ACCOUNT}.blob.core.windows.net" +DEFAULT_DATA_DIR = "data" +MANIFEST_FILE = Path(__file__).parent / "nhlbi_data.json" + + +def calc_sha256(filepath): + sha256 = hashlib.sha256() + with open(filepath, 'rb') as f: + for chunk in iter(lambda: f.read(65536), b''): + sha256.update(chunk) + return sha256.hexdigest() + + +def is_valid(filepath, expected_sha256): + if not os.path.isfile(filepath): + return False + return expected_sha256 == calc_sha256(filepath) + + +def get_container_client(): + sas_token = os.environ.get("NHLBI_AZURE_SAS_TOKEN") + if sas_token: + from azure.storage.blob import ContainerClient + return ContainerClient( + account_url=ACCOUNT_URL, + container_name=CONTAINER_NAME, + credential=sas_token, + ) + else: + from azure.identity import DefaultAzureCredential + from azure.storage.blob import ContainerClient + credential = DefaultAzureCredential() + return ContainerClient( + account_url=ACCOUNT_URL, + container_name=CONTAINER_NAME, + credential=credential, + ) + + +def download_blob_public(blob_name, destination, retries=3): + """Download a blob via public URL (no auth required).""" + import urllib.request + import urllib.error + import socket + + url = f"{ACCOUNT_URL}/{CONTAINER_NAME}/{blob_name}" + os.makedirs(os.path.dirname(destination), exist_ok=True) + for attempt in range(retries): + try: + with urllib.request.urlopen(url, timeout=60) as response: + with open(destination, 'wb') as f: + for chunk in iter(lambda: response.read(1024 * 1024), b''): + f.write(chunk) + return + except (urllib.error.URLError, ConnectionResetError, socket.timeout) as e: + if attempt == retries - 1: + raise RuntimeError(f"Failed to download {blob_name} after {retries} attempts: {e}") + print(f"Retry {attempt + 1} for {blob_name}: {e}") + + +def upload_blob(container_client, local_path, blob_name): + blob_client = container_client.get_blob_client(blob_name) + print(f"Uploading {local_path} -> {blob_name}") + with open(local_path, 'rb') as f: + blob_client.upload_blob(f, overwrite=True) + print(f"Upload complete: {blob_name}") + + +def load_manifest(): + with open(MANIFEST_FILE, 'r') as f: + return json.load(f) + + +def save_manifest(entries): + with open(MANIFEST_FILE, 'w') as f: + json.dump(entries, f, indent=2) + f.write('\n') + + +def download_data(args): + entries = load_manifest() + + if args.test: + entries = [e for e in entries if e.get('test') == args.test] + if not entries: + print(f"No data entries found for test '{args.test}'") + sys.exit(1) + + data_dir = args.destination + + def download_entry(entry): + destination = os.path.join(data_dir, entry['file']) + if is_valid(destination, entry['sha256']): + print(f"Verified: {destination}") + return + print(f"Downloading: {entry['file']}") + download_blob_public(entry['file'], destination) + if not is_valid(destination, entry['sha256']): + actual = calc_sha256(destination) + raise RuntimeError( + f"Downloaded file {destination} failed validation. " + f"Expected SHA256 {entry['sha256']}. Actual SHA256 {actual}" + ) + print(f"Saved: {destination}") + + with ThreadPoolExecutor(max_workers=4) as executor: + list(executor.map(download_entry, entries)) + + +def upload_data(args): + container_client = get_container_client() + upload_blob(container_client, args.local_path, args.remote_path) + sha256 = calc_sha256(args.local_path) + print(f"SHA256: {sha256}") + return sha256 + + +def main(): + parser = argparse.ArgumentParser( + description="NHLBI Integration Test Data Manager", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + subparsers = parser.add_subparsers(dest='command', help='Command to run') + + dl_parser = subparsers.add_parser('download', help='Download test data') + dl_parser.add_argument('-d', '--destination', type=str, + default=os.environ.get('NHLBI_DATA_CACHE', DEFAULT_DATA_DIR), + help="Local folder for downloaded data") + dl_parser.add_argument('-t', '--test', type=str, default=None, + help="Download data only for the specified test case") + dl_parser.add_argument('-l', '--list', type=str, default=str(MANIFEST_FILE), + help="Path to data manifest file") + + ul_parser = subparsers.add_parser('upload', help='Upload data to Azure') + ul_parser.add_argument('local_path', type=str, help="Local file to upload") + ul_parser.add_argument('remote_path', type=str, help="Blob path in container") + + args = parser.parse_args() + + if args.command == 'download': + download_data(args) + elif args.command == 'upload': + upload_data(args) + else: + parser.print_help() + sys.exit(1) + + +if __name__ == '__main__': + main() diff --git a/test/nhlbi_integration_tests/list_tests.py b/test/nhlbi_integration_tests/list_tests.py new file mode 100644 index 0000000..39c02aa --- /dev/null +++ b/test/nhlbi_integration_tests/list_tests.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +"""List all registered NHLBI integration tests and their status.""" + +import configparser +import json +import sys +from pathlib import Path + +from get_nhlbi_data import MANIFEST_FILE + +CASES_DIR = Path(__file__).parent / "cases" +BASELINES_DIR = Path(__file__).parent / "baselines" + +_codes = { + 'red': '\033[91m', + 'green': '\033[92m', + 'cyan': '\033[96m', + 'yellow': '\033[93m', + 'bold': '\033[1m', + 'end': '\033[0m', +} + + +def color(text, c): + return f"{_codes.get(c, '')}{text}{_codes.get('end', '')}" + + +def main(): + cfg_files = sorted(CASES_DIR.glob("*.cfg")) + + if not cfg_files: + print("No tests registered. Use submit_test.py to add one.") + sys.exit(0) + + # Load manifest + try: + with open(MANIFEST_FILE, 'r') as f: + manifest = json.load(f) + except Exception: + manifest = [] + + def get_manifest_entry(test_name, entry_type): + return next((e for e in manifest if e.get('test') == test_name and e.get('type') == entry_type), None) + + print(color(f"{'Test':<25} {'Config':<40} {'Baseline':<15} {'Recon Time':<12} {'Tags'}", 'bold')) + print("-" * 110) + + for cfg_file in cfg_files: + name = cfg_file.stem + config = configparser.ConfigParser() + config.read(cfg_file) + + # Config XML + recon_config = config.get('reconstruction.client', 'configuration', fallback='?') + + # Tags + tags = config.get('tags', 'tags', fallback='') + + # Baseline status + baseline_entry = get_manifest_entry(name, 'baseline') + has_local_output = (BASELINES_DIR / name / 'output.h5').exists() + + if baseline_entry: + validated_by = baseline_entry.get('validated_by', '?') + validated_date = baseline_entry.get('validated_date', '?') + baseline_status = color(f"yes ({validated_date})", 'green') + elif has_local_output: + baseline_status = color("local only", 'yellow') + else: + baseline_status = color("missing", 'red') + + # Recon time + recon_time = config.get('nhlbi', 'baseline_recon_time', fallback=None) + if recon_time: + recon_time_str = f"{float(recon_time):.1f}s" + elif baseline_entry and baseline_entry.get('recon_time_seconds'): + recon_time_str = f"{baseline_entry['recon_time_seconds']:.1f}s" + else: + recon_time_str = "-" + + # Description + description = config.get('nhlbi', 'description', fallback='') + + print(f"{name:<25} {recon_config:<40} {baseline_status:<27} {recon_time_str:<12} {tags}") + if description: + print(f" {color(description, 'cyan')}") + + # Data summary + noise_count = sum(1 for e in manifest if e.get('type') == 'noise') + input_count = sum(1 for e in manifest if e.get('type') == 'input') + baseline_count = sum(1 for e in manifest if e.get('type') == 'baseline') + print(f"\n{len(cfg_files)} tests, {baseline_count} baselines, " + f"{noise_count + input_count} data files in manifest") + + +if __name__ == '__main__': + main() diff --git a/test/nhlbi_integration_tests/nhlbi_data.json b/test/nhlbi_integration_tests/nhlbi_data.json new file mode 100644 index 0000000..3910383 --- /dev/null +++ b/test/nhlbi_integration_tests/nhlbi_data.json @@ -0,0 +1,80 @@ +[ + { + "file": "imoco_vds/noise_data.h5", + "sha256": "89396f6de4bf248aa9a37335ac086b73dd499773ec1f57c9aba83879b1cc144d", + "type": "noise", + "test": "imoco_vds" + }, + { + "file": "imoco_vds/recon_data.h5", + "sha256": "4d1888c4714467cbc4a0f7ccb454275d10596ffcd64b003099fdf6e110e79abf", + "type": "input", + "test": "imoco_vds" + }, + { + "file": "imoco_vds/baseline_output.h5", + "sha256": "a3080275eb0b509c917bf8398cad3a19cbff9f5ad5555e99cd2bf5ec4eec410d", + "type": "baseline", + "test": "imoco_vds", + "validated_by": "unknown", + "validated_date": "2026-03-19", + "git_sha": "d11d8915", + "recon_time_seconds": 386.0 + }, + { + "file": "csm_prescan/noise_data.h5", + "sha256": "8afbf9f5b988b110bdc164fba28205dbd1caf062ce8b7681dbc6034a2e8818fb", + "type": "noise", + "test": "csm_prescan" + }, + { + "file": "csm_prescan/recon_data.h5", + "sha256": "41dba6a2aff581500c8266bbc24a692d487f3991ca86cce67f68a54bf8a2f18f", + "type": "input", + "test": "csm_prescan" + }, + { + "file": "csm_prescan/baseline_output.h5", + "sha256": "20358f20e037ee1310dfeb097904e45ac9b6af686649027e67d55d349d5c4917", + "type": "baseline", + "test": "csm_prescan", + "validated_by": "unknown", + "validated_date": "2026-08-12", + "git_sha": "f6e3df1", + "recon_time_seconds": 56.2 + }, + { + "file": "3drt/noise_data.h5", + "sha256": "b82d79de3cd8e2f489d24b87f6bbcc0f66db3ea63f75a2f83bc2c974d7e14458", + "type": "noise", + "test": "3drt" + }, + { + "file": "3drt/recon_data.h5", + "sha256": "6d7125a55eb4c3d9b0b4e8bfd9d2ad793c18d2728afe3a887c6f7d579b723d80", + "type": "input", + "test": "3drt" + }, + { + "file": "3drt/prescan-reference_noise.h5", + "sha256": "8afbf9f5b988b110bdc164fba28205dbd1caf062ce8b7681dbc6034a2e8818fb", + "type": "additional", + "test": "3drt" + }, + { + "file": "3drt/prescan-reference.h5", + "sha256": "41dba6a2aff581500c8266bbc24a692d487f3991ca86cce67f68a54bf8a2f18f", + "type": "additional", + "test": "3drt" + }, + { + "file": "3drt/baseline_output.h5", + "sha256": "130da444d2d368a9f679191d28932a8f023ff684f5dd6319f1d4a06d9600977a", + "type": "baseline", + "test": "3drt", + "validated_by": "unknown", + "validated_date": "2026-08-28", + "git_sha": "87e65a6", + "recon_time_seconds": 42.3 + } +] diff --git a/test/nhlbi_integration_tests/requirements.txt b/test/nhlbi_integration_tests/requirements.txt new file mode 100644 index 0000000..309a806 --- /dev/null +++ b/test/nhlbi_integration_tests/requirements.txt @@ -0,0 +1,6 @@ +azure-storage-blob +azure-identity +h5py +ismrmrd +numpy +matplotlib diff --git a/test/nhlbi_integration_tests/run_nhlbi_tests.py b/test/nhlbi_integration_tests/run_nhlbi_tests.py new file mode 100644 index 0000000..38c4907 --- /dev/null +++ b/test/nhlbi_integration_tests/run_nhlbi_tests.py @@ -0,0 +1,327 @@ +#!/usr/bin/env python3 +"""NHLBI Integration Test Orchestrator. + +Thin wrapper around the existing run_gadgetron_test.py that handles +NHLBI-specific concerns: downloading test data from private Azure Blob +Storage and filtering by NHLBI tags. + +Usage: + python run_nhlbi_tests.py cases/*.cfg + python run_nhlbi_tests.py cases/imoco_vds.cfg + python run_nhlbi_tests.py cases/*.cfg --only fast +""" + +import argparse +import configparser +import csv +import glob +import itertools +import json +import os +import subprocess +import sys +from pathlib import Path +from test_utils import get_gadgetron_bin_path +# Ensure gadgetron binaries are on PATH +_gadgetron_bin = get_gadgetron_bin_path() +if _gadgetron_bin not in os.environ.get("PATH", ""): + os.environ["PATH"] = _gadgetron_bin + ":" + os.environ.get("PATH", "") + +# Reuse tag/requirement parsing from the existing test runner. +# In dev: test/nhlbi_integration_tests/../integration i.e. test/integration +# In RT container: /opt/nhlbi-integration-test -> sibling is /opt/integration-test +_dev_integration = Path(__file__).parent.parent / 'integration' +_rt_integration = Path('/opt/integration-test') +INTEGRATION_DIR = _dev_integration if _dev_integration.is_dir() else _rt_integration +sys.path.insert(0, str(INTEGRATION_DIR)) +from run_tests import ( + _colors_disabled, + _colors_enabled, + output_csv, + output_log_file, + query_gadgetron_capabilities, + ignore_gadgetron_capabilities, + read_test_details, + should_skip_test, + split_tag_list, +) + +SCRIPT_DIR = Path(__file__).parent +RUN_TEST_SCRIPT = INTEGRATION_DIR / 'run_gadgetron_test.py' + + +def get_data_dir(): + return os.environ.get('NHLBI_DATA_CACHE', str(SCRIPT_DIR / 'data')) + + +def download_test_data(test_names): + """Download data for the specified tests from Azure Blob Storage.""" + from get_nhlbi_data import download_data + import argparse as _argparse + + data_dir = get_data_dir() + for name in test_names: + dl_args = _argparse.Namespace( + destination=data_dir, + test=name, + list=str(SCRIPT_DIR / 'nhlbi_data.json'), + ) + try: + download_data(dl_args) + except Exception as e: + print(f"Warning: Could not download data for test '{name}': {e}") + + +def get_test_name_from_cfg(cfg_path): + """Extract the test name from a .cfg file path.""" + return Path(cfg_path).stem + + +def check_baseline_exists(cfg_path): + """Check if the test has a baseline registered in the manifest.""" + test_name = get_test_name_from_cfg(cfg_path) + manifest_path = SCRIPT_DIR / 'nhlbi_data.json' + try: + with open(manifest_path, 'r') as f: + manifest = json.load(f) + return any( + e.get('test') == test_name and e.get('type') == 'baseline' + for e in manifest + ) + except Exception: + return False + + +def get_baseline_recon_time(cfg_path): + """Read the baseline reconstruction time from the .cfg file.""" + config = configparser.ConfigParser() + config.read(cfg_path) + try: + return float(config['nhlbi']['baseline_recon_time']) + except (KeyError, ValueError): + return None + + +def check_speed_regression(test_file, actual_time, color_handler, speed_threshold): + """Compare actual reconstruction time against baseline and report.""" + baseline_time = get_baseline_recon_time(test_file) + if baseline_time is None: + return None + + ratio = actual_time / baseline_time + pct_change = (ratio - 1.0) * 100 + + if ratio > speed_threshold: + print(color_handler( + f" SPEED REGRESSION: {actual_time:.1f}s vs baseline {baseline_time:.1f}s " + f"({pct_change:+.1f}%, threshold {(speed_threshold - 1) * 100:.0f}%)", + 'red', + )) + return False + elif pct_change < -5: + print(color_handler( + f" Speed improved: {actual_time:.1f}s vs baseline {baseline_time:.1f}s ({pct_change:+.1f}%)", + 'green', + )) + else: + print(f" Speed OK: {actual_time:.1f}s vs baseline {baseline_time:.1f}s ({pct_change:+.1f}%)") + + return True + + +def main(): + parser = argparse.ArgumentParser( + description="NHLBI Integration Test Runner", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + + parser.add_argument('-p', '--port', type=int, default=9003, help="Port for Gadgetron instance") + parser.add_argument('-a', '--host', type=str, default="localhost", help="Address of Gadgetron host") + + parser.add_argument('-e', '--external', action='store_const', const=['-e'], default=[], + help="Use external Gadgetron; don't start a new instance each test.") + + parser.add_argument('-d', '--data-folder', type=str, default=None, + help="Look for test data in the specified folder (default: NHLBI_DATA_CACHE or ./data)") + parser.add_argument('-t', '--test-folder', type=str, default='test', + help="Save Gadgetron and Client output to specified folder") + + parser.add_argument('-F', '--ignore-failures', action='store_true', default=False, + help="Continue running tests after failures") + parser.add_argument('-s', '--stats', type=str, default=None, + help="Output individual test stats to CSV file") + + parser.add_argument('--timeout', type=int, default=None, + help="Fail test if it runs longer than timeout seconds") + + parser.add_argument('--echo-log-on-failure', action='store_true', default=False, + help="Send test logs to stdout on failure") + + parser.add_argument('--disable-color', dest='color_handler', action='store_const', + const=_colors_disabled, default=_colors_enabled, + help="Disable colors in output") + + parser.add_argument('--disable-capability-query', action='store_const', + dest='capability_query_function', + const=ignore_gadgetron_capabilities, + default=query_gadgetron_capabilities, + help="Disable querying Gadgetron capabilities") + + parser.add_argument('--ignore-requirements', type=split_tag_list, default='none', metavar='tags', + help="Run tests with specified tags regardless of capabilities") + parser.add_argument('--only', type=split_tag_list, default='all', metavar='tags', + help="Only run tests with the specified tags") + parser.add_argument('--exclude', type=split_tag_list, default='none', metavar='tags', + help="Do not run tests with the specified tags") + + parser.add_argument('--skip-download', action='store_true', default=False, + help="Skip automatic data download from Azure") + + parser.add_argument('--speed-threshold', type=float, default=1.5, + help="Fail if reconstruction takes longer than this multiple of baseline time (e.g., 1.5 = 50%% slower)") + parser.add_argument('--no-speed-check', action='store_true', default=False, + help="Disable speed regression checking") + + parser.add_argument('tests', type=str, nargs='+', help="Test case .cfg files or glob patterns") + + args = parser.parse_args() + + data_dir = args.data_folder or get_data_dir() + + # Resolve test files + files = sorted(set(itertools.chain(*[glob.glob(pattern) for pattern in args.tests]))) + if not files: + print("No test files found matching the specified patterns.") + sys.exit(1) + + # Check for missing baselines + missing_baselines = [] + for f in files: + if not check_baseline_exists(f): + missing_baselines.append(f) + print(args.color_handler( + f"Warning: No baseline for {f} — test will be skipped", + 'cyan', + )) + + # Filter out tests without baselines + files = [f for f in files if f not in missing_baselines] + if not files: + print("No tests with baselines to run.") + sys.exit(0) + + # Download test data + if not args.skip_download: + test_names = [get_test_name_from_cfg(f) for f in files] + print("Downloading test data...") + download_test_data(test_names) + + # Read test details and filter by capabilities/tags + tests = [read_test_details(f) for f in files] + capabilities = args.capability_query_function(args) + + stats = [] + passed = [] + failed = [] + skipped = [] + speed_regressions = [] + + def skip_handler(test, message): + skipped.append((test, message)) + + tests = [t for t in tests if not should_skip_test(t, capabilities, args, skip_handler)] + + if skipped: + print("\nSkipped tests:") + for test, message in skipped: + print(f"\t{test.get('file')} ({message})") + + # Run each test + for i, test in enumerate(tests, start=1): + print(args.color_handler(f"\nTest {i} of {len(tests)}: {test.get('file')}\n", 'bold')) + + disable_color = ['--disable-colors'] if args.color_handler == _colors_disabled else [] + + command = [ + sys.executable, str(RUN_TEST_SCRIPT), + '-a', str(args.host), + '-d', str(data_dir), + '-t', str(args.test_folder), + '-p', str(args.port), + ] + args.external + disable_color + [test.get('file')] + + with subprocess.Popen(command) as proc: + try: + import time as _time + test_start = _time.time() + proc.wait(timeout=args.timeout) + test_elapsed = _time.time() - test_start + + if proc.returncode == 0: + passed.append(test) + try: + with open('test/stats.json') as sf: + stat = json.loads(sf.read()) + stats.append(stat) + actual_time = stat.get('processing_time', test_elapsed) + except FileNotFoundError: + actual_time = test_elapsed + + # Check speed regression + if not args.no_speed_check: + speed_ok = check_speed_regression( + test.get('file'), actual_time, + args.color_handler, args.speed_threshold, + ) + if speed_ok is False: + speed_regressions.append(test) + else: + if args.echo_log_on_failure: + for log in glob.glob(os.path.join(args.test_folder, '*.log')): + output_log_file(log) + failed.append(test) + try: + with open('test/stats.json') as sf: + stats.append(json.loads(sf.read())) + except FileNotFoundError: + pass + if not args.ignore_failures: + break + except subprocess.TimeoutExpired: + print(f"Timeout during test: {test.get('file')}") + proc.kill() + failed.append(test) + if not args.ignore_failures: + break + + if args.stats and stats: + output_csv(stats, args.stats) + + # Summary + if failed: + print("\nFailed tests:") + for test in failed: + print(f"\t{test.get('file')}") + + if speed_regressions: + print("\nSpeed regressions:") + for test in speed_regressions: + print(f"\t{test.get('file')}") + + if missing_baselines: + print("\nTests skipped (no baseline):") + for f in missing_baselines: + print(f"\t{f}") + + print(f"\n{len(passed)} tests passed. {len(failed)} tests failed. " + f"{len(skipped)} tests skipped. {len(missing_baselines)} missing baselines. " + f"{len(speed_regressions)} speed regressions.") + + if stats: + print(f"Total processing time: {sum(s['processing_time'] for s in stats):.2f} seconds.") + + sys.exit(bool(failed)) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/test/nhlbi_integration_tests/submit_test.py b/test/nhlbi_integration_tests/submit_test.py new file mode 100644 index 0000000..4ba4e53 --- /dev/null +++ b/test/nhlbi_integration_tests/submit_test.py @@ -0,0 +1,295 @@ +#!/usr/bin/env python3 +"""Register a new NHLBI integration test case. + +Uploads data files to private Azure Blob Storage, computes SHA256 checksums, +creates a .cfg test case file, and updates the data manifest. + +Usage: + python submit_test.py \\ + --name imoco_vds \\ + --config imoco_recon_vds.xml \\ + --noise-file /path/to/noise.h5 \\ + --data-file /path/to/recon_data.h5 \\ + --noise-config default_measurement_dependencies.xml \\ + --description "iMOCO VDS 3D lung reconstruction" \\ + --gpu-memory 8192 +""" + +import argparse +import configparser +import os +import sys +from pathlib import Path +import os.path as op + +import h5py +from test_utils import get_gadgetron_config_path +from get_nhlbi_data import ( + calc_sha256, + get_container_client, + load_manifest, + save_manifest, + upload_blob, +) + +CASES_DIR = Path(__file__).parent / "cases" + + +def validate_hdf5(filepath): + try: + with h5py.File(filepath, 'r') as f: + pass + return True + except Exception as e: + print(f"Error: Cannot read HDF5 file {filepath}: {e}") + return False + + +def validate_config_exists(config_name): + """Check if the XML config exists in common gadgetron config locations.""" + search_paths = [ + Path(get_gadgetron_config_path()) / config_name, + Path("config") / config_name, + Path("config/config") / config_name, + ] + for p in search_paths: + if p.exists(): + return True + print(f"Warning: Config '{config_name}' not found in standard locations. " + f"Ensure it is installed before running the test.") + return True # Warning only, don't block submission + + +def create_cfg(name, config, noise_config, description, gpu_memory, system_memory, + value_threshold, scale_threshold, tags,optional_additional_datasets=[],optional_additional_dependency_datasets=[]): + cfg = configparser.ConfigParser() + + cfg['dependency.siemens'] = { + 'data_file': f'{name}/noise_data.h5', + 'measurement': '0', + 'additional_arguments': 'skip_converstion', + } + cfg['dependency.client'] = { + 'configuration': noise_config, + } + for i in range (len(optional_additional_dependency_datasets)//2): + cfg[f'dependency.siemens.{i+1}'] = { + 'data_file': f'{name}/{op.basename(optional_additional_dependency_datasets[2*i])}', + 'measurement': '0', + 'additional_arguments': 'skip_converstion', + } + cfg[f'dependency.client.{i+1}'] = { + 'configuration': optional_additional_dependency_datasets[2*i+1], + } + + cfg['reconstruction.siemens'] = { + 'data_file': f'{name}/recon_data.h5', + 'measurement': '0', + 'additional_arguments': 'skip_converstion', + } + cfg['reconstruction.client'] = { + 'configuration': config, + } + cfg['reconstruction.test'] = { + 'reference_file': f'{name}/baseline_output.h5', + 'reference_images': f'{config}/image_0', + 'output_images': f'{config}/image_0', + 'value_comparison_threshold': str(value_threshold), + 'scale_comparison_threshold': str(scale_threshold), + } + cfg['requirements'] = { + 'system_memory': str(system_memory), + 'gpu_support': '1', + 'gpu_memory': str(gpu_memory), + } + + tag_list = ['nhlbi'] + [t.strip() for t in tags.split(',') if t.strip()] + cfg['tags'] = { + 'tags': ','.join(tag_list), + } + cfg['nhlbi'] = { + 'description': description, + 'noise_file': f'{name}/noise_data.h5', + } + + if optional_additional_datasets: + for i, dataset in enumerate(optional_additional_datasets): + cfg['nhlbi'].update({ + f'additional_dataset_{i}': f"{name}/{op.basename(dataset)}"}) + + cfg_path = CASES_DIR / f"{name}.cfg" + with open(cfg_path, 'w') as f: + cfg.write(f) + + return cfg_path + + +def main(): + parser = argparse.ArgumentParser( + description="Register a new NHLBI integration test", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument('--name', required=True, help="Test case name (e.g., imoco_vds)") + parser.add_argument('--config', required=True, help="Gadgetron XML config for reconstruction") + parser.add_argument('--noise-file', required=True, help="Path to noise calibration HDF5 file") + parser.add_argument('--data-file', required=True, help="Path to reconstruction input HDF5 file") + parser.add_argument('--noise-config', default='default_measurement_dependencies.xml', + help="Gadgetron XML config for noise dependency") + parser.add_argument('--description', default='', help="Human-readable test description") + parser.add_argument('--gpu-memory', type=int, default=8192, help="Required GPU memory in MB") + parser.add_argument('--system-memory', type=int, default=8192, help="Required system memory in MB") + parser.add_argument('--value-threshold', type=float, default=0.01, + help="Value comparison threshold for baseline validation") + parser.add_argument('--scale-threshold', type=float, default=0.01, + help="Scale comparison threshold for baseline validation") + parser.add_argument('--tags', type=str, default='', + help="Comma-separated additional tags (nhlbi is always included)") + parser.add_argument('--skip-upload', action='store_true', + help="Skip uploading to Azure (for local-only testing)") + parser.add_argument('--additional-files',nargs='+',type=str,default=[], + help='List of additional files for testing (e.g traj_bSTAR.seq traj_bSTAR.h5)') + parser.add_argument('--additional-dependencies',nargs='+',type=str,default=[], + help='List of additional dependency data.h5 .xml data2 .xml2') + + args = parser.parse_args() + print(args) + print(type(args.additional_files)) + print(len(args.additional_files)) + print(type(args.additional_dependencies)) + print(len(args.additional_dependencies)) + print(args.additional_dependencies) + # Check if test already exists + cfg_path = CASES_DIR / f"{args.name}.cfg" + if cfg_path.exists(): + print(f"Error: Test '{args.name}' already exists at {cfg_path}") + print("Use update_test.py to modify existing tests.") + sys.exit(1) + + + # Validate input files + print("Validating input files...") + if not op.isfile(args.noise_file): + print(f"Error: Noise file not found: {args.noise_file}") + sys.exit(1) + if not op.isfile(args.data_file): + print(f"Error: Data file not found: {args.data_file}") + sys.exit(1) + if not validate_hdf5(args.noise_file): + sys.exit(1) + if not validate_hdf5(args.data_file): + sys.exit(1) + + validate_config_exists(args.config) + + # Validate additional files + for additional_file in args.additional_files: + if not op.isfile(additional_file): + print(f"Error: Additional file not found: {additional_file}") + sys.exit(1) + if additional_file.endswith('.h5'): + if not validate_hdf5(additional_file): + sys.exit(1) + + # Additional dependencies + if len(args.additional_dependencies) % 2 !=0 : + print(f"Error additional dependencies required data.h5 and config.xml") + sys.exit(1) + for i in range(len(args.additional_dependencies)//2): + dependency_file=args.additional_dependencies[2*i] + print(dependency_file) + if not op.isfile(dependency_file): + print(f"Error: Dependency file not found: {dependency_file}") + sys.exit(1) + if dependency_file.endswith('.h5'): + if not validate_hdf5(dependency_file): + sys.exit(1) + + # Compute checksums + print("Computing checksums...") + noise_sha256 = calc_sha256(args.noise_file) + data_sha256 = calc_sha256(args.data_file) + print(f" Noise SHA256: {noise_sha256}") + print(f" Data SHA256: {data_sha256}") + + additional_files_sha256 = [] + for additional_file in args.additional_files: + sha256 = calc_sha256(additional_file) + additional_files_sha256.append((additional_file, sha256)) + print(f" Additional file {additional_file} SHA256: {sha256}") + + additional_dependency_file_sha256 = [] + for i in range(len(args.additional_dependencies)//2): + dependency_file=args.additional_dependencies[2*i] + sha256 = calc_sha256(dependency_file) + additional_dependency_file_sha256.append((dependency_file, sha256)) + print(f" Dependency file {dependency_file} SHA256: {sha256}") + + # Upload to Azure + if not args.skip_upload: + print("Uploading to Azure Blob Storage...") + container_client = get_container_client() + upload_blob(container_client, args.noise_file, f"{args.name}/noise_data.h5") + upload_blob(container_client, args.data_file, f"{args.name}/recon_data.h5") + for additional_file, sha256 in additional_files_sha256: + upload_blob(container_client, additional_file, f"{args.name}/{op.basename(additional_file)}") + for additional_dependency_file, sha256 in additional_dependency_file_sha256: + upload_blob(container_client, additional_dependency_file, f"{args.name}/{op.basename(additional_dependency_file)}") + else: + print("Skipping Azure upload (--skip-upload)") + + # Update manifest + manifest = load_manifest() + manifest.append({ + 'file': f'{args.name}/noise_data.h5', + 'sha256': noise_sha256, + 'type': 'noise', + 'test': args.name, + }) + manifest.append({ + 'file': f'{args.name}/recon_data.h5', + 'sha256': data_sha256, + 'type': 'input', + 'test': args.name, + }) + + for additional_file, sha256 in additional_files_sha256: + manifest.append({ + 'file': f"{args.name}/{op.basename(additional_file)}", + 'sha256': sha256, + 'type': 'additional', + 'test': args.name, + }) + + for additional_dependency_file, sha256 in additional_dependency_file_sha256: + manifest.append({ + 'file': f"{args.name}/{op.basename(additional_dependency_file)}", + 'sha256': sha256, + 'type': 'additional', + 'test': args.name, + }) + + + save_manifest(manifest) + print(f"Updated manifest: {len(manifest)} entries") + + # Create .cfg file + CASES_DIR.mkdir(parents=True, exist_ok=True) + cfg_path = create_cfg( + name=args.name, + config=args.config, + noise_config=args.noise_config, + description=args.description, + gpu_memory=args.gpu_memory, + system_memory=args.system_memory, + value_threshold=args.value_threshold, + scale_threshold=args.scale_threshold, + tags=args.tags, + optional_additional_datasets=args.additional_files, + optional_additional_dependency_datasets=args.additional_dependencies + ) + print(f"Created test case: {cfg_path}") + print(f"\nNext step: Run 'python generate_baseline.py --test {args.name}' to create baseline") + + +if __name__ == '__main__': + main() diff --git a/test/nhlbi_integration_tests/test_utils.py b/test/nhlbi_integration_tests/test_utils.py new file mode 100644 index 0000000..898396c --- /dev/null +++ b/test/nhlbi_integration_tests/test_utils.py @@ -0,0 +1,204 @@ +import json +from pathlib import Path +import os +import re +from pathlib import Path +import os.path as op +import ismrmrd +from typing import List, Union, Tuple +import numpy as np + + +def sort_by_indexes(lst:List, indexes:Union[List[str],List[int]], reverse:bool=False) -> List: + """ + Sort a list based on a list of indexes + + Parameters + ---------- + + lst : List, + List + + indexes : Union[List[str],List[int]], + List of index + + reverse : bool (optional, default : False), + flag to reverse the sorting + + Returns + ------- + + sorted_lst : List, + Sorted List + + """ + return [val for (_, val) in sorted(zip(indexes, lst), key=lambda x: x[0], reverse=reverse)] + +def read_images_h5(filename:str)->List[np.ndarray]: + """ + Getting the images from ISMRMD file + + Parameters + ---------- + + filename : str, + ISMRMD file path + + Returns + ------- + + list_img : List[np.ndarray], + List of images + + """ + list_img=[] + with ismrmrd.File(filename,'r') as mrd: + for key_img in list(mrd.find_images()): + img=np.array(mrd[key_img].images.data) #np.array(mrd[key_img].images.data).T + #Complex dtype: dtype([('real', 'List[dict]: + """ + Getting the headers from ISMRMD file + + Parameters + ---------- + + filename : str, + ISMRMD file path + + Returns + ------- + + list_headers : List[dict], + List of ISMRMRD headers transformed in dictionnary + + """ + list_headers=[] + with ismrmrd.File(filename,'r') as mrd: + for key_img in list(mrd.find_images()): + list_headers.append([dict(zip(mrd[key_img].images.headers[i].dtype.names,mrd[key_img].images.headers[i])) for i in range(mrd[key_img].images.headers.shape[0])]) + return list_headers + +def read_h5(filename:str, ordered:str='image_series_index')-> Tuple[np.ndarray,dict]: + """ + Getting the headers from ISMRMD file + + Parameters + ---------- + + filename : str, + ISMRMD file path + + ordered : str (optional, default :'image_series_index'), #'image_index' + string used for ordering the images and headers + + Returns + ------- + + images : List[np.ndarray], + List of images + + headers : List[dict], + List of ISMRMRD headers transformed in dictionnary + + """ + headers = read_headers_h5(filename) + images = read_images_h5(filename) + if not(ordered == ""): + indexes=[header[0][ordered] for header in headers] + images=sort_by_indexes(images,indexes) + headers=sort_by_indexes(headers,indexes) + return images,headers + +def resolve_env_path(value:str): + """ + value: str + A string that may contain an environment variable reference in the format ${env:VAR_NAME}. + """ + match = re.fullmatch(r"\$\{env:([A-Za-z_][A-Za-z0-9_]*)\}", value) + if match: + var_name = match.group(1) + env_value = os.environ.get(var_name) + if env_value: + return Path(env_value) + return None + +def get_cmake_install_prefix(settings_path:str="/opt/code/gadgetron/.vscode/settings.json"): + """ + settings_path: str + Path to the settings.json file. Default is .vscode/settings.json + Returns: + str: The CMAKE_INSTALL_PREFIX value from settings.json or a default path if not found. + Raises: + ValueError: If the path specified in CMAKE_INSTALL_PREFIX does not exist. + Notes: + - The function checks for both "CMAKE_INSTALL_PREFIX" and "cmake.configureSettings.CMAKE_INSTALL_PREFIX" keys in the settings.json file. + - If neither key is found, it defaults to "/opt/package/". + """ + # Check well-known install prefixes first (RT container, dev container) + well_known_prefixes = [ + os.environ.get("GADGETRON_HOME"), + "/opt/conda/envs/gadgetron", + "/opt/package", + ] + for prefix in well_known_prefixes: + if prefix and op.isdir(op.join(prefix, "bin")): + return prefix + + settings_file = Path(settings_path) + cmake_settings ="" + if settings_file.exists(): + with open(settings_file) as f: + raw = f.read() + cleaned_str = re.sub(r'//.*', '', raw) + # Remove /* ... */ comments + cleaned_str = re.sub(r'/\*.*?\*/', '', cleaned_str, flags=re.DOTALL) + + cleaned_str=re.sub(r',\s*([\]}])', r'\1', cleaned_str) + settings = json.loads(cleaned_str) + if "CMAKE_INSTALL_PREFIX" in settings: + cmake_settings= settings["CMAKE_INSTALL_PREFIX"] + elif "cmake.configureSettings" in settings: + cmake_settings = settings["cmake.configureSettings"].get("CMAKE_INSTALL_PREFIX") + else: + cmake_settings="/opt/package/" + else: + cmake_settings="/opt/package/" + if not cmake_settings or not op.exists(cmake_settings): + cmake_settings = resolve_env_path(cmake_settings) if cmake_settings else None + if not cmake_settings or not op.exists(cmake_settings): + raise ValueError(f"Path {cmake_settings} does not exist") + return cmake_settings + +def get_gadgetron_bin_path(): + """ + Retrieves the path to the Gadgetron binary from the CMAKE_INSTALL_PREFIX setting. + Returns: + str: The path to the Gadgetron binary. + Raises: + ValueError: If the CMAKE_INSTALL_PREFIX path does not exist. + """ + cmake_install_prefix = get_cmake_install_prefix() + gadgetron_bin_path = op.join(cmake_install_prefix, "bin") + if not op.exists(gadgetron_bin_path): + raise ValueError(f"Gadgetron binary path {gadgetron_bin_path} does not exist.") + return gadgetron_bin_path + +def get_gadgetron_config_path(): + """ + Retrieves the path to the Gadgetron config from the CMAKE_INSTALL_PREFIX setting. + Returns: + str: The path to the Gadgetron config. + Raises: + ValueError: If the CMAKE_INSTALL_PREFIX path does not exist. + """ + cmake_install_prefix = get_cmake_install_prefix() + gadgetron_config_path = op.join(cmake_install_prefix, "share","gadgetron","config") + if not op.exists(gadgetron_config_path): + raise ValueError(f"Gadgetron config path {gadgetron_config_path} does not exist.") + return gadgetron_config_path \ No newline at end of file diff --git a/test/nhlbi_integration_tests/update_test.py b/test/nhlbi_integration_tests/update_test.py new file mode 100644 index 0000000..e7ae1c7 --- /dev/null +++ b/test/nhlbi_integration_tests/update_test.py @@ -0,0 +1,185 @@ +#!/usr/bin/env python3 +"""Update an existing NHLBI integration test case. + +Supports updating data files, config, thresholds, and GPU/memory requirements. + +Usage: + python update_test.py imoco_vds --data-file /path/to/new_recon_data.h5 + python update_test.py imoco_vds --noise-file /path/to/new_noise.h5 + python update_test.py imoco_vds --config new_imoco_config.xml + python update_test.py imoco_vds --value-threshold 0.05 --scale-threshold 0.05 + python update_test.py imoco_vds --gpu-memory 16384 + python update_test.py imoco_vds --regenerate-baseline +""" + +import argparse +import configparser +import os +import sys +from pathlib import Path + +import h5py + +from get_nhlbi_data import ( + calc_sha256, + get_container_client, + load_manifest, + save_manifest, + upload_blob, +) + +CASES_DIR = Path(__file__).parent / "cases" + + +def main(): + parser = argparse.ArgumentParser( + description="Update an existing NHLBI integration test", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument('name', help="Test case name to update") + parser.add_argument('--data-file', help="New reconstruction input HDF5 file") + parser.add_argument('--noise-file', help="New noise calibration HDF5 file") + parser.add_argument('--config', help="New Gadgetron XML config for reconstruction") + parser.add_argument('--noise-config', help="New Gadgetron XML config for noise dependency") + parser.add_argument('--description', help="Updated test description") + parser.add_argument('--value-threshold', type=float, help="New value comparison threshold") + parser.add_argument('--scale-threshold', type=float, help="New scale comparison threshold") + parser.add_argument('--gpu-memory', type=int, help="New GPU memory requirement in MB") + parser.add_argument('--system-memory', type=int, help="New system memory requirement in MB") + parser.add_argument('--tags', help="New comma-separated tags (nhlbi is always included)") + parser.add_argument('--regenerate-baseline', action='store_true', + help="Invalidate current baseline and prompt for regeneration") + parser.add_argument('--skip-upload', action='store_true', + help="Skip uploading new files to Azure") + + args = parser.parse_args() + name = args.name + + # Validate test exists + cfg_path = CASES_DIR / f"{name}.cfg" + if not cfg_path.exists(): + print(f"Error: Test '{name}' not found at {cfg_path}") + sys.exit(1) + + config = configparser.ConfigParser() + config.read(cfg_path) + manifest = load_manifest() + changes = [] + + # Update data file + if args.data_file: + if not os.path.isfile(args.data_file): + print(f"Error: Data file not found: {args.data_file}") + sys.exit(1) + try: + h5py.File(args.data_file, 'r').close() + except Exception as e: + print(f"Error: Cannot read HDF5 file: {e}") + sys.exit(1) + + sha256 = calc_sha256(args.data_file) + if not args.skip_upload: + container_client = get_container_client() + upload_blob(container_client, args.data_file, f"{name}/recon_data.h5") + + # Update manifest + for entry in manifest: + if entry.get('test') == name and entry.get('type') == 'input': + entry['sha256'] = sha256 + break + else: + manifest.append({'file': f'{name}/recon_data.h5', 'sha256': sha256, 'type': 'input', 'test': name}) + + changes.append("Updated reconstruction data file") + + # Update noise file + if args.noise_file: + if not os.path.isfile(args.noise_file): + print(f"Error: Noise file not found: {args.noise_file}") + sys.exit(1) + try: + h5py.File(args.noise_file, 'r').close() + except Exception as e: + print(f"Error: Cannot read HDF5 file: {e}") + sys.exit(1) + + sha256 = calc_sha256(args.noise_file) + if not args.skip_upload: + container_client = get_container_client() + upload_blob(container_client, args.noise_file, f"{name}/noise_data.h5") + + for entry in manifest: + if entry.get('test') == name and entry.get('type') == 'noise': + entry['sha256'] = sha256 + break + else: + manifest.append({'file': f'{name}/noise_data.h5', 'sha256': sha256, 'type': 'noise', 'test': name}) + + changes.append("Updated noise data file") + + # Update reconstruction config + if args.config: + config['reconstruction.client']['configuration'] = args.config + config['reconstruction.test']['reference_images'] = f'{args.config}/image_0' + config['reconstruction.test']['output_images'] = f'{args.config}/image_0' + changes.append(f"Updated reconstruction config to {args.config}") + + # Update noise config + if args.noise_config: + config['dependency.client']['configuration'] = args.noise_config + changes.append(f"Updated noise config to {args.noise_config}") + + # Update thresholds + if args.value_threshold is not None: + config['reconstruction.test']['value_comparison_threshold'] = str(args.value_threshold) + changes.append(f"Updated value threshold to {args.value_threshold}") + + if args.scale_threshold is not None: + config['reconstruction.test']['scale_comparison_threshold'] = str(args.scale_threshold) + changes.append(f"Updated scale threshold to {args.scale_threshold}") + + # Update requirements + if args.gpu_memory is not None: + config['requirements']['gpu_memory'] = str(args.gpu_memory) + changes.append(f"Updated GPU memory requirement to {args.gpu_memory} MB") + + if args.system_memory is not None: + config['requirements']['system_memory'] = str(args.system_memory) + changes.append(f"Updated system memory requirement to {args.system_memory} MB") + + # Update tags + if args.tags is not None: + tag_list = ['nhlbi'] + [t.strip() for t in args.tags.split(',') if t.strip()] + config['tags']['tags'] = ','.join(tag_list) + changes.append(f"Updated tags to {','.join(tag_list)}") + + # Update description + if args.description is not None: + config['nhlbi']['description'] = args.description + changes.append("Updated description") + + # Regenerate baseline + if args.regenerate_baseline: + # Remove baseline entry from manifest + manifest = [e for e in manifest if not (e.get('test') == name and e.get('type') == 'baseline')] + changes.append("Invalidated baseline") + + if not changes: + print("No changes specified. Use --help for options.") + sys.exit(0) + + # Write changes + with open(cfg_path, 'w') as f: + config.write(f) + save_manifest(manifest) + + print(f"Updated test '{name}':") + for change in changes: + print(f" - {change}") + + if args.regenerate_baseline or args.data_file or args.noise_file or args.config: + print(f"\nReminder: Run 'python generate_baseline.py --test {name}' to regenerate baseline") + + +if __name__ == '__main__': + main() diff --git a/toolboxes/core/gpu/cuSparseMatrix.h b/toolboxes/core/gpu/cuSparseMatrix.h index 1b1a1f5..081808b 100644 --- a/toolboxes/core/gpu/cuSparseMatrix.h +++ b/toolboxes/core/gpu/cuSparseMatrix.h @@ -39,18 +39,26 @@ namespace Gadgetron cuCsrMatrix &operator=(cuCsrMatrix &&other) { + if (this == &other) + return *this; + if (this->descr) + cusparseDestroySpMat(this->descr); this->descr = other.descr; other.descr = nullptr; + this->rows = other.rows; + this->cols = other.cols; + other.rows = 0; + other.cols = 0; this->csrColdnd = std::move(other.csrColdnd); this->csrRow = std::move(other.csrRow); - this->data = std::move(this->data); + this->data = std::move(other.data); return *this; } - size_t rows, cols; + size_t rows = 0, cols = 0; thrust::device_vector csrRow, csrColdnd; thrust::device_vector data; - cusparseSpMatDescr_t descr; + cusparseSpMatDescr_t descr = nullptr; }; /** diff --git a/toolboxes/nfft/gpu/ConvolverNC2C_sparse.cuh b/toolboxes/nfft/gpu/ConvolverNC2C_sparse.cuh index 132be11..066bb75 100644 --- a/toolboxes/nfft/gpu/ConvolverNC2C_sparse.cuh +++ b/toolboxes/nfft/gpu/ConvolverNC2C_sparse.cuh @@ -149,19 +149,19 @@ void check_csrMatrix(cuCsrMatrix &matrix) template class K> cuCsrMatrix make_conv_matrix( - const thrust::device_vector,D>> &points, - const vector_td& image_dims, - const ConvolutionKernel, D, K>* kernel) + const thrust::device_vector,D>> &points, + const vector_td& image_dims, + const ConvolutionKernel, D, K>* d_kernel, + realType_t radius) { auto csrRow = thrust::device_vector(points.size()+1); csrRow[0] = 0; CHECK_FOR_CUDA_ERROR(); - realType_t radius = kernel->get_radius(); { thrust::device_vector c_p_s(points.size()); thrust::transform(points.begin(), points.end(), c_p_s.begin(), - compute_num_cells_per_sample,D>(kernel->get_radius())); + compute_num_cells_per_sample,D>(radius)); thrust::inclusive_scan( c_p_s.begin(), c_p_s.end(), csrRow.begin()+1, thrust::plus()); // prefix sum @@ -183,11 +183,11 @@ cuCsrMatrix make_conv_matrix( thrust::raw_pointer_cast(csrRow.data()), thrust::raw_pointer_cast(data.data()), thrust::raw_pointer_cast(csrColdnd.data()), - vector_td(image_dims), points.size(),kernel); + vector_td(image_dims), points.size(),d_kernel); cudaDeviceSynchronize(); CHECK_FOR_CUDA_ERROR(); - cuCsrMatrix matrix(prod(image_dims), points.size(),std::move(csrRow),std::move(csrColdnd),std::move(data)); + cuCsrMatrix matrix(points.size(), prod(image_dims),std::move(csrRow),std::move(csrColdnd),std::move(data)); return matrix; } diff --git a/toolboxes/nfft/gpu/cuGriddingConvolution.cu b/toolboxes/nfft/gpu/cuGriddingConvolution.cu index 2253c43..7a9ba25 100644 --- a/toolboxes/nfft/gpu/cuGriddingConvolution.cu +++ b/toolboxes/nfft/gpu/cuGriddingConvolution.cu @@ -15,7 +15,18 @@ #include "ConvolverNC2C_standard.cuh" #define CUDA_CONV_MAX_COILS (16) -#define CUDA_CONV_THREADS_PER_KERNEL (512) // Optimized for Blackwell/Hopper (testing 512) +#if defined(__CUDA_ARCH__) + #if __CUDA_ARCH__ >= 1200 + // Optimized for Blackwell (compute capability 9.0+) + #define CUDA_CONV_THREADS_PER_KERNEL 512 + #else + // Default for other architectures + #define CUDA_CONV_THREADS_PER_KERNEL 192 + #endif +#else + // Host code or unknown arch + #define CUDA_CONV_THREADS_PER_KERNEL 192 +#endif namespace Gadgetron { @@ -316,7 +327,7 @@ namespace Gadgetron auto view_dims = to_std_vector(this->plan_.matrix_size_os_); view_dims.push_back(this->plan_.num_frames_); view_dims.push_back(0); // Placeholder for num_coils. - + //GDEBUG_STREAM("domain_size_coils_desired = " << domain_size_coils_desired << " num_repetitions = " << num_repetitions); for (unsigned int repetition = 0; repetition < num_repetitions; repetition++) { // Number of coils in this repetition. @@ -372,7 +383,7 @@ namespace Gadgetron REAL radius = this->plan_.kernel_.get_radius(); transform(trajectory.begin(), trajectory.end(), c_p_s.begin(), compute_num_cells_per_sample(radius)); - inclusive_scan(c_p_s.begin(), c_p_s.end(), c_p_s_ps.begin(), + thrust::inclusive_scan(c_p_s.begin(), c_p_s.end(), c_p_s_ps.begin(), thrust::plus()); // Prefix sum. // Build the vector of (grid_idx, sample_idx) tuples. Actually kept in @@ -668,7 +679,8 @@ namespace Gadgetron { this->conv_matrix_ = std::make_unique>( make_conv_matrix( - trajectory, this->plan_.matrix_size_os_, this->plan_.d_kernel_)); + trajectory, this->plan_.matrix_size_os_, this->plan_.d_kernel_, + this->plan_.kernel_.get_radius())); } diff --git a/toolboxes/nhlbi_gt_toolbox/CMakeLists.txt b/toolboxes/nhlbi_gt_toolbox/CMakeLists.txt index 4c5ba7e..11d61d6 100644 --- a/toolboxes/nhlbi_gt_toolbox/CMakeLists.txt +++ b/toolboxes/nhlbi_gt_toolbox/CMakeLists.txt @@ -66,6 +66,11 @@ set(gadgetron_nhlbi_gt_toolbox_config_files config/pulmonary_MOCOLR.xml config/cardiovascular_iMOCO.xml config/cardiopulmonary_recon.xml + config/spiral_3DRT.xml + config/spiral_3DRT_CSM.xml + config/spiral_3DRT_CSM_test.xml + config/spiral_3DRT_test.xml + config/imoco_recon_vds.xml ) diff --git a/toolboxes/nhlbi_gt_toolbox/config/imoco_recon_vds.xml b/toolboxes/nhlbi_gt_toolbox/config/imoco_recon_vds.xml new file mode 100644 index 0000000..1719a5a --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/config/imoco_recon_vds.xml @@ -0,0 +1,160 @@ + + + 2 + + + + gadgetron_core_readers + AcquisitionReader + + + gadgetron_core_readers + WaveformReader + + + + + + gadgetron_core_writers + ImageWriter + + + + + + + + WaveformToTrajectory + nhlbi_gt_gadgets + WaveformToTrajectory + perform_GIRFtrue + GIRF_folder/opt/GIRF/ + generateTrajtrue + attachWaveformfalse + + + + + gadgetron_mricore + NoiseAdjustGadget + + + + RemoveSpiralOversampling + nhlbi_gt_gadgets + RemoveSpiralOversampling + + + + RemoveNavsGadget + nhlbi_gt_gadgets + RemoveNavsGadget + + + + + + + + + + + + + + + + ImagetoVector + nhlbi_gt_gadgets + ImagetoVector + + + + PrepreconParams + nhlbi_gt_gadgets + PrepreconParams + matOSP_vector1.5 1.5 1.5 + downsampling_vector1 1 1 + warpCUDA_vectortrue true true + is3Dtrue + kernel_width3 + oversampling_factor1.5 + kernel_width_dcf3 + iterations_dcf10 + oversampling_factor_dcf2.1 + useIterativeDCWEstimatedfalse + lambda_spatial0.001 + lambda_spatial_imoco0.001 + lambda_time0.1 + lambda_time20.0 + iterations5 + iterations_imoco5 + iterations_inner2 + tolerance100 + norm2 + use_gccfalse + gcc_coils6 + doMC_iterfalse + iteration_count_moco3 + + + + + Noncart_recon_gadget + nhlbi_gt_gadgets + Noncart_recon_gadget + Debug0 + doConcomitantFieldCorrectiontrue + referencePhase0 + estimateCSM_perc50 + reconType0 + processingType0 + + + + + ImageArraySplit + gadgetron_mricore + ImageArraySplitGadget + + + + ComplexToFloatAttrib + gadgetron_mricore + ComplexToFloatGadget + + + + AutoScaleGadget + gadgetron_mricore + AutoScaleGadget + + + + FloatToShort + gadgetron_mricore + FloatToUShortGadget + + + + ImageFinish + gadgetron_mricore + ImageFinishGadget + + + + \ No newline at end of file diff --git a/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT.xml b/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT.xml new file mode 100644 index 0000000..43471d1 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT.xml @@ -0,0 +1,122 @@ + + + 2 + + + + gadgetron_core_readers + AcquisitionReader + + + gadgetron_core_readers + WaveformReader + + + + + + gadgetron_core_writers + ImageWriter + + + + + + + + NoiseAdjustgadgetron_mricoreNoiseAdjustGadget + + + + WaveformToTrajectory + + nhlbi_gt_gadgets + WaveformToTrajectory + perform_GIRFtrue + GIRF_folder/opt/GIRF/GIRF_20250225/ + attachWaveformfalse + generateTrajtrue + generateTrajtrue + + + + RemoveNavsGadget + + nhlbi_gt_gadgets + RemoveNavsGadget + + + + RemoveSpiralOversampling + + nhlbi_gt_gadgets + RemoveSpiralOversampling + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ImageArraySplit + gadgetron_mricore + ImageArraySplitGadget + + + + ComplexToFloatAttrib + gadgetron_mricore + ComplexToFloatGadget + + + + + AutoScaleGadget + gadgetron_mricore + AutoScaleGadget + + + + FloatToShort + gadgetron_mricore + FloatToUShortGadget + + + + + ImageFinish + gadgetron_mricore + ImageFinishGadget + + + + + + diff --git a/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_CSM.xml b/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_CSM.xml new file mode 100644 index 0000000..2158b36 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_CSM.xml @@ -0,0 +1,135 @@ + + + 2 + + + + gadgetron_core_readers + AcquisitionReader + + + gadgetron_core_readers + WaveformReader + + + + + + gadgetron_core_writers + ImageWriter + + + + + + NoiseAdjustgadgetron_mricoreNoiseAdjustGadget + + + + + WaveformToTrajectory + + nhlbi_gt_gadgets + WaveformToTrajectory + perform_GIRFtrue + GIRF_folder/opt/GIRF/GIRF_20250225/ + attachWaveformfalse + generateTrajtrue + + + + RemoveNavsGadget + + nhlbi_gt_gadgets + RemoveNavsGadget + + + + RemoveSpiralOversampling + + nhlbi_gt_gadgets + RemoveSpiralOversampling + + + + Spiral_Rovir + + nhlbi_gt_gadgets + Spiral_Rovir + combine_csm1 + send_csm0 + + + + + + + + + + + + + + + + + + + + + ImageArraySplit + gadgetron_mricore + ImageArraySplitGadget + + + + ComplexToFloatAttrib + gadgetron_mricore + ComplexToFloatGadget + + + + AutoScaleGadget + gadgetron_mricore + AutoScaleGadget + + + + + + + FloatToShort + gadgetron_mricore + FloatToUShortGadget + + + + ImageFinish + gadgetron_mricore + ImageFinishGadget + + + + + diff --git a/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_CSM_test.xml b/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_CSM_test.xml new file mode 100644 index 0000000..9cbeab6 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_CSM_test.xml @@ -0,0 +1,136 @@ + + + 2 + + + + gadgetron_core_readers + AcquisitionReader + + + gadgetron_core_readers + WaveformReader + + + + + + gadgetron_core_writers + ImageWriter + + + + + + NoiseAdjustgadgetron_mricoreNoiseAdjustGadget + + + + + WaveformToTrajectory + + nhlbi_gt_gadgets + WaveformToTrajectory + perform_GIRFtrue + GIRF_folder/opt/GIRF/GIRF_20250225/ + attachWaveformfalse + generateTrajtrue + + + + RemoveNavsGadget + + nhlbi_gt_gadgets + RemoveNavsGadget + + + + RemoveSpiralOversampling + + nhlbi_gt_gadgets + RemoveSpiralOversampling + + + + Spiral_Rovir + + nhlbi_gt_gadgets + Spiral_Rovir + combine_csm1 + send_csm0 + + + + + + + + + + + + + + + + + + + + + + ImageArraySplit + gadgetron_mricore + ImageArraySplitGadget + + + + ComplexToFloatAttrib + gadgetron_mricore + ComplexToFloatGadget + + + + AutoScaleGadget + gadgetron_mricore + AutoScaleGadget + + + + + + + FloatToShort + gadgetron_mricore + FloatToUShortGadget + + + + ImageFinish + gadgetron_mricore + ImageFinishGadget + + + + + diff --git a/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_test.xml b/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_test.xml new file mode 100644 index 0000000..d870521 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/config/spiral_3DRT_test.xml @@ -0,0 +1,116 @@ + + + 2 + + + + gadgetron_core_readers + AcquisitionReader + + + gadgetron_core_readers + WaveformReader + + + + + + gadgetron_core_writers + ImageWriter + + + + + + + + NoiseAdjustgadgetron_mricoreNoiseAdjustGadget + + + + WaveformToTrajectory + + nhlbi_gt_gadgets + WaveformToTrajectory + perform_GIRFtrue + GIRF_folder/opt/GIRF/GIRF_20250225/ + attachWaveformfalse + generateTrajtrue + generateTrajtrue + + + + RemoveNavsGadget + + nhlbi_gt_gadgets + RemoveNavsGadget + + + + RemoveSpiralOversampling + + nhlbi_gt_gadgets + RemoveSpiralOversampling + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ImageArraySplit + gadgetron_mricore + ImageArraySplitGadget + + + + ComplexToFloatAttrib + gadgetron_mricore + ComplexToFloatGadget + + + + + AutoScaleGadget + gadgetron_mricore + AutoScaleGadget + + + + FloatToShort + gadgetron_mricore + FloatToUShortGadget + + + + + ImageFinish + gadgetron_mricore + ImageFinishGadget + + + + + + diff --git a/toolboxes/nhlbi_gt_toolbox/doc/igt.jpg b/toolboxes/nhlbi_gt_toolbox/doc/igt.jpg new file mode 100644 index 0000000..6836cdc Binary files /dev/null and b/toolboxes/nhlbi_gt_toolbox/doc/igt.jpg differ diff --git a/toolboxes/nhlbi_gt_toolbox/doc/installation.rst b/toolboxes/nhlbi_gt_toolbox/doc/installation.rst index 51fd32d..913adb7 100644 --- a/toolboxes/nhlbi_gt_toolbox/doc/installation.rst +++ b/toolboxes/nhlbi_gt_toolbox/doc/installation.rst @@ -10,14 +10,14 @@ First of all, you will install Gadgetron : .. code-block:: console - git clone -b cardiopulmonary_bstar git@github.com:NHLBI/lit_gadgetron.git + git clone -b volumetric_rt_mri git@github.com:NHLBI/lit_gadgetron.git cd gadgetron conda env create -f environment.yml conda activate gadgetron mkdir build && cd build && cmake ../ -GNinja -DCMAKE_INSTALL_PREFIX=${CONDA_PREFIX} -DC_PREFIX_PATH=${CONDA_PREFIX} -DUSE_CUDA=ON -DUSE_MKL=ON ninja && ninja install -Once built, the package can be used with gadgetron using the config xml files provided with this repository (`config files repository `_). +Once built, the package can be used with gadgetron using the config xml files provided with this repository (`config files repository `_). Validate installation +++++++++++++++++++++ @@ -42,21 +42,117 @@ about your installed version of the Gadgetron and it would look something like t The output may vary on your specific setup, but you will see error messages if the Gadgetron is not installed or not installed correctly. +Validate image reconstruction pipelines ++++++++++++++++++++++++++++++++++++++++ +To validate that the Gadgetron is working correctly with the NHLBI toolbox, you can run the following command to test the prescan reference image reconstruction pipeline for example: + +.. code-block:: console + + conda activate gadgetron + cd test/nhlbi_integration_tests/ + python run_nhlbi_tests.py cases/csm_prescan.cfg -F + +The expected output of the test should look like this:: + + Downloading test data... + Downloading: csm_prescan/noise_data.h5 + Downloading: csm_prescan/recon_data.h5 + Downloading: csm_prescan/baseline_output.h5 + Saved: /opt/code/gadgetron/test/nhlbi_integration_tests/data/csm_prescan/baseline_output.h5 + Saved: /opt/code/gadgetron/test/nhlbi_integration_tests/data/csm_prescan/noise_data.h5 + Saved: /opt/code/gadgetron/test/nhlbi_integration_tests/data/csm_prescan/recon_data.h5 + Querying Gadgetron capabilities... + + Test 1 of 1: cases/csm_prescan.cfg + + Running Gadgetron test cases/csm_prescan.cfg with: + -- ISMRMRD_HOME : None + -- GADGETRON_HOME : None + -- TEST CASE : cases/csm_prescan.cfg + Starting MRD Storage Server on port 9113 + Starting Gadgetron instance on port 9003 + Copying prepared ISMRMRD data: /opt/code/gadgetron/test/nhlbi_integration_tests/data/csm_prescan/noise_data.h5 -> test/dependency.siemens.copied.mrd + Passing data to Gadgetron: test/dependency.siemens.copied.mrd -> test/dependency.client.output.mrd + Gadgetron processing time: 0.12 s + Copying prepared ISMRMRD data: /opt/code/gadgetron/test/nhlbi_integration_tests/data/csm_prescan/recon_data.h5 -> test/reconstruction.siemens.copied.mrd + Passing data to Gadgetron: test/reconstruction.siemens.copied.mrd -> test/reconstruction.client.output.mrd + Gadgetron processing time: 24.14 s + reconstruction.test [OK] (Norm: 9.2e-06 [0.01] Scale: 0.0e+00 [0.01]) + reconstruction.test [OK] (Output headers matched reference) + Test status: Passed + Speed improved: 24.3s vs baseline 56.2s (-56.8%) + + 1 tests passed. 0 tests failed. 0 tests skipped. 0 missing baselines. 0 speed regressions. + Total processing time: 24.26 seconds. + +Or the 3D real-time image reconstruction pipeline: + +.. code-block:: console + + cd test/nhlbi_integration_tests/ + python run_nhlbi_tests.py cases/3drt.cfg -F + +The expected output of the test should look like this:: + + Downloading test data... + Downloading: 3drt/noise_data.h5 + Downloading: 3drt/recon_data.h5 + Downloading: 3drt/prescan-reference_noise.h5 + Downloading: 3drt/prescan-reference.h5 + Saved: /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/prescan-reference_noise.h5 + Downloading: 3drt/baseline_output.h5 + Saved: /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/noise_data.h5 + Saved: /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/recon_data.h5 + Saved: /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/prescan-reference.h5 + Warning: Could not download data for test '3drt': Downloaded file /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/baseline_output.h5 failed validation. Expected SHA256 04e6d61654c5040e9589cb1ef02019e445990ef72c4381e151f2342e95ba4542. Actual SHA256 48011834c0279b9f4ef82a92b967015dff65ead9c3eba26b8ed081540a88b47c + Querying Gadgetron capabilities... + + Test 1 of 1: cases/3drt.cfg + + Running Gadgetron test cases/3drt.cfg with: + -- ISMRMRD_HOME : None + -- GADGETRON_HOME : None + -- TEST CASE : cases/3drt.cfg + Starting MRD Storage Server on port 9113 + Starting Gadgetron instance on port 9003 + Copying prepared ISMRMRD data: /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/noise_data.h5 -> test/dependency.siemens.copied.mrd + Passing data to Gadgetron: test/dependency.siemens.copied.mrd -> test/dependency.client.output.mrd + Gadgetron processing time: 0.12 s + Copying prepared ISMRMRD data: /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/prescan-reference_noise.h5 -> test/dependency.siemens.1.copied.mrd + Passing data to Gadgetron: test/dependency.siemens.1.copied.mrd -> test/dependency.client.1.output.mrd + Gadgetron processing time: 0.12 s + Copying prepared ISMRMRD data: /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/prescan-reference.h5 -> test/dependency.siemens.2.copied.mrd + Passing data to Gadgetron: test/dependency.siemens.2.copied.mrd -> test/dependency.client.2.output.mrd + Gadgetron processing time: 31.04 s + Copying prepared ISMRMRD data: /opt/code/gadgetron/test/nhlbi_integration_tests/data/3drt/recon_data.h5 -> test/reconstruction.siemens.copied.mrd + Passing data to Gadgetron: test/reconstruction.siemens.copied.mrd -> test/reconstruction.client.output.mrd + Gadgetron processing time: 41.94 s + reconstruction.test [OK] (Norm: 3.0e-03 [0.01] Scale: 7.3e-06 [0.01]) + reconstruction.test [OK] (Output headers matched reference) + Test status: Passed + SPEED REGRESSION: 73.2s vs baseline 42.3s (+73.1%, threshold 50%) + + Speed regressions: + cases/3drt.cfg + + 1 tests passed. 0 tests failed. 0 tests skipped. 0 missing baselines. 1 speed regressions. + Total processing time: 73.22 seconds. + Docker container ---------------- -Alternatively, you can test the code by pulling the provided docker image located in packages using the following command: +Alternatively, you can test the code by pulling the provided docker image located in `packages repository `_ using the following command: .. code-block:: console - docker pull ghcr.io/nhlbi/litgt_cardio_pulmonary_bstar_rt:20260205 + docker pull ghcr.io/nhlbi/litgt_volumetric_rt_mri_rt:20260828 This image can be deployed with: .. code-block:: console - docker run --gpus all --name=cardio_pulmonary_bstar_rt -ti -p 9063:9002 --volume=[LOCAL_DATA_FOLDER]:/opt/data --restart unless-stopped --detach ghcr.io/nhlbi/litgt_cardio_pulmonary_bstar_rt:20260205` + docker run --gpus all --name=volumetric_rt_mri -ti -p 9063:9002 -p 9072:9004 --volume=[LOCAL_DATA_FOLDER]:/opt/data --restart unless-stopped --detach ghcr.io/nhlbi/litgt_volumetric_rt_mri_rt:20260828` where **LOCAL_DATA_FOLDER** is the path to a folder containing raw data that can be used for testing the reconstruction. @@ -67,21 +163,22 @@ Once the docker container is running, you can start a bash terminal inside the c .. code-block:: console - docker exec -ti cardio_pulmonary_bstar_rt bash + docker exec -ti volumetric_rt_mri bash -and you can simply ou can simply navigate to `/opt/data/` and test the code : +and you can simply validate the image reconstruction pipeline using our integration tests (See precedent paragraph) or you can navigate to `/opt/data/` and test the code for the prescan reference image reconstruction pipeline for example: .. code-block:: console cd /opt/data - gadgetron_ismrmrd_client -p 9002 -f DATA_FILE -c XXX.xml -o OUTPUT_FILENAME.h5` + gadgetron_ismrmrd_client -p 9002 -f prescan-reference_noise.h5 -c default_measurement_dependencies.xml + gadgetron_ismrmrd_client -p 9002 -f prescan-reference.h5 -c spiral_3DRT_CSM_test.xml -o OUTPUT_FILENAME.h5 In another terminal session you can monitor the logs from the container .. code-block:: console - docker logs -f cardio_pulmonary_bstar_rt` + docker logs -f volumetric_rt_mri Please note that if you are using the gadgetron_ismrmrd_client from outside the container then you may need to specify the server address with **-a SERVER_ADDRESS** and the port **-p 9063** @@ -89,12 +186,22 @@ Please note that if you are using the gadgetron_ismrmrd_client from outside the .. code-block:: console cd LOCAL_DATA_FOLDER - gadgetron_ismrmrd_client -a SERVER_ADDRESS -p 9063 -f DATA_FILE -c XXX.xml -o OUTPUT_FILENAME.h5` + gadgetron_ismrmrd_client -a SERVER_ADDRESS -p 9063 -f prescan-reference_noise.h5 -c default_measurement_dependencies.xml + gadgetron_ismrmrd_client -a SERVER_ADDRESS -p 9063 -f prescan-reference.h5 -c spiral_3DRT_CSM_test.xml -o OUTPUT_FILENAME.h5 + + +3D Slicer Integration +--------------------- + +If the SlicerGadget is enabled, `3D Slicer `_ with the `SlicerIGT `_ extension can be used for live data streaming. In order to use, setup SlicerIGT on the host device with the corresponding forwarded port to the container (ex. 9072) in the OpenIGTLinkIF module, as shown in the picture below. + +.. image:: igt.jpg + :width: 50% Dataset ------- -The test data can be downloaded from zenodo: `18461603 `_ +The test data can be downloaded from zenodo: `19005977 `_ .. note:: More Information on Gadgetron are available over here : diff --git a/toolboxes/nhlbi_gt_toolbox/environment.yml b/toolboxes/nhlbi_gt_toolbox/environment.yml index 3014222..21db772 100644 --- a/toolboxes/nhlbi_gt_toolbox/environment.yml +++ b/toolboxes/nhlbi_gt_toolbox/environment.yml @@ -45,5 +45,6 @@ dependencies: - einops - wandb - onnxruntime-gpu + - pyigtl diff --git a/toolboxes/nhlbi_gt_toolbox/gadgets/CMakeLists.txt b/toolboxes/nhlbi_gt_toolbox/gadgets/CMakeLists.txt index 0be981d..b5495cc 100644 --- a/toolboxes/nhlbi_gt_toolbox/gadgets/CMakeLists.txt +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/CMakeLists.txt @@ -17,6 +17,8 @@ utility_gadgets/AcquisitionReconParamsFanout.cpp utility_gadgets/AcquisitionStopTimingGadget.cpp waveforms/AcquisitionWaveformFanout.cpp #waveforms/RemoveAcqsGadget.cpp +noncart_recon/Spiral_3DRT.cpp +noncart_recon/Spiral_Rovir.cpp ) add_library(nhlbi_gt_gadgets SHARED diff --git a/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Noncart_recon_gadget.cpp b/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Noncart_recon_gadget.cpp index ab53ac2..e072ff9 100644 --- a/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Noncart_recon_gadget.cpp +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Noncart_recon_gadget.cpp @@ -205,6 +205,8 @@ class Noncart_recon_gadget allAcq[idx] = std::move(Core::get(message)); if ((idx >= int((estimateCSM_perc / 100.0) * maxAcq)) && (!csm_calculated_ && recon_params_received)) { + auto& [headAcq_0, dataAcq_0, trajAcq_0] = allAcq[0]; + acqhdr = headAcq_0; GadgetronTimer timer_CSM("Calculating CSM"); GadgetronTimer timer_Average("Calculating Average Image"); cudaSetDevice(recon_params.selectedDevices[0]); @@ -225,6 +227,10 @@ class Noncart_recon_gadget } timer_CSM.stop(); csm_calculated_ = true; + if (save_csm){ + process_and_send_images(*csm, acqhdr, out, series_counter, "CSM", recon_params); + series_counter++; + } if (save_avg) { *channel_images *= *conj(csm.get()); auto combined = sum(channel_images.get(), channel_images->get_number_of_dimensions() - 1); @@ -251,6 +257,8 @@ class Noncart_recon_gadget if (!csm_calculated_){ + auto& [headAcq_1, dataAcq_1, trajAcq_1] = allAcq[0]; + acqhdr = headAcq_1; GadgetronTimer timer_CSM("Calculating CSM At the end"); GadgetronTimer timer_Average("Calculating Average Image"); cudaSetDevice(recon_params.selectedDevices[0]); @@ -271,6 +279,10 @@ class Noncart_recon_gadget } timer_CSM.stop(); csm_calculated_ = true; + if (save_csm){ + process_and_send_images(*csm, acqhdr, out, series_counter, "CSM", recon_params); + series_counter++; + } if (save_avg) { *channel_images *= *conj(csm.get()); auto combined = sum(channel_images.get(), channel_images->get_number_of_dimensions() - 1); @@ -465,6 +477,8 @@ class Noncart_recon_gadget } break; case 4: { + + GadgetronTimer timer_4D_respi("4D Respiratory Recon :"); std::vector binning_order_respi = {binning_order[1], binning_order[2], binning_order[0]}; auto output_collapsed = nhlbi_toolbox::utils::sort_idx_phases(idx_phases_vec, binning_order_respi, true,start_idx_nc); std::vector> idx_phases_respiratory = std::get<0>(output_collapsed); @@ -487,9 +501,13 @@ class Noncart_recon_gadget std::vector> trajVec_respi =reconstruction4D.arraytovector(&traj_respi, number_elements_respi); std::vector> dcwVec_respi = reconstruction4D.estimate_dcf(&trajVec_respi); auto ave_cuIimages = reconstruction4D.reconstruct(&cuData_respi, &trajVec_respi, &dcwVec_respi, csm,false); - process_and_send_images(ave_cuIimages, acqhdr, out, series_counter, + // Save respiratory-resolved images + if(save_intermediate_images){ + process_and_send_images(ave_cuIimages, acqhdr, out, series_counter, std::string("4DTresolved") + img_parameters_name, recon_params); - series_counter++; + series_counter++; + } + cuData_respi.clear(); trajVec_respi.clear(); dcwVec_respi.clear(); @@ -498,12 +516,13 @@ class Noncart_recon_gadget recon_params_adv.shots_per_time = shots_per_time; reconstruction5D.set_recon_params(recon_params_adv); reconstruction->set_recon_params(recon_params); - + /* Verbose for (size_t it = 0; it < shots_per_time.get_number_of_elements(); it++) { GDEBUG_STREAM("it " << it << "SHOTs " << *(shots_per_time.begin() + it)); size_t size_phase = (idx_phases[it].size()); GDEBUG_STREAM("it " << it << "Phase size " << size_phase); } + */ std::vector> trajVec_respi_cardiac =reconstruction5D.arraytovector(&traj_rc, number_elements_rc); std::vector> dcwVec_respi_cardiac =reconstruction5D.estimate_dcf(&trajVec_respi_cardiac); cuIimages = reconstruction5D.reconstructiMOCO_avg_image(&cuData_All, &trajVec_respi_cardiac, &dcwVec_respi_cardiac, ave_cuIimages, csm, referencePhase); @@ -573,8 +592,7 @@ class Noncart_recon_gadget std::vector> trajVec_respi =reconstruction4D.arraytovector(&traj_respi, number_elements_respi); std::vector> dcwVec_respi = reconstruction4D.estimate_dcf(&trajVec_respi); auto ave_cuIimages = reconstruction4D.reconstructMOCOLR(&cuData_respi, &trajVec_respi, &dcwVec_respi, csm); - process_and_send_images(ave_cuIimages, acqhdr, out, series_counter, - std::string("4DMOCOLR") + img_parameters_name, recon_params); + process_and_send_images(ave_cuIimages, acqhdr, out, series_counter,std::string("4DMOCOLR") + img_parameters_name, recon_params); series_counter++; cuData_respi.clear(); trajVec_respi.clear(); @@ -638,16 +656,25 @@ class Noncart_recon_gadget Gadgetron::reconParams& recon_params) { size_t NDim = cuImages.get_number_of_dimensions(); size_t CHA = 1; + size_t E0 = cuImages.get_size(0); + size_t E1 = cuImages.get_size(1); + size_t E2 = cuImages.get_size(2); + auto rmsize = recon_params.rmatrixSize_scanner; + if (E0 != rmsize.x || E1!=rmsize.y || E2!=rmsize.z){ + GDEBUG_STREAM("Cropping Images [E0 E1 E2] =[" << E0 << " " << E1 << " " << E2 <<"] != recon matrix [x y z] =[" << rmsize.x << " " << rmsize.y << " " << rmsize.z <<"]") + E0=rmsize.x;E1=rmsize.y;E2=rmsize.z; + } + cuNDArray cuimages_all =nhlbi_toolbox::utils::crop_to_recon_params_dims(cuImages,recon_params); size_t N = NDim > 3 ? cuImages.get_size(3) : 1; size_t S = NDim > 4 ? cuImages.get_size(4) : 1; size_t SLC = NDim > 5 ? cuImages.get_size(5) : 1; - GDEBUG_STREAM("CuImage SIZE " << NDim << " [RO E1 E2 CHA N S SLC] = [" << cuImages.get_size(0) << " " - << cuImages.get_size(1) << " " << cuImages.get_size(2) << " " << CHA << " " << N + GDEBUG_STREAM("CuImage SIZE " << NDim << " [RO E1 E2 CHA N S SLC] = [" << E0 << " " + << E1 << " " << E2 << " " << CHA << " " << N << " " << S << " " << SLC << "] "); IsmrmrdImageArray imarray_sense; auto images = hoNDArray>( - std::move(*boost::reinterpret_pointer_cast>>(cuImages.to_host()))); + std::move(*boost::reinterpret_pointer_cast>>(cuimages_all.to_host()))); auto tmp = hoNDArray>(images); tmp.reshape(tmp.get_size(0), tmp.get_size(1), tmp.get_size(2), 1, N, S, SLC); imarray_sense.data_ = tmp; @@ -663,6 +690,23 @@ class Noncart_recon_gadget imarray_sense.headers_(n, s, loc).image_index = offset + 1; imarray_sense.meta_[offset].append(GADGETRON_IMAGECOMMENT, image_comment.c_str()); imarray_sense.meta_[offset].append(GADGETRON_SEQUENCEDESCRIPTION, image_comment.c_str()); + imarray_sense.meta_[offset].append("ImageRowDir", imarray_sense.headers_(n, s, loc).read_dir[0]); + imarray_sense.meta_[offset].append("ImageRowDir", imarray_sense.headers_(n, s, loc).read_dir[1]); + imarray_sense.meta_[offset].append("ImageRowDir", imarray_sense.headers_(n, s, loc).read_dir[2]); + imarray_sense.meta_[offset].append("ImageColumnDir", imarray_sense.headers_(n, s, loc).phase_dir[0]); + imarray_sense.meta_[offset].append("ImageColumnDir", imarray_sense.headers_(n, s, loc).phase_dir[1]); + imarray_sense.meta_[offset].append("ImageColumnDir", imarray_sense.headers_(n, s, loc).phase_dir[2]); + /* + if (N >1){ + imarray_sense.meta_[offset].append("SiemensDicom_NumberInSeries", "long"); + imarray_sense.meta_[offset].append("SiemensDicom_NumberInSeries", long(N)); + imarray_sense.meta_[offset].append("SiemensDicom_ImageGroup", "long"); + imarray_sense.meta_[offset].append("SiemensDicom_ImageGroup", long(N)); + imarray_sense.meta_[offset].append("SiemensControl_CardiacRRInterval", "double"); + imarray_sense.meta_[offset].append("SiemensControl_CardiacRRInterval", double(N*25)); + } + */ + } } } @@ -738,6 +782,8 @@ class Noncart_recon_gadget NODE_PROPERTY(start_idx_nc, float, "With Collapse binning, only subsample NC", 0); NODE_PROPERTY(series_counter_initial, int, "series_counter_initial", 0); NODE_PROPERTY(save_avg, bool, "Saving Average image", true); + NODE_PROPERTY(save_csm, bool, "Saving CSM", false); + NODE_PROPERTY(save_intermediate_images, bool, "Saving intermediates image", true); NODE_PROPERTY(calculateKPRECOND, bool, "GT DCF of Kspace preconditioning", false); }; diff --git a/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Spiral_3DRT.cpp b/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Spiral_3DRT.cpp new file mode 100644 index 0000000..c1e6033 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Spiral_3DRT.cpp @@ -0,0 +1,359 @@ +/* + * SpiralMocoRecon.cpp + * + * Created on: September 17th, 2021 + * Author: Ahsan Javed + */ + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "../spiral/SpiralBuffer.h" +#include "../utils/gpu/cuda_utils.h" +//#include +#include "../../utils/util_functions.h" +#include "noncartesian_reconstruction.h" + +#include "reconParams.h" + +#include +#include +#include + +using namespace Gadgetron; +using namespace Gadgetron::Core; +using namespace Gadgetron::Indexing; +using namespace nhlbi_toolbox::reconstruction; + +class Spiral_3DRT : public ChannelGadget, + public ImageArraySendMixin +{ +public: + std::vector image_dims_; + uint64d3 image_dims_os_; + float oversampling_factor_; + float kernel_width_; + std::shared_ptr session_storage; + bool verbose; + + // "global heap variables".... do they need to be "public"? + boost::shared_ptr> reconstruction_; + boost::shared_ptr> csm_; + boost::shared_ptr> dcf_; + + Spiral_3DRT(const Core::Context &context, const Core::GadgetProperties &props) : ChannelGadget(context, props), + session_storage(context.storage.session) + { + kernel_width_ = 3; + oversampling_factor_ = oversampling_factor; + verbose = false; + } + + void process(InputChannel &in, + OutputChannel &out) override + { + + int selectedDevice = nhlbi_toolbox::utils::selectCudaDevice(); + unsigned int warp_size = cudaDeviceManager::Instance()->warp_size(); + size_t RO, E1, E2, CHA, N, S, SLC; + //nhlbi_toolbox::utils::enable_peeraccess(); + + ISMRMRD::AcquisitionHeader acqhdr; + + Gadgetron::reconParams recon_params; + + size_t recon_idx = 0; + size_t acq_count = 0; + + size_t out_counter = 0; + + auto n_tr_bin = repetitions; + if (n_tr_bin == 0) + { + n_tr_bin = ((header.encoding.at(0).trajectoryDescription.get().userParameterLong[0].value) * + (header.encoding[0].encodingLimits.kspace_encoding_step_2.get().maximum + 1)); + } + GDEBUG_STREAM("kspace_encode_1_max" << (header.encoding.at(0).encodingLimits.kspace_encoding_step_1.get().maximum + 1)); + GDEBUG_STREAM("kspace_encode_2_max" << (header.encoding.at(0).encodingLimits.kspace_encoding_step_2.get().maximum + 1)); + GDEBUG_STREAM("number of interleaves" << header.encoding.at(0).trajectoryDescription.get().userParameterLong[0].value); + GDEBUG_STREAM("n_tr_bin" << n_tr_bin); + + std::vector allAcq(n_tr_bin); + Gadgetron::hoNDArray headers(n_tr_bin); + + for (auto message : in) + { + auto &[head, data, traj] = message; + acqhdr = head; + + allAcq[recon_idx] = std::move(message); + headers[recon_idx] = std::move(head); + recon_idx ++; + acq_count ++; + + if (acq_count == n_tr_bin) + { + // this is the first time we are fully sampled, so we should get reconstruction, dcf, and csm. + // try to get them from the storage server, but if it fails.... generate anew. + + auto acq_toRecon = allAcq; + acq_toRecon.resize(n_tr_bin); + cudaSetDevice(selectedDevice); + auto &[headAcq, dataAcq, trajAcq] = acq_toRecon[0]; + RO = dataAcq.get_size(0); + CHA = dataAcq.get_size(1); + E2 = this->header.encoding.front().encodedSpace.matrixSize.z; + N = dataAcq.get_size(3); + S = 1; + SLC = 1; + + recon_params.numberChannels = CHA; + recon_params.RO = RO; + recon_params.ematrixSize = this->header.encoding.front().encodedSpace.matrixSize; + recon_params.rmatrixSize = this->header.encoding.front().reconSpace.matrixSize; + + std::ostringstream encoded_matrix_ss; + encoded_matrix_ss << recon_params.ematrixSize.x << "_" << recon_params.ematrixSize.y << + "_" << recon_params.ematrixSize.z << "_" << CHA; + std::string str_ematrix = encoded_matrix_ss.str(); + GDEBUG_STREAM("encoding matrix" << str_ematrix); + + // for now, set the rmatrix size to be a multiple of 32. + auto rmatrix = (int) (round(((float) recon_params.rmatrixSize.x) / 32 ) * 32); + recon_params.ematrixSize.x = rmatrix; + recon_params.ematrixSize.y = rmatrix; + recon_params.rmatrixSize.x = rmatrix; + recon_params.rmatrixSize.y = rmatrix; + + + GDEBUG_STREAM("recon_params.ematrixSize.x" << recon_params.ematrixSize.x); + GDEBUG_STREAM("recon_params.ematrixSize.y" << recon_params.ematrixSize.y); + GDEBUG_STREAM("recon_params.ematrixSize.z" << recon_params.ematrixSize.z); + GDEBUG_STREAM("recon_params.rmatrixSize.x" << recon_params.rmatrixSize.x); + GDEBUG_STREAM("recon_params.rmatrixSize.y" << recon_params.rmatrixSize.y); + GDEBUG_STREAM("recon_params.rmatrixSize.z" << recon_params.rmatrixSize.z); + + // TODO: PK revisit + // for now, set the rmatrix z size the same. This is because we want the full FOV for toeplitz (?) + recon_params.rmatrixSize.z = recon_params.ematrixSize.z; + recon_params.fov = this->header.encoding.front().encodedSpace.fieldOfView_mm; + recon_params.oversampling_factor_ = oversampling_factor_; + recon_params.kernel_width_ = kernel_width_; + recon_params.selectedDevice = selectedDevice; + recon_params.norm = 2; + recon_params.useIterativeDCWEstimated = true; + recon_params.oversampling_factor_dcf_ = oversampling_factor_; + recon_params.kernel_width_dcf_ = 5.5; + recon_params.iterations_dcf = 20; + this->initialize_encoding_space_limits(this->header); + + reconstruction_ = boost::make_shared>(recon_params); + + auto reconstruction = *reconstruction_; + + // reconstruct data. + auto [cuData, traj_csm, dcf_in] = reconstruction.organize_data(&acq_toRecon); + + //if (out_counter == 0) + // nhlbi_toolbox::utils::write_gpu_nd_array(cuData, "/opt/data/gt_data/kspace_0.complex"); + + std::ostringstream dcf_ss; + //dcf_ss << "dcf_" << recon_params.ematrixSize.x << "_" << recon_params.ematrixSize.y; + dcf_ss << "dcf_" << str_ematrix; + std::string dcf_str = dcf_ss.str(); + + std::ostringstream csm_ss; + //csm_ss << "csm_" << recon_params.ematrixSize.x << "_" << recon_params.ematrixSize.y; + csm_ss << "csm_" << str_ematrix; + std::string csm_str = csm_ss.str(); + + // Load the DCF. + try + { + auto dcf_cpu = this->session_storage->get_latest>(dcf_str); + + //PK funny: actually this debug message makes the try block fail gracefully + // TODO: figure out how to identify a 404 error and skip out gracefully, without + // having to resort to checking the dimensions to fail. + GDEBUG_STREAM("DCF did it fail? size: " << dcf_cpu->get_dimensions()->at(0)); + + nhlbi_toolbox::utils::write_cpu_nd_array(*dcf_cpu, "/opt/data/gt_data/dcf_scaled.real"); + + auto dcf_resize = hoNDArray(dcf_in.get_dimensions()); + for (size_t dcf_i = 0; dcf_i < dcf_resize.get_dimensions()->at(0); dcf_i++){ + dcf_resize.at(dcf_i) = dcf_cpu->at(dcf_i); + } + dcf_ = boost::make_shared>(cuNDArray(dcf_resize)); + } + catch (...) + { + //gen dcf if we can't load -> safety to regenerate! + GWARN_STREAM("Warning: DCF is being estimated from this scan, instead of loaded separately...."); + auto dcf = reconstruction.estimate_dcf(&traj_csm, &dcf_in); + dcf_ = boost::make_shared>(cuNDArray(dcf)); + } + + // Load the CSM. + try + { + auto csm_cpu = this->session_storage->get_latest>(csm_str); + csm_ = boost::make_shared>(cuNDArray(*csm_cpu)); + } + catch (...) + { + //gen csm if we can't load -> safety. But it could be undersampled. + GWARN_STREAM("Warning: CSM is being estimated from this scan, instead of loaded separately...."); + cuNDArray channel_images(reconstruction.get_recon_dims()); + { + reconstruction.reconstruct(&cuData, &channel_images, &traj_csm, dcf_.get()); + } + csm_ = reconstruction.generateRoemerCSM(&channel_images); + cuData.clear(); + traj_csm.clear(); + } + + // sqrt DCF if needed. + if (sqrt_dcf) + Gadgetron::sqrt_inplace(dcf_.get()); + + } + + if (recon_idx == n_tr_bin) + { + //re-set the recon_idx. We want to generate images with n_tr_bin over and over until acq_end. + recon_idx = 0; + out.push(allAcq); + GadgetronTimer timer_all("recon pipeline timer", false); + timer_all.start(); + + auto acq_toRecon = allAcq; + acq_toRecon.resize(n_tr_bin); + cudaSetDevice(selectedDevice); + auto &[headAcq, dataAcq, trajAcq] = acq_toRecon[0]; + + auto reconstruction = *reconstruction_; + + GadgetronTimer timer_organize("organize_data", false); + timer_organize.start(); + auto [cuData, traj, dcf_in] = reconstruction.organize_data(&acq_toRecon); + timer_organize.stop(); + + GadgetronTimer timer_dcf("dcf", false); + timer_dcf.start(); + auto dcf = *dcf_; + timer_dcf.stop(); + + cuNDArray channel_images(reconstruction.get_recon_dims()); + + { + GadgetronTimer recon_timer("reconstruction", false); + recon_timer.start(); + reconstruction.reconstruct(&cuData, &channel_images, &traj, &dcf); + cudaDeviceSynchronize(); + recon_timer.stop(); + } + + // save the data -- there must be a better way (I hope) + if (out_counter == 0) + { + std::ostringstream oss; + oss << "/opt/data/gt_data/kspace_" << out_counter << ".complex"; + std::string buf = oss.str(); + nhlbi_toolbox::utils::write_gpu_nd_array(cuData, buf); + } + + out_counter = out_counter + 1; + + GadgetronTimer coil_comb("coil_comb", false); + coil_comb.start(); + cuNDArray ci_cropped; + if (combine_csm == 0) + ci_cropped = channel_images; + else + { + channel_images *= *conj(csm_.get()); + + auto combined = sum(&channel_images, channel_images.get_number_of_dimensions() - 1); + + ci_cropped = reconstruction.crop_to_recondims(*combined); + + (*combined).clear(); + + } + channel_images.clear(); + coil_comb.stop(); + + GadgetronTimer finish("finish", false); + finish.start(); + auto images = hoNDArray>(std::move(*boost::reinterpret_pointer_cast>>(ci_cropped.to_host()))); + + using namespace Gadgetron::Indexing; + IsmrmrdImageArray imarray_sense; + imarray_sense.data_ = images; + nhlbi_toolbox::utils::attachHeadertoImageArray(imarray_sense, acqhdr, this->header); + prepare_image_array(imarray_sense, (size_t)0, ((int) series_counter), GADGETRON_IMAGE_REGULAR); + + timer_all.stop(); + finish.stop(); + + imarray_sense.headers_(0, 0, 0).user_int[0]=head.idx.repetition; + imarray_sense.headers_(0, 0, 0).data_type = ISMRMRD::ISMRMRD_CXFLOAT; + imarray_sense.headers_(0, 0, 0).image_type = ISMRMRD::ISMRMRD_IMTYPE_COMPLEX; + imarray_sense.acq_headers_ = headers; + + + out.push(imarray_sense); + } + } + } + +protected: + NODE_PROPERTY(oversampling_factor, float, "oversampling_factor", 1.25); + NODE_PROPERTY(crop_begin, size_t, "crop_begin", 0); + NODE_PROPERTY(crop_end, size_t, "crop_end", 0); + NODE_PROPERTY(repetitions, size_t, "repetitions", 0); + NODE_PROPERTY(combine_csm, size_t, "combine_csm", 1); + NODE_PROPERTY(send_csm, size_t, "send_csm", 0); + NODE_PROPERTY(sqrt_dcf, size_t, "sqrt_dcf", 0); + int series_counter = 0; +}; + +GADGETRON_GADGET_EXPORT(Spiral_3DRT) \ No newline at end of file diff --git a/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Spiral_Rovir.cpp b/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Spiral_Rovir.cpp new file mode 100644 index 0000000..d4d0f2c --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/noncart_recon/Spiral_Rovir.cpp @@ -0,0 +1,334 @@ +/* + * Spiral_Rovir.cpp + * + * Author: Prakash Kumar + */ + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "../spiral/SpiralBuffer.h" +#include "../utils/gpu/cuda_utils.h" +//#include +#include "../../utils/util_functions.h" +#include "noncartesian_reconstruction.h" + +#include "reconParams.h" + +#include +#include +#include + +using namespace Gadgetron; +using namespace Gadgetron::Core; +using namespace Gadgetron::Indexing; +using namespace nhlbi_toolbox::reconstruction; + +class Spiral_Rovir : public ChannelGadget, + public ImageArraySendMixin +{ +public: + std::vector image_dims_; + uint64d3 image_dims_os_; + float oversampling_factor_; + float kernel_width_; + std::shared_ptr session_storage; + bool verbose; + + // "global heap variables".... do they need to be "public"? + //boost::shared_ptr> reconstruction_; + //boost::shared_ptr> csm_; + //boost::shared_ptr> dcf_; + + Spiral_Rovir(const Core::Context &context, const Core::GadgetProperties &props) : ChannelGadget(context, props), + session_storage(context.storage.session) + { + kernel_width_ = 3; + oversampling_factor_ = oversampling_factor; + verbose = false; + } + + void process(InputChannel &in, + OutputChannel &out) override + { + + int selectedDevice = nhlbi_toolbox::utils::selectCudaDevice(); + unsigned int warp_size = cudaDeviceManager::Instance()->warp_size(); + size_t RO, E1, E2, CHA, N, S, SLC; + //nhlbi_toolbox::utils::enable_peeraccess(); + + ISMRMRD::AcquisitionHeader acqhdr; + + Gadgetron::reconParams recon_params; + + size_t recon_idx = 0; + size_t acq_count = 0; + + size_t out_counter = 0; + + auto n_tr_bin = ((header.encoding.at(0).encodingLimits.kspace_encoding_step_1.get().maximum + 1) * + (header.encoding[0].encodingLimits.kspace_encoding_step_2.get().maximum + 1)) * + (header.encoding[0].encodingLimits.repetition.get().maximum + 1); + (header.encoding[0].encodingLimits.average.get().maximum + 1); + + GDEBUG_STREAM("kspace_encode_1_max" << (header.encoding.at(0).encodingLimits.kspace_encoding_step_1.get().maximum + 1)); + GDEBUG_STREAM("kspace_encode_2_max" << (header.encoding.at(0).encodingLimits.kspace_encoding_step_2.get().maximum + 1)); + //GDEBUG_STREAM("number of interleaves" << header.encoding.at(0).trajectoryDescription.get().userParameterLong[0].value); + GDEBUG_STREAM("n_tr_bin" << n_tr_bin); + + std::vector allAcq(n_tr_bin); + Gadgetron::hoNDArray headers(n_tr_bin); + + std::vector csmAcq(n_tr_bin); + Gadgetron::hoNDArray csmHeaders(n_tr_bin); + + for (auto message : in) + { + if (acq_count < n_tr_bin) + { + auto &[head, data, traj] = message; + acqhdr = head; + csmAcq[acq_count] = std::move(message); + csmHeaders[acq_count] = std::move(head); + + if (acq_count == (n_tr_bin - 1)) + { + auto acq_toRecon = csmAcq; + acq_toRecon.resize(n_tr_bin); + cudaSetDevice(selectedDevice); + auto &[headAcq, dataAcq, trajAcq] = acq_toRecon[0]; + RO = dataAcq.get_size(0); + CHA = dataAcq.get_size(1); + E2 = this->header.encoding.front().encodedSpace.matrixSize.z; + N = dataAcq.get_size(3); + S = 1; + SLC = 1; + + recon_params.numberChannels = CHA; + recon_params.RO = RO; + recon_params.ematrixSize = this->header.encoding.front().encodedSpace.matrixSize; + recon_params.rmatrixSize = this->header.encoding.front().encodedSpace.matrixSize; //this->header.encoding.front().reconSpace.matrixSize; + + std::ostringstream encoded_matrix_ss; + encoded_matrix_ss << recon_params.ematrixSize.x << "_" << recon_params.ematrixSize.y << + "_" << recon_params.ematrixSize.z << "_" << CHA; + std::string str_ematrix = encoded_matrix_ss.str(); + GDEBUG_STREAM("encoding matrix" << str_ematrix); + + // for now, set the ematrix size to be a multiple of 32. + float matrix_overgrid = 2; + auto rmatrix = (int) (ceil(((float) recon_params.rmatrixSize.x * matrix_overgrid) / 32 ) * 32); + recon_params.ematrixSize.x = rmatrix; + recon_params.ematrixSize.y = rmatrix; + //recon_params.rmatrixSize.x = rmatrix; + //recon_params.rmatrixSize.y = rmatrix; + + // TODO: PK revisit + // for now, set the rmatrix z size the same. This is because we want the full FOV for toeplitz (?) + recon_params.rmatrixSize.z = recon_params.ematrixSize.z; + recon_params.rmatrixSize.z = recon_params.ematrixSize.z; + recon_params.fov = this->header.encoding.front().encodedSpace.fieldOfView_mm; + recon_params.oversampling_factor_ = oversampling_factor_; + recon_params.kernel_width_ = kernel_width_; + recon_params.selectedDevice = selectedDevice; + recon_params.norm = 2; + recon_params.useIterativeDCWEstimated = false; + recon_params.oversampling_factor_dcf_ = oversampling_factor_; + recon_params.kernel_width_dcf_ = 3; //3.2 + recon_params.iterations_dcf = 10; // 10 + this->initialize_encoding_space_limits(this->header); + + auto reconstruction = noncartesian_reconstruction<3>(recon_params); + + // reconstruct data. + auto [cuData, traj_csm, dcf_in] = reconstruction.organize_data(&acq_toRecon); + + out_counter = 1; + + std::ostringstream dcf_ss; + dcf_ss << "dcf_" << str_ematrix; + std::string dcf_str = dcf_ss.str(); + + std::ostringstream csm_ss; + csm_ss << "csm_" << str_ematrix; + std::string csm_str = csm_ss.str(); + + std::ostringstream csm_over_ss; + csm_over_ss << "csm_over_" << str_ematrix; + std::string csm_over_str = csm_over_ss.str(); + + std::ostringstream csm_image_ss; + csm_image_ss << "csm_image_" << str_ematrix; + std::string csm_image_str = csm_image_ss.str(); + + std::ostringstream csm_image_mc_ss; + csm_image_mc_ss << "csm_image_mc_" << str_ematrix; + std::string csm_image_mc_str = csm_image_mc_ss.str(); + + auto dcf = reconstruction.estimate_dcf(&traj_csm, &dcf_in); + + square_inplace(&dcf); + + auto dcf_cpu = dcf.to_host(); + //auto dcf_max = dcf_cpu->at((size_t) std::floor((float)RO * 0.8)); + //Gadgetron::clamp(dcf_cpu.get(), 0, dcf_max); + + //auto dcf_max_array = Gadgetron::hoNDArray(*dcf_cpu->get_dimensions().get()); + //Gadgetron::fill(&dcf_max_array, Gadgetron::max(dcf_cpu.get())); + //Gadgetron::divide(dcf_cpu.get(), &dcf_max_array, dcf_cpu.get()); + //auto scale = Gadgetron::max(dcf_cpu.get()); + + this->session_storage->store(dcf_str, *dcf_cpu); + + cuNDArray channel_images(reconstruction.get_recon_dims()); + { + + reconstruction.reconstruct(&cuData, &channel_images, &traj_csm, &dcf); + } + + + cuData.clear(); + traj_csm.clear(); + + auto csm_ = reconstruction.generateRoemerCSM(&channel_images); + + // GDEBUG_STREAM("SAVING CSM to File....") + // // save the CSM for PYTHON + // std::ostringstream csm_ss_python; + // //csm_ss_python << "/opt/data/gt_data/csm_" << csm_->get_dimensions()->at(0) << "_" << csm_->get_dimensions()->at(1) << ".complex"; + // csm_ss_python << "/opt/data/gt_data/csm_" << str_ematrix << ".complex"; + // std::string buf_csm = csm_ss_python.str(); + // nhlbi_toolbox::utils::write_gpu_nd_array(*csm_, buf_csm); + + + + if (send_csm == 1) + { + auto csm_send =hoNDArray>(std::move(*boost::reinterpret_pointer_cast>>(csm_->to_host()))); + using namespace Gadgetron::Indexing; + IsmrmrdImageArray imarray_sense; + imarray_sense.data_ = csm_send; + nhlbi_toolbox::utils::attachHeadertoImageArray(imarray_sense, acqhdr, this->header); + prepare_image_array(imarray_sense, (size_t)0, ((int)series_counter), GADGETRON_IMAGE_GFACTOR); + imarray_sense.headers_(0, 0, 0).user_int[0]=head.idx.repetition; + imarray_sense.headers_(0, 0, 0).data_type = ISMRMRD::ISMRMRD_CXFLOAT; + imarray_sense.headers_(0, 0, 0).image_type = ISMRMRD::ISMRMRD_IMTYPE_COMPLEX; + + out.push(imarray_sense); + series_counter++; + } + + // save the multi-coil CSM image to the storage server. This is useful for ROVIR. + // we also save it before cropping the dimensions. + auto csm_image_mc_cpu = std::move((channel_images).to_host()); + this->session_storage->store(csm_image_mc_str, *csm_image_mc_cpu); + + GDEBUG_STREAM("CSM estimation done. Sending images to storage server and/or next gadget....") + cuNDArray ci_cropped; + if (combine_csm == 0) + ci_cropped = channel_images; + else + { + channel_images *= *conj(csm_.get()); + + auto combined = sum(&channel_images, channel_images.get_number_of_dimensions() - 1); + + //HACK: update the supposedly private image_dims_ var (but it's public) so cropping actually does something. + reconstruction.image_dims_.at(0) = (int) recon_params.rmatrixSize.x * 1.5; + reconstruction.image_dims_.at(1) = (int) recon_params.rmatrixSize.y * 1.5; + + ci_cropped = reconstruction.crop_to_recondims_centered(*combined); + (*combined).clear(); + + } + + // save the CSM for the storage server + auto csm_cpu = std::move((*csm_).to_host()); + this->session_storage->store(csm_over_str, *csm_cpu); + + auto csm_cropped = reconstruction.crop_to_recondims_centered(*csm_); + auto csm_cpu_cropped = std::move((csm_cropped).to_host()); + this->session_storage->store(csm_str, *csm_cpu_cropped); + + // // SAVE THE CSM IMAGE for python + // std::ostringstream csm_image_ss_python; + // //csm_image_ss_python << "/opt/data/gt_data/csm_image_" << csm_->get_dimensions()->at(0) << "_" << csm_->get_dimensions()->at(1) << ".complex"; + // csm_image_ss_python << "/opt/data/gt_data/csm_image_" << str_ematrix << ".complex"; + // std::string buf_csm_image = csm_image_ss_python.str(); + // nhlbi_toolbox::utils::write_gpu_nd_array(ci_cropped, buf_csm_image); + + // save the CSM image to the storage server + auto csm_image_cpu = std::move((ci_cropped).to_host()); + this->session_storage->store(csm_image_str, *csm_image_cpu); + + channel_images.clear(); + + auto images = hoNDArray>(std::move(*boost::reinterpret_pointer_cast>>(ci_cropped.to_host()))); + + using namespace Gadgetron::Indexing; + IsmrmrdImageArray image_out; + image_out.data_ = images; + nhlbi_toolbox::utils::attachHeadertoImageArray(image_out, acqhdr, this->header); + prepare_image_array(image_out, (size_t)0, ((int) series_counter), GADGETRON_IMAGE_REGULAR); + + image_out.headers_(0, 0, 0).user_int[0]=head.idx.repetition; + image_out.headers_(0, 0, 0).data_type = ISMRMRD::ISMRMRD_CXFLOAT; + image_out.headers_(0, 0, 0).image_type = ISMRMRD::ISMRMRD_IMTYPE_COMPLEX; + image_out.acq_headers_ = csmHeaders; + GDEBUG_STREAM("ACQ HEADER EXPORT: " << image_out.acq_headers_->get_dimensions()->at(0)); + + out.push(image_out); + } + } + acq_count ++; + } + } +protected: + //NODE_PROPERTY(oversampling_factor, float, "oversampling_factor", 2.1); + NODE_PROPERTY(oversampling_factor, float, "oversampling_factor", 2.1); + NODE_PROPERTY(crop_begin, size_t, "crop_begin", 0); + NODE_PROPERTY(crop_end, size_t, "crop_end", 0); + NODE_PROPERTY(combine_csm, size_t, "combine_csm", 1); + NODE_PROPERTY(send_csm, size_t, "send_csm", 0); + int series_counter = 0; +}; + +GADGETRON_GADGET_EXPORT(Spiral_Rovir) diff --git a/toolboxes/nhlbi_gt_toolbox/gadgets/utility_gadgets/PrepreconParams.cpp b/toolboxes/nhlbi_gt_toolbox/gadgets/utility_gadgets/PrepreconParams.cpp index ed7727f..87088ba 100644 --- a/toolboxes/nhlbi_gt_toolbox/gadgets/utility_gadgets/PrepreconParams.cpp +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/utility_gadgets/PrepreconParams.cpp @@ -117,8 +117,12 @@ class PrepreconParams : public ChannelGadget recon_params.gcc_coils = gcc_coils; recon_params.selectedDevice = selectedGPUs[0]; - recon_params.selectedDevices = selectedGPUs; - + recon_params.selectedDevices_solver = selectedGPUs; + if(minGPU_utilization){ + recon_params.selectedDevices = {selectedGPUs[0]}; + }else{ + recon_params.selectedDevices = selectedGPUs; + } recon_params.try_channel_gridding=try_channel_gridding; @@ -162,9 +166,14 @@ class PrepreconParams : public ChannelGadget auto mr_y=size_t(ceil((matOSP_vector[1]*mr_dy)/warp_vector[1]))*warp_vector[1]; auto mr_z=size_t(ceil((matOSP_vector[2]*mr_dz)/warp_vector[2]))*warp_vector[2]; + auto mr_x_scanner = size_t(ceil(scannerOSP_vector[0]*mr_dx)); + auto mr_y_scanner = size_t(ceil(scannerOSP_vector[1]*mr_dy)); + auto mr_z_scanner = size_t(ceil(scannerOSP_vector[2]*mr_dz)); + if (recon_params.ematrixSize.z ==1 && recon_params.rmatrixSize.z ==1){ auto mr_z=1; auto me_z=1; + auto mr_z_scanner=1; } recon_params.ematrixSize.x = me_x; @@ -174,14 +183,21 @@ class PrepreconParams : public ChannelGadget recon_params.rmatrixSize.x = mr_x; recon_params.rmatrixSize.y = mr_y; recon_params.rmatrixSize.z = mr_z; + + recon_params.rmatrixSize_scanner.x = mr_x_scanner; + recon_params.rmatrixSize_scanner.y = mr_y_scanner; + recon_params.rmatrixSize_scanner.z = mr_z_scanner; + + recon_params.fov = this->header.encoding.front().encodedSpace.fieldOfView_mm; - recon_params.fov.x=recon_params.fov.x*(mr_x/mr_dx); - recon_params.fov.y=recon_params.fov.y*(mr_y/mr_dy); - recon_params.fov.z=recon_params.fov.z*(mr_z/mr_dz); + recon_params.fov.x=recon_params.fov.x*(mr_x_scanner/mr_dx); + recon_params.fov.y=recon_params.fov.y*(mr_y_scanner/mr_dy); + recon_params.fov.z=recon_params.fov.z*(mr_z_scanner/mr_dz); GDEBUG_STREAM("Encoded Matrix: X " << recon_params.ematrixSize.x << " Y " << recon_params.ematrixSize.y << " Z " << recon_params.ematrixSize.z ); GDEBUG_STREAM("Recon Matrix: X " << recon_params.rmatrixSize.x << " Y " << recon_params.rmatrixSize.y << " Z " << recon_params.rmatrixSize.z); + GDEBUG_STREAM("Recon Matrix scanner: X " << recon_params.rmatrixSize_scanner.x << " Y " << recon_params.rmatrixSize_scanner.y << " Z " << recon_params.rmatrixSize_scanner.z); GDEBUG_STREAM("Recon FOV: X " << recon_params.fov.x << " Y " << recon_params.fov.y << " Z " << recon_params.fov.z); @@ -201,6 +217,11 @@ class PrepreconParams : public ChannelGadget recon_params_avg.omatrixSize.y =recon_params.omatrixSize.y; recon_params_avg.omatrixSize.z =recon_params.omatrixSize.z; + recon_params_avg.rmatrixSize_scanner.x =recon_params.rmatrixSize_scanner.x; + recon_params_avg.rmatrixSize_scanner.y =recon_params.rmatrixSize_scanner.y; + recon_params_avg.rmatrixSize_scanner.z =recon_params.rmatrixSize_scanner.z; + + //FOV recon_params_avg.fov.x =recon_params.fov.x; recon_params_avg.fov.y =recon_params.fov.y; @@ -209,10 +230,18 @@ class PrepreconParams : public ChannelGadget GDEBUG_STREAM("AVERAGE RECON PARAMS" ) GDEBUG_STREAM("Encoded Matrix: X " << recon_params_avg.ematrixSize.x << " Y " << recon_params_avg.ematrixSize.y << " Z " << recon_params_avg.ematrixSize.z ); GDEBUG_STREAM("Recon Matrix: X " << recon_params_avg.rmatrixSize.x << " Y " << recon_params_avg.rmatrixSize.y << " Z " << recon_params_avg.rmatrixSize.z); + GDEBUG_STREAM("Recon Matrix scanner: X " << recon_params_avg.rmatrixSize_scanner.x << " Y " << recon_params_avg.rmatrixSize_scanner.y << " Z " << recon_params_avg.rmatrixSize_scanner.z); GDEBUG_STREAM("Recon FOV: X " << recon_params_avg.fov.x << " Y " << recon_params_avg.fov.y << " Z " << recon_params_avg.fov.z); - - - + /* + std::ostringstream str_lambda_spatial,str_lambda_time; + str_lambda_spatial << std::scientific << std::setprecision(2) << recon_params.lambda_spatial; + str_lambda_time << std::scientific << std::setprecision(2) << recon_params.lambda_time; + std::string img_parameters_name = std::string("r") + std::string("_ite_") + + std::to_string(recon_params.iterations) + std::string("_ls_") + + str_lambda_spatial.str() + std::string("_lt") + + str_lambda_time.str(); + GDEBUG_STREAM("IMAGE_NAME"< recon_params_avg.gcc_coils = gcc_coils; recon_params_avg.selectedDevice = selectedGPUs[0]; - recon_params_avg.selectedDevices = selectedGPUs; - + recon_params_avg.selectedDevices_solver = selectedGPUs; + if(minGPU_utilization){ + recon_params_avg.selectedDevices = {selectedGPUs[0]}; + }else{ + recon_params_avg.selectedDevices = selectedGPUs; + } recon_params_avg.try_channel_gridding=try_channel_gridding; GDEBUG_STREAM("CHANNEL GRIDDING " << recon_params_avg.try_channel_gridding << " " << recon_params.try_channel_gridding); size_t RO = 0; @@ -379,6 +412,7 @@ class PrepreconParams : public ChannelGadget NODE_PROPERTY(matOSP_vector, std::vector, "Large FOV factor",(std::vector{ 1, 1, 1})); // Vector of scaling factors for large Field Of View (FOV) NODE_PROPERTY(downsampling_vector, std::vector, "Downsampling factor plane(x,y) and z)",(std::vector{ 1, 1})); // Downsampling factors for plane (x, y) and z dimension NODE_PROPERTY(warpCUDA_vector, std::vector, "Warp CUDA (32)",(std::vector{ true, true, false})); // Flags for respecting CUDA size of warp ( matrix x,y,z should be a multiple of 32) + NODE_PROPERTY(scannerOSP_vector, std::vector, "FOV factor for reconstruction on the scanner",(std::vector{ 1, 1, 1})); // Vector of scaling factors for Field Of View (FOV) reconstructed on the scanner NODE_PROPERTY(is3D, bool, "is 3D not stack of 2D", false); // Flag indicating if data is 3D non cartesian (not a stack of stars, spirals) //NUFFT parameters @@ -415,6 +449,7 @@ class PrepreconParams : public ChannelGadget NODE_PROPERTY(use_gcc, bool, "use_gcc", false); // Flag to use GCC calibration NODE_PROPERTY(gcc_coils, size_t, "gcc_coils", 6); // Number of coils for GCC calibration NODE_PROPERTY(selectedDevices_STR, std::string, "String list of GPU device (0-N:device i, -1 : let GT choose, -2: No Device)", "-1 -2"); // String for selecting GPU devices + NODE_PROPERTY(minGPU_utilization, bool, "Only use multiple GPUs for solver",false); // Flag to use multiple GPUs only for the solver part of the reconstruction //NODE_PROPERTY(repeated_GPUs, unsigned int, "Repeat eligible GPUs x times",1); NODE_PROPERTY(maxIteRegistration, int, "Number of Iterations with estimation registration", 0); // Number of iterations for registration with estimation NODE_PROPERTY(try_channel_gridding, bool, "try_gridding over all channels", true); // Flag to enable gridding over all channels diff --git a/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.cpp b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.cpp index 6005e8a..562d6b6 100644 --- a/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.cpp +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.cpp @@ -20,7 +20,7 @@ using namespace Gadgetron; using namespace Gadgetron::Core; using namespace nhlbi_toolbox::utils; WaveformToTrajectory::WaveformToTrajectory(const Core::Context &context, const Core::GadgetProperties &props) - : ChannelGadget(context, props), header{context.header}, trajParams{context.header} {} + : ChannelGadget(context, props), header{context.header}, trajParams{context.header}, measurement{context.storage.measurement} {} namespace { @@ -131,7 +131,14 @@ void WaveformToTrajectory ::process( this->girf_kernel = nhlbi_toolbox::corrections::readGIRFKernel(GIRF_folder + "GIRF_fmax_"); // AJ fix for now else this->girf_kernel = nhlbi_toolbox::corrections::readGIRFKernel(GIRF_folder + "GIRF"); // Read GIRF Kernel from file - + + + // Set clock shift for GIRF correction + trajParams.set_clock_shift(this->clock_shift_s); + + // Set debug folder for writing out waveforms and trajectories if set + trajParams.set_debug_folder(this->debug_folder); + // Extract sampling time from the sequence ISMRMRD::TrajectoryDescription traj_desc; @@ -161,6 +168,9 @@ void WaveformToTrajectory ::process( } { + GDEBUG_STREAM("WaveformToTrajectory: GIRF parameters: perform_GIRF:" << perform_GIRF << " GIRF_folder:" << GIRF_folder << " GIRF_samplingtime:" << GIRF_samplingtime << "clock shift" << clock_shift_s); + GDEBUG_STREAM("Trajectory generation parameters: generateTraj:" << generateTraj << " attachWaveform" << attachWaveform << " realTime:" << realTime <<" acceleration_factor:" << acceleration_factor); + GDEBUG_STREAM(" setPre:" << setPre << " pre_cutoff_manual:" << pre_cutoff_manual << " crop_index_st:" << crop_index_st); GadgetronTimer timer("WaveformToTrajectory"); // #pragma omp parallel // #pragma omp for @@ -237,7 +247,8 @@ void WaveformToTrajectory ::process( GDEBUG_STREAM("rotations:" <(tw_gen); dcwgen = std::get<1>(tw_gen); + this->measurement->store("trajectory", trajgen); nhlbi_toolbox::utils::normalize_trajectory(&trajgen); Tsamp_us = trajParams.get_Tsampling_us(); } @@ -278,6 +300,7 @@ void WaveformToTrajectory ::process( traj_not_generated = false; trajgen = std::get<0>(tw_gen); dcwgen = std::get<1>(tw_gen); + this->measurement->store("trajectory", trajgen); nhlbi_toolbox::utils::normalize_trajectory(&trajgen); Tsamp_us = trajParams.get_Tsampling_us(); } @@ -299,39 +322,12 @@ void WaveformToTrajectory ::process( trajectory_and_weights = traj_dcw; for (int ii = 0; ii < trajectory_and_weights.get_size(1); ii++) { - trajectory_and_weights(0, ii) = traj_dcw(0, ii); - trajectory_and_weights(1, ii) = traj_dcw(1, ii); + trajectory_and_weights(0, ii) = (traj_dcw(0, ii)> 0.5f) ? 0.5f : ((traj_dcw(0, ii) < -0.5f) ? -0.5f : traj_dcw(0, ii)); + trajectory_and_weights(1, ii) = (traj_dcw(1, ii)> 0.5f) ? 0.5f : ((traj_dcw(1, ii) < -0.5f) ? -0.5f : traj_dcw(1, ii)); trajectory_and_weights(2, ii) = traj_dcw(2, ii); - size_t num = 0; - if (abs(trajectory_and_weights(0, ii)) > 0.5f || abs(trajectory_and_weights(1, ii)) > 0.5f) + if ((this->header.encoding.front().encodedSpace.matrixSize.z > 1)) // is 3D { - if (ii == 0) - GERROR("To Prevent recon failure setting to ±0.5 \n"); - - if (trajectory_and_weights(0, ii) > 0.5) - { - // GDEBUG_STREAM(" trajectory_and_weights(0, ii):" << trajectory_and_weights(0, ii)); - - trajectory_and_weights(0, ii) = 0.5; - } - else if (trajectory_and_weights(0, ii) < -0.5) - { - // GDEBUG_STREAM(" trajectory_and_weights(0, ii):" << trajectory_and_weights(0, ii)); - - trajectory_and_weights(0, ii) = -0.5; - } - if (trajectory_and_weights(1, ii) > 0.5) - { - // GDEBUG_STREAM(" trajectory_and_weights(1, ii):" << trajectory_and_weights(1, ii)); - trajectory_and_weights(1, ii) = 0.5; - } - else if (trajectory_and_weights(1, ii) < -0.5) - { - // GDEBUG_STREAM(" trajectory_and_weights(1, ii):" << trajectory_and_weights(1, ii)); - - trajectory_and_weights(1, ii) = -0.5; - } - num++; + trajectory_and_weights(3, ii) = traj_dcw(3, ii); } } if (head.discard_pre == 0 && setPre) @@ -386,38 +382,14 @@ void WaveformToTrajectory ::process( trajectory_and_weights.fill(0.0); for (int ii = 0; ii < trajectory_and_weights.get_size(1); ii++) { - trajectory_and_weights(0, ii) = temp(0, ii); - trajectory_and_weights(1, ii) = temp(1, ii); + trajectory_and_weights(0, ii) = (temp(0, ii)> 0.5f) ? 0.5f : ((temp(0, ii) < -0.5f) ? -0.5f : temp(0, ii)); // need to clip to 0.5 to prevent recon failure + trajectory_and_weights(1, ii) = (temp(1, ii)> 0.5f) ? 0.5f : ((temp(1, ii) < -0.5f) ? -0.5f : temp(1, ii)); // need to clip to 0.5 to prevent recon failure trajectory_and_weights(2, ii) = temp(2, ii); - if (!perform_GIRF) // only do this if not doing apply girf else apply girf takes care of this - { - size_t num = 0; - if (abs(trajectory_and_weights(0, ii)) > 0.5f || abs(trajectory_and_weights(1, ii)) > 0.5f) - { - if (ii == 0) - GERROR("To Prevent recon failure setting to ±0.5 \n"); - - if (trajectory_and_weights(0, ii) > 0.5f) - { - trajectory_and_weights(0, ii) = 0.5f; - } - else if (trajectory_and_weights(0, ii) < -0.5f) - { - trajectory_and_weights(0, ii) = -0.5f; - } - if (trajectory_and_weights(1, ii) > 0.5f) - { - trajectory_and_weights(1, ii) = 0.5f; - } - else if (trajectory_and_weights(1, ii) < -0.5f) - { - trajectory_and_weights(1, ii) = -0.5f; - } - num++; - } - } + if (header.encoding.front().encodedSpace.matrixSize.z > 1) trajectory_and_weights(3, ii) = temp(3, ii); + auto zencoding = float(-0.5 + head.idx.kspace_encode_step_2 * 1 / ((float)header.encoding.front().encodedSpace.matrixSize.z)); + trajectory_and_weights(2,ii) = zencoding; } if (perform_GIRF) // do_girf @@ -526,7 +498,7 @@ void WaveformToTrajectory::prepare_trajectory_from_waveforms(Core::Waveform &gra auto gradients_interpolated = zeroHoldInterpolation(gradients, upsampleFactor); if (perform_GIRF) - gradients_interpolated = nhlbi_toolbox::corrections::girf_correct(gradients_interpolated, this->girf_kernel, rotation_matrix, 2e-6, 10e-6, 0.85e-6); + gradients_interpolated = nhlbi_toolbox::corrections::girf_correct(gradients_interpolated, this->girf_kernel, rotation_matrix, 2e-6, this->GIRF_samplingtime, this->clock_shift_s); auto zencoding = float(-0.5 + head.idx.kspace_encode_step_2 * 1 / ((float)this->header.encoding.front().encodedSpace.matrixSize.z)); trajectory_and_weights(0, 0) = (gradients_interpolated(0)[0]) * GAMMA * 10 * head.sample_time_us * 1e-6 * kspace_scaling; @@ -675,7 +647,7 @@ void WaveformToTrajectory::applyGIRF(hoNDArray &trajectory_and_weights, I auto dcw_sep = std::move(*std::get<1>(traj_dcw).get()); auto gradients = nhlbi_toolbox::utils::traj2grad_3D2D(traj_sep, kspace_scaling, head); - gradients = nhlbi_toolbox::corrections::girf_correct(gradients, girf_kernel, rotation_matrix, head.sample_time_us * 1e-6, 10e-6, 0.85e-6); + gradients = nhlbi_toolbox::corrections::girf_correct(gradients, girf_kernel, rotation_matrix, head.sample_time_us * 1e-6, this->GIRF_samplingtime, this->clock_shift_s); auto zencoding = float(-0.5 + head.idx.kspace_encode_step_2 * 1 / ((float)header.encoding.front().encodedSpace.matrixSize.z)); trajectory_and_weights(0, 0) = (gradients(0)[0]) * GAMMA * 10 * head.sample_time_us * 1e-6 * kspace_scaling; @@ -724,7 +696,7 @@ void WaveformToTrajectory::applyGIRF(hoNDArray &trajectory_and_weights, I auto dcw_sep = std::move(*std::get<1>(traj_dcw).get()); auto gradients = nhlbi_toolbox::utils::traj2grad(traj_sep, kspace_scaling, head); - gradients = nhlbi_toolbox::corrections::girf_correct(gradients, girf_kernel, rotation_matrix, head.sample_time_us * 1e-6, 10e-6, 0.85e-6); + gradients = nhlbi_toolbox::corrections::girf_correct(gradients, girf_kernel, rotation_matrix, head.sample_time_us * 1e-6, this->GIRF_samplingtime, this->clock_shift_s); trajectory_and_weights(0, 0) = (gradients(0)[0]) * GAMMA * 10 * head.sample_time_us * 1e-6 * kspace_scaling; trajectory_and_weights(1, 0) = (gradients(0)[1]) * GAMMA * 10 * head.sample_time_us * 1e-6 * kspace_scaling; diff --git a/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.h b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.h index b3f37db..bff6a33 100644 --- a/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.h +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.h @@ -37,6 +37,8 @@ using namespace Gadgetron; std::map> trajectory_map; std::map gradient_wave_store; size_t curAvg=0; + + std::shared_ptr measurement; protected: @@ -44,13 +46,14 @@ using namespace Gadgetron; NODE_PROPERTY(perform_GIRF, bool, " Perform GIRF", false); NODE_PROPERTY(GIRF_folder, std::string, "Path where GIRF Data is stored", "/opt/GIRF/"); - NODE_PROPERTY(generateTraj, bool, "generate trajectories", false); NODE_PROPERTY(GIRF_samplingtime, float, "girf sampling time", 10e-6); + NODE_PROPERTY(clock_shift_s, float, "CLOCK SHIFT", 0.85e-6); NODE_PROPERTY(crop_index_st, size_t, "start index to crop acquisition data", 20); + NODE_PROPERTY(generateTraj, bool, "generate trajectories", false); NODE_PROPERTY(attachWaveform, bool, "attachWaveforms", true); NODE_PROPERTY(setPre, bool, "setPre", false); NODE_PROPERTY(realTime, bool, "realTime", false); - + NODE_PROPERTY(debug_folder, std::string, "If set, the debug output will be written out", ""); // debug folder for waveforms and trajectories ("/opt/data/gt_data/") NODE_PROPERTY(pre_cutoff_manual, size_t, "pre_cutoff_manual", 20); NODE_PROPERTY(acceleration_factor, size_t, "acceleration_factor", 1); // bug fix for a sequence bug with acc diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/TrajectoryParameters_lit.cpp b/toolboxes/nhlbi_gt_toolbox/spiral/TrajectoryParameters_lit.cpp index 91d6603..2b5afda 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/TrajectoryParameters_lit.cpp +++ b/toolboxes/nhlbi_gt_toolbox/spiral/TrajectoryParameters_lit.cpp @@ -13,6 +13,12 @@ namespace Gadgetron std::pair, hoNDArray> TrajectoryParameters_lit::calculate_trajectories_and_weight(const ISMRMRD::AcquisitionHeader &acq_header) { + + + bool debug_flag = !(this->debug_folder_.empty()); + + + // Two-fov percentage definition for variable density design if (strstr(systemModel.c_str(),"MAGNETOM eMeRge-XL") || strstr(systemModel.c_str(),"MAGNETOM Sola")) @@ -56,8 +62,7 @@ namespace Gadgetron fov_vds_temp[1] = std::round((-1 * fov_ * (1.0 - 1.0 * (vds_factor_ / 100.0)))*1000.0f)/1000.0f; // fov_vds_ = fov_vds_temp; - GDEBUG_STREAM("fov_vds_temp[0]:" << fov_vds_temp[0]); - + GDEBUG_STREAM("fov_vds_temp[0]:" << fov_vds_temp[0]); GDEBUG_STREAM("fov_vds_temp[1]:" << fov_vds_temp[1]); } @@ -70,18 +75,23 @@ namespace Gadgetron double sample_time = (1.0f * Tsamp_ns_) * 1.0e-9; // auto base_gradients = calculate_vds(smax_, gmax_, sample_time, sample_time, Nints_, &fov_, nfov, krmax_, ngmax, acq_header.number_of_samples); auto base_gradients = nhlbi_toolbox::Spiral::calculate_vds(smax_, gmax_, sample_time, sample_time, Nints_, fov_vds_, nfov, krmax_, ngmax, acq_header.number_of_samples); - auto filename = "/opt/data/gt_data/base_gradients.real2"; - nhlbi_toolbox::utils::write_cpu_nd_array(base_gradients, filename); + if (debug_flag){ + nhlbi_toolbox::utils::write_cpu_nd_array(base_gradients, this->debug_folder_ + std::string("base_gradients.real2")); + } + + int samples_per_interleave_ = base_gradients.get_number_of_elements(); if (spiral_rotations_ == 0) { // this is a hack which requires this parameter.. // normal operation + GDEBUG_STREAM("Using default spiral rotations: " << Nints_); base_gradients = nhlbi_toolbox::Spiral::create_rotations(base_gradients, Nints_); } else { + GDEBUG_STREAM("Using custom spiral rotations: " << spiral_rotations_ * this->acc); // Custom spiral rotations base_gradients = nhlbi_toolbox::Spiral::create_rotations(base_gradients, spiral_rotations_ * this->acc); } @@ -89,25 +99,27 @@ namespace Gadgetron auto trajectories = nhlbi_toolbox::Spiral::calculate_trajectories(base_gradients, sample_time, krmax_); auto weights = nhlbi_toolbox::Spiral::calculate_weights_Hoge(base_gradients, trajectories); - - filename = "/opt/data/gt_data/trajectories.real2"; - nhlbi_toolbox::utils::write_cpu_nd_array(trajectories, filename); - filename = "/opt/data/gt_data/weights.real"; + + if (debug_flag){ + nhlbi_toolbox::utils::write_cpu_nd_array(trajectories, this->debug_folder_ + std::string("trajectories.real2")); + nhlbi_toolbox::utils::write_cpu_nd_array(weights, this->debug_folder_ + std::string("weights.real")); - nhlbi_toolbox::utils::write_cpu_nd_array(weights, filename); + } if (this->girf_kernel) { - // base_gradients=Gadgetron::GIRF::girf_correct(base_gradients, this->girf_kernel, rotation_matrix, 2e-6, 10e-6, 0.85e-6); + // base_gradients=Gadgetron::GIRF::girf_correct(base_gradients, this->girf_kernel, rotation_matrix, 2e-6, 10e-6, this->clock_shift_s); base_gradients = correct_gradients(base_gradients, sample_time, this->girf_sampling_time_us, acq_header.read_dir, acq_header.phase_dir, acq_header.slice_dir); - auto filename = "/opt/data/gt_data/base_gradients_correct.real2"; - nhlbi_toolbox::utils::write_cpu_nd_array(base_gradients, filename); + if (debug_flag){ + nhlbi_toolbox::utils::write_cpu_nd_array(base_gradients, this->debug_folder_ + std::string("base_gradients_correct.real2")); + } // Weights should be calculated without GIRF corrections according to Hoge et al 2005 trajectories = nhlbi_toolbox::Spiral::calculate_trajectories(base_gradients, sample_time, krmax_); - - - filename = "/opt/data/gt_data/trajectories_correct.real2"; - nhlbi_toolbox::utils::write_cpu_nd_array(trajectories, filename); + + if (debug_flag){ + nhlbi_toolbox::utils::write_cpu_nd_array(trajectories, this->debug_folder_ + std::string("trajectories_correct.real2")); + } + weights = nhlbi_toolbox::Spiral::calculate_weights_Hoge(base_gradients, trajectories); } @@ -140,6 +152,15 @@ namespace Gadgetron this->acc = acc; } + void TrajectoryParameters_lit::set_debug_folder(std::string debug_folder) + { + this->debug_folder_ = debug_folder; + } + void TrajectoryParameters_lit::set_clock_shift(float shift_s) + { + this->clock_shift_s = shift_s; + } + TrajectoryParameters_lit::TrajectoryParameters_lit(const ISMRMRD::IsmrmrdHeader &h) { ISMRMRD::TrajectoryDescription traj_desc; @@ -202,6 +223,7 @@ namespace Gadgetron GDEBUG("gmax: %f\n", gmax_); GDEBUG("Tsamp_ns: %d\n", Tsamp_ns_); GDEBUG("Nints: %d\n", Nints_); + GDEBUG("spiral_rotation: %d\n", spiral_rotations_); GDEBUG("fov: %f\n", fov_); GDEBUG("krmax: %f\n", krmax_); GDEBUG("GIRF kernel: %d\n", bool(this->girf_kernel)); @@ -225,7 +247,7 @@ namespace Gadgetron rotation_matrix(1, 2) = slice_dir[1]; rotation_matrix(2, 2) = slice_dir[2]; - return nhlbi_toolbox::corrections::girf_correct(gradients, *girf_kernel, rotation_matrix, grad_samp_us, girf_samp_us, 0.85e-6); + return nhlbi_toolbox::corrections::girf_correct(gradients, *girf_kernel, rotation_matrix, grad_samp_us, girf_samp_us, this->clock_shift_s); } } // namespace Spiral } // namespace Gadgetron \ No newline at end of file diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/TrajectoryParameters_lit.h b/toolboxes/nhlbi_gt_toolbox/spiral/TrajectoryParameters_lit.h index a7a4cdb..13c9261 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/TrajectoryParameters_lit.h +++ b/toolboxes/nhlbi_gt_toolbox/spiral/TrajectoryParameters_lit.h @@ -30,7 +30,8 @@ namespace Gadgetron calculate_trajectories_and_weight(const ISMRMRD::AcquisitionHeader &acq_header); void set_girf_sampling_time(float time); void set_acceleration_factor(size_t acc); - + void set_debug_folder(std::string debug_folder); + void set_clock_shift(float shift_s); void read_girf_kernel(std::string girf_folder); hoNDArray> get_girf_kernel(); @@ -39,6 +40,8 @@ namespace Gadgetron double vds_factor_; // custom rotation number long spiral_rotations_; + + private: Core::optional>> girf_kernel; @@ -52,6 +55,8 @@ namespace Gadgetron float TE_; size_t acc; std::string systemModel; + std::string debug_folder_; + float clock_shift_s; hoNDArray correct_gradients(const hoNDArray &gradients, float grad_samp_us, float girf_samp_us, const float *read_dir, const float *phase_dir, diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.cpp b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.cpp index 2bd4d33..4359b13 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.cpp +++ b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.cpp @@ -293,6 +293,31 @@ cuNDArray noncartesian_reconstruction::crop_to_recondims(cuNDArray& inp return output; } +template +template +cuNDArray noncartesian_reconstruction::crop_to_recondims_centered(cuNDArray &input) +{ + cuNDArray output; //(this->recon_dims_reconSpace); + + if (input.get_number_of_dimensions() > 3){ + output.create(image_dims_[0], image_dims_[1], this->recon_dims_reconSpace[2], input.get_size(3)); + crop(uint64d4((input.get_size(0) - image_dims_[0]) / 2, (input.get_size(1) - image_dims_[1]) / 2, (image_dims_[2] - this->recon_dims_reconSpace[2]) / 2, 0), + uint64d4(image_dims_[0], image_dims_[1], this->recon_dims_reconSpace[2], input.get_size(3)), + input, + output); + } + else{ + output.create(image_dims_[0], image_dims_[1], this->recon_dims_reconSpace[2]); + crop(uint64d3((input.get_size(0) - image_dims_[0]) / 2, (input.get_size(1) - image_dims_[1]) / 2, (image_dims_[2] - this->recon_dims_reconSpace[2]) / 2), + //uint64d3(), + uint64d3(output.get_size(0), output.get_size(1), output.get_size(2)), + input, + output); + } + return output; +} + + template boost::shared_ptr> noncartesian_reconstruction::generateCSM(cuNDArray* channel_images) { @@ -1267,8 +1292,12 @@ template class noncartesian_reconstruction<2>; template class noncartesian_reconstruction<3>; template cuNDArray noncartesian_reconstruction<2>::crop_to_recondims(cuNDArray& input); template cuNDArray noncartesian_reconstruction<3>::crop_to_recondims(cuNDArray& input); +template cuNDArray noncartesian_reconstruction<2>::crop_to_recondims_centered(cuNDArray &input); +template cuNDArray noncartesian_reconstruction<3>::crop_to_recondims_centered(cuNDArray &input); template cuNDArray noncartesian_reconstruction<2>::crop_to_recondims(cuNDArray& input); template cuNDArray noncartesian_reconstruction<3>::crop_to_recondims(cuNDArray& input); +template cuNDArray noncartesian_reconstruction<2>::crop_to_recondims_centered(cuNDArray &input); +template cuNDArray noncartesian_reconstruction<3>::crop_to_recondims_centered(cuNDArray &input); template std::vector> noncartesian_reconstruction<2>::arraytovector(cuNDArray* inputArray, std::vector number_elements); template std::vector> diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.h b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.h index c8302dc..0d8997c 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.h +++ b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.h @@ -147,6 +147,9 @@ namespace nhlbi_toolbox template cuNDArray crop_to_recondims(cuNDArray &input); + + template + cuNDArray crop_to_recondims_centered(cuNDArray &input); boost::shared_ptr> nfft_plan_; std::vector image_dims_; diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/reconParams.h b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/reconParams.h index aa0d590..4b191b9 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/reconParams.h +++ b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/reconParams.h @@ -9,6 +9,7 @@ namespace Gadgetron ISMRMRD::MatrixSize ematrixSize; ISMRMRD::MatrixSize rmatrixSize; ISMRMRD::MatrixSize omatrixSize; + ISMRMRD::MatrixSize rmatrixSize_scanner; ISMRMRD::FieldOfView_mm fov; hoNDArray shots_per_time; size_t numberChannels; @@ -24,6 +25,7 @@ namespace Gadgetron float oversampling_factor_dcf_ = 2.1; int selectedDevice = 0; std::vector selectedDevices ; + std::vector selectedDevices_solver; float lambda_spatial = 1e-1; float lambda_spatial_imoco = 1e-1; float lambda_time = 1e-1; diff --git a/toolboxes/nhlbi_gt_toolbox/utils/CMakeLists.txt b/toolboxes/nhlbi_gt_toolbox/utils/CMakeLists.txt index 80b582c..c6df115 100644 --- a/toolboxes/nhlbi_gt_toolbox/utils/CMakeLists.txt +++ b/toolboxes/nhlbi_gt_toolbox/utils/CMakeLists.txt @@ -72,6 +72,12 @@ python/utils_function.py python/Pulseq_WaveformToTrajectory.py python/operators/low_rank.py python/kspace_preconditioning.py +python/GriddingTCRGadget.py +python/ROVirGadget.py +python/storage_server.py +python/tcr_utils.py +python/mrdhelper.py +python/SlicerGadget.py ) install(FILES ${gadgetron_python_gadgets_files} diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/GriddingTCRGadget.py b/toolboxes/nhlbi_gt_toolbox/utils/python/GriddingTCRGadget.py new file mode 100644 index 0000000..aa3a0e1 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/GriddingTCRGadget.py @@ -0,0 +1,538 @@ +import gadgetron +import numpy as np +import sigpy as sp +from sigpy.linop import NUFFT +import cupy as cp +import time +import ismrmrd as mrd +from sigpy.app import MaxEig +from sigpy.linop import FiniteDifference, Wavelet +import ctypes + +import matplotlib.pyplot as plt + +from storage_server import Storage +from tcr_utils import draw_number_indicators, compute_preemphasis_order, gram_schmidt, replace_zero, update_nonant, resize_with_crop_or_pad, analyticaldcf, remove_zero_padding, soft_threshold, crop_half_FOV, online_STCR_ISTA_2_timed, online_TCR_POGM_2, modified_gram_schmidt +from tcr_utils import * + +from skimage.util import montage +from skimage.transform import resize +import ismrmrd as mrd +from utils_function import eprint, parse_params, read_params + + +def pad_center_cupy(arr, new_shape, constant_values=0): + """ + Resizes a CuPy array by padding it from the center. + + This function takes a smaller CuPy array and places it in the center of a + new, larger array of a specified shape. The extra space is filled with a + constant value (zero by default). + + Args: + arr (cp.ndarray): The input CuPy array to be padded. + new_shape (tuple): The desired output shape for the new array. Each + dimension of new_shape must be greater than or equal + to the corresponding dimension in the input array's shape. + constant_values (scalar, optional): The value to use for padding. + Defaults to 0. + + Returns: + cp.ndarray: A new CuPy array with the specified shape, containing the + original array centered within it. + + Raises: + ValueError: If any dimension in `new_shape` is smaller than the + corresponding dimension of the input `arr`. + """ + # Get the shape of the input array + old_shape = arr.shape + + # Validate that the new shape is not smaller than the old shape + if any(new_dim < old_dim for new_dim, old_dim in zip(new_shape, old_shape)): + raise ValueError("Each dimension of new_shape must be greater than or equal to the old shape.") + + # Calculate the padding widths for each dimensfion. + # The padding is distributed as evenly as possible on both sides. + # If the total padding for a dimension is odd, the extra padding + # element is added to the end. + pad_width = [] + for new_dim, old_dim in zip(new_shape, old_shape): + # Total padding needed for the current dimension + delta = new_dim - old_dim + # Padding before the array + pad_before = delta // 2 + # Padding after the array + pad_after = delta - pad_before + pad_width.append((pad_before, pad_after)) + + # Use cupy.pad to apply the calculated padding + padded_arr = cp.pad(arr, pad_width, mode='constant', constant_values=constant_values) + + return padded_arr + + +def create_ismrmrd_image(data, field_of_view, index, acquisition): + + data = 4000 * data / np.max(data) + + im = mrd.image.Image.from_array( + np.float32(data), + image_series_index=index, + acquisition=acquisition, + image_type=mrd.IMTYPE_MAGNITUDE, + field_of_view=(field_of_view.x, field_of_view.y, field_of_view.z), + transpose=True + ) + + return im + +# RR add in for iMRI +def send_imri_image(data, mrdHeader,connection,acq_data): + # imagesOut = [None] * data.shape[-1] + + # Determine max value (12 or 16 bit) + BitsStored = 12 + # fix this: + # if (mrdhelper.get_userParameterLong_value(mrdHeader, "BitsStored") is not None): + # BitsStored = mrdhelper.get_userParameterLong_value(mrdHeader, "BitsStored") + maxVal = 2**BitsStored - 1 + + #data = abs(data) + #data = np.flip(data,axis=2) # flip for scanner? + + + + data = np.abs(data) + data *= maxVal / data.max() + data = np.flip(data,axis=2) # flip for scanner? + data = np.around(data) + + eprint(np.max(np.abs(data))) + eprint(np.min(data)) + + # add logic for over-sampled slices + #mrdHeader.encoding[0].encodedSpace.matrixSize.z + #mrdHeader.encoding[0].encodreconSpaceedSpace.matrixSize.z + + for partition in range(data.shape[-1]): + # Create new MRD instance for the processed image + #imagesOut = mrd.Image.from_array(data[...,partition], transpose=False) + imagesOut = mrd.image.Image.from_array( + np.float32(np.abs(data[:,:,partition])), + acquisition=acq_data, + image_series_index=0, + repetition=acq_data._head.idx.repetition, #=0, + image_type=mrd.IMTYPE_MAGNITUDE, + transpose=True + ) + + #data.data = data.data[:,partition,:,:] + + # Set the header information + # imagesOut.setHead(mrdhelper.update_img_header_from_raw(imagesOut.getHead(), rawHead[partition])) + # check everything has been assigned + + res = mrdHeader.encoding[0].reconSpace.fieldOfView_mm.x/mrdHeader.encoding[0].reconSpace.matrixSize.x + fovx = data.shape[0] * res + + imagesOut.field_of_view = (ctypes.c_float(fovx), + ctypes.c_float(fovx), + ctypes.c_float(mrdHeader.encoding[0].reconSpace.fieldOfView_mm.z/int(mrdHeader.encoding[0].encodedSpace.matrixSize.z))) + + imagesOut.slice = partition + + # Positions need to be set for the scanner to be happy + # ------------------------------------- + + #par_thickness = mrdHeader.encoding[0].encodedSpace.fieldOfView_mm.z/mrdHeader.encoding[0].encodedSpace.matrixSize.z + + par_thickness = mrdHeader.encoding[0].reconSpace.fieldOfView_mm.z/mrdHeader.encoding[0].reconSpace.matrixSize.z + + + # determine rotation matrix in [PRS] - stacking per lit-python + R_matrix = np.hstack((np.array(imagesOut.phase_dir).reshape(3,1),np.array(imagesOut.read_dir).reshape(3,1),np.array(imagesOut.slice_dir).reshape(3,1))) + + # partition is "third" dimension in matric coordinates + #partition_vector = np.array([0,0,(partition-data.shape[-1]/2.0)*par_thickness]) + partition_vector = np.array([0,0,(partition-data.shape[-1]/2.0)*par_thickness + par_thickness*3.5]) # testing slice offset, adding fudge for slice matching to Siemens + + + # multiply to determine appropriate offset position + position_offset = np.matmul(R_matrix,partition_vector) + for ii in range(3): # mrd wants c-float-array-3, and I dont know how to do a vector cast.. + # test casting # imagesOut.position[ii] = ctypes.c_float(imagesOut.position[ii] + position_offset[ii]) # added to test for imri + imagesOut.position[ii] = imagesOut.position[ii] + position_offset[ii] # added to test for imri + + #eprint(np.array(imagesOut.position)) # you will regret this + # ------------------------------------- + + + # Set ISMRMRD Meta Attributes + # ------------------------------------- + # An update to FIRE might fix a few things, but fixing some meta data manually + # note this only works with Fang's custom iMRI which helps us out + # Make overlay comment (need to get SOP instance?) + meta = mrd.Meta.deserialize(imagesOut.attribute_string) + #meta['GADGETRON_ImageComment'] = 'Overlay' + meta['DataRole'] = 'Image' + meta['ImageProcessingHistory'] = ['FIRE', 'PYTHON'] + meta['WindowCenter'] = str((maxVal + 1) / 2) + meta['DataRole'] = str((maxVal + 1)) + + # meta = ismrmrd.Meta({'DataRole': 'Image', + # 'ImageProcessingHistory': ['FIRE', 'PYTHON'], + # 'WindowCenter': str((maxVal + 1) / 2), + # 'WindowCenter': str((maxVal + 1))}) + + # Add image orientation directions to MetaAttributes if not already present + if meta.get('ImageRowDir') is None: + meta['ImageRowDir'] = ["{:.18f}".format(imagesOut.getHead().read_dir[0]), + "{:.18f}".format(imagesOut.getHead().read_dir[1]), + "{:.18f}".format(imagesOut.getHead().read_dir[2])] + + if meta.get('ImageColumnDir') is None: + meta['ImageColumnDir'] = ["{:.18f}".format(imagesOut.getHead().phase_dir[0]), + "{:.18f}".format(imagesOut.getHead().phase_dir[1]), + "{:.18f}".format(imagesOut.getHead().phase_dir[2])] + + meta['DimBoundariesIdc'] = 1 + meta['DimBoundariesIdd'] = 1 + meta['DimBoundariesIde'] = 1 + meta['DimBoundariesSeg'] = 1 + + metaXml = meta.serialize() + imagesOut.attribute_string = metaXml + + # Send device overlay + connection.send(imagesOut) + + +def _parse_params(xml): + return {p.get('name'): p.get('value') for p in xml.iter('property')} + +def GriddingTCRGadget(connection): + if (connection.header.subjectInformation is not None) and (connection.header.subjectInformation.patientID is not None): + subj_str = connection.header.subjectInformation.patientID + eprint(subj_str) + else: + subj_str = connection.header.measurementInformation.measurementID.split("_")[1] + eprint(subj_str) + + #params = _parse_params(connection.config) + params_init = parse_params(connection.config) + params={'n_tr_bin':int(216), ## GRIDDING parameters + 'montage':False, ## MONTAGE parameters + 'crop_percent': 0.6, + 'num_slices':int(1), + 'use_storage_server':False, + 'alg_type':"gridding", + 'niter':0, ## STCR Parameters + "lambdat":0, ## STCR Parameters + "lambdas":0, ## STCR Parameters + "n_rovir_coils":int(0),## ROVir Parameters + "recompute_dcf": False,## Analytical DCF + "recompute_L":False, + "time_per_frame_scale":0.7, + "max_frame_discard":int(6), + 'storage_port':int(9112), + } + + boolean_keys=['montage','use_storage_server','recompute_dcf','recompute_L'] + str_keys=['alg_type'] + int_keys=['n_tr_bin','num_slices','niter','n_rovir_coils','max_frame_discard','storage_port'] + float_keys=['crop_percent','lambdat','lambdas',"time_per_frame_scale"] + #eprint(f"storage port {params['storage_port']} {type(params['storage_port'])}") + #eprint(f"storage port {params_init['storage_port']} {type(params_init['storage_port'])}") + params=read_params(params_init,params_ref=params,boolean_keys=boolean_keys,str_keys=str_keys,int_keys=int_keys,float_keys=float_keys) + #eprint(f"storage port {params['storage_port']} {type(params['storage_port'])}") + + n_tr_bin = params['n_tr_bin'] + bool_montage=params["montage"] + crop_percent = params["crop_percent"] + n_slices_extract =params["num_slices"] + use_storage_server=params["use_storage_server"] + alg_type = params["alg_type"] + niter = params["niter"] + lambdat = params["lambdat"] + lambdas = params["lambdas"] + n_rovir_coils=params["n_rovir_coils"] + recompute_dcf=params["recompute_dcf"] + recompute_L=params["recompute_L"] + time_per_frame_scale=params["time_per_frame_scale"] + MAX_FRAME_DISCARD=params["max_frame_discard"] + storage_port=int(params['storage_port']) + eprint(f"storage port {storage_port} {type(storage_port)}") + # simply write frame_discard to a location for the SlicerGadget + np.save("frame_discard.npy", np.array(MAX_FRAME_DISCARD)) + + center_shift = 0 + + storage = Storage("localhost", storage_port, subject=subj_str) + storage_vars = Storage("localhost", storage_port) + + # variable initializations + L = None + n_frame_discard = 0 + n_frame_discard = MAX_FRAME_DISCARD + csm_loaded = False + n_connections = 0 + kspace_buff = [] + head_buff = [] + first_frame = True + n_images_sent = 0 + was_lip_computed = False + recomputed_dcf = False + + maximum_scale = None + + for data in connection: + st = time.time() + + # for now, very dummy buffer for acquiring the n_tr_bin. + # this can be optimized to not have to wait for an entire volume, but a simple case for now! + if n_tr_bin == 0: + if data.is_flag_set(mrd.ACQ_LAST_IN_SLICE): + n_tr_bin = data.scan_counter + eprint(f"n_tr_bin: {n_tr_bin}") + else: + continue + + if csm_loaded is False: + # setup all the beginning variables! + csm_loaded = True + csm_str = f"{connection.header.encoding[0].encodedSpace.matrixSize.x}_{connection.header.encoding[0].encodedSpace.matrixSize.y}_{connection.header.encoding[0].encodedSpace.matrixSize.z}_{data._head.active_channels}" + + eprint(f"CSM STR: {csm_str}") + + # PK HACK OF EPIC PROPORTIONS + #csm_str = "112_112_48_21" + + csm = storage.fetch_latest(name=f"csm_{csm_str}") + if n_rovir_coils > 0: + W_d = storage.fetch_latest(f"ROVIR_W_{csm_str}") + W = W_d['W'] + #csm_rovir = W_d['csm_rovir'][0:n_rovir_coils,:,:,:] + W, th = modified_gram_schmidt(cp.array(W[:,0:n_rovir_coils])) + csm_rovir = rovir_apply(cp.array(csm.transpose((1,2,3,0))), cp.array(W), n_rovir_coils).transpose((3,0,1,2)) + csm = crop_half_FOV(csm_rovir, dims=(2,3), size=np.flip(csm.shape[1:])) + # grahm schmidt orthogonalization + #Q, R = np.linalg.qr(W[:,0:n_rovir_coils]) + #W = cp.array(Q) + #csm = (csm.transpose((1,2,3,0)) @ W).transpose((3,0,1,2)) + + #csm = remove_zero_padding(csm, 1) # gt gridding artifact + #csm = cp.flip(csm, axis=1) + csm = cp.asarray(csm).transpose((0,3,2,1)) + n_coils = csm.shape[0] + + traj_complex = storage.fetch_latest("trajectory") + traj = np.stack((traj_complex.real, traj_complex.imag), axis=-1) + dcf = cp.array(storage.fetch_latest(name=f"dcf_{csm_str}")) + + nread = traj.shape[1]*n_tr_bin + dcf_sqrt = cp.sqrt(dcf[0:nread]) + What = sp.linop.Multiply([csm.shape[0], nread], dcf_sqrt) + S = sp.linop.Multiply(csm.shape, csm, conj=True) + R = sp.linop.Sum(csm.shape, [0]) + kspace_buff = cp.zeros((data._head.active_channels, n_tr_bin, traj.shape[1]), dtype=cp.complex64) + + csm_image = storage.fetch_latest(name=f"csm_image_{csm_str}") + # csm_image = remove_zero_padding(csm_image, 1) # gt gridding artifact + csm_image = cp.array(csm_image).transpose((2,1,0)) + #csm_image = crop_half_FOV(csm_image, dims=[0,1]) + + G = FiniteDifference(csm_image.shape) + + # try to pull the pre-computed L + L_scale_csm_scale_str = connection.header.measurementInformation.protocolName + "_len_" + str(n_tr_bin) + try: + eprint(f"L_scale_csm: {L_scale_csm_scale_str}") + storage_scale_fetch = storage.fetch_latest(custom_tags={"scale_factor_python":L_scale_csm_scale_str}) + L = storage_scale_fetch['L'] + eprint(f"FETCHED CSM L: {L}") + except: + eprint("failed to fetch storage server. Will re-generate....") + if recompute_L: + L = None + + eprint(f"TIME SPENT IN FIRST FRAME: {time.time() - st}") + + TR = connection.header.sequenceParameters.TR + time_per_frame = n_tr_bin * TR[0] + + kspace_buff[:,n_connections,:] = cp.array(data.data) + head_buff.append(data._head) + + # bin up the data for gridding + if n_connections == (n_tr_bin-1): + n_connections = 0 + encode_step_1 = [head_buff[i].idx.kspace_encode_step_1 for i in range(n_tr_bin)] + encode_step_2 = [head_buff[i].idx.kspace_encode_step_2 for i in range(n_tr_bin)] + head_buff = [] + n_slice = csm.shape[3] + + trajectory_frame = np.zeros((n_tr_bin, traj.shape[1], 3)) + trajectory_frame[:,:,0:2] = traj[encode_step_1] * csm.shape[1] + + trajectory_frame[:,:,2] = ((np.repeat(encode_step_2, traj.shape[1]).reshape((n_tr_bin,traj.shape[1])) / n_slice) - 0.5) * n_slice + trajectory_frame = trajectory_frame.reshape((n_tr_bin * traj.shape[1], 3)) + trajectory_frame = cp.array(trajectory_frame) + + if recompute_dcf and not recomputed_dcf: + eprint("recomputing the DCF...") + dcf = analyticaldcf(trajectory_frame.get(), ns=nread) + dcf_sqrt = cp.sqrt(cp.array(dcf)) + What = sp.linop.Multiply([csm.shape[0], nread], dcf_sqrt) + recomputed_dcf = True + + F = NUFFT(csm.shape, trajectory_frame, toeplitz=True) + #Ah = R * S * F.H * What + Aframe = What * F * S.H * R.H + + # apply rovir if necessary + if n_rovir_coils > 0: + kspace_rov = rovir_apply(kspace_buff.transpose(1,2,0), cp.array(W), n_rovir_coils).transpose(2,0,1) + y = (kspace_rov.reshape(kspace_rov.shape[0], -1) * dcf_sqrt) + image = Aframe.H * y + n_coils = n_rovir_coils + else: + y = (kspace_buff.reshape(kspace_buff.shape[0], -1) * dcf_sqrt) + image = Aframe.H * y + # precompute L if necessary + if L is None: + L = MaxEig(Aframe.N, max_iter=40, dtype=image.dtype, device=image.device, show_pbar=True).run() + eprint("Lipschitz Constant: {}".format(L)) + n_frame_discard = MAX_FRAME_DISCARD + was_lip_computed = True + storage.store({"L": L}, custom_tags={"scale_factor_python": L_scale_csm_scale_str}) + cp.cuda.stream.get_current_stream().synchronize() + eprint(f" time spent GRIDDING {time.time() - st}") + + def cost_stcr(x): + return 0.5 * cp.square(cp.linalg.norm((Aframe * (x + xn_1) - y).flatten())) + \ + (lambdat * cp.linalg.norm(x.flatten(), ord=1)) + (lambdas * cp.linalg.norm((G * x).flatten(), ord=1)) + + def cost_tcr(x): + return 0.5 * cp.square(cp.linalg.norm((Aframe * (x + xn_1) - y).flatten())) + \ + (lambdat * cp.linalg.norm(x.flatten(), ord=1)) + + stcrt = time.time() + # do STCR, or skip a few frames if needed to make up for L computation + if niter > 0 and n_frame_discard == 0: + if first_frame or (was_lip_computed and n_images_sent == (MAX_FRAME_DISCARD)): + max_image = cp.abs(image).max() + scale_factor = cp.abs(csm_image).max() / max_image + del_0 = csm_image / scale_factor + xn = csm_image / scale_factor + else: + max_image = cp.abs(image).max() + scale_factor = cp.abs(csm_image).max() / max_image + + lamt = lambdat * max_image + lams = lambdas * max_image + + if alg_type == "cg": + cg_alg = sp.alg.ConjugateGradient(Aframe.H*Aframe, image, xn_1, max_iter=niter) + while not cg_alg.done(): + cg_alg.update() + xn = cg_alg.x + elif alg_type == "stcr": + del_0 = (image - xn_1) / 4 # soft initialization + eprint(time_per_frame) + del_0 = online_STCR_ISTA_2_timed(Aframe, G, xn_1, image, lamt, lams, 1/L, mu=0.1, yn=(kspace_buff.reshape(kspace_buff.shape[0], -1) * dcf_sqrt), time_recon=time_per_frame*time_per_frame_scale, deln=del_0) + xn = xn_1 + del_0 + elif alg_type == "tcr": + del_0 = cp.zeros(image.shape) + del_0 = online_TCR_POGM_2(Aframe, xn_1, image, lamt, (1/L), n_iter=niter, del_0=del_0) + xn = xn_1 + del_0 + elif alg_type == "stcr_l1_approx": + xn = xn_1 + nesterov_dict = None + for it in range(niter): + delta = 1 + #xn = gradient_descent_iteration(xn, lambda d : grad_cost(d, xn_1, Aframe, y, G, delta, lamt, lams), 1/(L+(2/delta))) + xn, nesterov_dict = gradient_descent_iteration_nesterov(it, xn, lambda d : grad_cost(d, xn_1, Aframe, y, G, delta, lamt, lams), 1/(L+(2/delta)), max_iter=niter, nesterov_dict=nesterov_dict) + elif alg_type == "gridding": + xn = image + + eprint(f" ALG TYPE: {alg_type}") + else: + n_frame_discard = max(0, n_frame_discard-1) + xn = image + + eprint(f" RECON MATRIX SIZE: {xn.shape}") + # setup image out + xn_1 = cp.copy(xn) + + # do a fourier interpolation to 2x the matrix size + #xn = cp.fft.fftshift(cp.fft.fftn(xn, axes=(0,1,2)), axes=(0,1,2)) + #xn = pad_center_cupy(xn, [xn.shape[0]*2,xn.shape[1]*2,xn.shape[2]*2]) + #xn = cp.fft.ifftn(cp.fft.ifftshift(xn, axes=(0,1,2)), axes=(0,1,2)) + + # do we want to montage? + #if bool_montage is True: + if False: + if use_storage_server: + montage_info = storage_vars.fetch_latest(custom_tags={"montage":"1"}) + if 'center_shift' in montage_info: + center_shift = int(montage_info['center_shift']) + n_slices_extract = int(montage_info['n_slices']) + # we have images x, y, f. + image_for_m = crop_middle(extract_middle_slices(cp.abs(xn_1).squeeze(), n_slices_extract, 2, center_shift), [0,1], crop_percent) + image_for_m = draw_number_indicators(image_for_m, indicator_size=8) + image_for_m = image_for_m[:,:,compute_preemphasis_order(image_for_m.shape[2])] + image_m = montage(image_for_m.transpose((2,0,1)).get()) + elif 'oreintation_list' in montage_info: + # do something else + # list is of type ['a21', 'a22', 'a23', 'a3', 'a4', 'a5', 'a6', 'a7', 'a8'] + # number is index, and letter is a for axial, s for sagittal, c for coronal. + orientations = montage_info['oreintation_list'] + image_for_m = [] + for nonant_str in orientations: + image_for_m.append(update_nonant(xn, nonant_str).get()) + image_for_m = np.array(image_for_m) + image_for_m[:,:,compute_preemphasis_order(image_for_m.shape[0])] + image_m = montage(image_for_m) + eprint(f"{montage_info}") + else: + image_m = xn_1.get() + + # if it's the first image we output, we need the size to always be the same. So do a re-size if needed to match the first frame size. + if first_frame: + first_shape = image_m.shape + first_frame = False + else: + if not image_m.shape == first_shape: + eprint("resizing...") + image_m = resize(image_m, first_shape) + + # setup image and send. + n_images_sent = n_images_sent + 1 + if n_images_sent > MAX_FRAME_DISCARD and maximum_scale is None: + maximum_scale = np.abs(image_m).max() * 0.7 + if maximum_scale is not None: + # everywhere the image > maximum_scale, set to maximum_scale (but keep the phase) + #image_m = np.where(np.abs(image_m) > maximum_scale, maximum_scale * np.exp(1j * np.angle(image_m)), image_m) + pass + #out_im = create_ismrmrd_image(np.abs(image_m), connection.header.encoding[0].encodedSpace.fieldOfView_mm, 0, data) + cp.cuda.stream.get_current_stream().synchronize() + eprint(f"PYTHON TOTAL ELAPSED TIME: {time.time() - st}") + + if bool_montage: + out_im = create_ismrmrd_image(np.abs(image_m), connection.header.encoding[0].encodedSpace.fieldOfView_mm, 0, data) + connection.send(out_im) + else: + # we could, instead of sending gridding images, just ignore them. + if n_images_sent > MAX_FRAME_DISCARD: + eprint("send_imri_image called") + send_imri_image(image_m, connection.header,connection, data) + else: + n_connections = n_connections + 1 + + # cleaup + storage.close() + storage_vars.close() + +if __name__ == '__main__': + gadgetron.external.listen(2000, GriddingTCRGadget) \ No newline at end of file diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/ROVirGadget.py b/toolboxes/nhlbi_gt_toolbox/utils/python/ROVirGadget.py new file mode 100644 index 0000000..00601b4 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/ROVirGadget.py @@ -0,0 +1,410 @@ +import gadgetron +import ismrmrd as mrd + +# non-debugging +#from utils.python.storage_server import Storage +#from utils.python.tcr_utils import rovir_apply, crop_half_FOV +#from utils.python.coil_maps import calculate_csm_walsh_gpu + +#debugging +from storage_server import Storage +from tcr_utils import * +from coil_maps import calculate_csm_walsh_gpu + +import numpy as np +from scipy.ndimage import gaussian_filter +from scipy.linalg import eig +import cupy as cp + +from sigpy import fft +from sigpy.mri.app import EspiritCalib + +from skimage.morphology import convex_hull_image + +from utils_function import eprint, load_json, parse_params, read_params + +def automatic_mask_3D(image, threshold_ratio=0.1): + """ + Generates a 3D mask by thresholding the image and computing the convex hull + for each slice along the z-axis (axis 2). + + Parameters: + ----------- + image : numpy.ndarray + A 3D numpy array (x, y, z) representing the volume. + threshold_ratio : float + The fraction of the maximum intensity to use as a threshold + (default is 0.1, i.e., 10%). + + Returns: + -------- + mask : numpy.ndarray + A boolean 3D numpy array of the same shape as input image, + where True indicates the mask region. + """ + + # 1. Basic validation + if image.ndim != 3: + raise ValueError("Input image must be a 3D array.") + + # 2. Determine Signal Region (Thresholding) + # We use the max value of the entire 3D volume + max_val = np.max(image) + threshold = max_val * threshold_ratio + + # Create initial binary mask + # We treat any value > threshold as signal + signal_mask = image > threshold + + # 3. Convex Hull per Slice + # We assume 'z' is the last axis (axis 2). + # If your data is (z, x, y), change the iteration accordingly. + + final_mask = np.zeros_like(signal_mask, dtype=bool) + + depth = image.shape[2] # z-axis + + for z in range(depth): + slice_2d = signal_mask[:, :, z] + + # convex_hull_image requires at least one True value to work. + # If the slice is empty (no signal), the mask for this slice remains False. + if np.any(slice_2d): + # Compute convex hull for this slice + hull_2d = convex_hull_image(slice_2d) + final_mask[:, :, z] = hull_2d + else: + # Keep as zeros/False + final_mask[:, :, z] = slice_2d + + # Optional: Fill holes in 3D if needed, though slice-by-slice hull usually suffices + # for "wrapping" the object. + # To make it a solid 3D block if there are gaps between slices, + # one might perform binary_closing or binary_fill_holes in 3D, + # but the requirement was specifically "convex hull in each slice". + return split_mask_by_extent(final_mask, central_ratio=0.7, axis=0) + +def split_mask_by_extent(mask, central_ratio=0.8, axis=2): + """ + Splits a 3D mask into two regions: a central region and an off-region + (upper and lower tips), based on the extent of the mask along the specified axis. + + Parameters: + ----------- + mask : numpy.ndarray + A 3D boolean numpy array. + central_ratio : float + The proportion of the extent to keep in the central mask. + axis : int + The dimension along which to calculate extent and split. + + Returns: + -------- + central_mask : numpy.ndarray + off_mask : numpy.ndarray + """ + if mask.ndim != 3: + raise ValueError("Input mask must be a 3D array.") + + # Move target axis to 0 for easy processing + # If shape is (X, Y, Z) and axis is 2, shape becomes (Z, X, Y) + mask_moved = np.moveaxis(mask, axis, 0) + + # 1. Find indices along axis 0 where mask has positive values + # mask_moved has shape (Depth, H, W). We collapse (1, 2) + presence = np.any(mask_moved, axis=(1, 2)) + indices = np.where(presence)[0] + + central_moved = np.zeros_like(mask_moved) + off_moved = np.zeros_like(mask_moved) + + # Handle empty mask case + if len(indices) == 0: + return np.zeros_like(mask), np.zeros_like(mask) + + # 2. Determine extent + idx_min = indices[0] + idx_max = indices[-1] + + length = idx_max - idx_min + margin = length * (1.0 - central_ratio) / 2.0 + + cutoff_bottom = idx_min + margin + cutoff_top = idx_max - margin + + # 3. Create the masks + depth = mask_moved.shape[0] + + for i in range(depth): + if presence[i]: + if cutoff_bottom <= i <= cutoff_top: + central_moved[i, ...] = mask_moved[i, ...] + else: + off_moved[i, ...] = mask_moved[i, ...] + + # 4. Move axis back to original position + central_mask = np.moveaxis(central_moved, 0, axis) + off_mask = np.moveaxis(off_moved, 0, axis) + + return central_mask, off_mask + +def get_ES_mask_3D(im3d, mode='simple', splitud=False, scale_amt=2, margin=10): + """ + Generates a mask for a 3D image. + + Parameters: + - im3d: 3D numpy array (the input image). + - mode: 'simple' will not use thresholding (large rectangles), 'advanced' will use thresholding (smaller regions). + - splitud: bool, whether to create two masks (default False). + - scale_amt: float, intensity scaling factor (default 2). + + Returns: + - mask: 3D numpy array with the mask applied. + """ + matrix_size = im3d.shape + ny, nx, n_slice = matrix_size + + mask = np.zeros(matrix_size) + + if mode == 'simple': + mask_top = np.zeros(matrix_size) + mask_bot = np.zeros(matrix_size) + mask_fov = np.zeros(matrix_size) + + mask_top[:int(ny // 2.5), :, :] = 1 + mask_bot[3*ny // 4:, :, :] = 1 + + mask = mask_top + mask_bot + + # mask_fov is the center 40% of image + #mask_fov[ny // 3:2 * ny // 3, nx // 3:2 * nx // 3, (n_slice//6):(n_slice*5//6)] = 1 + mask_fov[ny // 5:4 * ny // 5, nx // 5:4 * nx // 5, (n_slice//6):(n_slice*5//6)] = 1 + + return (mask, mask_fov) + else: + # Gaussian filter + im_filt = gaussian_filter(im3d, sigma=2) + + # Split the filtered image if needed + if splitud: + im_filt_top = im_filt.copy() + im_filt_bot = im_filt.copy() + + im_filt_top[ny // 2:] = 0 + im_filt_bot[:ny // 2] = 0 + + im_split = [im_filt_top, im_filt_bot] + else: + im_split = [im_filt] + + for im_tmp in im_split: + M = np.zeros(n_slice) + row = np.zeros(n_slice, dtype=int) + col = np.zeros(n_slice, dtype=int) + + for slc in range(n_slice): + # Find the max value and its index for each slice + max_val = im_tmp[:, :, slc].max() + M[slc] = max_val + row[slc], col[slc] = np.unravel_index(np.argmax(im_tmp[:, :, slc]), im_tmp[:, :, slc].shape) + + win_size = 40 + + for slc in range(n_slice): + row_bounds = [np.clip(row[slc] - win_size, 0, nx), + np.clip(row[slc] + win_size, 0, nx)] + col_bounds = [np.clip(col[slc] - win_size, 0, ny), + np.clip(col[slc] + win_size, 0, ny)] + + mask1 = np.zeros(im_tmp.shape[:2]) + mask1[row_bounds[0]:row_bounds[1], col_bounds[0]:col_bounds[1]] = 1 + image_slice = im_tmp[:, :, slc] + mask2 = np.zeros_like(mask1) + mask2[(image_slice * mask1) > (M[slc] / scale_amt)] = 1 + mask[:, :, slc] += mask2 + + # Final FOV check: don't select inside FOV + margin = margin + fov_bound2dy = [ny // 2 - (ny // margin), ny // 2 + (ny // margin)] + fov_bound2dx = [nx // 2 - (nx // margin), nx // 2 + (nx // margin)] + + fov_mask = np.zeros((nx, ny)) + fov_mask[fov_bound2dy[0]:fov_bound2dy[1], fov_bound2dx[0]:fov_bound2dx[1]] = 1 + + fov_mask = fov_mask[:,:,np.newaxis] + mask = np.multiply(mask, (1-fov_mask)) + + fov_mask = np.ones(mask.shape) * fov_mask + + return (mask, fov_mask) + +def rovir_automatic_3d_pk_2(img_cc, csm_image_mc_, margin=5): + [rovir_sig_, rovir_int_] = automatic_mask_3D(np.abs(img_cc)) + nx, ny, nz = img_cc.shape + n_coil = csm_image_mc_.shape[-1] + rovir_sig_mask = rovir_sig_[:,:,:,np.newaxis] + rovir_int_mask = rovir_int_[:,:,:,np.newaxis] + signals = (rovir_sig_mask * csm_image_mc_).reshape(-1, n_coil) + interferences = (rovir_int_mask * csm_image_mc_).reshape(-1, n_coil) + + A = np.conj(signals.T) @ signals + B = np.conj(interferences.T) @ interferences + + # solve for generalized eigenvalues + eigvals, eigvecs = eig(B, A) + + # Sort by eigenvalues in descending order + idx = np.argsort(np.abs(eigvals))[::-1] + + # Normalize eigenvectors + V = eigvecs[:, np.flip(idx)] + + #V = V / np.linalg.norm(V, axis=0) + #np.save("W_rovir.npy", V) + + return [V, rovir_sig_, rovir_int_] + +def rovir_automatic_3d_pk(img, img_coils, margin=10): + """ + Automatically selects the signal region from a 3D image. + + Parameters: + - img: 3D numpy array (the main image data). + - img_coils: 4D numpy array (coil images). + + Returns: + - W: The ordered eigenvector matrix. + """ + [rovir_int_mask, rovir_sig_mask] = get_ES_mask_3D(np.abs(img), splitud=True, scale_amt=3, margin=margin) + + # Automatically select signal region + nx, ny, nz = img.shape + margin = 5 + margin_z = 3 + + """ + fov_bound2dy = [int(np.floor((ny / 2) - (ny / margin))), int(np.floor((ny / 2) + (ny / margin)))] + fov_bound2dx = [int(np.floor((nx / 2) - (nx / margin))), int(np.floor((nx / 2) + (nx / margin)))] + fov_bound2dz = [int(np.floor((nz / 2) - (nz / margin_z))), int(np.floor((nz / 2) + (nz / margin_z)))] + + rovir_sig_mask[fov_bound2dy[0]:fov_bound2dy[1], + fov_bound2dx[0]:fov_bound2dx[1], + fov_bound2dz[0]:fov_bound2dz[1]] = 1 + """ + + n_coil = img_coils.shape[3] + rovir_sig_mask = rovir_sig_mask[:,:,:,np.newaxis] + rovir_int_mask = rovir_int_mask[:,:,:,np.newaxis] + signals = (rovir_sig_mask * img_coils).reshape(-1, n_coil) + interferences = (rovir_int_mask * img_coils).reshape(-1, n_coil) + + A = signals.T @ signals + B = interferences.T @ interferences + + # Solve for generalized eigenvalues + eigvals, eigvecs = eig(B, A) + # np.linalg.eig(np.linalg.pinv(A) @ B) + + # Sort by eigenvalues in descending order + idx = np.argsort(np.abs(eigvals))[::-1] + + # Normalize eigenvectors + V = eigvecs[:, np.flip(idx)] + V = V / np.linalg.norm(V, axis=0) + + # Reorder W + W = V + + #np.save("W_rovir.npy", W) + + return W + +def remove_zero_padding(arr, axis): + """ + Remove zero padding along a specified axis in a NumPy array. + + Parameters: + arr (np.ndarray): Input array with potential zero padding. + axis (int): Axis along which to remove zero padding. + + Returns: + np.ndarray: Array with zero padding removed along the specified axis. + """ + if axis < 0 or axis >= arr.ndim: + raise ValueError(f"Invalid axis {axis} for array with {arr.ndim} dimensions.") + + # Create a boolean mask where True represents non-zero slices + non_zero_mask = np.any(arr != 0, axis=tuple(i for i in range(arr.ndim) if i != axis)) + + # Use the boolean mask to slice the array + slices = [slice(None)] * arr.ndim + slices[axis] = slice(None) # Keep the full range for the axis we are slicing + arr_sliced = arr[tuple(slices)] + + # Now use the boolean mask to filter out zero-only slices + arr_trimmed = arr_sliced.take(np.where(non_zero_mask)[0], axis=axis) + + return arr_trimmed + + +def ROVirGadget(connection): + if connection.header.subjectInformation.patientID is not None: + subj_str = connection.header.subjectInformation.patientID + else: + subj_str = connection.header.measurementInformation.measurementID.split("_")[1] + + params_init = parse_params(connection.config) + params={'storage_port':9112, + } + int_keys=['storage_port'] + + params=read_params(params_init,params_ref=params,int_keys=int_keys) + storage = Storage("localhost", params["storage_port"], subject=subj_str) + + setup = None + bool_espirit = True + + for data in connection: + if setup is None: + setup = True + size_path = f"{connection.header.encoding[0].encodedSpace.matrixSize.x}_{connection.header.encoding[0].encodedSpace.matrixSize.y}_{connection.header.encoding[0].encodedSpace.matrixSize.z}_{data.acq_headers[0].active_channels}" + csm_image_mc = storage.fetch_latest(name=f"csm_image_mc_{size_path}") + csm_over = storage.fetch_latest(name=f"csm_over_{size_path}") + img_cc = remove_zero_padding(np.sum(np.conj(csm_over) * csm_image_mc,axis=0),0).transpose((2, 1, 0)) + csm_image_mc_ = remove_zero_padding(csm_image_mc, 1).transpose((3,2,1,0)) + + if bool_espirit: + ksp = fft(remove_zero_padding(csm_image_mc, 1), axes=[-1, -2, -3]) + ksp = cp.array(ksp) + mps_rec = EspiritCalib(ksp, show_pbar=True, thresh=0.02, crop=0, device=cp.cuda.Device(device=0)).run() + + matrix_keep = int(connection.header.encoding[0].encodedSpace.matrixSize.x*1.5) + mps_rec = crop_half_FOV(mps_rec, dims=(-1, -2), size=[matrix_keep,matrix_keep]) + storage.store(mps_rec, name=f"csm_{size_path}") + + [W, rovir_sig_, rovir_int_] = rovir_automatic_3d_pk_2(img_cc, csm_image_mc_) + + if True: + rov_ims = rovir_apply(cp.array(csm_image_mc.transpose(1,2,3,0)), cp.array(W), 10).transpose(3,0,1,2) + # calculate sensitivity maps slice by slice. + rov_ims = cp.array(remove_zero_padding(rov_ims, 1)) + csm_rovir = np.zeros(rov_ims.shape, dtype=np.complex64) + + if bool_espirit: + ksp = fft(remove_zero_padding(rov_ims, 1), axes=[-1, -2, -3]) + ksp = cp.array(ksp) + csm_rovir = EspiritCalib(ksp, show_pbar=True, device=cp.cuda.Device(device=0)).run() + else: + for sl in range(rov_ims.shape[1]): + eprint(f"{sl} out of ") + [csm_, rho] = calculate_csm_walsh_gpu(rov_ims[:,sl,:,:].squeeze()) + csm_rovir[:,sl:sl+1,:,:] = csm_ + #storage.store({f"W":W, "csm_rovir": csm_rovir}, name=f"ROVIR_W_{size_path}") + storage.store({f"W":W, "csm_rovir": csm_rovir, "rovir_sig": rovir_sig_, "rovir_int": rovir_int_}, name=f"ROVIR_W_{size_path}") + else: + storage.store({f"W":W}, name=f"ROVIR_W_{size_path}") + + connection.send(data) + +if __name__ == '__main__': + gadgetron.external.listen(2020, ROVirGadget) \ No newline at end of file diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/SlicerGadget.py b/toolboxes/nhlbi_gt_toolbox/utils/python/SlicerGadget.py new file mode 100644 index 0000000..2819518 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/SlicerGadget.py @@ -0,0 +1,109 @@ +import gadgetron +import pyigtl +import numpy as np + +import pyigtl.comm +import socket +import struct +import fcntl +import sys +import signal +import time + + +def _parse_params(xml): + return {p.get('name'): p.get('value') for p in xml.iter('property')} + +def SlicerGadget(connection): + print("STARTING SLICER GADGET") + + params = _parse_params(connection.config) + + if "local_port" in params: + local_port = int(params["local_port"]) + else: + local_port = 9004 + + # try to load frame_discard from a file, if it doesn't exist default to 6 + try: + frame_discard = int(np.load("frame_discard.npy")) + except: + frame_discard = 6 + print(f"Frame discard set to: {frame_discard}") + + rep_count = frame_discard+1 + rep_buffer = [] + pyigtl.comm.OpenIGTLinkServer.__init__ = patched_init + igt_server = pyigtl.OpenIGTLinkServer(local_port, local_server=False, iface="0.0.0.0") + for data in connection: + st = time.time() + if data.repetition == rep_count: + rep_buffer.append(data.data) + else: + im = pyigtl.ImageMessage(np.flip((np.stack(rep_buffer,axis=0).squeeze()).transpose(2,1,0), axis=0)) + igt_server.send_message(im) + #print(f"Slicer Gadget Sent to 3D Slicer in: {time.time() - st}s") + rep_buffer = [data.data] + rep_count = data.repetition + #im = pyigtl.ImageMessage(np.fliplr(data.data.squeeze().transpose(1,2,0))) + #igt_server.send_message(im) + #print(f"Slicer Gadget Sent to 3D Slicer in: {time.time() - st}s") + connection.send(data) + igt_server.server_close() + + + +# Monkey patch the OpenIGTLinkServer class BEFORE using it +def patched_init(self, port=None, local_server=True, iface=None, start_now=True): + pyigtl.comm.OpenIGTLinkBase.__init__(self) + self.port = port + + if iface is None: + iface = 'eth0' + + if local_server: + self.host = "127.0.0.1" + else: + # Check if iface is already an IP address + try: + socket.inet_aton(iface) # This will succeed if it's a valid IP + self.host = iface # Use the IP directly + except socket.error: + # It's an interface name, so resolve it + if sys.platform.startswith('win32'): + self.host = socket.gethostbyname(socket.gethostname()) + elif sys.platform.startswith('linux'): + soc = socket.socket(socket.AF_INET, socket.SOCK_DGRAM) + try: + ifname = iface + # Fix: encode the string to bytes for Python 3 + self.host = socket.inet_ntoa(fcntl.ioctl( + soc.fileno(), 0x8915, + struct.pack('256s', ifname.encode('utf-8')[:15]) + )[20:24]) + except: + ifname = 'lo' + self.host = socket.inet_ntoa(fcntl.ioctl( + soc.fileno(), 0x8915, + struct.pack('256s', ifname.encode('utf-8')[:15]) + )[20:24]) + else: + self.host = iface + + # Continue with the rest of the original initialization + pyigtl.comm.SocketServer.TCPServer.allow_reuse_address = True + pyigtl.comm.SocketServer.TCPServer.__init__(self, (self.host, self.port), pyigtl.comm.TCPRequestHandler) + + # Signal handlers + self._previous_signal_handlers = {} + self._previous_signal_handlers[signal.SIGTERM] = signal.signal(signal.SIGTERM, self._signal_handler) + self._previous_signal_handlers[signal.SIGINT] = signal.signal(signal.SIGINT, self._signal_handler) + + if start_now: + self.start() + + +if __name__ == '__main__': + # start with a monkey patched igt_server __init__ + pyigtl.comm.OpenIGTLinkServer.__init__ = patched_init + gadgetron.external.listen(2020, SlicerGadget) \ No newline at end of file diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/mrdhelper.py b/toolboxes/nhlbi_gt_toolbox/utils/python/mrdhelper.py new file mode 100644 index 0000000..77f5638 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/mrdhelper.py @@ -0,0 +1,293 @@ +# MRD Helper functions + +# Taken from https://raw.githubusercontent.com/kspaceKelvin/python-ismrmrd-server/refs/heads/master/mrdhelper.py +# kelvin chow + +import ismrmrd +import re +import base64 + +def update_img_header_from_raw(imgHead, rawHead): + """Populate ImageHeader fields from AcquisitionHeader""" + + if rawHead is None: + return imgHead + + # # These fields are not translated from the raw header, but filled in + # # during image creation by from_array + # imgHead.data_type = + # imgHead.matrix_size = + # imgHead.channels = + + # # This is mandatory, but must be filled in from the XML header, + # # not from the acquisition header + # imgHead.field_of_view = + + imgHead.version = rawHead.version + imgHead.flags = rawHead.flags + imgHead.measurement_uid = rawHead.measurement_uid + + imgHead.position = rawHead.position + imgHead.read_dir = rawHead.read_dir + imgHead.phase_dir = rawHead.phase_dir + imgHead.slice_dir = rawHead.slice_dir + imgHead.patient_table_position = rawHead.patient_table_position + + imgHead.average = rawHead.idx.average + imgHead.slice = rawHead.idx.slice + imgHead.contrast = rawHead.idx.contrast + imgHead.phase = rawHead.idx.phase + imgHead.repetition = rawHead.idx.repetition + imgHead.set = rawHead.idx.set + + imgHead.acquisition_time_stamp = rawHead.acquisition_time_stamp + imgHead.physiology_time_stamp = rawHead.physiology_time_stamp + + # Defaults, to be updated by the user + imgHead.image_type = ismrmrd.IMTYPE_MAGNITUDE + imgHead.image_index = 1 + imgHead.image_series_index = 0 + + imgHead.user_float = rawHead.user_float + imgHead.user_int = rawHead.user_int + + return imgHead + +def get_userParameterLong_value(metadata, name): + """Get a value from MRD Header userParameterLong (returns None if key not found)""" + if metadata.userParameters is not None: + for param in metadata.userParameters.userParameterLong: + if param.name == name: + return int(param.value) + return None + +def get_userParameterDouble_value(metadata, name): + """Get a value from MRD Header userParameterDouble (returns None if key not found)""" + if metadata.userParameters is not None: + for param in metadata.userParameters.userParameterDouble: + if param.name == name: + return float(param.value) + return None + +def get_userParameterString_value(metadata, name): + """Get a value from MRD Header userParameterString (returns None if key not found)""" + if metadata.userParameters is not None: + for param in metadata.userParameters.userParameterString: + if param.name == name: + return param.value + return None + +def get_userParameterBase64_value(metadata, name): + """Get a value from MRD Header userParameterBase64 (returns None if key not found)""" + if metadata.userParameters is not None: + for param in metadata.userParameters.userParameterBase64: + if param.name == name: + return base64.b64decode(param.value).decode('utf-8') + return None + +def get_meta_value(meta, key): + """Get a value from MRD Meta Attributes (returns None if key not found)""" + if key in meta.keys(): + return meta[key] + else: + return None + +def extract_minihead_bool_param(miniHead, name): + """Extract a bool parameter from the serialized text of the ICE MiniHeader""" + val = extract_minihead_param(miniHead, name, 'ParamBool') + + if val is None: + return False + elif val.strip('" ').lower() == 'true'.lower(): + return True + else: + return False + +def extract_minihead_long_param(miniHead, name): + """Extract a long parameter from the serialized text of the ICE MiniHeader""" + val = extract_minihead_param(miniHead, name, 'ParamLong') + + if val is None: + return int(0) + else: + return int(val) + +def extract_minihead_double_param(miniHead, name): + """Extract a double parameter from the serialized text of the ICE MiniHeader""" + val = extract_minihead_param(miniHead, name, 'ParamDouble') + + if val is None: + return float(0) + else: + return float(val) + +def extract_minihead_string_param(miniHead, name): + """Extract a string parameter from the serialized text of the ICE MiniHeader""" + val = extract_minihead_param(miniHead, name, 'ParamString') + + return val.strip(' "') + +def extract_minihead_param(miniHead, name, strType): + """Extract a string parameter from the serialized text of the ICE MiniHeader""" + expr = r'(?<=<' + strType + r'."' + name + r'">)\s*[^}]*\s*' + res = re.search(expr, miniHead) + + if res is None: + return None + + # Strip off beginning '{' and whitespace, then split on newlines + values = res.group(0).strip('{\n ').split('\n') + + # Lines beginning with <> are properties -- ignore them + values = [val for val in values if bool(re.search(r'^\s*<\w+>', val)) is False] + + if len(values) != 1: + return None + else: + return values[0] + +def get_json_config_param(config, key, default=None, type='str'): + """ + Read a parameter from JSON config + Input: + - config : dict of parameters + - key : name (key) of parameter + - default : value if key is not present or config is invalid + - type : type casting of the parameter (int, float, string, bool) + Output: + - value of parameter, or default if absent + """ + if not isinstance(config, dict): + return default + + if not 'parameters' in config: + return default + + if not key in config['parameters']: + return default + + value = config['parameters'][key] + + if type == 'int': + return int(value) + elif (type == 'float') or (type == 'double'): + return float(value) + elif (type == 'string') or (type == 'str') or (type == 'choice'): + return str(value) + elif (type == 'bool') or (type == 'boolean'): + if isinstance(value, bool): + return value + elif 'true' in value.lower(): + return True + elif 'false' in value.lower(): + return False + else: + return default + else: + raise Exception("'type' must be int, float, string, or bool") + +def create_roi(x, y, rgb = (1, 0, 0), thickness = 1, style: int = 0, visibility: int = 1): + """ + Create an MRD-formatted ROI + Parameters: + - x (1D ndarray) : x coordinates in units of pixels, with (0,0) at the top left + - y (1D ndarray) : y coordinates in units of pixels, matching the length of x + - rgb (3 item tuple) : Colour as an (red, green, blue) tuple normalized to 1 + - thickness (float) : Line thickness + - style (int) : Line style (0 = solid, 1 = dashed) + - visibility (int) : Line visibility (0 = false, 1 = true) + Returns: + - roi (string list) : MRD-formatted ROI, intended to be stored as a MetaAttribute + with field name starting with "ROI_" + """ + xy = [(x[i], y[i]) for i in range(0, len(x))] # List of (x,y) tuples + + roi = [] + roi.append('%f' % rgb[0]) + roi.append('%f' % rgb[1]) + roi.append('%f' % rgb[2]) + roi.append('%f' % thickness) + roi.append('%d' % style) + roi.append('%d' % visibility) + + for i in range(0, len(xy)): + roi.append('%f' % xy[i][0]) + roi.append('%f' % xy[i][1]) + + return roi + +def parse_roi(roi): + """ + Parse an MRD-formatted ROI + Input: + - roi (string list) : MRD-formatted ROI from a MetaAttribute + Output: + - x (1D ndarray) : x coordinates in units of pixels, with (0,0) at the top left + - y (1D ndarray) : y coordinates in units of pixels, matching the length of x + - rgb (3 item tuple) : Colour as an (red, green, blue) tuple normalized to 1 + - thickness (float) : Line thickness + - style (int) : Line style (0 = solid, 1 = dashed) + - visibility (int) : Line visibility (0 = false, 1 = true) + """ + if (not isinstance(roi, list)) or (len(roi) < 8) or (len(roi) % 2): + raise Exception("ROI must be a list, have 6 metadata values, at least one coordinate, and an even number of values (x,y pairs)") + + fRoi = [float(x) for x in roi] + + rgb = tuple(fRoi[0:3]) + thickness = fRoi[3] + style = int(fRoi[4]) + visibility = int(fRoi[5]) + + x = fRoi[6::2] + y = fRoi[7::2] + + return x, y, rgb, thickness, style, visibility + +def create_text(x, y, rgb = (1, 0, 0), visibility: int = 1, string = ''): + """ + Create an MRD-formatted text object + Parameters: + - x (float) : x coordinate in units of pixels, with (0,0) at the top left + - y (float) : y coordinate in units of pixels + - rgb (3 item tuple) : Colour as an (red, green, blue) tuple normalized to 1 + - visibility (int) : Line visibility (0 = false, 1 = true) + - string (string) : Text string + Returns: + - txt (string list) : MRD-formatted text, intended to be stored as a MetaAttribute + with field name starting with "Text_" + """ + txt = [] + txt.append('%f' % rgb[0]) + txt.append('%f' % rgb[1]) + txt.append('%f' % rgb[2]) + txt.append('%f' % x) + txt.append('%f' % y) + txt.append('%d' % visibility) + txt.append('%s' % string) + + return txt + +def parse_text(txt): + """ + Parse an MRD-formatted text object + Input: + - txt (string list) : MRD-formatted text from a MetaAttribute + Output: + - x (float) : x coordinate in units of pixels, with (0,0) at the top left + - y (float) : y coordinate in units of pixels + - rgb (3 item tuple) : Colour as an (red, green, blue) tuple normalized to 1 + - visibility (int) : Line visibility (0 = false, 1 = true) + - string (string) : Text string + """ + if (not isinstance(txt, list)) or (len(txt) != 7): + raise Exception("txt must be a list that has exactly 7 metadata values") + + rgb = tuple([float(x) for x in txt[0:3]]) + x = float(txt[3]) + y = float(txt[4]) + visibility = int(float(txt[5])) + + string = txt[6] + + return x, y, rgb, visibility, string diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/registration/registration_gadget_call.py b/toolboxes/nhlbi_gt_toolbox/utils/python/registration/registration_gadget_call.py index a219ed8..8a932d1 100644 --- a/toolboxes/nhlbi_gt_toolbox/utils/python/registration/registration_gadget_call.py +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/registration/registration_gadget_call.py @@ -45,8 +45,8 @@ def eprint(*args, **kwargs): def registration_one_image(mov_image_np, ref_image_np,gpu_list=[]): GPU_freeM,dev_num=get_GPU_most_free(gpu_list) - print(f"Free Memory GPU {GPU_freeM} GPU num {dev_num}") - print("Registration of image 1: ", mov_image_np.shape) + eprint(f"Free Memory GPU {GPU_freeM} GPU num {dev_num}") + eprint("Registration of image 1: ", mov_image_np.shape) deformation_fields = reg.register_one_image_only_deformation(mov_image_np,ref_image_np,gpu_id=dev_num).transpose(1,2,3,0) return deformation_fields.astype(np.float32) @@ -69,7 +69,7 @@ def registration_images(images,ref_index=0,gpu_list=[]): eprint(f"----------------Running registration : Ref index {ref_index} Nbins {images.shape[0]}--------------") #with redirect_stdout(fnull) and redirect_stderr(fnull): deformation_fields = reg.register_images_only_deformation(images,ref_index,gpu_id=dev_num) - print(deformation_fields.shape) + eprint(deformation_fields.shape) deformation_fields = deformation_fields.transpose(2,3,4,1,0) np.nan_to_num(deformation_fields) eprint("Registration Time: ", time.time()-st) @@ -78,8 +78,8 @@ def registration_images(images,ref_index=0,gpu_list=[]): def registration_images_old(images,bidirectional=False,ref_index=0): - - print("Registration of images: ", images.shape) + + eprint("Registration of images: ", images.shape) #images t,nx,ny,nz #[0,2,3,1] #images = images.transpose(3,0,1,2) @@ -97,7 +97,7 @@ def registration_images_old(images,bidirectional=False,ref_index=0): eprint("----------------Running registration--------------") with redirect_stdout(fnull) and redirect_stderr(fnull): deformation_fields = reg.register_images_only_deformation(images,ref_index) - print(deformation_fields.shape) + eprint(deformation_fields.shape) deformation_fields = deformation_fields.transpose(2,3,4,1,0) np.nan_to_num(deformation_fields) eprint("Registration Time: ", time.time()-st) @@ -105,8 +105,8 @@ def registration_images_old(images,bidirectional=False,ref_index=0): return deformation_fields.astype(np.float32) def registration_images_back(images,bidirectional=False,ref_index=0): - - print("Registration of images: ", images.shape) + + eprint("Registration of images: ", images.shape) #images t,nx,ny,nz #[0,2,3,1] #images = images.transpose(3,0,1,2) @@ -127,7 +127,7 @@ def registration_images_back(images,bidirectional=False,ref_index=0): for idx in range(images.shape[0]): idxs = [ref_index,idx] deformation_fields[idx,...] = reg.register_images_only_deformation(images[idxs,...],1)[[0],...] - print(deformation_fields.shape) + eprint(deformation_fields.shape) deformation_fields = deformation_fields.transpose(2,3,4,1,0) np.nan_to_num(deformation_fields) eprint("Registration Time: ", time.time()-st) diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/registration/registration_oflow3D.py b/toolboxes/nhlbi_gt_toolbox/utils/python/registration/registration_oflow3D.py index 6619d8d..6a0e8ea 100644 --- a/toolboxes/nhlbi_gt_toolbox/utils/python/registration/registration_oflow3D.py +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/registration/registration_oflow3D.py @@ -1,6 +1,9 @@ import warnings warnings.filterwarnings("ignore") #import cv2 +import contextlib +import io +import os import opticalflow3D import cupy as cp from cupyx.scipy.ndimage import map_coordinates, zoom @@ -10,6 +13,13 @@ import gc import torch + +def _quiet_calculate_flow(farneback, *args, **kwargs): + """Call farneback.calculate_flow with stdout suppressed.""" + with open(os.devnull, 'w') as devnull: + with contextlib.redirect_stdout(devnull): + return farneback.calculate_flow(*args, **kwargs) + # from numba import config # config.CUDA_ENABLE_MINOR_VERSION_COMPATIBILITY = True def get_GPU_most_free(): @@ -46,9 +56,10 @@ def register_images_only_deformation(input_images, ref_index, filter_size=9, gpu filter_size=filter_size, presmoothing=0, # Default, none filter_type="gaussian", - sigma_k=0.05) + sigma_k=0.05, + device_id=gpu_id) - output_vx, output_vy, output_vz, x = farneback.calculate_flow( + output_vx, output_vy, output_vz, x = _quiet_calculate_flow(farneback, 0.05 * cp.abs(ref_image / cp.max(ref_image.ravel())), 0.05 * cp.abs(mov_image / cp.max(mov_image.ravel())), total_vol=(ref_image.shape[0], ref_image.shape[1], ref_image.shape[2]), @@ -84,9 +95,10 @@ def register_one_image_only_deformation(mov_image_np, ref_image_np, filter_size= filter_size=filter_size, presmoothing=0, # Default, none filter_type="gaussian", - sigma_k=0.05) - - output_vx, output_vy, output_vz, x = farneback.calculate_flow( + sigma_k=0.05, + device_id=gpu_id) + + output_vx, output_vy, output_vz, x = _quiet_calculate_flow(farneback, 0.05 * cp.abs(ref_image / cp.max(ref_image.ravel())), 0.05 * cp.abs(mov_image / cp.max(mov_image.ravel())), total_vol=(ref_image.shape[0], ref_image.shape[1], ref_image.shape[2]), @@ -125,12 +137,13 @@ def register_images(input_images, ref_index, filter_size=9, gpu_id=0): filter_size=filter_size, presmoothing=2, # Default, none filter_type="gaussian", - sigma_k=0.05) + sigma_k=0.05, + device_id=gpu_id) for ind in range(0, nimages): mov_image = images[ind, ...].squeeze() - output_vx, output_vy, output_vz, x = farneback.calculate_flow( + output_vx, output_vy, output_vz, x = _quiet_calculate_flow(farneback, 0.05 * cp.abs(ref_image / cp.max(ref_image.ravel())), 0.05 * cp.abs(mov_image / cp.max(mov_image.ravel())), total_vol=(ref_image.shape[0], ref_image.shape[1], ref_image.shape[2]), @@ -294,7 +307,6 @@ def findGPUs(): f,t = torch.cuda.mem_get_info() #memcap.append(float(torch.cuda.get_device_properties(devno).total_memory)/float(1024**3)) memcap.append(float(f)/float(1024**3)) - print(f'Memory: {memcap}') return np.argsort(np.array(memcap)) diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/storage_server.py b/toolboxes/nhlbi_gt_toolbox/utils/python/storage_server.py new file mode 100644 index 0000000..11e7df4 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/storage_server.py @@ -0,0 +1,259 @@ +# copied from https://github.com/fyrdahl/mrd-storage-client/tree/main +# modified slightly to utilize the GT storage server calls. + +import pickle +from typing import Any, Dict, Iterator, List, Optional, Union + +import requests +from requests.adapters import HTTPAdapter +from requests_toolbelt.sessions import BaseUrlSession +from urllib3.util.retry import Retry +import numpy as np + +DEFAULT_TIMEOUT = 3 + +# API documentation: https://github.com/ismrmrd/mrd-storage-server/blob/main/README.md + + +class StorageException(Exception): + """Base exception for Storage""" + + +class ConnectionStorageException(StorageException): + """Raised when a connection error occured""" + + +class SerializeStorageException(StorageException): + """Raised when a (de-)serialization error occured""" + + +class HealthcheckStorageException(StorageException): + """Raised when an healthcheck failed""" + + +class Blob: + def __init__(self, **kwargs) -> None: + """Initializes a Blob instance, dynamically assigning the given keyword arguments + as attributes to the instance. + """ + self.__dict__.update(**kwargs) + + def get_data(self) -> bytes: + """Fetch data from blob""" + response = requests.get(self.data) + return response.content + + def get(self, key: str, default=None): + """Provides a dict-like get method for the Blob instance.""" + return getattr(self, key, default) + + +class Storage: + def __init__( + self, + address: str, + port: int, + subject: str = "$null", + device: Optional[str] = None, + session: Optional[str] = None, + ): + self.address = address + self.port = port + + base_url = f"http://{address}:{port}/" + self.http = init_http(base_url) + + # Mandatory tags + self.base_params = {"subject": subject} + + # Optional tags + if device: + self.base_params["device"] = device + + if session: + self.base_params["session"] = session + + self.subject = subject + self.device = device + self.session = session + + def healthcheck(self) -> None: + """Healthcheck can be used to verify that that the server is functioning. + Raises HealthcheckStorageException if healthcheck failed or + ConnectionStorageException if the connection failed. + + Example: + >>> storage = Storage("localhost", 3333) + >>> storage.healthcheck() + """ + try: + response = self.http.get("healthcheck") + if not response.ok: + raise HealthcheckStorageException + except requests.exceptions.ConnectionError: + raise ConnectionStorageException + + def store(self, obj, name=None, ttl=None, custom_tags=None): + """Serialize and store an object as a blob. + Raises SerializeStorageException if the serialization failed. + + Example: + >>> import numpy as np + >>> storage = Storage("localhost", 3333) + >>> obj = {"key": "value"} + >>> storage.store(obj) + + """ + if custom_tags is None: + custom_tags = {} + try: + payload = pickle.dumps(obj) + except pickle.PicklingError as e: + raise SerializeStorageException(e) from e + + new_params = dict(self.base_params) + + if name: + new_params["name"] = name + if ttl: + new_params["_ttl"] = ttl + + new_params = {**new_params, **custom_tags} + self.http.post("v1/blobs/data", payload, params=new_params) + + def fetch(self, name=None, at=None, custom_tags=None): + """Return a list containing the data from all matching blobs + + Example: + >>> storage = Storage("localhost", 3333) + >>> data = storage.fetch() + + """ + if custom_tags is None: + custom_tags = {} + + blobs = self._search(name, at, custom_tags) + return [self._load_object(blob.get_data()) for blob in blobs] + + def fetch_blobs(self, name=None, at=None, custom_tags=None): + """Yield iterator for all matching blob objects + + Example: + >>> storage = Storage("localhost", 3333) + >>> for blob in storage.fetch_blobs(): + >>> data = blob.get_data() + + """ + if custom_tags is None: + custom_tags = {} + + yield from self._search(name, at, custom_tags) + + def fetch_latest(self, name=None, at=None, custom_tags=None): + """Shortcut to get the data from the latest blob matching a search query + + Example: + >>> storage = Storage("localhost", 3333) + >>> data = storage.fetch_latest() + + """ + if custom_tags is None: + custom_tags = {} + + new_params = dict(self.base_params) + + if name: + new_params["name"] = name + if at: + new_params["_at"] = at + + new_params = {**new_params, **custom_tags} + response = self.http.get("v1/blobs/data/latest", params=new_params) + return self._load_object(response.content) + + def close(self): + """Close the underlying HTTP session""" + self.http.close() + + def _search(self, name=None, at=None, custom_tags=None): + if custom_tags is None: + custom_tags = {} + + new_params = dict(self.base_params) + + if name: + new_params["name"] = name + if at: + new_params["_at"] = at + + new_params = {**new_params, **custom_tags} + + response = self.http.get("v1/blobs", params=new_params) + return self._create_blob_obj(response.json()) + + def _load_object(self, data): + try: + return pickle.loads(data) + except pickle.UnpicklingError as e: + try: + # get the dimensionality of the array. + ints = np.frombuffer(data, np.int64) + num_dims = ints[0] + dims = ints[1:num_dims+1] + + # PK HACK: for now, we assume float32... but this could be incorrect. + n_discard = (num_dims+1) # discard 2x amount because in64 and float32 mismatch. + np_data = np.frombuffer(data, dtype=' + # xk = del_n + # yk = new variable + # tk = new variable + costs = [] + if del_0 is None: + del_n = cp.zeros(xn_1.shape) + else: + del_n = cp.copy(del_0) + + w_n_1 = del_n + theta_n_1 = 1 + + z_n_1 = del_n + del_n_1 = del_n + + gamma_n_1 = 1 + + for i in range(n_iter): + if (i+1) < n_iter: + theta_n = 0.5 * (1 + cp.sqrt((4 * cp.square(theta_n_1)) + 1)) + else: + theta_n = 0.5 * (1 + cp.sqrt((8 * cp.square(theta_n_1)) + 1)) + + gamma_n = step * (((2 * theta_n_1) + theta_n - 1) / (theta_n)) + w_n = del_n_1 - (step*A.N*(del_n_1+xn_1)) + (step*ATyn) + + z_n = w_n + (((theta_n_1 - 1) / (theta_n)) * (w_n - w_n_1)) + (((theta_n_1) / (theta_n)) * (w_n - del_n_1)) + \ + (((theta_n_1 - 1) / ((1/step)*gamma_n_1*theta_n)) * (z_n_1 - del_n_1)) + del_n = soft_threshold(z_n, lam) + del_n = crop_FOV_edge(del_n, 0.8) + + # update theta, del, gamma, w, and z. + theta_n_1 = theta_n + del_n_1 = del_n + gamma_n_1 = gamma_n + w_n_1 = w_n + z_n_1 = z_n + + if cost_fn is not None: + costs.append(cost_fn(del_n+xn_1) + (lam * cp.abs(cp.sum(del_n)))) + + if cost_fn is not None: + return [del_n, costs] + return del_n + +def create_circular_mask(h, w, center=None, radius=None): + + if center is None: # use the middle of the image + center = (int(w/2), int(h/2)) + if radius is None: # use the smallest distance between the center and image walls + radius = min(center[0], center[1], w-center[0], h-center[1]) + + Y, X = cp.ogrid[:h, :w] + dist_from_center = cp.sqrt((X - center[0])**2 + (Y-center[1])**2) + + mask = dist_from_center <= radius + return mask + + +# crop edges of the FOV +def crop_FOV_edge(x, amt=1): + if amt > 1: + eprint("AMT NEEDS TO BE LESS THAN 1") + amt = 1 + x_fft = cp.fft.fftshift(cp.fft.fft2(x, axes=[0,1]), axes=[0,1]) + mask = create_circular_mask(x.shape[0], x.shape[1], radius=(x.shape[0]*amt/2)) + mask = cp.expand_dims(mask, 2) + x_fft = cp.multiply(x_fft , mask) + x_recomp = cp.fft.ifft2(cp.fft.ifftshift(x_fft, axes=[0,1]), axes=[0,1]) + return x_recomp + +def norm(x): + # L2 norm of vectorized x. + # flatten x + norm = (cp.sum(cp.square(cp.abs(x.flatten())))) + return norm + + +def cg(Af, b, x0, niter=20, tol=1e-6): + x = cp.copy(x0) + r = Af(x) - b + p = cp.copy(r) + for i in range(niter): + eprint(f"iter {i}") + Ap = Af(p) + rsold = cp.sum(cp.square(cp.abs(r))) + alpha = rsold / np.vdot(p, Ap) + x = x + alpha * p + r = r - alpha * Ap + rsnew = cp.sum(cp.square(cp.abs(r))) + if rsnew < tol: + break + p = r + (rsnew / rsold) * p + rsold = rsnew + return x + +def FISTA_iteration(iter, xk_1, gradf, step_size, fista_dict=None, threshold=0, max_iter=10): + if iter == 0: + fista_dict = { + 'yk_1': xk_1, + 'tk_1': 1 + } + elif iter >= max_iter: + eprint("not iterating, max_iter reached") + return [xk_1, fista_dict] + tk_1 = fista_dict['tk_1'] + yk_1 = fista_dict['yk_1'] + + tk = 0.5 * (1 + cp.sqrt((4 * cp.square(tk_1)) + 1)) + xk = yk_1 - (step_size * gradf(yk_1)) + xk = soft_threshold(xk, threshold) + yk = xk + (((tk_1 - 1) / tk) * (xk - xk_1)) + return [xk, {'yk_1': yk, 'tk_1': tk}] + +def POGM_iteration(iter, xk_1, gradf, step_size, pogm_dict=None, threshold=0, max_iter=10): + if iter == 0: + pogm_dict = { + 'zk_1': xk_1, + 'thetak_1': 1, + 'wk_1': xk_1, + 'gammak_1': 1 + } + elif iter >= max_iter: + eprint("not iterating, max_iter reached") + return [xk_1, pogm_dict] + thetak_1 = pogm_dict['thetak_1'] + gammak_1 = pogm_dict['gammak_1'] + zk_1 = pogm_dict['zk_1'] + wk_1 = pogm_dict['wk_1'] + + if iter < max_iter: + thetak = 0.5 * (1 + cp.sqrt((4 * cp.square(thetak_1)) + 1)) + else: + thetak = 0.5 * (1 + cp.sqrt((8 * cp.square(thetak_1)) + 1)) + + gammak = step_size * (((2 * thetak_1) + thetak - 1) / (thetak)) + wk = xk_1 - (step_size * gradf(xk_1)) + + zk = wk + (((thetak_1 - 1) / thetak) * (wk - wk_1)) + (((thetak_1) / thetak) * (wk - xk_1)) + \ + + (((thetak_1 - 1) * step_size / (gammak_1 * thetak)) * (zk_1 - xk_1)) + + xk = soft_threshold(zk, threshold) + return [xk, {'zk_1': zk, 'thetak_1': thetak, 'wk_1': wk, 'gammak_1': gammak}] + +def online_TCR_POGM_2(A, xn_1, ATyn, lam, step, n_iter=40, del_0=None, cost_fn=None): + + # construct gradf function + gradf = lambda x: A.N * (x + xn_1) - ATyn + + if del_0 is None: + del_n = cp.zeros(xn_1.shape) + else: + del_n = cp.copy(del_0) + + #initialize variables + pogm_dict = None + costs = [] + + for i in range(n_iter): + [del_n, pogm_dict] = POGM_iteration(i, del_n, gradf, step, pogm_dict, lam*step, n_iter) + del_n = crop_FOV_edge(del_n, 0.6) + if cost_fn is not None: + costs.append(cost_fn(del_n+xn_1)) + + if cost_fn is not None: + return [del_n, costs] + + return del_n + +# set up an online STCR reconstruction +def online_STCR_ISTA_timed(E, G, xn_1, Ahyn, lambdat, lambdas, step_size, mu=1 , yn=None, deln=None, time_recon=200, cost_fn=None): + + start_time = time.time() + running = True + + costs = [] + if deln is None: + deln = cp.zeros(xn_1.shape, Ahyn.dtype) + # initialize etan + zn = G*(xn_1 + deln) + etan = cp.zeros(zn.shape, Ahyn.dtype) + + niter = 0 + + while running: + zn = soft_threshold(G*(deln+xn_1) + etan, lambdas/mu) + inner_iter = 1 + + for i_i in range(inner_iter): + deln = soft_threshold(deln - step_size *(E.H * E * (deln+xn_1) - Ahyn + (G.H*(G*(deln+xn_1) - zn + etan)*mu)), lambdat*step_size) + #cost = norm(E * deln - yn + E*xn_1) + norm((cp.sqrt(mu) * (G * deln)) - (cp.sqrt(mu) * (zn-etan-(G*xn_1)))) + (lambdat * (cp.sum(cp.abs(deln)))) + etan = etan - zn + (G*(deln+xn_1)) + if cost_fn is not None: + costs.append(cost_fn(deln)) + + cp.cuda.stream.get_current_stream().synchronize() + elapsed = time.time() - start_time + + niter = niter+1 + if elapsed*1000 > time_recon: + running = False + + + eprint(f"niter {niter}") + if cost_fn is not None: + return [deln, costs] + return deln + + +# set up an online STCR reconstruction +def online_STCR_ISTA(E, G, xn_1, Ahyn, lambdat, lambdas, step_size, mu=1 , yn=None, deln=None, niter=20, cost_fn=None): + costs = [] + if deln is None: + deln = cp.zeros(xn_1.shape, Ahyn.dtype) + # initialize etan + zn = G*(xn_1 + deln) + etan = cp.zeros(zn.shape, Ahyn.dtype) + for n_i in range(niter): + zn = soft_threshold(G*(deln+xn_1) + etan, lambdas/mu) + inner_iter = 1 + """ + # cg instead? + Af = lambda x: (E.H * E * x + (mu * (G.H * G * x))) + y = (E.H * (yn + E*xn_1)) + (mu * (G.H * (zn - etan - (G*xn_1)))) + deln = cg(Af, y, deln, niter=inner_iter) + """ + for i_i in range(inner_iter): + deln = soft_threshold(deln - step_size *(E.H * E * (deln+xn_1) - Ahyn + (G.H*(G*(deln+xn_1) - zn + etan)*mu)), lambdat*step_size) + #cost = norm(E * deln - yn + E*xn_1) + norm((cp.sqrt(mu) * (G * deln)) - (cp.sqrt(mu) * (zn-etan-(G*xn_1)))) + (lambdat * (cp.sum(cp.abs(deln)))) + etan = etan - zn + (G*(deln+xn_1)) + if cost_fn is not None: + costs.append(cost_fn(deln)) + if cost_fn is not None: + return [deln, costs] + return deln + + +# set up an online STCR reconstruction +def online_STCR_ISTA_2(E, G, xn_1, Ahyn, lambdat, lambdas, step_size, mu=1 , yn=None, deln=None, niter=20, cost_fn=None): + + costs = [] + if deln is None: + deln = cp.zeros(xn_1.shape, Ahyn.dtype) + # initialize etan + zn = G*(xn_1 + deln) + etan = cp.zeros(zn.shape, Ahyn.dtype) + pogm_dict_outer = None + for n_i in range(niter): + gradout = lambda x: (zn - etan - (G*(deln+xn_1))) + step_size_outer = 1 + + inner_iter = 2 + pogm_dict_inner = None + for i_i in range(inner_iter): + gradf = lambda x:(E.H * E * (deln+xn_1) - Ahyn + (G.H*(G*(deln+xn_1) - zn + etan)*mu) ) + [deln, pogm_dict_inner] = FISTA_iteration(i_i, deln, gradf, step_size, pogm_dict_inner, lambdat*step_size, inner_iter) + + #zn = soft_threshold(G*(deln+xn_1) + etan, lambdas/mu) + [zn, pogm_dict_outer] = FISTA_iteration(n_i, zn, gradout, step_size_outer, pogm_dict_outer, lambdas/mu, niter) + + etan = etan - zn + (G*(deln+xn_1)) + if cost_fn is not None: + costs.append(cost_fn(deln)) + if cost_fn is not None: + return [deln, costs] + return deln + +# set up an online STCR reconstruction +def online_STCR_ISTA_2_timed(E, G, xn_1, Ahyn, lambdat, lambdas, step_size, mu=1 , yn=None, deln=None, time_recon=200, cost_fn=None): + + start_time = time.time() + running = True + n_i = 0 + + costs = [] + if deln is None: + deln = cp.zeros(xn_1.shape, Ahyn.dtype) + # initialize etan + zn = G*(xn_1 + deln) + etan = cp.zeros(zn.shape, Ahyn.dtype) + pogm_dict_outer = None + while running: + gradout = lambda x: (zn - etan - (G*(deln+xn_1))) + step_size_outer = 1 + + inner_iter = 1 + pogm_dict_inner = None + for i_i in range(inner_iter): + gradf = lambda x:(E.H * E * (deln+xn_1) - Ahyn + (G.H*(G*(deln+xn_1) - zn + etan)*mu) ) + [deln, pogm_dict_inner] = FISTA_iteration(i_i, deln, gradf, step_size, pogm_dict_inner, lambdat*step_size, inner_iter) + + #zn = soft_threshold(G*(deln+xn_1) + etan, lambdas/mu) + [zn, pogm_dict_outer] = FISTA_iteration(n_i, zn, gradout, step_size_outer, pogm_dict_outer, lambdas/mu, 100) + + etan = etan - zn + (G*(deln+xn_1)) + if cost_fn is not None: + costs.append(cost_fn(deln)) + + cp.cuda.stream.get_current_stream().synchronize() + elapsed = time.time() - start_time + eprint(elapsed) + + n_i = n_i + 1 + if elapsed*1000 > time_recon: + running = False + eprint(n_i) + if cost_fn is not None: + return [deln, costs] + return deln + + + +def crop_half_FOV(image, dims=(0,1), size=None): + h,w = (image.shape[dims[0]], image.shape[dims[1]]) + if size is None: + crop_h = h // 2 + crop_w = w // 2 + size = [crop_h, crop_w] + + start_indices = [ + (h - size[0]) // 2, + (w - size[1]) // 2 + ] + + eprint(start_indices) + eprint(size) + + # Create a slice object for each dimension + slices = [slice(None) for i in range(len(image.shape))] + slices[dims[0]] = slice(start_indices[0], start_indices[0] + size[0]) + slices[dims[1]] = slice(start_indices[1], start_indices[1] + size[1]) + + # slices = [slice(start_indices[0], start_indices[0] + size[0]), + # slice(start_indices[1], start_indices[1] + size[1])] + + # Use numpy's advanced indexing to crop the image + return image[tuple(slices)] + + +def soft_threshold(x, threshold): + x_phase = cp.angle(x) + # do I need cp.sign(x)? Probably only if it is a real numpber. maybe, 'x_phase' covers both. + return cp.maximum(cp.abs(x) - threshold, 0) * cp.exp(1j*x_phase) + + +def remove_zero_padding(arr, axis): + """ + Remove zero padding along a specified axis in a NumPy array. + + Parameters: + arr (np.ndarray): Input array with potential zero padding. + axis (int): Axis along which to remove zero padding. + + Returns: + np.ndarray: Array with zero padding removed along the specified axis. + """ + if axis < 0 or axis >= arr.ndim: + raise ValueError(f"Invalid axis {axis} for array with {arr.ndim} dimensions.") + + # Create a boolean mask where True represents non-zero slices + non_zero_mask = np.any(arr != 0, axis=tuple(i for i in range(arr.ndim) if i != axis)) + + # Use the boolean mask to slice the array + slices = [slice(None)] * arr.ndim + slices[axis] = slice(None) # Keep the full range for the axis we are slicing + arr_sliced = arr[tuple(slices)] + + # Now use the boolean mask to filter out zero-only slices + arr_trimmed = arr_sliced.take(np.where(non_zero_mask)[0], axis=axis) + + return arr_trimmed + +def analyticaldcf(trajectory, adc_dwell=1e-6, ns=1): + kxx = np.array(trajectory[:,0]) + kyy = np.array(trajectory[:,1]) + kzz = np.array(trajectory[:,2]) + + # kx = kx[ndiscard:,0] + # ky = ky[ndiscard:,0] + gx = np.diff(np.concatenate(([0], kxx)), axis=0)/adc_dwell/42.58e6 + gy = np.diff(np.concatenate(([0], kyy)), axis=0)/adc_dwell/42.58e6 + + # Analytical DCF formula + # 1. Hoge RD, Kwan RKS, Bruce Pike G. Density compensation functions for spiral MRI. + # Magnetic Resonance in Medicine. 1997;38(1):117-128. doi:10.1002/mrm.1910380117 + cosgk = np.cos(np.arctan2(kxx, kyy) - np.arctan2(gx, gy)) + w = np.sqrt(kxx*kxx+kyy*kyy)*np.sqrt(gx*gx+gy*gy)*np.abs(cosgk) + w[-int(ns//2):] = w[-int(ns//2)] # need this to correct weird jump at the end and improve SNR + w = w/np.max(w) + return w + +def resize_with_crop_or_pad(image, target_shape): + """ + Resize a 2D numpy array by center cropping and/or zero padding. + + Args: + image: Input numpy array of shape (nx, ny) + target_shape: Tuple (target_x, target_y) for desired output shape + + Returns: + Resized numpy array of shape target_shape + """ + current_shape = image.shape + target_x, target_y = target_shape + current_x, current_y = current_shape + + # Handle x dimension + if current_x > target_x: + # Need to crop in x dimension + start_x = (current_x - target_x) // 2 + end_x = start_x + target_x + image = image[start_x:end_x, :] + elif current_x < target_x: + # Need to pad in x dimension + pad_x = target_x - current_x + pad_before_x = pad_x // 2 + pad_after_x = pad_x - pad_before_x + image = np.pad(image, ((pad_before_x, pad_after_x), (0, 0)), mode='constant', constant_values=0) + + # Handle y dimension + current_y = image.shape[1] # Update current_y after potential x operations + if current_y > target_y: + # Need to crop in y dimension + start_y = (current_y - target_y) // 2 + end_y = start_y + target_y + image = image[:, start_y:end_y] + elif current_y < target_y: + # Need to pad in y dimension + pad_y = target_y - current_y + pad_before_y = pad_y // 2 + pad_after_y = pad_y - pad_before_y + image = np.pad(image, ((0, 0), (pad_before_y, pad_after_y)), mode='constant', constant_values=0) + + return image + + +def update_nonant(image, str_update, crop_percent=0.75): + # first, extract the image from xn based on the str_update + if str_update[0] == 'c': + # coronal + image_nonant = image[:, int(str_update[1:]), :] + if str_update[0] == 's': + # sagittal, so extract the image directly + image_nonant = image[:,:,int(str_update[1:])] + if str_update[0] == 'a': + # axial + image_nonant = image[int(str_update[1:]), :, :] + + # since each view is of different sizes (xy, xz, yz) we want to make sure the final image is zero-padded or cropped appropriately + image_size = (int(image.shape[0] * crop_percent), int(image.shape[1] * crop_percent)) + image_nonant = resize_with_crop_or_pad(image_nonant, image_size) + return image_nonant + + +def compute_preemphasis_order(n_slices): + """ + Compute the 1-based ordering needed to pre-compensate for a 90-degree CLOCKWISE rotation + after column-major montage arrangement. + + Args: + n_slices: Total number of image slices (must be a perfect square) + + Returns: + List of 1-based indices to reorder input array before montage + """ + # Verify it's a perfect square + grid_size = int(np.ceil(np.sqrt(n_slices))) + + desired_final = np.arange(1, n_slices + 1).reshape(grid_size, grid_size) + + before_rotation = np.rot90(desired_final, k=1) # k=1 is counterclockwise + + preemphasis_values = before_rotation.T.flatten() - 1 + + # take any preemphasis_values that are > n_slices and remove them (prevent wrong slicing) + preemphasis_values = preemphasis_values[preemphasis_values < n_slices] + + return preemphasis_values.tolist() + + +def draw_number_indicators(image_array: np.ndarray, + indicator_size: int = 20, + position: str = 'top-left', + border_width: int = 2, + fill_color: Union[int, float] = 255, + border_color: Union[int, float] = 0) -> np.ndarray: + """ + Draw geometric indicators (squares) on images to represent their index numbers. + Each image gets a number of squares equal to its 1-based index. + + Args: + image_array: Input array of shape (height, width, num_images) + indicator_size: Size of each square indicator + position: Where to place indicators ('top-left', 'top-right', 'bottom-left', 'bottom-right') + border_width: Width of the border around each square + fill_color: Color/intensity for the filled part of squares + border_color: Color/intensity for the border of squares + + Returns: + Modified image array with number indicators drawn on each image + """ + + # Create a copy to avoid modifying the original array + result_array = image_array.copy() + + # Get dimensions + height, width, num_images = image_array.shape + + for i in range(num_images): + # Number of squares to draw (1-indexed) + num_squares = i + 1 + + # Calculate starting position based on position parameter + if position == 'top-left': + start_x, start_y = 5, 5 + dx, dy = indicator_size + 5, 0 # Horizontal arrangement + elif position == 'top-right': + start_x = width - 5 - num_squares * (indicator_size + 5) + 5 + start_y = 5 + dx, dy = indicator_size + 5, 0 + elif position == 'bottom-left': + start_x, start_y = 5, height - indicator_size - 5 + dx, dy = indicator_size + 5, 0 + elif position == 'bottom-right': + start_x = width - 5 - num_squares * (indicator_size + 5) + 5 + start_y = height - indicator_size - 5 + dx, dy = indicator_size + 5, 0 + else: + # Default to top-left + start_x, start_y = 5, 5 + dx, dy = indicator_size + 5, 0 + + # Draw the squares + for square_idx in range(num_squares): + x = int(start_x + square_idx * dx) + y = int(start_y + square_idx * dy) + + # Make sure we don't go out of bounds + if x + indicator_size >= width or y + indicator_size >= height: + break + + # Draw filled square + result_array[y:y+indicator_size, x:x+indicator_size, i] = fill_color + + # Draw border + if border_width > 0: + # Top and bottom borders + result_array[y:y+border_width, x:x+indicator_size, i] = border_color + result_array[y+indicator_size-border_width:y+indicator_size, x:x+indicator_size, i] = border_color + + # Left and right borders + result_array[y:y+indicator_size, x:x+border_width, i] = border_color + result_array[y:y+indicator_size, x+indicator_size-border_width:x+indicator_size, i] = border_color + + return result_array + +def replace_zero(array): + + for i in range(len(array)) : + if array[i] == 0 : + array[i] = 1 + return array + +def ram_schmidt(A,norm=True,row_vect=False): + """Orthonormalizes vectors by gram-schmidt process + + Parameters + ----------- + A : ndarray, + Matrix having vectors in its columns + + norm : bool, + Do you need Normalized vectors? + + row_vect: bool, + Does Matrix A has vectors in its rows? + + Returns + ------- + G : ndarray, + Matrix of orthogonal vectors + + """ + if row_vect : + # if true, transpose it to make column vector matrix + A = A.T + + no_of_vectors = A.shape[1] + G = A[:,0:1].copy() # copy the first vector in matrix + # 0:1 is done to to be consistent with dimensions - [[1,2,3]] + + # iterate from 2nd vector to number of vectors + for i in range(1,no_of_vectors): + + # calculates weights(coefficents) for every vector in G + numerator = A[:,i].dot(G) + denominator = np.diag(np.dot(G.T,G)) #to get elements in diagonal + weights = np.squeeze(numerator/denominator) + + # projected vector onto subspace G + projected_vector = np.sum(weights * G, + axis=1, + keepdims=True) + + # orthogonal vector to subspace G + orthogonalized_vector = A[:,i:i+1] - projected_vector + + # now add the orthogonal vector to our set + G = np.hstack((G,orthogonalized_vector)) + + if norm : + # to get orthoNORMAL vectors (unit orthogonal vectors) + # replace zero to 1 to deal with division by 0 if matrix has 0 vector + G = G/replace_zero(np.linalg.norm(G,axis=0)) + + if row_vect: + return G.T + + return G + +def fair(x, delta): + return delta**2 * (cp.abs(x / delta) - cp.log(1 + cp.abs(x / delta))) + +def d_fair(x, delta): + return x / (1 + cp.abs(x / delta)) + +def grad_cost(x, x_prev, A, y, D, delta, lambdat, lambdas): + AhAx = A.H * (A * x - y) + spatial = lambdas * (D.H * d_fair(D * x, delta)) + temporal = lambdat * d_fair(x - x_prev, delta) + return AhAx + spatial + temporal + +def gradient_descent_iteration(xk_1, gradf, step_size): + xk = xk_1 - (step_size * gradf(xk_1)) + return xk + +def gradient_descent_iteration_nesterov(iter, xk_1, gradf, step_size, nesterov_dict=None, max_iter=10): + if iter == 0: + nesterov_dict = { + 'yk_1': xk_1, + 'tk_1': 0 + } + elif iter >= max_iter: + eprint("not iterating, max_iter reached") + return [xk_1, nesterov_dict] + tk_1 = nesterov_dict['tk_1'] + yk_1 = nesterov_dict['yk_1'] + + tk = 0.5 * (1 + cp.sqrt((4 * cp.square(tk_1)) + 1)) + xk = yk_1 - (step_size * gradf(yk_1)) + yk = xk + (((tk_1 - 1) / tk) * (xk - xk_1)) + return [xk, {'yk_1': yk, 'tk_1': tk}] \ No newline at end of file diff --git a/toolboxes/nhlbi_gt_toolbox/utils/util_functions.cpp b/toolboxes/nhlbi_gt_toolbox/utils/util_functions.cpp index 50080f2..7924c2c 100644 --- a/toolboxes/nhlbi_gt_toolbox/utils/util_functions.cpp +++ b/toolboxes/nhlbi_gt_toolbox/utils/util_functions.cpp @@ -486,7 +486,10 @@ namespace nhlbi_toolbox imarray.meta_.resize(N*S*LOC); auto fov = recon_params.fov; - auto rmsize = recon_params.rmatrixSize; + auto rmsize = recon_params.rmatrixSize_scanner; + if (E0 != rmsize.x || E1!=rmsize.y || E2!=rmsize.z){ + GDEBUG_STREAM("WARNING Images don't have the expected dimension [E0 E1 E2] =[" << E0 << " " << E1 << " " << E2 <<"] != recon matrix [x y z] =[" << rmsize.x << " " << rmsize.y << " " << rmsize.z <<"]") + } for (size_t loc = 0; loc < LOC; loc++) { for (size_t s = 0; s < S; s++) { for (size_t n = 0; n < N; n++) { @@ -500,10 +503,12 @@ namespace nhlbi_toolbox imarray.headers_(n, s, loc).average = acqhdr.idx.average; imarray.headers_(n, s, loc).slice = acqhdr.idx.slice; imarray.headers_(n, s, loc).contrast = acqhdr.idx.contrast; - imarray.headers_(n, s, loc).phase = acqhdr.idx.phase; + imarray.headers_(n, s, loc).phase = n; imarray.headers_(n, s, loc).repetition = acqhdr.idx.repetition; imarray.headers_(n, s, loc).set = acqhdr.idx.set; imarray.headers_(n, s, loc).acquisition_time_stamp = acqhdr.acquisition_time_stamp; + //imarray.headers_(n, s, loc).physiology_time_stamp = acqhdr.physiology_time_stamp; + imarray.headers_(n, s, loc).physiology_time_stamp[0]=(uint32_t)(n*25); imarray.headers_(n, s, loc).position[0] = acqhdr.position[0]; imarray.headers_(n, s, loc).position[1] = acqhdr.position[1]; imarray.headers_(n, s, loc).position[2] = acqhdr.position[2]; @@ -1209,6 +1214,87 @@ namespace nhlbi_toolbox return (permute(tempDef, {0, 1, 3, 2})); } + template + cuNDArray crop_to_recon_params_dims(cuNDArray& input,reconParams recon_params) + { + cuNDArray output; + size_t NDim = input.get_number_of_dimensions(); + GDEBUG_STREAM("Input dimensions: "); + for (size_t i = 0; i < input.get_number_of_dimensions(); ++i) { + GDEBUG_STREAM("Dim " << i << ": " << input.get_size(i)); + } + GDEBUG_STREAM("Number of dimensions: " << NDim); + + auto mr_x=recon_params.rmatrixSize_scanner.x; + auto mr_y=recon_params.rmatrixSize_scanner.y; + auto mr_z=recon_params.rmatrixSize_scanner.z; + if (recon_params.rmatrixSize_scanner.z == 1) { + + switch (NDim) { + case 2: + output.create({mr_x, mr_y}); + crop(uint64d2((input.get_size(0) - mr_x) / 2, + (input.get_size(1) - mr_y) / 2), + uint64d2(mr_x, mr_y), input, output); + break; + + case 3: + output.create({mr_x, mr_y, input.get_size(2)}); + crop(uint64d3((input.get_size(0) - mr_x) / 2, + (input.get_size(1) - mr_y) / 2, + 0), + uint64d3(mr_x, mr_y, input.get_size(2)),input, output); + break; + + case 4: + output.create( + {mr_x, mr_y, input.get_size(2), input.get_size(3)}); + crop(uint64d4((input.get_size(0) - mr_x) / 2, + (input.get_size(1) - mr_y) / 2, + 0, 0), + uint64d4(mr_x, mr_y,input.get_size(2), input.get_size(3)),input, output); + break; + + default: + GDEBUG_STREAM("crop_to_recondims is not working, unknow number of dimensions " << NDim); + } + }else{ + switch (NDim) { + case 3: + output.create( + {mr_x, mr_y, mr_z}); + crop(uint64d3((input.get_size(0) - mr_x) / 2, + (input.get_size(1) - mr_y) / 2, + (input.get_size(2) - mr_z) / 2), + uint64d3(mr_x, mr_y, mr_z),input, output); + break; + + case 4: + output.create({mr_x, mr_y, + mr_z, input.get_size(3)}); + crop(uint64d4((input.get_size(0) - mr_x) / 2, + (input.get_size(1) - mr_y) / 2, + (input.get_size(2) - mr_z) / 2, 0), + uint64d4(mr_x, mr_y,mr_z, input.get_size(3)),input, output); + break; + + case 5: + output.create({mr_x, mr_y, + mr_z, input.get_size(3), input.get_size(4)}); + crop(uint64d5((input.get_size(0) - mr_x) / 2, + (input.get_size(1) - mr_y) / 2, + (input.get_size(2) - mr_z) / 2, 0, 0), + uint64d5(mr_x, mr_y,mr_z, input.get_size(3), input.get_size(4)),input, output); + break; + + default: + GDEBUG_STREAM("crop_to_recondims is not working, unknow number of dimensions " << NDim); + } + } + return output; + } + + template hoNDArray std_real(hoNDArray input, unsigned int dim); template hoNDArray std_complex(hoNDArray input, unsigned int dim); template hoNDArray> std_complex(hoNDArray> input, unsigned int dim); @@ -1235,6 +1321,12 @@ namespace nhlbi_toolbox template cuNDArray concat(std::vector> &arrays); template cuNDArray concat(std::vector> &arrays); + + + template cuNDArray crop_to_recon_params_dims(cuNDArray & input,reconParams recon_params); + template cuNDArray crop_to_recon_params_dims(cuNDArray & input,reconParams recon_params); + + // template hoNDArray concat(std::vector> &arrays); // template hoNDArray concat(std::vector> &arrays); // template hoNDArray concat(std::vector> &arrays); diff --git a/toolboxes/nhlbi_gt_toolbox/utils/util_functions.h b/toolboxes/nhlbi_gt_toolbox/utils/util_functions.h index 807ffaf..5549a85 100644 --- a/toolboxes/nhlbi_gt_toolbox/utils/util_functions.h +++ b/toolboxes/nhlbi_gt_toolbox/utils/util_functions.h @@ -137,7 +137,9 @@ namespace nhlbi_toolbox cuNDArray padDeformations(cuNDArray deformation, std::vector size_deformation); - + template + cuNDArraycrop_to_recon_params_dims(cuNDArray& input,reconParams recon_params); + constexpr double GAMMA = 4258.0; /* Hz/G */ void enable_peeraccess();