%% THIS IS THE MOST RECENT VERSION 04/05/2018
clear all
close all
clc
warning('OFF');
%% Set File Select Location
selected_file = ('');
%% Set Parameters
% The following are the default parameters.  The user will adjust if
% desired in the input box (Line 19)
GthreshDef = 4.5;  %Default hreshold to isolate outline and centerline.  Typically between 5 and 8.  Higher threshold leads to more points, higher computation time.
scaleDef = 2.912; %Default scale in pixels/nm (200nm scale bar, VetMed scope)
ROIminDef = 30; %Identify Region of Interest (30-105nm is standard for perinexus, 10-90 standard for gap junction width)
ROImaxDef = 105;
agl = 0; %Image rotate
Parameters_prompt = {'Gthresh:','Scale (pixels/unit):','ROImin:','ROImax:','Start Point (0 Auto, 1 Manual):'};
dlg_title = 'Set Parameters';
num_lines = 1;
defaultans = {num2str(GthreshDef),num2str(scaleDef),num2str(ROIminDef),num2str(ROImaxDef),num2str(0)};
Parameters = inputdlg(Parameters_prompt,dlg_title,num_lines,defaultans);
if isempty(Parameters);
    return
end
Gthresh = str2double(Parameters(1));
scale = str2double(Parameters(2));
ROImin = str2double(Parameters(3));
ROImax = str2double(Parameters(4));
manstart = str2double(Parameters(5)); % Set to 1 if you want to choose start point of end line, set to 0 to autofind start point
%% Select File
[filename, pathname, filterindex] = uigetfile(selected_file, 'File Selector');
if pathname == 0
    return
end
addpath(pathname)
img = imread(filename);
img = img(:,:,1);
[pathstr,file_name,extens] = fileparts(filename);

%% Save Data Files
%If you want to save the data in a mat file, copy and past the following
%line into the command window and save it.  Then uncomment the next load
%line (line 38).
%WpData = [cellstr(filename ), cellstr( Ave150 ), cellstr( AveROI ), cellstr( Total Elapsed Time ), cellstr( User Input Time)];
% load WpData  %Load data file for storing Wp and time data

% Set the directory to save the list of Wp by distance from GJ edge
savefile = ['',file_name,'.mat'];

%% Figure Save Locations
savestringDC = ['',file_name]; %Directory for final working image figure:
file_name1 = [file_name,'A'];
savestringOCL = ['',file_name1]; %Directory for original image w/ centerline overlay:
savestringWPMP = ['',file_name]; %Directory for Wp/distance-from-GJ plot:

%% Select Perinexus of Interest
img1 = imrotate(img,agl);
figure(2)
%Crop image to isolate desired perinexus.  Cutt off bottom of perinexus, so thre is no part of the bottom line to dilate upwards.
%DOUBLE-CLICK WITHIN PERINEXUS OR REGION OF INTEREST
[~,~,cropim4,RECT] = imcrop(img1);  % Crop image to isolate desired perinexus.  Cut off "open" edge of perinexus, so thre is no part of the line to dilate back into the perinexal space.
if isempty(cropim4)
    return
end
mousecoor = get(gca,'CurrentPoint');
fillcol = floor(mousecoor(1)-floor(RECT(1)));
fillrow = floor(mousecoor(3)-floor(RECT(2)));
figure(1)
imshow(cropim4)
title([file_name, ' cropped'])
pause(0.1)
close figure 2
biimg4 = cropim4 == 255;

userinputtime = toc;   %Elapsed time of user input - everything after this is computation time.
ZZ = strel('disk',2);
biimg3 = imdilate(biimg4,ZZ);

[Ay,Ax] = find(biimg4);

filledimage4 = imfill(biimg3,[fillrow,fillcol]); %Results in filled perinexus, which acts as the guide to tell the program to stop dilating once the perinexus is filled in.
filledimage4 = imfill(filledimage4,'holes');
cropim = cropim4;
biimg1 = biimg4;
filledimage = imerode(filledimage4,ZZ);

%% Serial Dilations Count Pixels Between Membranes
figure(2)
imshow(biimg1)
finalimg = (imcomplement(biimg1)); %Initialize working image - inverse of the binary image
figure(3)
imshow(finalimg)
biimg = biimg1;
Z = biimg.*filledimage;  %Z is the "checking" image, to see if the perinexus is filled in.
matchcheck = Z == filledimage;

se = strel('disk',1); %Dilation kernel - cross-shaped, 1 pixel in size

%Perform 1-pixel dilations, adding the inverse of each dilation to the
%working image.  The process continues until the perinexus is filled in.
while mean(mean(matchcheck)) ~= 1

    dilated = imdilate(biimg,se);
    invertimg = imcomplement(dilated);
    biimg = dilated;
    finalimg = finalimg + invertimg;
    Z = biimg.*filledimage;
    matchcheck = Z == filledimage;
