"""
BuildILMGraph: Module for building graph of ILM layer
-----------------------------------------------------
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
w_min = 1e-5
graph_list = {}
endnodes = {}
[docs]def getILMGraph(_slice, scaling, shortest_path , y_offset):
"""
Construct the undirected, weighted ILM graph from the weights.
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
"""
start__slice = 0
#initial slice
if shortest_path is None:
_slice_center_y = _slice.shape[0]//2
#restrict area of search
start__slice = int((_slice_center_y-130)*scaling)
end__slice = int((_slice_center_y-20)*scaling)
_slice_buffer = _slice[start__slice:end__slice,:]
#add border to left and right
_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_1D = _slice_buffer.flatten()
#slice has predecessor
else:
_slice_buffer = _slice.copy()
#add border to left and right
_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)
#get predecessors shortest path
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 predecessors shortest path
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)
#set invalid pixels to nan
buffer_result = np.where(buffer_result == 1.0, _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])
# set borders to -2/+2
buffer_result[:,0] = -2
buffer_result[:,-1] = 2
_slice_1D = buffer_result.flatten()
# graph construction
g = Graph(directed=True)
g.add_vertices(_slice_1D.size)
sx = _slice_buffer.shape[1]
edge_list = []
weight_list = []
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(calculateWeight(_slice_1D[v_i],_slice_1D[v_j]))
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(calculateWeight(_slice_1D[v_i],_slice_1D[v_j]))
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(calculateWeight(_slice_1D[v_i],_slice_1D[v_j]))
#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(calculateWeight(_slice_1D[v_i],_slice_1D[v_j]))
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(calculateWeight(_slice_1D[v_i],_slice_1D[v_j]))
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 calculateWeight(value_v_i, value_v_j):
"""
Exponential weight calculation
When a pixel's intensity is zero, 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
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(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)