Repository: dicompyler/dicompyler-core

Issue:
Decubitus orientation
I'm hoping for some guidance on DVHs in decubitus situation - I have CT and RTDOSE files with ImageOrientationPatient of [0, -1, 0, 1, 0, 0] and I get empty DVHs returned.

On calling `get_dvh` with debugging on, I see many lines like:
`Dose plane not found for -100.00. Using -140.25 to calculate contour volume`.

I see that there have been some issues and PR related to non-HFS orientation (e.g. #183), and in reviewing the code, I don't see any particular limitations for decubitus.  I'm using the latest commit from master.

I haven't pursued it in further detail hoping to get some pointers in the right direction for further investigation.  For example I haven't tried changing any of the optional parameters to `get_dvh`.

Thanks



Discussion:
pep8speaks: Hello @darcymason! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:

* In the file [`dicompylercore/dicomparser.py`](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/dicompylercore/dicomparser.py):

> [Line 59:27](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/dicompylercore/dicomparser.py#L59): [E272](https://duckduckgo.com/?q=pep8%20E272) multiple spaces before keyword
> [Line 94:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/dicompylercore/dicomparser.py#L94): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (86 > 79 characters)
> [Line 592:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/dicompylercore/dicomparser.py#L592): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (80 > 79 characters)



* In the file [`dicompylercore/dvhcalc.py`](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/dicompylercore/dvhcalc.py):

> [Line 239:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/dicompylercore/dvhcalc.py#L239): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (89 > 79 characters)
> [Line 497:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/dicompylercore/dvhcalc.py#L497): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (83 > 79 characters)
> [Line 498:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/dicompylercore/dvhcalc.py#L498): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (83 > 79 characters)

* In the file [`tests/test_decubitus.py`](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py):

> [Line 25:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L25): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (88 > 79 characters)
> [Line 94:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L94): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (83 > 79 characters)
> [Line 97:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L97): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (88 > 79 characters)
> [Line 107:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L107): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (80 > 79 characters)
> [Line 244:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L244): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (80 > 79 characters)
> [Line 249:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L249): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (81 > 79 characters)
> [Line 284:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L284): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (82 > 79 characters)
> [Line 285:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L285): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (80 > 79 characters)
> [Line 326:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L326): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (80 > 79 characters)
> [Line 328:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L328): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (82 > 79 characters)
> [Line 339:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L339): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (87 > 79 characters)
> [Line 343:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L343): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (80 > 79 characters)
> [Line 345:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L345): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (82 > 79 characters)
> [Line 366:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L366): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (83 > 79 characters)
> [Line 421:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L421): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (82 > 79 characters)
> [Line 458:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L458): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (83 > 79 characters)
> [Line 469:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L469): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (89 > 79 characters)
> [Line 472:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L472): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (82 > 79 characters)
> [Line 476:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L476): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (83 > 79 characters)
> [Line 496:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L496): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (84 > 79 characters)
> [Line 500:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L500): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (81 > 79 characters)
> [Line 552:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L552): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (84 > 79 characters)
> [Line 555:80](https://github.com/dicompyler/dicompyler-core/blob/b32d78df1406e5e4e3188f1e8ad6e801c33bcfaa/tests/test_decubitus.py#L555): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (90 > 79 characters)




bastula: @darcymason Thanks for the comprehensive PR. Is this ready for review? Do you have any more commits to add?
darcymason: > Good PR! I think there are a few small discussion points, but if we can resolve/discuss them, it should be good to merge in.

Thanks for the review - I've had a quick look and can address these points, might just take me a few days due to other commitments.
lgtm-com[bot]: This pull request **introduces 1 alert** when merging 023090ab737f506c7548aedd57ce7cda2fb72408 into 24d9a9bb28a04280a9c4c8d2a0b6e051810eb7cf - [view on LGTM.com](https://lgtm.com/projects/g/dicompyler/dicompyler-core/rev/pr-555884c75ade14df118de7e74b9d69c0a9b504a8)

**new alerts:**

* 1 for Unused local variable
darcymason: @bastula, sorry for the delay in addressing this review - I think everything is covered now, please have a look and let me know what you think.  

Still a few codacy style issues, for me the lines complained about look standard, and my local style checkers pass (pycodestyle and pydocstyle).
bastula: Thanks @darcymason for addressing the changes. I'll try to take a look in the next day or two and see if we can get this merged in.
darcymason: Hi @bastula, can I make a friendly nudge on this?  I would like this in for a deployment that points to master.  Thanks.

Patch:
diff --git a/dicompylercore/dicomparser.py b/dicompylercore/dicomparser.py
index e204f97..39f9d62 100755
--- a/dicompylercore/dicomparser.py
+++ b/dicompylercore/dicomparser.py
@@ -487,27 +487,37 @@ class DicomParser:
 
         Referenced matrix can be found in Part 3 Section C.7.6.2.1.1
         """
-        di = self.ds.PixelSpacing[0]
-        dj = self.ds.PixelSpacing[1]
-        orientation = self.ds.ImageOrientationPatient
-        position = self.ds.ImagePositionPatient
+        drow, dcol = self.ds.PixelSpacing
 
-        m = np.matrix(
-            [[orientation[0]*di, orientation[3]*dj, 0, position[0]],
-             [orientation[1]*di, orientation[4]*dj, 0, position[1]],
-             [orientation[2]*di, orientation[5]*dj, 0, position[2]],
+        orientation = self.ds.ImageOrientationPatient
+        first_x, first_y, first_z = self.ds.ImagePositionPatient
+        num_cols = self.ds.Columns
+        num_rows = self.ds.Rows
+
+        # Determine which way X and Y real-world coords run
+        # X runs across columns if x_lut_index is 0
+        # limits to head-first/feet-first and prone/supine/decubitus
+        x_index = self.x_lut_index()
+
+        m = np.array(
+            [[orientation[0]*dcol, orientation[3]*drow, 0, first_x],
+             [orientation[1]*dcol, orientation[4]*drow, 0, first_y],
+             [orientation[2]*dcol, orientation[5]*drow, 0, first_z],
              [0, 0, 0, 1]])
 
-        x = []
-        y = []
-        for i in range(0, self.ds.Columns):
-            imat = m * np.matrix([[i], [0], [0], [1]])
-            x.append(float(imat[0]))
-        for j in range(0, self.ds.Rows):
-            jmat = m * np.matrix([[0], [j], [0], [1]])
-            y.append(float(jmat[1]))
+        last_xy = np.matmul(
+            m, np.array([[num_cols-1], [num_rows-1], [0], [1]])
+        )
+        last_x, last_y = last_xy[0][0], last_xy[1][0]
 
-        return (np.array(x), np.array(y))
+        if x_index == 0:
+            col_lut = np.linspace(first_x, last_x, num_cols)
+            row_lut = np.linspace(first_y, last_y, num_rows)
+        else:
+            col_lut = np.linspace(first_y, last_y, num_cols)
+            row_lut = np.linspace(first_x, last_x, num_rows)
+
+        return col_lut, row_lut
 
 # ========================= RT Structure Set Methods =========================
 
@@ -784,11 +794,13 @@ class DicomParser:
             z = float(z)
             # Get the initial dose grid position (z) in patient coordinates
             ipp = self.ds.ImagePositionPatient
-            iop = self.ds.ImageOrientationPatient
             gfov = self.ds.GridFrameOffsetVector
             # Add the position to the offset vector to determine the
             # z coordinate of each dose plane
-            planes = (iop[0] * iop[4] * np.array(gfov)) + ipp[2]
+
+            z_sign = 1 if self.is_head_first_orientation() else -1
+
+            planes = (z_sign * np.array(gfov)) + ipp[2]
             frame = -1
             # Check to see if the requested plane exists in the array
             if (np.amin(np.fabs(planes - z)) < threshold):
@@ -841,6 +853,89 @@ class DicomParser:
 
         return doseplane
 
+    def is_head_first_orientation(self):
+        """Return True if self.orientation is head-first.
+
+        Raises
+        ------
+        NotImplementedError
+            Raised if orientation is not one of head/feet-first
+            and supine/prone/decubitus
+
+        Returns
+        -------
+        bool
+            True if orientation is head-first, else False
+
+
+        """
+        orientation = self.ds.ImageOrientationPatient
+        if any(
+            all(np.isclose(orientation, hf_orientation))
+            for hf_orientation in (  # noqa
+                [1,  0,  0,  0,  1,  0],  # Head First Supine
+                [-1,  0,  0,  0, -1,  0],  # Head First Prone
+                [0, -1,  0,  1,  0,  0],  # Head First Decubitus Left
+                [0,  1,  0, -1,  0,  0]   # Head First Decubitus Right
+            )
+        ):
+            return True
+        elif any(
+            all(np.isclose(orientation, ff_orientation))
+            for ff_orientation in (
+                [0,  1,  0,  1,  0,  0],  # Feet First Decubitus Left
+                [0, -1,  0, -1,  0,  0],  # Feet First Decubitus Right
+                [1,  0,  0,  0, -1,  0],  # Feet First Prone
+                [-1,  0,  0,  0,  1,  0]   # Feet First Supine
+            )
+        ):
+            return False
+        else:
+            raise NotImplementedError(
+                "Cannot calculate dose plane sign for non-standard orientation"
+            )
+
+    def x_lut_index(self):
+        """Return LUT index for real-world X direction.
+
+        Raises
+        ------
+        NotImplementedError
+            Raised if orientation is not one of head/feet-first
+            and supine/prone/decubitus
+
+        Returns
+        -------
+        X direction LUT index, matching 'lut' from GetDoseData
+            0 if real-world X across columns
+            1 if real-world X along rows
+        """
+        orientation = self.ds.ImageOrientationPatient
+        if any(
+            all(np.isclose(orientation, non_decub))
+            for non_decub in (
+                [1,  0,  0,  0,  1,  0],  # Head First Supine
+                [-1,  0,  0,  0, -1,  0],  # Head First Prone
+                [-1,  0,  0,  0,  1,  0],  # Feet First Supine
+                [1,  0,  0,  0, -1,  0]   # Feet First Prone
+            )
+        ):
+            return 0
+        elif any(
+            all(np.isclose(orientation, decub))
+            for decub in (
+                [0, -1,  0,  1,  0,  0],  # Head First Decubitus Left
+                [0,  1,  0, -1,  0,  0],  # Head First Decubitus Right
+                [0,  1,  0,  1,  0,  0],  # Feet First Decubitus Left
+                [0, -1,  0, -1,  0,  0]   # Feet First Decubitus Right
+            )
+        ):
+            return 1
+        else:
+            raise NotImplementedError(
+                "Cannot calculate X direction for non-standard orientation"
+            )
+
     def GetIsodosePoints(self, z=0, level=100, threshold=0.5):
         """Return dose grid points from an isodose level & slice position.
 
@@ -882,6 +977,7 @@ class DicomParser:
                 del pixel_array
         data['dosemax'] = float(dosemax)
         data['lut'] = self.GetPatientToPixelLUT()
+        data['x_lut_index'] = self.x_lut_index()
         data['fraction'] = ''
         if "ReferencedRTPlanSequence" in self.ds:
             plan = self.ds.ReferencedRTPlanSequence[0]
diff --git a/dicompylercore/dvh.py b/dicompylercore/dvh.py
index 14c2777..4d2786f 100644
--- a/dicompylercore/dvh.py
+++ b/dicompylercore/dvh.py
@@ -61,7 +61,7 @@ class DVH(object):
 
     @classmethod
     def from_dicom_dvh(cls, dataset, roi_num, rx_dose=None, name=None,
-            color=None):
+                       color=None):
         """Initialization for a DVH from a pydicom RT Dose DVH sequence."""
         sequence_num = -1
         for i, d in enumerate(dataset.DVHSequence):
diff --git a/dicompylercore/dvhcalc.py b/dicompylercore/dvhcalc.py
index 2004d1d..31a77c5 100755
--- a/dicompylercore/dvhcalc.py
+++ b/dicompylercore/dvhcalc.py
@@ -76,6 +76,7 @@ def get_dvh(structure,
         different formats using the attributes and properties of the DVH class.
     """
     from dicompylercore import dicomparser
+
     rtss = dicomparser.DicomParser(structure)
     rtdose = dicomparser.DicomParser(dose, memmap_pixel_array=memmap_rtdose)
     structures = rtss.GetStructures()
@@ -178,7 +179,10 @@ def _calculate_dvh(structure,
         # Code taken from Stack Overflow Answer from Joe Kington:
         # https://stackoverflow.com/q/3654289/74123
         # Create vertex coordinates for each grid cell
-        x, y = np.meshgrid(np.array(dd['lut'][0]), np.array(dd['lut'][1]))
+        x_index = dd['x_lut_index']
+        x, y = np.meshgrid(
+            np.array(dd['lut'][x_index]), np.array(dd['lut'][1-x_index])
+        )
         x, y = x.flatten(), y.flatten()
         dosegridpoints = np.vstack((x, y)).T
 
@@ -228,8 +232,18 @@ def _calculate_dvh(structure,
                 logger.warning('Dose plane not found for %s.' +
                                ' Using %s to calculate contour volume.',
                                z, origin_z)
+                # Create a dummy dose grid with the correct size.
+                # calculate_plane_histogram() and its methods provide the
+                # volume calc needed, but do so as part of DVH calc,
+                # which requires a dose grid
+                dummy_dose = dose.GetDoseGrid(origin_z)
+                if use_structure_extents:
+                    extents = dgindexextents
+                    dummy_dose = dummy_dose[
+                        extents[1]:extents[3], extents[0]:extents[2]
+                    ]
                 _, vol = calculate_plane_histogram(
-                    plane, dose.GetDoseGrid(origin_z), dosegridpoints, maxdose,
+                    plane, dummy_dose, dosegridpoints, maxdose,
                     dd, id, structure, hist)
                 planedata[z] = (np.array([0]), vol)
                 notes = 'Dose grid does not encompass every contour.' + \
@@ -288,7 +302,10 @@ def get_contour_mask(dd, id, dosegridpoints, contour):
     # return mask.reshape((len(doselut[1]), len(doselut[0])))
 
     grid = c.contains_points(dosegridpoints)
-    grid = grid.reshape((len(doselut[1]), len(doselut[0])))
+    if dd['x_lut_index'] == 0:  # X values across columns
+        grid = grid.reshape((len(doselut[1]), len(doselut[0])))
+    else:  # decubitus
+        grid = grid.reshape((len(doselut[0]), len(doselut[1]))).T
 
     return grid
 
@@ -341,7 +358,7 @@ def dosegrid_extents_indices(extents, dd, padding=1):
     Parameters
     ----------
     extents : list
-        Structure extents in patient coordintes: [xmin, ymin, xmax, ymax].
+        Structure extents in patient coordinates: [xmin, ymin, xmax, ymax].
         If an empty list, no structure extents will be used in the calculation.
     dd : dict
         Dose data from dicomparser.GetDoseData.
@@ -351,24 +368,41 @@ def dosegrid_extents_indices(extents, dd, padding=1):
     Returns
     -------
     list
-        Dose grid extents in pixel coordintes as array indices:
-        [xmin, ymin, xmax, ymax].
+        Dose grid extents in pixel coordinates as array indices:
+        [col_min, row_min, col_max, row_max].
     """
+    col_lut, row_lut = dd['lut']
+    num_cols = len(col_lut)
+    num_rows = len(row_lut)
     if not len(extents):
-        return [0, 0, dd['lut'][0].shape[0] - 1, dd['lut'][1].shape[0] - 1]
-    dgxmin = np.argmin(np.fabs(dd['lut'][0] - extents[0])) - padding
-    if dd['lut'][0][dgxmin] > extents[0]:
-        dgxmin -= 1
-    dgxmax = np.argmin(np.fabs(dd['lut'][0] - extents[2])) + padding
-    dgymin = np.argmin(np.fabs(dd['lut'][1] - extents[1])) - padding
-    dgymax = np.argmin(np.fabs(dd['lut'][1] - extents[3])) + padding
-    dgxmin = 0 if dgxmin < 0 else dgxmin
-    dgymin = 0 if dgymin < 0 else dgymin
-    if dgxmax == dd['lut'][0].shape[0]:
-        dgxmax = dd['lut'][0].shape[0] - 1
-    if dgymax == dd['lut'][1].shape[0]:
-        dgymax = dd['lut'][1].shape[0] - 1
-    return [dgxmin, dgymin, dgxmax, dgymax]
+        return [0, 0, num_cols - 1, num_rows - 1]
+
+    if dd['x_lut_index'] == 0:  # X is across rows
+        strx_col_min, strx_col_max = extents[0], extents[2]
+        strx_row_min, strx_row_max = extents[1], extents[3]
+    else:  # decubitus case, X down rows
+        strx_col_min, strx_col_max = extents[1], extents[3]
+        strx_row_min, strx_row_max = extents[0], extents[2]
+
+    dg_col_min = np.argmin(np.fabs(col_lut - strx_col_min))
+    dg_col_max = np.argmin(np.fabs(col_lut - strx_col_max))
+    if dg_col_min > dg_col_max:
+        dg_col_min, dg_col_max = dg_col_max, dg_col_min
+    if col_lut[dg_col_min] > strx_col_min:
+        dg_col_min -= 1
+
+    dg_row_min = np.argmin(np.fabs(row_lut - strx_row_min))
+    dg_row_max = np.argmin(np.fabs(row_lut - strx_row_max))
+    if dg_row_min > dg_row_max:
+        dg_row_min, dg_row_max = dg_row_max, dg_row_min
+
+    # Ensure indices within array limits regardless of padding
+    dg_col_min = max(0, dg_col_min-padding)
+    dg_row_min = max(0, dg_row_min-padding)
+    dg_col_max = min(num_cols - 1, dg_col_max+padding)
+    dg_row_max = min(num_rows - 1, dg_row_max+padding)
+
+    return [dg_col_min, dg_row_min, dg_col_max, dg_row_max]
 
 
 def dosegrid_extents_positions(extents, dd):
@@ -377,7 +411,8 @@ def dosegrid_extents_positions(extents, dd):
     Parameters
     ----------
     extents : list
-        Dose grid extents in pixel coordintes: [xmin, ymin, xmax, ymax].
+        Dose grid extents in pixel coordintes:
+        [col_pos_min, row_pos_min, col_pos_max, row_pos_max].
     dd : dict
         Dose data from dicomparser.GetDoseData.
 
@@ -403,7 +438,8 @@ def get_resampled_lut(index_extents,
     index_extents : list
         Dose grid extents as array indices.
     extents : list
-        Dose grid extents in patient coordinates.
+        Dose grid extents in patient coordinates:
+        [col_pos_min, row_pos_min, col_pos_max, row_pos_max].
     new_pixel_spacing : tuple or float
         New pixel spacing in mm (row, column).
         If float is provided, original dose grid pixel spacing must be square.
@@ -413,7 +449,8 @@ def get_resampled_lut(index_extents,
     Returns
     -------
     tuple
-        A tuple of lists (x, y) of patient to pixel coordinate mappings.
+        A tuple of lut lists (col lut, row lut) with the coordinates of the
+        dose grid in patient coordinates
 
     Raises
     ------
@@ -455,15 +492,21 @@ def get_resampled_lut(index_extents,
             % min_pixel_spacing[1] +
             " where n is an integer. Value provided was %s."
             % new_pixel_spacing[1])
-    sampling_rate = np.array([
-        abs(index_extents[0] - index_extents[2]),
-        abs(index_extents[1] - index_extents[3])
-    ])
-    xsamples = sampling_rate[0] * min_pixel_spacing[1] / new_pixel_spacing[1]
-    ysamples = sampling_rate[1] * min_pixel_spacing[0] / new_pixel_spacing[0]
-    x = np.linspace(extents[0], extents[2], int(xsamples), dtype=np.float)
-    y = np.linspace(extents[1], extents[3], int(ysamples), dtype=np.float)
-    return x, y
+
+    # Existing number of cols, rows
+    num_cols = abs(index_extents[0] - index_extents[2])
+    num_rows = abs(index_extents[1] - index_extents[3])
+
+    col_samples = num_cols * min_pixel_spacing[1] / new_pixel_spacing[1]
+    row_samples = num_rows * min_pixel_spacing[0] / new_pixel_spacing[0]
+
+    col_lut = np.linspace(
+        extents[0], extents[2], int(col_samples), dtype=np.float
+    )
+    row_lut = np.linspace(
+        extents[1], extents[3], int(row_samples), dtype=np.float
+    )
+    return col_lut, row_lut
 
 
 def get_interpolated_dose(dose, z, resolution, extents):
@@ -488,6 +531,8 @@ def get_interpolated_dose(dose, z, resolution, extents):
     """
     # Return the dose bounded by extents if interpolation is not required
     d = dose.GetDoseGrid(z)
+    if not d.size:
+        return d  # cannot take 2d index below if empty
     extent_dose = d[extents[1]:extents[3],
                     extents[0]:extents[2]] if len(extents) else d
     if not resolution:


Test patch:
diff --git a/tests/test_dicomparser.py b/tests/test_dicomparser.py
index 03b8cfe..0733813 100644
--- a/tests/test_dicomparser.py
+++ b/tests/test_dicomparser.py
@@ -143,8 +143,8 @@ class TestImage(unittest.TestCase):
     def test_patient_to_pixel_lut(self):
         """Test if the image transformation matrix (LUT) can be generated."""
         lutvalue = 273.925909
-        self.assertAlmostEqual(
-            self.dp.GetPatientToPixelLUT()[0][-1], lutvalue)
+        doselut = self.dp.GetPatientToPixelLUT()
+        self.assertAlmostEqual(doselut[0][-1], lutvalue)
 
 
 class TestRTStructureSet(unittest.TestCase):
@@ -341,7 +341,8 @@ class TestRTDose(unittest.TestCase):
             'dosegridscaling': 1.4e-05,
             'dosemax': 1048626.0,
             'lut': 253.8458085,
-            'fraction': ''
+            'fraction': '',
+            'x_lut_index': 0,
         }
         dosedata = self.dp.GetDoseData()
         # Pop the LUT numpy array
diff --git a/tests/test_dvhcalc.py b/tests/test_dvhcalc.py
index 4dc6885..383f0d9 100644
--- a/tests/test_dvhcalc.py
+++ b/tests/test_dvhcalc.py
@@ -11,6 +11,7 @@ import os
 from dicompylercore import dicomparser, dvhcalc
 from dicompylercore.config import skimage_available
 from dicompylercore.dvh import DVH
+from dicompylercore.dvhcalc import get_dvh
 try:
     from pydicom.dataset import Dataset
     from pydicom.sequence import Sequence
@@ -18,6 +19,9 @@ except ImportError:
     from dicom.dataset import Dataset
     from dicom.sequence import Sequence
 from numpy import arange
+from numpy.testing import assert_allclose
+from .util import fake_rtdose, fake_ss
+
 
 basedata_dir = "tests/testdata"
 example_data = os.path.join(basedata_dir, "example_data")
@@ -27,7 +31,7 @@ class TestDVHCalc(unittest.TestCase):
     """Unit tests for DVH calculation."""
 
     def setUp(self):
-        """Setup files for common case testing."""
+        """Set up files for common case testing."""
         rtss_dcm = os.path.join(example_data, "rtss.dcm")
         rtdose_dcm = os.path.join(example_data, "rtdose.dcm")
         self.rtss = dicomparser.DicomParser(rtss_dcm)
@@ -53,7 +57,6 @@ class TestDVHCalc(unittest.TestCase):
 
     def create_new_contour(self, roi_id, extents, z):
         """Create a new contour sequence for the given ROI id."""
-
         roic = self.rtss.ds.ROIContourSequence[roi_id - 1]
         new_contour = Dataset()
         # Create a ContourImageSequence for the referenced Image
@@ -259,6 +262,316 @@ class TestDVHCalc(unittest.TestCase):
         self.assertAlmostEqual(dvh.mean, 6.4767105)
 
 
+class TestDVHCalcDecubitus(unittest.TestCase):
+    """Unit tests for DVH calculation in decubitus orientations."""
+
+    def setUp(self):
+        """Set up fake DICOM datasets used in various tests."""
+        self.ss = fake_ss()
+        self.dose = fake_rtdose()
+
+    def test_nondecub(self):
+        """Test that DVH is calculated correctly for standard orientation."""
+        self.dose.ImageOrientationPatient = [1, 0, 0, 0, 1, 0]
+        dvh = get_dvh(self.ss, self.dose, 1)
+        diffl = dvh.differential
+        # Counts are normalized to total, and to volume,
+        # So undo that here for test dose grid.
+        # 18=num dose voxels inside struct; 0.36=volume
+        got_counts = diffl.counts * 18 / 0.36
+        expected_counts = [0]*13 + [2, 2, 2, 0, 0, 0, 0, 0, 0, 0,
+                                    2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2]
+        assert_allclose(got_counts, expected_counts)
+
+    def test_HF_decubitus_left(self):
+        """Test DVH for head-first decubitus left orientation."""
+        # Keep same dose grid as std orientation but pixel-spacing in X, Y same
+        # For this case, use iop=[0, -1, 0, 1, 0, 0] Head first decubitus left
+        # Then X = r * dr + ipp[0]
+        #  and Y = -c * dc + ipp[1]
+        # (https://nipy.org/nibabel/dicom/dicom_orientation.html
+        # #dicom-affine-formula)
+        # Change ipp y of y to new max of 19 for similar y range
+        # Below show contours box of (3, 14.5) - (7, 17.5) on dose grid
+        #       Y=19 18  17                  12
+        # X=2   [10, 10, 10, 13, 14, 15, 16, 17],
+        #               |-----------|
+        #   4   [10, 10, 10, 13, 14, 15, 16, 17]
+        #   6   [10, 10, 10, 13, 14, 15, 16, 17]
+        #               |-----------|
+        #   8   [13, 13, 13, 16, 17, 18, 19, 20]
+        #  10   [14, 14, 14, 17, 18, 19, 20, 21]
+        #  12   [15, 15, 15, 18, 19, 20, 21, 22]
+        #  14   [16, 16, 16, 19, 20, 21, 22, 23]]
+
+        #       Y=19 18  17                  12
+        # X=2   [20, 20, 20, 23, 24, 25, 26, 27]
+        #               |-----------|
+        #   4   [20, 20, 20, 23, 24, 25, 26, 27]
+        #   6   [20, 20, 20, 23, 24, 25, 26, 27]
+        #               |-----------|
+        #   8   [23, 23, 23, 26, 27, 28, 29, 30]
+        #  10   [24, 24, 24, 27, 28, 29, 30, 31]
+        #  12   [25, 25, 25, 28, 29, 30, 31, 32]
+        #  14   [...]
+
+        #       Y=19 18  17                  12
+        # X=2   [30, 30, 30, 33, 34, 35, 36, 37]
+        #               |-----------|
+        #   4   [30, 30, 30, 33, 34, 35, 36, 37]
+        #   6   [30, 30, 30, 33, 34, 35, 36, 37]
+        #               |-----------|
+        #   8   [33, 33, 33, 36, 37, 38, 39, 40]
+        #   10  [34, 34, 34, 37, 38, 39, 40, 41]
+        #   12  [35, 35, 35, 38, 39, 40, 41, 42]
+        # X=14  [36, 36, 36, 39, 40, 41, 42, 43]
+
+        #                          10       13 14                20
+        expected_counts = [0]*10 + [2, 0, 0, 2, 2, 0, 0, 0, 0, 0, 2, 0, 0,
+                                    2, 2, 0, 0, 0, 0, 0, 2, 0, 0, 2, 2]
+        #                          23 24                30       33 34
+        self.dose.ImagePositionPatient = [2, 19, -20]  # X Y Z top left
+        self.dose.PixelSpacing = [2.0, 1.0]  # between Rows, Columns
+        dvh = get_dvh(self.ss, self.dose, 1)
+        diffl = dvh.differential
+        # Counts are normalized to total, and to volume,
+        # So undo that here for test dose grid.
+        # 18=num dose voxels inside struct; 0.36=volume
+        got_counts = diffl.counts * 18 / 0.36
+        assert_allclose(got_counts, expected_counts)
+
+    def test_HF_decubitus_left_structure_extents(self):
+        """Test DVH for HF decubitus Lt orientation structure_extents used."""
+        # Repeat test_HF_decubitus_left but with use_structure_extents
+        #                          10       13 14                20
+        expected_counts = [0]*10 + [2, 0, 0, 2, 2, 0, 0, 0, 0, 0, 2, 0, 0,
+                                    2, 2, 0, 0, 0, 0, 0, 2, 0, 0, 2, 2]
+        #                          23 24                30       33 34
+        self.dose.ImagePositionPatient = [2, 19, -20]  # X Y Z top left
+        self.dose.PixelSpacing = [2.0, 1.0]  # between Rows, Columns
+        dvh = get_dvh(self.ss, self.dose, 1, use_structure_extents=True)
+        diffl = dvh.differential
+        # Counts are normalized to total, and to volume,
+        # So undo that here for test dose grid.
+        # 18=num dose voxels inside struct; 0.36=volume
+        got_counts = diffl.counts * 18 / 0.36
+        assert_allclose(got_counts, expected_counts)
+
+    def test_HF_decubitus_right(self):
+        """Test DVH for head-first decubitus right orientation."""
+        # Keep same dose grid as std orientation
+
+        self.dose.ImageOrientationPatient = [0, 1, 0, -1, 0, 0]
+        self.dose.PixelSpacing = [2.0, 1.0]  # between Rows, Columns
+        # original ipp = [2, 12, -20]
+        # Then X = -r * dr + ipp[0], X decreases down the rows
+        #  and Y = c * dc + ipp[1], Y increases across cols
+        # (https://nipy.org/nibabel/dicom/dicom_orientation.html
+        # #dicom-affine-formula)
+        # Change ipp y of X to new max of 14 for similar y range
+        self.dose.ImagePositionPatient = [14, 12, -20]  # X Y Z top left
+        # Below show contours box of (3, 14.5) - (7, 17.5) on dose grid
+        #       Y=12 13  14  15  16  17  18  19
+        # X=14  [10, 10, 10, 13, 14, 15, 16, 17],
+        #   12  [10, 10, 10, 13, 14, 15, 16, 17]
+        #   10  [10, 10, 10, 13, 14, 15, 16, 17]
+        #    8  [13, 13, 13, 16, 17, 18, 19, 20]
+        #                   | ----------|
+        #    6  [14, 14, 14, 17, 18, 19, 20, 21]
+        #    4  [15, 15, 15, 18, 19, 20, 21, 22]
+        #                   | ----------|
+        #    2  [16, 16, 16, 19, 20, 21, 22, 23]]
+
+        #       Y=12 13  14                  19
+        # X=14  [20, 20, 20, 23, 24, 25, 26, 27]
+        #   12  [20, 20, 20, 23, 24, 25, 26, 27]
+        #   10  [20, 20, 20, 23, 24, 25, 26, 27]
+        #    8  [23, 23, 23, 26, 27, 28, 29, 30]
+        #                   | ----------|
+        #    6  [24, 24, 24, 27, 28, 29, 30, 31]
+        #    4  [25, 25, 25, 28, 29, 30, 31, 32]
+        #                   | ----------|
+        #    2  [...]
+
+        #       Y=12 13  14                  19
+        # X=14  [30, 30, 30, 33, 34, 35, 36, 37]
+        #   12  [30, 30, 30, 33, 34, 35, 36, 37]
+        #   10  [30, 30, 30, 33, 34, 35, 36, 37]
+        #    8  [33, 33, 33, 36, 37, 38, 39, 40]
+        #                   | ----------|
+        #    6  [34, 34, 34, 37, 38, 39, 40, 41]
+        #    4  [35, 35, 35, 38, 39, 40, 41, 42]
+        #                   | ----------|
+        # X= 2  [36, 36, 36, 39, 40, 41, 42, 43]
+
+        #                           17       20
+        expected_counts = [0]*17 + [1, 2, 2, 1, 0, 0, 0, 0, 0, 0,
+                                    1, 2, 2, 1, 0, 0, 0, 0, 0, 0, 1, 2, 2, 1]
+        #                          27 28 29 30                   37
+        dvh = get_dvh(self.ss, self.dose, 1)
+        diffl = dvh.differential
+        # Counts are normalized to total, and to volume,
+        # So undo that here for test dose grid.
+        # 18=num dose voxels inside struct; 0.36=volume
+        got_counts = diffl.counts * 18 / 0.36
+        assert_allclose(got_counts, expected_counts)
+
+    def test_FF_decubitus_right(self):
+        """Test DVH for feet-first decubitus right orientation."""
+        self.dose.ImageOrientationPatient = [0, -1, 0, -1, 0, 0]
+        self.dose.PixelSpacing = [2.0, 1.0]  # between Rows, Columns
+        # original ipp = [2, 12, -20]
+        # Then X = -r * dr + ipp[0], X decreases down the rows
+        #  and Y = -c * dc + ipp[1], Y decreases across cols
+        # (https://nipy.org/nibabel/dicom/dicom_orientation.html
+        # #dicom-affine-formula)
+        self.dose.ImagePositionPatient = [14, 19, 20]  # X Y Z top left
+        # Below show contours box of (3, 14.5) - (7, 17.5) on dose grid
+        #       Y=19 18  17  16  15  14  13  12
+        # X=14  [10, 10, 10, 13, 14, 15, 16, 17],
+        #   12  [10, 10, 10, 13, 14, 15, 16, 17]
+        #   10  [10, 10, 10, 13, 14, 15, 16, 17]
+        #    8  [13, 13, 13, 16, 17, 18, 19, 20]
+        #               | ----------|
+        #    6  [14, 14, 14, 17, 18, 19, 20, 21]
+        #    4  [15, 15, 15, 18, 19, 20, 21, 22]
+        #               | ----------|
+        #    2  [16, 16, 16, 19, 20, 21, 22, 23]]
+
+        #       Y=19 18  17  16  15  14  13  12
+        # X=14  [20, 20, 20, 23, 24, 25, 26, 27]
+        #   12  [20, 20, 20, 23, 24, 25, 26, 27]
+        #   10  [20, 20, 20, 23, 24, 25, 26, 27]
+        #    8  [23, 23, 23, 26, 27, 28, 29, 30]
+        #               | ----------|
+        #    6  [24, 24, 24, 27, 28, 29, 30, 31]
+        #    4  [25, 25, 25, 28, 29, 30, 31, 32]
+        #               | ----------|
+        #    2  [...]
+
+        #       Y=19 18  17  16  15  14  13  12
+        # X=14  [30, 30, 30, 33, 34, 35, 36, 37]
+        #   12  [30, 30, 30, 33, 34, 35, 36, 37]
+        #   10  [30, 30, 30, 33, 34, 35, 36, 37]
+        #    8  [33, 33, 33, 36, 37, 38, 39, 40]
+        #               | ----------|
+        #    6  [34, 34, 34, 37, 38, 39, 40, 41]
+        #    4  [35, 35, 35, 38, 39, 40, 41, 42]
+        #               | ----------|
+        # X= 2  [36, 36, 36, 39, 40, 41, 42, 43]
+
+        #                          14 15 16       19             24
+        expected_counts = [0]*14 + [1, 1, 0, 1, 2, 1, 0, 0, 0, 0, 1, 1, 0, 1,
+                                    2, 1, 0, 0, 0, 0, 1, 1, 0, 1, 2, 1]
+        #                                            34
+        dvh = get_dvh(self.ss, self.dose, 1)
+        diffl = dvh.differential
+        # Counts are normalized to total, and to volume,
+        # So undo that here for test dose grid.
+        # 18=num dose voxels inside struct; 0.36=volume
+        got_counts = diffl.counts * 18 / 0.36
+        assert_allclose(got_counts, expected_counts)
+
+    def test_FF_decubitus_right_structure_extents(self):
+        """Test DVH for FF decubitus Rt orientation using structure extents."""
+        self.dose.ImageOrientationPatient = [0, -1, 0, -1, 0, 0]
+        self.dose.PixelSpacing = [2.0, 1.0]  # between Rows, Columns
+        self.dose.ImagePositionPatient = [14, 19, 20]  # X Y Z top left
+        # see grid from test_FF_decubitus_right
+        #                          14 15 16       19             24
+        expected_counts = [0]*14 + [1, 1, 0, 1, 2, 1, 0, 0, 0, 0, 1, 1, 0,
+                                    1, 2, 1, 0, 0, 0, 0, 1, 1, 0, 1, 2, 1]
+        #                                               34
+        dvh = get_dvh(self.ss, self.dose, 1, use_structure_extents=True)
+        diffl = dvh.differential
+        # Counts are normalized to total, and to volume,
+        # So undo that here for test dose grid.
+        # 18=num dose voxels inside struct; 0.36=volume
+        got_counts = diffl.counts * 18 / 0.36
+        assert_allclose(got_counts, expected_counts)
+
+    def test_FF_decubitus_left(self):
+        """Test DVH for feet-first decubitus left orientation."""
+        self.dose.ImageOrientationPatient = [0, 1, 0, 1, 0, 0]
+        self.dose.PixelSpacing = [2.0, 1.0]  # between Rows, Columns
+        # original ipp = [2, 12, -20]
+        # Then X = r * dr + ipp[0], X increases down the rows
+        #  and Y = c * dc + ipp[1], Y increases across cols
+        # (https://nipy.org/nibabel/dicom/dicom_orientation.html
+        # #dicom-affine-formula)
+
+        # In this test, we also shift Z so three structure planes use the
+        #    first three dose planes rather than the middle three,
+        #    just to ensure asymmetry in z direction is checked.
+        #    Note, planes should really be reversed in pixel array, but doesn't
+        #    matter since contour is identical on each slice.
+        self.dose.ImagePositionPatient = [2, 12, 10]  # X Y Z top left
+        # Below show contours box of (3, 14.5) - (7, 17.5) on dose grid
+        #      Y=12  13  14  15  16  17      19
+        # X=2   [ 0,  0,  0,  3,  4,  5,  6,  7],
+        #                   |-----------|
+        #   4   [ 0,  0,  0,  3,  4,  5,  6,  7]
+        #   6   [ 0,  0,  0,  3,  4,  5,  6,  7]
+        #                   |-----------|
+        #   8   [ 3,  3,  3,  6,  7,  8,  9, 10]
+        #  10   [ 4,  4,  4,  7,  8,  9, 10, 11]
+        #  12   [ 5,  5,  5,  8,  9, 10, 11, 12]
+        #  14   [ 6,  6,  6,  9, 10, 11, 12, 13]]
+
+        #      Y=12  13  14                  19
+        # X=2   [10, 10, 10, 13, 14, 15, 16, 17],
+        #                   |-----------|
+        #   4   [10, 10, 10, 13, 14, 15, 16, 17]
+        #   6   [10, 10, 10, 13, 14, 15, 16, 17]
+        #                   |-----------|
+        #   8   [13, 13, 13, 16, 17, 18, 19, 20]
+        #  10   [14, 14, 14, 17, 18, 19, 20, 21]
+        #  12   [15, 15, 15, 18, 19, 20, 21, 22]
+        #  14   [16, 16, 16, 19, 20, 21, 22, 23]]
+
+        #      Y=12  13  14                  19
+        # X=2   [20, 20, 20, 23, 24, 25, 26, 27]
+        #                   |-----------|
+        #   4   [20, 20, 20, 23, 24, 25, 26, 27]
+        #   6   [20, 20, 20, 23, 24, 25, 26, 27]
+        #                   |-----------|
+        #   8   [23, 23, 23, 26, 27, 28, 29, 30]
+        #  10   [24, 24, 24, 27, 28, 29, 30, 31]
+        #  12   [25, 25, 25, 28, 29, 30, 31, 32]
+        #  14   [...]
+
+        #                          3
+        expected_counts = [0]*3 + [2, 2, 2, 0, 0, 0, 0, 0, 0, 0,
+                                   2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2]
+        #                         13                            23
+        dvh = get_dvh(self.ss, self.dose, 1)
+        diffl = dvh.differential
+        # Counts are normalized to total, and to volume,
+        # So undo that here for test dose grid.
+        # 18=num dose voxels inside struct; 0.36=volume
+        got_counts = diffl.counts * 18 / 0.36
+        assert_allclose(got_counts, expected_counts)
+
+    def test_empty_dose_grid(self):
+        """Test empty dose grid handled correctly."""
+        # See #274, prior to fixes this raised IndexError from
+        #  get_interpolated_dose() getting empty array from GetDoseGrid()
+        # Use z value to force no dose grid at that value
+        #  Otherwise make like decub example
+        self.dose.ImagePositionPatient = [2, 19, -1020]  # X Y Z top left
+        self.dose.PixelSpacing = [2.0, 1.0]  # between Rows, Columns
+
+        # 1 = roi number
+        dvh = get_dvh(self.ss, self.dose, 1, use_structure_extents=True)
+        self.assertTrue('Empty DVH' in dvh.notes)
+
+    def test_not_implemented_orientations(self):
+        """Test unhandled orientations raise NotImplementedError."""
+        self.dose.ImageOrientationPatient = [0.7071, 0.7071, 0, 1, 0, 0]
+        with self.assertRaises(NotImplementedError):
+            _ = get_dvh(self.ss, self.dose, 1)
+
+
 if __name__ == '__main__':
     import sys
     sys.exit(unittest.main())
diff --git a/tests/util.py b/tests/util.py
new file mode 100644
index 0000000..b252ec1
--- /dev/null
+++ b/tests/util.py
@@ -0,0 +1,265 @@
+"""Helper functions for test modules."""
+from pydicom.dataset import Dataset, FileMetaDataset
+from pydicom.sequence import Sequence
+from pydicom.uid import generate_uid
+import numpy
+
+
+# Constants for helper methods below
+SLICE_Z = [-20, -10, 0, 10, 20]
+NUM_SLICES = len(SLICE_Z)
+STUDY_iUID = generate_uid()
+CT_SERIES_iUID = generate_uid()
+CT_iUID = [generate_uid() for i in range(NUM_SLICES)]
+DOSE_SERIES_iUID = generate_uid()
+PLAN_SERIES_iUID = generate_uid()
+PLAN_iUID = generate_uid()
+SS_SERIES_iUID = generate_uid()
+SS_iUID = generate_uid()
+FoR_UID = generate_uid()
+DOSE_iUID = generate_uid()
+
+# Position, orientation
+CT_ipp = [2, 12, -20]
+CT_iop = [0, -1, 0, 1, 0, 0]  # head-first decubitus left
+DOSE_ipp = CT_ipp
+DOSE_iop = CT_iop
+
+RTStructureSetStorage = "1.2.840.10008.5.1.4.1.1.481.3"
+CTImageStorage = "1.2.840.10008.5.1.4.1.1.2"
+RTDoseStorage = "1.2.840.10008.5.1.4.1.1.481.2"
+
+
+# Helper methods for testing DVH calc, particularly decubitus
+def _basic_file_meta(class_UID):
+    """Create basic file meta for DICOM dataset."""
+    # File meta info data elements
+    file_meta = FileMetaDataset()
+    file_meta.FileMetaInformationVersion = b'\x00\x01'
+    file_meta.MediaStorageSOPClassUID = class_UID
+    file_meta.MediaStorageSOPInstanceUID = generate_uid()
+    file_meta.TransferSyntaxUID = '1.2.840.10008.1.2'
+    return file_meta
+
+
+def fake_rtdose():
+    """Create a fake RT Dose DICOM dataset."""
+    # Main data elements
+    file_meta = _basic_file_meta(RTDoseStorage)
+    ds = Dataset()
+    ds.SOPClassUID = RTDoseStorage
+    ds.SOPInstanceUID = DOSE_iUID
+    ds.StudyDate = '20220101'
+    ds.Modality = 'RTDOSE'
+    ds.PatientName = 'Decubitus'
+    ds.PatientID = 'Decubitus'
+    ds.PatientBirthDate = '18000101'
+    ds.PatientSex = 'O'
+
+    ds.SliceThickness = None
+    ds.StudyInstanceUID = STUDY_iUID
+    ds.SeriesInstanceUID = DOSE_SERIES_iUID
+    ds.StudyID = 'none'
+    ds.SeriesNumber = 1
+    ds.ImagePositionPatient = DOSE_ipp
+    ds.ImageOrientationPatient = DOSE_iop
+    ds.FrameOfReferenceUID = FoR_UID
+    ds.PositionReferenceIndicator = ''
+
+    ds.SamplesPerPixel = 1
+    ds.PhotometricInterpretation = 'MONOCHROME2'
+    ds.NumberOfFrames = 5
+    ds.FrameIncrementPointer = (0x3004, 0x000c)  # GridFrameOffsetVector
+    ds.Rows = 7
+    ds.Columns = 8
+    ds.PixelSpacing = [1.0, 2.0]  # (between rows, between cols)
+    ds.BitsAllocated = 32
+    ds.BitsStored = 32
+    ds.HighBit = 31
+    ds.PixelRepresentation = 0
+    ds.DoseUnits = 'GY'
+    ds.DoseType = 'PHYSICAL'
+    ds.DoseSummationType = 'PLAN'
+    ds.GridFrameOffsetVector = [0.0, 10.0, 20.0, 30.0, 40.0]
+
+    # Referenced RT Plan Sequence
+    refd_rt_plan_sequence = Sequence()
+    ds.ReferencedRTPlanSequence = refd_rt_plan_sequence
+
+    # Referenced RT Plan Sequence: Referenced RT Plan 1
+    refd_rt_plan1 = Dataset()
+    refd_rt_plan1.ReferencedSOPClassUID = '1.2.840.10008.5.1.4.1.1.481.5'
+    refd_rt_plan1.ReferencedSOPInstanceUID = PLAN_iUID
+    refd_rt_plan_sequence.append(refd_rt_plan1)
+
+    arr = numpy.zeros(
+        (ds.NumberOfFrames, ds.Rows, ds.Columns), dtype=numpy.uint32
+    )
+
+    ds.DoseGridScaling = 0.00001  # take to near integer cGy in range of < 100
+    # Add small amount so values stay above their int value when dose scaled
+    rounding_guard = 1
+    for zindex in range(len(arr)):
+        for row in range(len(arr[0])):
+            for col in range(len(arr[0][0])):
+                arr[zindex, row, col] = 1000 * (
+                    zindex*10 + (row > 2)*row + (col > 2) * (col)
+                ) + rounding_guard
+
+    # Three middle slices after above math:
+    # # Shown with location of contours:
+    # (3, 14.5), (7, 14.5), (7, 17.5), (3, 17.5)  (2 voxels across x, 3 down y)
+    #     2 mm x 1 mm pixel, across (cols, X) center of pixels from 2 to 16
+    #                         down (rows, Y) centers from 12 to 18
+    # volume = (2*2) x (3*1) x (3*10) = 360 mm^3 = 0.36 cm^3
+    #       X=2   4   6   8  10  12  14  16
+    # Y=12  [10, 10, 10, 13, 14, 15, 16, 17]
+    #       [10, 10, 10, 13, 14, 15, 16, 17]
+    #       [10, 10, 10, 13, 14, 15, 16, 17]
+    #           |-------|
+    #   15  [13, 13, 13, 16, 17, 18, 19, 20]
+    #       [14, 14, 14, 17, 18, 19, 20, 21]
+    #       [15, 15, 15, 18, 19, 20, 21, 22]
+    #           |-------|
+    # Y=18  [16, 16, 16, 19, 20, 21, 22, 23]
+
+    #       X=2   4   6   8  10  12  14  16
+    # Y=12  [20, 20, 20, 23, 24, 25, 26, 27]
+    #       [20, 20, 20, 23, 24, 25, 26, 27]
+    #       [20, 20, 20, 23, 24, 25, 26, 27]
+    #           |-------|
+    #       [23, 23, 23, 26, 27, 28, 29, 30]
+    #       [24, 24, 24, 27, 28, 29, 30, 31]
+    #       [25, 25, 25, 28, 29, 30, 31, 32]
+    #           |-------|
+    # Y=18  [26, 26, 26, 29, 30, 31, 32, 33]
+
+    #       X=2   4   6   8  10  12  14  16
+    # Y=12  [30, 30, 30, 33, 34, 35, 36, 37]
+    #       [30, 30, 30, 33, 34, 35, 36, 37]
+    #       [30, 30, 30, 33, 34, 35, 36, 37]
+    #           |-------|
+    #       [33, 33, 33, 36, 37, 38, 39, 40]
+    #       [34, 34, 34, 37, 38, 39, 40, 41]
+    #       [35, 35, 35, 38, 39, 40, 41, 42]
+    #           |-------|
+    # Y=18  [36, 36, 36, 39, 40, 41, 42, 43]
+
+    ds.PixelData = arr.tobytes()
+    ds.file_meta = file_meta
+    ds.is_implicit_VR = True
+    ds.is_little_endian = True
+
+    return ds
+
+
+def fake_ss():
+    """Create a fake RT Structure Set DICOM dataset."""
+    file_meta = _basic_file_meta(RTStructureSetStorage)
+    # Main data elements
+    ds = Dataset()
+    ds.SpecificCharacterSet = 'ISO_IR 192'
+    ds.InstanceCreationDate = '20220101'
+    ds.SOPClassUID = RTStructureSetStorage
+    ds.SOPInstanceUID = SS_iUID
+    ds.StudyDate = '20220101'
+    ds.Modality = 'RTSTRUCT'
+
+    ds.PatientName = 'Decubitus'
+    ds.PatientID = 'Decubitus'
+    ds.PatientBirthDate = '18000101'
+    ds.PatientSex = 'O'
+    ds.StudyInstanceUID = STUDY_iUID
+    ds.SeriesInstanceUID = SS_SERIES_iUID
+    ds.StudyID = '1'
+    ds.SeriesNumber = '10'
+    ds.StructureSetLabel = 'CT_1'
+    ds.StructureSetDate = '20220101'
+
+    # Referenced Frame of Reference Sequence
+    refd_frame_of_ref_sequence = Sequence()
+    ds.ReferencedFrameOfReferenceSequence = refd_frame_of_ref_sequence
+
+    # Referenced Frame of Reference Sequence: Referenced Frame of Reference 1
+    refd_frame_of_ref1 = Dataset()
+    refd_frame_of_ref1.FrameOfReferenceUID = FoR_UID
+
+    # RT Referenced Study Sequence
+    rt_refd_study_sequence = Sequence()
+    refd_frame_of_ref1.RTReferencedStudySequence = rt_refd_study_sequence
+
+    # RT Referenced Study Sequence: RT Referenced Study 1
+    rt_refd_study1 = Dataset()
+    rt_refd_study1.ReferencedSOPClassUID = '1.2.840.10008.5.1.4.1.1.481.3'
+    rt_refd_study1.ReferencedSOPInstanceUID = STUDY_iUID
+
+    # RT Referenced Series Sequence
+    rt_refd_series_sequence = Sequence()
+    rt_refd_study1.RTReferencedSeriesSequence = rt_refd_series_sequence
+
+    # RT Referenced Series Sequence: RT Referenced Series 1
+    rt_refd_series1 = Dataset()
+    rt_refd_series1.SeriesInstanceUID = CT_SERIES_iUID
+
+    # Contour Image Sequence
+    contour_image_sequence = Sequence()
+    rt_refd_series1.ContourImageSequence = contour_image_sequence
+
+    # Contour Image Sequence: Contour Image 1
+    contour_image1 = Dataset()
+    contour_image1.ReferencedSOPClassUID = CTImageStorage
+    contour_image1.ReferencedSOPInstanceUID = CT_iUID
+    contour_image_sequence.append(contour_image1)
+
+    rt_refd_series_sequence.append(rt_refd_series1)
+    rt_refd_study_sequence.append(rt_refd_study1)
+    refd_frame_of_ref_sequence.append(refd_frame_of_ref1)
+
+    # Structure Set ROI Sequence
+    structure_set_roi_sequence = Sequence()
+    ds.StructureSetROISequence = structure_set_roi_sequence
+
+    # Structure Set ROI Sequence: Structure Set ROI 1
+    structure_set_roi1 = Dataset()
+    structure_set_roi1.ROINumber = 1
+    structure_set_roi1.ReferencedFrameOfReferenceUID = FoR_UID
+    structure_set_roi1.ROIName = 'PTV'
+    structure_set_roi1.ROIGenerationAlgorithm = 'MANUAL'
+    structure_set_roi_sequence.append(structure_set_roi1)
+
+    # ROI Contour Sequence
+    #    Contour Sequence
+    #        Contour Image Sequence (optional, references the CT slices)
+    #        Contour Geometric Type (e.g. CLOSED_PLANAR)
+    #        Number of Contour Points
+    #        Contour Data
+    #    Referenced ROI Number
+
+    roi_contour_sequence = Sequence()
+    ds.ROIContourSequence = roi_contour_sequence
+
+    # ROI Contour Sequence: ROI Contour 1
+    roi_contour1 = Dataset()
+    roi_contour1.ROIDisplayColor = [255, 128, 0]
+
+    # Contour Sequence
+    contour_sequence = Sequence()
+    roi_contour1.ContourSequence = contour_sequence
+
+    for i in range(1, NUM_SLICES-1):  # contours on all but first, last slice
+        contour = Dataset()
+        contour.ContourGeometricType = 'CLOSED_PLANAR'
+
+        z = SLICE_Z[i]
+        contour.ContourData = [3, 14.5, z, 7, 14.5, z, 7, 17.5, z, 3, 17.5, z]
+        contour.NumberOfContourPoints = len(contour.ContourData) / 3
+        contour_sequence.append(contour)
+
+    roi_contour1.ReferencedROINumber = 1
+    roi_contour_sequence.append(roi_contour1)
+
+    ds.file_meta = file_meta
+    ds.is_implicit_VR = True
+    ds.is_little_endian = True
+
+    return ds