dicom_pixeldata/
transform.rs

1//! Private module for pixel sample value transformation functions.
2
3use snafu::Snafu;
4
5use crate::attribute::VoiLut;
6
7/// Description of a modality rescale function,
8/// defined by a _rescale slope_ and _rescale intercept_.
9#[derive(Debug, Copy, Clone, PartialEq)]
10pub struct Rescale {
11    /// the rescale slope
12    pub slope: f64,
13    /// the rescale intercept
14    pub intercept: f64,
15}
16
17impl Rescale {
18    /// Create a new rescale function.
19    #[inline]
20    pub fn new(slope: f64, intercept: f64) -> Self {
21        Rescale { slope, intercept }
22    }
23
24    /// Apply the rescale function to a value.
25    #[inline]
26    pub fn apply(&self, value: f64) -> f64 {
27        self.slope * value + self.intercept
28    }
29}
30
31/// A known DICOM Value of Interest (VOI) LUT function descriptor.
32#[derive(Debug, Default, Copy, Clone, Eq, Hash, PartialEq)]
33pub enum VoiLutFunction {
34    /// LINEAR
35    #[default]
36    Linear,
37    /// LINEAR_EXACT
38    LinearExact,
39    /// SIGMOID
40    Sigmoid,
41}
42
43/// Unrecognized VOI LUT function name
44#[derive(Debug, Copy, Clone, PartialEq, Snafu)]
45pub struct FromVoiLutFunctionError {
46    _private: (),
47}
48
49impl std::convert::TryFrom<&str> for VoiLutFunction {
50    type Error = FromVoiLutFunctionError;
51
52    fn try_from(s: &str) -> Result<Self, Self::Error> {
53        match s {
54            "LINEAR" => Ok(Self::Linear),
55            "LINEAR_EXACT" => Ok(Self::LinearExact),
56            "SIGMOID" => Ok(Self::Sigmoid),
57            _ => Err(FromVoiLutFunctionError { _private: () }),
58        }
59    }
60}
61
62/// The parameters of a single window level
63/// for a VOI LUT transformation,
64/// comprising the window center and the window width.
65#[derive(Debug, Copy, Clone, PartialEq)]
66pub struct WindowLevel {
67    /// The _Window Width_.
68    ///
69    /// Should be greater than 0
70    pub width: f64,
71    /// The _Window Center_.
72    pub center: f64,
73}
74
75/// A full description of a VOI LUT function transformation
76/// based on a window level.
77#[derive(Debug, PartialEq)]
78pub struct WindowLevelTransform {
79    voi_lut_function: VoiLutFunction,
80    window_level: WindowLevel,
81}
82
83impl WindowLevelTransform {
84    /// Create a new window level transformation.
85    ///
86    /// The width of the given `window_level` is automatically clamped
87    /// if it is incompatible with the given LUT function:
88    /// it muse be `>= 0` if the function is [`LinearExact`](VoiLutFunction::LinearExact),
89    /// and `>= 1` in other functions.
90    #[inline]
91    pub fn new(voi_lut_function: VoiLutFunction, window_level: WindowLevel) -> Self {
92        WindowLevelTransform {
93            voi_lut_function,
94            window_level: WindowLevel {
95                center: window_level.center,
96                width: match voi_lut_function {
97                    VoiLutFunction::LinearExact => window_level.width.max(0.),
98                    VoiLutFunction::Linear | VoiLutFunction::Sigmoid => window_level.width.max(1.),
99                },
100            },
101        }
102    }
103
104    /// Create a new window level transformation
105    /// with the `LINEAR` function.
106    ///
107    /// The width of the given `window_level` is automatically clamped
108    /// to 1 if it is lower than 1.
109    #[inline]
110    pub fn linear(window_level: WindowLevel) -> Self {
111        Self::new(VoiLutFunction::Linear, window_level)
112    }
113
114    /// Apply the window level transformation on a rescaled value,
115    /// into a number between `0` and `y_max`.
116    pub fn apply(&self, value: f64, y_max: f64) -> f64 {
117        match self.voi_lut_function {
118            VoiLutFunction::Linear => window_level_linear(
119                value,
120                self.window_level.width,
121                self.window_level.center,
122                y_max,
123            ),
124            VoiLutFunction::LinearExact => window_level_linear_exact(
125                value,
126                self.window_level.width,
127                self.window_level.center,
128                y_max,
129            ),
130            VoiLutFunction::Sigmoid => window_level_sigmoid(
131                value,
132                self.window_level.width,
133                self.window_level.center,
134                y_max,
135            ),
136        }
137    }
138}
139
140fn window_level_linear(value: f64, window_width: f64, window_center: f64, y_max: f64) -> f64 {
141    let ww = window_width;
142    let wc = window_center;
143    debug_assert!(ww >= 1.);
144
145    // C.11.2.1.2.1
146
147    let min = wc - (ww - 1.) / 2.;
148    let max = wc - 0.5 + (ww - 1.) / 2.;
149
150    if value <= min {
151        // if (x <= c - (w-1) / 2), then y = ymin
152        0.
153    } else if value > max {
154        // else if (x > c - 0.5 + (w-1) /2), then y = ymax
155        y_max
156    } else {
157        // else y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax- ymin) + ymin
158        ((value - (wc - 0.5)) / (ww - 1.) + 0.5) * y_max
159    }
160}
161
162fn window_level_linear_exact(value: f64, window_width: f64, window_center: f64, y_max: f64) -> f64 {
163    let ww = window_width;
164    let wc = window_center;
165    debug_assert!(ww >= 0.);
166
167    // C.11.2.1.3.2
168
169    let min = wc - ww / 2.;
170    let max = wc + ww / 2.;
171
172    if value <= min {
173        // if (x <= c - w/2), then y = ymin
174        0.
175    } else if value > max {
176        // else if (x > c + w/2), then y = ymax
177        y_max
178    } else {
179        // else y = ((x - c) / w + 0.5) * (ymax - ymin) + ymin
180        ((value - wc) / ww + 0.5) * y_max
181    }
182}
183
184fn window_level_sigmoid(value: f64, window_width: f64, window_center: f64, y_max: f64) -> f64 {
185    let ww = window_width;
186    let wc = window_center;
187    assert!(ww >= 1.);
188
189    // C.11.2.1.3.1
190
191    y_max / (1. + f64::exp(-4. * (value - wc) / ww))
192}
193
194/// A full description of a VOI LUT function transformation
195/// based on an explicit LUT.
196pub struct VoiLutTransform<'a> {
197    lut: &'a VoiLut,
198    shift: u32,
199}
200
201impl<'a> VoiLutTransform<'a> {
202    /// Create a new LUT transformation.
203    ///
204    /// # Panics
205    ///
206    /// Panics if `lut.data` is empty, or if `bits_stored` is smaller than
207    /// `lut.bits_stored`.
208    pub fn new(lut: &'a VoiLut, bits_stored: u16) -> Self {
209        if lut.data.is_empty() {
210            // we'll do unchecked accesses to the first and last items in apply()
211            panic!("LUT data is empty");
212        }
213
214        let next_pow_bits = bits_stored.next_power_of_two();
215
216        if next_pow_bits < (lut.bits_stored as u16) {
217            panic!(
218                "LUT with BitsStored {} cannot be used for an image with BitsStored {}",
219                lut.bits_stored, bits_stored
220            );
221        }
222
223        let shift = (next_pow_bits - (lut.bits_stored as u16)) as u32;
224
225        VoiLutTransform { lut, shift }
226    }
227
228    /// Apply the LUT transformation on a rescaled value.
229    pub fn apply(&self, value: f64) -> f64 {
230        let value_rounded = value.round().clamp(i32::MIN as f64, i32::MAX as f64) as i32;
231        let y = (*(if value_rounded <= self.lut.min_pixel_value {
232            self.lut.data.first().unwrap()
233        } else {
234            self.lut
235                .data
236                .get((value_rounded - self.lut.min_pixel_value) as usize)
237                .unwrap_or(self.lut.data.last().unwrap())
238        })) << self.shift;
239        y as f64
240    }
241}
242
243#[cfg(test)]
244mod tests {
245    use crate::Lut;
246
247    use super::*;
248
249    /// Applying a common rescale function to a value
250    /// gives the expected output.
251    #[test]
252    fn modality_lut_baseline() {
253        let rescale = Rescale::new(1., -1024.);
254
255        assert_eq!(rescale.apply(0.), -1024.);
256        assert_eq!(rescale.apply(1.), -1023.);
257        assert_eq!(rescale.apply(2.), -1022.);
258        assert_eq!(rescale.apply(1024.), 0.);
259    }
260
261    /// Applying a linear window level
262    /// as per the example described in the standard
263    /// (C.11.2.1.2.1)
264    /// gives us the expected outcome.
265    #[test]
266    fn window_level_linear_example() {
267        let window_level = WindowLevel {
268            width: 4096.,
269            center: 2048.,
270        };
271        let window_level_transform = WindowLevelTransform::linear(window_level);
272        let y_max = 255.;
273
274        // x <= 0
275        assert_eq!(window_level_transform.apply(-2., y_max), 0.);
276        assert_eq!(window_level_transform.apply(-1., y_max), 0.);
277        assert_eq!(window_level_transform.apply(0., y_max), 0.);
278
279        // x > 4095
280        assert_eq!(window_level_transform.apply(4095.5, y_max), y_max);
281        assert_eq!(window_level_transform.apply(4096., y_max), y_max);
282        assert_eq!(window_level_transform.apply(4097., y_max), y_max);
283
284        // inbetween:  y = ((x - 2047.5) / 4095 + 0.5) * 255
285
286        let x = 1024.;
287        let y = window_level_transform.apply(x, y_max);
288        let expected_y = ((x - 2047.5) / 4095. + 0.5) * 255.;
289
290        assert!((y - expected_y).abs() < 1e-3);
291    }
292
293    /// Applying a linear window level gives us the expected outcome.
294    #[test]
295    fn window_level_linear_1() {
296        let window_level = WindowLevel {
297            width: 300.,
298            center: 50.,
299        };
300        let window_level_transform = WindowLevelTransform::linear(window_level);
301        let y_max = 255.;
302
303        // x <= (50 - 150)
304        let y = window_level_transform.apply(-120., y_max);
305        assert_eq!(y, 0.);
306        let y = window_level_transform.apply(-100.5, y_max);
307        assert_eq!(y, 0.);
308        let y = window_level_transform.apply(-100., y_max);
309        assert_eq!(y, 0.);
310
311        // x >= (50 + 150)
312        let y = window_level_transform.apply(201., y_max);
313        assert_eq!(y, 255.);
314        let y = window_level_transform.apply(200., y_max);
315        assert_eq!(y, 255.);
316
317        // x inbetween
318        let y = window_level_transform.apply(50., y_max);
319        assert!(y > 127. && y < 129.);
320    }
321
322    #[test]
323    fn voi_lut_transform() {
324        let voi_data = (0..256).map(|x| x * 4).collect();
325        let voi = VoiLut {
326            min_pixel_value: 0,
327            bits_stored: 8,
328            data: voi_data,
329            explanation: None,
330        };
331        let rescale = Rescale::new(1., 0.);
332        let lut: Lut<u16> =
333            Lut::new_rescale_and_lut(8, false, rescale, VoiLutTransform::new(&voi, 8)).unwrap();
334
335        // test some values
336
337        assert_eq!(lut.get(0_u16), 0);
338        assert_eq!(lut.get(1_u16), 4);
339        assert_eq!(lut.get(2_u16), 8);
340        assert_eq!(lut.get(255_u16), 1020);
341    }
342}
dicom_pixeldata/transform.rs

dicom_pixeldata/
transform.rs