end
figure(4)
imagesc(finalimg)

%% Isolate Centerline With Spatial Derivative
[Gmag, Gdir] = imgradient(finalimg);
figure(5)
imagesc(Gmag)  %Magnitude of the spatial derivative isolates areas of discontinuity
colormap jet
img2 = Gmag <= Gthresh;
img3 = img2.*Z;
SE = strel('disk',1); %Secondary dilation kernel for outline removal.
biimg2 = imdilate(biimg1,SE);  %Dilate the original outline
Centerline1 = img3-biimg2;  %Subtract out original outline
Centerline = Centerline1 ==1;

%Another dilation is used on the centerline to fill in gaps for pathfinding
%algorithm.  SE2 size (NSE2) can be adjusted - higher values
%result in better fill, but can potentially create issues with pathfinding.
NSE2 = 3; % Secondary dilation radius
SE2 = strel('disk',NSE2); % Create secondary dilation kernel
NSE3 = NSE2; % Secondary erosion radius - initiated at one pixel less than dilation radius
SE3 = strel('disk',NSE3); % Create secondary erosion kernel
CLDilate = imdilate(Centerline,SE2); %Dilate to fill holes in centerline
CC1 = bwconncomp(CLDilate); %bwconncomp counts number of objects in BW image
n = 1;
while CC1.NumObjects >1 % Iteratively increase dilation radius until there are no more holes in the centerline
    NSE2 = NSE2 + 1;
    SE2 = strel('disk',NSE2);
    CLDilate = imdilate(Centerline,SE2);
    CC1 = bwconncomp(CLDilate);
end
CLIso = imerode(CLDilate,SE3); % Perform secondary erosion
CC2 = bwconncomp(CLIso); % Count number of objects in isolated centerline
while CC2.NumObjects >1 || CC2.NumObjects == 0 % Iteratively increase erosion radius to give thinnest line without re-introducing discontinuities
    NSE3 = NSE2-n;
    SE3 = strel('disk',NSE3);
    CLIso = imerode(CLDilate,SE3);
    CC2 = bwconncomp(CLIso);
    n = n+1;
end
CLBorderClearInv = imclearborder(CLIso);  %Identifies objects that don't touch an edge
CLBorderClear = CLIso - CLBorderClearInv; %Removes all objects that don't touch an edge
CLFinal1 = bwareaopen(CLBorderClear,50);  %Removes all connected objects with <50 pixels
CLThin = bwmorph(CLFinal1,'skel',Inf);  %Isolates one-pixel centerline
CLFinal2 = bwareaopen(CLThin,50);  %Final clear of non-centerline objects
SE4 = strel('disk',20); % Final dilation kernel ensures final centerline touches an edge
SE5 = strel('disk',20); % Final erosoin kernel reduces final centerline without re-introducing discontinuities
CLFinal3 = imdilate(CLFinal2,SE4); % Perform final dilation/erosion
CLFinal = imerode(CLFinal3,SE5);
CLFinal = CLFinal.*filledimage;
%% Identify Start Point
medfilter = medfilt2(CLFinal,[4 4]);
endpoints = medfilter-CLFinal;
[SR,SC] = find(endpoints);
for i = 1:length(SR)
CLFinal(SR(i),SC(i)) = 1;
end
figure(10)
imshow(CLFinal)
pause(0.1)
busymsg = msgbox('Identifying Centerline. Please Wait.','Busy');
%If centerline is not continuous at this point, Matlab will throw an error,
%typically "undefined function or variable 'xy'".
%% Pathfinding Algorithm Map Centerline
% Copyright (c) 2012, Wasit Limprasert 
% All rights reserved.
finalimgiso = (finalimg+1).*CLFinal;  %Combine final working image and dilated centerline.
map = finalimgiso;
map=double(map(:,:,1))/255; % normalizing to 1
g=Graph(map); %Generating a graph from a mobility map.
[X,Y] = find(finalimgiso);
A = [X,Y];

%% Set end point for pathfinding.
if max(finalimgiso(1,:))>0 %If CL ends on upper edge
    [~, IND] = max(finalimgiso(1,:));
    xy = [IND,1];
end
if max(finalimgiso(end,:))>0 %If CL ends on bottom edge
    [~, IND] = max(finalimgiso(end,:));
    xy = [IND,length(finalimgiso(:,1))];
end
if max(finalimgiso(:,1))>0  %If CL ends on left edge
    [~, IND] = max(finalimgiso(:,1));
    xy = [1,IND];
end
if max(finalimgiso(:,end))>0 %If CL ends on right edge
    [~, IND] = max(finalimgiso(:,end));
    xy = [length(finalimgiso(1,:)), IND];
