Testing DICOMizing Croaker acquisitions

raw_data/rebuild.py

import sys
import numpy as np
from operator import itemgetter, attrgetter
import Image
from math import *
import math

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.mlab as mlab

# DICOM
import dicom, dicom.UID
from dicom.dataset import Dataset, FileDataset
import datetime, time

Tableau = []
with open("raw_data/20161217-222737-total.csv", 'r') as echOpenLog:
    for line in echOpenLog:
        if(line.startswith("#")):
        #print "commentaire"        
        ikl = 0
        else:
        line = line.strip().split(';')    
        Tableau.append(line)
Tableau = np.asarray(Tableau, dtype=np.float32)
NbOfLines,PointsPerLine = (np.shape(Tableau))
plt.imshow((np.transpose(Tableau) ), cmap="gray")
plt.show()
Raw = Tableau*65535/np.amax(Tableau)

maxAngle = 30.0
ImagePoints = Tableau
step = 2.0*maxAngle/(NbOfLines+1)
size = (NbOfLines,PointsPerLine)

X=np.zeros(shape=(size[0],size[1]))
Y=np.zeros(shape=(size[0],size[1]))
for i in range(size[0]):
    for j in range(size[1]):
    X[i][j] = j*math.cos(math.radians(step*i-maxAngle)) 
    Y[i][j] = (size[1]+1)/2.0+j*math.sin(math.radians(step*(i)-maxAngle)) # same
print int(size[1]*math.cos(math.radians(maxAngle)))


MaxDepth = int(size[1]*math.cos(math.radians(30)))

sizeSC = (size[1],size[1])
ScanConverted=np.zeros(shape=(size[1],size[1]))
ScanConverted = ScanConverted + 40

pix = np.zeros(shape=(size[1],size[1]))
print sizeSC
for i in range(size[1]):
    for j in range(size[1]):
        #value = np.amax(ImagePoints)
        pix[i,j] = 0

for i in range(size[1]):
    if (i<MaxDepth):
        if (i>=0 & i<=(size[1])):
            sweep = int(i*(254/2)/220)
            for j in range((size[1]/2-sweep),(size[1]/2+sweep)):
            D = (X-i)**2 + (Y-j)**2
            resul = np.unravel_index(D.argmin(), D.shape)
            # here is a basic NN, not even a 2-tap
            ScanConverted[i][j] = ImagePoints[resul[0]][resul[1]]
            #print value
            value = ScanConverted[i][j]
            pix[j,i] = value
        if not(i%50):print i 
    else:
        for j in range(size[1]):
            if ( (i**2) + (j-(size[1]/2))**2 ) < ((size[1])**2 - 1):
                D = (X-i)**2 + (Y-j)**2
                resul = np.unravel_index(D.argmin(), D.shape)
                # here is a basic NN, not even a 2-tap
                ScanConverted[i][j] = ImagePoints[resul[0]][resul[1]]
                value = ScanConverted[i][j]
                #print value
                pix[j,i] = value
        if not(i%50):print i 
pix = np.asarray(pix, dtype=np.float32)
pix = pix*65535/np.amax(pix)

plt.imshow((np.transpose(pix) ), cmap="gray")
plt.show()

220
(255, 255)
0
50
100
150
200
250

x = range(sizeSC[0])
y = range(sizeSC[1])

hf = plt.figure()
ha = hf.add_subplot(111, projection='3d')

X, Y = np.meshgrid(x, y)  # `plot_surface` expects `x` and `y` data to be 2D
ha.plot_surface(X, Y, np.sqrt(pix))

plt.show()

def write_dicom(pixel_array,filename):
    """
    INPUTS:
    pixel_array: 2D numpy ndarray.  If pixel_array is larger than 2D, errors.
    filename: string name for the output file.
    """

    ## This code block was taken from the output of a MATLAB secondary
    ## capture.
    file_meta = Dataset()
    file_meta.MediaStorageSOPClassUID = 'Ultrasound_Image_Storage'
    file_meta.MediaStorageSOPInstanceUID = '1.2.840.10008.5.1.4.1.1.6.1' # Ultrasound
    file_meta.ImplementationClassUID = '1.3.6.1.4.1.9590.100.1.0.100.4.0'
    ds = FileDataset(filename, {},file_meta = file_meta,preamble="\0"*128)
    ds.Modality = 'WSD' # Workstation
    ds.ContentDate = str(datetime.date.today()).replace('-','')
    ds.ContentTime = str(time.time()) #milliseconds since the epoch
    #ds.StudyInstanceUID =  '1.3.6.1.4.1.9590.100.1.1.12431397741236017523427128747280487209311'
    #ds.SeriesInstanceUID = '1.3.6.1.4.1.9590.100.1.1.3692311180110610034034218591726431436492'
    #ds.SOPInstanceUID =    '1.3.6.1.4.1.9590.100.1.1.1111656844110176690217683857207368737803'
    ds.SOPClassUID = 'Ultrasound_Image_Storage'
    ds.SecondaryCaptureDeviceManufctur = 'Croaker - https://github.com/kelu124/echomods/blob/master/croaker/data/20161217/20161217-TestingArduinoAndPhantom.md'
    ## These are the necessary imaging components of the FileDataset object.
    ds.SamplesPerPixel = 1
    ds.PhotometricInterpretation = "MONOCHROME2"
    ds.PixelRepresentation = 0
    ds.HighBit = 15
    ds.BitsStored = 16
    # spatial resolution: see http://dicom.nema.org/medical/dicom/2014c/output/chtml/part03/sect_10.7.html#sect_10.7.1.3
    ds.PixelSpacing = '\\0.44\\0.44' # mm/pixel
    # F acq = 1.715MHz see https://github.com/kelu124/bomanz/blob/master/ADC08200/20170430-ADC08200-FirstAcqs.ipynb
    # then with v = 1500 then 1 sample is (1500/2)*1.715*10^6 ( 437,317784257µm by pixel )
    ds.BitsAllocated = 16
    ds.SmallestImagePixelValue = '\\x00\\x00'
    ds.LargestImagePixelValue = '\\xff\\xff'
    # Creating the image itself
    pixel_array = np.array(np.transpose(pixel_array))
    ds.Columns = pixel_array.shape[1]
    ds.Rows = pixel_array.shape[0]
    if pixel_array.dtype != np.uint16:
        pixel_array = pixel_array.astype(np.uint16)
    ds.PixelData = pixel_array.tostring()
    ds.save_as(filename)
    return filename

write_dicom(pix,"20170502-DICOMizingCroakerData.dcm")

'20170502-DICOMizingCroakerData.dcm'

print repr('\\x0.44\\x0.44')

'\\x0.44\\x0.44'