Image denoising: Threshold calculation for oraclesoft

Senthilkumar July 30, 2011 Coded in Matlab
function thr = oracleshrink1(CD,T)
%function used to calculate threshold for oracleshrink method
CD = CD(:)';
n = length(CD);
sx2 = (T-CD).^2;
b  = cumsum(sx2);
[risk,best] = min(b);
hr = sqrt(sx2(best));

Image denoising - Oracleshrink method

Senthilkumar July 30, 20111 comment Coded in Matlab
%Using oracle threshold alone
%subband dependent threshold
%Soft and Hard threshold

function [soft_X1,SOFT_PSNR]  = Oracle_soft(X,Y)
%One -level decomposition
[CA,CH,CV,CD] = dwt2(Y,'haar');
[ca,ch,cv,cd] = dwt2(X,'haar');
%Call the function to calculate the threshold
T_CH  = func_threshold(CH);
T1_CH = oracleshrink1(ca,T_CH);
T_CV  = func_threshold(CV);
T1_CV = oracleshrink1(cv,T_CV);
T_CD  = func_threshold(CD);
T1_CD = oracleshrink1(cd,T_CD);
HH = uint8(T1_CD);
HL = uint8(T1_CV);
LH = uint8(T1_CH);
% Call the function to perfom soft shrinkage
de_CH = soft1(CH,T1_CH);
de_CV = soft1(CV,T1_CV);
de_CD = soft1(CD,T1_CD);
%
%Two -level decomposition
[CA1,CH1,CV1,CD1] = dwt2(CA,'haar');
[ca1,ch1,cv1,cd1] = dwt2(ca,'haar');
%Call the function to calculate the threshold
T_CH1  = func_threshold(CH1);
T1_CH1 = oracleshrink1(ch1,T_CH1);
T_CV1  = func_threshold(CV1);
T1_CV1 = oracleshrink1(cv1,T_CV1);
T_CD1  = func_threshold(CD1);
T1_CD1 = oracleshrink1(cd1,T_CD1);
HH1 = uint8(T1_CD1);
HL1 = uint8(T1_CV1);
LH1 = uint8(T1_CH1);
% Call the function to perfom soft shrinkage
de_CH1 = soft1(CH1,T1_CH1);
de_CV1 = soft1(CV1,T1_CV1);
de_CD1 = soft1(CD1,T1_CD1);
%
%
%Reconstruction for soft shrinkage
X2 = idwt2(CA1,de_CH1,de_CV1,de_CD1,'haar');
X1 = idwt2(X2,de_CH,de_CV,de_CD,'haar');

SOFT_PSNR = PSNR(X,X1);
soft_X1 = uint8(X1);

Image Denoising -threshold calculation using modified bivariate method

Senthilkumar July 30, 2011 Coded in Matlab
function [T,sigma1,sigma2] = model3_soft(CD,CH,CH1)
%function used to calculate the threshold and variance using modified 
%bivariate model
%CD and CH : obtained after first stage wavelet decomposition
%CH1: obtained after second stage 2D-wavelet decomposition
sigman = (median(median(abs(CD))))/0.6745;
sigmay1 = 0;
sigmay2 = 0;
[m,n] = size(CD);
[m1,n1] = size(CH1);
%To calculate sigma1 value
for i =1:m
    for j = 1:n
        sigmay1 = sigmay1+((CH(i,j))^2);
    end
end
sigmay1 = sqrt(2)*sigmay1/(m*n);
sigma1 = sqrt(max((((sigmay1))-((sigman)^2)),0));
if sigma1~=0
    T = sqrt(3)*(sigman^2);
else
    T = max(max(abs(CH)));
end
%To calculate sigma2 value
for i =1:m1
    for j = 1:n1
        sigmay2 = sigmay2+((CH1(i,j))^2);
    end
end
sigmay2 = sqrt(2)*sigmay2/(m1*n1);
sigma2 = sqrt(max((((sigmay2))-((sigman)^2)),0));

Image denoising - Bishrink method

Senthilkumar July 30, 2011 Coded in Matlab
% function [X1_denoise,DENOISE_PSNR] = Bivariate_soft(X1,Y)
clear all;
close all;
clc;
x = imread('lena.tif');
X =double(x);
[M,N] = size(X);
W =wgn(M,N,25);
Y = X+W;
    
[CA,CH,CV,CD] = dwt2(Y,'haar');
[CA1,CH1,CV1,CD1] = dwt2(CA,'haar');
[CA2,CH2,CV2,CD2] = dwt2(CA1,'haar');
[m,n] = size(CA);
[m2,n2] = size(CA1);
[m1,n1] = size(CA2);

