Carotid artery 2D demo: Max-tree Filtering

MainPage

Max-tree: artery wall segmentation

  1 import numpy as np
  2 import iamxt
  3 import ia636
  4 import cv2
  5 from skimage.transform import hough_circle
  6 from skimage.draw import circle_perimeter
  7 from skimage.feature import peak_local_max
  8 import ia870
  9 import time
 10 
 11 #Filtering thresholds
 12 a_min,a_max = 100,400
 13 dx_max,dy_max = 35,35
 14 rr_max = 0.5
 15 ar_min = 0.75
 16 circ_min = 1.0
 17 circ_max = 0.18
 18 ecc_max = 0.6
 19 
 20 t = time.time()
 21 #Loading the image
 22 img = adreadgray(find_attachment_file('homeroberto/carotid_demo/carotid_slice.png'))
 23 H,W = img.shape
 24 print "Image shape: %d x %d"%(H,W)
 25 print
 26 adshow(img,"Original image")
 27 
 28 #Structuring element: C4 connectivity
 29 Bc = np.zeros((3,3), dtype = bool)
 30 Bc[1,:] = True
 31 Bc[:,1] = True
 32 
 33 #Building the max-tree
 34 mxt = iamxt.MaxTreeAlpha(img,Bc)
 35 print "Original"
 36 print "Max-tree"
 37 print "Number of nodes: %d" %mxt.node_array.shape[1]
 38 print "Number of leaves: %d" %((mxt.node_array[1,:] == 0).sum())
 39 
 40 
 41 #Area filtering
 42 area = mxt.node_array[3,:]
 43 mxt.contractDR(np.logical_and(area>a_min,area< a_max))
 44 img2 = mxt.getImage()
 45 adshow(img2,"Only nodes with %d < area < %d" %(a_min,a_max))
 46 
 47 print "Area filtering"
 48 print "Max-tree"
 49 print "Number of nodes: %d" %mxt.node_array.shape[1]
 50 print "Number of leaves: %d" %((mxt.node_array[1,:] == 0).sum())
 51 print
 52 
 53 #Bounding-box filtering
 54 xmin,xmax = mxt.node_array[6,:],mxt.node_array[7,:]
 55 ymin,ymax = mxt.node_array[9,:],mxt.node_array[10,:]
 56 dx = xmax - xmin
 57 dy = ymax - ymin
 58 mxt.contractDR(np.logical_and(dy<dy_max,dx< dx_max))
 59 img3 = mxt.getImage()
 60 adshow(img3,"Only nodes with dx < %d and dy < %d" %(dx_max,dy_max))
 61 
 62 print "Bounding_box filtering"
 63 print "Max-tree"
 64 print "Number of nodes: %d" %mxt.node_array.shape[1]
 65 print "Number of leaves: %d" %((mxt.node_array[1,:] == 0).sum())
 66 print
 67 
 68 
 69 #Aspect ratio filtering
 70 xmin,xmax = mxt.node_array[6,:],mxt.node_array[7,:]
 71 ymin,ymax = mxt.node_array[9,:],mxt.node_array[10,:]
 72 AR = np.asarray([xmax-xmin,ymax - ymin], dtype = float)
 73 AR = 1.0*AR.min(axis = 0)/AR.max(axis = 0)
 74 mxt.contractDR(AR>ar_min)
 75 img4 = mxt.getImage()
 76 adshow(img4,"Only nodes with AR > %f" %ar_min)
 77 
 78 print "Aspect ratio filtering"
 79 print "Max-tree"
 80 print "Number of nodes: %d" %mxt.node_array.shape[1]
 81 print "Number of leaves: %d" %((mxt.node_array[1,:] == 0).sum())
 82 print
 83 
 84 
 85 #Filtering the max-tree by the circularity
 86 n = mxt.node_array.shape[1]
 87 P=np.zeros(n, dtype=np.int)
 88 P[0] = H*2 + W*2 - 4
 89 eroded_img = ia870.iaero (img4, Bc)
 90 node_index = mxt.node_index.copy()
 91 h = mxt.node_array[2,:]
 92 pindex = mxt.node_array[0,:]
 93 for i in xrange(n-1,0,-1):
 94     xmin,xmax = mxt.node_array[6,i], mxt.node_array[7,i] + 1
 95     ymin,ymax = mxt.node_array[9,i], mxt.node_array[10,i] + 1
 96     cc_bb = node_index[xmin:xmax,ymin:ymax]
 97     cc = cc_bb == i
 98     cc_bb[cc] = pindex[i]
 99     grad = cc - (eroded_img[xmin:xmax,ymin:ymax]>=h[i])*cc
