Repository: mhe/pynrrd

Issue:
Confusing representation of "space direction"
Hi,
I've recently been working with MONAI that optionally uses pynrrd to read .nrrd and has an error when handling the space direction. 
The way the variable space direction is read by pynrrd (in the 3d case as a 3x3 matrix) makes it easy to misinterpret as the upper left 3x3 part of the affine matrix (this is what is done in MONAI). 
However, the matrix loaded by pynrrd is the transpose of the upper left 3x3 part of the affine matrix.

Loading the space direction variable to a matrix might cause this confusion, an alternative would be to instead load it as a list of vectors.

Discussion:

Patch:
diff --git a/docs/source/background/datatypes.rst b/docs/source/background/datatypes.rst
index 75e08b6..0b13dfa 100644
--- a/docs/source/background/datatypes.rst
+++ b/docs/source/background/datatypes.rst
@@ -88,4 +88,22 @@ double matrix
 :Python Datatype: (M,N) :class:`numpy.ndarray` of :class:`float`
 :Python Example: np.array([[2.54, 1.3, 0.0], [3.14, 0.3, 3.3], [np.nan, np.nan, np.nan], [0.0, -12.3, -3.3]])
 
-This datatype has the added feature where rows can be defined as empty by setting the vector as :code:`none`. In the NRRD specification, instead of the row, the :code:`none` keyword is used in it's place. This is represented in the Python NumPy array as a row of all NaN's. An example use case for this optional row matrix is for the 'space directions' field where one row may be empty because it is not a domain type.
\ No newline at end of file
+This datatype has the added feature where rows can be defined as empty by setting the vector as :code:`none`. In the NRRD specification, instead of the row, the :code:`none` keyword is used in it's place. This is represented in the Python NumPy array as a row of all NaN's. An example use case for this optional row matrix is for the 'space directions' field where one row may be empty because it is not a domain type.
+
+int vector list
+~~~~~~~~~~~~~~~~~~
+:NRRD Syntax: (<i>,<i>,...,<i>) (<i>,<i>,...,<i>) ... (<i>,<i>,...,<i>)
+:NRRD Example: (1,0,0) (0,1,0) none (0,0,1)
+:Python Datatype: (M,N) :class:`list` of (N,) :class:`numpy.ndarray` of :class:`int`
+:Python Example: [np.array([1, 0, 0]), np.array([0, 1, 0]), None, np.array([0, 0, 1])]
+
+This datatype is similar to `int matrix`_ except instead of returning a (M,N) :class:`numpy.ndarray`, it returns a list of (N,) :class:`numpy.ndarray`. Each row is optional and designated by :code:`none` in the NRRD specification and represented as :obj:`None` in this library.
+
+double vector list
+~~~~~~~~~~~~~~~~~~
+:NRRD Syntax: (<d>,<d>,...,<d>) (<d>,<d>,...,<d>) ... (<d>,<d>,...,<d>)
+:NRRD Example: (2.54, 1.3, 0.0) (3.14, 0.3, 3.3) none (0.05, -12.3, -3.3)
+:Python Datatype: (M,N) :class:`list` of (N,) :class:`numpy.ndarray` of :class:`float`
+:Python Example: [np.array([2.54, 1.3, 0.0]), np.array([3.14, 0.3, 3.3]), None, np.array([0.0, -12.3, -3.3])]
+
+This datatype is similar to `double matrix`_ except instead of returning a (M,N) :class:`numpy.ndarray`, it returns a list of (N,) :class:`numpy.ndarray`. Each row is optional and designated by :code:`none` in the NRRD specification and represented as :obj:`None` in this library.
\ No newline at end of file
diff --git a/docs/source/background/fields.rst b/docs/source/background/fields.rst
index c0ac23e..03209ee 100644
--- a/docs/source/background/fields.rst
+++ b/docs/source/background/fields.rst
@@ -33,7 +33,7 @@ centerings_               :ref:`background/datatypes:string list`
 space_                    :ref:`background/datatypes:string`
 `space dimension`_        :ref:`background/datatypes:int`
 `space units`_            :ref:`background/datatypes:quoted string list`
-`space directions`_       :ref:`background/datatypes:double matrix`
+`space directions`_       :ref:`background/datatypes:double matrix` or :ref:`background/datatypes:double vector list` depending on :data:`nrrd.SPACE_DIRECTIONS_TYPE`
 `space origin`_           :ref:`background/datatypes:double vector`
 `measurement frame`_      :ref:`background/datatypes:int matrix`
 ========================  ==============================================
