Convert an ECG in dicom waveform format to PDF. Obsolete, see instead: https://github.com/marcodebe/dicomecg_convert

dicom-ecg2pdf

########################################################################
#                                                                      #
# Obsolete, see instead: https://github.com/marcodebe/dicomecg_convert #
#                                                                      #
########################################################################










from __future__ import division
from datetime import datetime
import dicom
import numpy as np
import matplotlib
matplotlib.use('Agg')

import matplotlib.pyplot as plt
import locale
locale.setlocale(locale.LC_TIME, "it_IT.utf8")

from scipy.signal import butter, lfilter


def butter_bandpass(lowcut, highcut, fs, order=5):
    nyquist_freq = 0.5 * fs
    low = lowcut / nyquist_freq
    high = highcut / nyquist_freq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y


def norm(x, max_value):
    sgn = lambda x: lambda y: x-y > 0 and 1 or -1
    norm = lambda y: lambda x: x % (sgn(y)(x) * y)
    return norm(max_value)(x)


def _signals(sequence_item):
    """
    sequence_item := dicom.dataset.FileDataset.WaveformData[n]
    Return a list of signals.
    """

    channel_definitions = sequence_item.ChannelDefinitionSequence
    wavewform_data = sequence_item.WaveformData
    channels_no = sequence_item.NumberOfWaveformChannels

    signals = []
    factor = np.zeros(channels_no) + 1
    for idx in range(channels_no):
        signals.append([])

        definition = channel_definitions[idx]

        if definition.get('ChannelSensitivity'):
            factor[idx] = (float(definition.ChannelSensitivity) *
                           float(
                               definition.ChannelSensitivityCorrectionFactor))

    for idx in xrange(0, len(wavewform_data), 2):
        channel = int((idx/2) % channels_no)
        signals[channel].append(factor[channel] *
                                norm(ord(wavewform_data[idx]) +
                                     256 * ord(wavewform_data[idx+1]), 2**15))

    for i, s in enumerate(signals):
        low = .05
        high = 40.0

        # micro to milli volt
        signals[i] = butter_bandpass_filter(np.asarray(s),
                                            low, high, 1000, order=1)/1000

    return signals


def graphpaper():
    """
    Draw a "graph paper"-like grid on plt figure
    """

    hl01 = np.arange(0, frame_height)
    hl05 = np.arange(0, frame_height, 5)
    hl10 = np.arange(0, frame_height, 10)
    vl01 = np.arange(0, samples*(1+1/frame_width), samples/frame_width)
    vl05 = np.arange(0, samples*(1+1/frame_width), 5*samples/frame_width)
    vl10 = np.arange(0, samples*(1+1/frame_width), 10*samples/frame_width)

    plt.vlines(vl01, 0, frame_height, color='red', linewidth=0.2, alpha=0.3)
    plt.vlines(vl05, 0, frame_height, color='red', linewidth=0.4, alpha=0.6)
    plt.vlines(vl10, 0, frame_height, color='red', linewidth=0.6)
    plt.hlines(hl01, 0, samples, color='red', linewidth=0.2, alpha=0.3)
    plt.hlines(hl05, 0, samples, color='red', linewidth=0.4, alpha=0.6)
    plt.hlines(hl10, 0, samples, color='red', linewidth=0.6)


if __name__ == '__main__':

    ecg_dicom = dicom.read_file('ecg.dcm')
    sequence = ecg_dicom.WaveformSequence
    sequence_item = sequence[0]
    samples = sequence_item.NumberOfWaveformSamples
    frequency = int(sequence_item.SamplingFrequency)
    time = 1. / frequency * np.arange(samples)
    start_time = 0
    end_time = 1./frequency*samples

    time = 1. / frequency * np.arange(samples)
    start_time = 0
    end_time = 1./frequency*samples

    # mm
    a4 = 297, 210
    margin_left = 27
    margin_right = 20
    margin_top = 30
    margin_bottom = 10

    frame_height = a4[1] - margin_bottom - margin_top
    frame_width = a4[0] - margin_right - margin_left

    pat_name = ecg_dicom.PatientName.replace('^', ' ')
    pat_id = ecg_dicom.PatientID
    ecg_date_str = (ecg_dicom.InstanceCreationDate +
                    ecg_dicom.InstanceCreationTime)
    ecg_date = datetime.strftime(datetime.strptime(ecg_date_str,
                                                   '%Y%m%d%H%M%S'),
                                 '%d %b %Y %H:%M')

    text_y = frame_height+15

    plt.text(0, text_y-5, "id: " + pat_id, fontsize=10)
    plt.suptitle(ecg_date)

    left = margin_left / a4[0]
    right = (a4[0] - margin_right) / a4[0]
    bottom = margin_bottom / a4[1]
    top = (a4[1] - margin_top) / a4[1]

    print "left: %s, right: %s" % (left, right)
    #plt.subplots_adjust(left=left, right=right, top=top, bottom=bottom)

    graphpaper()

    plt.ylim(0, frame_height)
    plt.xlim(0, samples-1)

    frame = plt.gca()
    frame.axes.get_xaxis().set_visible(False)
    frame.axes.get_yaxis().set_visible(False)

    fig = plt.gcf()
    # A4R 297 x 210 mm -> 11.69 x 8.27 "
    fig.set_size_inches(11.69, 8.27)
    premax = delta = premin = 0

    signals = _signals(sequence_item)
    signals.reverse()

    delta = 3
    for num, signal in enumerate(signals):
        delta += 1 + 10 * (premax - signal.min())
        plt.plot(10.0*signal+delta, linewidth=0.6, color='black',
                 antialiased=True)
        premax = signal.max()
        premin = signal.min()

    plt.savefig('ecg.pdf')