% Improvements for global thresholding

clc;
close all; clear all;
%% Part1
f=imread('Fig1015(a)[noiseless].tif'); % bivalued image
fn=imnoise(f,'gaussian',0,0.038); % zero mean and 50 variance
% J=imnoise(I,Type,M,V); add noise of a given type to the image
% Type: 'gaussian','poisson',... and normalized numerical parameters
% M and V are 0 and 0.01 by default, V here = 50^2/255^2
figure; subplot(2,3,1); imshow(fn); title('noisy image');

subplot(2,3,2); imhist(fn); title('histogram'); % unimodal histogram
gn=im2bw(fn,graythresh(fn)); % Otsu's thresholding
subplot(2,3,3); imshow(gn); title('thresholded');
% segmentation is unsuccessful

h=fspecial('average',5); %  square size filtering, [3,3] by default
% h=fspecial('average',hsize); return an averaging filter h of size hsize
fa=imfilter(fn,h,'replicate'); % B=imfilter(A,H, Options)
% A and B are the same in the size and class
% 'replicate', input values outside bounds are computed by replication
subplot(2,3,4); imshow(fa); title('smoothed image');
subplot(2,3,5); imhist(fa); title('histogram'); % bimodal histogram
ga=im2bw(fa,graythresh(fa)); subplot(2,3,6); imshow(ga); title('segmentation');
% segmentation is significantly improved after averaging the raw image

%% Part2
% Improving thresholding by gradient-based edge detection
f=imread('Fig1016(a).tif');
figure; subplot(2,3,1); imshow(f); title('raw image');
subplot(2,3,2); imhist(f); title('raw histogram');
f=tofloat(f); % conversion to single class, scaled to [0,1];
sx=fspecial('sobel'); % H=fspecial(Type); create a filter mask
% Sobel detector computes the gradient using the mask of size 3*3
% By default,Sobel emphasizes on horizontal edges
sy=sx'; % verticcal edge, H=H';

% B=imfilter(A,B, options), B is computed using floating point
gx=imfilter(f,sx,'replicate'); % smoothing by a vertical gradient
gy=imfilter(f,sy,'replicate'); % smoothing by a horizontal gradient
grad=sqrt(gx.*gx+gy.*gy); % first-order derivative or gradient(magnitude)
grad=grad/max(grad(:)); % normalized to [0,1], floating point
h=imhist(grad);

% identify edge pixel,those which have large values of gradient
Q=percentile2i(h,0.999);% I=percentile2i(h,P), Q intensity level [0,1]
% using high percentile (0.999) return the intensity level near boders

% locate the indexes of edge  pixels
markerimage=grad>Q; % logical comparison, with edge pixels of value 1
subplot(2,3,3); imshow(markerimage); title('gradient thresholded binary image');
% some pixels selected are located in backgrounds

fp=f.*markerimage; % fp contains edge pixels and also some background pixels
subplot(2,3,4); imshow(fp); title('gradient thresholed gray-scale image');

% fp is dominated by pixels of 0 values
% the frequency for 0 should be removed for subsequent Otsu's thresholding
hp=imhist(fp); % frequency for all intensity values, column vector
hp(1)=0; % set the frequencye to 0 for the 0 intensity pixel
subplot(2,3,5); bar(hp,0); % bar(Y,width), specifies the bar width
title('gradient thresholded histogram');

% Otsu's thresholding given the histogram, otsuthresh
[T,SM]=otsuthresh(hp); % Again, thresholding on left pixels
% resulting separability = 0.89
% actual threshold level T=T*(numel(hp)-1); % T*255;

g=im2bw(f,T); % segmented with the resultant threshold
subplot(2,3,6); imshow(g); title('improved segmentation');

%% Part3
% Laplacian detector to obtain edge information for thresholding
f=imread('Fig1017(a).tif');
figure; subplot(2,3,1) ;imshow(f); title('raw image');
% bright spots in the objects immersed in background

f=tofloat(f); % single floating point, scaled to [0,1]
subplot(2,3,2); imhist(f); title('raw histogram');
[T_f,SM_f]=graythresh(f); % Otsu's thresholding on raw image
% SM_f=0.636; T=T_f*255=42
gf=im2bw(f,T_f);
subplot(2,3,3); imshow(gf); title('thresholded image');
% bright spots are not identified by directly thresholding on raw image

% Linear spatial filter: Laplacian filter
W=[-1 -1 -1; -1 8 -1; -1 -1 -1]; % with center 8 surrounded by -1
% Thefilter is manually set, since it is beyond those available in fspeciall
% W=fspecial('laplacian',alpha), a 3*3 laplacian fiter
% the filter shape is specified by alpha [0,1], 0.2 be default

lap=abs(imfilter(f,W,'replicate')); % using absolute values of lap
lap=lap/max(lap(:));

% identify the pixels that are near the egdes of bright spots
h=imhist(lap);
Q=percentile2i(h,0.995);
markerimage=lap>Q; % locate pixels on and within the borders

% show the identified pixels
fp=f.*markerimage; % in the raw image
subplot(2,3,4); imshow(fp); title('Lap-thresholded image');

hp=imhist(fp); % histogram of selected pixels
hp(1)=0; % remove the dominant pixels of value 0
subplot(2,3,5); bar(hp,0); title('Lap thresholded histogram');

% Again, Otsu's thresholding
[T,SM]=otsuthresh(hp);
% SM=0.762; T=0.451*255=115
g=im2bw(f,T);
subplot(2,3,6); imshow(g); title('improvement segmentation');