%
k =1;
for i = 1:1:m2
    for j = 1:1:n2
        Y2_CH([2*i-1:2*i+k-1],[2*j-1:2*j+k-1])=CH1(i,j);
        Y2_CV([2*i-1:2*i+k-1],[2*j-1:2*j+k-1]) =CV1(i,j);
        Y2_CD([2*i-1:2*i+k-1],[2*j-1:2*j+k-1]) =CD1(i,j) ;
        
    end
end

for i = 1:1:m1
    for j = 1:1:n1
        Y1_CH([2*i-1:2*i+k-1],[2*j-1:2*j+k-1])=CH2(i,j);
        Y1_CV([2*i-1:2*i+k-1],[2*j-1:2*j+k-1]) =CV2(i,j);
        Y1_CD([2*i-1:2*i+k-1],[2*j-1:2*j+k-1]) =CD2(i,j) ;
        
    end
end
T1 = bishrink_threshold(CD,CH);
T2 = bishrink_threshold(CD,CV);
T3 = bishrink_threshold(CD,CD);
% 
w1_ch = bishrink(CH,Y2_CH,T1);
w1_cv = bishrink(CV,Y2_CV,T2);
w1_cd = bishrink(CD,Y2_CD,T3);

T4 = bishrink_threshold(CD,CH);
T5 = bishrink_threshold(CD,CV);
T6 = bishrink_threshold(CD,CD);
% 
w1_ch1 = bishrink(CH1,Y1_CH,T4);
w1_cv1 = bishrink(CV1,Y1_CV,T5);
w1_cd1 = bishrink(CD1,Y1_CD,T6);
% 
X2 = idwt2(CA1,w1_ch1,w1_cv1,w1_cd1,'haar');
X1_denoise = idwt2(X2,w1_ch,w1_cv,w1_cd,'haar');

%modified estimate is used

T1 = bishrink_threshold1(CD,CH);
T2 = bishrink_threshold1(CD,CV);
T3 = bishrink_threshold1(CD,CD);
% 
w1_ch = bishrink(CH,Y2_CH,T1);
w1_cv = bishrink(CV,Y2_CV,T2);
w1_cd = bishrink(CD,Y2_CD,T3);

T4 = bishrink_threshold1(CD,CH);
T5 = bishrink_threshold1(CD,CV);
T6 = bishrink_threshold1(CD,CD);
% 
w1_ch1 = bishrink(CH1,Y1_CH,T4);
w1_cv1 = bishrink(CV1,Y1_CV,T5);
w1_cd1 = bishrink(CD1,Y1_CD,T6);
% 
X4 = idwt2(CA1,w1_ch1,w1_cv1,w1_cd1,'haar');
X3_denoise = idwt2(X4,w1_ch,w1_cv,w1_cd,'haar');

noisy_PSNR = PSNR(X,Y)
% DENOISE_PSNR = PSNR(X,X1_denoise)
X1_PSNR = PSNR(X,X1_denoise)
X3_PSNR = PSNR(X,X3_denoise)
X1_denoise = uint8(X1_denoise);
X3_denoise = uint8(X3_denoise);
imshow(x)
figure;
imshow(uint8(Y))
figure;
imshow(X1_denoise)
figure;
imshow(X3_denoise)

Image denoising - Threshold calculation for Bishrink method

Senthilkumar July 30, 2011 Coded in Matlab
function T = bishrink2(CD,CY)
sigman = (median(median(abs(CD))))/0.6745;
[m,n] = size(CY);
sigmay = 0;
for i =1:m
    for j = 1:n
        sigmay = sigmay+((CY(i,j))^2);
    end
end
sigmay = sqrt(2)*sigmay/(m*n);
sigma = sqrt(max((((sigmay))-((sigman)^2)),0)); % Variance calculation
if sigma~=0
    T = sqrt(3)*(sigman^2)/sigma;   %Threshold
else
    T = max(max(abs(CY)));
end

Image denoising - Bishrink threshold calculation

Senthilkumar July 30, 2011 Coded in Matlab
function [w1] = bishrink1(y1,y2,T)
% Bivariate Shrinkage Function
%      y1 - a noisy coefficient value
%      y2 - the corresponding parent value
%      T  - threshold value
%      w1 - the denoised coefficient

R  = sqrt(abs(y1).^2 + abs(y2).^2);
R = R - T;
R  = R .* (R > 0);
w1 = y1 .* R./(R+T);

Image denoising -Bayeshrink threshold calculation

