Merge c601ace into 9a1671d

.gitignore

@@ -5,6 +5,7 @@
 .pytest_cache/
 .ipynb_checkpoints/
 .vscode/
+build/
 
 # Project's data
 testing_data

installer/windows_install.bat

@@ -42,7 +42,7 @@ REM Installing Shimming Toolbox
 copy "%ST_SOURCE_FILES%\config\dcm2bids.json" "%ST_DIR%\dcm2bids.json" || goto error
 
 cd "%ST_SOURCE_FILES%"
-"%ST_DIR%\%PYTHON_DIR%\python.exe" -m pip install -e ".[docs,dev]" --no-warn-script-location || goto error
+"%ST_DIR%\%PYTHON_DIR%\python.exe" -m pip install . --no-warn-script-location || goto error
 
 REM Create launchers for Shimming Toolbox
 set "BIN_DIR=bin"
@@ -52,6 +52,17 @@ for %%f in ("%ST_DIR%\%PYTHON_DIR%\Scripts\st_*.*") do (
 	copy "%%f" "%ST_DIR%\%BIN_DIR%" || goto error
 )
 
+REM Copy dcm2niix to the bin directory, there are 2 places it might be
+set "PATH_DCM2NIIX=%ST_DIR%\%PYTHON_DIR%\Library\bin\dcm2niix.exe"
+if exist "%PATH_DCM2NIIX%" (copy "%PATH_DCM2NIIX%" "%ST_DIR%\%BIN_DIR%" || goto error) else (
+    set "PATH_DCM2NIIX=%ST_DIR%\%PYTHON_DIR%\Scripts\dcm2niix.exe"
+    if exist "%PATH_DCM2NIIX%" (copy "%PATH_DCM2NIIX%" "%ST_DIR%\%BIN_DIR%" || goto error) else (
+        echo "dcm2niix.exe not found"
+        goto error
+    )
+)
+
+
 REM Add scripts to the User's path
 for /F "skip=2 tokens=2,*" %%A in ('reg.exe query "HKEY_CURRENT_USER\Environment" /v path') do set "OLD_PATH=%%B"
 REM If OLD_PATH is not an empty string

shimmingtoolbox/coils/coil.py

@@ -5,9 +5,8 @@
 import numpy as np
 from typing import Tuple
 
-from shimmingtoolbox.coils.spher_harm_basis import siemens_basis, ge_basis
+from shimmingtoolbox.coils.spher_harm_basis import siemens_basis, ge_basis, philips_basis, SHIM_CS
 from shimmingtoolbox.coils.coordinates import generate_meshgrid
-from shimmingtoolbox.shim.shim_utils import shim_cs
 
 logger = logging.getLogger(__name__)
 
@@ -123,8 +122,8 @@ def __init__(self, dim_volume, affine, constraints, order, manufacturer=""):
         self.order = order
 
         manufacturer = manufacturer.upper()
-        if manufacturer in shim_cs:
-            self.coord_system = shim_cs[manufacturer.upper()]
+        if manufacturer in SHIM_CS:
+            self.coord_system = SHIM_CS[manufacturer.upper()]
         else:
             logger.warning(f"Unknown manufacturer {manufacturer}, assuming RAS")
             self.coord_system = 'RAS'
@@ -155,6 +154,9 @@ def _create_coil_profile(self, dim, manufacturer=None):
             elif manufacturer == 'GE':
                 profile_orders = ge_basis(mesh1, mesh2, mesh3, orders=tuple(range(1, self.order + 1)),
                                           shim_cs=self.coord_system)
+            elif manufacturer == 'PHILIPS':
+                profile_orders = philips_basis(mesh1, mesh2, mesh3, orders=tuple(range(1, self.order + 1)),
+                                               shim_cs=self.coord_system)
             else:
                 logger.warning(f"{manufacturer} manufacturer not implemented. Outputting in Hz, uT/m, uT/m^2 for order "
                                f"0, 1 and 2 respectively")

shimmingtoolbox/coils/spher_harm_basis.py

@@ -8,10 +8,13 @@
 
 logger = logging.getLogger(__name__)
 