diff --git a/docs/source/reference/formatting.rst b/docs/source/reference/formatting.rst
index 1e83557..259b1e9 100644
--- a/docs/source/reference/formatting.rst
+++ b/docs/source/reference/formatting.rst
@@ -9,8 +9,10 @@ Formatting NRRD fields
     nrrd.format_optional_vector
     nrrd.format_matrix
     nrrd.format_optional_matrix
+    nrrd.format_vector_list
+    nrrd.format_optional_vector_list
 
 .. automodule:: nrrd
-    :members: format_number, format_number_list, format_vector, format_optional_vector, format_matrix, format_optional_matrix
+    :members: format_number, format_number_list, format_vector, format_optional_vector, format_matrix, format_optional_matrix, format_vector_list, format_optional_vector_list
     :undoc-members:
     :show-inheritance:
diff --git a/docs/source/reference/parsing.rst b/docs/source/reference/parsing.rst
index e324cd6..c89ba5b 100644
--- a/docs/source/reference/parsing.rst
+++ b/docs/source/reference/parsing.rst
@@ -9,8 +9,10 @@ Parsing NRRD fields
     nrrd.parse_optional_vector
     nrrd.parse_matrix
     nrrd.parse_optional_matrix
+    nrrd.parse_vector_list
+    nrrd.parse_optional_vector_list
 
 .. automodule:: nrrd
-    :members: parse_number_auto_dtype, parse_number_list, parse_vector, parse_optional_vector, parse_matrix, parse_optional_matrix
+    :members: parse_number_auto_dtype, parse_number_list, parse_vector, parse_optional_vector, parse_matrix, parse_optional_matrix, parse_vector_list, parse_optional_vector_list
     :undoc-members:
     :show-inheritance:
diff --git a/docs/source/reference/reading.rst b/docs/source/reference/reading.rst
index 9aca987..277a1ab 100644
--- a/docs/source/reference/reading.rst
+++ b/docs/source/reference/reading.rst
@@ -7,6 +7,7 @@ Reading NRRD files
     nrrd.read_header
     nrrd.read_data
     nrrd.reader.ALLOW_DUPLICATE_FIELD
+    nrrd.SPACE_DIRECTIONS_TYPE
 
 .. automodule:: nrrd
     :members: read, read_header, read_data
@@ -14,3 +15,4 @@ Reading NRRD files
     :show-inheritance:
 
 .. autodata:: nrrd.reader.ALLOW_DUPLICATE_FIELD
+.. autodata:: nrrd.SPACE_DIRECTIONS_TYPE
diff --git a/nrrd/__init__.py b/nrrd/__init__.py
index 6bf8b3e..d7f8dcf 100644
--- a/nrrd/__init__.py
+++ b/nrrd/__init__.py
@@ -1,3 +1,5 @@
+from typing_extensions import Literal
+
 from nrrd._version import __version__
 from nrrd.formatters import *
 from nrrd.parsers import *
@@ -5,7 +7,40 @@ from nrrd.reader import read, read_data, read_header
 from nrrd.types import NRRDFieldMap, NRRDFieldType, NRRDHeader
 from nrrd.writer import write
 
+# TODO Change to 'double vector list' in next major release
+SPACE_DIRECTIONS_TYPE: Literal['double matrix', 'double vector list'] = 'double matrix'
+"""Datatype to use for 'space directions' field when reading/writing NRRD files
+
+The 'space directions' field can be represented in two different ways: as a matrix or as a list of vectors. Per the
+NRRD specification, the 'space directions' field is a per-axis definition that represents the direction and spacing of
+each axis. Non-spatial axes are represented as 'none'.
+
+The current default is to return a matrix, where each non-spatial axis is represented as a row of `NaN` in the matrix.
+In the next major release, this default option will change to return a list of optional vectors, where each non
+spatial axis is represented as `None`.
+
+Example:
+    Reading a NRRD file with space directions type set to 'double matrix' (the default).
+
+    >>> nrrd.SPACE_DIRECTIONS_TYPE = 'double matrix'
+    >>> data, header = nrrd.read('file.nrrd')
+    >>> print(header['space directions'])
+    [[1.5 0.  0. ]
+     [0.  1.5 0. ]
+     [0.  0.  1. ]
+     [nan nan nan]]
+
+    Reading a NRRD file with space directions type set to 'double vector list'.
+
+    >>> nrrd.SPACE_DIRECTIONS_TYPE = 'double vector list'
+    >>> data, header = nrrd.read('file.nrrd')
+    >>> print(header['space directions'])
+    [array([1.5, 0. , 0. ]), array([0. , 1.5, 0. ]), array([0., 0., 1.]), None]
+"""
+
 __all__ = ['read', 'read_data', 'read_header', 'write', 'format_number_list', 'format_number', 'format_matrix',
-           'format_optional_matrix', 'format_optional_vector', 'format_vector', 'parse_matrix',
-           'parse_number_auto_dtype', 'parse_number_list', 'parse_optional_matrix', 'parse_optional_vector',
-           'parse_vector', 'NRRDFieldType', 'NRRDFieldMap', 'NRRDHeader', '__version__']
+           'format_optional_matrix', 'format_optional_vector', 'format_vector', 'format_vector_list',
+           'format_optional_vector_list', 'parse_matrix', 'parse_number_auto_dtype', 'parse_number_list',
+           'parse_optional_matrix',
+           'parse_optional_vector', 'parse_vector', 'parse_vector_list', 'parse_optional_vector_list', 'NRRDFieldType',
+           'NRRDFieldMap', 'NRRDHeader', 'SPACE_DIRECTIONS_TYPE', '__version__']
diff --git a/nrrd/formatters.py b/nrrd/formatters.py
index 53a6d3d..3a1064d 100644
--- a/nrrd/formatters.py
+++ b/nrrd/formatters.py
@@ -1,4 +1,4 @@
-from typing import Optional, Union
+from typing import List, Optional, Union
 
 import numpy as np
 import numpy.typing as npt
