Skip to content

Instantly share code, notes, and snippets.

@YaronBlinder

YaronBlinder/calculateFVD.m

Created August 10, 2016 17:27
Show Gist options
  • Save YaronBlinder/4805958e7f1c86c1edc1b0f1ff5b9005 to your computer and use it in GitHub Desktop.
Save YaronBlinder/4805958e7f1c86c1edc1b0f1ff5b9005 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
calculateFVD.m
function FVD = calculateFVD(filename, black, intensity_threshold, size_threshold)
% Prompt user for filename.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%prompt={'File name (including extension'}; % Bring up prompt to ask for file name for processing
%def={'example.tif'}; % Predefined answer
%dlgTitle='Input Parameters (all files must be in current directory)'; % Title of prompt
%lines=1; % Lines in prompt
%Resize='on'; % Allow resizing
% answer=inputdlg(prompt,dlgTitle,lines,def,Resize); % Record response
% filename = answer{1}; % Take file name and call save as variable 'filename'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert image to double precision and convert to grayscale if RGB.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I = double(imread(filename)); % Read in file for processing ans convert to double
if size(I,3) == 3 % Convert to grayscale if image is RGB
I = 0.2990.*I(:,:,1) + 0.5870.*I(:,:,2) + 0.1140.*I(:,:,3); % Convert to grayscale by scaling R,G, and B components.
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ask user if vessels are white on black background or vice versa (done for
% proper thresholding).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure
% imagesc(I);axis image;colormap gray % Display image
%S = {'Black','White'}; % Bring up prompt asking user whether vessels are white or black aginst the opposing color background (i.e. white vessels on black background)
% [Selection,ok] = listdlg('ListString',S,'SelectionMode','single','PromptString','What color are the vessels?','ListSize',[160 80]);
% close all
Selection = black;
if Selection == 1; % If vessels are black against white background, invert image
I = I.*(-1)+(16.^2-1); % Invert image
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Enhance edges of image for ease in selecting ROI and better thresholding
% results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I_FT = fft2(I); % Compute fft of image
I_FT_SHIFT = fftshift(I_FT); % fftshifimage (DC now at u,v = 0,0)
num_rows = size(I_FT_SHIFT,1); % Find number of rows of image in freq space
num_cols = size(I_FT_SHIFT,2); % Find number of columns of image in freq space
u = -num_rows./2+1:num_rows./2; % Define bounds of u and v
v = -num_cols./2+1:num_cols./2;
% Use a low pass filter on image, call the result 'Lowpass'
sigma = 7; % Std dev parameter for gaussian (chosen empirically)
for i = 1:num_rows;
for j = 1:num_cols;
Lowpass(i,j)= exp(-2.*(pi.*sigma).^2.*((u(i)./num_rows).^2+(v(j)./num_cols).^2)); % Create low pass filter
end
end
I_FT_SHIFT_LP_FILTERED = I_FT_SHIFT.*Lowpass; % Applying low pass filter
I_FT_LP_FILTERED = ifftshift(I_FT_SHIFT_LP_FILTERED); % Inverse fftshift
I_LP_FILTERED = ifft2(I_FT_LP_FILTERED); % Inverse FT
I_edges = I - I_LP_FILTERED; % Find image with accentuated edges by subtracting LP filtered image from original
I2 = 10.*real(I_edges)+I; % Create new image composed of oringal image with addition of accentuated edges (multiplicative factor chosen empirically)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Allow user to select ROI.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display('Choose ROI on edge-enhanced image.');
% display('Left click to select points.'); % Display instructions
% display('Right click to complete selection.');
% display('Double-click ROI when finished.');
% roi = roipoly(I2./max(I2(:))); % Allow user to select ROI
% area_pixel_space = sum(roi(:)); % Find area of ROI (for later calculation of FVD)
% I2 = I2.*roi; % Use of the portion of I2 within the selected ROI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Allow user to threshold image using intensity criterion. Threshold value can be
% adjusted.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%close all % Close all open images
reply = 1; % Set reply variable which will be used for looping if intensity threshold value is to be changed
while reply == 1; % While user wants to change intensity threshold (or during first run)...
%prompt={'Intensity threshold value'}; % Bring ip prompt to ask for intensity threshold level
%def={'230'};
%dlgTitle='Input Parameter';
%lines=1;
%Resize='on';
% answer1=inputdlg(prompt,dlgTitle,lines,def,Resize);
answer1 = {intensity_threshold};
I_thresh_level = str2double(answer1{1}); % Input is double precision number
fin = I2<I_thresh_level; % Threshold I2 by taking only values less than user defined threshold
fin = not(fin); % Invert image
% figure
% imagesc(fin);colormap gray;axis off;axis image;title(['Image thresholded using intensity >' num2str(I_thresh_level)]) % Display thresholded image for user to accept/decline result
% reply = input('Accept this threshold value? Yes: enter. No: 1,enter '); % Ask user if resulting image is accecptable at command prompt. Hit enter for yes and 1, enter for no
reply = '';
if isempty(reply) % If user hits enter, no longer cycle through intensity threshold level selection
reply = 'done';
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Allow user to threshold image using size criterion. Threshold value can be
% adjusted.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
reply1= 1; % Set reply variable which will be used for looping if intensity threshold value is to be changed
while reply1 == 1; % While user wants to change size threshold (or during first run)...
prompt={'Size threshold value'}; % Bring ip prompt to ask for size threshold level
def={'100'};
dlgTitle='Input Parameter';
lines=1;
Resize='on';
% answer2=inputdlg(prompt,dlgTitle,lines,def,Resize);
answer2 = {size_threshold}
size_thresh_level = str2double(answer2{1}); % Input is double precision number
fin_label = bwlabel(fin); % Label all features with value 1 using bwlabel
stats = regionprops(fin_label,'all'); % Acquire stats for feature labeled
idx = find([stats.Area] > size_thresh_level); % Select features with area less than selected area threshold
fin2 = ismember(fin_label,idx); % New is image is previous image with features with area less than size threshold removed
% figure
% imagesc(fin2);colormap gray;axis off;axis image;title(['Image thresholded using area >' num2str(size_thresh_level)]) % Display thresholded image for user to accept/decline result
% reply1 = input('Accept this threshold value? Yes: enter. No: 1,enter '); % Ask user if resulting image is accecptable at command prompt. Hit enter for yes and 1, enter for no
reply1 = '';
if isempty(reply) % If user hits enter, no longer cycle through intensity threshold level selection
reply1 = 'done';
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use a median filter on thresholded image to reduce edge artifacts during
% skeletonization and then skeletonize image using Zhang Suen method.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fin2 = medfilt2(fin2,[3 3]); % Apply 3x3 median filter to image before skeletonizing to avoid edge artifacts
skel = Zhang_Suen(fin2); % Skeletonize image using modified Zhang Suen method
% figure
% imagesc(skel);colormap gray;axis image;title('Skeletonzied image') % Display skeletonized image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate FVD and display result.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
total_length = length_calculator(skel); % Find length of all vessels using point by point accumulation (see length_calculator.m)
FVD = total_length./area_pixel_space; % Calculate FVD
% disp(['Total vessel length [pixels] = ' num2str(total_length)]); % Display total vessel length in pixels
% disp(['FVD [pixels^(-1)] = ' num2str(FVD)]); % Display computed FVD in pixels^(-1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% close all
Raw
calculateFVD_ported.py
# Original user input:
# 1. Filename
# 2. Black/White vessels (default = White)
# 3. Intensity threshold value (default = 230)
# 4. Size threshold value (default = 100)
# CLI - calculateFVD [-filename] filename [--black] [-i] i [-s] s
from helpers_calculateFVD import str2double_m, ismember_m, disp_m
from roipoly import roipoly
import numpy as np
from cv2 import imread, threshold, THRESH_BINARY_INV
from scipy.signal import medfilt
from scipy.ndimage import label
from skimage.measure import regionprops
import matplotlib.pyplot as plt
from Zhang_Suen_ported import Zhang_Suen_ported
from length_calculator_ported import length_calculator
import argparse
from rt_insert import mimic
GLOBAL_filename = None
GLOBAL_black = None
GLOBAL_intensity_threshold = None
GLOBAL_size_threshold = None
def inputdlg_m1(*args, **kwargs):
answer = {1 : GLOBAL_filename}
return(answer)
def inputdlg_m2(*args, **kwargs):
answer = {1 : str(GLOBAL_intensity_threshold)}
return(answer)
def inputdlg_m3(*args, **kwargs):
answer = {1 : str(GLOBAL_size_threshold)}
return(answer)
def listdlg_m(*args, **kwargs):
return(GLOBAL_black)
@mimic('calculateFVD')
def calculateFVD_ported(filename, black, intensity_threshold, size_threshold):
####################################################################################################################################################################
# Code written to compute functional vascular density (FVD) from images of vasculature.
# FVD is defined as the length of functional (perfused) vessels divided by
# the area of interest.
# Images require only contrast between vessels and background and can be collected
# using any appropriate methodology. The following m-files are required for proper
# functionality:
#
# 1) calculateFVD.m (this m-file)
# 2) length_calculator.m
# 3) Zhang_Suen.m
# 4) Zhang_Suen_Even.m
# 5) Zhang_Suen_Odd.m
#
# Code written by Sean M. White
# Cardiopulmonary Transport and Tissue Remodeling Lab &
# Microvascular Therapeutics and Imaging Lab
# Biomedical Engineering Department
# University of California, Irvine
# E-mail address: seanmw@uci.edu
#
#
# References:
# 1. T.Y. Zhang and C.Y. Suen, �A Fast Parallel Algorithm for Thinning Digital Patterns,� Image Processing and Computer Vision. 27(3), 236-239 (1984).
# 2. H. Lu and P. Wang, �A Comment on A Fast Parallel Algorithm for Thinning Digital Patterns,� Image Processing and Computer Vision. 29(3), 239-242 (1986).
# 3. Alasdair McAndrew, Introduction to digital image processing. Thomson. 2004
#
#
# Version 1.0
####################################################################################################################################################################
# # Prompt user for filename.
# ####################################################################################################################################################################
# prompt={'File name (including extension'}; # Bring up prompt to ask for file name for processing
# def={'example.tif'}; # Predefined answer
# dlgTitle='Input Parameters (all files must be in current directory)'; # Title of prompt
# lines=1; # Lines in prompt
# Resize='on'; # Allow resizing
# answer=inputdlg(prompt,dlgTitle,lines,def,Resize); # Record response
# filename = answer{1}; # Take file name and call save as variable 'filename'
# ####################################################################################################################################################################
# Prompt user for filename.
####################################################################################################################################################################
# prompt = 'File name (including extension' # Bring up prompt to ask for file name for processing
# _def = 'example.tif' # Predefined answer
# dlgTitle = 'Input Parameters (all files must be in current directory)' # Title of prompt
# lines = 1 # Lines in prompt
# Resize = 'on' # Allow resizing
# answer=inputdlg_m1(prompt,dlgTitle,lines,_def,Resize) # Record response
# filename = answer[1] # Take file name and call save as variable 'filename'
####################################################################################################################################################################
# # Convert image to double precision and convert to grayscale if RGB.
# ####################################################################################################################################################################
# I = double(imread(filename)); # Read in file for processing ans convert to double
# if size(I,3) == 3 # Convert to grayscale if image is RGB
# I = 0.2990.*I(:,:,1) + 0.5870.*I(:,:,2) + 0.1140.*I(:,:,3); # Convert to grayscale by scaling R,G, and B components.
# end
# ####################################################################################################################################################################
# Convert image to double precision and convert to grayscale if RGB.
####################################################################################################################################################################
I = np.double(imread(filename)) # Read in file for processing ans convert to double
if I.shape[2] == 3: # Convert to grayscale if image is RGB
I = 0.2990*I[:,:,2] + 0.5870*I[:,:,1] + 0.1140*I[:,:,0]; # Convert to grayscale by scaling R,G, and B components. (corrected for BGR)
####################################################################################################################################################################
# # Ask user if vessels are white on black background or vice versa (done for
# # proper thresholding).
# ####################################################################################################################################################################
# figure
# imagesc(I);axis image;colormap gray # Display image
# S = {'Black','White'}; # Bring up prompt asking user whether vessels are white or black aginst the opposing color background (i.e. white vessels on black background)
# [Selection,ok] = listdlg('ListString',S,'SelectionMode','single','PromptString','What color are the vessels?','ListSize',[160 80]);
# close all
# if Selection == 1; # If vessels are black against white background, invert image
# I = I.*(-1)+(16.^2-1); # Invert image
# end
# ####################################################################################################################################################################
# Ask user if vessels are white on black background or vice versa (done for
# proper thresholding).
####################################################################################################################################################################
# plt.imshow(I, cmap = 'gray') # Display image
# plt.axis('image')
# plt.show()
# S = ['Black','White']; # Bring up prompt asking user whether vessels are white or black aginst the opposing color background (i.e. white vessels on black background)
# Selection = listdlg_m('ListString',S,'SelectionMode','single','PromptString','What color are the vessels?','ListSize',[160, 80]);
# plt.close()
Selection = black
if Selection == 1: # If vessels are black against white background, invert image
I = I*(-1)+(16 **2-1); # Invert image
####################################################################################################################################################################
# # Enhance edges of image for ease in selecting ROI and better thresholding
# # results.
# ####################################################################################################################################################################
# I_FT = fft2(I); # Compute fft of image
# I_FT_SHIFT = fftshift(I_FT); # fftshifimage (DC now at u,v = 0,0)
# num_rows = size(I_FT_SHIFT,1); # Find number of rows of image in freq space
# num_cols = size(I_FT_SHIFT,2); # Find number of columns of image in freq space
# u = -num_rows./2+1:num_rows./2; # Define bounds of u and v
# v = -num_cols./2+1:num_cols./2;
# # Use a low pass filter on image, call the result 'Lowpass'
# sigma = 7; # Std dev parameter for gaussian (chosen empirically)
# for i = 1:num_rows;
# for j = 1:num_cols;
# Lowpass(i,j)= exp(-2.*(pi.*sigma).^2.*((u(i)./num_rows).^2+(v(j)./num_cols).^2)); # Create low pass filter
# end
# end
# I_FT_SHIFT_LP_FILTERED = I_FT_SHIFT.*Lowpass; # Applying low pass filter
# I_FT_LP_FILTERED = ifftshift(I_FT_SHIFT_LP_FILTERED); # Inverse fftshift
# I_LP_FILTERED = ifft2(I_FT_LP_FILTERED); # Inverse FT
# I_edges = I - I_LP_FILTERED; # Find image with accentuated edges by subtracting LP filtered image from original
# I2 = 10.*real(I_edges)+I; # Create new image composed of oringal image with addition of accentuated edges (multiplicative factor chosen empirically)
# ####################################################################################################################################################################
# Enhance edges of image for ease in selecting ROI and better thresholding
# results.
####################################################################################################################################################################
I_FT = np.fft.fft2(I) # Compute fft of image
I_FT_SHIFT = np.fft.fftshift(I_FT) # fftshifimage (DC now at u,v = 0,0)
num_rows = I_FT_SHIFT.shape[0] # Find number of rows of image in freq space
num_cols = I_FT_SHIFT.shape[1] # Find number of columns of image in freq space
u = np.arange(-num_rows/2+1, num_rows/2+1) # Define bounds of u and v
v = np.arange(-num_cols/2+1, num_cols/2+1)
# Use a low pass filter on image, call the result 'Lowpass'
sigma = 7 # Std dev parameter for gaussian (chosen empirically)
Lowpass = np.zeros([num_rows, num_cols])
for i in range(num_rows):
for j in range(num_cols):
Lowpass[i,j]= np.exp(-2*(np.pi*sigma)**2*((u[i]/num_rows)**2+(v[j]/num_cols)**2)) # Create low pass filter *** check element-wise operations!
I_FT_SHIFT_LP_FILTERED = I_FT_SHIFT*Lowpass # Applying low pass filter
I_FT_LP_FILTERED = np.fft.ifftshift(I_FT_SHIFT_LP_FILTERED) # Inverse fftshift
I_LP_FILTERED = np.fft.ifft2(I_FT_LP_FILTERED) # Inverse FT
I_edges = I - I_LP_FILTERED # Find image with accentuated edges by subtracting LP filtered image from original
I2 = 10*np.real(I_edges)+I # Create new image composed of oringal image with addition of accentuated edges (multiplicative factor chosen empirically)
####################################################################################################################################################################
# # Allow user to select ROI.
# ####################################################################################################################################################################
# display('Choose ROI on edge-enhanced image.');
# display('Left click to select points.'); # Display instructions
# display('Right click to complete selection.');
# display('Double-click ROI when finished.');
# roi = roipoly(I2./max(I2(:))); # Allow user to select ROI
# area_pixel_space = sum(roi(:)); # Find area of ROI (for later calculation of FVD)
# I2 = I2.*roi; # Use of the portion of I2 within the selected ROI
# ####################################################################################################################################################################
# Allow user to select ROI.
####################################################################################################################################################################
# print('Choose ROI on edge-enhanced image.')
# print('Left click to select points.') # Display instructions
# print('Right click to complete selection.')
# print('Double-click ROI when finished.')
# plt.imshow(I2/I2.max(), cmap = 'gray')
# plt.title('Choose ROI on edge-enhanced image.')
# roi = roipoly(roicolor='r') # Allow user to select ROI
# area_pixel_space = roi.getMask(I2).sum() # Find area of ROI (for later calculation of FVD)
# I2 = I2*roi.getMask(I2) # Use of the portion of I2 within the selected ROI
####################################################################################################################################################################
# # Allow user to threshold image using intensity criterion. Threshold value can be
# # adjusted.
# ####################################################################################################################################################################
# close all # Close all open images
# reply = 1; # Set reply variable which will be used for looping if intensity threshold value is to be changed
# while reply == 1; # While user wants to change intensity threshold (or during first run)...
# prompt={'Intensity threshold value'}; # Bring ip prompt to ask for intensity threshold level
# def={'230'};
# dlgTitle='Input Parameter';
# lines=1;
# Resize='on';
# answer1=inputdlg(prompt,dlgTitle,lines,def,Resize);
# I_thresh_level = str2double(answer1{1}); # Input is double precision number
# fin = I2<I_thresh_level; # Threshold I2 by taking only values less than user defined threshold
# fin = not(fin); # Invert image
# figure
# imagesc(fin);colormap gray;axis off;axis image;title(['Image thresholded using intensity >' num2str(I_thresh_level)]) # Display thresholded image for user to accept/decline result
# reply = input('Accept this threshold value? Yes: enter. No: 1,enter '); # Ask user if resulting image is accecptable at command prompt. Hit enter for yes and 1, enter for no
# if isempty(reply) # If user hits enter, no longer cycle through intensity threshold level selection
# reply = 'done';
# end
# end
# ####################################################################################################################################################################
# Allow user to threshold image using intensity criterion. Threshold value can be
# adjusted.
####################################################################################################################################################################
# plt.close() # Close all open images
reply = 1 # Set reply variable which will be used for looping if intensity threshold value is to be changed
while reply == 1: # While user wants to change intensity threshold (or during first run)...
prompt='Intensity threshold value' #Deprecated, moved to command line flag. Bring ip prompt to ask for intensity threshold level
_def='230' #Deprecated, moved to command line flag.
dlgTitle='Input Parameter' #Deprecated, moved to command line flag.
lines=1 #Deprecated, moved to command line flag.
Resize='on' #Deprecated, moved to command line flag.
answer1=inputdlg_m2(prompt,dlgTitle,lines,_def,Resize)
I_thresh_level = str2double_m(answer1[1]) # Input is double precision number
fin = (I2<I_thresh_level) # Threshold I2 by taking only values less than user defined threshold
# fin = threshold(I2, 0, I_thresh_level, THRESH_BINARY_INV) # Threshold I2 by taking only values less than user defined threshold
fin = np.logical_not(fin) # Invert image
# plt.imshow(fin, cmap = 'gray')
# plt.axis('off')
# plt.axis('image')
# plt.title('Image thresholded using intensity > ' + str(I_thresh_level))
# plt.show()
# reply = input('Accept this threshold value? Yes: enter. No: 1,enter ') # Ask user if resulting image is accecptable at command prompt. Hit enter for yes and 1, enter for no
reply = ''
if reply == '': # If user hits enter, no longer cycle through intensity threshold level selection
reply = 'done'
####################################################################################################################################################################
#
# # Allow user to threshold image using size criterion. Threshold value can be
# # adjusted.
# ####################################################################################################################################################################
# reply1= 1; # Set reply variable which will be used for looping if intensity threshold value is to be changed
# while reply1 == 1; # While user wants to change size threshold (or during first run)...
# prompt={'Size threshold value'}; # Bring ip prompt to ask for size threshold level
# def={'100'};
# dlgTitle='Input Parameter';
# lines=1;
# Resize='on';
# answer2=inputdlg(prompt,dlgTitle,lines,def,Resize);
# size_thresh_level = str2double(answer2{1}); # Input is double precision number
# fin_label = bwlabel(fin); # Label all features with value 1 using bwlabel
# stats = regionprops(fin_label,'all'); # Acquire stats for feature labeled
# idx = find([stats.Area] > size_thresh_level); # Select features with area less than selected area threshold
# fin2 = ismember(fin_label,idx); # New is image is previous image with features with area less than size threshold removed
# figure
# imagesc(fin2);colormap gray;axis off;axis image;title(['Image thresholded using area >' num2str(size_thresh_level)]) # Display thresholded image for user to accept/decline result
# reply1 = input('Accept this threshold value? Yes: enter. No: 1,enter '); # Ask user if resulting image is accecptable at command prompt. Hit enter for yes and 1, enter for no
# if isempty(reply) # If user hits enter, no longer cycle through intensity threshold level selection
# reply1 = 'done';
# end
# end
# ####################################################################################################################################################################
#
# Allow user to threshold image using size criterion. Threshold value can be
# adjusted.
####################################################################################################################################################################
reply1 = 1 # Set reply variable which will be used for looping if intensity threshold value is to be changed
while reply1 == 1: # While user wants to change size threshold (or during first run)...
prompt = 'Size threshold value' # Bring ip prompt to ask for size threshold level
_def = '100' #Deprecated, moved to command line flag.
dlgTitle = 'Input Parameter' #Deprecated, moved to command line flag.
lines = 1 #Deprecated, moved to command line flag.
Resize = 'on' #Deprecated, moved to command line flag.
answer2 = inputdlg_m3(prompt,dlgTitle,lines,_def,Resize) #Deprecated, moved to command line flag.
size_thresh_level = str2double_m(answer2[1]) # Input is double precision number
fin_label = label(fin)[0] # Label all features with value 1 using bwlabel
stats = regionprops(fin_label) # Acquire stats for feature labeled
idx = [element.label for element in stats if element.area > size_thresh_level] # Select features with area less than selected area threshold
fin2 = ismember_m(fin_label,idx) # New is image is previous image with features with area less than size threshold removed
# plt.imshow(fin2, cmap = 'gray')
# plt.axis('off')
# plt.axis('image')
# plt.title('Image thresholded using area > ' + str(size_thresh_level))
# plt.show()
# reply1 = input('Accept this threshold value? Yes: enter. No: 1,enter '); # Ask user if resulting image is accecptable at command prompt. Hit enter for yes and 1, enter for no
reply1 = ''
if reply1 == '': # If user hits enter, no longer cycle through intensity threshold level selection
reply1 = 'done'
####################################################################################################################################################################
# # Use a median filter on thresholded image to reduce edge artifacts during
# # skeletonization and then skeletonize image using Zhang Suen method.
# ####################################################################################################################################################################
# fin2 = medfilt2(fin2,[3 3]); # Apply 3x3 median filter to image before skeletonizing to avoid edge artifacts
# skel = Zhang_Suen(fin2); # Skeletonize image using modified Zhang Suen method
# figure
# imagesc(skel);colormap gray;axis image;title('Skeletonzied image') # Display skeletonized image
# ####################################################################################################################################################################
# Use a median filter on thresholded image to reduce edge artifacts during
# skeletonization and then skeletonize image using Zhang Suen method.
####################################################################################################################################################################
fin2 = medfilt(fin2, kernel_size = [3, 3]) # Apply 3x3 median filter to image before skeletonizing to avoid edge artifacts
skel = Zhang_Suen(fin2) # Skeletonize image using modified Zhang Suen method
# plt.imshow(skel, cmap = 'gray') # Display skeletonized image
# plt.axis('image')
# plt.title('Skeletonzied image')
# plt.show()
####################################################################################################################################################################
# # Calculate FVD and display result.
# ####################################################################################################################################################################
# total_length = length_calculator(skel); # Find length of all vessels using point by point accumulation (see length_calculator.m)
# FVD = total_length./area_pixel_space; # Calculate FVD
# disp(['Total vessel length [pixels] = ' num2str(total_length)]); # Display total vessel length in pixels
# disp(['FVD [pixels^(-1)] = ' num2str(FVD)]); # Display computed FVD in pixels^(-1)
# ####################################################################################################################################################################
# close all
# Calculate FVD and display result.
####################################################################################################################################################################
total_length = length_calculator(skel) # Find length of all vessels using point by point accumulation (see length_calculator.m)
FVD = total_length/area_pixel_space # Calculate FVD
# disp_m('Total vessel length [pixels] = ' + str(total_length)) # Display total vessel length in pixels
# disp_m('FVD [pixels^(-1)] = ' + str(FVD)) # Display computed FVD in pixels^(-1)
return(FVD)
####################################################################################################################################################################
# plt.close()
def main():
global GLOBAL_filename
global GLOBAL_black
global GLOBAL_intensity_threshold
global GLOBAL_size_threshold
parser = argparse.ArgumentParser()
parser.add_argument('-filename', default='example.tif')
parser.add_argument('--black', default=False, action='store_true')
parser.add_argument('-i', '--intensity-threshold', default=230, type=int)
parser.add_argument('-s', '--size-threshold', default=100, type=int)
ns = parser.parse_args()
GLOBAL_filename = ns.filename
GLOBAL_black = ns.black
GLOBAL_intensity_threshold = ns.intensity_threshold
GLOBAL_size_threshold = ns.size_threshold
calculateFVD(filename=ns.filename, black=ns.black, intensity_threshold=ns.intensity_threshold, size_threshold=ns.size_threshold)
if __name__ == '__main__':
main()
Raw
test_calculateFVD.py
from rt_insert import MANAGER
from calculateFVD_ported import calculateFVD_ported
import numpy as np
from Zhang_Suen_ported import Zhang_Suen_ported
from helpers_calculateFVD import makelut_m, applylut_m
try:
import matlab.engine as engine
except ImportError:
try:
import oct2py
engine = oct2py.Oct2Py(convert_to_float=False)
except ImportError:
raise ImportError("Found neither 'matlab.engine' nor 'oct2py'")
# logging.basicConfig(level=logging.DEBUG)
# engine = Oct2Py(logger=logging.getLogger())
engine.addpath('C:/Program Files/octave-4.0.3/share/octave/packages/image-2.4.1')
engine.addpath('C:/Program Files/octave-4.0.3/lib/octave/packages/image-2.4.1/i686-w64-mingw32-api-v50+')
def callback(func, func_args, func_dargs, dec_args, dec_dargs):
mfunc_name = dec_args[0]
print(func_args)
mfunc_result = getattr(engine, mfunc_name)(*func_args)
local_result = func(*func_args, **func_dargs)
assert type(mfunc_result) == type(local_result), \
"octave type={} != numpy type={}".format(type(mfunc_result), type(local_result))
if type(mfunc_result) == type(local_result) == type(np.array([])):
assert mfunc_result.dtype == local_result.dtype, \
"octave dtype={} != numpy dtype={}".format(mfunc_result.dtype, local_result.dtype)
assert mfunc_result.shape == local_result.shape, \
"octave shape={} != numpy shape={}".format(mfunc_result.shape, local_result.shape)
print(mfunc_name)
print('octave:')
print(mfunc_result)
print('python:')
print(local_result)
print('octave==python:')
print(mfunc_result==local_result)
truth = (mfunc_result==local_result)
print(np.where(~(truth)))
assert np.allclose(mfunc_result, local_result), mfunc_name
return local_result
MANAGER.callback = callback
TEST_CASES = [('example.tif',0,230,100),('example.tif',1,230,100)]
def test_calculateFVD_ported():
for filename, black, intensity_threshold, size_threshold in TEST_CASES:
yield calculateFVD_ported, filename, black, intensity_threshold, size_threshold
def test_Zhang_Suen_ported():
I = np.array([[1,0,1,0,1],[0,1,0,1,0],[1,0,1,0,1],[0,1,0,1,0],[1,0,1,0,1]])
yield Zhang_Suen_ported, I
def test_makelut_m():
TEST_LUT = ['Zhang_Suen_Odd','Zhang_Suen_Even']
for lut in TEST_LUT:
yield makelut_m, lut, 3
def test_applylut_m():
I = np.array([[1,0,1,0,1],[0,1,0,1,0],[1,0,1,0,1],[0,1,0,1,0],[1,0,1,0,1]])
lut = makelut_m('Zhang_Suen_Odd',3)
yield applylut_m, I, lut
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment