-
Notifications
You must be signed in to change notification settings - Fork 0
/
RetinalFeatureExtractor.py
164 lines (128 loc) · 5.43 KB
/
RetinalFeatureExtractor.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
import numpy as np
from os import path
import cv2
from matplotlib import pyplot as plt
from PIL import Image
from skimage.morphology import disk, opening, erosion
from scipy.ndimage import median_filter, binary_opening, binary_erosion, label, gaussian_filter
class RetinalFeatureExtractor:
"""
This class helps extract features from a retinal image for the purpose of later
using these feature points for 2D registration.
"""
@staticmethod
def segment_blood_vessel(image: np.ndarray) -> np.ndarray:
"""
Receives an image of a retina as a ndarray and tries to segment the blood vessels from within the retina image.
Parameters
----------
image
Returns
segmentation image.
-------
"""
# worked according to this method: https://iopscience.iop.org/article/10.1088/1742-6596/1376/1/012023/pdf
# crop bottom label
image = image[:image.shape[0] - 100, :, :]
# Get green channel
segmented_image = image[:, :, 1]
# take out cancer
# segmented_image[450:1000, 500:1200] = 255
# get complement of green channel
segmented_image = np.max(segmented_image) - segmented_image
plt.imshow(segmented_image, cmap='gray')
plt.show()
# use contrast limited adaptive histogram equalization
clahe = cv2.createCLAHE(clipLimit=20, tileGridSize=(8, 8))
segmented_image = clahe.apply(segmented_image)
# # Otsu thresholding
segmented_image = median_filter(segmented_image, 7)
ret3, segmented_image = cv2.threshold(segmented_image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
# Morphological operations
segmented_image = binary_opening(segmented_image, disk(3), iterations=2)
segmented_image = binary_opening(segmented_image, disk(1), iterations=5)
segmented_image = binary_erosion(segmented_image, disk(3), iterations=2)
segmented_image = median_filter(segmented_image, footprint=disk(3))
segmented_image, components = label(segmented_image)
segmented_image = segmented_image == np.argmax(np.bincount(segmented_image.flat)[1:]) + 1
plt.imshow(segmented_image, cmap='gray')
plt.show()
return segmented_image
@staticmethod
def segment_blood_vessel_2(image: np.ndarray) -> np.ndarray:
"""
Receives an image of a retina as a ndarray and tries to segment the blood vessels from within the retina image.
Parameters
----------
image
Returns
segmentation image.
-------
"""
# worked according to this method: https://iopscience.iop.org/article/10.1088/1742-6596/1376/1/012023/pdf
# crop bottom label
image = image[:image.shape[0] - 100, :, :]
# Get green channel
segmented_image = image[:, :, 1]
# get complement of green channel
segmented_image = np.max(segmented_image) - segmented_image
plt.imshow(segmented_image, cmap='gray')
plt.show()
# use contrast limited adaptive histogram equalization
clahe = cv2.createCLAHE(clipLimit=100, tileGridSize=(8, 8))
segmented_image = clahe.apply(segmented_image)
plt.imshow(segmented_image, cmap='gray')
plt.show()
segmented_image = median_filter(segmented_image, 2)
segmented_image = median_filter(segmented_image, 3)
segmented_image = median_filter(segmented_image, 4)
segmented_image = median_filter(segmented_image, 5)
plt.imshow(segmented_image, cmap='gray')
plt.show()
# take out cancer
segmented_image[450:1000, 500:1200] = 0
segmented_image[opening(segmented_image, disk(8)) < 50] = 0
segmented_image = erosion(segmented_image, disk(3))
segmented_image = erosion(segmented_image, disk(3))
segmented_image = erosion(segmented_image, disk(3))
segmented_image = erosion(segmented_image, disk(3))
plt.imshow(segmented_image, cmap='gray')
plt.show()
segmented_image = median_filter(segmented_image, 2)
segmented_image = median_filter(segmented_image, 3)
segmented_image = median_filter(segmented_image, 4)
segmented_image = median_filter(segmented_image, 5)
segmented_image = gaussian_filter(segmented_image, sigma=3)
plt.imshow(segmented_image, cmap='gray')
plt.show()
return segmented_image
@staticmethod
def find_retina_features(image: np.ndarray) -> list:
"""
Uses cv2 to extract SIFT keypoints
Parameters
----------
image
Returns
-------
"""
sift = cv2.SIFT_create(100)
kp = sift.detect(image, None)
sift_image = cv2.drawKeypoints(image, kp, image)
# show the image
plt.imshow(sift_image)
plt.show()
points = np.array([point.pt for point in kp])
# n = 2000
# index = np.random.choice(points.shape[0], n, replace=False)
# return points[index]
return points
if __name__ == '__main__':
dataset_path = "/home/edan/Desktop/HighRad/Exercises/data/Targil2_data_2018-20230315T115832Z-001/Targil2_data_2018"
# bl = "BL01.tif"
bl = "BL04.bmp"
bl_path = path.join(dataset_path, bl)
image = Image.open(bl_path)
im_arr = np.array(image)
RetinalFeatureExtractor.segment_blood_vessel(im_arr)
# RetinalFeatureExtractor.find_retina_features(im_arr)