%% Setup - run this first

clear all; close all; clc; %#ok<CLALL> 

% MATLAB has a few built-in images that you can load to test with
rgb_pears = imread('pears.png');
rgb_chips = imread('coloredChips.png');

% Different formats can be read by the same function, but the data type
% that results can vary. 
I_rice = imread('rice.png');
I_pout = imread('pout.tif'); % Note that imread takes a variety of different extensions, but the data type isn't necessarily the same. 

% make a double image from a uint8 image.
d_rice = double(I_rice); d_rice=d_rice./max(d_rice(:));


%
%% Whole image quantifiers (single-image)

% Many functions can apply to entire images or pairs of images. Here we
% will look at single-image quantification and later will touch on
% referenced and referenceless image quality quantification.

% Gray level co-occurrence matrix is a matrix-based quantifier of image
% texture:
glcm_rice = graycomatrix(I_rice);
tiledlayout(1,2); nexttile; 
imagesc(I_rice); axis equal tight off; colorbar('south');
title('Raw image');
nexttile;
imagesc(log10(glcm_rice)); axis equal tight off ; colorbar('south');
title('log_{10}(glcm)');
colormap gray;
graycoprops(glcm_rice)

% Calculate the entropy of the rice image vs the pout image.
entropy_rice = entropy(I_rice)
entropy_pout = entropy(I_pout);

% Calculate the histograms of the rice image using multiple bins.
figure;
imhist(I_rice,256); xlabel('Intensity values'); ylabel('Counts'); title('256 bins');

figure;
imhist(I_rice,64); xlabel('Intensity values'); ylabel('Counts'); title('64 bins');


%
%% Whole image operations

% There are many built-in functions that filter entire images. The dropdown
% below is not an exhaustive list. 
filteredImage = imgaussfilt(I_rice);


tiledlayout(1,2,'TileSpacing','compact'); 
nexttile; imagesc(I_rice); title('Raw image'); axis equal tight off; colorbar('south');
nexttile; imagesc(filteredImage); title('Filtered image'); axis equal tight off; colorbar('south');


%% Running many different filtering functions
basicFilters = {'medfilt2','imboxfilt','imgaussfilt','rangefilt','entropyfilt','stdfilt'}

tiledlayout(2,3,'TileSpacing','compact');
N_filt = numel(basicFilters);
filtered_I_rice = cell(N_filt,1); 

% Call the built-in functions all in a row using `feval`
% This is the same as calling medfilt2(I_rice), imboxfilt(I_rice), etc.
% but saves coding.
for k = 1:N_filt
    nexttile;
    filtered_I_rice{k} = feval(basicFilters{k},I_rice);
    imagesc(filtered_I_rice{k}); colorbar; axis equal off;
    title(basicFilters{k});
end

%

% There are also filter kernels, which are created using the fspecial
% function and applied with imfilter. The following list is non-exhaustive.

specialFilters = {'average','disk','gaussian','laplacian','log','sobel'}
N_filt = numel(specialFilters);

% Generating many kernels and applying them to the image.
filtered_I_rice = cell(N_filt,1);

for k = 1:N_filt
    figure;
    tiledlayout(1,2);
    
    nexttile;
    filt = fspecial(specialFilters{k});
    imagesc(filt); axis equal tight; title([specialFilters{k} ' kernel']);
    
    nexttile;
    filtered_I_rice{k} = feval('imfilter',I_rice,filt);
    imagesc(filtered_I_rice{k}); colorbar; axis equal off;
    title([specialFilters{k} ' applied']);
end



%
%% Edge Filtering

I_rc = d_rice(1:64,1:64);

figure;
imagesc(I_rc); axis square off; title('Base image'); 

% One can filter to smooth noise, but this can ruin resolution by smoothing edges
GAUSSIAN_SIGMA = 1.3;
I_rc_smooth = imgaussfilt(I_rc,GAUSSIAN_SIGMA); % Large standard deviation of Gaussian.
imagesc(I_rc_smooth); axis equal tight; title(['Gaussian \sigma=',num2str(GAUSSIAN_SIGMA)]); axis off;


% One can filter for edges 'manually', by using classical filters.
% Formally, this is the application of a kernel (made by 'fspecial') to a
% target image through an operation known as convolution.

sobelY = fspecial('sobel');
sobelX = sobelY.'; 

I_edge_X = imfilter(I_rc,sobelX);
I_edge_Y = imfilter(I_rc,sobelY);
tiledlayout(1,2); 
nexttile; imagesc(I_edge_X); axis equal tight off; title({'X edge intensity';'sobel imfilter'}); 
nexttile; imagesc(I_edge_Y); axis equal tight off; title({'Y edge intensity';'sobel imfilter'}); 

figure;
sobelMag_hypot = hypot(I_edge_X,I_edge_Y);
imagesc(sobelMag_hypot); axis equal tight off; title({'Edge magnitude';'sobel imfilter'});

