import pydicom as dicom
import matplotlib.pyplot as plt
import numpy as np


def grouped(iterable, n):
    "s -> (s0,s1,s2,...sn-1), (sn,sn+1,sn+2,...s2n-1), (s2n,s2n+1,s2n+2,...s3n-1), ..."
    return zip(*[iter(iterable)] * n)


def bytes2int(bytes_arr):
    return [y << 8 | x for x, y in grouped(bytes_arr.__iter__(), 2)]


def mammo_convert(dcm, which_LUT=0):
    arr = dcm.pixel_array
    lut = dcm.VOILUTSequence[which_LUT]
    lut_data = bytes2int(lut.LUTData)
    first_mapped = bytes2int(lut.LUTDescriptor)[1]
    last_mapped = first_mapped + len(lut_data) - 1

    ret = arr.copy()
    ret[ret > last_mapped] = last_mapped
    ret[ret < first_mapped] = first_mapped
    ret -= first_mapped

    return np.array(lut_data)[ret]


def mammo_convert2(dcm, which_window=0):
    arr = dcm.pixel_array
    wc = dcm.WindowCenter[which_window]
    ww = dcm.WindowWidth[which_window]
    out_range = 2 ** (dcm.BitsStored)

    ret = out_range / (1 + np.exp((-4) * (arr.astype(np.float) - wc) / ww))

    return ret


dc = dicom.read_file(r'C:\temp\1.2.840.113619.2.255.224244498451.3934190211141740.243.dcm', force=True)
dc.file_meta.TransferSyntaxUID = dicom.uid.ImplicitVRLittleEndian
ret = mammo_convert(dc)
ret2 = mammo_convert2(dc)

plt.imshow(ret - np.round(ret2), cmap='gray')
plt.imshow(ret - ret2, cmap='gray')