Use and add conversions of units to the basis functions

shimmingtoolbox/coils/spher_harm_basis.py

@@ -5,10 +5,12 @@
 import numpy as np
 
 from shimmingtoolbox.coils.spherical_harmonics import spherical_harmonics
+from shimmingtoolbox.conversion import (hz_per_cm_to_micro_tesla_per_m, hz_per_cm2_to_micro_tesla_per_m2,
+                                        metric_unit_to_metric_unit, unit_per_metric_unit_to_unit_per_metric_unit,
+                                        tesla_to_hz)
 
 logger = logging.getLogger(__name__)
 
-GYROMAGNETIC_RATIO = 42.5774785178325552  # [MHz/T] or equivalently [Hz/uT]
 SHIM_CS = {'SIEMENS': 'LAI',
            'GE': 'LPI',
            'PHILIPS': 'RPI'}
@@ -164,12 +166,10 @@ def _reorder_scaling_to_shim(coefs):
 
             # Scale
             orders_to_order2_ut = np.zeros_like(orders_to_order2)
-            # Hz/cm2/A, -> uT/m2/A = order2_to_order2 * 1e6 * (100 ** 2) / (GYROMAGNETIC_RATIO * 1e6)
-            # = order2_to_order2 * (100 ** 2) / GYROMAGNETIC_RATIO
-            orders_to_order2_ut[4:] = orders_to_order2[0:5] * (100 ** 2) / GYROMAGNETIC_RATIO
-            # Hz/cm/A, -> uT/m/A = order1_to_order2 * 1e6 * 100 / (GYROMAGNETIC_RATIO * 1e6)
-            # = order2_to_order2 * 100 / GYROMAGNETIC_RATIO
-            orders_to_order2_ut[1:4] = orders_to_order2[5:8] * 100 / GYROMAGNETIC_RATIO
+            # Hz/cm2/A, -> uT/m2/A = order2_to_order2 * 1e6 * (100 ** 2) / GYROMAGNETIC_RATIO
+            orders_to_order2_ut[4:] = hz_per_cm2_to_micro_tesla_per_m2(orders_to_order2[:5])
+            # Hz/cm/A, -> uT/m/A = order1_to_order2 * 1e6 * 100 / GYROMAGNETIC_RATIO
+            orders_to_order2_ut[1:4] = hz_per_cm_to_micro_tesla_per_m(orders_to_order2[5:8])
             # Hz/A
             orders_to_order2_ut[0] = orders_to_order2[8]
 
@@ -187,8 +187,10 @@ def _reorder_scaling_to_shim(coefs):
                 reordered_spher_necessary = {key: reordered_spher[key] for key in (0, 1, 2)}
                 reordered_spher_array = convert_spher_harm_to_array(reordered_spher_necessary)
 
-                scaled[2][..., i_channel] = np.matmul(reordered_spher_array,
-                                                      orders_to_order2_ut[:, i_channel]) / 1000
+                scaled[2][..., i_channel] = np.matmul(reordered_spher_array, orders_to_order2_ut[:, i_channel])
+                scaled[2][..., i_channel] = unit_per_metric_unit_to_unit_per_metric_unit(scaled[2][..., i_channel],
+                                                                                         '',
+                                                                                         'm')
                 # Todo: We need a /2 between expected zx, zy, xy results and calculated results
                 if i_channel in [0, 1, 2]:
                     scaled[2][..., i_channel] /= 2
@@ -253,7 +255,8 @@ def philips_basis(x, y, z, orders=(1, 2), shim_cs=SHIM_CS['PHILIPS']):
 
     # 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
+    # We currently have 1uT/m scaled to Hz/m. We want the equivalent Hz/m for 1mT/m
+    # order1: *1e3, order2: *1e3
     for order in orders:
         if order == 0:
             continue