end

x=round(xy(2));
y=round(xy(1));

%% Identify Starting Node
for i = 1:length(g.node_list)
    posy(i) = g.node_list(i).posx;
    posx(i) = g.node_list(i).posy;
end
posxy = [posx;posy];

for h = 1:length(SC)
    startpointdist(h) = sqrt((SR(h)-x).^2+(SC(h)-y).^2);
end
startindx = find(startpointdist == max(startpointdist));
if manstart == 1
figure(12)
imshow(cropim)
hold on
plot(Y,X,'b','linewidth',2)
    manstartindx = ginput(1);
    close(12);
for s = 1:length(posx)
startdist(s) = sqrt((posx(s)-manstartindx(1))^2+(posy(s)-manstartindx(2))^2);
end
startpoint = find(startdist==min(startdist));
else
B = find(posx==SC(startindx));
C = find(posy==SR(startindx));

startpoint = intersect(B,C);
end

%% Display found path
initv = startpoint;
shortpathv = initv;
finalv = sub2ind(size(finalimgiso),x,y);
[xmax, ymax]=size(map); 
ShowPath(g,map,shortpathv);

finalv=g.node_index_map(x,y);
[shortpathv, g]=GetPath(g,initv,finalv);
shortpathxy=ShowPath(g,map,shortpathv);   
Pathxy = [shortpathxy(2,:);shortpathxy(1,:)];  %Index path

figure(6)
imagesc(finalimg)
colormap jet
hold on
plot(Pathxy(1,:),Pathxy(2,:),'Linewidth',4, 'Color', 'k')
hold off
title(['Dilation Count & CL - ',file_name])

%% Combine Centerline and Final Working Image
CL_map = zeros(size(finalimg));  %initialize space for isolated centerline
finalimg = finalimg./scale*2; %Divide by scale.  finalimg counts dilations to the centerline, must multiply by 2 to account for Cl to opposite membrane.
for i = 1:length(Pathxy)
CL_map(Pathxy(2,i),Pathxy(1,i))=1;  %create binary image locating centerline
end
CL_Dist = CL_map.*finalimg;    %Centerline with dilation count at each point

figure(7) %**ORIGINAL IMAGE = FIGURE 7**
imshow(cropim)
title(['TEM Image + Centerline: ',file_name])
hold on
plot(Pathxy(1,:),Pathxy(2,:),'r', 'linewidth',3)
hold off
figure(8)
imagesc(CL_Dist)
colormap(jet)
title(['W_p Measurements - ',file_name])
coloaxis = colorbar();
ylabel(coloaxis,'nm')
caxis([0 round(max(max(CL_Dist)))]);
%% Calculate average Wp
avg_Wp = mean(CL_Dist(CL_Dist~=0));
sumdist = zeros(1,length(Pathxy));
dist = zeros(1,length(Pathxy));
for j = 2:length(Pathxy);
    dist(j) = sqrt((Pathxy(1,j)-Pathxy(1,j-1))^2+(Pathxy(2,j)-Pathxy(2,j-1))^2);
    dist(j) = dist(j)./scale;
    sumdist(j) = sum(dist);
end
I = zeros(1,length(Pathxy));
for k = 1:length(Pathxy)
I(k) = CL_Dist(Pathxy(2,k),Pathxy(1,k));
end

distthresh = 8;  %Initial distance threshold - if first measurement is greater than this, flip the matrix.
ylimit = max(I) + 3;  %Max y value for plot y-axis

%% Combine Distance from GJ values and corresponding Wp measurements
zerodist = find(sumdist==0);
if I(zerodist) > distthresh
    distflipped = [fliplr(sumdist); fliplr(I)];
    figure(9)
    subplot(2,1,1)
    whitebg(9,[0.1,0.1,0.1])
    plot(distflipped(1,:),distflipped(2,:),'g','linewidth',2)
    WpList = [distflipped];
    ylim([0,ylimit])
    ROI150 = (distflipped(:,distflipped(1,:)<150));
    Ave150 = mean(ROI150(2,:));
    maxROI = distflipped(:,distflipped(1,:)<ROImax);
    ROI = maxROI(:,maxROI(1,:)>ROImin);
    ylimitROI = max(ROI(2,:)) + 3;
else
    combined = [sumdist;I];
    figure(9)
    subplot(2,1,1)
    whitebg(9,[0.1,0.1,0.1])
    plot(sumdist, I,'g','linewidth',2)
    WpList = [sumdist; I];
    ylim([0,ylimit])
    maxROI = combined(:,combined(1,:)<ROImax);
    ROI = maxROI(:,maxROI(1,:)>ROImin);
    ylimitROI = max(ROI(2,:)) +3;
    ROI150 = (combined(:,combined(1,:)<150));
    Ave150 = mean(ROI150(2,:));  
