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/apps/gadgetron/connection/HeaderConnection.cpp b/apps/gadgetron/connection/HeaderConnection.cpp index 2dc1b89..0af9097 100644 --- a/apps/gadgetron/connection/HeaderConnection.cpp +++ b/apps/gadgetron/connection/HeaderConnection.cpp @@ -41,7 +41,7 @@ namespace { std::string raw_header(read_string_from_stream(stream)); ISMRMRD::IsmrmrdHeader header{}; - + //GDEBUG_STREAM("HEADER" << raw_header.c_str()); ISMRMRD::deserialize(raw_header.c_str(), header); //PD @@ -57,10 +57,17 @@ namespace { if (it->name == "tSequenceVariant") { GDEBUG("Hash bstar found %s found\n", it->value.c_str()); std::string hash_bstar=it->value.c_str(); - std::string bstar_folder_traj="/opt/data/bstar_traj/"; - // Loop through all .seq files in bstar_folder_traj - for (const auto& entry : std::filesystem::directory_iterator(bstar_folder_traj)) { - + // Search for .seq trajectory files in known locations + std::vector traj_search_paths = { + "/opt/data/bstar_traj/", + "/opt/nhlbi-integration-test/data/", + "/opt/code/gadgetron_lit/test/nhlbi_integration_tests/data/", + }; + bool seq_found = false; + for (const auto& bstar_folder_traj : traj_search_paths) { + if (!std::filesystem::exists(bstar_folder_traj)) continue; + for (const auto& entry : std::filesystem::recursive_directory_iterator(bstar_folder_traj)) { + if (entry.path().extension()==".seq"){ //Find the hash (Last line of the .seq file) std::ifstream fin; @@ -85,23 +92,23 @@ namespace { } } - std::string lastLine; + std::string lastLine; std::getline(fin,lastLine); // Read the current line fin.close(); // Comparison hash of the data with hash of the .seq file if(lastLine.find(hash_bstar) != std::string::npos){ GDEBUG_STREAM("Seq file found with hash"<< lastLine); - - std::string traj_h5=(bstar_folder_traj +std::string("traj_") + entry.path().stem().string() + std::string(".h5")); + + std::string traj_h5=(entry.path().parent_path().string() + std::string("/traj_") + entry.path().stem().string() + std::string(".h5")); GDEBUG_STREAM("Pulseq trajectory filepath " << traj_h5); - + std::string xml_config; std::string hdf5_in_group="/dataset"; std::shared_ptr ismrmrd_dataset= std::shared_ptr(new ISMRMRD::Dataset(traj_h5.c_str(), hdf5_in_group.c_str(), false)); ismrmrd_dataset->readHeader(xml_config); ISMRMRD::IsmrmrdHeader h_traj; - + ISMRMRD::deserialize(xml_config.c_str(), h_traj); // MODIFying Header matrixSize and FOV : Only reconSpace.matrixSize/FOV are correct in the trajectory file auto factor_0r = float((size_t(round(float(h_traj.encoding.front().reconSpace.matrixSize.x) / 32.0))) * 32.0) / float(size_t(h_traj.encoding.front().reconSpace.matrixSize.x)); @@ -109,10 +116,11 @@ namespace { auto factor_2r = float((size_t(round(float(h_traj.encoding.front().reconSpace.matrixSize.z) / 32.0))) * 32.0) / float(size_t(h_traj.encoding.front().reconSpace.matrixSize.z)); GDEBUG_STREAM("Recon Matrix traj : X " << h_traj.encoding.front().reconSpace.matrixSize.x << " Y " << h_traj.encoding.front().reconSpace.matrixSize.y << " Z " << h_traj.encoding.front().reconSpace.matrixSize.z); GDEBUG_STREAM("Factor traj : X " << factor_0r << " Y " << factor_1r << " Z " << factor_2r); - + factor_0r = factor_0r == 0 ? 1 : factor_0r; factor_1r = factor_1r == 0 ? 1 : factor_1r; factor_2r = factor_2r == 0 ? 1 : factor_2r; + GDEBUG_STREAM("Raw Header information :"); GDEBUG_STREAM("Encoded Matrix: X " << header.encoding.front().encodedSpace.matrixSize.x << " Y " << header.encoding.front().encodedSpace.matrixSize.y << " Z " << header.encoding.front().encodedSpace.matrixSize.z); GDEBUG_STREAM("Recon Matrix: X " << header.encoding.front().reconSpace.matrixSize.x << " Y " << header.encoding.front().reconSpace.matrixSize.y << " Z " << header.encoding.front().reconSpace.matrixSize.z); @@ -120,13 +128,13 @@ namespace { GDEBUG_STREAM("Recon FOV: X " << header.encoding.front().reconSpace.fieldOfView_mm.x << " Y " << header.encoding.front().reconSpace.fieldOfView_mm.y << " Z " << header.encoding.front().reconSpace.fieldOfView_mm.z); GDEBUG_STREAM("Encoding Limits: Encoded step 1 max " << header.encoding.at(0).encodingLimits.kspace_encoding_step_1.get().maximum << " Encoded step 2 max " << header.encoding.at(0).encodingLimits.kspace_encoding_step_2.get().maximum ); - header.encoding.front().encodedSpace.matrixSize.x=size_t(h_traj.encoding.front().reconSpace.matrixSize.x* factor_0r); - header.encoding.front().encodedSpace.matrixSize.y=size_t(h_traj.encoding.front().reconSpace.matrixSize.y* factor_1r); - header.encoding.front().encodedSpace.matrixSize.z=size_t(h_traj.encoding.front().reconSpace.matrixSize.z* factor_2r); + header.encoding.front().encodedSpace.matrixSize.x=size_t(h_traj.encoding.front().reconSpace.matrixSize.x * factor_0r); + header.encoding.front().encodedSpace.matrixSize.y=size_t(h_traj.encoding.front().reconSpace.matrixSize.y * factor_1r); + header.encoding.front().encodedSpace.matrixSize.z=size_t(h_traj.encoding.front().reconSpace.matrixSize.z * factor_2r); - header.encoding.front().reconSpace.matrixSize.x=size_t(h_traj.encoding.front().reconSpace.matrixSize.x* factor_0r); - header.encoding.front().reconSpace.matrixSize.y=size_t(h_traj.encoding.front().reconSpace.matrixSize.y* factor_1r); - header.encoding.front().reconSpace.matrixSize.z=size_t(h_traj.encoding.front().reconSpace.matrixSize.z* factor_2r); + header.encoding.front().reconSpace.matrixSize.x=size_t(h_traj.encoding.front().reconSpace.matrixSize.x * factor_0r); + header.encoding.front().reconSpace.matrixSize.y=size_t(h_traj.encoding.front().reconSpace.matrixSize.y * factor_1r); + header.encoding.front().reconSpace.matrixSize.z=size_t(h_traj.encoding.front().reconSpace.matrixSize.z * factor_2r); header.encoding.front().encodedSpace.fieldOfView_mm.x=h_traj.encoding.front().reconSpace.fieldOfView_mm.x; header.encoding.front().encodedSpace.fieldOfView_mm.y=h_traj.encoding.front().reconSpace.fieldOfView_mm.y; @@ -138,7 +146,7 @@ namespace { // Modifying the header encodingLimits : - + header.encoding.at(0).encodingLimits.kspace_encoding_step_1.get().maximum =h_traj.encoding.front().encodingLimits.kspace_encoding_step_1.get().maximum; header.encoding.at(0).encodingLimits.kspace_encoding_step_2.get().maximum =h_traj.encoding.front().encodingLimits.segment.get().maximum; @@ -149,13 +157,20 @@ namespace { GDEBUG_STREAM("Recon FOV: X " << header.encoding.front().reconSpace.fieldOfView_mm.x << " Y " << header.encoding.front().reconSpace.fieldOfView_mm.y << " Z " << header.encoding.front().reconSpace.fieldOfView_mm.z); GDEBUG_STREAM("Encoding Limits: Encoded step 1 max " << header.encoding.at(0).encodingLimits.kspace_encoding_step_1.get().maximum << " Encoded step 2 max " << header.encoding.at(0).encodingLimits.kspace_encoding_step_2.get().maximum ); + seq_found = true; break; } } - + } + if (seq_found) break; } + if (seq_found) break; + } + if (!seq_found) { + GDEBUG_STREAM("Seq file not found in any search path for hash " << hash_bstar); + } } } } 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/cases/mocolr_bSTAR.cfg b/test/nhlbi_integration_tests/cases/mocolr_bSTAR.cfg new file mode 100644 index 0000000..97acdd2 --- /dev/null +++ b/test/nhlbi_integration_tests/cases/mocolr_bSTAR.cfg @@ -0,0 +1,44 @@ +[dependency.siemens] +data_file = mocolr_bSTAR/noise_data.h5 +measurement = 0 +additional_arguments = skip_converstion + +[dependency.client] +configuration = default_measurement_dependencies.xml + +[reconstruction.siemens] +data_file = mocolr_bSTAR/recon_data.h5 +measurement = 0 +additional_arguments = skip_converstion + +[reconstruction.client] +configuration = pulmonary_echo0.xml + +[reconstruction.test.1] +reference_file = mocolr_bSTAR/baseline_output.h5 +reference_images = pulmonary_echo0.xml/image_6 +output_images = pulmonary_echo0.xml/image_6 +value_comparison_threshold = 0.01 +scale_comparison_threshold = 0.01 + +[reconstruction.test.2] +reference_file = mocolr_bSTAR/baseline_output.h5 +reference_images = pulmonary_echo0.xml/image_7 +output_images = pulmonary_echo0.xml/image_7 +value_comparison_threshold = 0.01 +scale_comparison_threshold = 0.01 + +[requirements] +system_memory = 8192 +gpu_support = 1 +gpu_memory = 8192 + +[tags] +tags = nhlbi,mocolr,bSTAR,lung + +[nhlbi] +description = MOCO Low Rank respiratory resolved pulmonary reconstruction +noise_file = mocolr_bSTAR/noise_data.h5 +additional_dataset_0 = mocolr_bSTAR/bstar_450mm_1.20mm_true0_TR2.32ms_rf200_i110_67k_FA40_WASP_self0_fid0_noise0.seq +additional_dataset_1 = mocolr_bSTAR/traj_bstar_450mm_1.20mm_true0_TR2.32ms_rf200_i110_67k_FA40_WASP_self0_fid0_noise0.h5 +baseline_recon_time = 179.6 \ No newline at end of file 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..4bfdf01 --- /dev/null +++ b/test/nhlbi_integration_tests/generate_baseline.py @@ -0,0 +1,527 @@ +#!/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) + + 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) + + # 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 = [ + 'dependency.siemens', + 'dependency.client', + '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..f6be8a3 --- /dev/null +++ b/test/nhlbi_integration_tests/nhlbi_data.json @@ -0,0 +1,114 @@ +[ + { + "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": "spiral_vds_dlwd/noise_data.h5", + "sha256": "89396f6de4bf248aa9a37335ac086b73dd499773ec1f57c9aba83879b1cc144d", + "type": "noise", + "test": "spiral_vds_dlwd" + }, + { + "file": "spiral_vds_dlwd/recon_data.h5", + "sha256": "4d1888c4714467cbc4a0f7ccb454275d10596ffcd64b003099fdf6e110e79abf", + "type": "input", + "test": "spiral_vds_dlwd" + }, + { + "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": "spiral_vds_dlwd/baseline_output.h5", + "sha256": "6251769fc250e16436d5274599dc859cbc8f0e936a9892ee218a81c593cc44e8", + "type": "baseline", + "test": "spiral_vds_dlwd", + "validated_by": "unknown", + "validated_date": "2026-03-19", + "git_sha": "d11d8915", + "recon_time_seconds": 236.2 + }, + { + "file": "mocolr_bSTAR/noise_data.h5", + "sha256": "0e635ec41daeae49f01788aeddc6e751e21ba7b52a12276929307acb26eef615", + "type": "noise", + "test": "mocolr_bSTAR" + }, + { + "file": "mocolr_bSTAR/recon_data.h5", + "sha256": "e2d2b684eba9c78c7c9bc1252a47bf080117a99f9705e9032e6ba200ca8b680c", + "type": "input", + "test": "mocolr_bSTAR" + }, + { + "file": "mocolr_bSTAR/bstar_450mm_1.20mm_true0_TR2.32ms_rf200_i110_67k_FA40_WASP_self0_fid0_noise0.seq", + "sha256": "9fed494c5103abd5e4abe7dceda46eb22eaf256a7decc9172807369ed599546d", + "type": "additional", + "test": "mocolr_bSTAR" + }, + { + "file": "mocolr_bSTAR/traj_bstar_450mm_1.20mm_true0_TR2.32ms_rf200_i110_67k_FA40_WASP_self0_fid0_noise0.h5", + "sha256": "25c3528406a8a2c04a452955c9ddb1f4483c79f29f276826f0a6f3aff5fab3ca", + "type": "additional", + "test": "mocolr_bSTAR" + }, + { + "file": "cardiovascular_bSTAR/noise_data.h5", + "sha256": "0e635ec41daeae49f01788aeddc6e751e21ba7b52a12276929307acb26eef615", + "type": "noise", + "test": "cardiovascular_bSTAR" + }, + { + "file": "cardiovascular_bSTAR/recon_data.h5", + "sha256": "e2d2b684eba9c78c7c9bc1252a47bf080117a99f9705e9032e6ba200ca8b680c", + "type": "input", + "test": "cardiovascular_bSTAR" + }, + { + "file": "cardiovascular_bSTAR/bstar_450mm_1.20mm_true0_TR2.32ms_rf200_i110_67k_FA40_WASP_self0_fid0_noise0.seq", + "sha256": "9fed494c5103abd5e4abe7dceda46eb22eaf256a7decc9172807369ed599546d", + "type": "additional", + "test": "cardiovascular_bSTAR" + }, + { + "file": "cardiovascular_bSTAR/traj_bstar_450mm_1.20mm_true0_TR2.32ms_rf200_i110_67k_FA40_WASP_self0_fid0_noise0.h5", + "sha256": "25c3528406a8a2c04a452955c9ddb1f4483c79f29f276826f0a6f3aff5fab3ca", + "type": "additional", + "test": "cardiovascular_bSTAR" + }, + { + "file": "cardiovascular_bSTAR/baseline_output.h5", + "sha256": "0e2dd48076fda2baee8faa98ebd716dace76f5c7f1d38aa8b9a700ad3c6eb2b9", + "type": "baseline", + "test": "cardiovascular_bSTAR", + "validated_by": "unknown", + "validated_date": "2026-04-08", + "git_sha": "315e7878", + "recon_time_seconds": 294.9 + }, + { + "file": "mocolr_bSTAR/baseline_output.h5", + "sha256": "d4aa67155c57272b8fca249d4b43f010a52929b4ab9b39e345d6afa577a43c0a", + "type": "baseline", + "test": "mocolr_bSTAR", + "validated_by": "unknown", + "validated_date": "2026-04-08", + "git_sha": "315e7878", + "recon_time_seconds": 179.6 + } +] 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..c1a5dab --- /dev/null +++ b/test/nhlbi_integration_tests/submit_test.py @@ -0,0 +1,249 @@ +#!/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=[]): + 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, + } + 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)') + + args = parser.parse_args() + print(args) + print(type(args.additional_files)) + print(len(args.additional_files)) + + # 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) + + # 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}") + + + + # 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)}") + 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, + }) + 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 + ) + 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..0ffcbcd 100644 --- a/toolboxes/nhlbi_gt_toolbox/CMakeLists.txt +++ b/toolboxes/nhlbi_gt_toolbox/CMakeLists.txt @@ -66,6 +66,8 @@ set(gadgetron_nhlbi_gt_toolbox_config_files config/pulmonary_MOCOLR.xml config/cardiovascular_iMOCO.xml config/cardiopulmonary_recon.xml + config/imoco_recon_vds.xml + config/pulmonary_echo0.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/pulmonary_echo0.xml b/toolboxes/nhlbi_gt_toolbox/config/pulmonary_echo0.xml new file mode 100644 index 0000000..30810d3 --- /dev/null +++ b/toolboxes/nhlbi_gt_toolbox/config/pulmonary_echo0.xml @@ -0,0 +1,168 @@ + + + + 2 + + + + gadgetron_core_readers + AcquisitionReader + + + gadgetron_core_readers + WaveformReader + + + gadgetron_core_readers + ImageReader + + + + + + gadgetron_core_writers + ImageWriter + + + + NoiseAdjustgadgetron_mricoreNoiseAdjustGadget + + + + + + + + + + + + + SelectEchoes + nhlbi_gt_gadgets + SelectEchoes + setnum0 + + + + + + + + + + + + + + + + + + + + + + + + + + + ImagetoVector + nhlbi_gt_gadgets + ImagetoVector + + + + gadgetron_mricore + PCACoilGadget + + + + gadgetron_mricore + CoilReductionGadget + + + + + + PrepreconParams + nhlbi_gt_gadgets + PrepreconParams + matOSP_vector1.0 1.0 0.8 + scannerOSP_vector1 1 0.8 + warpCUDA_vectortrue true true + downsampling_vector1 1 + is3Dtrue + kernel_width_dcf4 + oversampling_factor_dcf2.1 + kernel_width_dcf_avg3 + oversampling_factor_dcf_avg3 + lambda_spatial0.05 + lambda_spatial_imoco0.01 + lambda_time0.1 + lambda_time20 + lambda_LR0.1 + iteration_count_moco1 + iterations3 + iterations_imoco0 + doMC_iterfalse + tolerance1e-3 + selectedDevices_STR3 + try_channel_griddingfalse + + + + + Noncart_recon_gadget + nhlbi_gt_gadgets + Noncart_recon_gadget + Debug0 + doConcomitantFieldCorrectionfalse + referencePhase0.49 + estimateCSM_perc100 + reconType6 + binning_order0 1 2 + binning_collapse_to_lastfalse + series_counter_initial6 + save_avgtrue + save_intermediate_imagesfalse + + + + + + 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/doc/installation.rst b/toolboxes/nhlbi_gt_toolbox/doc/installation.rst index 51fd32d..2b0d196 100644 --- a/toolboxes/nhlbi_gt_toolbox/doc/installation.rst +++ b/toolboxes/nhlbi_gt_toolbox/doc/installation.rst @@ -19,8 +19,8 @@ First of all, you will install Gadgetron : Once built, the package can be used with gadgetron using the config xml files provided with this repository (`config files repository `_). -Validate installation -+++++++++++++++++++++ +Validate Gadgetron installation ++++++++++++++++++++++++++++++++ First, validate that the Gadgetron is installed and working.After activating the environment (with ``conda activate gadgetron``), the command ``gadgetron --info`` should give you information about your installed version of the Gadgetron and it would look something like this:: @@ -42,10 +42,64 @@ 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 first need to download the test data using the following command: + +.. code-block:: console + + conda activate gadgetron + python test/nhlbi_integration_tests/get_nhlbi_data.py download + +Then, you can run the following command to test the bSTAR pulmonary image reconstruction pipeline for example: + +.. code-block:: console + + cd test/nhlbi_integration_tests/ + python run_nhlbi_tests.py cases/mocolr_bSTAR.cfg -F + +The expected output of the test should look like this:: + + Downloading test data... + Verified: /opt/code/gadgetron/test/nhlbi_integration_tests/data/mocolr_bSTAR/bstar_450mm_1.20mm_true0_TR2.32ms_rf200_i110_67k_FA40_WASP_self0_fid0_noise0.seq + Verified: /opt/code/gadgetron/test/nhlbi_integration_tests/data/mocolr_bSTAR/noise_data.h5 + Verified: /opt/code/gadgetron/test/nhlbi_integration_tests/data/mocolr_bSTAR/baseline_output.h5 + Verified: /opt/code/gadgetron/test/nhlbi_integration_tests/data/mocolr_bSTAR/traj_bstar_450mm_1.20mm_true0_TR2.32ms_rf200_i110_67k_FA40_WASP_self0_fid0_noise0.h5 + Verified: /opt/code/gadgetron/test/nhlbi_integration_tests/data/mocolr_bSTAR/recon_data.h5 + Querying Gadgetron capabilities... + + Test 1 of 1: cases/mocolr_bSTAR.cfg + + Running Gadgetron test cases/mocolr_bSTAR.cfg with: + -- ISMRMRD_HOME : None + -- GADGETRON_HOME : None + -- TEST CASE : cases/mocolr_bSTAR.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/mocolr_bSTAR/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.13 s + Copying prepared ISMRMRD data: /opt/code/gadgetron/test/nhlbi_integration_tests/data/mocolr_bSTAR/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: 378.75 s + reconstruction.test.1 [OK] (Norm: 3.5e-05 [0.01] Scale: 6.0e-08 [0.01]) + reconstruction.test.1 [OK] (Output headers matched reference) + reconstruction.test.2 [OK] (Norm: 7.3e-04 [0.01] Scale: 3.0e-05 [0.01]) + reconstruction.test.2 [OK] (Output headers matched reference) + Test status: Passed + SPEED REGRESSION: 378.9s vs baseline 179.6s (+111.0%, threshold 50%) + + Speed regressions: + cases/mocolr_bSTAR.cfg + + 1 tests passed. 0 tests failed. 0 tests skipped. 0 missing baselines. 1 speed regressions. + Total processing time: 378.88 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 @@ -69,19 +123,20 @@ Once the docker container is running, you can start a bash terminal inside the c docker exec -ti cardio_pulmonary_bstar_rt 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 using the following command: .. 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 noise/noise_Freemax_XL_NIH_2025-03-06-112829_FID016823_bstar_1_20mm_FA40_FOV450.h5 -c default_measurement_dependencies.xml + gadgetron_ismrmrd_client -p 9002 -f h5/Freemax_XL_NIH_2025-03-06-112829_FID016823_bstar_1_20mm_FA40_FOV450.h5 -c pulmonary_echo0.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 cardio_pulmonary_bstar_rt 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,13 +144,15 @@ 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 noise/noise_Freemax_XL_NIH_2025-03-06-112829_FID016823_bstar_1_20mm_FA40_FOV450.h5 -c default_measurement_dependencies.xml + gadgetron_ismrmrd_client -a SERVER_ADDRESS -p 9063 -f h5/Freemax_XL_NIH_2025-03-06-112829_FID016823_bstar_1_20mm_FA40_FOV450.h5 -c pulmonary_echo0.xml -o OUTPUT_FILENAME.h5 Dataset ------- -The test data can be downloaded from zenodo: `18461603 `_ +The test data can also be downloaded from zenodo: `18461603 `_ .. note:: + More Information on Gadgetron are available over here : `Gadgetron repository `_ and `Gadgetron documentation `_ 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/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/AcquisitionWaveformFanout.h b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/AcquisitionWaveformFanout.h index 4c0b180..e5af8f2 100644 --- a/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/AcquisitionWaveformFanout.h +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/AcquisitionWaveformFanout.h @@ -9,4 +9,4 @@ using namespace Gadgetron; using namespace Gadgetron::Core; -using AcquisitionWaveformFanout = Gadgetron::Core::Parallel::Fanout>; +using AcquisitionWaveformFanout = Gadgetron::Core::Parallel::Fanout>>>>>; diff --git a/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.cpp b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.cpp index 6005e8a..08e51a5 100644 --- a/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.cpp +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.cpp @@ -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:" < 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 +380,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 +496,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 +645,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 +694,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..c525720 100644 --- a/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.h +++ b/toolboxes/nhlbi_gt_toolbox/gadgets/waveforms/WaveformToTrajectory.h @@ -37,20 +37,21 @@ using namespace Gadgetron; std::map> trajectory_map; std::map gradient_wave_store; size_t curAvg=0; - + protected: ISMRMRD::IsmrmrdHeader header; 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..5b4ae7c 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.cpp +++ b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction.cpp @@ -450,7 +450,6 @@ void noncartesian_reconstruction::reconstruct(cuNDArray* data // 3) * std::pow(recon_params.oversampling_factor_, D) * 4 + (stride_results * 4) * (2 * CHA) * // std::pow(recon_params.oversampling_factor_, D) * 4) / float(std::pow(1024, 3)); GDEBUG_STREAM("data and image // space: " << float(data_and_imageSize)); // this is not working - if(!recon_params.try_channel_gridding) { for (int iCHA = 0; iCHA < CHA; iCHA++) { @@ -459,12 +458,10 @@ void noncartesian_reconstruction::reconstruct(cuNDArray* data this->nfft_plan_->compute(data_view, results_view, dcw, NFFT_comp_mode::BACKWARDS_NC2C); } - cudaDeviceSynchronize(); + }else{ - try { this->nfft_plan_->compute(*data, *image, dcw, NFFT_comp_mode::BACKWARDS_NC2C); - cudaDeviceSynchronize(); } catch (const std::exception& e) { cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { @@ -476,10 +473,8 @@ void noncartesian_reconstruction::reconstruct(cuNDArray* data this->nfft_plan_->compute(data_view, results_view, dcw, NFFT_comp_mode::BACKWARDS_NC2C); } - cudaDeviceSynchronize(); this->recon_params.try_channel_gridding=false; GDEBUG_STREAM("Try Channel gridding to false "); - cudaDeviceSynchronize(); } catch (...) { cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { @@ -490,10 +485,11 @@ void noncartesian_reconstruction::reconstruct(cuNDArray* data auto results_view = cuNDArray>(image_dimensions, image->data() + stride_results * iCHA); this->nfft_plan_->compute(data_view, results_view, dcw, NFFT_comp_mode::BACKWARDS_NC2C); } - cudaDeviceSynchronize(); this->recon_params.try_channel_gridding=false; + GDEBUG_STREAM("Try Channel gridding to false "); } } + cudaDeviceSynchronize(); } @@ -516,6 +512,7 @@ void noncartesian_reconstruction::reconstruct(cuNDArray* data auto out_dimensions = *image->get_dimensions(); auto in_dimensions = *data->get_dimensions(); + auto channel_dimensions = *data->get_dimensions(); if (CHA != 1 || csm->get_size(csm->get_number_of_dimensions() - 1) == CHA) { in_dimensions.pop_back(); // remove CHA @@ -536,19 +533,68 @@ void noncartesian_reconstruction::reconstruct(cuNDArray* data //#pragma omp target teams num_teams(numteams) //#pragma omp distribute parallel for reduction(complex_add: test) - for (size_t ich = 0; ich < CHA; ich++) { + if(!recon_params.try_channel_gridding) + { + for (size_t ich = 0; ich < CHA; ich++) { + + auto slice_view=cuNDArray(in_dimensions, data->get_data_ptr() + stride_ch * ich); + auto tmpview = cuNDArray(out_dimensions); - auto slice_view=cuNDArray(in_dimensions, data->get_data_ptr() + stride_ch * ich); - auto tmpview = cuNDArray(out_dimensions); + this->nfft_plan_->compute(&slice_view, tmpview, dcw, NFFT_comp_mode::BACKWARDS_NC2C); - this->nfft_plan_->compute(&slice_view, tmpview, dcw, NFFT_comp_mode::BACKWARDS_NC2C); + auto csm_view = cuNDArray(out_dimensions, csm->get_data_ptr() + stride_out * ich); + tmpview *= *conj(&csm_view); + out_view_ch += tmpview; + } + + }else{ + try { + out_dimensions.push_back(CHA); + auto channel_images=cuNDArray(out_dimensions); + this->nfft_plan_->compute(*data, channel_images, dcw, NFFT_comp_mode::BACKWARDS_NC2C); + channel_images *= *conj(csm); + out_view_ch += *sum(&channel_images, channel_images.get_number_of_dimensions() - 1); - auto csm_view = cuNDArray(out_dimensions, csm->get_data_ptr() + stride_out * ich); - tmpview *= *conj(&csm_view); - out_view_ch += tmpview; + } catch (const std::exception& e) { + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + GERROR_STREAM("CUDA error in deconstruct: " << cudaGetErrorString(err)); + } + for (size_t ich = 0; ich < CHA; ich++) { + auto slice_view=cuNDArray(in_dimensions, data->get_data_ptr() + stride_ch * ich); + auto tmpview = cuNDArray(out_dimensions); + + this->nfft_plan_->compute(&slice_view, tmpview, dcw, NFFT_comp_mode::BACKWARDS_NC2C); + + auto csm_view = cuNDArray(out_dimensions, csm->get_data_ptr() + stride_out * ich); + tmpview *= *conj(&csm_view); + out_view_ch += tmpview; + } + this->recon_params.try_channel_gridding=false; + GDEBUG_STREAM("Try Channel gridding to false "); + } catch (...) { + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + GERROR_STREAM("Unknown CUDA error in deconstruct: " << cudaGetErrorString(err)); + } + for (size_t ich = 0; ich < CHA; ich++) { + + + auto slice_view=cuNDArray(in_dimensions, data->get_data_ptr() + stride_ch * ich); + auto tmpview = cuNDArray(out_dimensions); + + this->nfft_plan_->compute(&slice_view, tmpview, dcw, NFFT_comp_mode::BACKWARDS_NC2C); + + auto csm_view = cuNDArray(out_dimensions, csm->get_data_ptr() + stride_out * ich); + tmpview *= *conj(&csm_view); + out_view_ch += tmpview; + } + this->recon_params.try_channel_gridding=false; + GDEBUG_STREAM("Try Channel gridding to false "); + } } - cudaDeviceSynchronize(); + cudaDeviceSynchronize(); } template @@ -567,7 +613,6 @@ void noncartesian_reconstruction::deconstruct(cuNDArray* imag // GDEBUG_STREAM("Y: " << images->get_size(1)); // GDEBUG_STREAM("Z: " << images->get_size(2)); // GDEBUG_STREAM("C: " << images->get_size(3)); - // if (!this->isprocessed) { this->nfft_plan_->preprocess(*traj, NFFT_prep_mode::C2NC); @@ -612,7 +657,6 @@ void noncartesian_reconstruction::deconstruct(cuNDArray* imag cudaGetDeviceProperties(&properties, data->get_device()); // cudaSetDevice(data->get_device()); - // GDEBUG_STREAM("Failed: now running in slower channel by channel mode"); if(!recon_params.try_channel_gridding) { for (int iCHA = 0; iCHA < CHA; iCHA++) { @@ -623,15 +667,13 @@ void noncartesian_reconstruction::deconstruct(cuNDArray* imag tmp_view *= csm_view; this->nfft_plan_->compute(&tmp_view, data_view, dcw, NFFT_comp_mode::FORWARDS_C2NC); } - cudaDeviceSynchronize(); + }else{ - try { auto dims_csm = csm->get_dimensions(); cuNDArray images_mult_csm(dims_csm); csm_mult_M(images, &images_mult_csm, csm); this->nfft_plan_->compute(images_mult_csm, *data, dcw, NFFT_comp_mode::FORWARDS_C2NC); - cudaDeviceSynchronize(); } catch (const std::exception& e) { cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { @@ -647,7 +689,6 @@ void noncartesian_reconstruction::deconstruct(cuNDArray* imag } this->recon_params.try_channel_gridding=false; GDEBUG_STREAM("Try Channel gridding to false "); - cudaDeviceSynchronize(); } catch (...) { cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { @@ -661,10 +702,10 @@ void noncartesian_reconstruction::deconstruct(cuNDArray* imag tmp_view *= csm_view; this->nfft_plan_->compute(&tmp_view, data_view, dcw, NFFT_comp_mode::FORWARDS_C2NC); } - cudaDeviceSynchronize(); this->recon_params.try_channel_gridding=false; } } + cudaDeviceSynchronize(); // this->nfft_plan_->compute(*images, *data, dcw, NFFT_comp_mode::FORWARDS_C2NC); } @@ -891,7 +932,7 @@ cuNDArray noncartesian_reconstruction::estimate_dcf(cuNDArrayget_dimensions()); std::vector flat_dims = {traj->get_number_of_elements()}; auto hoTraj = hoNDArray>( @@ -920,7 +961,7 @@ cuNDArray noncartesian_reconstruction::estimate_dcf(cuNDArray* dcf_in) { float kw_dcf = recon_params.kernel_width_dcf_; // 1e-2 float osf_dcf = recon_params.oversampling_factor_dcf_; // 1.5 - GDEBUG_STREAM("DCF parameters: kw " << kw_dcf << " os " << osf_dcf); + //GDEBUG_STREAM("DCF parameters: kw " << kw_dcf << " os " << osf_dcf); auto dims_traj = *(traj->get_dimensions()); auto hoTraj = hoNDArray>( diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_2Dt.cpp b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_2Dt.cpp index d47bcb6..677f4f4 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_2Dt.cpp +++ b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_2Dt.cpp @@ -91,8 +91,8 @@ namespace nhlbi_toolbox solver_.add_regularization_operator(Ry, recon_params.norm); GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); reg_image = *solver_.solve(data); auto reg_image_dims = *reg_image.get_dimensions(); @@ -175,8 +175,8 @@ namespace nhlbi_toolbox GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); reg_image = *solver_.solve(data); cuNDArray images_cropped = this->crop_to_recondims(reg_image); diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_3D.cpp b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_3D.cpp index f93e8d3..47d1b2b 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_3D.cpp +++ b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_3D.cpp @@ -86,8 +86,8 @@ namespace nhlbi_toolbox solver_.add_regularization_operator(Rz, recon_params.norm); GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); reg_image = *solver_.solve(data); cuNDArray images_cropped = this->crop_to_recondims(reg_image); @@ -236,7 +236,8 @@ namespace nhlbi_toolbox std::replace(ho_prereconi.begin(),ho_prereconi.end(),INFINITY,0.0f); auto ho_prereconri = *real_imag_to_complex(&ho_prereconr,&ho_prereconi); precon_weights = boost::make_shared>(hoNDArray(ho_prereconri)); - + _precon_weights->clear(); + _precon_weights_cropped.clear(); D_->set_weights(precon_weights); // setup solver spit-bergman @@ -265,9 +266,18 @@ namespace nhlbi_toolbox boost::shared_ptr> csm) { auto data_dims = *data->get_dimensions(); + auto stride = std::accumulate(data_dims.begin(), data_dims.end() - 1, size_t(1), std::multiplies()); + // prep data and dcw - doing this data save in memory to prevent data from being affected by recon. + cudaSetDevice(data->get_device()); + hoNDArray hodata(*data->get_dimensions()); + cudaMemcpy(hodata.get_data_ptr(), data->get_data_ptr(), data->get_number_of_elements() * sizeof(float_complext), cudaMemcpyDeviceToHost); + auto dcwPtr = boost::make_shared>(*dcw); - // need to multiply by the weights to correctly to the FWD transform because we did sqrt of dcw - *data *= *dcw; + for (auto iCHA = 0; iCHA < recon_params.numberChannels; iCHA++) + { + auto dataview = cuNDArray>((*dcw).get_dimensions(), data->data() + stride * iCHA); + dataview *= (*dcw); + } auto E_ = boost::shared_ptr>(new cuNonCartesianSenseOperator(ConvolutionType::ATOMIC)); auto D_ = boost::shared_ptr>(new cuCgPreconditioner()); @@ -319,8 +329,7 @@ namespace nhlbi_toolbox reg_image = *solver_.solve(data); cuNDArray images_cropped = this->crop_to_recondims(reg_image); - // de-prep data - *data /= *dcw; + cudaMemcpy(data->get_data_ptr(), hodata.get_data_ptr(), data->get_number_of_elements() * sizeof(float_complext), cudaMemcpyHostToDevice); return images_cropped; } diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_4D.cpp b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_4D.cpp index 3a45590..e63856d 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_4D.cpp +++ b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_4D.cpp @@ -145,8 +145,8 @@ namespace nhlbi_toolbox // gpus_input_possible.erase(std::remove(gpus_input_possible.begin(), gpus_input_possible.end(), data->get_device()), gpus_input_possible.end()); GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); auto reg_image = *solver_.solve(data); @@ -306,8 +306,8 @@ namespace nhlbi_toolbox // gpus_input_possible.erase(std::remove(gpus_input_possible.begin(), gpus_input_possible.end(), data->get_device()), gpus_input_possible.end()); GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); auto reg_image = *solver_.solve(data); @@ -607,8 +607,8 @@ namespace nhlbi_toolbox solver_.add_regularization_operator(Rz, recon_params.norm); GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); reg_image = *solver_.solve(data); @@ -855,7 +855,7 @@ namespace nhlbi_toolbox */ GDEBUG_STREAM("Data_device:" << data->get_device()); GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); @@ -995,8 +995,8 @@ namespace nhlbi_toolbox solver_.add_regularization_operator(Rz, recon_params.norm); */ GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); auto reg_image = *solver_.solve(data); @@ -1139,8 +1139,8 @@ namespace nhlbi_toolbox solver_.add_regularization_operator(Rz, recon_params.norm); */ GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); auto reg_image = *solver_.solve(data); @@ -1341,8 +1341,8 @@ namespace nhlbi_toolbox // solver_.add_group(recon_params.norm); GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); reg_image = *solver_.solve(data); @@ -1551,8 +1551,8 @@ namespace nhlbi_toolbox // solver_.add_regularization_group_operator(Rz); // solver_.add_group(recon_params.norm); GDEBUG_STREAM("Data_device:" << data->get_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); reg_image = *solver_.solve(data); diff --git a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_5D.cpp b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_5D.cpp index f4a326e..6a33786 100644 --- a/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_5D.cpp +++ b/toolboxes/nhlbi_gt_toolbox/spiral/reconstruction/noncartesian_reconstruction_5D.cpp @@ -124,8 +124,8 @@ precon_weights = boost::make_shared>(padget_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); reg_image = *solver_.solve(data); @@ -463,8 +463,8 @@ precon_weights = boost::make_shared>(padget_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); reg_image = *solver_.solve(data); @@ -701,8 +701,8 @@ precon_weights = boost::make_shared>(padget_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); reg_image = *solver_.solve(data); @@ -930,8 +930,8 @@ precon_weights = boost::make_shared>(padget_device()); - GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices[0]<< " [1] if exist" << recon_params.selectedDevices[1]); - solver_.set_gpus(recon_params.selectedDevices); + GDEBUG_STREAM("gpus_input_possible[0]:" << recon_params.selectedDevices_solver[0]<< " [1] if exist" << recon_params.selectedDevices_solver[1]); + solver_.set_gpus(recon_params.selectedDevices_solver); cudaSetDevice(data->get_device()); reg_image = *solver_.solve(data); @@ -964,4 +964,4 @@ precon_weights = boost::make_shared>(pad 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/gpu/gpuSVD.cu b/toolboxes/nhlbi_gt_toolbox/utils/gpu/gpuSVD.cu index a3b20e0..5932228 100644 --- a/toolboxes/nhlbi_gt_toolbox/utils/gpu/gpuSVD.cu +++ b/toolboxes/nhlbi_gt_toolbox/utils/gpu/gpuSVD.cu @@ -136,8 +136,8 @@ std::tuple,cuNDArray,cuNDArray> gpuSVD::cuda_DNSg float ms=0.f; CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop)); int info_h=0; CUDA_CHECK(cudaMemcpy(&info_h, d_info, sizeof(int), cudaMemcpyDeviceToHost)); - if (info_h != 0) std::cerr << "cuSOLVER Sgesvd info=" << info_h << std::endl; - std::cout << "GPU cuSOLVER gesvd (float) time: " << ms << " ms" << std::endl; + if (info_h != 0) GERROR_STREAM("cuSOLVER Sgesvd info=" << info_h); + GDEBUG_STREAM("GPU cuSOLVER gesvd (float) time: " << ms << " ms" ); CUDA_CHECK(cudaFree(d_work)); CUDA_CHECK(cudaFree(d_info)); CUSOLVER_CHECK(cusolverDnDestroy(handle)); @@ -196,8 +196,8 @@ std::tuple,cuNDArray,cuNDArray> float ms=0.f; CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop)); int info_h=0; CUDA_CHECK(cudaMemcpy(&info_h, d_info, sizeof(int), cudaMemcpyDeviceToHost)); - if (info_h != 0) std::cerr << "cuSOLVER Sgesvd info=" << info_h << std::endl; - std::cout << "GPU cuSOLVER gesvd (float) time: " << ms << " ms" << std::endl; + if (info_h != 0) GERROR_STREAM("cuSOLVER Sgesvd info=" << info_h); + GDEBUG_STREAM("GPU cuSOLVER gesvd (float) time: " << ms << " ms" ); d_work.clear(); CUDA_CHECK(cudaFree(d_info)); CUSOLVER_CHECK(cusolverDnDestroy(handle)); @@ -265,8 +265,8 @@ std::tuple,cuNDArray,cuNDArray> gpuSVD::cuda_DNSg float ms=0.f; CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop)); int info_h=0; CUDA_CHECK(cudaMemcpy(&info_h, d_info, sizeof(int), cudaMemcpyDeviceToHost)); - if (info_h != 0) std::cerr << "cuSOLVER Sgesvd info=" << info_h << std::endl; - std::cout << "GPU cuSOLVER gesvd (float) time: " << ms << " ms" << std::endl; + if (info_h != 0) GERROR_STREAM("cuSOLVER Sgesvd info=" << info_h); + GDEBUG_STREAM("GPU cuSOLVER gesvd (float) time: " << ms << " ms" ); CUSOLVER_CHECK(cusolverDnDestroyGesvdjInfo(params)); d_work.clear(); CUDA_CHECK(cudaFree(d_info)); @@ -406,8 +406,8 @@ std::tuple,cuNDArray,cuNDArray> float ms=0.f; CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop)); int info_h=0; CUDA_CHECK(cudaMemcpy(&info_h, d_info, sizeof(int), cudaMemcpyDeviceToHost)); - if (info_h != 0) std::cerr << "cuSOLVER Sgesvd info=" << info_h << std::endl; - std::cout << "GPU cuSOLVER gesvd (float) time: " << ms << " ms" << std::endl; + if (info_h != 0) GERROR_STREAM("cuSOLVER Sgesvd info=" << info_h); + GDEBUG_STREAM("GPU cuSOLVER gesvd (float) time: " << ms << " ms" ); CUSOLVER_CHECK(cusolverDnDestroyGesvdjInfo(params)); d_work.clear(); CUDA_CHECK(cudaFree(d_info)); diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/CardiacRespiGadget.py b/toolboxes/nhlbi_gt_toolbox/utils/python/CardiacRespiGadget.py index f93c3d4..6b5d736 100644 --- a/toolboxes/nhlbi_gt_toolbox/utils/python/CardiacRespiGadget.py +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/CardiacRespiGadget.py @@ -71,17 +71,18 @@ def CardiacRespiGadget(connection): "C_smoothing": True, "C_numBins_to_ms":False, "C_waveforms":False, + "RC_enforced_1_echo":False, } BPfilter_freqs= [0.08,0.1,0.45,0.50] boolean_keys=['phantom','gaussian','bstar','useDC','C_PHYSIO','R_stableBinning','R_evenbins','R_bidirectional','R_angular_filteration','R_binningPercent', - 'C_evenbins','C_arrythmia_detection','C_angular_filteration','C_HRinfo','C_even_timing','C_smoothing','C_numBins_to_ms','C_stableBinning','C_waveforms'] + 'C_evenbins','C_arrythmia_detection','C_angular_filteration','C_HRinfo','C_even_timing','C_smoothing','C_numBins_to_ms','C_stableBinning','C_waveforms','RC_enforced_1_echo'] str_keys=[] int_keys=['R_numBins','C_numBins','R_binningPercent','C_binningPercent','samples'] params=read_params(params_init,params_ref=params,boolean_keys=boolean_keys,str_keys=str_keys,int_keys=int_keys) if params['C_PHYSIO'] and params['C_waveforms']: connection.filter(lambda input: type(input)==mrd.Acquisition or type(input)==mrd.Waveform) - print("Receiving Waveforms") + eprint("Receiving Waveforms") else: connection.filter(lambda input: type(input)==mrd.Acquisition) gaussian_flag=params['gaussian'] @@ -103,6 +104,10 @@ def CardiacRespiGadget(connection): encoding_limits = mrd_header.encoding[0].encodingLimits number_of_sets=encoding_limits.set.maximum+1 + # Enforcing the binning one echoe time + if params['RC_enforced_1_echo']: + number_of_sets=1 + mz=mrd_header.encoding[0].encodedSpace.matrixSize.z ecg_data = [] ecg_tstamp = [] @@ -181,7 +186,7 @@ def CardiacRespiGadget(connection): if gaussian_flag: - print("3D cine: 1 sample every 3 samples") + eprint("3D cine: 1 sample every 3 samples") nav_data=nav_data[:,1::3,:] nav_tstamp=nav_tstamp[1::3] @@ -211,11 +216,11 @@ def CardiacRespiGadget(connection): #Respiratory binning - print (f"Respiratory Binning Nbins {numRBins} stable {params['R_stableBinning']} bidirectionnal {params['R_bidirectional']}") - print(f'DATA {nav_data_copy.shape}') - print(f'Tstamp {nav_tstamp_copy.shape}') + eprint (f"Respiratory Binning Nbins {numRBins} stable {params['R_stableBinning']} bidirectionnal {params['R_bidirectional']}") + eprint(f'DATA {nav_data_copy.shape}') + eprint(f'Tstamp {nav_tstamp_copy.shape}') if numRBins==1 and params["R_stableBinning"]==False: - print("No respiratory binning") + eprint("No respiratory binning") idx_acceptedTimes=[np.arange(len(nav_tstamp_copy))] acceptedTimes=[nav_tstamp_copy.squeeze()] else: @@ -228,7 +233,7 @@ def CardiacRespiGadget(connection): eprint('Execution time (GatingSignal):', elapsed_time, 'seconds') acceptedTimes,idx_acceptedTimes = binning(respiratory_waveform,nav_tstamp_copy.get(),params['R_binningPercent'],params['R_bidirectional'], params['R_stableBinning'], params['R_evenbins'], numRBins) - if bstar_flag and params['samples']==1: + if (bstar_flag and params['samples']==1) or (bstar_flag and params['samples']==2): resp_bins_index=[] for idx_a in idx_acceptedTimes: idx_a.sort() @@ -237,9 +242,13 @@ def CardiacRespiGadget(connection): resp_bins_index,maxSize=get_idx_to_send(acq_tstamp,acceptedTimes, samplingTime) # Cardiac binning - print (f"Cardiac Binning Nbins {numCBins} PHYSIO {params['C_PHYSIO']} Nbinorms {params['C_numBins_to_ms']}") + eprint (f"Cardiac Binning Nbins {numCBins} PHYSIO {params['C_PHYSIO']} Nbinorms {params['C_numBins_to_ms']}") if params['C_PHYSIO']: nav_tstamp=cp.array(ecg_tstamp) + eprint("nav_tstamp before",nav_tstamp.shape) + if bstar_flag: + nav_tstamp=nav_tstamp[::params['samples']] + eprint("nav_tstamp after",nav_tstamp.shape) if params["C_waveforms"]: waveform_np=np.concatenate(waveform_data,1) waveform_t_np = np.concatenate(waveform_timestamp,0) @@ -261,14 +270,19 @@ def CardiacRespiGadget(connection): cardiac_waveform_smooth = np.interp(2.5*nav_tstamp.get(), waveform_t_np, ecgtrigger)[None,:] else: cardiac_waveform_smooth=np.array(ecg_data)[:,:,0] + eprint("cardiac_waveform_smooth before",cardiac_waveform_smooth.shape) + if bstar_flag: + cardiac_waveform_smooth=cardiac_waveform_smooth[::params['samples'],:] + eprint("cardiac_waveform_smooth after",cardiac_waveform_smooth.shape) + if ((numCBins==1 and params['C_numBins_to_ms']==False) and params["C_stableBinning"]==False): - print("No Cardiac binning") + eprint("No Cardiac binning") bins_index=[np.arange(len(nav_tstamp.squeeze()))] C_acceptedTimes=[nav_tstamp.squeeze()] else : if (params["C_stableBinning"]==True and numCBins>1): - print("Stable Binning : Modifying number of cardiac bins") + eprint("Stable Binning : Modifying number of cardiac bins") numCBins=1 params['C_numBins_to_ms']=False @@ -295,7 +309,7 @@ def CardiacRespiGadget(connection): cardiac_waveform_smooth = cardiac_waveform bins_index,ecg_freq_final=cardiacbinning(cardiac_waveform_smooth,samplingTime,numCBins,ecg_freq,evenbins=params['C_evenbins'],phantomflag=params['phantom'],arrythmia_detection=params['C_arrythmia_detection'],even_timing=params["C_even_timing"],stable_binning=params['C_stableBinning'],stable_perc=params['C_binningPercent']) C_Index=np.arange(len(nav_tstamp.squeeze())) - print(len(C_Index)) + eprint(len(C_Index)) C_acceptedTimes=[] for set in range(number_of_sets): for nbin in range(len(bins_index)): @@ -329,15 +343,15 @@ def CardiacRespiGadget(connection): imageSize = maxSize #pow(2,math.ceil(math.log2(math.sqrt(maxSize)))) max_nchannels=(np.power(2,16)-1) #Bug nchannels in Image Header is a np.uint16 (max value =65535) if np.prod([imageSize, numRBins,numCBins,number_of_sets])==imageSize and imageSize>max_nchannels : - print("Carefull too much data in each bins, required to collapse Respiratory dimension !!!") + eprint("Carefull too much data in each bins, required to collapse Respiratory dimension !!!") numRBins=int(np.ceil(imageSize/(max_nchannels-1))) maxSize=max_nchannels-1 idx_0=idx_to_send_2D[0][1:] idx_to_send_2D = [] for idx_r in range(numRBins): - print(idx_r) - print(idx_r*max_nchannels) - print(np.min([max_nchannels*(idx_r+1),len(idx_0)])) + eprint(idx_r) + eprint(idx_r*max_nchannels) + eprint(np.min([max_nchannels*(idx_r+1),len(idx_0)])) tmp_idx=idx_0[idx_r*maxSize:np.min([maxSize*(idx_r+1),len(idx_0)])] idx_to_send_2D.append( np.concatenate( ([tmp_idx.shape[0]],tmp_idx))) imageSize=max_nchannels @@ -356,6 +370,4 @@ def CardiacRespiGadget(connection): if __name__ == '__main__': - gadgetron.external.listen(2000,CardiacRespiGadget) - - \ No newline at end of file + gadgetron.external.listen(2000,CardiacRespiGadget) \ No newline at end of file diff --git a/toolboxes/nhlbi_gt_toolbox/utils/python/Pulseq_WaveformToTrajectory.py b/toolboxes/nhlbi_gt_toolbox/utils/python/Pulseq_WaveformToTrajectory.py index c934b36..46be454 100644 --- a/toolboxes/nhlbi_gt_toolbox/utils/python/Pulseq_WaveformToTrajectory.py +++ b/toolboxes/nhlbi_gt_toolbox/utils/python/Pulseq_WaveformToTrajectory.py @@ -209,7 +209,7 @@ def readGIRFKernel(girf_path): girfz_data=np.reshape(girfz_data,[sGIRF,2]) girfz=girfz_data[:,0]+1j*girfz_data[:,1] - print('GIRF not swap') + eprint('GIRF not swap') #GIRF=np.column_stack((girfy,girfx,girfz)) # Swap x,y Needs investigation GIRF=np.column_stack((girfx,girfy,girfz)) # Swap x,y Needs investigation return GIRF,dtGIRF @@ -431,23 +431,38 @@ def Pulseq_WaveformtoTrajectoryGadget(connection): reading_time=time() trj_file=params['traj_folder_or_file'] - print(trj_file) + # Search for .seq files in known test data locations + trj_folder_test_candidates=[ + '/opt/nhlbi-integration-test/data', # RT container + '/opt/code/gadgetron/test/nhlbi_integration_tests/data', # Docker build + '/opt/code/gadgetron_lit/test/nhlbi_integration_tests/data', # Dev container + ] + seq_files_test=[] + for trj_folder_test in trj_folder_test_candidates: + if op.isdir(trj_folder_test): + seq_files_test=glob.glob(op.join(trj_folder_test,'*','*.seq')) + break + seq_files_test.sort() + if not op.exists(trj_file): + eprint(f"This trajectory file or folder does not exist : {trj_file}") + trj_file=trj_folder_test # Detect Traj_file.h5 based on the hash if op.isdir(trj_file): hash_traj='' for k in range(len(mrd_header.userParameters.userParameterString)): if mrd_header.userParameters.userParameterString[k].name =='tSequenceVariant': hash_traj=mrd_header.userParameters.userParameterString[k].value - print(hash_traj) + break + eprint(f"Hash trajectory :{hash_traj}") if hash_traj: - seq_files=glob.glob(op.join(trj_file,'*.seq')) + seq_files=glob.glob(op.join(trj_file,'*.seq'))+seq_files_test seq_files.sort() for seq_file in seq_files: hash_seq=read_n_to_last_line(seq_file).split(' ')[-1][:-1] if hash_traj==hash_seq: - print(f"Seq file found : {op.basename(seq_file)[:-4]}") + eprint(f"Seq file found : {op.basename(seq_file)[:-4]}") break - trj_file=glob.glob(op.join(trj_file,f"*{op.basename(seq_file)[:-4]}*.h5"))[0] + trj_file=glob.glob(op.join(op.dirname(seq_file),f"*{op.basename(seq_file)[:-4]}*.h5"))[0] with mrd.File(trj_file,'r') as mrd_file: traj_header=mrd_file['dataset'].header @@ -455,7 +470,7 @@ def Pulseq_WaveformtoTrajectoryGadget(connection): reading_time_traj=time() - print(' Reading traj running time %f s'%(reading_time_traj-reading_time)) + eprint(' Reading traj running time %f s'%(reading_time_traj-reading_time)) traj_unscaled = rearrange(np.array([tr.traj for tr in traj_acq if tr.flags==0]).squeeze(),'INT RO DIM -> RO INT DIM') ## ANGLES INFORMATION @@ -480,7 +495,7 @@ def Pulseq_WaveformtoTrajectoryGadget(connection): discard_pre=traj_acq[idx_ref].discard_pre discard_post=traj_acq[idx_ref].discard_post real_dwell_time=traj_acq[idx_ref].sample_time_us*1e-6 - print(real_dwell_time) + eprint(real_dwell_time) nr_readouts=traj_header.encoding[0].encodingLimits.kspace_encoding_step_1.maximum+1 nr_interleaves=traj_header.encoding[0].encodingLimits.segment.maximum+1 @@ -521,7 +536,10 @@ def Pulseq_WaveformtoTrajectoryGadget(connection): traj_unscaled[...,ax] *= 0.5 / max_xyz[ax] eprint(f"shape traj {traj_unscaled.shape}") - traj_pulseq=traj_unscaled + if not params['applyGIRF']: + traj_pulseq=traj_unscaled + traj_echo1=traj_unscaled[:int(traj_unscaled.shape[0]/2),...] + traj_echo2=traj_unscaled[-int(traj_unscaled.shape[0]/2):,...] FOV=traj_header.encoding[0].reconSpace.fieldOfView_mm matrixR=traj_header.encoding[0].reconSpace.matrixSize @@ -540,7 +558,7 @@ def Pulseq_WaveformtoTrajectoryGadget(connection): # Fix GIRF correction but it is incorrect trajectory_GIRF, grad_GIRF = apply_GIRF_plus_cp(all_gradients, real_dwell_time, sR,0,batchsize=32*5) #RO INT DIM #tRR=1.25 GIRF_time_end=time() - print(' GIRF time %f s'%(GIRF_time_end-GIRF_time)) + eprint(' GIRF time %f s'%(GIRF_time_end-GIRF_time)) traj_c = trajectory_GIRF[:,:,:].real/scale_factor_before_GIRF/scale_factor_inside_GIRF for ax in range(len(max_xyz)): traj_c[...,ax] *= 0.5 / max_xyz[ax] 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/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();