100     P[i]= grad.sum()
101 
102 A = mxt.node_array[3,:]
103 circ = 4*np.pi*A/(P**2)
104 mxt.contractDR(circ<circ_max)
105 img5 = mxt.getImage()
106 adshow(img5,"Only nodes with circularity < %f" %circ_max)
107 print "Circularity filtering"
108 print "Max-tree"
109 print "Number of nodes: %d" %mxt.node_array.shape[1]
110 print "Number of leaves: %d" %((mxt.node_array[1,:] == 0).sum())
111 print
112 
113 #Filtering by the eccentricity
114 ecc,L1,L2 = mxt.computeEccentricity()
115 mxt.contractDR(ecc<ecc_max)
116 img6 = mxt.getImage()
117 adshow(img6,"Only nodes with eccentricity < %f" %ecc_max)
118 
119 #Extinction filter
120 A = mxt.node_array[3,:]
121 Aext = mxt.computeExtinctionValues(A,"area")
122 Height = mxt.computeHeight()
123 Hext = mxt.computeExtinctionValues(Height,"height")
124 ext = .25*Hext/Hext.max() + .75*Aext/Aext.max()
125 
126 mxt.extinctionFilter(ext,2)
127 img7 = mxt.getImage()
128 adshow(img7,"Area + Height EF")
129 print "Time: %f s" %(time.time()-t)
Image shape: 512 x 512

Original
Max-tree
Number of nodes: 43368
Number of leaves: 15571
Area filtering
Max-tree
Number of nodes: 977
Number of leaves: 157

Bounding_box filtering
Max-tree
Number of nodes: 720
Number of leaves: 143

Aspect ratio filtering
Max-tree
Number of nodes: 423
Number of leaves: 76

Circularity filtering
Max-tree
Number of nodes: 69
Number of leaves: 35

Time: 0.351501 s

Original image

Only nodes with 100 < area < 400

Only nodes with dx < 35 and dy < 35

Only nodes with AR > 0.750000

Only nodes with circularity < 0.180000

Only nodes with eccentricity < 0.600000

Area + Height EF

Min-tree: inside the carotid wall segmentation

 1 img_neg = 255-img
 2 adshow(img_neg,"Negative of the image")
 3 
 4 
 5 t = time.time()
 6 #Building the min-tree
 7 mxt2 = iamxt.MaxTreeAlpha(img_neg,Bc)
 8 print "Min-tree"
 9 print "Number of nodes: %d" %mxt2.node_array.shape[1]
