Pydicom: ValueError when running apply_voi_lut function

[aisosalo]

We are working with mammography data and came across a set of files for which an error is produced when using the pydicom function apply_voi_lut in pixel_data_handlers referring to this line here:

lut_data = np.asarray(item.LUTData, dtype=dtype)

The following error message is given:

Error using _asarray>asarray (line 85)
ValueError: invalid literal for int() with base 10: 
b'\x00\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x00\x01\x

Dataset information does not seem to give indications what might be different with those specific files producing the above error. We are not at the moment able to provide more details about the contents of (0028, 3006) LUT Data in those cases when the error is issued, but we can look this piece of information up later if it helps. (We cannot provide a sample DICOM file, as it is patient data.) As a side note, for majority of our files, apply_voi_lut function works just fine.

Our system configuration is as follows:

Windows-10-10.0.17763-SP0
Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)]
numpy 1.18.1
pydicom 1.4.2

We are running the Python function in question within a Matlab script:

ds = py.pydicom.dcmread(filepath);
dicom_utils = py.importlib.import_module('pydicom.pixel_data_handlers.util');
img = dicom_utils.apply_voi_lut(ds.pixel_array.copy(), ds);
img = img.uint16;

Can you provide guidance?

There has already been a brief discussion about this issue with @darcymason.

[scaramallion]

Please post the full traceback and the group 0x0028 elements.

print(ds[0x00280000:0x00300000])

Is the dataset explicit VR and the VR for (0028,3006) LUT Data OW?

If so, try doing this before calling apply_voi_lut():

from struct import unpack

item = ds.VOILUTSequence[0]
nr_entries = item.LUTDescriptor[0] or 2**16
item.LUTData = unpack('<{}H'.format(nr_entries), item.LUTData)

[aisosalo]

The requested group 0x0028 elements:

(0028, 0002) Samples per Pixel                   US: 1
(0028, 0004) Photometric Interpretation          CS: 'MONOCHROME2'
(0028, 0010) Rows                                US: 2294
(0028, 0011) Columns                             US: 1914
(0028, 0100) Bits Allocated                      US: 16
(0028, 0101) Bits Stored                         US: 12
(0028, 0102) High Bit                            US: 11
(0028, 0103) Pixel Representation                US: 0
(0028, 0120) Pixel Padding Value                 US: 0
(0028, 0121) Pixel Padding Range Limit           US: 670
(0028, 0300) Quality Control Image               CS: 'NO'
(0028, 0301) Burned In Annotation                CS: 'NO'
(0028, 1040) Pixel Intensity Relationship        CS: 'LOG'
(0028, 1041) Pixel Intensity Relationship Sign   SS: -1
(0028, 1050) Window Center                       DS: [2875, 2923, 2797]
(0028, 1051) Window Width                        DS: [900, 750, 1050]
(0028, 1052) Rescale Intercept                   DS: "0.0"
(0028, 1053) Rescale Slope                       DS: "1.0"
(0028, 1054) Rescale Type                        LO: 'US'
(0028, 1055) Window Center & Width Explanation   LO: ['NORMAL', 'HARDER', 'SOFTER']
(0028, 1056) VOI LUT Function                    CS: 'SIGMOID'
(0028, 1300) Breast Implant Present              CS: 'NO'
(0028, 2110) Lossy Image Compression             CS: '00'
(0028, 3010)  VOI LUT Sequence   3 item(s) ---- 
   (0028, 3002) LUT Descriptor                      US: [3243, 847, 12]
   (0028, 3003) LUT Explanation                     LO: 'NORMAL'
   (0028, 3006) LUT Data                            OW: Array of 6486 elements
   ---------
   (0028, 3002) LUT Descriptor                      US: [2849, 1233, 12]
   (0028, 3003) LUT Explanation                     LO: 'HARDER'
   (0028, 3006) LUT Data                            OW: Array of 5698 elements
   ---------
   (0028, 3002) LUT Descriptor                      US: [3661, 431, 12]
   (0028, 3003) LUT Explanation                     LO: 'SOFTER'
   (0028, 3006) LUT Data                            OW: Array of 7322 elements
   ---------

Doing

item = ds.VOILUTSequence[0]
nr_entries = item.LUTDescriptor[0] or 2**16
item.LUTData = unpack('<{}H'.format(nr_entries), item.LUTData)

before calling apply_voi_lut() makes the issue go away.

[scaramallion]

OK thanks for the confirmation, it should be fixed soon

[scaramallion]

Fixed in current master, do you mind testing it to confirm?

[aisosalo]

I tested the new implementation with few files and the solution is working as it should. And it did not break anything at least on our corner either regarding the files for which the function was already working well. Thank you!

I can run a more covering test using hundreds of files within a week or so, if you wish?

I suppose the fix could earn its own function:

def get_lut_data(item, nr_entries, ds)
    # Ambiguous VR, US or OW
    if item['LUTData'].VR == 'OW':
        endianness = '<' if ds.is_little_endian else '>'
        unpack_fmt = '{}{}H'.format(endianness, nr_entries)
        lut_data = unpack(unpack_fmt, item.LUTData)
    else:
        lut_data = item.LUTData

    return lut_data

[aisosalo]

Since we were working with Matlab I wrote a Matlab version of the apply_voi_lut as an intermediate solution for our problem by following as closely as possible the pydicom apply_voi_lut implementation. It worked as such without any significant additions receiving as an input the outcomes of Matlab dicominfo and dicomread.

I was wondering if we could make the Matlab version / conversion also public? To my knowledge, Matlab is lacking a good function for this kind of task. @scaramallion could be named there as one of the authors if you wish, and we can surely reference the pydicom licence according to the requests of the pydicom team and provide a link to the documentation of this repository? How does this sound?

[scaramallion]

I can run more covering test using hundreds of files within a week or so, if you wish?

If you don't mind doing it, it might be useful to see if there are any other cases I've missed.

I was wondering if we could make the Matlab version / conversion also public?

It's fine by me, everything is licensed under MIT anyway so as long as you include the copyright notice you're all good, but a link back to pydicom would be appreciated

[aisosalo]

If you don't mind doing it, it might be useful to see if there are any other cases I've missed.

I run a test with 10000 mammograms. Everything worked fine regarding the apply_voi_lut function.

The issue can now be closed on my behalf. I'll open a new one if we come across additional exceptions to the current rules.

[aisosalo]

It's fine by me, everything is licensed under MIT anyway so as long as you include the copyright notice you're all good, but a link back to pydicom would be appreciated

How does this look?

One line came to me as little bit surprising:

            elif width <= 0:
                raise ValueError(
                    "The (0028,1051) Window Width must be greater than 0 "
                    "for a 'LINEAR_EXACT' windowing operation"
                )

Should there be somewhere and voi_func == 'LINEAR_EXACT' or something like that?

Matlab translation:

function arr_out = dicom_apply_voi_lut(arr_in, nfo, index)
%DICOM_APPLY_VOI_LUT Apply a VOI lookup table or windowing operation to
% image data from a compliant DICOM file
%
% INPUTS
%    arr_in     image data from a compliant DICOM file
%    nfo        metadata from a compliant DICOM file
%    index      the index of the view
%
% OUTPUTS:
%    out_arr    image data with applied VOI LUT or windowing operation
%
% See also DICOMINFO, DICOMREAD

% Author: @aisosalo, 2020-

% Method adapted from the Pydicom library function apply_voi_lut, 
% which is originally licensed under the following license:
% 
% Copyright (c) 2008-2020 Darcy Mason and pydicom contributors
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
%
% See: https://pydicom.github.io/
% See: https://github.com/pydicom/pydicom
% See: https://github.com/pydicom/pydicom/blob/b6709beba54f534418ad621e1832f20f4f5e0af2/pydicom/pixel_data_handlers/util.py#L263

if ~exist('index', 'var') ||  isempty(index)
    index = 1;
end

if sum(strcmp(fieldnames(nfo), 'VOILUTSequence')) == 1
    items = nfo.VOILUTSequence;
    items = struct2cell(items);
    item1 = items{index};

    if ~isempty(item1.LUTDescriptor(1)) && any(item1.LUTDescriptor(1)) % empty or zero
        nr_entries = item1.LUTDescriptor(1);
    else
        nr_entries = 2^16;
    end

    first_map = item1.LUTDescriptor(2);

    nominal_depth = item1.LUTDescriptor(3);

    if any([10, 11, 12, 13, 14, 15, 16] == nominal_depth)
        dtype = 'uint16';
    elseif nominal_depth == 8
        dtype = 'uint8';
    else
        error('NotImplementedError:entryNotSupported',...
            'NotImplementedError. \n%s bits per LUT entry is not supported.',num2str(nominal_depth))
    end

    if strcmp(dtype,'uint16')
        lut_data = item1.LUTData;
        lut_data = uint16(lut_data);
    elseif strcmp(dtype,'uint8')
        lut_data = item1.LUTData;
        lut_data = uint8(lut_data);
    end

    if strcmp(class(arr_in), 'uint16')  %#ok<STISA>
        clipped_iv = uint16(zeros(size(arr_in)));
    elseif strcmp(class(arr_in), 'uint8') %#ok<STISA>
        clipped_iv = uint8(zeros(size(arr_in)));
    else
        error('NotImplementedError:inputArrayNotSupported',...
            'NotImplementedError. \nInput array type %s not supported.',class(arr_in))
    end

    mapped_pixels = arr_in >= first_map;

    clipped_iv(mapped_pixels) = arr_in(mapped_pixels) - first_map;

    clipped_iv(clipped_iv > (nr_entries - 1)) = (nr_entries - 1);

    [~, ~, ind] = unique(clipped_iv);
    clipped_iv_mapped = lut_data(ind);

    arr_out = reshape(clipped_iv_mapped, size(clipped_iv));