% Notice that the previous approach, though technically correct and
% appropriate for most of the image, has the problem that it introduces
% edge artifacts. 
% imgradient can do the same in one step, and handles boundaries correctly.
sobelMag = imgradient(I_rc); 
imagesc(sobelMag); axis equal tight off; title({'Edge magnitude';'imgradient'});

[I_edge_X_grad ,I_edge_Y_grad] = imgradientxy(I_rc);
tiledlayout(1,2); 
nexttile; imagesc(I_edge_X_grad); axis equal tight off; title({'X edge intensity';'imgradientxy'}); 
nexttile; imagesc(I_edge_Y_grad); axis equal tight off; title({'Y edge intensity';'imgradientxy'}); 


%

% Noise can make edge filtering unstable, magnifying the noise.
figure;
tiledlayout(1,2);
nexttile;
I_rc_noise = I_rc + 50*imnoise(I_rc);
sobelMag_noise = imgradient(I_rc_noise);
imagesc(sobelMag_noise); axis equal tight off; title({'Edge magnitude';'noise'}); colorbar('south');


% Median filtering can be a good way to remove some types of noise. Notice corner
% artifacts.
nexttile;
sobelMag_medfilt = imgradient(medfilt2(I_rc_noise));
imagesc(sobelMag_medfilt); axis equal tight off; title({'Edge magnitude';'noise and medfilt2'}); colorbar('south');

% It is also possible to filter images using edge-aware methods:
smoothing = 50;
I_bilat = imbilatfilt(I_rc,smoothing);
I_bilat_noise = imbilatfilt(I_rc_noise,smoothing);

tiledlayout(2,2,'TileSpacing','compact');
nexttile;
imagesc(I_rc); axis equal tight off; title('Raw image');
nexttile;
imagesc(I_rc_noise); axis equal tight off; title('Noisy image');
nexttile;
imagesc(I_bilat);  axis equal tight off; title({'Bilateral filt of raw';['Smoothing: ' num2str(smoothing)]});
nexttile;
imagesc(I_bilat_noise); axis equal tight off; title({'Bilateral filt of noisy';['Smoothing: ' num2str(smoothing)]});

% The wiener2 function does spatially-adaptive filtering to smooth out
% noise while retaining edge quality, by solving a particular inverse
% problem related to the local variance of a given pixel.
tiledlayout(1,2);
nexttile;
imagesc(I_rc_noise); axis equal tight off; title('Noisy image');
nexttile;
I_rc_noise_wiener = wiener2(I_rc_noise);
imagesc(I_rc_noise_wiener); axis equal tight off; title('wiener2 filtered image');


%
%% Deep Learning Denoising

% Recent MATLAB versions include pretrained networks which can be used for
% a variety of applications. Their performance can vary heavily depending
% on the data input to them. 

net = denoisingNetwork('DnCNN');
I_dnn_denoised = denoiseImage(I_rc_noise,net);

tiledlayout(1,2);
nexttile;
imagesc(I_rc_noise); axis equal tight off; title('Noisy image');
nexttile;
imagesc(I_dnn_denoised); axis equal tight off; title('DNN Denoised image');


%
%% Image comparison and referenced quantification

% If you have a known ground-truth image, you can quantify image quality
% using reference-based image metrics. Different metrics convey different
% information - some more structural, some more noise-related.

% Normalized MSE - a smaller value means the image is closer to the ground
% truth. This is symmetric, so the arguments to immse can be in either
% order. Also note that the 'perfect' image has an MSE of 0.
immse_noise_raw = immse(I_rc_noise,I_rc)./numel(I_rc)
immse_wiener_raw = immse(I_rc_noise_wiener,I_rc)./numel(I_rc)
immse_raw_raw = immse(I_rc,I_rc)./numel(I_rc);

% Structural similarity (ssim). This is asymmetric, and the reference image
% must be the second argument. Note that the 'perfect' image has an SSIM of
% 1, and values of SSIM are between 0 and 1.
ssim_noise_raw = ssim(I_rc_noise,I_rc)
ssim_wiener_raw = ssim(I_rc_noise_wiener,I_rc)
ssim_raw_raw = ssim(I_rc,I_rc);

%
%% Reference-less quantification

% It is also possible to use reference-less image quality metrics. These
% are often trained on particular image models and so may not be
% appropriate for your application. As always, read the documentation. 

piqe_raw = piqe(I_rc)
piqe_noise = piqe(I_rc_noise)
piqe_noise_wiener = piqe(I_rc_noise_wiener)
piqe_bilat = piqe(I_bilat)

[score,activity,artifacts,noise] = piqe(I_rc);


brisque_raw = brisque(I_rc)
brisque_noise = brisque(I_rc_noise)
brisque_noise_wiener = brisque(I_rc_noise_wiener)
brisque_bilat = brisque(I_bilat)