@@ -423,7 +426,7 @@ def _reorder_order3(sph, manuf):
 def _get_scaling_factors(orders):
     """
     Get scaling factors for the 1st/2nd/3rd order spherical harmonic
-    fields for rescaling them to 1 uT/unit-shim in units of Hz/mm:
+    fields for rescaling them to 1 uT/unit-shim in units of Hz/mm^x:
 
     Gx, Gy, and Gz should yield 1 micro-T of field shift per metre: equivalently, 0.042576 Hz/mm
 
@@ -488,8 +491,10 @@ def _get_scaling_factors(orders):
         for i in range(channels_per_order(order)):
             field = sh[:, :, :, i_ch]
             if order != 0:
-                scaling_factors[i_ch] = (GYROMAGNETIC_RATIO * ((r[order][i] * 0.001) ** order) /
-                                         field[iref[order][i]][0])
+                # 1uT in Hz
+                hz = tesla_to_hz(metric_unit_to_metric_unit(1, 'u', ''))
+                # 1 mm^order in m^order
+                scaling_factors[i_ch] = hz * ((r[order][i] * 0.001) ** order) / field[iref[order][i]][0]
             else:
                 # B0 term needs to be flipped. Producing a 300 Hz in f means that the coil profile should be negative.
                 # (f0 - f)

shimmingtoolbox/conversion.py

@@ -3,7 +3,56 @@
 
 import math
 
-GYROMAGNETIC_RATIO = 42.577478518e+6  # [Hz/T]
+GYROMAGNETIC_RATIO = 42.5774785178325552e+6  # [Hz/T]
+METRIC_PREFIXES = {
+    'n': 1e-9,  # 'nano'
+    'u': 1e-6,  # 'micro'
+    'm': 1e-3,  # 'milli'
+    'c': 1e-2,  # 'centi'
+    '': 1,  # 'base'
+    'b': 1,  # 'base'
+    'h': 1e2,  # 'hecto'
+    'k': 1e3,  # 'kilo'
+    'M': 1e6,  # 'mega'
+    'G': 1e9  # 'giga'
+}
+
+
+def metric_unit_to_metric_unit(x, prefix_in, prefix_out, power=1):
+    """ Convert units with metric prefixes to other metric prefixes (i.e. T to uT)
+
+    Args:
+        x: Float or array of floats
+        prefix_in (str): Prefix of the input unit
+        prefix_out (str): Prefix of the output unit
+        power (int): Power of the unit, use 2 for squared units (i.e. m^2)
+
+    Returns:
+        Float or array of floats converted
+    """
+    if prefix_in not in METRIC_PREFIXES or prefix_out not in METRIC_PREFIXES:
+        raise ValueError('Invalid prefix')
+    factor = (METRIC_PREFIXES[prefix_in] / METRIC_PREFIXES[prefix_out]) ** power
+    return x * factor
+
+
+def unit_per_metric_unit_to_unit_per_metric_unit(x, prefix_in, prefix_out, power=1):
+    """ Convert units per metric prefixes to other metric prefixes (i.e. unit/cm to unit/m)
+
+    Args:
+        x: Float or array of floats
+        prefix_in (str): Prefix of the input unit
+        prefix_out (str): Prefix of the output unit
+        power (int): Power of the unit, use 2 for squared units (i.e. m^2)
+
+    Returns:
+        Float or array of floats converted
+    """
+    if prefix_in not in METRIC_PREFIXES or prefix_out not in METRIC_PREFIXES:
+        raise ValueError('Invalid prefix')
+
+    factor = (METRIC_PREFIXES[prefix_out] / METRIC_PREFIXES[prefix_in]) ** power
+    return x * factor
 
 
 def hz_to_rad_per_sec(hz):
@@ -26,6 +75,20 @@ def tesla_to_hz(tesla):
     return tesla * GYROMAGNETIC_RATIO
 
 
+def hz_to_micro_tesla(hz):
+    """ Convert Hz to microTesla """
+    tesla = hz_to_tesla(hz)
+    micro_tesla = metric_unit_to_metric_unit(tesla, '', 'u')
+    return micro_tesla
+
+
+def micro_tesla_to_hz(micro_tesla):
+    """ Convert microTesla to Hz """
+    tesla = metric_unit_to_metric_unit(micro_tesla, 'u', '')
+    hz = tesla_to_hz(tesla)
+    return hz
+
+
 def rad_per_sec_to_rad(rad_per_sec, dt):
     """ Convert rad/s to rad """
     return rad_per_sec * dt
@@ -64,27 +127,17 @@ def rad_to_tesla(rad, dt):
     return tesla
 
 
-def milli_tesla_to_tesla(milli_tesla):
-    """ Convert milliTesla to Tesla """
-    return milli_tesla * 1e-3
-
-
-def tesla_to_milli_tesla(tesla):
-    """ Convert Tesla to milliTesla """
-    return tesla * 1e3
-
-
 def milli_tesla_to_rad(milli_tesla, dt):
     """ Convert milliTesla to rad """
-    tesla = milli_tesla_to_tesla(milli_tesla)
+    tesla = metric_unit_to_metric_unit(milli_tesla, 'm', '')
     rad = tesla_to_rad(tesla, dt)
     return rad
 
 
 def rad_to_milli_tesla(rad, dt):
     """ Convert rad to milliTesla """
     tesla = rad_to_tesla(rad, dt)
-    milli_tesla = tesla_to_milli_tesla(tesla)
+    milli_tesla = metric_unit_to_metric_unit(tesla, '', 'm')
     return milli_tesla
 
 
@@ -110,3 +163,31 @@ def rad_to_gauss(rad, dt):
     tesla = rad_to_tesla(rad, dt)
     gauss = tesla_to_gauss(tesla)
     return gauss
+
+
+def hz_per_cm_to_micro_tesla_per_m(hz_per_cm):
+    """ Convert Hz/cm to microTesla/m """
+    micro_tesla_per_cm = hz_to_micro_tesla(hz_per_cm)
+    micro_tesla_per_m = unit_per_metric_unit_to_unit_per_metric_unit(micro_tesla_per_cm, 'c', '')
+    return micro_tesla_per_m
+
+
+def micro_tesla_per_m_to_hz_per_cm(micro_tesla_per_m):
+    """ Convert microTesla/m to Hz/cm """
+    micro_tesla_per_cm = unit_per_metric_unit_to_unit_per_metric_unit(micro_tesla_per_m, '', 'c')
+    hz_per_cm = micro_tesla_to_hz(micro_tesla_per_cm)
+    return hz_per_cm
+
+
+def hz_per_cm2_to_micro_tesla_per_m2(hz_per_cm2):
+    """ Convert Hz/cm^2 to microTesla/m^2 """
+    micro_tesla_per_cm2 = hz_to_micro_tesla(hz_per_cm2)
+    micro_tesla_per_m2 = unit_per_metric_unit_to_unit_per_metric_unit(micro_tesla_per_cm2, 'c', '', 2)
+    return micro_tesla_per_m2
+
+
+def micro_tesla_per_m2_to_hz_per_cm2(micro_tesla_per_m2):
+    """ Convert microTesla/m^2 to Hz/cm^2 """
+    micro_tesla_per_cm2 = unit_per_metric_unit_to_unit_per_metric_unit(micro_tesla_per_m2, '', 'c', 2)
+    hz_per_cm2 = micro_tesla_to_hz(micro_tesla_per_cm2)
+    return hz_per_cm2

test/test_conversion.py

@@ -8,7 +8,10 @@
 from shimmingtoolbox.conversion import (hz_to_rad_per_sec, rad_per_sec_to_hz, hz_to_tesla, tesla_to_hz,
                                         rad_per_sec_to_rad, rad_to_rad_per_sec, hz_to_rad, rad_to_hz, tesla_to_rad,
                                         rad_to_tesla, milli_tesla_to_rad, rad_to_milli_tesla, gauss_to_rad,
-                                        rad_to_gauss)
+                                        rad_to_gauss, hz_to_micro_tesla, micro_tesla_to_hz,
+                                        hz_per_cm2_to_micro_tesla_per_m2, micro_tesla_per_m2_to_hz_per_cm2,
+                                        hz_per_cm_to_micro_tesla_per_m, micro_tesla_per_m_to_hz_per_cm,
+                                        metric_unit_to_metric_unit, unit_per_metric_unit_to_unit_per_metric_unit)
 
 
 def test_hz_to_rad_per_sec():
@@ -22,15 +25,23 @@ def test_rad_per_sec_to_hz():
 
 
 def test_hz_to_tesla():
-    assert hz_to_tesla(42.577478518e+6) == 1
+    assert math.isclose(hz_to_tesla(42.577478518e+6), 1)
     assert hz_to_tesla(0) == 0
 
 
 def test_tesla_to_hz():
-    assert tesla_to_hz(1) == 42.577478518e+6
+    assert math.isclose(tesla_to_hz(1), 42.577478518e+6)
     assert tesla_to_hz(0) == 0
 
 
+def test_hz_to_micro_tesla():
+    assert math.isclose(hz_to_micro_tesla(42.577478518), 1)
+
+
+def test_micro_tesla_to_hz():
+    assert math.isclose(micro_tesla_to_hz(1), 42.577478518)
+
+
 def test_rad_per_sec_to_rad():
     assert rad_per_sec_to_rad(1, 0.1) == 0.1
     assert rad_per_sec_to_rad(0, 0.1) == 0
@@ -52,12 +63,12 @@ def test_rad_to_hz():
 
 
 def test_tesla_to_rad():
-    assert tesla_to_rad(1, 0.1) == 0.2 * math.pi * 42.577478518e+6
+    assert math.isclose(tesla_to_rad(1, 0.1), 0.2 * math.pi * 42.577478518e+6)
     assert tesla_to_rad(0, 0.1) == 0
 
 
 def test_rad_to_tesla():
-    assert rad_to_tesla(2 * math.pi, 0.1) == 10 / 42.577478518e+6
+    assert math.isclose(rad_to_tesla(2 * math.pi, 0.1), 10 / 42.577478518e+6)
     assert rad_to_tesla(0, 0.1) == 0
 
 
@@ -67,15 +78,63 @@ def test_milli_tesla_to_rad():
 
 
 def test_rad_to_milli_tesla():
-    assert rad_to_milli_tesla(2 * math.pi, 0.1) == 10e3 / 42.577478518e+6
+    assert math.isclose(rad_to_milli_tesla(2 * math.pi, 0.1), 10e3 / 42.577478518e+6)
     assert rad_to_milli_tesla(0, 0.1) == 0
 
 
 def test_gauss_to_rad():
-    assert gauss_to_rad(1, 0.1) == 0.2e-4 * math.pi * 42.577478518e+6
+    assert math.isclose(gauss_to_rad(1, 0.1), 0.2e-4 * math.pi * 42.577478518e+6)
     assert gauss_to_rad(0, 0.1) == 0
 
 
 def test_rad_to_gauss():
-    assert rad_to_gauss(2 * math.pi, 0.1) == 10e4 / 42.577478518e+6
+    assert math.isclose(rad_to_gauss(2 * math.pi, 0.1), 10e4 / 42.577478518e+6)
     assert rad_to_gauss(0, 0.1) == 0
+
+
+def test_cm_to_m():
+    assert metric_unit_to_metric_unit(500, 'c', '') == 5
+
+
+def test_m_to_cm():
+    assert metric_unit_to_metric_unit(5, 'b', 'c') == 500
+
+
+def test_cm2_to_m2():
+    assert metric_unit_to_metric_unit(50000, 'c', '', 2) == 5
+
+
+def test_m2_to_cm2():
+    assert metric_unit_to_metric_unit(5, '', 'c', 2) == 50000
+
+
+def test_unit_per_cm_to_unit_per_m():
+    assert unit_per_metric_unit_to_unit_per_metric_unit(1, 'c', '') == 1e2
+
+
+def test_unit_per_m_to_unit_per_cm():
+    assert unit_per_metric_unit_to_unit_per_metric_unit(1e2, '', 'c') == 1
+
+
+def test_hz_per_cm_to_micro_tesla_per_m():
+    assert math.isclose(hz_per_cm_to_micro_tesla_per_m(0.42577478518), 1)
+
+
+def test_micro_tesla_per_m_to_hz_per_cm():
+    assert math.isclose(micro_tesla_per_m_to_hz_per_cm(1), 0.42577478518)
+
+
+def test_unit_per_cm2_to_unit_per_m2():
+    assert unit_per_metric_unit_to_unit_per_metric_unit(1, 'c', '', 2) == 1e4
+
+
+def test_unit_per_m2_to_unit_per_cm2():
+    assert unit_per_metric_unit_to_unit_per_metric_unit(1e4, '', 'c', 2) == 1
+
+
+def test_hz_per_cm2_to_micro_tesla_per_m2():
+    assert math.isclose(hz_per_cm2_to_micro_tesla_per_m2(0.0042577478518), 1)
+
+
+def test_micro_tesla_per_m2_to_hz_per_cm2():
+    assert math.isclose(micro_tesla_per_m2_to_hz_per_cm2(1), 0.0042577478518)