function resp_raw = convolve_image_open(fname, outdir)
%--------------------------------------------------------------------------
% convolve_image.m
% main script for evaluating SWF for RPE cobblestone morphology
% (c) 2015-2016 Andrew R. Cohen, Walter Mankowski and Rohini Joshi
%
%
%       This is the source code for the paper:

%  R. Joshi, W. Mankowski, 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".
%
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     The software is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     in file "gnu gpl v3.txt". If not, see <http://www.gnu.org/licenses/>.
%
% 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 'at' coe.drexel.edu
% http://bioimage.coe.drexel.edu
%
%--------------------------------------------------------------------------

%% Initialization

% These folders from the Steerable Wavelet Filters code need to be in your path:
% addpath U
% addpath Wavelet
% addpath Wavelet/Radial

img = imread(fname);
if size(img, 3) > 1
    img = rgb2gray(img);
end
img = mat2gray(img);

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);
[~,~,hfun,U,harmonics] = Init(htype,Utype,Uparams);

%% Preprocessing

img_raw = img;

FS = 7;
MF = 3;
im = img_raw;
imfilt = imgaussfilt(im,FS);
imhfreq = max((im - imfilt), zeros(size(im)));
imf = medfilt2(imhfreq,[MF MF]);
imf = mat2gray(imf);
imwrite(imf, fullfile(outdir, 'michel.tif'), 'tif');
img = imf;

%% Convolving

[m, n] = size(img);
[loc, im_skel] = rpe_morph_open(mat2gray(img));
imwrite(im_skel, fullfile(outdir, 'skeleton.tif'), 'tif');

raw_skel = repmat(img_raw, 1, 1, 3);
skel_clr = repmat(im_skel, 1, 1, 3);
skel_clr(:,:,1) = zeros(size(im_skel));
skel_clr(:,:,3) = zeros(size(im_skel));
raw_skel = raw_skel + skel_clr;
imwrite(raw_skel, fullfile(outdir, 'orig_skel.tif'), 'tif');

ANGLE_DIVS = 120;
SPAN = 2*pi/3;
angles = SPAN/ANGLE_DIVS:SPAN/ANGLE_DIVS:SPAN;
abs_psi = cell(1,ANGLE_DIVS);
for t=1:ANGLE_DIVS
    steerangle = angles(t);
    psi = psiWav(hfun, U(9,:), harmonics, 1, 1, 25, steerangle);
    abs_psi{t} = abs(psi);
%    abs_psi{t} = psi;
end

winsize = 50;
winsize2 = winsize/2;
response = zeros(m, n);
[r,c] = find(loc);
ind = find(loc);
spmd
    for i=labindex:numlabs:length(r)
        if r(i) - winsize2 <= 0 || r(i) + winsize2 > m || c(i) - winsize2 <= 0 || c(i) + winsize2 > n
            continue;
        end
        
        temp = double(img(r(i) - winsize2:r(i) + winsize2-1, c(i) - winsize2:c(i) + winsize2-1));
        resp = zeros(ANGLE_DIVS,1);
        for t=1:ANGLE_DIVS
            resp(t) = sum(sum(temp.*abs_psi{t}));
        end
        response(ind(i)) = max(resp);
%        response(ind(i)) = abs(max(resp));
    end
end

% create r by combining everything in response, which is Composite
r = zeros(m,n);
for i=1:length(response)
    tmp = response(i);
    r = r + tmp{1};
end


%% create a number of images to show the results
r(r < 0) = 0;
resp_raw = r(r > 0);
r = mat2gray(r);
[I, J] = find(r>0.0);
r_255 = floor(r*255);
cmap = jet(256);
c_ind = zeros(length(I),3);
img_clr = repmat(img_raw, 1, 1, 3);
skel_clr = repmat(double(im_skel), 1, 1, 3);
skbr_clr = repmat(double(im_skel), 1, 1, 3);
for c=1:length(I)
    c_ind(c,:) = cmap(r_255(I(c),J(c))+1,:);
    r1 = I(c);
    c1 = J(c);
    img_clr(r1-3:r1+3,c1-3:c1+3,:) = repmat(reshape(c_ind(c,:),[1,1,3]),7,7,1);
    skel_clr(r1-3:r1+3,c1-3:c1+3,:) = repmat(reshape(c_ind(c,:),[1,1,3]),7,7,1);
    raw_skel(r1-3:r1+3,c1-3:c1+3,:) = repmat(reshape(c_ind(c,:),[1,1,3]),7,7,1);
    skbr_clr(r1-3:r1+3,c1-3:c1+3,:) = repmat(reshape([1 1 0],[1,1,3]),7,7,1);
end

% responses shown on top of image
imwrite(img_clr, fullfile(outdir, 'branch_points_resp.tif'), 'tif');
% responses shown on top of skeleton
imwrite(skel_clr, fullfile(outdir, 'skel_branch_points.tif'), 'tif');
% responses and skelton show on top of image
imwrite(raw_skel, fullfile(outdir, 'raw_skel_branch_points.tif'), 'tif');
% branch points (without response) show on top of image
imwrite(skbr_clr, fullfile(outdir, 'skel_branch_points_noresp.tif'), 'tif');
end