10 print "Number of leaves: %d" %((mxt2.node_array[1,:] == 0).sum())
11 print
12 
13 #Area filtering
14 area2 = mxt2.node_array[3,:]
15 mxt2.contractDR(np.logical_and(area2>a_min,area2< a_max))
16 img2_neg = mxt2.getImage()
17 adshow(img2_neg,"Only nodes with %d < area < %d" %(a_min,a_max))
18 print "Area filtering"
19 print "Min-tree"
20 print "Number of nodes: %d" %mxt2.node_array.shape[1]
21 print "Number of leaves: %d" %((mxt2.node_array[1,:] == 0).sum())
22 print
23 
24 #Bounding-box filtering
25 xmin2,xmax2 = mxt2.node_array[6,:],mxt2.node_array[7,:]
26 ymin2,ymax2 = mxt2.node_array[9,:],mxt2.node_array[10,:]
27 dx2 = xmax2 - xmin2
28 dy2 = ymax2 - ymin2
29 mxt2.contractDR(np.logical_and(dy2<dy_max,dx2< dx_max))
30 img3_neg = mxt2.getImage()
31 adshow(img3_neg,"Only nodes with dx < %d and dy < %d" %(dx_max,dy_max))
32 print "Bounding_box filtering"
33 print "Min-tree"
34 print "Number of nodes: %d" %mxt2.node_array.shape[1]
35 print "Number of leaves: %d" %((mxt2.node_array[1,:] == 0).sum())
36 print
37 
38 
39 #Aspect ratio filtering
40 xmin2,xmax2 = mxt2.node_array[6,:],mxt2.node_array[7,:]
41 ymin2,ymax2 = mxt2.node_array[9,:],mxt2.node_array[10,:]
42 AR = np.asarray([xmax2-xmin2,ymax2 - ymin2], dtype = float)
43 AR = 1.0*AR.min(axis = 0)/AR.max(axis = 0)
44 mxt2.contractDR(AR>ar_min)
45 img4_neg = mxt2.getImage()
46 adshow(img4_neg,"Only nodes with AR > %f" %ar_min)
47 
48 print "Aspect ratio filtering"
49 print "Min-tree"
50 print "Number of nodes: %d" %mxt2.node_array.shape[1]
51 print "Number of leaves: %d" %((mxt2.node_array[1,:] == 0).sum())
52 print
53 
54 
55 #Filtering the min-tree by the circularity
56 n = mxt2.node_array.shape[1]
57 P=np.zeros(n, dtype=np.int)
58 P[0] = H*2 + W*2 - 4
59 eroded_img = ia870.iaero (img4_neg, Bc)
60 node_index = mxt2.node_index.copy()
61 h = mxt2.node_array[2,:]
62 pindex = mxt2.node_array[0,:]
63 for i in xrange(n-1,0,-1):
64     xmin,xmax = mxt2.node_array[6,i], mxt2.node_array[7,i] + 1
65     ymin,ymax = mxt2.node_array[9,i], mxt2.node_array[10,i] + 1
66     cc_bb = node_index[xmin:xmax,ymin:ymax]
67     cc = cc_bb == i
68     cc_bb[cc] = pindex[i]
69     grad = cc - (eroded_img[xmin:xmax,ymin:ymax]>=h[i])*cc
70     P[i]= grad.sum()
71 
72 A = mxt2.node_array[3,:]
73 circ = 4*np.pi*A/(P**2)
74 mxt2.contractDR(circ>circ_min)
75 img5_neg = mxt2.getImage()
76 adshow(img5_neg,"Only nodes with circularity > %f" %circ_min)
77 
78 print "Circularity filtering"
79 print "Min-tree"
80 print "Number of nodes: %d" %mxt2.node_array.shape[1]
81 print "Number of leaves: %d" %((mxt2.node_array[1,:] == 0).sum())
82 print
83 print "Time: %f s" %(time.time()-t)
Min-tree
Number of nodes: 37003
Number of leaves: 16747

Area filtering
Min-tree
Number of nodes: 509
Number of leaves: 132

Bounding_box filtering
Min-tree
Number of nodes: 422
Number of leaves: 121

Aspect ratio filtering
Min-tree
Number of nodes: 220
Number of leaves: 60

Circularity filtering
Min-tree
Number of nodes: 68
Number of leaves: 2

Time: 0.210079 s

Negative of the image

Only nodes with 100 < area < 400

Only nodes with dx < 35 and dy < 35

Only nodes with AR > 0.750000

Only nodes with circularity > 1.000000

Min-tree + max-tree results

1 adshow(ia870.iagshow(img,img7>0,img5_neg>0),"Composition of the max-tree and min-tree results")#

Composition of the max-tree and min-tree results