convolve_image.m

%--------------------------------------------------------------------------
% convolve_image.m
% main script for evaluating SWF for RPE cobblestone morphology
% (c) 2015-2016 Andrew R. Cohen and Rohini Joshi
%
%   license - for review evaluation only
%       THIS IS UNRELEASED CODE - VIEWING AND DOWNLOAD ONLY BY MANUSCRIPT REVIEWERS
%   the software will be released under the GPL concurrent with publication
%
%       This is the source code for the paper:

%  R. Joshi, M. Winter, J. S. Saini, T. A. Blenkinsop, J. H. Stern, S. Temple and 
%       A. R. Cohen, "Automated measurement of cobblestone morphology for 
%       characterizing Stem cell-Derived Retinal Pigment Epithelial cell cultures".

% EXCEPTION: as described in the paper, the source code for the steerable 
%   wavelet filters must be downloaded separately from http://bigwww.epfl.ch/demo/circular-wavelets. 
%   All use of that component is governed by the license agreement described there.
% 
% Andrew Cohen
% acohen@coe.drexel.edu
% http://bioimage.coe.drexel.ed
%
%--------------------------------------------------------------------------

%% Initialization

clc;
% close all;
clear;
addpath U
addpath Wavelet
addpath Wavelet/Radial

tic;
fname = '..\Images\Temple\RPE\ZOI\Category 5\Category 5 ZO-1 32X C.jpg';
img = rgb2gray(imread(fname));

htype = 'Simoncelli'; % isotropic wavelet type (string), one of 'Simoncelli', 'Meyer', 'Papadakis' or 'Shannon'

Utype = 'Prolate';  %steerable construction type (string). Available constructions are 
%               'Gradient', 'Monogenic', 'Hessian', 'CircularHarmonic',
%               'Prolate', 'ProlateUniSided', 'Slepian', 'RieszEven',
%               'RieszOdd', 'RieszComplete', 'SimoncelliEven', 'SimoncelliOdd',
%               'SimoncelliComplete', 'UserDefined'.

Uparams.harmonics = 0:3:27; % -4:2:4; % For 'RieszEven', 'RieszOdd', 'RieszComplete', 'SimoncelliEven', 'SimoncelliOdd', 
                            %'SimoncelliComplete' - wavelet order
                            
Uparams.B = .1;     % For 'Slepian' - bandwidth 
num_scales = 4;     % number of scales
steerangle = pi/3;  % angle by which to steer the wavelet (optional, defaults to 0)

Mequiangular = 9;   % number of equiangular channels

% init wavelet transform
[hafun,hdfun,hfun,U,harmonics] = Init(htype,Utype,Uparams);

%% Preprocessing

% img = img(350:849, 500:999);
img_raw = img;
% junct = img(230:279, 450:499);
% img = img(1:399, 100:499);
% figure;imagesc(img_raw);colormap gray;
% img = mat2gray(stdfilt(imadjust(img), getnhood(strel('disk',3))));


% Michel contrast enhancement
h=fspecial('gaussian', 7,20);
smooth = imfilter(img,h);
sub=img-smooth;
out=medfilt2(sub, [3 3]);
figure, imagesc(out);colormap gray
img=out;

%% Convolving

[m n] = size(img);
% response = [];
response = zeros([50 50 13]);
loc = rpe_morph(img_raw);
for t=1:13
    psi = psiWav(hfun, U(9,:), harmonics, 0.7, 1, 25, t, 0, 0);
    abs_psi{t} = abs(psi);
end

for i=1:m-50
    for j=1:n-50
        if loc(i+26,j+26)
            for t=1:13
                temp = img(i:i+50-1,j:j+50-1);
                response(i,j,t) = sum(sum(double(temp).*abs_psi{t}));
            end
        else
            response(i,j,:) = 0;
            continue;
        end
    end

end

toc

%%
r = max(response, [], 3);
r = mat2gray(r);
for i=1:15:size(r,1)-15
    for j=1:15:size(r,2)-15
        sub = r(i:i+14,j:j+14);
        r_max(i:i+14,j:j+14) = zeros(15);
        [ii jj] = find(sub == max(sub(:)));
        r_max(i+ii,j+jj) = max(sub(:));
    end
end
rr=zeros(size(r));
rr(1:size(r_max,1),1:size(r_max,2)) = r_max;
r=rr;
[peaks locn] = sort(r(:),'descend');
[I,J] = ind2sub(size(r),locn);
[I J] = find(r>0.3*max(r(:)));
% [I J] = find(r);
r_plot = zeros(size(r));
r_plot(I,J) = r(I,J);
% img_marked(:,:,1:3) = cat(3,img_raw,img_raw,img_raw);
% img_marked(26:size(img_raw,1)-25,:,1) = r_plot;
% img_marked(26:size(img_raw,1)-25,26:size(img_raw,2)-25,1) = r_plot;
figure;
imagesc(img_raw);colormap gray;
axis image
set(gca, 'xtick', [], 'ytick', []);
hold on;
% plot(J(1:end)+26, I(1:end)+26, '+b','markersize',8,'linewidth',3);
% plot([26,26,m-25,m-25],[26,m-25,26,m-25],'+r');
display(['Number of junctions = ' num2str(length(I))]);
r_255 = floor(r*255);
cmap = jet(256);
% figure;
% Find cmap index
for c =1:length(I)
    c_ind(c,:) = cmap(r_255(I(c),J(c))+1,:);
    plot(J(c)+26, I(c)+26, 'h','markerfacecolor', c_ind(c,:),'markersize',7, 'markeredgecolor', 'none');
    hold on;
end

% colorbar
% r=r_max;
% [peaks locn] = sort(r(:),'descend');
% [I,J] = ind2sub(size(r),locn);
% [I J] = find(r>0.35*max(r(:)));
% % [I J] = find(r);
% figure, imagesc(img_raw);colormap gray;
% axis image
% set(gca, 'xtick', [], 'ytick', []);
% hold on;
% plot(J(1:end)+26, I(1:end)+26, '+b','markersize',8,'linewidth',3);
% % plot([26,26,m-25,m-25],[26,m-25,26,m-25],'+r');
% display(['Number of junctions = ' num2str(length(I))]);