MRI - 16 - Eerste leerresultaten

Het duurt even voordat de data de juiste vorm heeft om gebruikt te worden. Ergens in Keras worden 4 dimensies verwacht terwijl er maar 3 zijn gegeven.

'Wrong number of dimensions: expected 4, got 3 with shape (1, 96, 96)

Vaak gaat het om de verschillende kleuren in verschillende kanalen. Maar eens proberen om er 1 toe te voegen. Ik vind dit statement:

ModImaX = np.expand_dims(ModImaX, axis=1)

En ja! Het werkt! Nu kijken naar de resultaten. Met 1 epoch valt het best wel tegen. Het is toch tijd om even wat tv te gaan kijken dus ik zet in op 100 epochs. De 'loss' duikt van 23.8 naar 5.94. Lijkt een goede verbetering. De evaluatie score gaat naar 10.58. Ik weet niet precies wat dat bij Keras inhoud. Bij een classificatie netwerk kan je de juist gescoorde tegen de fout gescoorde berekenen. Bij een regressie score zoals hier zal je iets moeten met de gemiddelde kwardratische afwijking of zo-iets. Nog eens even uitzoeken. De afbeeldingen zijn soms aardig maar ook soms nog erg slecht. Morgen eens kijken of we daar verbetering in kunnen krijgen.

 Na 500 epochs lijkt het resultaat alweer beter hoewel er ook nog grote fouten in blijven. De 'loss' is van 5.94 naar 4.74 gedaald. Lijkt niet zo spectaculair maar de daling wordt steeds moeilijker bereikbaar. De evaluatiescore daalt van 10.58 naar 5.76. Deze score wordt berekend bij evaluatie van de 'test-data'. Dat lijkt een enorme verbetering. Ik denk dat in ieder geval zoveel mogelijk extra data aangeboden moet gaan worden. Wordt weer een uurtje 'clicken'.  

MRI - 15 - Betere puntselectie

De eerdere 2 punt-selectie is waarschijnlijk wat minder betrouwbaar. Ik bouw het om naar een 3 punts selectie:

Het verschil is dat deze twee aanhechtingspunten van de hartkleppen wat beter te selecteren zijn dan, zoals in de eerdere opzet 'het midden daartussen'. Uiteraard kan het midden alsnog, indien noodzakelijk, worden uitgerekend. Ik besluit ook zowel de 2 chamber view als de 4 chamber view te gebruiken. Inmiddels al 3 x 500 punten als zodanig opnieuw geselecteerd.  Ook bouw ik het tweede programma om dat deze afbeeldingen moet muteren naar een heleboel afbeeldingen. Moet uiteraard nu 3 x X,Y coordinaten mee muteren. Waarschijnlijk kan Python dat eenvoudig voor een onbepaald aantal punten maar daar wil ik nu geen tijd aan besteden.

Het wordt tijd om de data aan een neuraal netwerk aan te bieden. Het zijn nu ongeveer 2800 afbeeldingen. Bij de Kaggle competitie is er toevallig een Keras tutorial extra aangemaakt. Ik besluit dat platform maar weer eens van stal te halen. Na alweer de nodige updates van Keras en Theano (gaat ook wel erg snel - als die updates!!!) lijkt het programma te werken. Helaas krijg ik een vervelende melding:

ValueError: setting an array element with a sequence.

Ik snap er weinig van. Het lijkt erop dat het combineren van de afbeelding en de punt-coordinaten in een numpy array mij nu parten speelt. Ik besluit ze in aparte arrays aan te maken. De afbeeldingen maak ik 96 x 96 pixels. Tevens maak ik van de coordinaten float-types zodat ze zo'n groot mogelijke nauwkeurigheid behouden.  Nu nog maar eens proberen ....  

MRI - 14 - Veeeeel data

