"""
BuildRPEGraph (Refine): Refinement of RPE segmentation (9-neighborhood)
-------------------------------------------------------------------------
PRLEC Framework for OCT Processing and Visualization
"""
# This framework evolved from a collaboration of:
# - Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambdrige, MA, US
# - Pattern Recognition Lab, Friedrich-Alexander-Universitaet Erlangen-Nuernberg, Germany
# - Department of Biomedical Engineering, Peking University, Beijing, China
# - New England Eye Center, Tufts Medical Center, Boston, MA, US
# v1.0: Updated on Mar 20, 2019
# @author: Daniel Stromer - EMAIL:daniel.stromer@fau.de
# Copyright (C) 2018-2019 - Daniel Stromer
# PRLE is developed as an Open Source project under the GNU General Public License (GPL) v3.0.
from igraph import *
import numpy as np
import cv2
# added to weights to avoid division by zero (called system stability)
w_min = 1e-5
[docs]def getRPEGraph(_slice, shortest_path, scaling, y_offset):
"""
Construct the undirected, weighted RPE graph from the calculated.
For the initial slice (center), an extensive calculation is conducted.
Afterwards, the shortest path of the predecessor is dilated to
speed up the process. Invalid edges are set to np.nan
Parameters
----------
_slice: ndarray
weight slice
scaling: scalar
scaling factor for calculation
shortest_path: ndarray
the shortest path from the predecessor
y_offset: scalar
the top row extracted from the predecessor
Returns
---------
g: graph
the resulting segmentation
endpoint: scalar, int
the endpoint of the graph
y_offset:
the top row extracted from the predecessor
"""
#initial slice
if shortest_path is None:
# limit to certain region above/beneath flattened RPE
slice_center_y = _slice.shape[0]//2
start_slice = int((slice_center_y-50)*scaling)
end_slice = int((slice_center_y+20)*scaling)
slice_buffer = _slice[start_slice:end_slice,:]
#add a column at the left and right border for start and end points (adding value +2.0 to the right border
slice_buffer = np.where(slice_buffer == 0, 0.1, slice_buffer)
slice_buffer = cv2.copyMakeBorder(cv2.copyMakeBorder(slice_buffer, top=0,bottom=0, left=0, right=1, borderType= cv2.BORDER_CONSTANT, value = 2.0), top=0,bottom=0, left=1, right=0, borderType= cv2.BORDER_CONSTANT, value = -2.0)
slice_buffer = np.where(slice_buffer == -1, np.nan , slice_buffer)
slice_1D = slice_buffer.flatten()
#has predecessor
else:
slice_buffer = _slice.copy()
#add borders and get predecessors shortest path
slice_buffer = cv2.copyMakeBorder(cv2.copyMakeBorder(slice_buffer, top=0,bottom=0, left=0, right=1, borderType= cv2.BORDER_CONSTANT, value = 2.0), top=0,bottom=0, left=1, right=0, borderType= cv2.BORDER_CONSTANT, value = -2.0)
path_x = (shortest_path % slice_buffer.shape[1]).astype(np.int32)
path_y = (shortest_path / slice_buffer.shape[1]).astype(np.int32) + y_offset
start_x = np.where(path_x == 1)[0][-1]
end_x = np.where(path_x == slice_buffer.shape[1] - 3)[0][-1]
min_y = np.amin(path_y[int(start_x):int(end_x)])
max_y = np.amax(path_y[int(start_x):int(end_x)])
buffer_result = np.zeros((slice_buffer.shape)).astype(np.float32)
buffer_result[path_y, path_x] = 1.0
#dilate shortest path from predecessor
kernel = np.maximum(int(10*scaling),2)
y_offset = np.maximum(int(10*scaling),2)
buffer_result[:,1:-1] = cv2.dilate(buffer_result[:,1:-1], np.ones((kernel,kernel),np.uint8),iterations = 1)
buffer_result = np.where(buffer_result == 1, slice_buffer, np.nan)
buffer_result = np.where(slice_buffer[min_y-y_offset:max_y+y_offset] == 0, 0, buffer_result[min_y-y_offset:max_y+y_offset])
buffer_result[:,0] = -2
buffer_result[:,-1] = 2
buffer_result = np.where(buffer_result == -1, np.nan , buffer_result)
slice_1D = buffer_result.flatten()
# construct graph
g = Graph(directed=True)
g.add_vertices(slice_1D.size)
sx = slice_buffer.shape[1]
edge_list = []
weight_list = []
#construct with 9-neighborhood
for idx in range(slice_1D.size):
if np.isnan(slice_1D[idx]):
continue
v_i = idx
#pixel is not at the right boundary of the image -> has right neighbors
if((v_i+1) % sx != 0 ):
v_j = v_i +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i - sx +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i + sx +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i - 2*sx +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i + 2*sx +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i - 3*sx +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i + 3*sx +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i - 4*sx +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i + 4*sx +1
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
#pixel is located at the left/right boundary -> zero weights
if(slice_1D[v_i] > 1.0 or slice_1D[v_i] < -1.0):
slice_1D[v_i] = 0.0
v_j = v_i - sx
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(w_min)
v_j = v_i + sx
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(w_min)
else:
#not a boundary pixel
v_j = v_i - sx
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
v_j = v_i + sx
if (v_j >= 0 and v_j < slice_1D.shape[0] and np.isnan(slice_1D[v_j]) == False):
edge_list.append((v_i, v_j))
weight_list.append(calculateExpWeight(slice_1D[v_i],slice_1D[v_j],shortest_path))
g.add_edges(edge_list)
g.es['weight'] = weight_list
if shortest_path is None:
return g, slice_1D.shape[0]-1, start_slice
else:
return g, slice_1D.size-1, min_y - y_offset
[docs]def calculateExpWeight(value_v_i, value_v_j,shortest_path):
"""
Exponential weight calculation.
When a pixel's intensity is smaller than 0.4, the weight is set to a constant (exp(6))
Parameters
----------
value_v_i: scalar, float
first pixel's intensity
value_v_j: scalar, float
second's pixel's intensity
shortest_path: ndarray
when initial slice (shortest path = None), different behavior
Returns
----------
weight: scalar, float
the edge weight
"""
if(value_v_j > 1.0 or value_v_j < -1.0):
return np.exp(2.0 - (value_v_i) + w_min)
elif shortest_path is None and (value_v_j < 0.4 or value_v_i< 0.4):
return np.exp(6.0)
elif shortest_path is not None and (value_v_j == 0 or value_v_i == 0):
return np.exp(6.0)
else:
return np.exp(2.0 - (value_v_i + value_v_j) + w_min)