|
18 | 18 |
|
19 | 19 | from monai.deploy.utils.importutil import optional_import |
20 | 20 |
|
| 21 | +apply_presentation_lut, _ = optional_import("pydicom.pixels", name="apply_presentation_lut") |
21 | 22 | apply_rescale, _ = optional_import("pydicom.pixels", name="apply_rescale") |
22 | 23 |
|
23 | 24 | from monai.deploy.core import ConditionType, Fragment, Operator, OperatorSpec |
@@ -109,72 +110,82 @@ def generate_voxel_data(self, series): |
109 | 110 | Returns: |
110 | 111 | A 3D numpy tensor representing the volumetric data. |
111 | 112 | """ |
| 113 | + |
| 114 | + def _get_rescaled_pixel_array(sop_instance): |
| 115 | + # Use pydicom utility to apply a modality lookup table or rescale operator to the pixel array. |
| 116 | + # The pydicom Dataset is required which can be obtained from the first slice's native SOP instance. |
| 117 | + # If Modality LUT is present the return array is of np.uint8 or np.uint16, and if Rescale |
| 118 | + # Intercept and Rescale Slope are present, np.float64. |
| 119 | + # If the pixel array is already in the correct type, the return array is the same as the input array. |
| 120 | + |
| 121 | + if not sop_instance: |
| 122 | + return np.array([]) |
| 123 | + |
| 124 | + native_sop = None |
| 125 | + try: |
| 126 | + native_sop = sop_instance.get_native_sop_instance() |
| 127 | + rescaled_pixel_data = apply_rescale(sop_instance.get_pixel_array(), native_sop) |
| 128 | + # In our use cases, pixel data will be interpreted as if MONOCHROME2, hence need to |
| 129 | + # apply the presentation lut. |
| 130 | + rescaled_pixel_data = apply_presentation_lut(rescaled_pixel_data, native_sop) |
| 131 | + except Exception as e: |
| 132 | + logging.error(f"Failed to apply rescale to DICOM volume: {e}") |
| 133 | + raise RuntimeError("Failed to apply rescale to DICOM volume.") from e |
| 134 | + |
| 135 | + # The following tests are expecting the array is already of the Numpy type. |
| 136 | + if rescaled_pixel_data.dtype == np.uint8 or rescaled_pixel_data.dtype == np.uint16: |
| 137 | + logging.info("Rescaled pixel array is alreadydy of type uint8 or uint16.") |
| 138 | + # Check if casting to uint16 and back to float results in the same values. |
| 139 | + elif np.all(rescaled_pixel_data > 0) and np.array_equal( |
| 140 | + rescaled_pixel_data, rescaled_pixel_data.astype(np.uint16) |
| 141 | + ): |
| 142 | + logging.info("Rescaled pixel array can be safely cast to uint16 with equivalence test.") |
| 143 | + rescaled_pixel_data = rescaled_pixel_data.astype(dtype=np.uint16) |
| 144 | + # Check if casting to int16 and back to float results in the same values. |
| 145 | + elif np.array_equal(rescaled_pixel_data, rescaled_pixel_data.astype(np.int16)): |
| 146 | + logging.info("Rescaled pixel array can be safely cast to int16 with equivalence test.") |
| 147 | + rescaled_pixel_data = rescaled_pixel_data.astype(dtype=np.int16) |
| 148 | + # Check casting to float32 with equivalence test |
| 149 | + elif np.array_equal(rescaled_pixel_data, rescaled_pixel_data.astype(np.float32)): |
| 150 | + logging.info("Rescaled pixel array can be cast to float32 with equivalence test.") |
| 151 | + rescaled_pixel_data = rescaled_pixel_data.astype(np.float32) |
| 152 | + else: |
| 153 | + logging.info("Rescaled pixel data remains as of type float64.") |
| 154 | + |
| 155 | + return rescaled_pixel_data |
| 156 | + |
112 | 157 | slices = series.get_sop_instances() |
113 | 158 | # The sop_instance get_pixel_array() returns a 2D NumPy array with index order |
114 | 159 | # of `HW`. The pixel array of all instances will be stacked along the first axis, |
115 | 160 | # so the final 3D NumPy array will have index order of [DHW]. This is consistent |
116 | 161 | # with the NumPy array returned from the ITK GetArrayViewFromImage on the image |
117 | 162 | # loaded from the same DICOM series. |
118 | | - vol_data = np.stack([s.get_pixel_array() for s in slices], axis=0) |
| 163 | + # The below code loads all slice pixel data into a list of NumPy arrays in memory |
| 164 | + # before stacking them into a single 3D volume. This can be inefficient for series |
| 165 | + # with many slices. |
| 166 | + if not slices: |
| 167 | + return np.array([]) |
119 | 168 |
|
120 | | - # Use pydicom utility to apply a modality lookup table or rescale operator to the pixel array. |
121 | | - # The pydicom Dataset is required which can be obtained from the first slice's native SOP instance. |
122 | | - # If Modality LUT is present the return array is of np.uint8 or np.uint16, and if Rescale |
123 | | - # Intercept and Rescale Slope are present, np.float64. |
124 | | - # If the pixel array is already in the correct type, the return array is the same as the input array. |
| 169 | + # Get shape and dtype from the first slice to pre-allocate numpy array. |
125 | 170 | try: |
126 | | - native_sop = slices[0].get_native_sop_instance() |
127 | | - vol_data = apply_rescale(vol_data, native_sop) |
| 171 | + first_slice_pixel_array = _get_rescaled_pixel_array(slices[0]) |
| 172 | + vol_shape = (len(slices),) + first_slice_pixel_array.shape |
| 173 | + dtype = first_slice_pixel_array.dtype |
128 | 174 | except Exception as e: |
129 | | - logging.error(f"Failed to apply rescale to DICOM volume: {e}") |
130 | | - raise RuntimeError("Failed to apply rescale to DICOM volume.") from e |
131 | | - |
132 | | - # For now we support monochrome image only, for which DICOM Photometric Interpretation |
133 | | - # (0028,0004) has defined terms, MONOCHROME1 and MONOCHROME2, with the former being: |
134 | | - # Pixel data represent a single monochrome image plane. The minimum sample value is |
135 | | - # intended to be displayed as white after any VOI gray scale transformations have been |
136 | | - # performed. See PS3.4. This value may be used only when Samples per Pixel (0028,0002) |
137 | | - # has a value of 1. May be used for pixel data in a Native (uncompressed) or Encapsulated |
138 | | - # (compressed) format; see Section 8.2 in PS3.5. |
139 | | - # and for the latter "The minimum sample value is intended to be displayed as black" |
140 | | - # |
141 | | - # In this function, pixel data will be interpreted as if MONOCHROME2, hence inverting |
142 | | - # MONOCHROME1 for the final voxel data. |
143 | | - |
144 | | - photometric_interpretation = ( |
145 | | - slices[0].get_native_sop_instance().get("PhotometricInterpretation", "").strip().upper() |
146 | | - ) |
147 | | - presentation_lut_shape = slices[0].get_native_sop_instance().get("PresentationLUTShape", "").strip().upper() |
148 | | - |
149 | | - if not photometric_interpretation: |
150 | | - logging.warning("Cannot get value of attribute Photometric Interpretation.") |
151 | | - |
152 | | - if photometric_interpretation != "MONOCHROME2": |
153 | | - if photometric_interpretation == "MONOCHROME1" or presentation_lut_shape == "INVERSE": |
154 | | - logging.debug("Applying INVERSE transformation as required for MONOCHROME1 image.") |
155 | | - vol_data = np.amax(vol_data) - vol_data |
156 | | - else: |
157 | | - raise ValueError( |
158 | | - f"Cannot process pixel data with Photometric Interpretation of {photometric_interpretation}." |
159 | | - ) |
160 | | - |
161 | | - # The following tests are expecting the array is already of the Numpy type. |
162 | | - if vol_data.dtype == np.uint8 or vol_data.dtype == np.uint16: |
163 | | - logging.info("Rescaled pixel array is alrady of type uint8 or uint16.") |
164 | | - # Check if casting to uint16 and back to float results in the same values. |
165 | | - elif np.all(vol_data > 0) and np.array_equal(vol_data, vol_data.astype(np.uint16)): |
166 | | - logging.info("Rescaled pixel array can be safely cast to uint16 with equivalence test.") |
167 | | - vol_data = vol_data.astype(dtype=np.uint16) |
168 | | - # Check if casting to int16 and back to float results in the same values. |
169 | | - elif np.array_equal(vol_data, vol_data.astype(np.int16)): |
170 | | - logging.info("Rescaled pixel array can be safely cast to int16 with equivalence test.") |
171 | | - vol_data = vol_data.astype(dtype=np.int16) |
172 | | - # Check casting to float32 with equivalence test |
173 | | - elif np.array_equal(vol_data, vol_data.astype(np.float32)): |
174 | | - logging.info("Rescaled pixel array can be cast to float32 with equivalence test.") |
175 | | - vol_data = vol_data.astype(np.float32) |
176 | | - else: |
177 | | - logging.info("Rescaled pixel data remains as of type float64.") |
| 175 | + logging.error(f"Failed to get pixel array from the first slice: {e}") |
| 176 | + raise |
| 177 | + |
| 178 | + # Pre-allocate the volume data array. |
| 179 | + vol_data = np.empty(vol_shape, dtype=dtype) |
| 180 | + vol_data[0] = first_slice_pixel_array |
| 181 | + |
| 182 | + # Read subsequent slices directly into the pre-allocated array. |
| 183 | + for i, s in enumerate(slices[1:], 1): |
| 184 | + try: |
| 185 | + vol_data[i] = _get_rescaled_pixel_array(s) |
| 186 | + except Exception as e: |
| 187 | + logging.error(f"Failed to get pixel array from slice {i}: {e}") |
| 188 | + raise |
178 | 189 |
|
179 | 190 | return vol_data |
180 | 191 |
|
|
0 commit comments