-GYROMAGNETIC_RATIO = 42.5774785178325552  # [MHz/T]
+GYROMAGNETIC_RATIO = 42.5774785178325552  # [MHz/T] or equivalently [Hz/uT]
+SHIM_CS = {'SIEMENS': 'LAI',
+           'GE': 'LPI',
+           'PHILIPS': 'RPI'}
 
 
-def siemens_basis(x, y, z, orders=(1, 2), shim_cs='LAI'):
+def siemens_basis(x, y, z, orders=(1, 2), shim_cs=SHIM_CS['SIEMENS']):
     """
     The function first wraps ``shimmingtoolbox.coils.spherical_harmonics`` to generate 1st and 2nd order spherical
     harmonic ``basis`` fields at the grid positions given by arrays ``X,Y,Z``. *Following Siemens convention*,
@@ -53,7 +56,8 @@ def siemens_basis(x, y, z, orders=(1, 2), shim_cs='LAI'):
 
     # Create spherical harmonics from first to second order
     all_orders = np.array(range(1, 3))
-    spher_harm = scaled_spher_harm(x, y, z, all_orders, shim_cs=shim_cs)
+    flip = get_flip_matrix(shim_cs, manufacturer='SIEMENS', xyz=True)
+    spher_harm = scaled_spher_harm(x * flip[0], y * flip[1], z * flip[2], all_orders)
 
     # Reorder according to siemens convention: X, Y, Z, Z2, ZX, ZY, X2-Y2, XY
     reordered_spher = _reorder_to_siemens(spher_harm)
@@ -94,7 +98,7 @@ def _reorder_to_siemens(spher_harm):
     return reordered
 
 
-def ge_basis(x, y, z, orders=(1, 2), shim_cs='LPI'):
+def ge_basis(x, y, z, orders=(1, 2), shim_cs=SHIM_CS['GE']):
     """
     The function first wraps ``shimmingtoolbox.coils.spher_harm_basis.scaled_spher_harm`` to generate 1st and 2nd
     order spherical harmonic ``basis`` fields at the grid positions given by arrays ``X,Y,Z``.
@@ -133,7 +137,8 @@ def ge_basis(x, y, z, orders=(1, 2), shim_cs='LPI'):
 
     # Create spherical harmonics from first to second order
     all_orders = np.array(range(1, 3))
-    spher_harm = scaled_spher_harm(x, y, z, all_orders, shim_cs=shim_cs)
+    flip = get_flip_matrix(shim_cs, manufacturer='GE', xyz=True)
+    spher_harm = scaled_spher_harm(x * flip[0], y * flip[1], z * flip[2], all_orders)
 
     # The following matrix (8 x 5) refers to the following:
     #  \  xy, zy, zx, XY, z2
@@ -175,25 +180,30 @@ def ge_basis(x, y, z, orders=(1, 2), shim_cs='LPI'):
     for i_channel in range(reordered_spher.shape[3]):
         if i_channel in [0, 1, 2]:
 
-            # Rescale to unit-shims that are G/cm
+            # Rescale to unit-shims that are XXXXX
             # They are currently in uT/m
-            # 1G = 1e-4T, 1T = 1e4G
-            # uT/m --> G/cm = reordered_spher * (1/1e6) * 1e4  * 100 = reordered_sphere
-            scaled[..., i_channel] = reordered_spher[..., i_channel]
+            # Todo: This seems to be the appropriate scaling factor, we need to verify the units
+            scaled[..., i_channel] = reordered_spher[..., i_channel] / 10
+
         else:
             # Since reordered_spher contains the values of 1uT/m^2 in Hz/mm^2. We simply multiply by the amount of
             # uT/m^2 / A
             # This gives us a value in Hz/mm^2 / A which we need to modify to Hz/mm^2 / mA
             scaled[..., i_channel] = np.matmul(reordered_spher[..., 3:], orders_to_order2_uT[i_channel - 3, :]) / 1000
+            # Todo: We need a /2 between expected zx, zy, xy results and calculated results
+            if i_channel in [3, 4, 5]:
+                scaled[..., i_channel] /= 2
+
+    ordered = _reorder_to_ge2(scaled)
 
     # Output according to the specified order
     range_per_order = {1: list(range(3)), 2: list(range(3, 8))}
     length_dim3 = np.sum([len(values) for key, values in range_per_order.items() if key in orders])
-    output = np.zeros(scaled[..., 0].shape + (length_dim3,), dtype=scaled.dtype)
+    output = np.zeros(ordered[..., 0].shape + (length_dim3,), dtype=ordered.dtype)
     start_index = 0
     for order in orders:
         end_index = start_index + len(range_per_order[order])
-        output[..., start_index:end_index] = scaled[..., range_per_order[order]]
+        output[..., start_index:end_index] = ordered[..., range_per_order[order]]
         # prep next iteration
         start_index = end_index
 
