In [5]:

    

%matplotlib inline


    

In [6]:

    

from __future__ import division

__author__ = 'Diego'

import os
import nibabel as nib
from scipy import ndimage
import braviz
import numpy as np
import skimage
from skimage import data, filter, color
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage.draw import circle_perimeter
from skimage.util import img_as_float
import  matplotlib
#matplotlib.use("Qt4Agg")
import matplotlib.pyplot as plt

os.chdir(r"C:\Users\Diego\Documents\prueba_free_surf_400")

skull = nib.load("skull.nii.gz")
data = skull.get_data()
#===============

#data2 = data

data2 = ndimage.median_filter(data,size=3)
data2 = data2 > 50
data3 = data2 - ndimage.binary_erosion(data2,np.ones((3,3,3)))
slices = np.arange(130,165)

data3 = data3[:,slices,:]
data2 = data2[:,slices,:]

data2 = data2 - np.min(data2)
data2 = data2/np.max(data2)

data3 = data3 - np.min(data2)
data3 = data3/np.max(data2)

tot_scores = []


    

In [7]:

    

for s in xrange(data2.shape[1]):
    print s
    ax = plt.gca()
    ax.clear()
    slice = data2[:,s,:].squeeze()
    e_slice = data3[:,s,:].squeeze()
    image = img_as_float(e_slice)
    edges = img_as_float(e_slice)
    hough_radii = np.arange(4, 15, 1,dtype=np.intp)
    hough_res = hough_circle(edges, hough_radii)

    centers = []
    accums = []
    radii = []
    score = []
    for radius, h in zip(hough_radii, hough_res):
        # For each radius, extract two circles
        peaks = peak_local_max(h, num_peaks=2)
        centers.extend(peaks)
        accums.extend(h[peaks[:, 0], peaks[:, 1]])
        radii.extend([radius, radius])
        coords_scores = [h[tuple(p)] for p in peaks]
        score.extend(coords_scores)


    # Draw the most prominent 5 circles
    #image[image==0]=200
    image = color.gray2rgb(image)
    for idx in np.argsort(accums)[::-1][:2]:
        print accums[idx]
        center_x, center_y = centers[idx]
        radius = radii[idx]
        cx, cy = circle_perimeter(center_y, center_x, radius)
        image[cy, cx] = (20, 200, 20)
    ax.imshow(image, cmap=plt.cm.gray)
    plt.show()
    tot_scores.extend(score)


    

    




0
0.71875
0.6875






    













    




1
0.75
0.65625






    













    




2
0.6875
0.6875






    













    




3
0.8125
0.6875






    













    




4
0.875
0.8125






    













    




5
0.78125
0.78125






    













    




6
0.78125
0.78125






    













    




7
0.78125
0.71875






    













    




8
0.84375
0.75






    













    




9
0.8125
0.78125






    













    




10
0.875
0.84375






    













    




11
0.875
0.78125






    













    




12
0.78125
0.78125






    













    




13
0.9375
0.78125






    













    




14
0.8125
0.78125






    













    




15
1.0
0.96875






    













    




16
0.875
0.84375






    













    




17
0.96875
0.9375






    













    




18
0.875
0.8125






    













    




19
0.84375
0.84375






    













    




20
0.78125
0.78125






    













    




21
0.75
0.736111111111






    













    




22
0.763888888889
0.75






    













    




23
0.9375
0.84375






    













    




24
0.85
0.84375






    













    




25
0.71875
0.65625






    













    




26
0.71875
0.71875






    













    




27
0.84375
0.78125






    













    




28
0.71875
0.71875






    













    




29
0.71875
0.71875






    













    




30
0.71875
0.71875






    













    




31
0.6875
0.65625






    













    




32
0.78125
0.71875






    













    




33
0.75
0.71875






    













    




34
0.8125
0.65625






    

In [7]: