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

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: [email protected]
#
#
# 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()

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