diff --git a/shimmingtoolbox/coils/coordinates.py b/shimmingtoolbox/coils/coordinates.py
index 802ee951..cb8fcd44 100644
--- a/shimmingtoolbox/coils/coordinates.py
+++ b/shimmingtoolbox/coils/coordinates.py
@@ -2,9 +2,8 @@
 # -*- coding: utf-8 -*
 # Deals with coordinate systems, going from voxel-based to physical-based coordinates.
 
-# TODO: create a test for this API
-
 import numpy as np
+from nibabel.affines import apply_affine
 
 
 def generate_meshgrid(dim, affine):
@@ -19,15 +18,16 @@ def generate_meshgrid(dim, affine):
     """
 
     nx, ny, nz = dim
-    coord_vox = np.meshgrid(np.array(range(nx)), np.array(range(ny)), np.array(range(nz)), indexing='ij')
-    coord_phys = [np.zeros_like(coord_vox[0]), np.zeros_like(coord_vox[1]), np.zeros_like(coord_vox[2])]
+    coord_vox = np.meshgrid(np.arange(nx), np.arange(ny), np.arange(nz), indexing='ij')
+    coord_phys = [np.zeros_like(coord_vox[0]).astype(float),
+                  np.zeros_like(coord_vox[1]).astype(float),
+                  np.zeros_like(coord_vox[2]).astype(float)]
 
     # TODO: Better code
     for ix in range(nx):
         for iy in range(ny):
             for iz in range(nz):
-                coord_phys_list = \
-                    np.dot([coord_vox[i][ix, iy, iz] for i in range(3)], affine[0:3, 0:3]) + affine[0:3, 3]
+                coord_phys_list = apply_affine(affine, [coord_vox[i][ix, iy, iz] for i in range(3)])
                 for i in range(3):
                     coord_phys[i][ix, iy, iz] = coord_phys_list[i]
     return coord_phys
diff --git a/test/test_coordinates.py b/test/test_coordinates.py
new file mode 100644
index 00000000..d8d66507
--- /dev/null
+++ b/test/test_coordinates.py
@@ -0,0 +1,35 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*
+
+import numpy as np
+
+from shimmingtoolbox.coils.coordinates import generate_meshgrid
+
+
+def test_generate_meshgrid():
+    """Test to verify generate_meshgrid outputs the correct scanner coordinates from input voxels"""
+
+    affine = np.array([[0., 0.,    3., -3.61445999],
+                      [-2.91667008, 0., 0., 101.76699829],
+                      [0., 2.91667008, 0., -129.85464478],
+                      [0., 0., 0., 1.]])
+
+    nx = 2
+    ny = 2
+    nz = 2
+    coord = generate_meshgrid((nx, ny, nz), affine)
+
+    expected = [np.array([[[-3.61445999, -0.61445999],
+                           [-3.61445999, -0.61445999]],
+                          [[-3.61445999, -0.61445999],
+                           [-3.61445999, -0.61445999]]]),
+                np.array([[[101.76699829, 101.76699829],
+                           [101.76699829, 101.76699829]],
+                          [[98.85032821,  98.85032821],
+                           [98.85032821,  98.85032821]]]),
+                np.array([[[-129.85464478, -129.85464478],
+                           [-126.9379747, -126.9379747]],
+                          [[-129.85464478, -129.85464478],
+                           [-126.9379747, -126.9379747]]])]
+
+    assert(np.all(np.isclose(coord, expected)))
diff --git a/test/test_siemens_basis.py b/test/test_siemens_basis.py
index 858fea16..19fc7310 100644
--- a/test/test_siemens_basis.py
+++ b/test/test_siemens_basis.py
@@ -47,8 +47,8 @@ def test_siemens_basis_resample():
     basis = siemens_basis(coord_phys[0], coord_phys[1], coord_phys[2])
 
     # Hard-coded values corresponding to the mid-point of the FOV.
-    expected = np.array([5.21405621e-18, -8.51520000e-02,  1.02182400e+00,  2.44386240e-02,
-                         2.50274698e-19, -4.08729600e-03, -1.70304000e-04, -2.08562248e-20])
+    expected = np.array([5.31991104e-18, -8.68807405e-02,  1.03212742e+00,  2.49321889e-02,
+                         2.57930573e-19, -4.21232592e-03, -1.77289155e-04, -2.17116596e-20])
 
     nx, ny, nz = nii.get_fdata().shape
     assert(np.all(np.isclose(basis[int(nx/2), int(ny/2), int(nz/2), :], expected, rtol=1e-05)))