De vorige keer ben ik bezig geweest met het verrijken van de 2ch view met 2 punten. Een 'bovenin' de LV tussen de kleppen en een 'onderin', bij de apex. Ik heb er al ongeveer 300 gedaan en zal daar nog wel even mee door gaan. Het is echter ruim onvoldoende om een deep learning netwerk mee te 'voeden'. De truc is dan ook om er meer data van te maken. Door alle beschikbare afbeeldingen in verschillende hoeveelheden te schuiven, roteren, knippen, schalen kunnen er in principe heel veel afbeeldingen bij worden gemaakt. Zaak is dan wel om ook de 2 'points of interest (poi's)' op dezelfde wijze mee te muteren. Dat lijkt eenvoudig maar blijkt in praktijk (voor mij) erg lastig. Voor elke mutatie functie duurt het veel langer dan verwacht om deze 'bug-vrij' te krijgen.

Uiteindelijk lukt het vrijwel goed. In de 'crop-functie' lijkt nog wel in enkele gevallen een fout te zitten. Heeft vermoedelijk te maken met een crop grens buiten de afbeelding. Ook blijken af en toe de poi's buiten de afbeelding te vallen. Ik maak een extra functie die een laatste kwaliteits-test doet op de afbeelding.
Het lijkt nu goed te werken. Ik kan nu makkelijk van mijn 300 uitgangsafbeeldingen een veelvoud maken. Tevens zijn deze dan al 'op maat gemaakt' voor het DL programma.

'''

execfile('/Users/DWW/Documents/_EF Kaggle/Create multiple images.py')

'''

import os

import cv2

import dicom

import random

import numpy as np

from scipy import ndimage

from sklearn.externals import joblib

from math import cos, sin, radians

def reImage(ImIn):

    f = 100.0 # Bepaald afrondings klasses Hoe groter f Hoe minder klasses.

    return np.round(ImIn/f) * f

def NormIm(InIm): # Verwijder outliers en normalize naar waarde.

    Perc = 95  # alles boven Perc% wordt op maximum gezet

    Max = 500.

    Percentiel = np.percentile(InIm, Perc)

    h=np.clip(InIm, 0, Percentiel)

    return h * Max / h.max() # breng alles naar schaal 'Max'

def WhitenIm(ImIn):

    hx,wx = ImIn.shape

    #print ImIn[h-1][w-1]

    si = 20

    m = np.average(ImIn)

    ImIn[0:si,0:wx] = 0

    ImIn[hx-si:hx,0:wx] = 0

    ImIn[0:hx,0:si] = 0

    ImIn[0:hx,wx-si:wx] = 0

    

    for wi in range(0,wx-si,si):

        for ho in range(0,hx-si,si):

            if np.average(ImIn[ho:ho+si,wi:wi+si]) < m:

                ImIn[ho:ho+si,wi:wi+si] = 0.

    return ImIn

def SchrinkIm(InIm, size, kx,ky,px,py):

    # Function for schrinking or expanding an (square) image to width 'size'

    # Also the point coordinates are corrected simultaniously

    l, b = InIm.shape

    factor = 1.0 * size / b

    kx = factor * kx

    px = factor * px

    ky = factor * ky

    py = factor * py

    if factor < 1:

        Im = cv2.resize(InIm,None, fx = factor, fy = factor,interpolation = cv2.INTER_CUBIC)

    else:

        Im = cv2.resize(InIm,None, fx = factor, fy = factor,interpolation = cv2.INTER_AREA)

    return Im, int(kx),int(ky),int(px),int(py)

def CropIm(InIm, x,y , size, kx,ky,px,py):

    # Function for cutting-out (cropping) a part of the image. Corrected for borders.

    # Also the point coordinates are corrected simultaniously

    # InIm is image in, x,y is center coordinates, size = lengte / breedte in pixels

    # kx,ky, px,py zijn de puntcoordinaten. Deze moeten alle bewerkingen meekrijgen

    l, b = InIm.shape

    hsize = int(size /2)

    xn = x - hsize                          # top left point

    yn = y - hsize                          # top left point

    if xn < 0:

        xn = 0

    if (xn + size) > b:

        xn = b - size

    if yn < 0:

        yn = 0

    if (yn + size) > l:

        yn = l - size

    kx -= xn

    px -= xn

    ky -= yn

    py -= yn

    Im = InIm[ yn:yn+size, xn:xn+size]

    return Im, int(kx),int(ky),int(px),int(py)

def Rotate(x,y, angle, origin=(0, 0)):

    cos_theta, sin_theta = cos(radians(angle)), sin(radians(angle))

    x0, y0 = origin

    x, y = x - x0, y - y0

    return (x * cos_theta - y * sin_theta + x0,

            x * sin_theta + y * cos_theta + y0)