@@ -204,7 +214,7 @@ def _reorder_to_ge(spher_harm):
     """
     Reorder 1st - 2nd order coefficients along the last dim. From
     1. Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 (output by shimmingtoolbox.coils.spherical_harmonics.spherical_harmonics), to
-    2. x, y, z, xy, zy, zx, X2 - Y2, z2 (in line with GE shims)
+    2. x, y, z, xy, zy, zx, X2 - Y2, z2 (in line with GE scaling matrix)
 
     Args:
         spher_harm (numpy.ndarray): Coefficients with the last dimensions representing the different order channels.
@@ -222,7 +232,119 @@ def _reorder_to_ge(spher_harm):
     return reordered
 
 
-def scaled_spher_harm(x, y, z, orders=(1, 2), shim_cs='ras'):
+def _reorder_to_ge2(spher_harm):
+    """
+    Reorder 1st - 2nd order coefficients along the last dim. From
+    1. *x, y, z, xy, zy, zx, X2 - Y2, z2* (scaling matrix)
+    2. *X, Y, Z, Z2, ZX, ZY, X2-Y2, XY* (in line with GE shims)
+
+    Args:
+        spher_harm (numpy.ndarray): Coefficients with the last dimensions representing the different order channels.
+                                    ``spher_harm.shape[-1]`` must equal 8.
+
+    Returns:
+        numpy.ndarray: Coefficients ordered following GE convention
+    """
+
+    if spher_harm.shape[-1] != 8:
+        raise RuntimeError("Input arrays should have 4th dimension's shape equal to 8")
+
+    reordered = spher_harm[..., [0, 1, 2, 7, 5, 4, 6, 3]]
+
+    return reordered
+
+
+def philips_basis(x, y, z, orders=(1, 2), shim_cs=SHIM_CS['PHILIPS']):
+    """
+    The function first wraps ``shimmingtoolbox.coils.spherical_harmonics`` to generate 1st and 2nd order spherical
+    harmonic ``basis`` fields at the grid positions given by arrays ``X,Y,Z``. *Following Philips convention*,
+    ``basis`` is then:
+
+        - Rescaled to Hz/unit-shim, where "unit-shim" refers to:
+
+            - 1 milli-T/m for *X,Y,Z* gradients (= 42.576 Hz/mm)
+            - 1 milli-T/m^2 for 2nd order terms (= 0.042576 Hz/mm^2)
+
+        - Reordered along the 4th dimension as *X, Y, Z, Z2, ZX, ZY, X2-Y2, XY*
+
+    The returned ``basis`` is thereby in the form of ideal "shim reference maps", ready for optimization.
+
+    Args:
+        x (numpy.ndarray): 3-D arrays of grid coordinates, "Left->Right" grid coordinates in the patient coordinate
+                           system (i.e. NIfTI reference (RAS), units of mm)
+        y (numpy.ndarray): 3-D arrays of grid coordinates (same shape as x). "Posterior->Anterior" grid coordinates in
+                           the patient coordinate system (i.e. NIfTI reference (RAS), units of mm)
+        z (numpy.ndarray): 3-D arrays of grid coordinates (same shape as x). "Inferior->Superior" grid coordinates in
+                           the patient coordinate system (i.e. NIfTI reference, units of mm)
+        orders (tuple): Degrees of the desired terms in the series expansion, specified as a vector of non-negative
+                        integers (``(0:1:n)`` yields harmonics up to n-th order, implemented 1st and 2nd order)
+        shim_cs (str): Coordinate system of the shims. Letter 1 'R' or 'L', letter 2 'A' or 'P', letter 3 'S' or 'I'.
+
+    Returns:
+        numpy.ndarray: 4-D array of spherical harmonic basis fields
+
+    Note:
+        Philips coordinate system has its x in the AP direction and y axis in the RL direction. Therefore, channel 0 (x)
+        changes along axis 1 and channel 1 (y) changes along axis 0.
+    """
+    # Check inputs
+    _check_basis_inputs(x, y, z, orders)
+
+    # Create spherical harmonics from first to second order
+    all_orders = np.array(range(1, 3))
+    # Philips' y and x axis are flipped (x is AP, y is RL)
+    flip = get_flip_matrix(shim_cs, manufacturer='Philips', xyz=True)
+    spher_harm = scaled_spher_harm(y * flip[0], x * flip[1], z * flip[2], all_orders)
+
+    # Scale according to Philips convention
+    # milli-T/m for order 1, milli-T/m^2 for order 2
+    # uT/m * 1e3 = mT/m, uT/m^2 * 1e3 = mT/m^2
+    spher_harm *= 1000
+
+    # Reorder according to philips convention: X, Y, Z, Z2, ZX, ZY, X2-Y2, XY
+    reordered_spher = _reorder_to_philips(spher_harm)
+
+    # Todo: We need a /2 between expected zx, zy results and calculated results (Similar to GE but not with XY)
+    reordered_spher[..., 4] /= 2
+    reordered_spher[..., 5] /= 2
+
+    # Select order
+    range_per_order = {1: list(range(3)), 2: list(range(3, 8))}
+    length_dim3 = np.sum([len(values) for key, values in range_per_order.items() if key in orders])
+    output = np.zeros(reordered_spher[..., 0].shape + (length_dim3,), dtype=reordered_spher.dtype)
+    start_index = 0
+    for order in orders:
+        end_index = start_index + len(range_per_order[order])
+        output[..., start_index:end_index] = reordered_spher[..., range_per_order[order]]
+        # prep next iteration
+        start_index = end_index
+
+    return output
+
+
+def _reorder_to_philips(spher_harm):
+    """
+    Reorder 1st - 2nd order coefficients along the last dim. From
+    1. Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 (output by shimmingtoolbox.coils.spherical_harmonics.spherical_harmonics), to
+    2. X, Y, Z, Z2, ZX, ZY, X2 - Y2, XY (in line with Philips shims)
+
+    Args:
+        spher_harm (numpy.ndarray): Coefficients with the last dimensions representing the different order channels.
+                                    ``spher_harm.shape[-1]`` must equal 8.
+
+    Returns:
+        numpy.ndarray: Coefficients ordered following Philips convention
+    """
+
+    if spher_harm.shape[-1] != 8:
+        raise RuntimeError("Input arrays should have 4th dimension's shape equal to 8")
+
+    reordered = spher_harm[..., [2, 0, 1, 5, 6, 4, 7, 3]]
+
+    return reordered
+
+
+def scaled_spher_harm(x, y, z, orders=(1, 2)):
     """ The function first wraps ``shimmingtoolbox.coils.spherical_harmonics`` to generate 1st and 2nd order spherical
     harmonic ``basis`` fields at the grid positions given by arrays ``X,Y,Z``. It is then:
 