end
aveROI = mean(ROI(2,:));
%% Rotation Correction
% The orientation of the perinexus can influence how the dilation count is
% calculated, due to the use of a plus-shaped kernel.  This correction is
% calculates the approximate angle of the perinexus (end-to-end in a
% straight line) and performs a trigonomic calculation to quantify the
% perinexus as if it had been imaged horizontally, resulting in the minimal
% Wp value (shortest distance between two points is a straight line).

a = [Pathxy(1,1),Pathxy(1,end)];
b = [Pathxy(2,1),Pathxy(2,end)];
aa = Pathxy(1,:);
bb = Pathxy(2,:);
deltax = a(2)-a(1);
deltay = b(1)-b(2);
theta = atand(deltay/-deltax);
% figure(12)
% imshow(imrotate(cropim,theta))
% hold on
% plot(aa,bb,'b','linewidth',3)
% plot(a,b,'k','linewidth',3)
% plot([a(1),a(2)],[b(1),b(1)],'r','linewidth',3)
% plot([a(2),a(2)],[b(1),b(2)],'g','linewidth',3)
% hold off
% aveROICorrected = aveROI-(scale*abs(sind(2*theta)));
% Ave150Corrected = Ave150-(scale*abs(sind(2*theta)));
if abs(theta)<45
    aveROICorrected = aveROI*abs(cosd(theta));
    Ave150Corrected = Ave150*abs(cosd(theta));
else
    aveROICorrected = aveROI*abs(sind(theta));
    Ave150Corrected = Ave150*abs(sind(theta)); 
end

%% Display and Save Results

if round(150*scale) > length(aa)
    x150 = aa(1:end);
    y150 = bb(1:end);
    fprintf('WARNING: Selected perinexus may be shorter than 150nm\n')
else
x150 = aa(1:round(150*scale));
y150 = bb(1:round(150*scale));
end

if round(ROImax*scale) > length(aa)
    if round(ROImin*scale) > length(aa)
        fprintf('WARNING: You have a tiny perinexus.\n')
        ROIx = 1;
        ROIy = 1;
    else
    ROIx = aa(round(ROImin*scale):end);
    ROIy = bb(round(ROImin*scale):end);
    fprintf('WARNING: Selected region may be smaller than defined ROI\n')
    end
else
ROIx = (aa(round(ROImin*scale):round(ROImax*scale)));
ROIy = (bb(round(ROImin*scale):round(ROImax*scale)));
end
figure(7)
hold on
plot(x150,y150,'g','linewidth',3)
plot(ROIx,ROIy,'b','linewidth',3)
hold off

figure(9)
title(['W_p & Distance from Start: ',file_name])
xlabel('Distance from Start, nm')
ylabel('W_p, nm')
xlim([0,150])

subplot(2,1,2,'color',[0.1,0.1,0.1])
plot(ROI(1,:),ROI(2,:),'b','linewidth',2)
xlim([ROImin,ROImax])
ylim([0,ylimitROI])
xlabel('Distance from Start, nm')
ylabel('W_p, nm')

fprintf('File: %s\n',file_name)
fprintf('Average W_p (0-150nm): %2.2f\n', Ave150Corrected);
printformat = 'Average W_p(%2.0f-%3.0fnm): %2.2f\n';
fprintf(printformat, ROImin,ROImax,aveROICorrected)
eltime = toc;  %Total time elapsed - user input and computational.

% fprintf('User input time: %2.2f\n', userinputtime)
% fprintf('Total Elapsed Time: %2.2f\n',eltime)

figs2keep = [7,9];

% Uncomment the following to 
% include ALL windows, including those with hidden handles (e.g. GUIs)
% all_figs = findall(0, 'type', 'figure');

all_figs = findobj(0, 'type', 'figure');
delete(setdiff(all_figs, figs2keep));
close (busymsg)
figure(7)
% save(savefile,'WpList')
% 
% 
% % If you are using the WpData file from the top of the code, uncomment the
% %following lines to save relevant data and figures.
% 
% %Save data - [File name, 0-150nm Wp, ROIWp, total time, user input time]
% % WpData = [WpData;cellstr(file_name), num2cell(Ave150), num2cell(aveROI), num2cell(eltime), num2cell(userinputtime)];
% % save ('WpData','WpData');
% % %Save final working image figure
% % savestringDC = ['',file_name]; %Insert save location
% saveas(figure(6),savestringDC)
% % %Save original image w/ centerline overlay
% % savestringOCL = ['',file_name]; %Insert save location
% saveas(figure(7),savestringOCL)
% % %Save Wp/distance from GJ plot
% % savestringWPMP = ['',file_name]; %Insert save location
% saveas(figure(9),savestringWPMP)