def RotIm(InIm, x,y, angle, kx,ky,px,py):

    # Function to rotate image around x,y with angle 'angle'

    # Also the point coordinates are corrected simultaniously

    rows, cols = InIm.shape

    kx,ky = Rotate(kx,ky, -1 * angle, origin=(x,y))

    px,py = Rotate(px,py, -1 * angle, origin=(x,y))

    M = cv2.getRotationMatrix2D((x,y),angle,1)

    Im = cv2.warpAffine(InIm,M,(cols,rows))

    return Im, int(kx),int(ky),int(px),int(py)

def ShiftIm(InIm, xshift, yshift, kx,ky,px,py):

    # Function for shifting the image.

    # Also the point coordinates are corrected simultaniously

    le, br = InIm.shape

    if xshift < 0:

        fromx = 0

        tox = br + xshift

    else:

        fromx = xshift

        tox = br

        kx -= xshift

        px -= xshift

    if yshift < 0:

        fromy = 0

        toy = le + yshift

    else:

        fromy = yshift

        toy = le

        ky -= yshift

        py -= yshift

    Im = InIm[fromy: toy, fromx: tox]

    return Im, int(kx),int(ky),int(px),int(py)

def TestIm(img, size, kx,ky,px,py):

    # Test if image is size*size and points of interest are in range eof image

    test = True

    b,l = img.shape

    if b != size or l != size:

        test = False

    if kx < 0 or ky < 0 or px < 0 or py < 0:

        test = False

    if kx > size or ky > size or px > size or py > size:

        test = False

    return test

def ShIm(InIm):

    return InIm /500

def draw_circle(event,x,y,flags,param):

    global clickprint

    global Kleppen

    global MouseX

    global MouseY

    if event == cv2.EVENT_LBUTTONDOWN:

        clickprint = True

        if Kleppen:

            Kleppen = False

        else:

            Kleppen = True

        MouseX = x

        MouseY = y

        cv2.circle(img,(x,y),10,(255,0,0),1)

##############################################

# Create additional images for learning      #

##############################################

RootDir = '/Users/DWW/Data/validate/'

fname = '/Users/DWW/Data/Hearthpoints.pkl'

Wname = 'Hearth'

cv2.namedWindow(Wname)

cv2.destroyAllWindows()

Hearthpos = np.array(joblib.load(fname))

ModIma = np.array([])

c = 0

for i in Hearthpos:

    c += 1

    if c == 10:

        break

    ds = dicom.read_file(i[0])

    di = ds.pixel_array

    img = di.copy()

    b, l = img.shape

    x, y = (int(i[1]) + int(i[3]))/2, (int(i[2]) + int(i[4]))/2 # Center of interest points.

    for j in range(3):

        img = di.copy()

        b, l = img.shape

        # Shift

        shiftx = int(b / 6 - b / 3 * random.random())

        shifty = int(l / 6 - l / 3 * random.random())

        img, kx,ky,px,py = ShiftIm(img, shiftx,shifty, int(i[1]),int(i[2]),int(i[3]),int(i[4]) )

        print 'shift:',img.shape

        # Rotate

        angle = 360*random.random()

        x, y = (kx + px)/2, (ky + py)/2 # New center of point of interest

        img, kx,ky,px,py = RotIm(img,x,y,angle, kx,ky,px,py)

        print 'rotate:',img.shape

        # Resize

        b, l = img.shape

        resize = int(b * (.8 + 0.4 * random.random()))

        img, kx,ky,px,py = SchrinkIm(img,resize, kx,ky,px,py)

        print 'resize:',img.shape

        # Crop

        x, y = (kx + px)/2, (ky + py)/2 # New center of point of interest

        crop =  (0.5 + 0.5* random.random()) * b

        img, kx,ky,px,py = CropIm(img, x,y, crop, kx,ky,px,py)

        print 'crop',img.shape

        # Standard size

        resize = 300

        img, kx,ky,px,py = SchrinkIm(img,resize, kx,ky,px,py)

        xs = np.array([img, kx,ky,px,py])

        print 'standard',img.shape

        if TestIm(img, resize, kx,ky,px,py):

            ModIma = np.vstack([ModIma, xs]) if ModIma.size else xs