@@ -57,6 +57,7 @@ def format_vector(x: npt.NDArray) -> str:
     vector : :class:`str`
         String containing NRRD vector
     """
+    x = np.asarray(x)
 
     return '(' + ','.join([format_number(y) for y in x]) + ')'
 
@@ -80,10 +81,15 @@ def format_optional_vector(x: Optional[npt.NDArray]) -> str:
     vector : :class:`str`
         String containing NRRD vector
     """
+    # If vector is None, return none
+    if x is None:
+        return 'none'
+
+    x = np.asarray(x)
 
-    # If vector is None or all elements are NaN, then return none
+    # If all elements are None or NaN, then return none
     # Otherwise format the vector as normal
-    if x is None or np.all(np.isnan(x)):
+    if np.all(x == None) or np.all(np.isnan(x)):  # noqa: E711
         return 'none'
     else:
         return format_vector(x)
@@ -131,6 +137,8 @@ def format_optional_matrix(x: Optional[npt.NDArray]) -> str:
     matrix : :class:`str`
         String containing NRRD matrix
     """
+    # Convert to float dtype to convert None to NaN
+    x = np.asarray(x, dtype=float)
 
     return ' '.join([format_optional_vector(y) for y in x])
 
@@ -151,5 +159,49 @@ def format_number_list(x: npt.NDArray) -> str:
     list : :class:`str`
         String containing NRRD list
     """
+    x = np.asarray(x)
 
     return ' '.join([format_number(y) for y in x])
+
+
+def format_vector_list(x: List[npt.NDArray]) -> str:
+    """Format a :class:`list` of (N,) :class:`numpy.ndarray` into a NRRD vector list string
+
+    See :ref:`background/datatypes:int vector list` and :ref:`background/datatypes:double vector list` for more
+    information on the format.
+
+    Parameters
+    ----------
+    x : :class:`list` of (N,) :class:`numpy.ndarray`
+        Vector list to convert to NRRD vector list string
+
+    Returns
+    -------
+    vector_list : :class:`str`
+        String containing NRRD vector list
+    """
+
+    return ' '.join([format_vector(y) for y in x])
+
+
+def format_optional_vector_list(x: List[Optional[npt.NDArray]]) -> str:
+    """Format a :class:`list` of (N,) :class:`numpy.ndarray` or :obj:`None` into a NRRD optional vector list string
+
+    Function converts a :class:`list` of (N,) :class:`numpy.ndarray` or :obj:`None` into a string using
+    the NRRD vector list format.
+
+    See :ref:`background/datatypes:int vector list` and :ref:`background/datatypes:double vector list` for more
+    information on the format.
+
+    Parameters
+    ----------
+    x : :class:`list` of (N,) :class:`numpy.ndarray` or :obj:`None`
+        Vector list to convert to NRRD vector list string
+
+    Returns
+    -------
+    vector_list : :class:`str`
+        String containing NRRD vector list
+    """
+
+    return ' '.join([format_optional_vector(y) for y in x])
diff --git a/nrrd/parsers.py b/nrrd/parsers.py
index e164bb3..01d654d 100644
--- a/nrrd/parsers.py
+++ b/nrrd/parsers.py
@@ -1,4 +1,4 @@
-from typing import Optional, Type, Union
+from typing import List, Optional, Type, Union
 
 import numpy as np
 import numpy.typing as npt
