From afffacbca04ab0e631d375a8615b099948df322d Mon Sep 17 00:00:00 2001 From: Adrien Bustany Date: Thu, 17 Oct 2024 12:30:46 +0200 Subject: [PATCH] [pixeldata] Apply VOI LUT when rendering images Fixes: #232 --- pixeldata/src/attribute.rs | 150 +++++++++++++++++++++++++ pixeldata/src/gdcm.rs | 4 + pixeldata/src/lib.rs | 221 +++++++++++++++++++++++++++---------- pixeldata/src/lut.rs | 63 ++++++++++- pixeldata/src/transform.rs | 49 ++++++++ 5 files changed, 426 insertions(+), 61 deletions(-) diff --git a/pixeldata/src/attribute.rs b/pixeldata/src/attribute.rs index 825215562..f3c4706c3 100644 --- a/pixeldata/src/attribute.rs +++ b/pixeldata/src/attribute.rs @@ -30,6 +30,9 @@ pub enum AttributeName { VoiLutFunction, WindowCenter, WindowWidth, + LutDescriptor, + LutData, + LutExplanation, } impl std::fmt::Display for AttributeName { @@ -610,6 +613,153 @@ pub fn photometric_interpretation( .into()) } +/// A decoded representation of the +/// DICOM _VOI LUT Sequence_ attribute. +/// +/// See [section C.8.11.3.1.5][1] of the standard for more details. +/// +/// [1]: https://dicom.nema.org/dicom/2013/output/chtml/part03/sect_C.8.html#sect_C.8.11.3.1.5 +#[derive(Clone, Debug)] +pub struct VoiLut { + /// Minimum pixel value to be mapped. All values below this should be mapped + /// to the first entry of the table. + /// + /// The value can either be a signed or an unsigned 16-bit integer, an i32 + /// accomodates both. + pub min_pixel_value: i32, + /// Number of bits stored for each LUT entry + pub bits_stored: u8, + /// LUT data. Pixels with value min_pixel_value or below get mapped to the + /// first entry, pixels with value min_pixel_value+1 to the second entry, + /// and so forth. Pixels with value higher or equal to min_pixel_value + + /// data.len() get mapped to the last entry. + pub data: Vec, + /// Free form text explanation of the meaning of the LUT. + pub explanation: Option, +} + +fn parse_voi_lut_entry(entry: &InMemDicomObject) -> Result { + let descriptor_elements: Vec = entry + .element_opt(tags::LUT_DESCRIPTOR) + .context(RetrieveSnafu { + name: AttributeName::LutDescriptor, + })? + .context(MissingRequiredSnafu { + name: AttributeName::LutDescriptor, + })? + .to_multi_int() + .context(ConvertValueSnafu { + name: AttributeName::LutDescriptor, + })?; + ensure!( + descriptor_elements.len() == 3, + InvalidValueSnafu { + name: AttributeName::LutDescriptor, + value: format!("value with multiplicity {}", descriptor_elements.len()), + } + ); + ensure!( + descriptor_elements[0] > 0, + InvalidValueSnafu { + name: AttributeName::LutDescriptor, + value: format!("value with LUT length {}", descriptor_elements[0]), + } + ); + let expected_lut_len: usize = descriptor_elements[0] as usize; + + let min_pixel_value = descriptor_elements[1]; + + ensure!( + descriptor_elements[2] >= 0 && descriptor_elements[2] <= 16, + InvalidValueSnafu { + name: AttributeName::LutDescriptor, + value: format!("value with bits stored {}", descriptor_elements[2]) + } + ); + let bits_stored = descriptor_elements[2] as u8; + + let lut_data = entry + .element_opt(tags::LUT_DATA) + .context(RetrieveSnafu { + name: AttributeName::LutData, + })? + .context(MissingRequiredSnafu { + name: AttributeName::LutData, + })? + .uint16_slice() + .context(CastValueSnafu { + name: AttributeName::LutData, + })?; + ensure!( + expected_lut_len == lut_data.len(), + InvalidValueSnafu { + name: AttributeName::LutData, + value: format!( + "sequence with {} elements (expected {expected_lut_len})", + lut_data.len() + ), + } + ); + + let explanation = if let Some(val) = + entry + .element_opt(tags::LUT_EXPLANATION) + .context(RetrieveSnafu { + name: AttributeName::LutExplanation, + })? { + Some( + val.string() + .context(CastValueSnafu { + name: AttributeName::LutExplanation, + })? + .to_string(), + ) + } else { + None + }; + + Ok(VoiLut { + min_pixel_value, + bits_stored, + data: lut_data.to_vec(), + explanation, + }) +} + +pub fn voi_lut_sequence( + obj: &FileDicomObject>, +) -> Option> { + obj.get(tags::VOILUT_SEQUENCE) + .and_then(|e| { + e.items().and_then(|items| { + items + .iter() + .map(|e| parse_voi_lut_entry(e).ok()) + .collect::>>() + }) + }) + .or_else(|| { + get_from_per_frame(obj, [tags::FRAME_VOILUT_SEQUENCE, tags::VOILUT_SEQUENCE]).and_then( + |v| { + v.into_iter() + .flat_map(|el| el.items()) + .flat_map(|items| items.iter().map(|e| parse_voi_lut_entry(e).ok())) + .collect() + }, + ) + }) + .or_else(|| { + get_from_shared(obj, [tags::FRAME_VOILUT_SEQUENCE, tags::VOILUT_SEQUENCE]).and_then( + |v| { + v.into_iter() + .flat_map(|el| el.items()) + .flat_map(|items| items.iter().map(|e| parse_voi_lut_entry(e).ok())) + .collect() + }, + ) + }) +} + #[cfg(test)] mod tests { use super::rescale_intercept; diff --git a/pixeldata/src/gdcm.rs b/pixeldata/src/gdcm.rs index e24fb7377..917827114 100644 --- a/pixeldata/src/gdcm.rs +++ b/pixeldata/src/gdcm.rs @@ -71,6 +71,7 @@ where .map(|v| VoiLutFunction::try_from((*v).as_str()).ok()) .collect() }); + let voi_lut_sequence = voi_lut_sequence(self); ensure!( rescale_intercept.len() == rescale_slope.len(), @@ -178,6 +179,7 @@ where rescale, voi_lut_function, window, + voi_lut_sequence, enforce_frame_fg_vm_match: false, }) } @@ -236,6 +238,7 @@ where .map(|v| VoiLutFunction::try_from((*v).as_str()).ok()) .collect() }); + let voi_lut_sequence = voi_lut_sequence(self); let decoded_pixel_data = match pixel_data.value() { DicomValue::PixelSequence(v) => { @@ -356,6 +359,7 @@ where rescale: rescale, voi_lut_function, window, + voi_lut_sequence, enforce_frame_fg_vm_match: false, }) } diff --git a/pixeldata/src/lib.rs b/pixeldata/src/lib.rs index f9a8616c4..71b3e798a 100644 --- a/pixeldata/src/lib.rs +++ b/pixeldata/src/lib.rs @@ -102,6 +102,7 @@ //! including the default behavior for each method. //! +use attribute::VoiLut; use byteorder::{ByteOrder, NativeEndian}; #[cfg(not(feature = "gdcm"))] use dicom_core::{DataDictionary, DicomValue}; @@ -245,6 +246,16 @@ pub enum InnerError { ww_vm: u32, backtrace: Backtrace, }, + #[snafu(display( + "Value multiplicity of VOI LUT must match the number of frames. Expected `{:?}`, found `{:?}`", + nr_frames, + vm + ))] + LengthMismatchVoiLut { + vm: u32, + nr_frames: u32, + backtrace: Backtrace, + }, } pub type Result = std::result::Result; @@ -355,6 +366,9 @@ pub enum VoiLutOption { /// only when converting to an image; /// no VOI LUT function is performed /// when converting to an ndarray or to bare pixel values. + /// + /// If both a VOI LUT table and window values are present in the pixel data, + /// use the VOI LUT table. #[default] Default, /// Apply the first VOI LUT function transformation @@ -438,6 +452,8 @@ pub struct DecodedPixelData<'a> { voi_lut_function: Option>, /// the window level specified via width and center window: Option>, + /// the explicit VOI LUTs + voi_lut_sequence: Option>, /// Enforce frame functional groups VMs match `number_of_frames` enforce_frame_fg_vm_match: bool, @@ -642,6 +658,36 @@ impl DecodedPixelData<'_> { } } + pub fn voi_lut_sequence(&self) -> Result> { + if let Some(inner) = &self.voi_lut_sequence { + match &inner.len() { + 0 => Ok(None), + 1 => Ok(Some(inner.as_slice())), + len => { + if *len == self.number_of_frames as usize { + Ok(Some(inner.as_slice())) + } else { + if self.enforce_frame_fg_vm_match { + LengthMismatchVoiLutSnafu { + vm: *len as u32, + nr_frames: self.number_of_frames(), + } + .fail()? + } + tracing::warn!( + "Expected `{:?}` VOI LUTs, found `{:?}`, using first value for all", + self.number_of_frames, + len + ); + Ok(Some(&inner[0..0])) + } + } + } + } else { + Ok(None) + } + } + // converter methods /// Convert the decoded pixel data of a specific frame into a dynamic image. @@ -889,6 +935,8 @@ impl DecodedPixelData<'_> { #[cfg(feature = "image")] fn build_monochrome_image(&self, frame: u32, options: &ConvertOptions) -> Result { + use transform::VoiLutTransform; + let ConvertOptions { modality_lut, voi_lut, @@ -919,11 +967,30 @@ impl DecodedPixelData<'_> { let signed = self.pixel_representation == PixelRepresentation::Signed; - let lut: Lut = match (voi_lut, self.window()?) { - (VoiLutOption::Identity, _) => { + let lut: Lut = match (voi_lut, self.window()?, self.voi_lut_sequence()?) + { + (VoiLutOption::Identity, _, _) => { Lut::new_rescale(8, false, rescale).context(CreateLutSnafu)? } - (VoiLutOption::Default | VoiLutOption::First, Some(window)) => { + ( + VoiLutOption::Default | VoiLutOption::First, + _, + Some(voi_lut_sequence), + ) => Lut::new_rescale_and_lut( + 8, + signed, + rescale, + VoiLutTransform::new( + if voi_lut_sequence.len() > 1 { + &voi_lut_sequence[frame as usize] + } else { + &voi_lut_sequence[0] + }, + 8, + ), + ) + .context(CreateLutSnafu)?, + (VoiLutOption::Default | VoiLutOption::First, Some(window), _) => { Lut::new_rescale_and_window( 8, signed, @@ -948,8 +1015,10 @@ impl DecodedPixelData<'_> { ) .context(CreateLutSnafu)? } - (VoiLutOption::Default | VoiLutOption::First, None) => { - tracing::warn!("Could not find window level for object"); + (VoiLutOption::Default | VoiLutOption::First, None, None) => { + tracing::warn!( + "Could find neither VOI LUT nor window level for object" + ); Lut::new_rescale_and_normalize( 8, signed, @@ -958,7 +1027,7 @@ impl DecodedPixelData<'_> { ) .context(CreateLutSnafu)? } - (VoiLutOption::Custom(window), _) => Lut::new_rescale_and_window( + (VoiLutOption::Custom(window), _, _) => Lut::new_rescale_and_window( 8, signed, rescale, @@ -977,7 +1046,7 @@ impl DecodedPixelData<'_> { ), ) .context(CreateLutSnafu)?, - (VoiLutOption::Normalize, _) => Lut::new_rescale_and_normalize( + (VoiLutOption::Normalize, _, _) => Lut::new_rescale_and_normalize( 8, signed, rescale, @@ -1060,70 +1129,92 @@ impl DecodedPixelData<'_> { let samples = self.frame_data_ow(frame)?; // use 16-bit precision to prevent possible loss of precision in image - let lut: Lut = match (voi_lut, self.window()?) { - (VoiLutOption::Identity, _) => { - Lut::new_rescale(self.bits_stored, signed, rescale) - } - (VoiLutOption::Default | VoiLutOption::First, Some(window)) => { - Lut::new_rescale_and_window( + let lut: Lut = + match (voi_lut, self.window()?, self.voi_lut_sequence()?) { + (VoiLutOption::Identity, _, _) => { + Lut::new_rescale(self.bits_stored, signed, rescale) + } + ( + VoiLutOption::Default | VoiLutOption::First, + _, + Some(voi_lut_sequence), + ) => Lut::new_rescale_and_lut( self.bits_stored, signed, rescale, - WindowLevelTransform::new( - match self.voi_lut_function()? { - Some(lut) => { - if lut.len() > 1 { - lut[frame as usize] - } else { - lut[0] - } - } - None => VoiLutFunction::Linear, - }, - if window.len() > 1 { - window[frame as usize] + VoiLutTransform::new( + if voi_lut_sequence.len() > 1 { + &voi_lut_sequence[frame as usize] } else { - window[0] + &voi_lut_sequence[0] }, + self.bits_stored, ), - ) - } - (VoiLutOption::Default | VoiLutOption::First, None) => { - tracing::warn!("Could not find window level for object"); - - Lut::new_rescale_and_normalize( + ), + (VoiLutOption::Default | VoiLutOption::First, Some(window), _) => { + Lut::new_rescale_and_window( + self.bits_stored, + signed, + rescale, + WindowLevelTransform::new( + match self.voi_lut_function()? { + Some(lut) => { + if lut.len() > 1 { + lut[frame as usize] + } else { + lut[0] + } + } + None => VoiLutFunction::Linear, + }, + if window.len() > 1 { + window[frame as usize] + } else { + window[0] + }, + ), + ) + } + (VoiLutOption::Default | VoiLutOption::First, None, None) => { + tracing::warn!( + "Could find neither VOI LUT nor window level for object" + ); + + Lut::new_rescale_and_normalize( + self.bits_stored, + signed, + rescale, + samples.iter().copied(), + ) + } + (VoiLutOption::Custom(window), _, _) => { + Lut::new_rescale_and_window( + self.bits_stored, + signed, + rescale, + WindowLevelTransform::new( + match self.voi_lut_function()? { + Some(lut) => { + if lut.len() > 1 { + lut[frame as usize] + } else { + lut[0] + } + } + None => VoiLutFunction::Linear, + }, + *window, + ), + ) + } + (VoiLutOption::Normalize, _, _) => Lut::new_rescale_and_normalize( self.bits_stored, signed, rescale, samples.iter().copied(), - ) - } - (VoiLutOption::Custom(window), _) => Lut::new_rescale_and_window( - self.bits_stored, - signed, - rescale, - WindowLevelTransform::new( - match self.voi_lut_function()? { - Some(lut) => { - if lut.len() > 1 { - lut[frame as usize] - } else { - lut[0] - } - } - None => VoiLutFunction::Linear, - }, - *window, ), - ), - (VoiLutOption::Normalize, _) => Lut::new_rescale_and_normalize( - self.bits_stored, - signed, - rescale, - samples.iter().copied(), - ), - } - .context(CreateLutSnafu)?; + } + .context(CreateLutSnafu)?; #[cfg(feature = "rayon")] { @@ -1741,6 +1832,7 @@ impl DecodedPixelData<'_> { rescale: self.rescale.to_vec(), voi_lut_function: self.voi_lut_function.clone(), window: self.window.clone(), + voi_lut_sequence: self.voi_lut_sequence.clone(), enforce_frame_fg_vm_match: self.enforce_frame_fg_vm_match, } } @@ -1927,6 +2019,7 @@ pub(crate) struct ImagingProperties { pub(crate) number_of_frames: u32, pub(crate) voi_lut_function: Option>, pub(crate) window: Option>, + pub(crate) voi_lut_sequence: Option>, } #[cfg(not(feature = "gdcm"))] @@ -1957,6 +2050,7 @@ impl ImagingProperties { .map(|v| VoiLutFunction::try_from((*v).as_str()).ok()) .collect() }); + let voi_lut_sequence = voi_lut_sequence(obj); ensure!( rescale_intercept.len() == rescale_slope.len(), @@ -2006,6 +2100,7 @@ impl ImagingProperties { number_of_frames, voi_lut_function, window, + voi_lut_sequence, }) } } @@ -2033,6 +2128,7 @@ where number_of_frames, voi_lut_function, window, + voi_lut_sequence, } = ImagingProperties::from_obj(self)?; let transfer_syntax = &self.meta().transfer_syntax; @@ -2085,6 +2181,7 @@ where rescale, voi_lut_function, window, + voi_lut_sequence, enforce_frame_fg_vm_match: false, }); } @@ -2117,6 +2214,7 @@ where rescale, voi_lut_function, window, + voi_lut_sequence, enforce_frame_fg_vm_match: false, }) } @@ -2139,6 +2237,7 @@ where number_of_frames, voi_lut_function, window, + voi_lut_sequence, } = ImagingProperties::from_obj(self)?; let transfer_syntax = &self.meta().transfer_syntax; @@ -2191,6 +2290,7 @@ where rescale, voi_lut_function, window, + voi_lut_sequence, enforce_frame_fg_vm_match: false, }); } @@ -2234,6 +2334,7 @@ where rescale, voi_lut_function, window, + voi_lut_sequence, enforce_frame_fg_vm_match: false, }) } diff --git a/pixeldata/src/lut.rs b/pixeldata/src/lut.rs index 7f3adb2af..f383ed309 100644 --- a/pixeldata/src/lut.rs +++ b/pixeldata/src/lut.rs @@ -13,7 +13,7 @@ use num_traits::{NumCast, ToPrimitive}; use rayon::iter::{IntoParallelIterator, ParallelIterator}; use snafu::{OptionExt, Snafu}; -use crate::{Rescale, WindowLevelTransform}; +use crate::{transform::VoiLutTransform, Rescale, WindowLevelTransform}; /// The LUT could not be created: /// entry #{index} was mapped to {y_value}, @@ -251,6 +251,39 @@ where }) } + /// Create a new LUT containing the modality rescale transformation + /// and the VOI transformation defined by an explicit VOI LUT. + /// + /// The amplitude of the output values + /// goes from 0 to `2^n - 1`, where `n` is the power of two + /// which follows `bits_stored` (or itself if it is a power of two). + /// For instance, if `bits_stored` is 12, the output values + /// will go from 0 to 65535. + /// + /// - `bits_stored`: + /// the number of bits effectively used to represent the sample values + /// (the _Bits Stored_ DICOM attribute) + /// - `signed`: + /// whether the input sample values are expected to be signed + /// (_Pixel Representation_ = 1) + /// - `rescale`: the rescale parameters + /// - `voi`: the value of interest (VOI) function and parameters + /// + /// # Panics + /// + /// Panics if `bits_stored` is 0 or too large. + pub fn new_rescale_and_lut( + bits_stored: u16, + signed: bool, + rescale: Rescale, + voi: VoiLutTransform, + ) -> Result { + Self::new_with_fn(bits_stored, signed, |v| { + let v = rescale.apply(v); + voi.apply(v) + }) + } + /// Create a new LUT containing /// a VOI transformation defined by a window level. /// @@ -281,6 +314,34 @@ where Self::new_with_fn(bits_stored, signed, |v| voi.apply(v, y_max)) } + /// Create a new LUT containing a VOI transformation defined by an explicit + /// VOI LUT. + /// + /// The amplitude of the output values + /// goes from 0 to `2^n - 1`, where `n` is the power of two + /// which follows `bits_stored` (or itself if it is a power of two). + /// For instance, if `bits_stored` is 12, the output values + /// will go from 0 to 65535. + /// + /// - `bits_stored`: + /// the number of bits effectively used to represent the sample values + /// (the _Bits Stored_ DICOM attribute) + /// - `signed`: + /// whether the input sample values are expected to be signed + /// (_Pixel Representation_ = 1) + /// - `voi`: the value of interest (VOI) function and parameters + /// + /// # Panics + /// + /// Panics if `bits_stored` is 0 or too large. + pub fn new_lut( + bits_stored: u16, + signed: bool, + voi: VoiLutTransform, + ) -> Result { + Self::new_with_fn(bits_stored, signed, |v| voi.apply(v)) + } + /// Apply the transformation to a single pixel sample value. /// /// Although the input is expected to be one of `u8`, `u16`, or `u32`, diff --git a/pixeldata/src/transform.rs b/pixeldata/src/transform.rs index 2d04b2cb7..2ad657c79 100644 --- a/pixeldata/src/transform.rs +++ b/pixeldata/src/transform.rs @@ -2,6 +2,8 @@ use snafu::Snafu; +use crate::attribute::VoiLut; + /// Description of a modality rescale function, /// defined by a _rescale slope_ and _rescale intercept_. #[derive(Debug, Copy, Clone, PartialEq)] @@ -189,6 +191,53 @@ fn window_level_sigmoid(value: f64, window_width: f64, window_center: f64, y_max y_max / (1. + f64::exp(-4. * (value - wc) / ww)) } +/// A full description of a VOI LUT function transformation +/// based on an explicit LUT. +pub struct VoiLutTransform<'a> { + lut: &'a VoiLut, + shift: u32, +} + +impl<'a> VoiLutTransform<'a> { + /// Create a new LUT transformation. + /// + /// # Panics + /// + /// Panics if `lut.data` is empty, or if `bits_stored` is smaller than + /// `lut.bits_stored`. + pub fn new(lut: &'a VoiLut, bits_stored: u16) -> Self { + if lut.data.is_empty() { + // we'll do unchecked accesses to the first and last items in apply() + panic!("LUT data is empty"); + } + + if bits_stored < (lut.bits_stored as u16) { + panic!( + "LUT with BitsStored {} cannot be used for an image with BitsStored {}", + lut.bits_stored, bits_stored + ); + } + + let shift = (bits_stored.next_power_of_two() - (lut.bits_stored as u16)) as u32; + + VoiLutTransform { lut, shift } + } + + /// Apply the LUT transformation on a rescaled value. + pub fn apply(&self, value: f64) -> f64 { + let value_rounded = value.round().clamp(i32::MIN as f64, i32::MAX as f64) as i32; + let y = (*(if value_rounded <= self.lut.min_pixel_value { + self.lut.data.first().unwrap() + } else { + self.lut + .data + .get((value_rounded - self.lut.min_pixel_value) as usize) + .unwrap_or(self.lut.data.last().unwrap()) + })) << self.shift; + y as f64 + } +} + #[cfg(test)] mod tests { use super::*;