From a26e6ca23053d521cbea00d8d1581be6ec2168f8 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Fri, 1 Nov 2024 12:22:52 +0100 Subject: [PATCH 01/18] Split trajectory.maths.py into a submodule --- src/mrinufft/trajectories/maths.py | 343 ------------------- src/mrinufft/trajectories/maths/__init__.py | 12 + src/mrinufft/trajectories/maths/constants.py | 7 + src/mrinufft/trajectories/maths/fibonacci.py | 135 ++++++++ src/mrinufft/trajectories/maths/primes.py | 32 ++ src/mrinufft/trajectories/maths/rotations.py | 164 +++++++++ 6 files changed, 350 insertions(+), 343 deletions(-) delete mode 100644 src/mrinufft/trajectories/maths.py create mode 100644 src/mrinufft/trajectories/maths/__init__.py create mode 100644 src/mrinufft/trajectories/maths/constants.py create mode 100644 src/mrinufft/trajectories/maths/fibonacci.py create mode 100644 src/mrinufft/trajectories/maths/primes.py create mode 100644 src/mrinufft/trajectories/maths/rotations.py diff --git a/src/mrinufft/trajectories/maths.py b/src/mrinufft/trajectories/maths.py deleted file mode 100644 index 906b1e38..00000000 --- a/src/mrinufft/trajectories/maths.py +++ /dev/null @@ -1,343 +0,0 @@ -"""Utility functions for mathematical operations.""" - -import numpy as np -import numpy.linalg as nl - -CIRCLE_PACKING_DENSITY = np.pi / (2 * np.sqrt(3)) -EIGENVECTOR_2D_FIBONACCI = (0.4656, 0.6823, 1) - - -########## -# PRIMES # -########## - - -def compute_coprime_factors(Nc, length, start=1, update=1): - """Compute a list of coprime factors of Nc. - - Parameters - ---------- - Nc : int - Number to factorize. - length : int - Number of coprime factors to compute. - start : int, optional - First number to check. The default is 1. - update : int, optional - Increment between two numbers to check. The default is 1. - - Returns - ------- - list - List of coprime factors of Nc. - """ - count = start - coprimes = [] - while len(coprimes) < length: - # Check greatest common divider (gcd) - if np.gcd(Nc, count) == 1: - coprimes.append(count) - count += update - return coprimes - - -############# -# ROTATIONS # -############# - - -def R2D(theta): - """Initialize 2D rotation matrix. - - Parameters - ---------- - theta : float - Rotation angle in rad. - - Returns - ------- - np.ndarray - 2D rotation matrix. - """ - return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) - - -def Rx(theta): - """Initialize 3D rotation matrix around x axis. - - Parameters - ---------- - theta : float - Rotation angle in rad. - - Returns - ------- - np.ndarray - 3D rotation matrix. - """ - return np.array( - [ - [1, 0, 0], - [0, np.cos(theta), -np.sin(theta)], - [0, np.sin(theta), np.cos(theta)], - ] - ) - - -def Ry(theta): - """Initialize 3D rotation matrix around y axis. - - Parameters - ---------- - theta : float - Rotation angle in rad. - - Returns - ------- - np.ndarray - 3D rotation matrix. - """ - return np.array( - [ - [np.cos(theta), 0, np.sin(theta)], - [0, 1, 0], - [-np.sin(theta), 0, np.cos(theta)], - ] - ) - - -def Rz(theta): - """Initialize 3D rotation matrix around z axis. - - Parameters - ---------- - theta : float - Rotation angle in rad. - - Returns - ------- - np.ndarray - 3D rotation matrix. - """ - return np.array( - [ - [np.cos(theta), -np.sin(theta), 0], - [np.sin(theta), np.cos(theta), 0], - [0, 0, 1], - ] - ) - - -def Rv(v1, v2, normalize=True, eps=1e-8): - """Initialize 3D rotation matrix from two vectors. - - Initialize a 3D rotation matrix from two vectors using Rodrigues's rotation - formula. Note that the rotation is carried around the axis orthogonal to both - vectors from the origin, and therefore is undetermined when both vectors - are colinear. While this case is handled manually, close cases might result - in approximative behaviors. - - Parameters - ---------- - v1 : np.ndarray - Source vector. - v2 : np.ndarray - Target vector. - normalize : bool, optional - Normalize the vectors. The default is True. - - Returns - ------- - np.ndarray - 3D rotation matrix. - """ - # Check for colinearity, not handled by Rodrigues' coefficients - if nl.norm(np.cross(v1, v2)) < eps: - sign = np.sign(np.dot(v1, v2)) - return sign * np.identity(3) - - if normalize: - v1, v2 = v1 / np.linalg.norm(v1), v2 / np.linalg.norm(v2) - cos_theta = np.dot(v1, v2) - v3 = np.cross(v1, v2) - cross_matrix = np.cross(v3, np.identity(v3.shape[0]) * -1) - return np.identity(3) + cross_matrix + cross_matrix @ cross_matrix / (1 + cos_theta) - - -def Ra(vector, theta): - """Initialize 3D rotation matrix around an arbitrary vector. - - Initialize a 3D rotation matrix to rotate around `vector` by an angle `theta`. - It corresponds to a generalized formula with `Rx`, `Ry` and `Rz` as subcases. - - Parameters - ---------- - vector : np.ndarray - Vector defining the rotation axis, automatically normalized. - theta : float - Angle in radians defining the rotation around `vector`. - - Returns - ------- - np.ndarray - 3D rotation matrix. - """ - cos_t = np.cos(theta) - sin_t = np.sin(theta) - v_x, v_y, v_z = vector / np.linalg.norm(vector) - return np.array( - [ - [ - cos_t + v_x**2 * (1 - cos_t), - v_x * v_y * (1 - cos_t) + v_z * sin_t, - v_x * v_z * (1 - cos_t) - v_y * sin_t, - ], - [ - v_y * v_x * (1 - cos_t) - v_z * sin_t, - cos_t + v_y**2 * (1 - cos_t), - v_y * v_z * (1 - cos_t) + v_x * sin_t, - ], - [ - v_z * v_x * (1 - cos_t) + v_y * sin_t, - v_z * v_y * (1 - cos_t) - v_x * sin_t, - cos_t + v_z**2 * (1 - cos_t), - ], - ] - ) - - -############# -# FIBONACCI # -############# - - -def is_from_fibonacci_sequence(n): - """Check if an integer belongs to the Fibonacci sequence. - - An integer belongs to the Fibonacci sequence if either - :math:`5*n²+4` or :math:`5*n²-4` is a perfect square - (`Wikipedia `_). - - Parameters - ---------- - n : int - Integer to check. - - Returns - ------- - bool - Whether or not ``n`` belongs to the Fibonacci sequence. - """ - - def _is_perfect_square(n): - r = int(np.sqrt(n)) - return r * r == n - - return _is_perfect_square(5 * n**2 + 4) or _is_perfect_square(5 * n**2 - 4) - - -def get_closest_fibonacci_number(x): - """Provide the closest Fibonacci number. - - Parameters - ---------- - x : float - Value to match. - - Returns - ------- - int - Closest number from the Fibonacci sequence. - """ - # Find the power such that x ~= phi ** power - phi = (1 + np.sqrt(5)) / 2 - power = np.ceil(np.log(x) / np.log(phi)) + 1 - - # Check closest between the ones below and above n - lower_xf = int(np.around(phi ** (power) / np.sqrt(5))) - upper_xf = int(np.around(phi ** (power + 1) / np.sqrt(5))) - xf = lower_xf if (x - lower_xf) < (upper_xf - x) else upper_xf - return xf - - -def generate_fibonacci_lattice(nb_points, epsilon=0.25): - """Generate 2D Cartesian coordinates using the Fibonacci lattice. - - Place 2D points over a 1x1 square following the Fibonacci lattice. - - Parameters - ---------- - nb_points : int - Number of 2D points to generate. - epsilon : float - Continuous offset used to reduce initially wrong lattice behavior. - - Returns - ------- - np.ndarray - Array of 2D Cartesian coordinates covering a 1x1 square. - """ - angle = (1 + np.sqrt(5)) / 2 - fibonacci_square = np.zeros((nb_points, 2)) - fibonacci_square[:, 0] = (np.arange(nb_points) / angle) % 1 - fibonacci_square[:, 1] = (np.arange(nb_points) + epsilon) / ( - nb_points - 1 + 2 * epsilon - ) - return fibonacci_square - - -def generate_fibonacci_circle(nb_points, epsilon=0.25): - """Generate 2D Cartesian coordinates shaped as Fibonacci spirals. - - Place 2D points structured as Fibonacci spirals by distorting - a square Fibonacci lattice into a circle of radius 1. - - Parameters - ---------- - nb_points : int - Number of 2D points to generate. - epsilon : float - Continuous offset used to reduce initially wrong lattice behavior. - - Returns - ------- - np.ndarray - Array of 2D Cartesian coordinates covering a circle of radius 1. - """ - fibonacci_square = generate_fibonacci_lattice(nb_points, epsilon) - radius = np.sqrt(fibonacci_square[:, 1]) - angles = 2 * np.pi * fibonacci_square[:, 0] - - fibonacci_circle = np.zeros((nb_points, 2)) - fibonacci_circle[:, 0] = radius * np.cos(angles) - fibonacci_circle[:, 1] = radius * np.sin(angles) - return fibonacci_circle - - -def generate_fibonacci_sphere(nb_points, epsilon=0.25): - """Generate 3D Cartesian coordinates as a Fibonacci sphere. - - Place 3D points almost evenly over a sphere surface of radius - 1 by distorting a square Fibonacci lattice into a sphere. - - Parameters - ---------- - nb_points : int - Number of 3D points to generate. - epsilon : float - Continuous offset used to reduce initially wrong lattice behavior. - - Returns - ------- - np.ndarray - Array of 3D Cartesian coordinates covering a sphere of radius 1. - """ - fibonacci_square = generate_fibonacci_lattice(nb_points, epsilon) - theta = 2 * np.pi * fibonacci_square[:, 0] - phi = np.arccos(1 - 2 * fibonacci_square[:, 1]) - - fibonacci_sphere = np.zeros((nb_points, 3)) - fibonacci_sphere[:, 0] = np.cos(theta) * np.sin(phi) - fibonacci_sphere[:, 1] = np.sin(theta) * np.sin(phi) - fibonacci_sphere[:, 2] = np.cos(phi) - return fibonacci_sphere diff --git a/src/mrinufft/trajectories/maths/__init__.py b/src/mrinufft/trajectories/maths/__init__.py new file mode 100644 index 00000000..2f721cfc --- /dev/null +++ b/src/mrinufft/trajectories/maths/__init__.py @@ -0,0 +1,12 @@ +"""Utility module for mathematical operations.""" + +from .constants import CIRCLE_PACKING_DENSITY, EIGENVECTOR_2D_FIBONACCI +from .fibonacci import ( + is_from_fibonacci_sequence, + get_closest_fibonacci_number, + generate_fibonacci_lattice, + generate_fibonacci_circle, + generate_fibonacci_sphere, +) +from .primes import compute_coprime_factors +from .rotations import R2D, Rx, Ry, Rz, Rv, Ra diff --git a/src/mrinufft/trajectories/maths/constants.py b/src/mrinufft/trajectories/maths/constants.py new file mode 100644 index 00000000..adbd338c --- /dev/null +++ b/src/mrinufft/trajectories/maths/constants.py @@ -0,0 +1,7 @@ +"""Constant values for mathematical purposes.""" + +import numpy as np + + +CIRCLE_PACKING_DENSITY = np.pi / (2 * np.sqrt(3)) +EIGENVECTOR_2D_FIBONACCI = (0.4656, 0.6823, 1) diff --git a/src/mrinufft/trajectories/maths/fibonacci.py b/src/mrinufft/trajectories/maths/fibonacci.py new file mode 100644 index 00000000..d7cc20a8 --- /dev/null +++ b/src/mrinufft/trajectories/maths/fibonacci.py @@ -0,0 +1,135 @@ +"""Fibonacci-related functions.""" + +import numpy as np + + +def is_from_fibonacci_sequence(n): + """Check if an integer belongs to the Fibonacci sequence. + + An integer belongs to the Fibonacci sequence if either + :math:`5*n²+4` or :math:`5*n²-4` is a perfect square + (`Wikipedia `_). + + Parameters + ---------- + n : int + Integer to check. + + Returns + ------- + bool + Whether or not ``n`` belongs to the Fibonacci sequence. + """ + + def _is_perfect_square(n): + r = int(np.sqrt(n)) + return r * r == n + + return _is_perfect_square(5 * n**2 + 4) or _is_perfect_square(5 * n**2 - 4) + + +def get_closest_fibonacci_number(x): + """Provide the closest Fibonacci number. + + Parameters + ---------- + x : float + Value to match. + + Returns + ------- + int + Closest number from the Fibonacci sequence. + """ + # Find the power such that x ~= phi ** power + phi = (1 + np.sqrt(5)) / 2 + power = np.ceil(np.log(x) / np.log(phi)) + 1 + + # Check closest between the ones below and above n + lower_xf = int(np.around(phi ** (power) / np.sqrt(5))) + upper_xf = int(np.around(phi ** (power + 1) / np.sqrt(5))) + xf = lower_xf if (x - lower_xf) < (upper_xf - x) else upper_xf + return xf + + +def generate_fibonacci_lattice(nb_points, epsilon=0.25): + """Generate 2D Cartesian coordinates using the Fibonacci lattice. + + Place 2D points over a 1x1 square following the Fibonacci lattice. + + Parameters + ---------- + nb_points : int + Number of 2D points to generate. + epsilon : float + Continuous offset used to reduce initially wrong lattice behavior. + + Returns + ------- + np.ndarray + Array of 2D Cartesian coordinates covering a 1x1 square. + """ + angle = (1 + np.sqrt(5)) / 2 + fibonacci_square = np.zeros((nb_points, 2)) + fibonacci_square[:, 0] = (np.arange(nb_points) / angle) % 1 + fibonacci_square[:, 1] = (np.arange(nb_points) + epsilon) / ( + nb_points - 1 + 2 * epsilon + ) + return fibonacci_square + + +def generate_fibonacci_circle(nb_points, epsilon=0.25): + """Generate 2D Cartesian coordinates shaped as Fibonacci spirals. + + Place 2D points structured as Fibonacci spirals by distorting + a square Fibonacci lattice into a circle of radius 1. + + Parameters + ---------- + nb_points : int + Number of 2D points to generate. + epsilon : float + Continuous offset used to reduce initially wrong lattice behavior. + + Returns + ------- + np.ndarray + Array of 2D Cartesian coordinates covering a circle of radius 1. + """ + fibonacci_square = generate_fibonacci_lattice(nb_points, epsilon) + radius = np.sqrt(fibonacci_square[:, 1]) + angles = 2 * np.pi * fibonacci_square[:, 0] + + fibonacci_circle = np.zeros((nb_points, 2)) + fibonacci_circle[:, 0] = radius * np.cos(angles) + fibonacci_circle[:, 1] = radius * np.sin(angles) + return fibonacci_circle + + +def generate_fibonacci_sphere(nb_points, epsilon=0.25): + """Generate 3D Cartesian coordinates as a Fibonacci sphere. + + Place 3D points almost evenly over a sphere surface of radius + 1 by distorting a square Fibonacci lattice into a sphere. + + Parameters + ---------- + nb_points : int + Number of 3D points to generate. + epsilon : float + Continuous offset used to reduce initially wrong lattice behavior. + + Returns + ------- + np.ndarray + Array of 3D Cartesian coordinates covering a sphere of radius 1. + """ + fibonacci_square = generate_fibonacci_lattice(nb_points, epsilon) + theta = 2 * np.pi * fibonacci_square[:, 0] + phi = np.arccos(1 - 2 * fibonacci_square[:, 1]) + + fibonacci_sphere = np.zeros((nb_points, 3)) + fibonacci_sphere[:, 0] = np.cos(theta) * np.sin(phi) + fibonacci_sphere[:, 1] = np.sin(theta) * np.sin(phi) + fibonacci_sphere[:, 2] = np.cos(phi) + return fibonacci_sphere diff --git a/src/mrinufft/trajectories/maths/primes.py b/src/mrinufft/trajectories/maths/primes.py new file mode 100644 index 00000000..3a9df50c --- /dev/null +++ b/src/mrinufft/trajectories/maths/primes.py @@ -0,0 +1,32 @@ +"""Prime-related functions.""" + +import numpy as np + + +def compute_coprime_factors(Nc, length, start=1, update=1): + """Compute a list of coprime factors of Nc. + + Parameters + ---------- + Nc : int + Number to factorize. + length : int + Number of coprime factors to compute. + start : int, optional + First number to check. The default is 1. + update : int, optional + Increment between two numbers to check. The default is 1. + + Returns + ------- + list + List of coprime factors of Nc. + """ + count = start + coprimes = [] + while len(coprimes) < length: + # Check greatest common divider (gcd) + if np.gcd(Nc, count) == 1: + coprimes.append(count) + count += update + return coprimes diff --git a/src/mrinufft/trajectories/maths/rotations.py b/src/mrinufft/trajectories/maths/rotations.py new file mode 100644 index 00000000..56878078 --- /dev/null +++ b/src/mrinufft/trajectories/maths/rotations.py @@ -0,0 +1,164 @@ +"""Rotation functions in 2D & 3D spaces.""" + +import numpy as np +import numpy.linalg as nl + + +def R2D(theta): + """Initialize 2D rotation matrix. + + Parameters + ---------- + theta : float + Rotation angle in rad. + + Returns + ------- + np.ndarray + 2D rotation matrix. + """ + return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) + + +def Rx(theta): + """Initialize 3D rotation matrix around x axis. + + Parameters + ---------- + theta : float + Rotation angle in rad. + + Returns + ------- + np.ndarray + 3D rotation matrix. + """ + return np.array( + [ + [1, 0, 0], + [0, np.cos(theta), -np.sin(theta)], + [0, np.sin(theta), np.cos(theta)], + ] + ) + + +def Ry(theta): + """Initialize 3D rotation matrix around y axis. + + Parameters + ---------- + theta : float + Rotation angle in rad. + + Returns + ------- + np.ndarray + 3D rotation matrix. + """ + return np.array( + [ + [np.cos(theta), 0, np.sin(theta)], + [0, 1, 0], + [-np.sin(theta), 0, np.cos(theta)], + ] + ) + + +def Rz(theta): + """Initialize 3D rotation matrix around z axis. + + Parameters + ---------- + theta : float + Rotation angle in rad. + + Returns + ------- + np.ndarray + 3D rotation matrix. + """ + return np.array( + [ + [np.cos(theta), -np.sin(theta), 0], + [np.sin(theta), np.cos(theta), 0], + [0, 0, 1], + ] + ) + + +def Rv(v1, v2, normalize=True, eps=1e-8): + """Initialize 3D rotation matrix from two vectors. + + Initialize a 3D rotation matrix from two vectors using Rodrigues's rotation + formula. Note that the rotation is carried around the axis orthogonal to both + vectors from the origin, and therefore is undetermined when both vectors + are colinear. While this case is handled manually, close cases might result + in approximative behaviors. + + Parameters + ---------- + v1 : np.ndarray + Source vector. + v2 : np.ndarray + Target vector. + normalize : bool, optional + Normalize the vectors. The default is True. + + Returns + ------- + np.ndarray + 3D rotation matrix. + """ + # Check for colinearity, not handled by Rodrigues' coefficients + if nl.norm(np.cross(v1, v2)) < eps: + sign = np.sign(np.dot(v1, v2)) + return sign * np.identity(3) + + if normalize: + v1, v2 = v1 / np.linalg.norm(v1), v2 / np.linalg.norm(v2) + cos_theta = np.dot(v1, v2) + v3 = np.cross(v1, v2) + cross_matrix = np.cross(v3, np.identity(v3.shape[0]) * -1) + return np.identity(3) + cross_matrix + cross_matrix @ cross_matrix / (1 + cos_theta) + + +def Ra(vector, theta): + """Initialize 3D rotation matrix around an arbitrary vector. + + Initialize a 3D rotation matrix to rotate around `vector` by an angle `theta`. + It corresponds to a generalized formula with `Rx`, `Ry` and `Rz` as subcases. + + Parameters + ---------- + vector : np.ndarray + Vector defining the rotation axis, automatically normalized. + theta : float + Angle in radians defining the rotation around `vector`. + + Returns + ------- + np.ndarray + 3D rotation matrix. + """ + cos_t = np.cos(theta) + sin_t = np.sin(theta) + v_x, v_y, v_z = vector / np.linalg.norm(vector) + return np.array( + [ + [ + cos_t + v_x**2 * (1 - cos_t), + v_x * v_y * (1 - cos_t) + v_z * sin_t, + v_x * v_z * (1 - cos_t) - v_y * sin_t, + ], + [ + v_y * v_x * (1 - cos_t) - v_z * sin_t, + cos_t + v_y**2 * (1 - cos_t), + v_y * v_z * (1 - cos_t) + v_x * sin_t, + ], + [ + v_z * v_x * (1 - cos_t) + v_y * sin_t, + v_z * v_y * (1 - cos_t) - v_x * sin_t, + cos_t + v_z**2 * (1 - cos_t), + ], + ] + ) From 0f12b7e45a0b423a28374af26b4aa4f27bc359ba Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Fri, 1 Nov 2024 19:28:51 +0100 Subject: [PATCH 02/18] Add densities, random walk and travelling salesman trajectories --- src/mrinufft/trajectories/__init__.py | 85 ++++---- src/mrinufft/trajectories/densities.py | 132 ++++++++++++ src/mrinufft/trajectories/inits/__init__.py | 14 ++ .../trajectories/inits/random_walk.py | 95 +++++++++ .../trajectories/inits/travelling_salesman.py | 196 ++++++++++++++++++ src/mrinufft/trajectories/maths/__init__.py | 9 +- src/mrinufft/trajectories/maths/constants.py | 1 - src/mrinufft/trajectories/maths/tsp_solver.py | 33 +++ src/mrinufft/trajectories/trajectory3D.py | 2 +- 9 files changed, 526 insertions(+), 41 deletions(-) create mode 100644 src/mrinufft/trajectories/densities.py create mode 100644 src/mrinufft/trajectories/inits/__init__.py create mode 100644 src/mrinufft/trajectories/inits/random_walk.py create mode 100644 src/mrinufft/trajectories/inits/travelling_salesman.py create mode 100644 src/mrinufft/trajectories/maths/tsp_solver.py diff --git a/src/mrinufft/trajectories/__init__.py b/src/mrinufft/trajectories/__init__.py index bed80d3b..1e823c03 100644 --- a/src/mrinufft/trajectories/__init__.py +++ b/src/mrinufft/trajectories/__init__.py @@ -1,55 +1,62 @@ """Collection of trajectories and tools used for non-Cartesian MRI.""" +from .densities import ( + create_chauffert_density, + create_cutoff_decay_density, + create_energy_density, + create_fast_chauffert_density, + create_polynomial_density, + sample_from_density, +) +from .display import display_2D_trajectory, display_3D_trajectory, displayConfig +from .gradients import patch_center_anomaly +from .inits import ( + initialize_2D_random_walk, + initialize_2D_travelling_salesman, + initialize_3D_random_walk, + initialize_3D_travelling_salesman, +) +from .tools import ( + conify, + duplicate_along_axes, + precess, + radialize_center, + rotate, + shellify, + stack, + stack_spherically, +) from .trajectory2D import ( - initialize_2D_radial, - initialize_2D_spiral, - initialize_2D_fibonacci_spiral, initialize_2D_cones, - initialize_2D_sinusoide, + initialize_2D_fibonacci_spiral, + initialize_2D_lissajous, + initialize_2D_polar_lissajous, initialize_2D_propeller, - initialize_2D_rosette, + initialize_2D_radial, initialize_2D_rings, - initialize_2D_polar_lissajous, - initialize_2D_lissajous, + initialize_2D_rosette, + initialize_2D_sinusoide, + initialize_2D_spiral, initialize_2D_waves, ) - from .trajectory3D import ( - initialize_3D_phyllotaxis_radial, - initialize_3D_golden_means_radial, - initialize_3D_wong_radial, - initialize_3D_park_radial, + initialize_3D_annular_shells, initialize_3D_cones, initialize_3D_floret, - initialize_3D_wave_caipi, - initialize_3D_seiffert_spiral, + initialize_3D_golden_means_radial, initialize_3D_helical_shells, - initialize_3D_annular_shells, + initialize_3D_park_radial, + initialize_3D_phyllotaxis_radial, + initialize_3D_repi, initialize_3D_seiffert_shells, + initialize_3D_seiffert_spiral, initialize_3D_turbine, - initialize_3D_repi, -) - -from .tools import ( - stack, - rotate, - precess, - conify, - stack_spherically, - shellify, - duplicate_along_axes, - radialize_center, -) - -from .display import ( - displayConfig, - display_2D_trajectory, - display_3D_trajectory, + initialize_3D_wave_caipi, + initialize_3D_wong_radial, ) -from .gradients import patch_center_anomaly - __all__ = [ + # Trajectories "initialize_2D_radial", "initialize_2D_spiral", "initialize_2D_fibonacci_spiral", @@ -75,6 +82,7 @@ "initialize_3D_seiffert_shells", "initialize_3D_turbine", "initialize_3D_repi", + # Tools "stack", "rotate", "precess", @@ -87,4 +95,11 @@ "display_2D_trajectory", "display_3D_trajectory", "patch_center_anomaly", + # Densities + "sample_from_density", + "create_cutoff_decay_density", + "create_polynomial_density", + "create_energy_density", + "create_chauffert_density", + "create_fast_chauffert_density", ] diff --git a/src/mrinufft/trajectories/densities.py b/src/mrinufft/trajectories/densities.py new file mode 100644 index 00000000..6c0cb8ae --- /dev/null +++ b/src/mrinufft/trajectories/densities.py @@ -0,0 +1,132 @@ +import numpy as np +import numpy.fft as nf +import numpy.linalg as nl +import numpy.random as nr +import pywt as pw +from sklearn.cluster import BisectingKMeans, KMeans +from tqdm.auto import tqdm + +from .utils import KMAX + + +def sample_from_density(nb_samples, density, method="random"): + rng = nr.default_rng() + + shape = density.shape + nb_dims = len(shape) + max_nb_samples = np.prod(shape) + + if nb_samples > max_nb_samples: + raise ValueError("`nb_samples` must be lower than the size of `density`.") + + if method == "random": + choices = rng.choice( + np.arange(max_nb_samples), + size=nb_samples, + p=density.flatten(), + replace=False, + ) + locations = np.indices(shape).reshape((nb_dims, -1))[:, choices] + locations = locations.T + locations = 2 * KMAX * locations / np.max(locations) - KMAX + elif method == "lloyd": + kmeans = ( + KMeans(n_clusters=nb_samples) + if nb_dims <= 2 + else BisectingKMeans(n_clusters=nb_samples) + ) + kmeans.fit( + np.indices(density.shape).reshape((nb_dims, -1)).T, + sample_weight=density.flatten(), + ) + locations = kmeans.cluster_centers_ - np.array(density.shape) / 2 + locations = KMAX * locations / np.max(np.abs(locations)) + else: + raise ValueError(f"Unknown sampling method {method}.") + return locations + + +def create_cutoff_decay_density(shape, cutoff, decay, resolution=None): + shape = np.array(shape) + nb_dims = len(shape) + + if not resolution: + resolution = np.ones(nb_dims) + + distances = np.indices(shape).astype(float) + for i in range(nb_dims): + distances[i] = distances[i] + 0.5 - shape[i] / 2 + distances[i] = distances[i] / shape[i] * resolution[i] + distances = nl.norm(distances, axis=0) + + cutoff = cutoff * np.max(distances) if cutoff else np.min(distances) + density = np.ones(shape) + decay_mask = np.where(distances > cutoff, True, False) + + density[decay_mask] = (cutoff / distances[decay_mask]) ** decay + density = density / np.sum(density) + return density + + +def create_polynomial_density(shape, decay, resolution=None): + return create_cutoff_decay_density( + shape, cutoff=0, decay=decay, resolution=resolution + ) + + +def create_energy_density(dataset): + nb_dims = len(dataset.shape) - 1 + axes = range(1, nb_dims + 1) + + kspace = nf.fftshift(nf.fftn(nf.fftshift(dataset, axes=axes), axes=axes), axes=axes) + energy = np.mean(np.abs(kspace), axis=0) + density = energy / np.sum(energy) + return density + + +def create_chauffert_density(shape, wavelet_basis, nb_wavelet_scales, verbose=False): + nb_dims = len(shape) + indices = np.indices(shape).reshape((nb_dims, -1)).T + + density = np.zeros(shape) + unit_vector = np.zeros(shape) + for ids in tqdm(indices, disable=not verbose): + ids = tuple(ids) + unit_vector[ids] = 1 + fourier_vector = nf.ifftn(unit_vector) + coeffs = pw.wavedecn( + fourier_vector, wavelet=wavelet_basis, level=nb_wavelet_scales + ) + coeffs, _ = pw.coeffs_to_array(coeffs) + density[ids] = np.max(np.abs(coeffs)) ** 2 + unit_vector[ids] = 0 + + density = density / np.sum(density) + return nf.ifftshift(density) + + +def create_fast_chauffert_density(shape, wavelet_basis, nb_wavelet_scales): + nb_dims = len(shape) + + density = np.ones(shape) + for k, s in enumerate(shape): + unit_vector = np.zeros(s) + density_1d = np.zeros(s) + + for i in range(s): + unit_vector[i] = 1 + fourier_vector = nf.ifft(unit_vector) + coeffs = pw.wavedec( + fourier_vector, wavelet=wavelet_basis, level=nb_wavelet_scales + ) + coeffs, _ = pw.coeffs_to_array(coeffs) + density_1d[i] = np.max(np.abs(coeffs)) ** 2 + unit_vector[i] = 0 + + reshape = np.ones(nb_dims).astype(int) + reshape[k] = s + density_1d = density_1d.reshape(reshape) + density = density * density_1d + + density = density / np.sum(density) + return nf.ifftshift(density) diff --git a/src/mrinufft/trajectories/inits/__init__.py b/src/mrinufft/trajectories/inits/__init__.py new file mode 100644 index 00000000..31c9da93 --- /dev/null +++ b/src/mrinufft/trajectories/inits/__init__.py @@ -0,0 +1,14 @@ +"""Module containing all trajectories.""" + +from .random_walk import initialize_2D_random_walk, initialize_3D_random_walk +from .travelling_salesman import ( + initialize_2D_travelling_salesman, + initialize_3D_travelling_salesman, +) + +__all__ = [ + "initialize_2D_random_walk", + "initialize_3D_random_walk", + "initialize_2D_travelling_salesman", + "initialize_3D_travelling_salesman", +] diff --git a/src/mrinufft/trajectories/inits/random_walk.py b/src/mrinufft/trajectories/inits/random_walk.py new file mode 100644 index 00000000..01ab9962 --- /dev/null +++ b/src/mrinufft/trajectories/inits/random_walk.py @@ -0,0 +1,95 @@ +"""Trajectories based on random walks.""" + +import numpy as np + +from ..utils import KMAX + + +def _get_adjacent_neighbors_flat_offsets(shape): + nb_dims = len(shape) + neighborhood = np.indices([3] * nb_dims) - 1 + distances = np.sum(np.abs(neighborhood), axis=0) + + center = np.ravel_multi_index([1] * nb_dims, dims=shape) + neighbors = np.ravel_multi_index(np.where(distances == 1), dims=shape) - center + return neighbors + + +def _get_diagonal_neighbors_flat_offsets(shape): + nb_dims = len(shape) + neighborhood = np.indices([3] * nb_dims) - 1 + distances = np.sum(np.abs(neighborhood), axis=0) + + center = np.ravel_multi_index([1] * nb_dims, dims=shape) + neighbors = np.ravel_multi_index(np.where(distances > 1), dims=shape) - center + return neighbors + + +def _initialize_ND_random_walk(Nc, Ns, density, diagonals=True, pseudo_random=True): + flat_density = np.copy(density.flatten()) + max_id = np.prod(density.shape) + mask = np.ones_like(flat_density) + + # Prepare neighbor offsets once + offsets = _get_adjacent_neighbors_flat_offsets(density.shape) + if diagonals: + offsets = np.concatenate( + [offsets, _get_diagonal_neighbors_flat_offsets(density.shape)] + ) + + # Make all random draws at once for performance + draws = np.random.random((Ns, Nc)) # inverted shape for convenience + + # Initialize shot starting points + choices = np.random.choice(np.arange(len(flat_density)), size=Nc, p=flat_density) + routes = [choices] + + # Walk + for i in range(1, Ns): + neighbors = choices[:, None] + offsets[None] + + # Find out-of-bound neighbors and ignore them + invalids = (neighbors < 0) | (neighbors > max_id) + neighbors[invalids] = 0 + + # Set walk probabilities + walk_probs = flat_density[neighbors] + walk_probs[invalids] = 0 + walk_probs = walk_probs / np.sum(walk_probs, axis=-1, keepdims=True) + cum_walk_probs = np.cumsum(walk_probs, axis=-1) + + # Select next walk steps + indices = np.argmax(draws[i][:, None] < cum_walk_probs, axis=-1) + choices = neighbors[np.arange(Nc), indices] + routes.append(choices) + + # Update density to account for already drawed positions + if pseudo_random: + flat_density[choices] = ( + mask[choices] * flat_density[choices] / (mask[choices] + 1) + ) + mask[choices] += 1 + routes = np.array(routes).T + + # Create trajectory from routes + locations = np.indices(density.shape) + locations = locations.reshape((len(density.shape), -1)) + trajectory = np.array([locations[:, r].T for r in routes]) + trajectory = 2 * KMAX * trajectory / density.shape - KMAX + return trajectory + + +def initialize_2D_random_walk(Nc, Ns, density, diagonals=True, pseudo_random=True): + if len(density.shape) != 2: + raise ValueError("`density` is expected to be 2-dimensional.") + return _initialize_ND_random_walk( + Nc, Ns, density, diagonals=diagonals, pseudo_random=pseudo_random + ) + + +def initialize_3D_random_walk(Nc, Ns, density, diagonals=True, pseudo_random=True): + if len(density.shape) != 3: + raise ValueError("`density` is expected to be 3-dimensional.") + return _initialize_ND_random_walk( + Nc, Ns, density, diagonals=diagonals, pseudo_random=pseudo_random + ) diff --git a/src/mrinufft/trajectories/inits/travelling_salesman.py b/src/mrinufft/trajectories/inits/travelling_salesman.py new file mode 100644 index 00000000..0cb69ba2 --- /dev/null +++ b/src/mrinufft/trajectories/inits/travelling_salesman.py @@ -0,0 +1,196 @@ +"""Trajectories based on the Travelig Salesman Problem.""" + +import numpy as np +import numpy.linalg as nl +from scipy.interpolate import CubicSpline +from tqdm.auto import tqdm + +from ..densities import sample_from_density +from ..maths import solve_tsp_with_2opt + + +def _get_approx_cluster_sizes(nb_total, nb_clusters): + cluster_sizes = round(nb_total / nb_clusters) * np.ones(nb_clusters).astype(int) + delta_sum = nb_total - np.sum(cluster_sizes) + cluster_sizes[: int(np.abs(delta_sum))] += np.sign(delta_sum) + return cluster_sizes + + +def _sort_by_coordinate(array, coord): + if array.shape[-1] < 3 and coord.lower() in ["z", "theta"]: + raise ValueError( + f"Invalid `coord`='{coord}' for arrays with less than 3 dimensions." + ) + + match coord.lower(): + case "x": + coord = array[..., 0] + case "y": + coord = array[..., 1] + case "z": + coord = array[..., 2] + case "r": + coord = np.linalg.norm(array, axis=-1) + case "phi": + coord = np.sign(array[..., 1]) * np.arccos( + array[..., 0] / nl.norm(array[..., :2], axis=-1) + ) + case "theta": + coord = np.arccos(array[..., 2] / nl.norm(array, axis=-1)) + case _: + raise ValueError(f"Unknown coordinate `{coord}`") + order = np.argsort(coord) + return array[order] + + +def _cluster_by_coordinate( + locations, nb_clusters, cluster_by, second_cluster_by=None, sort_by=None +): + # Gather dimension variables + nb_dims = locations.shape[-1] + locations = locations.reshape((-1, nb_dims)) + nb_locations = locations.shape[0] + + # Check arguments validity + if nb_locations % nb_clusters: + raise ValueError("`nb_clusters` should divide the number of locations") + cluster_size = nb_locations // nb_clusters + + # Create chunks of cluters by a first coordinate + locations = _sort_by_coordinate(locations, cluster_by) + + if second_cluster_by: + # Cluster each location within the chunks of clusters by a second coordinate + chunk_sizes = _get_approx_cluster_sizes( + nb_clusters, round(np.sqrt(nb_clusters)) + ) + chunk_ranges = np.cumsum([0] + list(chunk_sizes)) + for i in range(len(chunk_sizes)): + i_s, i_e = ( + chunk_ranges[i] * cluster_size, + chunk_ranges[i + 1] * cluster_size, + ) + locations[i_s:i_e] = _sort_by_coordinate( + locations[i_s:i_e], second_cluster_by + ) + locations = locations.reshape((nb_clusters, cluster_size, nb_dims)) + + # Order locations within each cluster by another coordinate + if sort_by: + for i in range(nb_clusters): + locations[i] = _sort_by_coordinate(locations[i], sort_by) + return locations + + +def _initialize_ND_travelling_salesman( + Nc, + Ns, + density, + first_cluster_by=None, + second_cluster_by=None, + sort_by=None, + nb_tsp_points="auto", + sampling="random", + tsp_tol=1e-8, + verbose=False, +): + # Handle variable inputs + nb_tsp_points = Ns if nb_tsp_points == "auto" else nb_tsp_points + + # Check arguments validity + if Nc * nb_tsp_points > np.prod(density.shape): + raise ValueError( + "`density` array not large enough to peak `Nc` * `nb_tsp_points` points." + ) + Nd = len(density.shape) + + # Select k-space locations + density = density / np.sum(density) + locations = sample_from_density(Nc * Ns, density, method=sampling) + + # Re-organise locations into Nc clusters + if first_cluster_by: + locations = _cluster_by_coordinate( + locations, + Nc, + cluster_by=first_cluster_by, + second_cluster_by=second_cluster_by, + sort_by=sort_by, + ) + + # Compute TSP solution within each cluster/shot + for i in tqdm(range(Nc), disable=not verbose): + order = solve_tsp_with_2opt(locations[i], improvement_threshold=tsp_tol) + locations[i] = locations[i][order] + else: + locations = ( + _sort_by_coordinate(locations, coord=sort_by) if sort_by else locations + ) + + # Compute TSP solution over the whole cloud + order = solve_tsp_with_2opt(locations, improvement_threshold=tsp_tol) + locations = locations[order] + locations = locations.reshape((Nc, Ns, Nd)) + + # Interpolate shot points up to full length + trajectory = np.zeros((Nc, Ns, Nd)) + for i in range(Nc): + cbs = CubicSpline(np.linspace(0, 1, nb_tsp_points), locations[i]) + trajectory[i] = cbs(np.linspace(0, 1, Ns)) + return trajectory + + +def initialize_2D_travelling_salesman( + Nc, + Ns, + density, + first_cluster_by=None, + second_cluster_by=None, + sort_by=None, + nb_tsp_points="auto", + sampling="random", + tsp_tol=1e-8, + verbose=False, +): + if len(density.shape) != 2: + raise ValueError("`density` is expected to be 2-dimensional.") + return _initialize_ND_travelling_salesman( + Nc, + Ns, + density, + first_cluster_by=first_cluster_by, + second_cluster_by=second_cluster_by, + sort_by=sort_by, + nb_tsp_points=nb_tsp_points, + sampling=sampling, + tsp_tol=tsp_tol, + verbose=verbose, + ) + + +def initialize_3D_travelling_salesman( + Nc, + Ns, + density, + first_cluster_by=None, + second_cluster_by=None, + sort_by=None, + nb_tsp_points="auto", + sampling="random", + tsp_tol=1e-8, + verbose=False, +): + if len(density.shape) != 3: + raise ValueError("`density` is expected to be 3-dimensional.") + return _initialize_ND_travelling_salesman( + Nc, + Ns, + density, + first_cluster_by=first_cluster_by, + second_cluster_by=second_cluster_by, + sort_by=sort_by, + nb_tsp_points=nb_tsp_points, + sampling=sampling, + tsp_tol=tsp_tol, + verbose=verbose, + ) diff --git a/src/mrinufft/trajectories/maths/__init__.py b/src/mrinufft/trajectories/maths/__init__.py index 2f721cfc..36537a5e 100644 --- a/src/mrinufft/trajectories/maths/__init__.py +++ b/src/mrinufft/trajectories/maths/__init__.py @@ -2,11 +2,12 @@ from .constants import CIRCLE_PACKING_DENSITY, EIGENVECTOR_2D_FIBONACCI from .fibonacci import ( - is_from_fibonacci_sequence, - get_closest_fibonacci_number, - generate_fibonacci_lattice, generate_fibonacci_circle, + generate_fibonacci_lattice, generate_fibonacci_sphere, + get_closest_fibonacci_number, + is_from_fibonacci_sequence, ) from .primes import compute_coprime_factors -from .rotations import R2D, Rx, Ry, Rz, Rv, Ra +from .rotations import R2D, Ra, Rv, Rx, Ry, Rz +from .tsp_solver import solve_tsp_with_2opt diff --git a/src/mrinufft/trajectories/maths/constants.py b/src/mrinufft/trajectories/maths/constants.py index adbd338c..179dbf05 100644 --- a/src/mrinufft/trajectories/maths/constants.py +++ b/src/mrinufft/trajectories/maths/constants.py @@ -2,6 +2,5 @@ import numpy as np - CIRCLE_PACKING_DENSITY = np.pi / (2 * np.sqrt(3)) EIGENVECTOR_2D_FIBONACCI = (0.4656, 0.6823, 1) diff --git a/src/mrinufft/trajectories/maths/tsp_solver.py b/src/mrinufft/trajectories/maths/tsp_solver.py new file mode 100644 index 00000000..092b8b7a --- /dev/null +++ b/src/mrinufft/trajectories/maths/tsp_solver.py @@ -0,0 +1,33 @@ +"""Solver for the Travelling Salesman Problem.""" + +import numpy as np + + +def solve_tsp_with_2opt(locations, improvement_threshold=1e-8): + route = np.arange(locations.shape[0]) + distances = np.linalg.norm(locations[None] - locations[:, None], axis=-1) + route_length = np.sum(distances[0]) + + improvement_factor = 1 + while improvement_factor >= improvement_threshold: + old_route_length = route_length + for i in range(1, len(route) - 2): + # Check new distance by reversing chunks between i and j + for j in range(i + 1, len(route) - 1): + # Compute new route distance variation + delta_length = ( + distances[route[i - 1], route[j]] + + distances[route[i], route[j + 1]] + - distances[route[i - 1], route[i]] + - distances[route[j], route[j + 1]] + ) + + if delta_length < 0: + # Reverse route chunk + route = np.concatenate( + [route[:i], route[i : j + 1][::-1], route[j + 1 :]] + ) + route_length += delta_length + + improvement_factor = abs(1 - route_length / old_route_length) + return route diff --git a/src/mrinufft/trajectories/trajectory3D.py b/src/mrinufft/trajectories/trajectory3D.py index 12dd9241..cdcf4a0a 100644 --- a/src/mrinufft/trajectories/trajectory3D.py +++ b/src/mrinufft/trajectories/trajectory3D.py @@ -8,12 +8,12 @@ from .maths import ( CIRCLE_PACKING_DENSITY, + EIGENVECTOR_2D_FIBONACCI, R2D, Ra, Ry, Rz, generate_fibonacci_circle, - EIGENVECTOR_2D_FIBONACCI, ) from .tools import conify, duplicate_along_axes, epify, precess, stack from .trajectory2D import initialize_2D_radial, initialize_2D_spiral From e72f52e66ea8558b5033f634558a90db1e628bad Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Sun, 3 Nov 2024 23:19:47 +0100 Subject: [PATCH 03/18] fixup! Add densities, random walk and travelling salesman trajectories --- src/mrinufft/trajectories/inits/travelling_salesman.py | 4 ++-- src/mrinufft/trajectories/maths/__init__.py | 7 ++++++- src/mrinufft/trajectories/maths/constants.py | 6 ------ 3 files changed, 8 insertions(+), 9 deletions(-) delete mode 100644 src/mrinufft/trajectories/maths/constants.py diff --git a/src/mrinufft/trajectories/inits/travelling_salesman.py b/src/mrinufft/trajectories/inits/travelling_salesman.py index 0cb69ba2..1e14e1ad 100644 --- a/src/mrinufft/trajectories/inits/travelling_salesman.py +++ b/src/mrinufft/trajectories/inits/travelling_salesman.py @@ -106,7 +106,7 @@ def _initialize_ND_travelling_salesman( # Select k-space locations density = density / np.sum(density) - locations = sample_from_density(Nc * Ns, density, method=sampling) + locations = sample_from_density(Nc * nb_tsp_points, density, method=sampling) # Re-organise locations into Nc clusters if first_cluster_by: @@ -130,7 +130,7 @@ def _initialize_ND_travelling_salesman( # Compute TSP solution over the whole cloud order = solve_tsp_with_2opt(locations, improvement_threshold=tsp_tol) locations = locations[order] - locations = locations.reshape((Nc, Ns, Nd)) + locations = locations.reshape((Nc, nb_tsp_points, Nd)) # Interpolate shot points up to full length trajectory = np.zeros((Nc, Ns, Nd)) diff --git a/src/mrinufft/trajectories/maths/__init__.py b/src/mrinufft/trajectories/maths/__init__.py index 36537a5e..2d49c099 100644 --- a/src/mrinufft/trajectories/maths/__init__.py +++ b/src/mrinufft/trajectories/maths/__init__.py @@ -1,6 +1,8 @@ """Utility module for mathematical operations.""" -from .constants import CIRCLE_PACKING_DENSITY, EIGENVECTOR_2D_FIBONACCI +# Constants +import numpy as np + from .fibonacci import ( generate_fibonacci_circle, generate_fibonacci_lattice, @@ -11,3 +13,6 @@ from .primes import compute_coprime_factors from .rotations import R2D, Ra, Rv, Rx, Ry, Rz from .tsp_solver import solve_tsp_with_2opt + +CIRCLE_PACKING_DENSITY = np.pi / (2 * np.sqrt(3)) +EIGENVECTOR_2D_FIBONACCI = (0.4656, 0.6823, 1) diff --git a/src/mrinufft/trajectories/maths/constants.py b/src/mrinufft/trajectories/maths/constants.py deleted file mode 100644 index 179dbf05..00000000 --- a/src/mrinufft/trajectories/maths/constants.py +++ /dev/null @@ -1,6 +0,0 @@ -"""Constant values for mathematical purposes.""" - -import numpy as np - -CIRCLE_PACKING_DENSITY = np.pi / (2 * np.sqrt(3)) -EIGENVECTOR_2D_FIBONACCI = (0.4656, 0.6823, 1) From 166590678f3dbb1d3f47baf5c814dff18193ec54 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Sun, 3 Nov 2024 23:20:45 +0100 Subject: [PATCH 04/18] Add docstrings for random walk and TSP solver --- src/mrinufft/trajectories/densities.py | 2 + .../trajectories/inits/random_walk.py | 58 ++++++++++++++++++- .../trajectories/inits/travelling_salesman.py | 2 +- src/mrinufft/trajectories/maths/tsp_solver.py | 22 +++++++ src/mrinufft/trajectories/trajectory2D.py | 5 ++ 5 files changed, 85 insertions(+), 4 deletions(-) diff --git a/src/mrinufft/trajectories/densities.py b/src/mrinufft/trajectories/densities.py index 6c0cb8ae..9755d7d3 100644 --- a/src/mrinufft/trajectories/densities.py +++ b/src/mrinufft/trajectories/densities.py @@ -1,3 +1,5 @@ +"""Sampling densities and methods.""" + import numpy as np import numpy.fft as nf import numpy.linalg as nl diff --git a/src/mrinufft/trajectories/inits/random_walk.py b/src/mrinufft/trajectories/inits/random_walk.py index 01ab9962..f91765c5 100644 --- a/src/mrinufft/trajectories/inits/random_walk.py +++ b/src/mrinufft/trajectories/inits/random_walk.py @@ -25,7 +25,7 @@ def _get_diagonal_neighbors_flat_offsets(shape): return neighbors -def _initialize_ND_random_walk(Nc, Ns, density, diagonals=True, pseudo_random=True): +def _initialize_ND_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random=True): flat_density = np.copy(density.flatten()) max_id = np.prod(density.shape) mask = np.ones_like(flat_density) @@ -79,7 +79,33 @@ def _initialize_ND_random_walk(Nc, Ns, density, diagonals=True, pseudo_random=Tr return trajectory -def initialize_2D_random_walk(Nc, Ns, density, diagonals=True, pseudo_random=True): +def initialize_2D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random=True): + """Initialize a 2D random walk trajectory. + + It creates a trajectory by walking randomly to neighboring points + following a provided sampling density. + + Parameters + ---------- + Nc : int + Number of shots + Ns : int + Number of samples per shot + density : array_like + Sampling density used to determine the walk probabilities. + diagonals : bool, optional + Whether to draw the next walk step from the diagional neighbors + on top of the adjacent ones. Default to True. + pseudo_random : bool, optional + Whether to adapt the density dynamically to reduce areas + already covered. The density is still statistically followed + for undersampled acquisitions. Default to True. + + Returns + ------- + array_like + 2D random walk trajectory + """ if len(density.shape) != 2: raise ValueError("`density` is expected to be 2-dimensional.") return _initialize_ND_random_walk( @@ -87,7 +113,33 @@ def initialize_2D_random_walk(Nc, Ns, density, diagonals=True, pseudo_random=Tru ) -def initialize_3D_random_walk(Nc, Ns, density, diagonals=True, pseudo_random=True): +def initialize_3D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random=True): + """Initialize a 3D random walk trajectory. + + It creates a trajectory by walking randomly to neighboring points + following a provided sampling density. + + Parameters + ---------- + Nc : int + Number of shots + Ns : int + Number of samples per shot + density : array_like + Sampling density used to determine the walk probabilities. + diagonals : bool, optional + Whether to draw the next walk step from the diagional neighbors + on top of the adjacent ones. Default to True. + pseudo_random : bool, optional + Whether to adapt the density dynamically to reduce areas + already covered. The density is still statistically followed + for undersampled acquisitions. Default to True. + + Returns + ------- + array_like + 3D random walk trajectory + """ if len(density.shape) != 3: raise ValueError("`density` is expected to be 3-dimensional.") return _initialize_ND_random_walk( diff --git a/src/mrinufft/trajectories/inits/travelling_salesman.py b/src/mrinufft/trajectories/inits/travelling_salesman.py index 1e14e1ad..cda20659 100644 --- a/src/mrinufft/trajectories/inits/travelling_salesman.py +++ b/src/mrinufft/trajectories/inits/travelling_salesman.py @@ -100,7 +100,7 @@ def _initialize_ND_travelling_salesman( # Check arguments validity if Nc * nb_tsp_points > np.prod(density.shape): raise ValueError( - "`density` array not large enough to peak `Nc` * `nb_tsp_points` points." + "`density` array not large enough to pick `Nc` * `nb_tsp_points` points." ) Nd = len(density.shape) diff --git a/src/mrinufft/trajectories/maths/tsp_solver.py b/src/mrinufft/trajectories/maths/tsp_solver.py index 092b8b7a..9b254eb8 100644 --- a/src/mrinufft/trajectories/maths/tsp_solver.py +++ b/src/mrinufft/trajectories/maths/tsp_solver.py @@ -4,6 +4,28 @@ def solve_tsp_with_2opt(locations, improvement_threshold=1e-8): + """Solve the TSP problem using a 2-opt approach. + + A sub-optimal solution to the Travelling Salesman Problem (TSP) + is provided using the 2-opt approach in O(n²) where chunks of + an initially random route are reversed, and selected if the + total distance is reduced. As a result the route solution + does not cross its own path in 2D. + + Parameters + ---------- + locations : array_like + An array of N points with shape (N, D) with D + the space dimension. + improvement_threshold: float + Threshold used as progress criterion to stop the optimization + process. + + Returns + ------- + array_like + The new positions order of shape (N,). + """ route = np.arange(locations.shape[0]) distances = np.linalg.norm(locations[None] - locations[:, None], axis=-1) route_length = np.sum(distances[0]) diff --git a/src/mrinufft/trajectories/trajectory2D.py b/src/mrinufft/trajectories/trajectory2D.py index f87c6362..e930e90b 100644 --- a/src/mrinufft/trajectories/trajectory2D.py +++ b/src/mrinufft/trajectories/trajectory2D.py @@ -27,6 +27,11 @@ def initialize_2D_radial(Nc, Ns, tilt="uniform", in_out=False): Tilt of the shots, by default "uniform" in_out : bool, optional Whether to start from the center or not, by default False + + Returns + ------- + array_like + 2D radial trajectory """ # Initialize a first shot segment = np.linspace(-1 if (in_out) else 0, 1, Ns) From 55b3ba4cd1ca7ca22bd1cda2c63702dc97123d9b Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Thu, 7 Nov 2024 22:30:41 +0100 Subject: [PATCH 05/18] fixup! Add densities, random walk and travelling salesman trajectories --- src/mrinufft/trajectories/__init__.py | 16 ++++++------- .../trajectories/inits/travelling_salesman.py | 23 +++++++++++++------ .../{densities.py => sampling_densities.py} | 8 ++++--- src/mrinufft/trajectories/tools.py | 7 +++++- 4 files changed, 35 insertions(+), 19 deletions(-) rename src/mrinufft/trajectories/{densities.py => sampling_densities.py} (95%) diff --git a/src/mrinufft/trajectories/__init__.py b/src/mrinufft/trajectories/__init__.py index 1e823c03..7051c64f 100644 --- a/src/mrinufft/trajectories/__init__.py +++ b/src/mrinufft/trajectories/__init__.py @@ -1,13 +1,5 @@ """Collection of trajectories and tools used for non-Cartesian MRI.""" -from .densities import ( - create_chauffert_density, - create_cutoff_decay_density, - create_energy_density, - create_fast_chauffert_density, - create_polynomial_density, - sample_from_density, -) from .display import display_2D_trajectory, display_3D_trajectory, displayConfig from .gradients import patch_center_anomaly from .inits import ( @@ -16,6 +8,14 @@ initialize_3D_random_walk, initialize_3D_travelling_salesman, ) +from .sampling_densities import ( + create_chauffert_density, + create_cutoff_decay_density, + create_energy_density, + create_fast_chauffert_density, + create_polynomial_density, + sample_from_density, +) from .tools import ( conify, duplicate_along_axes, diff --git a/src/mrinufft/trajectories/inits/travelling_salesman.py b/src/mrinufft/trajectories/inits/travelling_salesman.py index cda20659..97667f31 100644 --- a/src/mrinufft/trajectories/inits/travelling_salesman.py +++ b/src/mrinufft/trajectories/inits/travelling_salesman.py @@ -1,15 +1,17 @@ -"""Trajectories based on the Travelig Salesman Problem.""" +"""Trajectories based on the Travelling Salesman Problem.""" import numpy as np import numpy.linalg as nl from scipy.interpolate import CubicSpline from tqdm.auto import tqdm -from ..densities import sample_from_density from ..maths import solve_tsp_with_2opt +from ..sampling_densities import sample_from_density +from ..tools import oversample def _get_approx_cluster_sizes(nb_total, nb_clusters): + # Give a list of cluster sizes close to sqrt(`nb_total`) cluster_sizes = round(nb_total / nb_clusters) * np.ones(nb_clusters).astype(int) delta_sum = nb_total - np.sum(cluster_sizes) cluster_sizes[: int(np.abs(delta_sum))] += np.sign(delta_sum) @@ -92,6 +94,8 @@ def _initialize_ND_travelling_salesman( nb_tsp_points="auto", sampling="random", tsp_tol=1e-8, + *, + mask_density=False, verbose=False, ): # Handle variable inputs @@ -105,7 +109,8 @@ def _initialize_ND_travelling_salesman( Nd = len(density.shape) # Select k-space locations - density = density / np.sum(density) + if mask_density: + density = density ** (Nd / (Nd - 1)) locations = sample_from_density(Nc * nb_tsp_points, density, method=sampling) # Re-organise locations into Nc clusters @@ -133,10 +138,8 @@ def _initialize_ND_travelling_salesman( locations = locations.reshape((Nc, nb_tsp_points, Nd)) # Interpolate shot points up to full length - trajectory = np.zeros((Nc, Ns, Nd)) - for i in range(Nc): - cbs = CubicSpline(np.linspace(0, 1, nb_tsp_points), locations[i]) - trajectory[i] = cbs(np.linspace(0, 1, Ns)) + trajectory = oversample(locations, Ns) + return trajectory @@ -150,6 +153,8 @@ def initialize_2D_travelling_salesman( nb_tsp_points="auto", sampling="random", tsp_tol=1e-8, + *, + mask_density=False, verbose=False, ): if len(density.shape) != 2: @@ -164,6 +169,7 @@ def initialize_2D_travelling_salesman( nb_tsp_points=nb_tsp_points, sampling=sampling, tsp_tol=tsp_tol, + mask_density=mask_density, verbose=verbose, ) @@ -178,6 +184,8 @@ def initialize_3D_travelling_salesman( nb_tsp_points="auto", sampling="random", tsp_tol=1e-8, + *, + mask_density=False, verbose=False, ): if len(density.shape) != 3: @@ -192,5 +200,6 @@ def initialize_3D_travelling_salesman( nb_tsp_points=nb_tsp_points, sampling=sampling, tsp_tol=tsp_tol, + mask_density=mask_density, verbose=verbose, ) diff --git a/src/mrinufft/trajectories/densities.py b/src/mrinufft/trajectories/sampling_densities.py similarity index 95% rename from src/mrinufft/trajectories/densities.py rename to src/mrinufft/trajectories/sampling_densities.py index 9755d7d3..0fb0ba30 100644 --- a/src/mrinufft/trajectories/densities.py +++ b/src/mrinufft/trajectories/sampling_densities.py @@ -14,7 +14,8 @@ def sample_from_density(nb_samples, density, method="random"): rng = nr.default_rng() - shape = density.shape + density = density / np.sum(density) + shape = np.array(density.shape) nb_dims = len(shape) max_nb_samples = np.prod(shape) @@ -29,8 +30,9 @@ def sample_from_density(nb_samples, density, method="random"): replace=False, ) locations = np.indices(shape).reshape((nb_dims, -1))[:, choices] - locations = locations.T - locations = 2 * KMAX * locations / np.max(locations) - KMAX + locations = locations.T + 0.5 + locations = locations / shape[None, :] + locations = 2 * KMAX * locations - KMAX elif method == "lloyd": kmeans = ( KMeans(n_clusters=nb_samples) diff --git a/src/mrinufft/trajectories/tools.py b/src/mrinufft/trajectories/tools.py index 5897b7f9..00a1f2c2 100644 --- a/src/mrinufft/trajectories/tools.py +++ b/src/mrinufft/trajectories/tools.py @@ -1,7 +1,7 @@ """Functions to manipulate/modify trajectories.""" import numpy as np -from scipy.interpolate import CubicSpline +from scipy.interpolate import CubicSpline, interp1d from .maths import Rv, Rx, Ry, Rz from .utils import KMAX, initialize_tilt @@ -429,6 +429,11 @@ def rewind(trajectory, Ns_transitions): return assembled_trajectory +def oversample(trajectory, new_Ns, kind="cubic"): + f = interp1d(np.linspace(0, 1, trajectory.shape[1]), trajectory, axis=1, kind=kind) + return f(np.linspace(0, 1, new_Ns)) + + #################### # FUNCTIONAL TOOLS # #################### From 74b08f86376818787531c107edad0a01e59b67e2 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Mon, 11 Nov 2024 15:20:03 +0100 Subject: [PATCH 06/18] Add docstrings to TSP & sampling densities --- src/mrinufft/trajectories/__init__.py | 2 +- .../trajectories/inits/random_walk.py | 34 +- .../trajectories/inits/travelling_salesman.py | 158 ++++++-- src/mrinufft/trajectories/sampling.py | 338 ++++++++++++++++++ .../trajectories/sampling_densities.py | 136 ------- src/mrinufft/trajectories/tools.py | 35 ++ 6 files changed, 524 insertions(+), 179 deletions(-) create mode 100644 src/mrinufft/trajectories/sampling.py delete mode 100644 src/mrinufft/trajectories/sampling_densities.py diff --git a/src/mrinufft/trajectories/__init__.py b/src/mrinufft/trajectories/__init__.py index 7051c64f..719660e3 100644 --- a/src/mrinufft/trajectories/__init__.py +++ b/src/mrinufft/trajectories/__init__.py @@ -8,7 +8,7 @@ initialize_3D_random_walk, initialize_3D_travelling_salesman, ) -from .sampling_densities import ( +from .sampling import ( create_chauffert_density, create_cutoff_decay_density, create_energy_density, diff --git a/src/mrinufft/trajectories/inits/random_walk.py b/src/mrinufft/trajectories/inits/random_walk.py index f91765c5..2671a790 100644 --- a/src/mrinufft/trajectories/inits/random_walk.py +++ b/src/mrinufft/trajectories/inits/random_walk.py @@ -82,9 +82,15 @@ def _initialize_ND_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random def initialize_2D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random=True): """Initialize a 2D random walk trajectory. + This is an adaptation of the proposition from [Cha+14]_. It creates a trajectory by walking randomly to neighboring points following a provided sampling density. + This implementation is different from the original proposition: + trajectories are continuous with a fixed length instead of + making random jumps to other locations, and an option + is provided to have pseudo-random walks to improve coverage. + Parameters ---------- Nc : int @@ -95,16 +101,23 @@ def initialize_2D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= Sampling density used to determine the walk probabilities. diagonals : bool, optional Whether to draw the next walk step from the diagional neighbors - on top of the adjacent ones. Default to True. + on top of the adjacent ones. Default to ``True``. pseudo_random : bool, optional Whether to adapt the density dynamically to reduce areas already covered. The density is still statistically followed - for undersampled acquisitions. Default to True. + for undersampled acquisitions. Default to ``True``. Returns ------- array_like 2D random walk trajectory + + References + ---------- + .. [Cha+14] Chauffert, Nicolas, Philippe Ciuciu, + Jonas Kahn, and Pierre Weiss. + "Variable density sampling with continuous trajectories." + SIAM Journal on Imaging Sciences 7, no. 4 (2014): 1962-1992. """ if len(density.shape) != 2: raise ValueError("`density` is expected to be 2-dimensional.") @@ -116,9 +129,15 @@ def initialize_2D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= def initialize_3D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random=True): """Initialize a 3D random walk trajectory. + This is an adaptation of the proposition from [Cha+14]_. It creates a trajectory by walking randomly to neighboring points following a provided sampling density. + This implementation is different from the original proposition: + trajectories are continuous with a fixed length instead of + making random jumps to other locations, and an option + is provided to have pseudo-random walks to improve coverage. + Parameters ---------- Nc : int @@ -129,16 +148,23 @@ def initialize_3D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= Sampling density used to determine the walk probabilities. diagonals : bool, optional Whether to draw the next walk step from the diagional neighbors - on top of the adjacent ones. Default to True. + on top of the adjacent ones. Default to ``True``. pseudo_random : bool, optional Whether to adapt the density dynamically to reduce areas already covered. The density is still statistically followed - for undersampled acquisitions. Default to True. + for undersampled acquisitions. Default to ``True``. Returns ------- array_like 3D random walk trajectory + + References + ---------- + .. [Cha+14] Chauffert, Nicolas, Philippe Ciuciu, + Jonas Kahn, and Pierre Weiss. + "Variable density sampling with continuous trajectories." + SIAM Journal on Imaging Sciences 7, no. 4 (2014): 1962-1992. """ if len(density.shape) != 3: raise ValueError("`density` is expected to be 3-dimensional.") diff --git a/src/mrinufft/trajectories/inits/travelling_salesman.py b/src/mrinufft/trajectories/inits/travelling_salesman.py index 97667f31..590ede8e 100644 --- a/src/mrinufft/trajectories/inits/travelling_salesman.py +++ b/src/mrinufft/trajectories/inits/travelling_salesman.py @@ -6,7 +6,7 @@ from tqdm.auto import tqdm from ..maths import solve_tsp_with_2opt -from ..sampling_densities import sample_from_density +from ..sampling import sample_from_density from ..tools import oversample @@ -19,6 +19,7 @@ def _get_approx_cluster_sizes(nb_total, nb_clusters): def _sort_by_coordinate(array, coord): + # Sort a list of N-D locations by a Cartesian/spherical coordinate if array.shape[-1] < 3 and coord.lower() in ["z", "theta"]: raise ValueError( f"Invalid `coord`='{coord}' for arrays with less than 3 dimensions." @@ -48,6 +49,7 @@ def _sort_by_coordinate(array, coord): def _cluster_by_coordinate( locations, nb_clusters, cluster_by, second_cluster_by=None, sort_by=None ): + # Cluster approximately a list of N-D locations by Cartesian/spherical coordinates # Gather dimension variables nb_dims = locations.shape[-1] locations = locations.reshape((-1, nb_dims)) @@ -91,32 +93,23 @@ def _initialize_ND_travelling_salesman( first_cluster_by=None, second_cluster_by=None, sort_by=None, - nb_tsp_points="auto", - sampling="random", tsp_tol=1e-8, *, - mask_density=False, verbose=False, + **sampling_kwargs, ): - # Handle variable inputs - nb_tsp_points = Ns if nb_tsp_points == "auto" else nb_tsp_points - # Check arguments validity - if Nc * nb_tsp_points > np.prod(density.shape): - raise ValueError( - "`density` array not large enough to pick `Nc` * `nb_tsp_points` points." - ) + if Nc * Ns > np.prod(density.shape): + raise ValueError("`density` array not large enough to pick `Nc` * `Ns` points.") Nd = len(density.shape) # Select k-space locations - if mask_density: - density = density ** (Nd / (Nd - 1)) - locations = sample_from_density(Nc * nb_tsp_points, density, method=sampling) + trajectory = sample_from_density(Nc * Ns, density, **sampling_kwargs) # Re-organise locations into Nc clusters if first_cluster_by: - locations = _cluster_by_coordinate( - locations, + trajectory = _cluster_by_coordinate( + trajectory, Nc, cluster_by=first_cluster_by, second_cluster_by=second_cluster_by, @@ -125,20 +118,17 @@ def _initialize_ND_travelling_salesman( # Compute TSP solution within each cluster/shot for i in tqdm(range(Nc), disable=not verbose): - order = solve_tsp_with_2opt(locations[i], improvement_threshold=tsp_tol) - locations[i] = locations[i][order] + order = solve_tsp_with_2opt(trajectory[i], improvement_threshold=tsp_tol) + trajectory[i] = trajectory[i][order] else: - locations = ( - _sort_by_coordinate(locations, coord=sort_by) if sort_by else locations + trajectory = ( + _sort_by_coordinate(trajectory, coord=sort_by) if sort_by else trajectory ) # Compute TSP solution over the whole cloud - order = solve_tsp_with_2opt(locations, improvement_threshold=tsp_tol) - locations = locations[order] - locations = locations.reshape((Nc, nb_tsp_points, Nd)) - - # Interpolate shot points up to full length - trajectory = oversample(locations, Ns) + order = solve_tsp_with_2opt(trajectory, improvement_threshold=tsp_tol) + trajectory = trajectory[order] + trajectory = trajectory.reshape((Nc, Ns, Nd)) return trajectory @@ -150,13 +140,61 @@ def initialize_2D_travelling_salesman( first_cluster_by=None, second_cluster_by=None, sort_by=None, - nb_tsp_points="auto", - sampling="random", tsp_tol=1e-8, *, - mask_density=False, verbose=False, + **sampling_kwargs, ): + """ + Initialize a 2D trajectory using a Travelling Salesman Problem (TSP)-based path. + + This is a reproduction of the work from [Cha+14]_. The TSP solution + is obtained using the 2-opt method in O(n²). An additional option + is provided to cluster shots before solving the TSP and thus + reduce drastically the computation time. The initial sampling method + can also be customized. + + Parameters + ---------- + Nc : int + The number of clusters (or shots) to divide the trajectory into. + Ns : int + The number of points per cluster. + density : np.ndarray + A 2-dimensional density array from which points are sampled. + first_cluster_by : str, optional + The coordinate used to cluster points initially, by default ``None``. + second_cluster_by : str, optional + A secondary coordinate used for clustering within primary clusters, + by default ``None``. + sort_by : str, optional + The coordinate by which to order points within each cluster, + by default ``None``. + tsp_tol : float, optional + Convergence tolerance for the TSP solution, by default ``1e-8``. + verbose : bool, optional + If ``True``, displays a progress bar, by default ``False``. + **sampling_kwargs : dict, optional + Additional arguments to pass to + ``mrinufft.trajectories.sampling.sample_from_density``. + + Returns + ------- + np.ndarray + A 2D array representing the TSP-ordered trajectory. + + Raises + ------ + ValueError + If ``density`` is not a 2-dimensional array. + + References + ---------- + .. [Cha+14] Chauffert, Nicolas, Philippe Ciuciu, + Jonas Kahn, and Pierre Weiss. + "Variable density sampling with continuous trajectories." + SIAM Journal on Imaging Sciences 7, no. 4 (2014): 1962-1992. + """ if len(density.shape) != 2: raise ValueError("`density` is expected to be 2-dimensional.") return _initialize_ND_travelling_salesman( @@ -166,11 +204,9 @@ def initialize_2D_travelling_salesman( first_cluster_by=first_cluster_by, second_cluster_by=second_cluster_by, sort_by=sort_by, - nb_tsp_points=nb_tsp_points, - sampling=sampling, tsp_tol=tsp_tol, - mask_density=mask_density, verbose=verbose, + **sampling_kwargs, ) @@ -181,13 +217,61 @@ def initialize_3D_travelling_salesman( first_cluster_by=None, second_cluster_by=None, sort_by=None, - nb_tsp_points="auto", - sampling="random", tsp_tol=1e-8, *, - mask_density=False, verbose=False, + **sampling_kwargs, ): + """ + Initialize a 3D trajectory using a Travelling Salesman Problem (TSP)-based path. + + This is a reproduction of the work from [Cha+14]_. The TSP solution + is obtained using the 2-opt method in O(n²). An additional option + is provided to cluster shots before solving the TSP and thus + reduce drastically the computation time. The initial sampling method + can also be customized. + + Parameters + ---------- + Nc : int + The number of clusters (or shots) to divide the trajectory into. + Ns : int + The number of points per cluster. + density : np.ndarray + A 3-dimensional density array from which points are sampled. + first_cluster_by : str, optional + The coordinate used to cluster points initially, by default ``None``. + second_cluster_by : str, optional + A secondary coordinate used for clustering within primary clusters, + by default ``None``. + sort_by : str, optional + The coordinate by which to order points within each cluster, + by default ``None``. + tsp_tol : float, optional + Convergence tolerance for the TSP solution, by default ``1e-8``. + verbose : bool, optional + If ``True``, displays a progress bar, by default ``False``. + **sampling_kwargs : dict, optional + Additional arguments to pass to + ``mrinufft.trajectories.sampling.sample_from_density``. + + Returns + ------- + np.ndarray + A 3D array representing the TSP-ordered trajectory. + + Raises + ------ + ValueError + If ``density`` is not a 3-dimensional array. + + References + ---------- + .. [Cha+14] Chauffert, Nicolas, Philippe Ciuciu, + Jonas Kahn, and Pierre Weiss. + "Variable density sampling with continuous trajectories." + SIAM Journal on Imaging Sciences 7, no. 4 (2014): 1962-1992. + """ if len(density.shape) != 3: raise ValueError("`density` is expected to be 3-dimensional.") return _initialize_ND_travelling_salesman( @@ -197,9 +281,7 @@ def initialize_3D_travelling_salesman( first_cluster_by=first_cluster_by, second_cluster_by=second_cluster_by, sort_by=sort_by, - nb_tsp_points=nb_tsp_points, - sampling=sampling, tsp_tol=tsp_tol, - mask_density=mask_density, verbose=verbose, + **sampling_kwargs, ) diff --git a/src/mrinufft/trajectories/sampling.py b/src/mrinufft/trajectories/sampling.py new file mode 100644 index 00000000..2e1378fb --- /dev/null +++ b/src/mrinufft/trajectories/sampling.py @@ -0,0 +1,338 @@ +"""Sampling densities and methods.""" + +import numpy as np +import numpy.fft as nf +import numpy.linalg as nl +import numpy.random as nr +import pywt as pw +from sklearn.cluster import BisectingKMeans, KMeans +from tqdm.auto import tqdm + +from .utils import KMAX + + +def sample_from_density( + nb_samples, density, method="random", *, binary_mask_density=False +): + """ + Sample points based on a given density distribution. + + Parameters + ---------- + nb_samples : int + The number of samples to draw. + density : np.ndarray + An array representing the density distribution from which samples are drawn. + method : str, optional + The sampling method to use, either 'random' for random sampling over + the discrete grid defined by the density or 'lloyd' for Lloyd's + algorithm over a continuous space, by default "random". + binary_mask_density : bool, optional + Whether to apply a specific density correction introduced + in [Cha+14]_. An exponent ``N/(N-1)`` with ``N`` the number of + dimensions in ``density`` is applied to fix the observed + density expectation when paths are traced between drawn samples + specifically over binary masks/grids. Otherwise it is demonstrated + that the empirical density will be ``density ** (N - 1) / N``. + + Returns + ------- + np.ndarray + An array of range-normalized sampled locations. + + Raises + ------ + ValueError + If ``nb_samples`` exceeds the total size of the density array or if the + specified ``method`` is unknown. + + + References + ---------- + .. [Cha+14] Chauffert, Nicolas, Philippe Ciuciu, + Jonas Kahn, and Pierre Weiss. + "Variable density sampling with continuous trajectories." + SIAM Journal on Imaging Sciences 7, no. 4 (2014): 1962-1992. + """ + rng = nr.default_rng() + + # Define dimension variables + shape = np.array(density.shape) + nb_dims = len(shape) + max_nb_samples = np.prod(shape) + + if nb_samples > max_nb_samples: + raise ValueError("`nb_samples` must be lower than the size of `density`.") + + # Adapt density to use binary mask use case + density = density ** (nb_dims / (nb_dims - 1)) if binary_mask_density else density + density = density / np.sum(density) + + # Sample using specified method + if method == "random": + choices = rng.choice( + np.arange(max_nb_samples), + size=nb_samples, + p=density.flatten(), + replace=False, + ) + locations = np.indices(shape).reshape((nb_dims, -1))[:, choices] + locations = locations.T + 0.5 + locations = locations / shape[None, :] + locations = 2 * KMAX * locations - KMAX + elif method == "lloyd": + kmeans = ( + KMeans(n_clusters=nb_samples) + if nb_dims <= 2 + else BisectingKMeans(n_clusters=nb_samples) + ) + kmeans.fit( + np.indices(density.shape).reshape((nb_dims, -1)).T, + sample_weight=density.flatten(), + ) + locations = kmeans.cluster_centers_ - np.array(density.shape) / 2 + locations = KMAX * locations / np.max(np.abs(locations)) + else: + raise ValueError(f"Unknown sampling method {method}.") + return locations + + +def create_cutoff_decay_density(shape, cutoff, decay, resolution=None): + """ + Create a density with central plateau and polynomial decay. + + Create a density composed of a central constant-value ellipsis + defined by a cutoff ratio, followed by a polynomial decay over + outer regions as defined in [Cha+22]_. + + Parameters + ---------- + shape : tuple of int + The shape of the density grid. + cutoff : float + The k-space radius cutoff ratio between 0 and 1 within + which density remains uniform and beyond which it decays. + decay : float + The rate of decay of density beyond the cutoff ratio. + resolution : np.ndarray, optional + Resolution scaling factors for each dimension of the density grid, + by default None. + + Returns + ------- + np.ndarray + A density array with values decaying based on the specified + cutoff ratio and decay rate. + + References + ---------- + .. [Cha+22] Chaithya, G. R., Pierre Weiss, Guillaume Daval-Frérot, + Aurélien Massire, Alexandre Vignaud, and Philippe Ciuciu. + "Optimizing full 3d sparkling trajectories for high-resolution + magnetic resonance imaging." + IEEE Transactions on Medical Imaging 41, no. 8 (2022): 2105-2117. + """ + shape = np.array(shape) + nb_dims = len(shape) + + if not resolution: + resolution = np.ones(nb_dims) + + distances = np.indices(shape).astype(float) + for i in range(nb_dims): + distances[i] = distances[i] + 0.5 - shape[i] / 2 + distances[i] = distances[i] / shape[i] * resolution[i] + distances = nl.norm(distances, axis=0) + + cutoff = cutoff * np.max(distances) if cutoff else np.min(distances) + density = np.ones(shape) + decay_mask = np.where(distances > cutoff, True, False) + + density[decay_mask] = (cutoff / distances[decay_mask]) ** decay + density = density / np.sum(density) + return density + + +def create_polynomial_density(shape, decay, resolution=None): + """ + Create a density with polynomial decay from the center. + + Parameters + ---------- + shape : tuple of int + The shape of the density grid. + decay : float + The exponent that controls the rate of decay for density. + resolution : np.ndarray, optional + Resolution scaling factors for each dimension of the density grid, + by default None. + + Returns + ------- + np.ndarray + A density array with polynomial decay. + """ + return create_cutoff_decay_density( + shape, cutoff=0, decay=decay, resolution=resolution + ) + + +def create_energy_density(dataset): + """ + Create a density based on energy in the Fourier spectrum. + + A density is made based on the average energy in the Fourier domain + of volumes from a target image dataset. + + Parameters + ---------- + dataset : np.ndarray + The dataset from which to calculate the density + based on its Fourier transform, with an expected + shape (nb_volumes, dim_1, ..., dim_N). + An N-dimensional Fourier transform is performed. + + Returns + ------- + np.ndarray + A density array derived from the mean energy in the Fourier + domain of the input dataset. + """ + nb_dims = len(dataset.shape) - 1 + axes = range(1, nb_dims + 1) + + kspace = nf.fftshift(nf.fftn(nf.fftshift(dataset, axes=axes), axes=axes), axes=axes) + energy = np.mean(np.abs(kspace), axis=0) + density = energy / np.sum(energy) + return density + + +def create_chauffert_density(shape, wavelet_basis, nb_wavelet_scales, verbose=False): + """Create a density based on Chauffert's method. + + This is a reproduction of the proposition from [CCW13]_. + A sampling density is derived from compressed sensing + equations to maximize guarantees of exact image recovery + for a specified sparse wavelet domain decomposition. + + Parameters + ---------- + shape : tuple of int + The shape of the density grid. + wavelet_basis : str, pywt.Wavelet + The wavelet basis to use for wavelet decomposition, either + as a built-in wavelet name from the PyWavelets package + or as a custom ``pywt.Wavelet`` object. + nb_wavelet_scales : int + The number of wavelet scales to use in decomposition. + verbose : bool, optional + If ``True``, displays a progress bar. Default to ``False``. + + Returns + ------- + np.ndarray + A density array created based on wavelet transform coefficients. + + See Also + -------- + pywt.wavelist : A list of wavelet decompositions available in the + PyWavelets package used inside the function. + pywt.Wavelet : A wavelet object accepted to generate Chauffert densities. + + References + ---------- + .. [CCW13] Chauffert, Nicolas, Philippe Ciuciu, and Pierre Weiss. + "Variable density compressed sensing in MRI. + Theoretical vs heuristic sampling strategies." + In 2013 IEEE 10th International Symposium on Biomedical Imaging, + pp. 298-301. IEEE, 2013. + """ + nb_dims = len(shape) + indices = np.indices(shape).reshape((nb_dims, -1)).T + + density = np.zeros(shape) + unit_vector = np.zeros(shape) + for ids in tqdm(indices, disable=not verbose): + ids = tuple(ids) + unit_vector[ids] = 1 + fourier_vector = nf.ifftn(unit_vector) + coeffs = pw.wavedecn( + fourier_vector, wavelet=wavelet_basis, level=nb_wavelet_scales + ) + coeffs, _ = pw.coeffs_to_array(coeffs) + density[ids] = np.max(np.abs(coeffs)) ** 2 + unit_vector[ids] = 0 + + density = density / np.sum(density) + return nf.ifftshift(density) + + +def create_fast_chauffert_density(shape, wavelet_basis, nb_wavelet_scales): + """Create a density based on an approximated Chauffert method. + + This implementation is based on this + tutorial: https://github.com/philouc/mri_acq_recon_tutorial. + It is a fast approximation of the proposition from [CCW13]_, + where a sampling density is derived from compressed sensing + equations to maximize guarantees of exact image recovery + for a specified sparse wavelet domain decomposition. + + In this approximation, the decomposition dimensions are + considered independent and computed separately to accelerate + the density generation. + + Parameters + ---------- + shape : tuple of int + The shape of the density grid. + wavelet_basis : str, pywt.Wavelet + The wavelet basis to use for wavelet decomposition, either + as a built-in wavelet name from the PyWavelets package + or as a custom ``pywt.Wavelet`` object. + nb_wavelet_scales : int + The number of wavelet scales to use in decomposition. + + Returns + ------- + np.ndarray + A density array created using a faster approximation + based on 1D projections of the wavelet transform. + + See Also + -------- + pywt.wavelist : A list of wavelet decompositions available in the + PyWavelets package used inside the function. + pywt.Wavelet : A wavelet object accepted to generate Chauffert densities. + + References + ---------- + .. [CCW13] Chauffert, Nicolas, Philippe Ciuciu, and Pierre Weiss. + "Variable density compressed sensing in MRI. + Theoretical vs heuristic sampling strategies." + In 2013 IEEE 10th International Symposium on Biomedical Imaging, + pp. 298-301. IEEE, 2013. + """ + nb_dims = len(shape) + + density = np.ones(shape) + for k, s in enumerate(shape): + unit_vector = np.zeros(s) + density_1d = np.zeros(s) + + for i in range(s): + unit_vector[i] = 1 + fourier_vector = nf.ifft(unit_vector) + coeffs = pw.wavedec( + fourier_vector, wavelet=wavelet_basis, level=nb_wavelet_scales + ) + coeffs, _ = pw.coeffs_to_array(coeffs) + density_1d[i] = np.max(np.abs(coeffs)) ** 2 + unit_vector[i] = 0 + + reshape = np.ones(nb_dims).astype(int) + reshape[k] = s + density_1d = density_1d.reshape(reshape) + density = density * density_1d + + density = density / np.sum(density) + return nf.ifftshift(density) diff --git a/src/mrinufft/trajectories/sampling_densities.py b/src/mrinufft/trajectories/sampling_densities.py deleted file mode 100644 index 0fb0ba30..00000000 --- a/src/mrinufft/trajectories/sampling_densities.py +++ /dev/null @@ -1,136 +0,0 @@ -"""Sampling densities and methods.""" - -import numpy as np -import numpy.fft as nf -import numpy.linalg as nl -import numpy.random as nr -import pywt as pw -from sklearn.cluster import BisectingKMeans, KMeans -from tqdm.auto import tqdm - -from .utils import KMAX - - -def sample_from_density(nb_samples, density, method="random"): - rng = nr.default_rng() - - density = density / np.sum(density) - shape = np.array(density.shape) - nb_dims = len(shape) - max_nb_samples = np.prod(shape) - - if nb_samples > max_nb_samples: - raise ValueError("`nb_samples` must be lower than the size of `density`.") - - if method == "random": - choices = rng.choice( - np.arange(max_nb_samples), - size=nb_samples, - p=density.flatten(), - replace=False, - ) - locations = np.indices(shape).reshape((nb_dims, -1))[:, choices] - locations = locations.T + 0.5 - locations = locations / shape[None, :] - locations = 2 * KMAX * locations - KMAX - elif method == "lloyd": - kmeans = ( - KMeans(n_clusters=nb_samples) - if nb_dims <= 2 - else BisectingKMeans(n_clusters=nb_samples) - ) - kmeans.fit( - np.indices(density.shape).reshape((nb_dims, -1)).T, - sample_weight=density.flatten(), - ) - locations = kmeans.cluster_centers_ - np.array(density.shape) / 2 - locations = KMAX * locations / np.max(np.abs(locations)) - else: - raise ValueError(f"Unknown sampling method {method}.") - return locations - - -def create_cutoff_decay_density(shape, cutoff, decay, resolution=None): - shape = np.array(shape) - nb_dims = len(shape) - - if not resolution: - resolution = np.ones(nb_dims) - - distances = np.indices(shape).astype(float) - for i in range(nb_dims): - distances[i] = distances[i] + 0.5 - shape[i] / 2 - distances[i] = distances[i] / shape[i] * resolution[i] - distances = nl.norm(distances, axis=0) - - cutoff = cutoff * np.max(distances) if cutoff else np.min(distances) - density = np.ones(shape) - decay_mask = np.where(distances > cutoff, True, False) - - density[decay_mask] = (cutoff / distances[decay_mask]) ** decay - density = density / np.sum(density) - return density - - -def create_polynomial_density(shape, decay, resolution=None): - return create_cutoff_decay_density( - shape, cutoff=0, decay=decay, resolution=resolution - ) - - -def create_energy_density(dataset): - nb_dims = len(dataset.shape) - 1 - axes = range(1, nb_dims + 1) - - kspace = nf.fftshift(nf.fftn(nf.fftshift(dataset, axes=axes), axes=axes), axes=axes) - energy = np.mean(np.abs(kspace), axis=0) - density = energy / np.sum(energy) - return density - - -def create_chauffert_density(shape, wavelet_basis, nb_wavelet_scales, verbose=False): - nb_dims = len(shape) - indices = np.indices(shape).reshape((nb_dims, -1)).T - - density = np.zeros(shape) - unit_vector = np.zeros(shape) - for ids in tqdm(indices, disable=not verbose): - ids = tuple(ids) - unit_vector[ids] = 1 - fourier_vector = nf.ifftn(unit_vector) - coeffs = pw.wavedecn( - fourier_vector, wavelet=wavelet_basis, level=nb_wavelet_scales - ) - coeffs, _ = pw.coeffs_to_array(coeffs) - density[ids] = np.max(np.abs(coeffs)) ** 2 - unit_vector[ids] = 0 - - density = density / np.sum(density) - return nf.ifftshift(density) - - -def create_fast_chauffert_density(shape, wavelet_basis, nb_wavelet_scales): - nb_dims = len(shape) - - density = np.ones(shape) - for k, s in enumerate(shape): - unit_vector = np.zeros(s) - density_1d = np.zeros(s) - - for i in range(s): - unit_vector[i] = 1 - fourier_vector = nf.ifft(unit_vector) - coeffs = pw.wavedec( - fourier_vector, wavelet=wavelet_basis, level=nb_wavelet_scales - ) - coeffs, _ = pw.coeffs_to_array(coeffs) - density_1d[i] = np.max(np.abs(coeffs)) ** 2 - unit_vector[i] = 0 - - reshape = np.ones(nb_dims).astype(int) - reshape[k] = s - density_1d = density_1d.reshape(reshape) - density = density * density_1d - - density = density / np.sum(density) - return nf.ifftshift(density) diff --git a/src/mrinufft/trajectories/tools.py b/src/mrinufft/trajectories/tools.py index 00a1f2c2..cc099c07 100644 --- a/src/mrinufft/trajectories/tools.py +++ b/src/mrinufft/trajectories/tools.py @@ -430,6 +430,41 @@ def rewind(trajectory, Ns_transitions): def oversample(trajectory, new_Ns, kind="cubic"): + """ + Resample a trajectory to increase the number of samples using interpolation. + + Parameters + ---------- + trajectory : np.ndarray + The original trajectory array, where interpolation + is applied along the second axis. + new_Ns : int + The desired number of samples in the resampled trajectory. + kind : str, optional + The type of interpolation to use, such as 'linear', + 'quadratic', or 'cubic', by default "cubic". + + Returns + ------- + np.ndarray + The resampled trajectory array with ``new_Ns`` points along the second axis. + + Notes + ----- + This function uses ``scipy.interpolate.interp1d`` to perform + the interpolation along the second axis of the input `trajectory` array. + + Warnings + -------- + Using 'quadratic' or 'cubic' interpolations is likely to generate + samples located slightly beyond the original k-space limits by + making smooth transitions. + + See Also + -------- + scipy.interpolate.interp1d : The underlying interpolation function + used for resampling. + """ f = interp1d(np.linspace(0, 1, trajectory.shape[1]), trajectory, axis=1, kind=kind) return f(np.linspace(0, 1, new_Ns)) From 362ca9fc11dd69b366cde661fc7235d388713552 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Sat, 23 Nov 2024 23:46:48 +0100 Subject: [PATCH 07/18] Add samplings to random walk initializations and fix boundaries --- src/mrinufft/__init__.py | 39 ++++++- src/mrinufft/trajectories/__init__.py | 16 ++- .../trajectories/inits/random_walk.py | 100 ++++++++++++------ src/mrinufft/trajectories/sampling.py | 29 +++-- 4 files changed, 132 insertions(+), 52 deletions(-) diff --git a/src/mrinufft/__init__.py b/src/mrinufft/__init__.py index 69bae1a2..041c454f 100644 --- a/src/mrinufft/__init__.py +++ b/src/mrinufft/__init__.py @@ -14,6 +14,7 @@ ) from .trajectories import ( + # trajectories initialize_2D_radial, initialize_2D_spiral, initialize_2D_fibonacci_spiral, @@ -25,6 +26,8 @@ initialize_2D_polar_lissajous, initialize_2D_lissajous, initialize_2D_waves, + initialize_2D_random_walk, + initialize_2D_travelling_salesman, initialize_3D_phyllotaxis_radial, initialize_3D_golden_means_radial, initialize_3D_wong_radial, @@ -38,6 +41,9 @@ initialize_3D_seiffert_shells, initialize_3D_turbine, initialize_3D_repi, + initialize_3D_random_walk, + initialize_3D_travelling_salesman, + # tools stack, rotate, precess, @@ -46,6 +52,15 @@ shellify, duplicate_along_axes, radialize_center, + oversample, + # densities + sample_from_density, + create_cutoff_decay_density, + create_polynomial_density, + create_energy_density, + create_chauffert_density, + create_fast_chauffert_density, + # display displayConfig, display_2D_trajectory, display_3D_trajectory, @@ -54,10 +69,16 @@ from .density import voronoi, cell_count, pipe, get_density __all__ = [ + # nufft "get_operator", "check_backend", "list_backends", "get_interpolators_from_fieldmap", + "voronoi", + "cell_count", + "pipe", + "get_density", + # trajectories "initialize_2D_radial", "initialize_2D_spiral", "initialize_2D_fibonacci_spiral", @@ -69,6 +90,8 @@ "initialize_2D_polar_lissajous", "initialize_2D_lissajous", "initialize_2D_waves", + "initialize_2D_random_walk", + "initialize_2D_travelling_salesman", "initialize_3D_phyllotaxis_radial", "initialize_3D_golden_means_radial", "initialize_3D_wong_radial", @@ -83,6 +106,9 @@ "initialize_3D_seiffert_shells", "initialize_3D_turbine", "initialize_3D_repi", + "initialize_3D_random_walk", + "initialize_3D_travelling_salesman", + # tools "stack", "rotate", "precess", @@ -92,13 +118,18 @@ "duplicate_along_axes", "radialize_center", "patch_center_anomaly", + "oversample", + # densities + "sample_from_density", + "create_cutoff_decay_density", + "create_polynomial_density", + "create_energy_density", + "create_chauffert_density", + "create_fast_chauffert_density", + # display "displayConfig", "display_2D_trajectory", "display_3D_trajectory", - "voronoi", - "cell_count", - "pipe", - "get_density", ] from importlib.metadata import version, PackageNotFoundError diff --git a/src/mrinufft/trajectories/__init__.py b/src/mrinufft/trajectories/__init__.py index 719660e3..77c588cd 100644 --- a/src/mrinufft/trajectories/__init__.py +++ b/src/mrinufft/trajectories/__init__.py @@ -19,6 +19,7 @@ from .tools import ( conify, duplicate_along_axes, + oversample, precess, radialize_center, rotate, @@ -56,7 +57,7 @@ ) __all__ = [ - # Trajectories + # trajectories "initialize_2D_radial", "initialize_2D_spiral", "initialize_2D_fibonacci_spiral", @@ -68,6 +69,8 @@ "initialize_2D_polar_lissajous", "initialize_2D_lissajous", "initialize_2D_waves", + "initialize_2D_random_walk", + "initialize_2D_travelling_salesman", "initialize_3D_phyllotaxis_radial", "initialize_3D_golden_means_radial", "initialize_3D_wong_radial", @@ -82,7 +85,9 @@ "initialize_3D_seiffert_shells", "initialize_3D_turbine", "initialize_3D_repi", - # Tools + "initialize_3D_random_walk", + "initialize_3D_travelling_salesman", + # tools "stack", "rotate", "precess", @@ -95,11 +100,16 @@ "display_2D_trajectory", "display_3D_trajectory", "patch_center_anomaly", - # Densities + "oversample", + # densities "sample_from_density", "create_cutoff_decay_density", "create_polynomial_density", "create_energy_density", "create_chauffert_density", "create_fast_chauffert_density", + # display + "displayConfig", + "display_2D_trajectory", + "display_3D_trajectory", ] diff --git a/src/mrinufft/trajectories/inits/random_walk.py b/src/mrinufft/trajectories/inits/random_walk.py index 2671a790..da68feb1 100644 --- a/src/mrinufft/trajectories/inits/random_walk.py +++ b/src/mrinufft/trajectories/inits/random_walk.py @@ -2,55 +2,63 @@ import numpy as np +from ..sampling import sample_from_density from ..utils import KMAX -def _get_adjacent_neighbors_flat_offsets(shape): - nb_dims = len(shape) - neighborhood = np.indices([3] * nb_dims) - 1 - distances = np.sum(np.abs(neighborhood), axis=0) - - center = np.ravel_multi_index([1] * nb_dims, dims=shape) - neighbors = np.ravel_multi_index(np.where(distances == 1), dims=shape) - center - return neighbors +def _get_adjacent_neighbors_offsets(shape): + return np.concatenate([np.eye(len(shape)), -np.eye(len(shape))], axis=0).astype(int) -def _get_diagonal_neighbors_flat_offsets(shape): +def _get_neighbors_offsets(shape): nb_dims = len(shape) - neighborhood = np.indices([3] * nb_dims) - 1 - distances = np.sum(np.abs(neighborhood), axis=0) - - center = np.ravel_multi_index([1] * nb_dims, dims=shape) - neighbors = np.ravel_multi_index(np.where(distances > 1), dims=shape) - center + neighbors = (np.indices([3] * nb_dims) - 1).reshape((nb_dims, -1)).T + nb_half = neighbors.shape[0] // 2 + # Remove full zero entry + neighbors = np.concatenate([neighbors[:nb_half], neighbors[-nb_half:]], axis=0) return neighbors -def _initialize_ND_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random=True): +def _initialize_ND_random_walk( + Nc, Ns, density, *, diagonals=True, pseudo_random=True, **sampling_kwargs +): + density = density / np.sum(density) flat_density = np.copy(density.flatten()) - max_id = np.prod(density.shape) + shape = np.array(density.shape) mask = np.ones_like(flat_density) # Prepare neighbor offsets once - offsets = _get_adjacent_neighbors_flat_offsets(density.shape) - if diagonals: - offsets = np.concatenate( - [offsets, _get_diagonal_neighbors_flat_offsets(density.shape)] - ) + offsets = ( + _get_neighbors_offsets(shape) + if diagonals + else _get_adjacent_neighbors_offsets(shape) + ) # Make all random draws at once for performance draws = np.random.random((Ns, Nc)) # inverted shape for convenience # Initialize shot starting points - choices = np.random.choice(np.arange(len(flat_density)), size=Nc, p=flat_density) + locations = sample_from_density(Nc, density, **sampling_kwargs) + choices = np.around((locations + KMAX) * (np.array(density.shape) - 1)).astype(int) + choices = np.ravel_multi_index(choices.T, density.shape) + # choices = np.random.choice(np.arange(len(flat_density)), size=Nc, p=flat_density) routes = [choices] # Walk for i in range(1, Ns): - neighbors = choices[:, None] + offsets[None] + # Express indices back to multi-dim coordinates to check bounds + neighbors = np.array(np.unravel_index(choices, shape)) + neighbors = neighbors[:, None] + offsets.T[..., None] # Find out-of-bound neighbors and ignore them - invalids = (neighbors < 0) | (neighbors > max_id) - neighbors[invalids] = 0 + invalids = (neighbors < 0).any(axis=0) | ( + neighbors >= shape[:, None, None] + ).any(axis=0) + neighbors[:, invalids] = 0 + invalids = invalids.T + + # Flatten indices to use np.random.choice + neighbors = np.ravel_multi_index(neighbors, shape).T # Set walk probabilities walk_probs = flat_density[neighbors] @@ -72,14 +80,16 @@ def _initialize_ND_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random routes = np.array(routes).T # Create trajectory from routes - locations = np.indices(density.shape) - locations = locations.reshape((len(density.shape), -1)) + locations = np.indices(shape) + locations = locations.reshape((len(shape), -1)) trajectory = np.array([locations[:, r].T for r in routes]) - trajectory = 2 * KMAX * trajectory / density.shape - KMAX + trajectory = 2 * KMAX * trajectory / (shape - 1) - KMAX return trajectory -def initialize_2D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random=True): +def initialize_2D_random_walk( + Nc, Ns, density, *, diagonals=True, pseudo_random=True, **sampling_kwargs +): """Initialize a 2D random walk trajectory. This is an adaptation of the proposition from [Cha+14]_. @@ -98,7 +108,8 @@ def initialize_2D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= Ns : int Number of samples per shot density : array_like - Sampling density used to determine the walk probabilities. + Sampling density used to determine the walk probabilities, + normalized automatically by its sum during the call for convenience. diagonals : bool, optional Whether to draw the next walk step from the diagional neighbors on top of the adjacent ones. Default to ``True``. @@ -106,6 +117,10 @@ def initialize_2D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= Whether to adapt the density dynamically to reduce areas already covered. The density is still statistically followed for undersampled acquisitions. Default to ``True``. + **sampling_kwargs + Sampling parameters in + ``mrinufft.trajectories.sampling.sample_from_density`` used for the + shot starting positions. Returns ------- @@ -122,11 +137,18 @@ def initialize_2D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= if len(density.shape) != 2: raise ValueError("`density` is expected to be 2-dimensional.") return _initialize_ND_random_walk( - Nc, Ns, density, diagonals=diagonals, pseudo_random=pseudo_random + Nc, + Ns, + density, + diagonals=diagonals, + pseudo_random=pseudo_random, + **sampling_kwargs, ) -def initialize_3D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random=True): +def initialize_3D_random_walk( + Nc, Ns, density, *, diagonals=True, pseudo_random=True, **sampling_kwargs +): """Initialize a 3D random walk trajectory. This is an adaptation of the proposition from [Cha+14]_. @@ -145,7 +167,8 @@ def initialize_3D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= Ns : int Number of samples per shot density : array_like - Sampling density used to determine the walk probabilities. + Sampling density used to determine the walk probabilities, + normalized automatically by its sum during the call for convenience. diagonals : bool, optional Whether to draw the next walk step from the diagional neighbors on top of the adjacent ones. Default to ``True``. @@ -153,6 +176,10 @@ def initialize_3D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= Whether to adapt the density dynamically to reduce areas already covered. The density is still statistically followed for undersampled acquisitions. Default to ``True``. + **sampling_kwargs + Sampling parameters in + ``mrinufft.trajectories.sampling.sample_from_density`` used for the + shot starting positions. Returns ------- @@ -169,5 +196,10 @@ def initialize_3D_random_walk(Nc, Ns, density, *, diagonals=True, pseudo_random= if len(density.shape) != 3: raise ValueError("`density` is expected to be 3-dimensional.") return _initialize_ND_random_walk( - Nc, Ns, density, diagonals=diagonals, pseudo_random=pseudo_random + Nc, + Ns, + density, + diagonals=diagonals, + pseudo_random=pseudo_random, + **sampling_kwargs, ) diff --git a/src/mrinufft/trajectories/sampling.py b/src/mrinufft/trajectories/sampling.py index 2e1378fb..ea730260 100644 --- a/src/mrinufft/trajectories/sampling.py +++ b/src/mrinufft/trajectories/sampling.py @@ -12,7 +12,7 @@ def sample_from_density( - nb_samples, density, method="random", *, binary_mask_density=False + nb_samples, density, method="random", *, dim_compensation="auto" ): """ Sample points based on a given density distribution. @@ -22,18 +22,20 @@ def sample_from_density( nb_samples : int The number of samples to draw. density : np.ndarray - An array representing the density distribution from which samples are drawn. + An array representing the density distribution from which samples are drawn, + normalized automatically by its sum during the call for convenience. method : str, optional The sampling method to use, either 'random' for random sampling over the discrete grid defined by the density or 'lloyd' for Lloyd's algorithm over a continuous space, by default "random". - binary_mask_density : bool, optional - Whether to apply a specific density correction introduced + dim_compensation : str, bool, optional + Whether to apply a specific dimensionality compensation introduced in [Cha+14]_. An exponent ``N/(N-1)`` with ``N`` the number of dimensions in ``density`` is applied to fix the observed - density expectation when paths are traced between drawn samples - specifically over binary masks/grids. Otherwise it is demonstrated - that the empirical density will be ``density ** (N - 1) / N``. + density expectation when set to ``"auto"`` and ``method="lloyd"``. + It is also relevant to set it to ``True`` when ``method="random"`` + and one wants to create binary masks with continuous paths between + drawn samples. Returns ------- @@ -46,7 +48,6 @@ def sample_from_density( If ``nb_samples`` exceeds the total size of the density array or if the specified ``method`` is unknown. - References ---------- .. [Cha+14] Chauffert, Nicolas, Philippe Ciuciu, @@ -60,13 +61,19 @@ def sample_from_density( shape = np.array(density.shape) nb_dims = len(shape) max_nb_samples = np.prod(shape) + density = density / np.sum(density) if nb_samples > max_nb_samples: raise ValueError("`nb_samples` must be lower than the size of `density`.") - # Adapt density to use binary mask use case - density = density ** (nb_dims / (nb_dims - 1)) if binary_mask_density else density - density = density / np.sum(density) + # Check for dimensionality compensation + if isinstance(dim_compensation, str) and dim_compensation != "auto": + raise ValueError(f"Unknown string {dim_compensation} for `dim_compensation`.") + if (dim_compensation == "auto" and method == "lloyd") or ( + isinstance(dim_compensation, bool) and dim_compensation + ): + density = density ** (nb_dims / (nb_dims - 1)) + density = density / np.sum(density) # Sample using specified method if method == "random": From fcf699b235996dc692581a5cd22da5b41dd4802f Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Sat, 23 Nov 2024 23:47:56 +0100 Subject: [PATCH 08/18] Add utils to examples to show densities --- examples/example_2D_trajectories.py | 44 +++++------ examples/example_3D_trajectories.py | 72 ++++++++--------- examples/example_trajectory_tools.py | 76 +++++++++--------- examples/utils.py | 112 ++++++++++++++++++++++----- 4 files changed, 189 insertions(+), 115 deletions(-) diff --git a/examples/example_2D_trajectories.py b/examples/example_2D_trajectories.py index ed50a9f7..eae7f9e1 100644 --- a/examples/example_2D_trajectories.py +++ b/examples/example_2D_trajectories.py @@ -17,7 +17,7 @@ # External import matplotlib.pyplot as plt import numpy as np -from utils import show_argument, show_trajectory +from utils import show_trajectories, show_trajectory # Internal import mrinufft as mn @@ -75,7 +75,7 @@ arguments = [8, 16, 32, 64] function = lambda x: mn.initialize_2D_radial(x, Ns, tilt=tilt, in_out=in_out) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -89,7 +89,7 @@ arguments = [8, 16, 32, 64] function = lambda x: mn.initialize_2D_radial(Nc, x, tilt=tilt, in_out=in_out) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -105,7 +105,7 @@ arguments = ["uniform", "golden", "mri-golden", np.pi / 17] function = lambda x: mn.initialize_2D_radial(Nc, Ns, tilt=x, in_out=in_out) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -128,7 +128,7 @@ arguments = [True, False] function = lambda x: mn.initialize_2D_radial(Nc, Ns, tilt=tilt, in_out=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -170,7 +170,7 @@ function = lambda x: mn.initialize_2D_spiral( Nc, Ns, tilt=tilt, nb_revolutions=x, in_out=in_out ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -194,7 +194,7 @@ arguments = ["galilean", "archimedes", "fermat", 1 / 4] function = lambda x: mn.initialize_2D_spiral(Nc, Ns, tilt=tilt, spiral=x, in_out=in_out) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -221,7 +221,7 @@ Ns, patch_center=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -267,7 +267,7 @@ Ns, spiral_reduction=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -293,7 +293,7 @@ Ns, patch_center=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -333,7 +333,7 @@ function = lambda x: mn.initialize_2D_cones( Nc, Ns, tilt=tilt, in_out=in_out, nb_zigzags=x ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -348,7 +348,7 @@ arguments = [0.2, 1, 2, 3] function = lambda x: mn.initialize_2D_cones(Nc, Ns, tilt=tilt, in_out=in_out, width=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -413,7 +413,7 @@ arguments = [2, 3, 4, 6] function = lambda x: mn.initialize_2D_propeller(Nc, Ns, nb_strips=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -450,7 +450,7 @@ arguments = [Nc, int(2 * Nc / 3), int(Nc / 3)] function = lambda x: mn.initialize_2D_rings(Nc=x, Ns=Ns, nb_rings=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% # @@ -461,7 +461,7 @@ arguments = [Nc, int(4 * Nc / 3), 2 * Nc] function = lambda x: mn.initialize_2D_rings(Nc=x, Ns=Ns, nb_rings=Nc) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -501,7 +501,7 @@ arguments = [0, 1, 5, 10] function = lambda x: mn.initialize_2D_rosette(Nc, Ns, in_out=in_out, coprime_index=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -541,7 +541,7 @@ function = lambda x: mn.initialize_2D_polar_lissajous( Nc, Ns, in_out=in_out, coprime_index=x ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -564,7 +564,7 @@ function = lambda x: mn.initialize_2D_polar_lissajous( Nc, Ns, in_out=in_out, nb_segments=x ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -588,7 +588,7 @@ function = lambda x: mn.initialize_2D_polar_lissajous( Nc, Ns, in_out=io, coprime_index=cpi, nb_segments=x ) - show_argument( + show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size ) @@ -625,7 +625,7 @@ arguments = [1, 2.5, 5, 10] function = lambda x: mn.initialize_2D_waves(Nc, Ns, nb_zigzags=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -644,7 +644,7 @@ arguments = [0, 1, 1.5, 3] function = lambda x: mn.initialize_2D_waves(Nc, Ns, width=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -679,7 +679,7 @@ arguments = [1, 1.5, 2, 3] function = lambda x: mn.initialize_2D_lissajous(Nc, Ns, density=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% diff --git a/examples/example_3D_trajectories.py b/examples/example_3D_trajectories.py index 21bab047..5224f343 100644 --- a/examples/example_3D_trajectories.py +++ b/examples/example_3D_trajectories.py @@ -26,7 +26,7 @@ # External import matplotlib.pyplot as plt import numpy as np -from utils import show_argument, show_trajectory +from utils import show_trajectories, show_trajectory # Internal import mrinufft as mn @@ -93,7 +93,7 @@ arguments = [Nc // 4, Nc // 2, Nc, Nc * 2] function = lambda x: mn.initialize_3D_phyllotaxis_radial(x, Ns, in_out=in_out) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -106,7 +106,7 @@ arguments = [10, 25, 40, 100] function = lambda x: mn.initialize_3D_phyllotaxis_radial(Nc, x, in_out=in_out) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -129,7 +129,7 @@ arguments = [True, False] function = lambda x: mn.initialize_3D_phyllotaxis_radial(Nc, Ns, in_out=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -265,7 +265,7 @@ arguments = ["uniform", "golden", "mri-golden", np.pi / 17] function = lambda x: mn.initialize_3D_cones(Nc, Ns, tilt=x, in_out=in_out) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -278,7 +278,7 @@ arguments = [0.5, 2, 5, 10] function = lambda x: mn.initialize_3D_cones(Nc, Ns, in_out=in_out, nb_zigzags=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -293,7 +293,7 @@ arguments = ["archimedes", "fermat", 0.5, 1.5] function = lambda x: mn.initialize_3D_cones(Nc, Ns, in_out=in_out, spiral=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -308,7 +308,7 @@ arguments = [0.2, 1, 2, 3] function = lambda x: mn.initialize_3D_cones(Nc, Ns, in_out=in_out, width=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -372,11 +372,11 @@ max_angle=np.pi / 4, axes=x, )[::-1] -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -418,7 +418,7 @@ arguments = [0.5, 2.5, 5, 10] function = lambda x: mn.initialize_3D_wave_caipi(Nc, Ns, nb_revolutions=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% # ``width (float)`` @@ -434,7 +434,7 @@ arguments = [0.2, 1, 2, 3] function = lambda x: mn.initialize_3D_wave_caipi(Nc, Ns, width=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% # ``packing (str)`` @@ -456,11 +456,11 @@ arguments = ["triangular", "square", "circular", "fibonacci", "random"] function = lambda x: mn.initialize_3D_wave_caipi(Nc, Ns, packing=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -482,11 +482,11 @@ arguments = ["circle", "square", "diamond", 0.5] function = lambda x: mn.initialize_3D_wave_caipi(Nc, Ns, shape=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -503,11 +503,11 @@ arguments = [(1, 1), (2, 1), (1, 2), (2.3, 1.8)] function = lambda x: mn.initialize_3D_wave_caipi(Nc, Ns, packing="square", spacing=x) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -561,7 +561,7 @@ function = lambda x: mn.initialize_3D_seiffert_spiral( Nc, Ns, in_out=in_out, curve_index=x ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -581,7 +581,7 @@ in_out=in_out, nb_revolutions=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -605,7 +605,7 @@ axis_tilt=x, spiral_tilt=0, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -631,7 +631,7 @@ axis_tilt="golden", spiral_tilt=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -680,7 +680,7 @@ function = lambda x: mn.initialize_3D_helical_shells( Nc=x, Ns=Ns, nb_shells=x, spiral_reduction=2 ) -show_argument(function, arguments, one_shot=False, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=False, subfig_size=subfigure_size) # %% @@ -699,7 +699,7 @@ function = lambda x: mn.initialize_3D_helical_shells( Nc=Nc, Ns=Ns, nb_shells=nb_repetitions, spiral_reduction=x ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -713,7 +713,7 @@ function = lambda x: mn.initialize_3D_helical_shells( Nc=Nc, Ns=Ns, nb_shells=nb_repetitions, spiral_reduction=2, shell_tilt=x ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -730,7 +730,7 @@ function = lambda x: mn.initialize_3D_helical_shells( Nc=Nc, Ns=Ns, nb_shells=nb_repetitions, spiral_reduction=2, shot_tilt=x ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -782,7 +782,7 @@ function = lambda x: mn.initialize_3D_annular_shells( Nc, Ns, nb_shells=nb_repetitions, ring_tilt=x ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -878,7 +878,7 @@ Ns_transitions=x, nb_blades=nb_blades, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -897,7 +897,7 @@ Ns_transitions=Ns // 10, nb_blades=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -915,11 +915,11 @@ nb_blades=nb_blades, blade_tilt=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -941,7 +941,7 @@ nb_blades=nb_blades, nb_trains=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -965,11 +965,11 @@ nb_blades=nb_blades, skip_factor=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -1054,7 +1054,7 @@ nb_blade_revolutions=x, nb_spiral_revolutions=0, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -1071,7 +1071,7 @@ nb_blade_revolutions=x, nb_spiral_revolutions=nb_revolutions, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% diff --git a/examples/example_trajectory_tools.py b/examples/example_trajectory_tools.py index e0a2efff..83afeb5d 100644 --- a/examples/example_trajectory_tools.py +++ b/examples/example_trajectory_tools.py @@ -24,7 +24,7 @@ # External import matplotlib.pyplot as plt import numpy as np -from utils import show_argument, show_trajectory +from utils import show_trajectories, show_trajectory # Internal import mrinufft as mn @@ -87,7 +87,9 @@ arguments = ["Radial", "Spiral", "2D Cones", "3D Cones"] function = lambda x: single_trajectories[x] -show_argument(function, arguments, one_shot=bool(one_shot), subfig_size=subfigure_size) +show_trajectories( + function, arguments, one_shot=bool(one_shot), subfig_size=subfigure_size +) # %% @@ -119,7 +121,9 @@ arguments = ["Radial", "Spiral", "2D Cones", "3D Cones"] function = lambda x: planar_trajectories[x] -show_argument(function, arguments, one_shot=bool(one_shot), subfig_size=subfigure_size) +show_trajectories( + function, arguments, one_shot=bool(one_shot), subfig_size=subfigure_size +) # %% @@ -160,9 +164,9 @@ arguments = ["Radial", "Spiral", "2D Cones", "3D Cones"] function = lambda x: tools.stack(planar_trajectories[x], nb_stacks=nb_repetitions) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -188,7 +192,7 @@ ), nb_stacks=nb_repetitions, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -206,7 +210,7 @@ trajectory = np.copy(planar_trajectories["3D Cones"]) trajectory[..., 2] *= 2 function = lambda x: tools.stack(trajectory, nb_stacks=nb_repetitions, hard_bounded=x) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -258,9 +262,9 @@ nb_rotations=nb_repetitions, x_tilt="uniform", ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -326,7 +330,7 @@ half_sphere=in_out, axis=2, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% # @@ -345,7 +349,7 @@ half_sphere=in_out, axis=0, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -369,7 +373,7 @@ half_sphere=x, axis=0, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -396,7 +400,7 @@ partition=x, axis=0, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -422,7 +426,7 @@ partition=x, axis=0, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -445,7 +449,7 @@ partition=x, axis=0, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% # ``axis (int, array)`` @@ -478,7 +482,7 @@ half_sphere=in_out, axis=x, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -501,7 +505,7 @@ half_sphere=in_out, axis=x, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% # @@ -556,9 +560,9 @@ function = lambda x: tools.conify( planar_trajectories[x], nb_cones=nb_repetitions, in_out=in_out ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -578,9 +582,9 @@ function = lambda x: tools.conify( single_trajectories[x], nb_cones=Nc, z_tilt="golden", in_out=in_out ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -607,7 +611,7 @@ in_out=in_out, max_angle=x, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -634,7 +638,7 @@ max_angle=np.pi / 2, borderless=x, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -691,11 +695,11 @@ nb_trains=Nc_planes // 2, reverse_odd_shots=True, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -716,7 +720,7 @@ nb_trains=Nc_planes // 2, reverse_odd_shots=True, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -735,7 +739,7 @@ nb_trains=x, reverse_odd_shots=True, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -755,7 +759,7 @@ nb_trains=Nc_planes // 2, reverse_odd_shots=x, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -799,11 +803,11 @@ tools.rewind(planar_trajectories[x], Ns_transitions=Ns // 10), Ns_transitions=Ns // 10, ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -821,7 +825,7 @@ tools.rewind(planar_trajectories["2D Cones"], Ns_transitions=x), Ns_transitions=x, ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, subfig_size=subfigure_size, dim="2D" ) @@ -898,9 +902,9 @@ function = lambda x: tools.stack_spherically( init_trajectories[x], Nc=Nc, nb_stacks=nb_repetitions ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, @@ -935,7 +939,7 @@ trajectories = {"Classic": traj_classic, "Normalized": traj_normal} arguments = ["Classic", "Normalized"] function = lambda x: trajectories[x] -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% @@ -988,7 +992,7 @@ function = lambda x: tools.shellify( init_trajectories[x], Nc=Nc, nb_shells=nb_repetitions ) -show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) # %% # ``hemisphere_mode (str)`` @@ -1005,7 +1009,7 @@ function = lambda x: tools.shellify( init_trajectories["Spiral"], Nc=Nc, nb_shells=nb_repetitions, hemisphere_mode=x ) -show_argument( +show_trajectories( function, arguments, one_shot=one_shot, diff --git a/examples/utils.py b/examples/utils.py index 7c70595b..5edaedc5 100644 --- a/examples/utils.py +++ b/examples/utils.py @@ -3,16 +3,39 @@ the examples. """ -import matplotlib.pyplot as plt - # External imports import numpy as np +from matplotlib import colors +import matplotlib.pyplot as plt # Internal imports from mrinufft import display_2D_trajectory, display_3D_trajectory, displayConfig +from mrinufft.trajectories.utils import KMAX + + +def show_trajectory(trajectory, one_shot, figure_size): + if trajectory.shape[-1] == 2: + ax = display_2D_trajectory( + trajectory, size=figure_size, one_shot=one_shot % trajectory.shape[0] + ) + ax.set_aspect("equal") + plt.tight_layout() + plt.show() + else: + ax = display_3D_trajectory( + trajectory, + size=figure_size, + one_shot=one_shot % trajectory.shape[0], + per_plane=False, + ) + plt.tight_layout() + plt.subplots_adjust(bottom=0.1) + plt.show() -def show_argument(function, arguments, one_shot, subfig_size, dim="3D", axes=(0, 1)): +def show_trajectories( + function, arguments, one_shot, subfig_size, dim="3D", axes=(0, 1) +): # Initialize trajectories with varying option trajectories = [function(arg) for arg in arguments] @@ -42,25 +65,72 @@ def show_argument(function, arguments, one_shot, subfig_size, dim="3D", axes=(0, ax.set_xlabel(labels[axes[0]], fontsize=displayConfig.fontsize) ax.set_ylabel(labels[axes[1]], fontsize=displayConfig.fontsize) ax.set_aspect("equal") - ax.set_title(str(arg), fontsize=4 * subfig_size) + ax.set_title(str(arg), fontsize=displayConfig.fontsize) plt.show() -def show_trajectory(trajectory, one_shot, figure_size): - if trajectory.shape[-1] == 2: - ax = display_2D_trajectory( - trajectory, size=figure_size, one_shot=one_shot % trajectory.shape[0] - ) - ax.set_aspect("equal") - plt.tight_layout() - plt.show() +def show_density(density, figure_size, *, log_scale=False): + density = density.T[::-1] + + plt.figure(figsize=(figure_size, figure_size)) + if log_scale: + plt.imshow(density, cmap="jet", norm=colors.LogNorm()) else: - ax = display_3D_trajectory( - trajectory, - size=figure_size, - one_shot=one_shot % trajectory.shape[0], - per_plane=False, - ) - plt.tight_layout() - plt.subplots_adjust(bottom=0.1) - plt.show() + plt.imshow(density, cmap="jet") + + ax = plt.gca() + ax.set_xlabel("kx", fontsize=displayConfig.fontsize) + ax.set_ylabel("ky", fontsize=displayConfig.fontsize) + ax.set_aspect("equal") + + plt.axis(False) + plt.colorbar() + plt.show() + + +def show_densities(function, arguments, subfig_size, *, log_scale=False): + # Initialize k-space densities with varying option + densities = [function(arg).T[::-1] for arg in arguments] + + # Plot the trajectories side by side + fig, axes = plt.subplots( + 1, + len(densities), + figsize=(len(densities) * subfig_size, subfig_size), + constrained_layout=True, + ) + + for ax, arg, density in zip(axes, arguments, densities): + ax.set_title(str(arg), fontsize=displayConfig.fontsize) + ax.set_xlabel("kx", fontsize=displayConfig.fontsize) + ax.set_ylabel("ky", fontsize=displayConfig.fontsize) + ax.set_aspect("equal") + if log_scale: + ax.imshow(density, cmap="jet", norm=colors.LogNorm()) + else: + ax.imshow(density, cmap="jet") + ax.axis(False) + plt.show() + + +def show_locations(function, arguments, subfig_size, *, log_scale=False): + # Initialize k-space locations with varying option + locations = [function(arg) for arg in arguments] + + # Plot the trajectories side by side + fig, axes = plt.subplots( + 1, + len(locations), + figsize=(len(locations) * subfig_size, subfig_size), + constrained_layout=True, + ) + + for ax, arg, location in zip(axes, arguments, locations): + ax.set_title(str(arg), fontsize=displayConfig.fontsize) + ax.set_xlim(-1.05 * KMAX, 1.05 * KMAX) + ax.set_ylim(-1.05 * KMAX, 1.05 * KMAX) + ax.set_xlabel("kx", fontsize=displayConfig.fontsize) + ax.set_ylabel("ky", fontsize=displayConfig.fontsize) + ax.set_aspect("equal") + ax.scatter(location[..., 0], location[..., 1], s=3) + plt.show() From fdcda5531949cec4352c2aa57634a5deb0bd1799 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Sun, 24 Nov 2024 22:33:42 +0100 Subject: [PATCH 09/18] Add sampling density example skeleton --- examples/example_sampling_densities.py | 387 ++++++++++++++++++ src/mrinufft/trajectories/maths/tsp_solver.py | 10 +- src/mrinufft/trajectories/sampling.py | 4 +- 3 files changed, 394 insertions(+), 7 deletions(-) create mode 100644 examples/example_sampling_densities.py diff --git a/examples/example_sampling_densities.py b/examples/example_sampling_densities.py new file mode 100644 index 00000000..c07fb32b --- /dev/null +++ b/examples/example_sampling_densities.py @@ -0,0 +1,387 @@ +""" +================== +Sampling densities +================== + +A collection of sampling densities and density-based non-Cartesian trajectories. + +""" + +# %% +# In this example we illustrate the use of different sampling densities, +# and show how to generate trajectories based on them, such as random +# walks and travelling-salesman trajectories. +# + +# External +import brainweb_dl as bwdl +import matplotlib.pyplot as plt +import numpy as np +from utils import ( + show_density, + show_densities, + show_locations, + show_trajectory, + show_trajectories, +) + +# Internal +import mrinufft as mn +from mrinufft import display_2D_trajectory, display_3D_trajectory + + +# %% +# Script options +# ============== +# +# These options are used in the examples below as default values for +# all densities and trajectories. + +# Density parameters +shape_2d = (100, 100) +shape_3d = (100, 100, 100) + +# Trajectory parameters +Nc = 10 # Number of shots +Ns = 50 # Number of samples per shot + +# Display parameters +figure_size = 5.5 # Figure size for trajectory plots +subfigure_size = 3 # Figure size for subplots +one_shot = 0 # Highlight one shot in particular + +# %% +# Densities +# ========= +# +# In this section we present different sampling densities +# with various properties. +# +# Cutoff/decay density +# -------------------- +# +# Create a density composed of a central constant-value ellipsis +# defined by a cutoff ratio, followed by a polynomial decay over +# outer regions as defined in [Cha+22]_. + +cutoff_density = mn.create_cutoff_decay_density(shape=shape_2d, cutoff=0.2, decay=2) +show_density(cutoff_density, figure_size=figure_size) + +# %% +# ``cutoff (float)`` +# ~~~~~~~~~~~~~~~~~~ +# +# The k-space radius cutoff ratio between 0 and 1 within +# which density remains uniform and beyond which it decays. +# It is modulated by ``resolution`` to create ellipsoids. +# +# The ``mrinufft.trajectories.sampling.create_polynomial_density`` +# simply calls this function with ``cutoff=0``. + +arguments = [0, 0.1, 0.2, 0.3] +function = lambda x: mn.create_cutoff_decay_density( + shape=shape_2d, + cutoff=x, + decay=2, +) +show_densities( + function, + arguments, + subfig_size=subfigure_size, +) + + +# %% +# ``decay (float)`` +# ~~~~~~~~~~~~~~~~~ +# +# The polynomial decay in density beyond the cutoff ratio. +# It can be zero or negative as shown below, but most applications +# are expected have decays in the positive range. + + +arguments = [-1, 0, 0.5, 2] +function = lambda x: mn.create_cutoff_decay_density( + shape=shape_2d, + cutoff=0.2, + decay=x, +) +show_densities( + function, + arguments, + subfig_size=subfigure_size, +) + +# %% +# ``resolution (tuple)`` +# ~~~~~~~~~~~~~~~~~~~~~~ +# +# Resolution scaling factors for each dimension of the density grid, +# by default ``None``. Note on the example below that the unit doesn't +# matter because ``cutoff`` is a ratio and ``decay`` is an exponent, +# so only the relative factor between the dimensions is important. +# +# This argument can be used to handle anisotropy but also to produce +# ellipsoidal densities. + + +arguments = [None, (1, 1), (1, 2), (1e-3, 0.5e-3)] +function = lambda x: mn.create_cutoff_decay_density( + shape=shape_2d, + cutoff=0.2, + decay=2, + resolution=x, +) +show_densities( + function, + arguments, + subfig_size=subfigure_size, +) + +# %% +# Energy-based density +# -------------------- +# +# A common intuition is to consider that the sampling density +# should be proportional to the k-space amplitude. It can be +# learned from existing datasets and used for new acquisitions. + +dataset = bwdl.get_mri(4, "T1")[:, ::2, ::2] +energy_density = mn.create_energy_density(dataset=dataset) +show_density(energy_density, figure_size=figure_size, log_scale=True) + +# %% +# ``dataset (np.ndarray)`` +# ~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The dataset from which to calculate the density +# based on its Fourier transform, with an expected +# shape (nb_volumes, dim_1, ..., dim_N). +# An N-dimensional Fourier transform is then performed. +# +# In the example below, we show the resulting densities +# from different slices of a single volume for convenience. +# More relevant use cases would be to learn densities for +# different organs and/or contrasts. + + +arguments = [50, 100, 150] +function = lambda x: mn.create_energy_density(dataset=bwdl.get_mri(4, "T1")[x : x + 20]) +show_densities( + function, + arguments, + subfig_size=subfigure_size, + log_scale=True, +) + + +# %% +# Chauffert's density +# ------------------- +# + + +chauffert_density = mn.create_chauffert_density( + shape=shape_2d, + wavelet_basis="haar", + nb_wavelet_scales=3, +) +show_density(chauffert_density, figure_size=figure_size) + + +# %% +# ``wavelet_basis (str)`` +# ~~~~~~~~~~~~~~~~~~~~~~~ + +arguments = ["haar", "rbio2.2", "coif4", "sym8"] +function = lambda x: mn.create_chauffert_density( + shape=shape_2d, + wavelet_basis=x, + nb_wavelet_scales=3, +) +show_densities( + function, + arguments, + subfig_size=subfigure_size, +) + +# %% +# ``nb_wavelet_scales (int)`` +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +arguments = [1, 2, 3, 4] +function = lambda x: mn.create_chauffert_density( + shape=shape_2d, + wavelet_basis="haar", + nb_wavelet_scales=x, +) +show_densities( + function, + arguments, + subfig_size=subfigure_size, +) + + +# %% +# Custom density +# -------------- + +# Linear gradient +density = np.tile(np.linspace(0, 1, shape_2d[1])[:, None], (1, shape_2d[0])) +# Square center +density[ + 3 * shape_2d[0] // 8 : 5 * shape_2d[0] // 8, + 3 * shape_2d[1] // 8 : 5 * shape_2d[1] // 8, +] = 2 +# Outer circle +density = np.where( + np.linalg.norm(np.indices(shape_2d) - shape_2d[0] / 2, axis=0) < shape_2d[0] / 2, + density, + 0, +) +# Normalization +custom_density = density / np.sum(density) + +show_density(custom_density, figure_size=figure_size) + + +# %% +# Sampling +# ======== +# +# In this section we present random, pseudo-random and +# algorithm-based sampling methods. The examples are based +# on a few densities picked from the ones presented above. +# + +densities = { + "Cutoff/Decay": cutoff_density, + "Energy": energy_density, + "Chauffert": chauffert_density, + "Custom": custom_density, +} + +arguments = densities.keys() +function = lambda x: densities[x] +show_densities(function, arguments, subfig_size=subfigure_size) + + +# %% +# Random sampling +# --------------- + +arguments = densities.keys() +function = lambda x: mn.sample_from_density(Nc * Ns, densities[x], method="random") +show_locations(function, arguments, subfig_size=subfigure_size) + + +# %% +# Lloyd's sampling +# ---------------- + +arguments = densities.keys() +function = lambda x: mn.sample_from_density(Nc * Ns, densities[x], method="lloyd") +show_locations(function, arguments, subfig_size=subfigure_size) + + +# %% +# Density-based trajectories +# ========================== +# +# In this section we present 2D and 3D trajectories based +# on arbitrary densities, and also sampling for some of them. +# +# Random walks +# ------------ +# + +arguments = densities.keys() +function = lambda x: mn.initialize_2D_random_walk( + Nc, Ns, density=densities[x][::4, ::4] +) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) + +# %% +# + +arguments = densities.keys() +function = lambda x: mn.initialize_2D_random_walk( + Nc, Ns, density=densities[x][::4, ::4], method="lloyd" +) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) + +# %% +# +# Oversampled + +arguments = densities.keys() +function = lambda x: mn.oversample( + mn.initialize_2D_random_walk( + Nc, Ns, density=densities[x][::4, ::4], method="lloyd" + ), + 4 * Ns, +) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) + + +# %% +# Travelling Salesman +# ------------------- +# + +arguments = densities.keys() +function = lambda x: mn.initialize_2D_travelling_salesman( + Nc, + Ns, + density=densities[x], +) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) + +# %% +# + +arguments = densities.keys() +function = lambda x: mn.initialize_2D_travelling_salesman( + Nc, + Ns, + density=densities[x], + method="lloyd", +) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) + +# %% +# +# Oversampled + +arguments = densities.keys() +function = lambda x: mn.oversample( + mn.initialize_2D_travelling_salesman(Nc, Ns, density=densities[x], method="lloyd"), + 4 * Ns, +) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) + +# %% +# + +arguments = ((None, None, None), ("y", None, "x"), ("phi", None, "r"), ("y", "x", "r")) +function = lambda x: mn.initialize_2D_travelling_salesman( + Nc, + Ns, + density=densities["Custom"], + first_cluster_by=x[0], + second_cluster_by=x[1], + sort_by=x[2], + method="lloyd", +) +show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) + + +# %% +# References +# ========== +# +# .. [Cha+22] Chaithya, G. R., Pierre Weiss, Guillaume Daval-Frérot, +# Aurélien Massire, Alexandre Vignaud, and Philippe Ciuciu. +# "Optimizing full 3d sparkling trajectories for high-resolution +# magnetic resonance imaging." +# IEEE Transactions on Medical Imaging 41, no. 8 (2022): 2105-2117. diff --git a/src/mrinufft/trajectories/maths/tsp_solver.py b/src/mrinufft/trajectories/maths/tsp_solver.py index 9b254eb8..b4bcf859 100644 --- a/src/mrinufft/trajectories/maths/tsp_solver.py +++ b/src/mrinufft/trajectories/maths/tsp_solver.py @@ -6,11 +6,11 @@ def solve_tsp_with_2opt(locations, improvement_threshold=1e-8): """Solve the TSP problem using a 2-opt approach. - A sub-optimal solution to the Travelling Salesman Problem (TSP) - is provided using the 2-opt approach in O(n²) where chunks of - an initially random route are reversed, and selected if the - total distance is reduced. As a result the route solution - does not cross its own path in 2D. + A sub-optimal solution to the Travelling Salesman Problem (TSP) + is provided using the 2-opt approach in O(n²) where chunks of + an initially random route are reversed, and selected if the + total distance is reduced. As a result the route solution + does not cross its own path in 2D. Parameters ---------- diff --git a/src/mrinufft/trajectories/sampling.py b/src/mrinufft/trajectories/sampling.py index ea730260..d58b3539 100644 --- a/src/mrinufft/trajectories/sampling.py +++ b/src/mrinufft/trajectories/sampling.py @@ -120,10 +120,10 @@ def create_cutoff_decay_density(shape, cutoff, decay, resolution=None): The k-space radius cutoff ratio between 0 and 1 within which density remains uniform and beyond which it decays. decay : float - The rate of decay of density beyond the cutoff ratio. + The polynomial decay in density beyond the cutoff ratio. resolution : np.ndarray, optional Resolution scaling factors for each dimension of the density grid, - by default None. + by default ``None``. Returns ------- From 2f631cacb8100e6da1d40adc4250363fd7b22a34 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Mon, 25 Nov 2024 23:44:48 +0100 Subject: [PATCH 10/18] Add descriptions to sampling density examples --- examples/example_sampling_densities.py | 79 ++++++++++++++++--- .../trajectories/inits/travelling_salesman.py | 10 ++- src/mrinufft/trajectories/sampling.py | 2 +- 3 files changed, 77 insertions(+), 14 deletions(-) diff --git a/examples/example_sampling_densities.py b/examples/example_sampling_densities.py index c07fb32b..cb6da78d 100644 --- a/examples/example_sampling_densities.py +++ b/examples/example_sampling_densities.py @@ -75,7 +75,7 @@ # which density remains uniform and beyond which it decays. # It is modulated by ``resolution`` to create ellipsoids. # -# The ``mrinufft.trajectories.sampling.create_polynomial_density`` +# The ``mrinufft.create_polynomial_density`` # simply calls this function with ``cutoff=0``. arguments = [0, 0.1, 0.2, 0.3] @@ -90,7 +90,6 @@ subfig_size=subfigure_size, ) - # %% # ``decay (float)`` # ~~~~~~~~~~~~~~~~~ @@ -99,7 +98,6 @@ # It can be zero or negative as shown below, but most applications # are expected have decays in the positive range. - arguments = [-1, 0, 0.5, 2] function = lambda x: mn.create_cutoff_decay_density( shape=shape_2d, @@ -138,6 +136,7 @@ subfig_size=subfigure_size, ) + # %% # Energy-based density # -------------------- @@ -164,7 +163,6 @@ # More relevant use cases would be to learn densities for # different organs and/or contrasts. - arguments = [50, 100, 150] function = lambda x: mn.create_energy_density(dataset=bwdl.get_mri(4, "T1")[x : x + 20]) show_densities( @@ -179,7 +177,14 @@ # Chauffert's density # ------------------- # - +# This is a reproduction of the proposition from [CCW13]_. +# A sampling density is derived from compressed sensing +# equations to maximize guarantees of exact image recovery +# for a specified sparse wavelet domain decomposition. +# +# This principle is valid for any linear transform but +# for convenience it was limited to wavelets as in the +# original implementation. chauffert_density = mn.create_chauffert_density( shape=shape_2d, @@ -192,6 +197,10 @@ # %% # ``wavelet_basis (str)`` # ~~~~~~~~~~~~~~~~~~~~~~~ +# +# The wavelet basis to use for wavelet decomposition, either +# as a built-in wavelet name from the PyWavelets package +# or as a custom ``pywt.Wavelet`` object. arguments = ["haar", "rbio2.2", "coif4", "sym8"] function = lambda x: mn.create_chauffert_density( @@ -208,6 +217,8 @@ # %% # ``nb_wavelet_scales (int)`` # ~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The number of wavelet scales to use in decomposition. arguments = [1, 2, 3, 4] function = lambda x: mn.create_chauffert_density( @@ -225,6 +236,9 @@ # %% # Custom density # -------------- +# +# Any density can be defined and later used for sampling and +# trajectories. # Linear gradient density = np.tile(np.linspace(0, 1, shape_2d[1])[:, None], (1, shape_2d[0])) @@ -269,6 +283,9 @@ # %% # Random sampling # --------------- +# +# This sampling simply consists of weighted-random selection from the +# density grid locations. arguments = densities.keys() function = lambda x: mn.sample_from_density(Nc * Ns, densities[x], method="random") @@ -278,6 +295,11 @@ # %% # Lloyd's sampling # ---------------- +# +# This sampling is based on a Voronoi/Dirichlet tesselation using Lloyd's +# weighted KMeans algorithm. The implementation is based on +# ``sklearn.cluster.KMeans`` in 2D and ``sklearn.cluster.BisectingKMeans`` +# in 3D, mostly to reduce computation times in the most demanding cases. arguments = densities.keys() function = lambda x: mn.sample_from_density(Nc * Ns, densities[x], method="lloyd") @@ -293,7 +315,15 @@ # # Random walks # ------------ +# +# This is an adaptation of the proposition from [Cha+14]_. +# It creates a trajectory by walking randomly to neighboring points +# following a provided sampling density. # +# This implementation is different from the original proposition: +# trajectories are continuous with a fixed length instead of +# making random jumps to other locations, and an option +# is provided to have pseudo-random walks to improve coverage. arguments = densities.keys() function = lambda x: mn.initialize_2D_random_walk( @@ -303,6 +333,11 @@ # %% # +# The starting shot positions can be modified to follow Lloyd's sampling +# method rather than the default random approach, resulting in more evenly +# spaced shots that still respect the prescribed density. +# Additional ``kwargs`` can provided to set the arguments in +# ``mrinufft.sample_from_density``. arguments = densities.keys() function = lambda x: mn.initialize_2D_random_walk( @@ -310,9 +345,10 @@ ) show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) -# %% +# %% # -# Oversampled +# The random paths can be made into a smooth and continuous +# trajectory by oversampling the shots with cubic splines. arguments = densities.keys() function = lambda x: mn.oversample( @@ -328,6 +364,9 @@ # Travelling Salesman # ------------------- # +# This is a reproduction of the work from [Cha+14]_. The Travelling +# Salesman Problem (TSP) solution is obtained using the 2-opt method +# with a complexity in O(n²) in time and memory. arguments = densities.keys() function = lambda x: mn.initialize_2D_travelling_salesman( @@ -339,6 +378,10 @@ # %% # +# It is possible to customize the sampling method using ``kwargs`` +# to provide arguments to ``mrinufft.sample_from_density``. +# For example, one can use Lloyd's sampling method to create evenly +# spaced point distributions and obtain a more deterministic coverage. arguments = densities.keys() function = lambda x: mn.initialize_2D_travelling_salesman( @@ -351,7 +394,10 @@ # %% # -# Oversampled +# Similarly to random walks, the travelling paths can be smoothed +# by oversampling the shots with cubic splines. Another use case +# is to reduce the number of TSP points to reduce the computation load +# and then oversample up to the desired shot length. arguments = densities.keys() function = lambda x: mn.oversample( @@ -362,6 +408,12 @@ # %% # +# An option is provided to cluster the points before calling the TSP solver, +# reducing drastically the computation time. +# Clusters are chosen by Cartesian (``"x"``, ``"y"``, ``"z"``) or spherical +# (``"r"``, ``"phi"``, ``"theta"``) coordinate with up to two coordinates. +# Then the points can be sorted within each cluster in order to define a general +# shot direction as shown below. arguments = ((None, None, None), ("y", None, "x"), ("phi", None, "r"), ("y", "x", "r")) function = lambda x: mn.initialize_2D_travelling_salesman( @@ -380,8 +432,17 @@ # References # ========== # +# .. [CCW13] Chauffert, Nicolas, Philippe Ciuciu, and Pierre Weiss. +# "Variable density compressed sensing in MRI. +# Theoretical vs heuristic sampling strategies." +# In 2013 IEEE 10th International Symposium on Biomedical Imaging, +# pp. 298-301. IEEE, 2013. +# .. [Cha+14] Chauffert, Nicolas, Philippe Ciuciu, +# Jonas Kahn, and Pierre Weiss. +# "Variable density sampling with continuous trajectories." +# SIAM Journal on Imaging Sciences 7, no. 4 (2014): 1962-1992. # .. [Cha+22] Chaithya, G. R., Pierre Weiss, Guillaume Daval-Frérot, # Aurélien Massire, Alexandre Vignaud, and Philippe Ciuciu. -# "Optimizing full 3d sparkling trajectories for high-resolution +# "Optimizing full 3D SPARKLING trajectories for high-resolution # magnetic resonance imaging." # IEEE Transactions on Medical Imaging 41, no. 8 (2022): 2105-2117. diff --git a/src/mrinufft/trajectories/inits/travelling_salesman.py b/src/mrinufft/trajectories/inits/travelling_salesman.py index 590ede8e..6b4f7764 100644 --- a/src/mrinufft/trajectories/inits/travelling_salesman.py +++ b/src/mrinufft/trajectories/inits/travelling_salesman.py @@ -226,10 +226,12 @@ def initialize_3D_travelling_salesman( Initialize a 3D trajectory using a Travelling Salesman Problem (TSP)-based path. This is a reproduction of the work from [Cha+14]_. The TSP solution - is obtained using the 2-opt method in O(n²). An additional option - is provided to cluster shots before solving the TSP and thus - reduce drastically the computation time. The initial sampling method - can also be customized. + is obtained using the 2-opt method with a complexity in O(n²) in time + and memory. + + An additional option is provided to cluster shots before solving the + TSP and thus reduce drastically the computation time. The initial + sampling method can also be customized. Parameters ---------- diff --git a/src/mrinufft/trajectories/sampling.py b/src/mrinufft/trajectories/sampling.py index d58b3539..0a3059a8 100644 --- a/src/mrinufft/trajectories/sampling.py +++ b/src/mrinufft/trajectories/sampling.py @@ -135,7 +135,7 @@ def create_cutoff_decay_density(shape, cutoff, decay, resolution=None): ---------- .. [Cha+22] Chaithya, G. R., Pierre Weiss, Guillaume Daval-Frérot, Aurélien Massire, Alexandre Vignaud, and Philippe Ciuciu. - "Optimizing full 3d sparkling trajectories for high-resolution + "Optimizing full 3D SPARKLING trajectories for high-resolution magnetic resonance imaging." IEEE Transactions on Medical Imaging 41, no. 8 (2022): 2105-2117. """ From 83b90815f84bb5a123444d780b870ab65066d84f Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Thu, 28 Nov 2024 23:40:33 +0100 Subject: [PATCH 11/18] Add dependency to pywavelet --- examples/example_sampling_densities.py | 4 ++-- pyproject.toml | 1 + src/mrinufft/trajectories/sampling.py | 18 +++++++++++++++++- 3 files changed, 20 insertions(+), 3 deletions(-) diff --git a/examples/example_sampling_densities.py b/examples/example_sampling_densities.py index cb6da78d..cc505a6f 100644 --- a/examples/example_sampling_densities.py +++ b/examples/example_sampling_densities.py @@ -315,7 +315,7 @@ # # Random walks # ------------ -# +# # This is an adaptation of the proposition from [Cha+14]_. # It creates a trajectory by walking randomly to neighboring points # following a provided sampling density. @@ -345,7 +345,7 @@ ) show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size) -# %% +# %% # # The random paths can be made into a smooth and continuous # trajectory by oversampling the shots with cubic splines. diff --git a/pyproject.toml b/pyproject.toml index bfe9bcd7..db2f0b32 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,6 +19,7 @@ pynfft = ["pynfft2>=1.4.3", "numpy>=2.0.0"] pynufft = ["pynufft"] io = ["pymapvbvd"] smaps = ["scikit-image"] +wavelets = ["pywavelets"] autodiff = ["torch"] diff --git a/src/mrinufft/trajectories/sampling.py b/src/mrinufft/trajectories/sampling.py index 0a3059a8..bca52cbb 100644 --- a/src/mrinufft/trajectories/sampling.py +++ b/src/mrinufft/trajectories/sampling.py @@ -4,12 +4,16 @@ import numpy.fft as nf import numpy.linalg as nl import numpy.random as nr -import pywt as pw from sklearn.cluster import BisectingKMeans, KMeans from tqdm.auto import tqdm from .utils import KMAX +try: + import pywt as pw +except ImportError: + pw = None + def sample_from_density( nb_samples, density, method="random", *, dim_compensation="auto" @@ -254,6 +258,12 @@ def create_chauffert_density(shape, wavelet_basis, nb_wavelet_scales, verbose=Fa In 2013 IEEE 10th International Symposium on Biomedical Imaging, pp. 298-301. IEEE, 2013. """ + if pw is None: + raise ImportError( + "The PyWavelets package must be installed " + "as an additional dependency for this function." + ) + nb_dims = len(shape) indices = np.indices(shape).reshape((nb_dims, -1)).T @@ -319,6 +329,12 @@ def create_fast_chauffert_density(shape, wavelet_basis, nb_wavelet_scales): In 2013 IEEE 10th International Symposium on Biomedical Imaging, pp. 298-301. IEEE, 2013. """ + if pw is None: + raise ImportError( + "The PyWavelets package must be installed " + "as an additional dependency for this function." + ) + nb_dims = len(shape) density = np.ones(shape) From 82172be11b92c30790d2169ca8c992205f72801b Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Mon, 2 Dec 2024 20:03:20 +0100 Subject: [PATCH 12/18] Improve cutoff/decay density, fix dependencies --- pyproject.toml | 2 +- src/mrinufft/trajectories/sampling.py | 17 +++++++++-------- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index fe1b1185..5e3cb71a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,7 +19,7 @@ pynfft = ["pynfft2>=1.4.3; python_version < '3.12'", "numpy>=2.0.0; python_versi pynufft = ["pynufft"] io = ["pymapvbvd"] smaps = ["scikit-image"] -wavelets = ["pywavelets"] +sampling = ["pywavelets", "scikit-learn"] autodiff = ["torch"] diff --git a/src/mrinufft/trajectories/sampling.py b/src/mrinufft/trajectories/sampling.py index bca52cbb..cbaa5877 100644 --- a/src/mrinufft/trajectories/sampling.py +++ b/src/mrinufft/trajectories/sampling.py @@ -119,10 +119,11 @@ def create_cutoff_decay_density(shape, cutoff, decay, resolution=None): Parameters ---------- shape : tuple of int - The shape of the density grid. + The shape of the density grid, analog to the field-of-view + as opposed to ``resolution`` below. cutoff : float - The k-space radius cutoff ratio between 0 and 1 within - which density remains uniform and beyond which it decays. + The ratio of the largest k-space dimension between 0 + and 1 within which density remains uniform and beyond which it decays. decay : float The polynomial decay in density beyond the cutoff ratio. resolution : np.ndarray, optional @@ -149,13 +150,13 @@ def create_cutoff_decay_density(shape, cutoff, decay, resolution=None): if not resolution: resolution = np.ones(nb_dims) - distances = np.indices(shape).astype(float) + differences = np.indices(shape).astype(float) for i in range(nb_dims): - distances[i] = distances[i] + 0.5 - shape[i] / 2 - distances[i] = distances[i] / shape[i] * resolution[i] - distances = nl.norm(distances, axis=0) + differences[i] = differences[i] + 0.5 - shape[i] / 2 + differences[i] = differences[i] / shape[i] / resolution[i] + distances = nl.norm(differences, axis=0) - cutoff = cutoff * np.max(distances) if cutoff else np.min(distances) + cutoff = cutoff * np.max(differences) if cutoff else np.min(differences) density = np.ones(shape) decay_mask = np.where(distances > cutoff, True, False) From 34144b80e62f1f7e7bd772d443e6e815af0f701f Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Tue, 10 Dec 2024 23:42:30 +0100 Subject: [PATCH 13/18] Gather extra dependencies, homogeneize extra imports --- pyproject.toml | 4 +-- src/mrinufft/extras/smaps.py | 11 ++++++-- src/mrinufft/io/siemens.py | 6 +++++ src/mrinufft/operators/stacked.py | 1 - src/mrinufft/trajectories/sampling.py | 39 ++++++++++++++++----------- 5 files changed, 40 insertions(+), 21 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 5e3cb71a..1d594e0c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,9 +17,7 @@ cufinufft = ["cufinufft<2.3", "cupy-cuda12x"] finufft = ["finufft"] pynfft = ["pynfft2>=1.4.3; python_version < '3.12'", "numpy>=2.0.0; python_version < '3.12'"] pynufft = ["pynufft"] -io = ["pymapvbvd"] -smaps = ["scikit-image"] -sampling = ["pywavelets", "scikit-learn"] +extra = ["pymapvbvd", "scikit-image", "scikit-learn", "pywavelets"] autodiff = ["torch"] diff --git a/src/mrinufft/extras/smaps.py b/src/mrinufft/extras/smaps.py index c3f7018f..5060c3c3 100644 --- a/src/mrinufft/extras/smaps.py +++ b/src/mrinufft/extras/smaps.py @@ -149,8 +149,15 @@ def low_frequency( """ # defer import to later to prevent circular import from mrinufft import get_operator - from skimage.filters import threshold_otsu, gaussian - from skimage.morphology import convex_hull_image + try: + from skimage.filters import threshold_otsu, gaussian + from skimage.morphology import convex_hull_image + except ImportError as err: + raise ImportError( + "The scikit-image module is not available. Please install " + "it along with the [extra] dependencies " + "or using `pip install scikit-image`." + ) from err k_space, samples, dc = _extract_kspace_center( kspace_data=kspace_data, diff --git a/src/mrinufft/io/siemens.py b/src/mrinufft/io/siemens.py index dd8c96d7..d172e355 100644 --- a/src/mrinufft/io/siemens.py +++ b/src/mrinufft/io/siemens.py @@ -60,6 +60,12 @@ def read_siemens_rawdat( "The mapVBVD module is not available. Please install it using " "the following command: pip install pymapVBVD" ) from err + except ImportError as err: + raise ImportError( + "The mapVBVD module is not available. Please install " + "it along with the [extra] dependencies " + "or using `pip install pymapVBVD`." + ) from err twixObj = mapVBVD(filename) if isinstance(twixObj, list): twixObj = twixObj[-1] diff --git a/src/mrinufft/operators/stacked.py b/src/mrinufft/operators/stacked.py index c3136bd1..cdb2bb62 100644 --- a/src/mrinufft/operators/stacked.py +++ b/src/mrinufft/operators/stacked.py @@ -23,7 +23,6 @@ try: import cupy as cp from cupyx.scipy import fft as cpfft - except ImportError: CUPY_AVAILABLE = False diff --git a/src/mrinufft/trajectories/sampling.py b/src/mrinufft/trajectories/sampling.py index cbaa5877..b5fb465d 100644 --- a/src/mrinufft/trajectories/sampling.py +++ b/src/mrinufft/trajectories/sampling.py @@ -4,16 +4,10 @@ import numpy.fft as nf import numpy.linalg as nl import numpy.random as nr -from sklearn.cluster import BisectingKMeans, KMeans from tqdm.auto import tqdm from .utils import KMAX -try: - import pywt as pw -except ImportError: - pw = None - def sample_from_density( nb_samples, density, method="random", *, dim_compensation="auto" @@ -59,7 +53,15 @@ def sample_from_density( "Variable density sampling with continuous trajectories." SIAM Journal on Imaging Sciences 7, no. 4 (2014): 1962-1992. """ - rng = nr.default_rng() + try: + from sklearn.cluster import BisectingKMeans, KMeans + except ImportError as err: + raise ImportError( + "The scikit-learn module is not available. Please install " + "it along with the [extra] dependencies " + "or using `pip install scikit-learn`." + ) from err + # Define dimension variables shape = np.array(density.shape) @@ -80,6 +82,7 @@ def sample_from_density( density = density / np.sum(density) # Sample using specified method + rng = nr.default_rng() if method == "random": choices = rng.choice( np.arange(max_nb_samples), @@ -259,11 +262,14 @@ def create_chauffert_density(shape, wavelet_basis, nb_wavelet_scales, verbose=Fa In 2013 IEEE 10th International Symposium on Biomedical Imaging, pp. 298-301. IEEE, 2013. """ - if pw is None: + try: + import pywt + except ImportError as err: raise ImportError( - "The PyWavelets package must be installed " - "as an additional dependency for this function." - ) + "The PyWavelets module is not available. Please install " + "it along with the [extra] dependencies " + "or using `pip install pywavelets`." + ) from err nb_dims = len(shape) indices = np.indices(shape).reshape((nb_dims, -1)).T @@ -330,11 +336,14 @@ def create_fast_chauffert_density(shape, wavelet_basis, nb_wavelet_scales): In 2013 IEEE 10th International Symposium on Biomedical Imaging, pp. 298-301. IEEE, 2013. """ - if pw is None: + try: + import pywt + except ImportError as err: raise ImportError( - "The PyWavelets package must be installed " - "as an additional dependency for this function." - ) + "The PyWavelets module is not available. Please install " + "it along with the [extra] dependencies " + "or using `pip install pywavelets`." + ) from err nb_dims = len(shape) From e7589f2ac08a64fc3b9bcc917dd2880c0108d5a9 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Tue, 10 Dec 2024 23:45:54 +0100 Subject: [PATCH 14/18] Fix imports to pywavelets --- src/mrinufft/trajectories/sampling.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/mrinufft/trajectories/sampling.py b/src/mrinufft/trajectories/sampling.py index b5fb465d..14ed18fc 100644 --- a/src/mrinufft/trajectories/sampling.py +++ b/src/mrinufft/trajectories/sampling.py @@ -62,7 +62,6 @@ def sample_from_density( "or using `pip install scikit-learn`." ) from err - # Define dimension variables shape = np.array(density.shape) nb_dims = len(shape) @@ -263,7 +262,7 @@ def create_chauffert_density(shape, wavelet_basis, nb_wavelet_scales, verbose=Fa pp. 298-301. IEEE, 2013. """ try: - import pywt + import pywt as pw except ImportError as err: raise ImportError( "The PyWavelets module is not available. Please install " @@ -337,7 +336,7 @@ def create_fast_chauffert_density(shape, wavelet_basis, nb_wavelet_scales): pp. 298-301. IEEE, 2013. """ try: - import pywt + import pywt as pw except ImportError as err: raise ImportError( "The PyWavelets module is not available. Please install " From 432df369673c3373dec37f7ef30f239b65925add Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Tue, 10 Dec 2024 23:48:33 +0100 Subject: [PATCH 15/18] Fix black format in smaps --- src/mrinufft/extras/smaps.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/mrinufft/extras/smaps.py b/src/mrinufft/extras/smaps.py index 5060c3c3..6bb0ba3b 100644 --- a/src/mrinufft/extras/smaps.py +++ b/src/mrinufft/extras/smaps.py @@ -149,6 +149,7 @@ def low_frequency( """ # defer import to later to prevent circular import from mrinufft import get_operator + try: from skimage.filters import threshold_otsu, gaussian from skimage.morphology import convex_hull_image From 4ff1bf0ff5cbe22804d9629a7afb5f2a8a443694 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Tue, 10 Dec 2024 23:53:52 +0100 Subject: [PATCH 16/18] Remove except duplicate in io/siemens --- src/mrinufft/io/siemens.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/mrinufft/io/siemens.py b/src/mrinufft/io/siemens.py index d172e355..e0ea009d 100644 --- a/src/mrinufft/io/siemens.py +++ b/src/mrinufft/io/siemens.py @@ -55,11 +55,6 @@ def read_siemens_rawdat( """ try: from mapvbvd import mapVBVD - except ImportError as err: - raise ImportError( - "The mapVBVD module is not available. Please install it using " - "the following command: pip install pymapVBVD" - ) from err except ImportError as err: raise ImportError( "The mapVBVD module is not available. Please install " From e60da32a1b8c56aca43dc7ad11a79a2ee8897933 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Thu, 12 Dec 2024 18:59:58 +0100 Subject: [PATCH 17/18] Update CI with new [extra] dependencies --- .github/workflows/test-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test-ci.yml b/.github/workflows/test-ci.yml index 5ae8cd2d..4720af6a 100644 --- a/.github/workflows/test-ci.yml +++ b/.github/workflows/test-ci.yml @@ -206,7 +206,7 @@ jobs: shell: bash run: | python -m pip install --upgrade pip - python -m pip install -e .[test,dev] + python -m pip install -e .[extra,test,dev] python -m pip install finufft pooch brainweb-dl torch fastmri - name: Install GPU related interfaces From 2f34957a22f73ad218932043f390ee989738695c Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Mon, 16 Dec 2024 18:33:25 +0100 Subject: [PATCH 18/18] [docs] Trigger documentation build