joblib.dump(ModIma, '/Users/DWW/Data/ModIma.pkl')

for i in ModIma:

    cv2.namedWindow('image')

    cv2.circle(i[0], (i[1],i[2]) ,5, (255,0,0), 2)

    cv2.circle(i[0], (i[3],i[4]) ,5, (255,0,0), 1)

    cv2.imshow("image", ShIm(NormIm(i[0])))

    cv2.waitKey()

cv2.destroyAllWindows()          

MRI - 13 - Informatie toevoegen

Of het toegestaan is in de Kaggle competitie weet ik niet. Maar gegeven de voorgaande analyse, waarbij blijkt dat het met de regelmaat van de locatie van de slices slecht is gesteld, heb ik behoefte om te weten welke slice ik als eerste moet gebruiken. De 2ch en de 4ch view lijken daar mogelijkheden voor te bieden. Ik wil dan een DL leren om deze startpunten en eindpunten zelfstandig te vinden. Op basis daarvan kan ik dan hopelijk de van toepassing zijnde slices selecteren. Het vraagt wel een herverdieping in de 3d meetkunde. Hoe bereken je de coordinaten ook alweer van een punt in een vlak waarvan 3 punten bekend zijn?  

Als ik het goed begrijp worden in de Dicom metadata punten O, R en S gegeven. Mijn tussen doel is om van punt P de coordinaten te berekenen op basis van de L en B welke ik uit kan rekenen door de pixel-coordinaten te vermenigvuldigen met de, ook gegeven, 'PixelSpacing' (Aantal mm per pixel).
Het 2e doel is dan op basis van de P coordinaten het dichtstbijzijnde sax - vlak te vinden. Hoe moet dat ook alweer? Vanuit daar kan dan het volume worden opgebouwd. Natuurlijk moet dan ook de oppervlaktes van het LV in de sax-vlakken nog worden bepaald.

Maar goed, dat komt later wel. Eerst maar eens de P locaties zien te vinden. Met muisklikken selecteer ik de locatie tussen uiteinden van de LV spieren en de punt van het LV. Ik sla het op in een bestand dat ik voor DL wil gaan gebruiken.

Kleine witte circeltjes om de selectiepunten.

Het is een hoop werk en het blijkt ook wel weer dat ofwel de harten zeer gevarieerd zijn ofwel (maar waarschijnlijk en) de opnames niet altijd even juist worden gemaakt. Zie bijvoorbeeld deze:

Er zijn ook nog veel onduidelijker. Ik hoop dat straks het DL algoritme er iets van weet te maken. Dit is de gebruikte code:

'''

execfile('/Users/DWW/Documents/Create 2ch shapes.py')

'''

import os

import cv2

import dicom

import random

import numpy as np

from scipy import ndimage

from sklearn.externals import joblib

def reImage(ImIn):

    f = 100.0 # Bepaald afrondings klasses Hoe groter f Hoe minder klasses.

    return np.round(ImIn/f) * f

def NormIm(InIm): # Verwijder outliers en normalize naar waarde.

    Perc = 95  # alles boven Perc% wordt op maximum gezet

    Max = 500.

    Percentiel = np.percentile(InIm, Perc)

    h=np.clip(InIm, 0, Percentiel)

    return h * Max / h.max() # breng alles naar schaal 'Max'

def WhitenIm(ImIn):

    hx,wx = ImIn.shape

    #print ImIn[h-1][w-1]

    si = 20

    m = np.average(ImIn)

    ImIn[0:si,0:wx] = 0

    ImIn[hx-si:hx,0:wx] = 0

    ImIn[0:hx,0:si] = 0

    ImIn[0:hx,wx-si:wx] = 0

    

    for wi in range(0,wx-si,si):

        for ho in range(0,hx-si,si):

            if np.average(ImIn[ho:ho+si,wi:wi+si]) < m:

                ImIn[ho:ho+si,wi:wi+si] = 0.

    return ImIn