Senthilkumar July 30, 20112 comments Coded in Matlab
function T = bayeshrink(CD,CV)
[m,n] = size(CD);
sigma  = (median(median(abs(CD))))/0.6745; %Standard deviation
var_sig = std_sig^2;            % Variance
mean_sig = mean(mean(CH));      %mean value
sigmax = var_sig-(mean_sig^2);
sigmas = sqrt(max((sigmax-(sigma^2)),0));
%Calculation of threshold
if sigmas ~=0
    T = ((sigma)^2)/(sigmas);
elseif sigmas == 0
    T = max(max(abs(CV)));
end

Image Denoising - Bayesshrink

Senthilkumar July 30, 20113 comments Coded in Matlab
%Using BAYESSHRINK
%subband dependent threshold
%Soft  threshold
function [soft_X1,SOFT_PSNR] = Bayes_soft(X,Y)
%One -level decomposition
[CA,CH,CV,CD] = dwt2(Y,'haar');
%Call the function to calculate the threshold
CH =CH/1.2;
CV = CV/1.2;
CD = CD/1.2;
T1_CH = bayeshrink(CD,CH);
T1_CV = bayeshrink(CD,CV);
T1_CD = bayeshrink(CD,CD);
% Call the function to perfom soft shrinkage
de_CH = soft1(CH,T1_CH);
de_CV = soft1(CV,T1_CV);
de_CD = soft1(CD,T1_CD);
%Two -level decomposition
[CA1,CH1,CV1,CD1] = dwt2(CA,'haar');
CH1 = CH1/1.2;
CV1 = CV1/1.2;
CD1 = CD1/1.2;
%Call the function to calculate the threshold
T1_CH1 = bayeshrink(CD1,CH1);
T1_CV1 = bayeshrink(CD1,CV1);
T1_CD1 = bayeshrink(CD1,CD1);

% Call the function to perfom soft shrinkage
de_CH1 = soft1(CH1,T1_CH1);
de_CV1 = soft1(CV1,T1_CV1);
de_CD1 = soft1(CD1,T1_CD1);
%Reconstruction for soft shrinkage
X2 = idwt2(CA1,de_CH1,de_CV1,de_CD1,'haar');
X1 = idwt2(X2,de_CH,de_CV,de_CD,'haar');
SOFT_PSNR = PSNR(X,X1);
soft_X1 = uint8(X1);

Phase changing Allpass IIR Filter

July 24, 2011 Coded in Matlab
function [b,a] = allpass(order, Fc, Fs, Q)
% Returns allpass filter coefficients.
% Currently only 1st and 2nd orders are supported.
%
% Usage: [b,a] = ALLPASS(N, FC, FS, Q);
%
%        N is the order of the allpass
%        FC is the frequency a the 90deg phase shift point
%        FS is the sampling rate
%        Q is quality factor describing the slope of phase shift
%
% Author: sparafucile17 01/07/2004

if(order > 2)
    error('Only 1st and 2nd orders are supported');
end

%Bilinear transform
g = tan(pi*(Fc/Fs));

if(order == 2)
    d  = 1/Q;
    K  = 1/(1 + d*g + g^2);
    b0 = (1 - g*d + g^2) * K;
    b1 = 2 * (g^2 - 1) * K;
    b2 = 1;
    a1 = b1;
    a2 = b0;
else
    b0 = (g-1)/(g+1);
    b1 = 1;
    b2 = 0;
    a1 = b0;
    a2 = 0;
end

b = [b0 b1 b2];
a = [1  a1 a2];

Batch WAV file processing

July 24, 2011 Coded in Matlab
%In this example we are going to apply a low-pass filter to all WAV files,
%but it could be a multitude of different "processing" methods.  This is
%only used to illustrate the batch processing example.
Fs = 44100; %Hz
Fc = 1000;  %Hz
[b,a] = butter(2, Fc/(Fs/2), 'low');

%These are the input and output directories relative to this M-file
input_directory = 'test_input_database\';
output_direcory = 'test_output\';

%What extensions are you looking for?  In our case, we want to batch
%process audio files that are store in the uncompressed *.WAV format
extension = '*.wav';

%Get the files
files=dir([input_directory '*.wav']);
N_files=numel(files);

%Loop through one file at a time
for i=1:N_files

    %This is a stereo file and for this example, 
    %I only care about channel number 1
    x = wavread([input_directory files(i).name ]);
    x = x(:,1);
    
    %Process the data by applying our filter
    y = filter(b,a,x);
    
    %Output the data as a unique filename based on input
    wavwrite(y, Fs, [output_directory 'processed' files(i).name]);
    disp(['Processed File: ' input_directory files(i).name 'as: ' ...
                             output_directory 'processed' files(i).name]);
    
end