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()