def CropIm(InIm, x,y,pix,size):

    pi = int(size/pix)

    #print 'pi=',pi

    x = x-pi/2

    if x < 0:

         x = 0

    y = y-pi/2

    if y<0:

        y = 0

    h,w = InIm.shape

    if x + pi > h:

        x = h - pi

    if y + pi > w:

        y = w - pi

    NormSize = 300#106

    if pi < NormSize:

        return cv2.resize(InIm[x:x+pi, y:y+pi],(NormSize, NormSize),interpolation = cv2.INTER_CUBIC)

    if pi > NormSize:

        return cv2.resize(InIm[x:x+pi, y:y+pi],(NormSize, NormSize),interpolation = cv2.INTER_AREA)

    if pi == NormSize:

        return InIm[x:x+pi, y:y+pi]

def ShIm(InIm):

    return InIm /500

def draw_circle(event,x,y,flags,param):

    global clickprint

    global Kleppen

    global MouseX

    global MouseY

    if event == cv2.EVENT_LBUTTONDOWN:

        clickprint = True

        if Kleppen:

            Kleppen = False

        else:

            Kleppen = True

        MouseX = x

        MouseY = y

        cv2.circle(img,(x,y),10,(255,0,0),1)

##############################################

# Select points of interest with mouseclicks #

##############################################

RootDir = '/Users/DWW/Data/validate/'

fname = '/Users/DWW/Data/Hearthpoints.pkl'

#joblib.dump([], fname)

cv2.destroyAllWindows()

Wname = 'Hearth'

cv2.namedWindow(Wname)

Hearthpos = np.array([]) # Array filled with point data

Wsize = 800 # Standard display width

Number = 0 # Number of tests

Stop = False

for root, dirs, files in os.walk(RootDir):

    ra = random.random()

    if not(Stop) and "2ch" in root:

        for name in files:

            ra = random.random()

            if not(Stop) and ra < .01 and ".dcm" in name:

                Number += 1

                ds = dicom.read_file(os.path.join(root, name))

                di = ds.pixel_array # afbeelding in numpy array laden

                

                img = ShIm(NormIm(di))

                height, width = img.shape[:2]

                

                Factor = 1.0 * Wsize / width

                img = cv2.resize(img,None,fx=Factor, fy=Factor, interpolation = cv2.INTER_CUBIC)

                Kleppen = True

                cv2.setMouseCallback(Wname,draw_circle)

                clickprint = True

                MouseX = 0

                MouseY = 0

                while(1):

                    if clickprint:

                        if Kleppen:

                            print 'Kleppen'

                        else:

                            print 'Punt'

                        clickprint = False

                    if not(Kleppen):

                        KlepX = int(MouseX/Factor)

                        KlepY = int(MouseY/Factor)

                    else:

                        PuntX = int(MouseX/Factor)

                        PuntY = int(MouseY/Factor)

                    cv2.imshow(Wname, img)

                    cv2.moveWindow(Wname, 100, 150)

                    k = cv2.waitKey(33)

                    if k == 27:

                        break

                    if k == ord('s'):

                        Stop = True

                        break

                xs = np.array([os.path.join(root, name), KlepX, KlepY, PuntX, PuntY])

                #print xs

                Hearthpos = np.vstack([Hearthpos,xs]) if Hearthpos.size else xs

                if Stop:

                    break

#cv2.waitKey()

cv2.destroyAllWindows()

Hearthpos2 = np.array(joblib.load(fname))

Hearthpos = np.vstack([Hearthpos2, Hearthpos]) if Hearthpos2.size else Hearthpos

joblib.dump(Hearthpos, fname)

Hearthpos = np.array(joblib.load(fname))

Hearthpos = Hearthpos[Hearthpos[:,4] !='0']

c = 0

for i in Hearthpos:

    c += 1

    '''

    ds = dicom.read_file(i[0])

    di = ds.pixel_array

    cv2.namedWindow('im')

    img = ShIm(NormIm(di))

    

    cv2.circle(img, (int(i[1]),int(i[2])) ,5, (255,0,0), 2)

    cv2.circle(img, (int(i[3]),int(i[4])) ,5, (255,0,0), 1)

    cv2.imshow("image", img)

    '''

    print c, i[0], i[1],i[2],i[3],i[4]

    #cv2.waitKey()

cv2.destroyAllWindows()