@@ -212,6 +212,103 @@ def parse_number_list(x: str, dtype: Optional[Type[Union[int, float]]] = None) -
     return number_list
 
 
+def parse_vector_list(x: str, dtype: Optional[Type[Union[int, float]]] = None) -> List[npt.NDArray]:
+    """Parse NRRD vector list from string into a :class:`list` of (N,) :class:`numpy.ndarray`.
+
+    Parses input string to convert it into a list of Numpy arrays using the NRRD vector list format.
+
+    See :ref:`background/datatypes:int vector list` and :ref:`background/datatypes:double vector list` for more
+    information on the format.
+
+    Parameters
+    ----------
+    x : :class:`str`
+        String containing NRRD vector list
+    dtype : data-type, optional
+        Datatype to use for the resulting Numpy arrays. Datatype can be :class:`float`, :class:`int` or :obj:`None`. If
+        :obj:`dtype` is :obj:`None`, it will be automatically determined by checking any of the vector elements
+        for fractional numbers. If found, the vectors will be converted to :class:`float`, otherwise :class:`int`.
+        Default is to automatically determine datatype.
+
+    Returns
+    -------
+    vector_list : :class:`list` of (N,) :class:`numpy.ndarray`
+        List of vectors that are parsed from the :obj:`x` string
+    """
+
+    # Split input by spaces, convert each row into a vector
+    vector_list = [parse_vector(x, dtype=float) for x in x.split()]
+
+    # Get the size of each row vector and then remove duplicate sizes
+    # There should be exactly one value in the matrix because all row sizes need to be the same
+    if len(np.unique([len(x) for x in vector_list])) != 1:
+        raise NRRDError('Vector list should have same number of elements in each row')
+
+    # If using automatic datatype detection, then start by converting to float and determining if the number is whole
+    # Truncate to integer if dtype is int also
+    if dtype is None:
+        vector_list_trunc = [x.astype(int) for x in vector_list]
+        if np.all([np.array_equal(x, y) for x, y in zip(vector_list, vector_list_trunc)]):
+            vector_list = vector_list_trunc
+    elif dtype == int:
+        vector_list = [x.astype(int) for x in vector_list]
+    elif dtype != float:
+        raise NRRDError('dtype should be None for automatic type detection, float or int')
+
+    return vector_list
+
+
+def parse_optional_vector_list(x: str, dtype: Optional[Type[Union[int, float]]] = None) -> List[Optional[npt.NDArray]]:
+    """Parse optional NRRD vector list from string into :class:`list` of (N,) :class:`numpy.ndarray` of :class:`float`.
+
+    Function parses optional NRRD vector list from string into a list of (N,) :class:`numpy.ndarray` or :obj:`None`.
+    This function works the same as :meth:`parse_vector_list` except if a row vector in the list is none, the resulting
+    row in the returned list will be :obj:`None`.
+
+    See :ref:`background/datatypes:int vector list` and :ref:`background/datatypes:double vector list` for more
+    information on the format.
+
+    Parameters
+    ----------
+    x : :class:`str`
+        String containing NRRD vector list
+
+    Returns
+    -------
+    vector_list : :class:`list` of (N,) :class:`numpy.ndarray` or :obj:`None`
+        List of vectors that is parsed from the :obj:`x` string
+    """
+
+    # Split input by spaces to get each row and convert into a vector. The row can be 'none', in which case it will
+    # return None
+    vector_list = [parse_optional_vector(x, dtype=float) for x in x.split()]
+
+    # Get the size of each row vector, 0 if None
+    sizes = np.array([0 if x is None else len(x) for x in vector_list])
+
+    # Get sizes of each row vector removing duplicate sizes
+    # Since each row vector should be same size, the unique sizes should return one value for the row size or it may
+    # return a second one (0) if there are None vectors
+    unique_sizes = np.unique(sizes)
+
+    if len(unique_sizes) != 1 and (len(unique_sizes) != 2 or unique_sizes.min() != 0):
+        raise NRRDError('Vector list should have same number of elements in each row')
+
+    # If using automatic datatype detection, then start by converting to float and determining if the number is whole
+    # Truncate to integer if dtype is int also
+    if dtype is None:
+        vector_list_trunc = [x.astype(int) if x is not None else None for x in vector_list]
+
+        if np.all([np.array_equal(x, y) for x, y in zip(vector_list, vector_list_trunc)]):
+            vector_list = vector_list_trunc
+    elif dtype == int:
+        vector_list = [x.astype(int) if x is not None else None for x in vector_list]
+    elif dtype != float:
+        raise NRRDError('dtype should be None for automatic type detection, float or int')
+
+    return vector_list
+
+
 def parse_number_auto_dtype(x: str) -> Union[int, float]:
     """Parse number from string with automatic type detection.
 
diff --git a/nrrd/reader.py b/nrrd/reader.py
index b6a69b9..1dc757e 100644
--- a/nrrd/reader.py
+++ b/nrrd/reader.py
@@ -8,6 +8,7 @@ import zlib
 from collections import OrderedDict
 from typing import IO, Any, AnyStr, Iterable, Tuple
 
+import nrrd
 from nrrd.parsers import *
 from nrrd.types import IndexOrder, NRRDFieldMap, NRRDFieldType, NRRDHeader
 
@@ -19,7 +20,7 @@ _READ_CHUNKSIZE: int = 2 ** 32
 
 _NRRD_REQUIRED_FIELDS = ['dimension', 'type', 'encoding', 'sizes']
 
-ALLOW_DUPLICATE_FIELD = False
+ALLOW_DUPLICATE_FIELD: bool = False
 """Allow duplicate header fields when reading NRRD files
 
 When there are duplicated fields in a NRRD file header, pynrrd throws an error by default. Setting this field as
