maandag 25 januari 2016

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

Geen opmerkingen:

Een reactie posten