@@ -253,15 +375,13 @@ def scaled_spher_harm(x, y, z, orders=(1, 2), shim_cs='ras'):
     spher_harm = spherical_harmonics(all_orders, x, y, z)
     # 1. Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 (output by shimmingtoolbox.coils.spherical_harmonics.spherical_harmonics)
 
-    spher_harm_cs = get_flip_matrix(shim_cs) * spher_harm
-
     # scale according to
     # - 1 micro-T/m for *X,Y,Z* gradients (= 0.042576 Hz/mm)
     # - 1 micro-T/m^2 for 2nd order terms (= 0.000042576 Hz/mm^2)
     scaling_factors = _get_scaling_factors()
-    scaled = np.zeros_like(spher_harm_cs)
-    for i_channel in range(0, spher_harm_cs.shape[3]):
-        scaled[:, :, :, i_channel] = scaling_factors[i_channel] * spher_harm_cs[:, :, :, i_channel]
+    scaled = np.zeros_like(spher_harm)
+    for i_channel in range(0, spher_harm.shape[3]):
+        scaled[:, :, :, i_channel] = scaling_factors[i_channel] * spher_harm[:, :, :, i_channel]
 
     # 1 uT/m, 1 uT/m2 in Hz/mm, Hz/mm2
     return scaled
@@ -302,7 +422,7 @@ def _get_scaling_factors():
     i_x1y1 = np.nonzero((x_iso == 1) & (y_iso == 1) & (z_iso == 0))
 
     # order the reference indices like the sh field terms
-    # 1. Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 (output by shimmingtoolbox.coils.spherical_harmonics.spherical_harmonics), to
+    # 1. Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 (output by shimmingtoolbox.coils.spherical_harmonics.spherical_harmonics)
     i_ref = [i_y1, i_z1, i_x1, i_x1y1, i_y1z1, i_z1, i_x1z1, i_x1]
 
     # distance from iso/origin to adopted reference point[units: mm]
@@ -374,16 +494,18 @@ def get_flip_matrix(shim_cs='ras', manufacturer=None, xyz=False):
     elif manufacturer == 'GE':
         out_matrix = _reorder_to_ge(out_matrix)
     elif manufacturer == 'PHILIPS':
-        logger.warning("Philips shim CS not implemented yet")
+        out_matrix = _reorder_to_philips(out_matrix)
     else:
         # Do not reorder if the manufacturer is not specified
         pass
 
     if xyz:
-        # X, Y, Z
+        # None: Y, Z, X
+        # GE/Siemens/Philips: X, Y, Z
         return out_matrix[:3]
     else:
         # None: Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2
         # GE: x, y, z, xy, zy, zx, X2 - Y2, z2
         # Siemens: X, Y, Z, Z2, ZX, ZY, X2 - Y2, XY
+        # Philips: X, Y, Z, Z2, ZX, ZY, X2 - Y2, XY
         return out_matrix