@@ -109,7 +110,7 @@ def _get_field_type(field: str, custom_field_map: Optional[NRRDFieldMap]) -> NRR
     elif field in ['measurement frame']:
         return 'double matrix'
     elif field in ['space directions']:
-        return 'double matrix'
+        return nrrd.SPACE_DIRECTIONS_TYPE
     else:
         if custom_field_map and field in custom_field_map:
             return custom_field_map[field]
@@ -144,6 +145,10 @@ def _parse_field_value(value: str, field_type: NRRDFieldType) -> Any:
         # This is only valid for double matrices because the matrix is represented with NaN in the entire row
         # for none rows. NaN is only valid for floating point numbers
         return parse_optional_matrix(value)
+    elif field_type == 'int vector list':
+        return parse_optional_vector_list(value, dtype=int)
+    elif field_type == 'double vector list':
+        return parse_optional_vector_list(value, dtype=float)
     else:
         raise NRRDError(f'Invalid field type given: {field_type}')
 
diff --git a/nrrd/writer.py b/nrrd/writer.py
index e439409..52da6d4 100644
--- a/nrrd/writer.py
+++ b/nrrd/writer.py
@@ -97,6 +97,10 @@ def _format_field_value(value: Any, field_type: NRRDFieldType) -> str:
         return format_matrix(value)
     elif field_type == 'double matrix':
         return format_optional_matrix(value)
+    elif field_type == 'int vector list':
+        return format_optional_vector_list(value)
+    elif field_type == 'double vector list':
+        return format_optional_vector_list(value)
     else:
         raise NRRDError(f'Invalid field type given: {field_type}')
 


Test patch:
diff --git a/nrrd/tests/test_formatting.py b/nrrd/tests/test_formatting.py
index de4461f..97dad70 100644
--- a/nrrd/tests/test_formatting.py
+++ b/nrrd/tests/test_formatting.py
@@ -44,6 +44,8 @@ class TestFieldFormatting(unittest.TestCase):
         self.assertEqual(nrrd.format_optional_vector(None), 'none')
         self.assertEqual(nrrd.format_optional_vector(np.array([np.nan, np.nan, np.nan])), 'none')
         self.assertEqual(nrrd.format_optional_vector([np.nan, np.nan, np.nan]), 'none')
+        self.assertEqual(nrrd.format_optional_vector([None, None, None]), 'none')
+        self.assertEqual(nrrd.format_optional_vector(np.array([None, None, None])), 'none')
 
     def test_format_matrix(self):
         self.assertEqual(nrrd.format_matrix(np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])), '(1,2,3) (4,5,6) (7,8,9)')
@@ -82,6 +84,19 @@ class TestFieldFormatting(unittest.TestCase):
         self.assertEqual(nrrd.format_optional_matrix(np.array([
             [1, 2, 3], [np.nan, np.nan, np.nan], [4, 5, 6], [7, 8, 9]])),
             '(1,2,3) none (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_matrix(np.array([
+            [None, None, None], [1, 2, 3], [4, 5, 6], [7, 8, 9]])),
+            'none (1,2,3) (4,5,6) (7,8,9)')
+
+        self.assertEqual(nrrd.format_optional_matrix([
+            [np.nan, np.nan, np.nan], [1, 2, 3], [4, 5, 6], [7, 8, 9]]),
+            'none (1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_matrix([
+            [1, 2, 3], [np.nan, np.nan, np.nan], [4, 5, 6], [7, 8, 9]]),
+            '(1,2,3) none (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_matrix([
+            [None, None, None], [1, 2, 3], [4, 5, 6], [7, 8, 9]]),
+            'none (1,2,3) (4,5,6) (7,8,9)')
 
     def test_format_number_list(self):
         self.assertEqual(nrrd.format_number_list([1, 2, 3]), '1 2 3')