elseif sum(strcmp(fieldnames(nfo), 'WindowCenter')) == 1 && sum(strcmp(fieldnames(nfo), 'WindowCenter')) == 1
    if sum(strcmp(nfo.PhotometricInterpretation, 'MONOCHROME1')) == 0 && sum(strcmp(nfo.PhotometricInterpretation, 'MONOCHROME2')) == 0
        error('ValueError:photometricInterpretationNotAllowed',...
            'ValueError. \nOnly MONOCHROME1 and MONOCHROME2 are allowed for (0028,0004) Photometric Interpretation.')
    end

    if sum(strcmp(fieldnames(nfo), 'VOILUTFunction')) == 0 % no such field as VOILUTFunction
        voi_func = 'LINEAR';
    else
        voi_func = upper(nfo.VOILUTFunction); % convert to uppecase
    end

    data_elem_wc = nfo.WindowCenter;

    if isstruct(data_elem_wc) && sum(strcmp(fieldnames(data_elem_wc), 'VM')) == 0 && data_elem_wc.VM > 1
        center = data_elem_wc(index);
    else
        center = data_elem_wc;
    end

    data_elem_ww = nfo.WindowWidth;

    if isstruct(data_elem_ww) && sum(strcmp(fieldnames(data_elem_ww), 'VM')) == 0 && data_elem_ww.VM > 1
        width = data_elem_ww(index);
    else
        width = data_elem_ww;
    end

    % The output range depends on whether or not a modality LUT or rescale
    % operation has been applied
    if sum(strcmp(fieldnames(nfo), 'ModalityLUTSequence')) == 1 % unsigned
        y_min = 0;
        modality_lut_seq = nfo.ModalityLUTSequence;
        lut_desc = modality_lut_seq{1};
        bit_depth = lut_desc{2};
        y_max = int32(2^bit_depth - 1);
    elseif sum(strcmp(fieldnames(nfo), 'PixelRepresentation')) == 1 && nfo.PixelRepresentation == 0 % unsigned
        y_min = 0;
        y_max = int32(2^nfo.BitsStored - 1);
    else % signed
        y_min = -int32(2^(nfo.BitsStored - 1));
        y_max = int32(2^(nfo.BitsStored - 1) - 1);
    end

    if sum(strcmp(fieldnames(nfo), 'RescaleSlope')) == 1 && sum(strcmp(fieldnames(nfo), 'RescaleIntercept')) == 1
        y_min = y_min * nfo.RescaleSlope + nfo.RescaleIntercept;
        y_max = y_max * nfo.RescaleSlope + nfo.RescaleIntercept;
    end

    y_range = y_max - y_min;
    y_range = double(y_range);

    arr_in = double(arr_in);

    if strcmp(voi_func,'LINEAR') || strcmp(voi_func,'LINEAR_EXACT')
        if strcmp(voi_func,'LINEAR')
            if width < 1
                error('ValueError:invalidWindowWidth',...
                    'ValueError. \nWindow Width must be greater than or equal to 1 for a %s windowing operation.', 'LINEAR')
            end    

            center = center - 0.5;
            center = double(center);

            width = width - 1;
            width = double(width);

        elseif width <= 0
            error('ValueError:invalidWindowWidth',...
                    'ValueError. \nWindow Width must be greater than 0 for a %s windowing operation.', 'LINEAR_EXACT')
        end        

        below = arr_in <= (center - width / 2);
        above = arr_in > (center + width / 2);

        between = and(~below,~above);    

        arr_in(below) = y_min;
        arr_in(above) = y_max;

        if any(between)
            arr_in(between) = ((arr_in(between) - center) / width + 0.5) * y_range + double(y_min);
        end

    elseif strcmp(voi_func,'SIGMOID')
        if width <= 0
            error('ValueError:invalidWindowWidth',...
                    'ValueError. \nWindow Width must be greater than 0 for a %s windowing operation.', 'SIGMOID')
        end
        arr_in = y_range ./ (1 + exp(-4 * (arr_in - center) / width)) + double(y_min);
    else
        error('ValueError:invalidWindowWidth',...
                    'ValueError. \nUnsupported (0028,1056) VOI LUT Function value %s.', voi_func)
    end

    arr_out = arr_in;
end

end

[mrbean-bremen]

One line came to me as little bit surprising:

Here is the context:

    if voi_func in ['LINEAR', 'LINEAR_EXACT']:
        if voi_func == 'LINEAR':
            ...      
        elif width <= 0:
            raise ValueError(...)

As you can see, this is only done for the LINEAR_EXACT VOI LUT function.

[scaramallion]

Looks good!