@@ -92,6 +107,49 @@ class TestFieldFormatting(unittest.TestCase):
         self.assertEqual(nrrd.format_number_list(np.array([1., 2., 3.])), '1 2 3')
         self.assertEqual(nrrd.format_number_list(np.array([1.2, 2., 3.2])), '1.2 2 3.2000000000000002')
 
+    def test_format_vector_list(self):
+        self.assertEqual(nrrd.format_vector_list([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
+                         '(1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_vector_list([[1., 2., 3.], [4., 5., 6.], [7., 8., 9.]]),
+                         '(1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_vector_list([[1, 2.2, 3.3], [4.4, 5.5, 6.6], [7.7, 8.8, 9.9]]),
+                         '(1,2.2000000000000002,3.2999999999999998) (4.4000000000000004,5.5,6.5999999999999996) '
+                         '(7.7000000000000002,8.8000000000000007,9.9000000000000004)')
+
+        self.assertEqual(nrrd.format_vector_list(np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])),
+                         '(1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_vector_list(np.array([[1., 2., 3.], [4., 5., 6.], [7., 8., 9.]])),
+                         '(1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_vector_list(np.array([[1, 2.2, 3.3], [4.4, 5.5, 6.6], [7.7, 8.8, 9.9]])),
+                         '(1,2.2000000000000002,3.2999999999999998) (4.4000000000000004,5.5,6.5999999999999996) '
+                         '(7.7000000000000002,8.8000000000000007,9.9000000000000004)')
+
+    def test_format_optional_vector_list(self):
+        self.assertEqual(nrrd.format_optional_vector_list([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
+                         '(1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_vector_list([[1., 2., 3.], [4., 5., 6.], [7., 8., 9.]]),
+                         '(1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_vector_list([[1, 2.2, 3.3], [4.4, 5.5, 6.6], [7.7, 8.8, 9.9]]),
+                         '(1,2.2000000000000002,3.2999999999999998) (4.4000000000000004,5.5,6.5999999999999996) '
+                         '(7.7000000000000002,8.8000000000000007,9.9000000000000004)')
+
+        self.assertEqual(nrrd.format_optional_vector_list([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
+                         '(1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_vector_list([[1., 2., 3.], [4., 5., 6.], [7., 8., 9.]]),
+                         '(1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_vector_list([[1, 2.2, 3.3], [4.4, 5.5, 6.6], [7.7, 8.8, 9.9]]),
+                         '(1,2.2000000000000002,3.2999999999999998) (4.4000000000000004,5.5,6.5999999999999996) '
+                         '(7.7000000000000002,8.8000000000000007,9.9000000000000004)')
+
+        self.assertEqual(nrrd.format_optional_vector_list([[np.nan, np.nan, np.nan], [1, 2, 3], [4, 5, 6], [7, 8, 9]]),
+                         'none (1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_vector_list([[1, 2, 3], [np.nan, np.nan, np.nan], [4, 5, 6], [7, 8, 9]]),
+                         '(1,2,3) none (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_vector_list([[None, None, None], [1, 2, 3], [4, 5, 6], [7, 8, 9]]),
+                         'none (1,2,3) (4,5,6) (7,8,9)')
+        self.assertEqual(nrrd.format_optional_vector_list([None, [1, 2, 3], [4, 5, 6], [7, 8, 9]]),
+                         'none (1,2,3) (4,5,6) (7,8,9)')
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/nrrd/tests/test_parsing.py b/nrrd/tests/test_parsing.py
index 2897bc7..ef53f93 100644
--- a/nrrd/tests/test_parsing.py
+++ b/nrrd/tests/test_parsing.py
@@ -10,7 +10,10 @@ class TestFieldParsing(unittest.TestCase):
         pass
 
     def assert_equal_with_datatype(self, desired, actual):
-        self.assertEqual(desired.dtype, np.array(actual[0]).dtype)
+        if isinstance(desired, list):
+            self.assertEqual(desired[0].dtype, np.array(actual[0]).dtype)
+        else:
+            self.assertEqual(desired.dtype, np.array(actual[0]).dtype)
         np.testing.assert_equal(desired, actual)
 
     def test_parse_vector(self):
@@ -154,6 +157,64 @@ class TestFieldParsing(unittest.TestCase):
         self.assertEqual(nrrd.parse_number_auto_dtype('25'), 25)
         self.assertEqual(nrrd.parse_number_auto_dtype('25.125'), 25.125)
 
+    def test_parse_vector_list(self):
+        self.assert_equal_with_datatype(
+            nrrd.parse_vector_list('(1.4726600000000003,-0,0) (-0,1.4726600000000003,-0) (0,-0,4.7619115092114601)'),
+            [[1.4726600000000003, 0, 0], [0, 1.4726600000000003, 0], [0, 0, 4.7619115092114601]])
+
+        self.assert_equal_with_datatype(
+            nrrd.parse_vector_list('(1.4726600000000003,-0,0) (-0,1.4726600000000003,-0) (0,-0,4.7619115092114601)',
+                                   dtype=float),
+            [[1.4726600000000003, 0, 0], [0, 1.4726600000000003, 0], [0, 0, 4.7619115092114601]])
+
+        self.assert_equal_with_datatype(
+            nrrd.parse_vector_list('(1.4726600000000003,-0,0) (-0,1.4726600000000003,-0) (0,-0,4.7619115092114601)',
+                                   dtype=int), [[1, 0, 0], [0, 1, 0], [0, 0, 4]])
+
+        self.assert_equal_with_datatype(nrrd.parse_vector_list('(1,0,0) (0,1,0) (0,0,1)'),
+                                        [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+        self.assert_equal_with_datatype(nrrd.parse_vector_list('(1,0,0) (0,1,0) (0,0,1)', dtype=float),
+                                        [[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
+        self.assert_equal_with_datatype(nrrd.parse_vector_list('(1,0,0) (0,1,0) (0,0,1)', dtype=int),
+                                        [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+
+        with self.assertRaisesRegex(nrrd.NRRDError, 'Vector list should have same number of elements in each row'):
+            nrrd.parse_vector_list('(1,0,0,0) (0,1,0) (0,0,1)')
+
+        with self.assertRaisesRegex(nrrd.NRRDError, 'dtype should be None for automatic type detection, float or int'):
+            nrrd.parse_vector_list('(1,0,0) (0,1,0) (0,0,1)', dtype=np.uint8)
+
+        vector_list = nrrd.parse_vector_list('(1,0,0) (0,1,0) (0,0,1)')
+        self.assertIsInstance(vector_list, list)
+        self.assertTrue(all(isinstance(vector, np.ndarray) for vector in vector_list))
+
+    def test_parse_optional_vector_list(self):
+        self.assert_equal_with_datatype(nrrd.parse_optional_vector_list(
+            '(1.4726600000000003,-0,0) (-0,1.4726600000000003,-0) (0,-0,4.7619115092114601)'),
+            [[1.4726600000000003, 0, 0], [0, 1.4726600000000003, 0], [0, 0, 4.7619115092114601]])
+
+        self.assert_equal_with_datatype(nrrd.parse_optional_vector_list('(1,0,0) (0,1,0) (0,0,1)'),
+                                        [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+
+        self.assert_equal_with_datatype(nrrd.parse_optional_vector_list(
+            '(1.4726600000000003,-0,0) none (-0,1.4726600000000003,-0) (0,-0,4.7619115092114601)'),
+            [[1.4726600000000003, 0, 0], None, [0, 1.4726600000000003, 0],
+             [0, 0, 4.7619115092114601]])
+
+        self.assert_equal_with_datatype(nrrd.parse_optional_vector_list(
+            '(1.4726600000000003,-0,0) none (-0,1.4726600000000003,-0) (0,-0,4.7619115092114601)'),
+            [[1.4726600000000003, 0, 0], None, [0, 1.4726600000000003, 0], [0, 0, 4.7619115092114601]])
+
+        with self.assertRaisesRegex(nrrd.NRRDError, 'Vector list should have same number of elements in each row'):
+            nrrd.parse_optional_vector_list('(1,0,0,0) (0,1,0) (0,0,1)')
+
+        with self.assertRaisesRegex(nrrd.NRRDError, 'Vector list should have same number of elements in each row'):
+            nrrd.parse_optional_vector_list('none (1,0,0,0) (0,1,0) (0,0,1)')
+
+        vector_list = nrrd.parse_optional_vector_list('(1,0,0) (0,1,0) none (0,0,1)')
+        self.assertIsInstance(vector_list, list)
+        self.assertTrue(all(vector is None or isinstance(vector, np.ndarray) for vector in vector_list))
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/nrrd/tests/test_reading.py b/nrrd/tests/test_reading.py
index c0a6b68..ba3b9da 100644
--- a/nrrd/tests/test_reading.py
+++ b/nrrd/tests/test_reading.py
@@ -501,6 +501,22 @@ class Abstract:
                 with self.subTest(filename):
                     test(filename)
 
+        def test_read_space_directions_list(self):
+            try:
+                nrrd.SPACE_DIRECTIONS_TYPE = 'double vector list'
+
+                _, header = nrrd.read(RAW_4D_NRRD_FILE_PATH, index_order=self.index_order)
+                self.assertIsInstance(header['space directions'], list)
+                self.assertTrue(
+                    all(vector is None or isinstance(vector, np.ndarray) for vector in header['space directions']))
+                np.testing.assert_equal(header['space directions'][0].dtype, np.float64)
+                np.testing.assert_equal(header['space directions'], [np.array([1.5, 0., 0.]),
+                                        np.array([0., 1.5, 0.]),
+                                        np.array([0., 0., 1.]),
+                                        None])
+            finally:
+                nrrd.SPACE_DIRECTIONS_TYPE = 'double matrix'
+
 
 class TestReadingFunctionsFortran(Abstract.TestReadingFunctions):
     index_order = 'F'
diff --git a/nrrd/tests/test_writing.py b/nrrd/tests/test_writing.py
index 9b3e842..ce7bb2e 100644
--- a/nrrd/tests/test_writing.py
+++ b/nrrd/tests/test_writing.py
@@ -478,6 +478,50 @@ class Abstract:
             self.assertEqual(header.pop('sizes').all(), memory_header.pop('sizes').all())
             self.assertSequenceEqual(header, memory_header)
 
+        def test_space_directions_header(self):
+            space_directions_list = [np.array([1.5, 0., 0.]),
+                                     np.array([0., 1.5, 0.]),
+                                     np.array([0., 0., 1.]),
+                                     None]
+            space_directions_matrix = np.array([[1.5, 0., 0.],
+                                               [0., 1.5, 0.],
+                                               [0., 0., 1.],
+                                               [np.nan, np.nan, np.nan]])
+
+            with self.subTest('double matrix -> double matrix'):
+                nrrd.SPACE_DIRECTIONS_TYPE = 'double matrix'
+                test_output_filename = os.path.join(self.temp_write_dir, 'testfile_space_directions_matrix.nrrd')
+                nrrd.write(test_output_filename, self.data_input, {'space directions': space_directions_matrix},
+                           index_order=self.index_order)
+                _, header = nrrd.read(test_output_filename, index_order=self.index_order)
+                np.testing.assert_equal(header['space directions'], space_directions_matrix)
+
+            with self.subTest('double matrix -> double vector list'):
+                nrrd.SPACE_DIRECTIONS_TYPE = 'double matrix'
+                test_output_filename = os.path.join(self.temp_write_dir, 'testfile_space_directions_matrix.nrrd')
+                nrrd.write(test_output_filename, self.data_input, {'space directions': space_directions_matrix},
+                           index_order=self.index_order)
+                nrrd.SPACE_DIRECTIONS_TYPE = 'double vector list'
+                _, header = nrrd.read(test_output_filename, index_order=self.index_order)
+                np.testing.assert_equal(header['space directions'], space_directions_list)
+
+            with self.subTest('double vector list'):
+                nrrd.SPACE_DIRECTIONS_TYPE = 'double vector list'
+                test_output_filename = os.path.join(self.temp_write_dir, 'testfile_space_directions_list.nrrd')
+                nrrd.write(test_output_filename, self.data_input, {'space directions': space_directions_list},
+                           index_order=self.index_order)
+                _, header = nrrd.read(test_output_filename, index_order=self.index_order)
+                np.testing.assert_equal(header['space directions'], space_directions_list)
+
+            with self.subTest('double vector list -> double matrix'):
+                nrrd.SPACE_DIRECTIONS_TYPE = 'double vector list'
+                test_output_filename = os.path.join(self.temp_write_dir, 'testfile_space_directions_list.nrrd')
+                nrrd.write(test_output_filename, self.data_input, {'space directions': space_directions_list},
+                           index_order=self.index_order)
+                nrrd.SPACE_DIRECTIONS_TYPE = 'double matrix'
+                _, header = nrrd.read(test_output_filename, index_order=self.index_order)
+                np.testing.assert_equal(header['space directions'], space_directions_matrix)
+
 
 class TestWritingFunctionsFortran(Abstract.TestWritingFunctions):
     index_order = 'F'