Matlab image processing functions VERSUS my image processing functions in Matlab

Main Code:

%% My BitGet

tstart=tic;

% Get bit what do you want.

mybitget = @(in,pos) abs( bitand(in,(2.^uint8(pos-1))) );

% Guartied to not cleared old webcam objected.

clear('cam')

cam = webcam;

F=rgb2gray(snapshot(cam));

%F = gpuArray(F);

figure('Name','My BitGet: High Contrast');

subplot(4,4,[1 8]);imshow(F);title('Original 8 bit image');  %subplot original Image

for i=1:8

    % subplot(4,4,8+i);imshow(logical(bitget(F,i))); title(['Bitplane ',num2str(i)]); %bitget

    subplot(4,4,8+i);imshow(logical(mybitget(F,i))); title(['Bitplane ',num2str(i)]); %mybitget

end

% Guartied to clear webcam object

clear('cam')

tElapsed = toc(tstart);

disp(tElapsed)

High Contrast

Low Contrast

%% My Contrast Stretching

% use help to read report and help

Stretching(imread('img.png'), 0, 255, 30, 4, 255, 4);

Stretching(imread('low_contrast.png'), 0, 255, 30, 4, 220, 4);

High Contrast

Low Contrast

%% Matlab Contrast Stretching

A=rgb2gray(imread('img.png'));

A_F=imadjust(A); %output default dynamic range is 0-255

A_EF=imadjust(A,[0.3 0.7],[]); %output default dynamic range is 0-255

figure('Name','Matlab Contrast Stretching: High Contrast');

subplot(3,2,1); imshow(A); title('Original');

ax1=subplot(3,2,2); imhist(A); title('Original Histogram');

subplot(3,2,3); imshow(A_F); title('imadjust(A)');

ax2=subplot(3,2,4); imhist(A_F); title('imadjust(A) Histogram');

subplot(3,2,5); imshow(A_EF); title('imadjust(A,[0.3 0.7],[])');

ax3=subplot(3,2,6); imhist(A_EF); title('imadjust(A,[0.3 0.7],[]) Histogram');

linkaxes([ax1 ax2 ax3],'xy');

B=rgb2gray(imread('low_contrast.png'));

B_F=imadjust(B); %output default dynamic range is 0-255

B_EF=imadjust(B,[0.3 0.7],[]); %output default dynamic range is 0-255

figure('Name','Matlab Contrast Stretching: Low Contrast');

subplot(3,2,1); imshow(B); title('Original');

ax1=subplot(3,2,2); imhist(B); title('Original Histogram');

subplot(3,2,3); imshow(B_F); title('imadjust(B)');

ax2=subplot(3,2,4); imhist(B_F); title('imadjust(B) Histogram');

subplot(3,2,5); imshow(B_EF); title('imadjust(B,[0.3 0.7],[])');

ax3=subplot(3,2,6); imhist(B_EF); title('imadjust(B,[0.3 0.7],[]) Histogram');

linkaxes([ax1 ax2 ax3],'xy');

High Contrast

Low Contrast

%% My Histogram Equalization

A=rgb2gray(imread('img.png')); % dynamic range is 0-255

A_F=hist_eq(A,100); %output default dynamic range is 0-100

A_EF=hist_eq(A); %output default dynamic range is 0-255

A_AEF=hist_eq(A,400); %output default dynamic range is 0-255 but corrupted

figure('Name','My Histogram Equalization: High Contrast');

nx=4;

ny=2;

subplot(nx,ny,1); imshow(A); title('Original');

ax1=subplot(nx,ny,2); imhist(A); title('Original Histogram');

subplot(nx,ny,3); imshow(A_F); title('hist_eq(A,100)');

ax2=subplot(nx,ny,4); imhist(A_F); title('hist_eq(A,100) Histogram');

subplot(nx,ny,5); imshow(A_EF); title('hist_eq(A)');

ax3=subplot(nx,ny,6); imhist(A_EF); title('hist_eq(A) Histogram');

subplot(nx,ny,7); imshow(A_AEF); title('hist_eq(A,400)');

ax4=subplot(nx,ny,8); imhist(A_AEF); title('hist_eq(A,400) Histogram');

linkaxes([ax1 ax2 ax3 ax4],'xy');

A=rgb2gray(imread('low_contrast.png')); % dynamic range is 0-255

A_F=hist_eq(A,100); %output default dynamic range is 0-100

A_EF=hist_eq(A); %output default dynamic range is 0-255

A_AEF=hist_eq(A,400); %output default dynamic range is 0-255 but corrupted

figure('Name','My Histogram Equalization: Low Contrast');

nx=4;

ny=2;

subplot(nx,ny,1); imshow(A); title('Original');

ax1=subplot(nx,ny,2); imhist(A); title('Original Histogram');

subplot(nx,ny,3); imshow(A_F); title('hist_eq(A,100)');

ax2=subplot(nx,ny,4); imhist(A_F); title('hist_eq(A,100) Histogram');

subplot(nx,ny,5); imshow(A_EF); title('hist_eq(A)');

ax3=subplot(nx,ny,6); imhist(A_EF); title('hist_eq(A) Histogram');

subplot(nx,ny,7); imshow(A_AEF); title('hist_eq(A,400)');

ax4=subplot(nx,ny,8); imhist(A_AEF); title('hist_eq(A,400) Histogram');

linkaxes([ax1 ax2 ax3 ax4],'xy');

High Contrast

Low Contrast

%% Matlab Histogram equalization

A=rgb2gray(imread('img.png')); % dynamic range is 0-255

A_F=histeq(A,100); %output default dynamic range is 0-100

A_EF=histeq(A); %output default dynamic range is 0-255

A_AEF=histeq(A,400); %output default dynamic range is 0-255 but corrupted

figure('Name','Matlab Histogram Equalization: High Contrast');

nx=4;

ny=2;

subplot(nx,ny,1); imshow(A); title('Original');

ax1=subplot(nx,ny,2); imhist(A); title('Original Histogram');

subplot(nx,ny,3); imshow(A_F); title('histeq(A,100)');

ax2=subplot(nx,ny,4); imhist(A_F); title('histeq(A,100) Histogram');

subplot(nx,ny,5); imshow(A_EF); title('histeq(A)');

ax3=subplot(nx,ny,6); imhist(A_EF); title('histeq(A) Histogram');

subplot(nx,ny,7); imshow(A_AEF); title('histeq(A,400)');

ax4=subplot(nx,ny,8); imhist(A_AEF); title('histeq(A,400) Histogram');

linkaxes([ax1 ax2 ax3 ax4],'xy');

A=rgb2gray(imread('low_contrast.png')); % dynamic range is 0-255

A_F=histeq(A,100); %output default dynamic range is 0-100

A_EF=histeq(A); %output default dynamic range is 0-255

A_AEF=histeq(A,400); %output default dynamic range is 0-255 but corrupted

figure('Name','Matlab Histogram Equalization: Low Contrast');

nx=4;

ny=2;

subplot(nx,ny,1); imshow(A); title('Original');

ax1=subplot(nx,ny,2); imhist(A); title('Original Histogram');

subplot(nx,ny,3); imshow(A_F); title('histeq(A,100)');

ax2=subplot(nx,ny,4); imhist(A_F); title('histeq(A,100) Histogram');

subplot(nx,ny,5); imshow(A_EF); title('histeq(A)');

ax3=subplot(nx,ny,6); imhist(A_EF); title('histeq(A) Histogram');

subplot(nx,ny,7); imshow(A_AEF); title('histeq(A,400)');

ax4=subplot(nx,ny,8); imhist(A_AEF); title('histeq(A,400) Histogram');

linkaxes([ax1 ax2 ax3 ax4],'xy');

 

%% My Laplacian Filtering

image = imread('cameraman.tif');

mylaplacian_filter(image,3,3,8);

%% Matlab Laplacian Filtering center 8

I = imread('cameraman.tif');

H = [1, 1, 1;

    1,  8, 1;

    1, 1, 1]; %laplace kernel

lapla = imfilter(I,H,'replicate'); %laplace filter

figure('Name','Matlab Laplacian Filtering');

nrow=2; ncol=2;

subplot(nrow,ncol,1); imshow(I); subplot(nrow,ncol,2); imhist(I);

subplot(nrow,ncol,3); imshow(lapla); subplot(nrow,ncol,4); imhist(lapla);

%% Matlab gauss_AND_laplace filter(with new function) gauss: 3x3, ?2 = 5 ; laplace: 3x3 center -4

% just gives black image

I = imread('cameraman.tif');%load

Iblur = imgaussfilt(I, sqrt(5), 'FilterSize',[3 3]); %gauss filter

figure('Name','Matlab gauss_AND_laplace filter(with new function) gauss: 3x3, ?2 = 5 ; laplace: 3x3 center -4');

nrow=3; ncol=2;

subplot(nrow,ncol,1); imshow(I); title('Original Image'); subplot(nrow,ncol,2); imhist(I)

subplot(nrow,ncol,3); imshow(Iblur); title('Gaussian filtered image, \sigma = sqrt(5)'); subplot(nrow,ncol,4); imhist(Iblur)

H = [0, 1, 0;

    1,  -4, 1;

    0, 1, 0]; %laplace kernel

gauss_AND_lapla = imfilter(Iblur,H,'replicate'); %laplace filter

subplot(nrow,ncol,5); imshow(gauss_AND_lapla); title('Gauss and laplace'); subplot(nrow,ncol,6); imhist(gauss_AND_lapla)

%% Matlab gauss_AND_laplace filter(with new function) gauss: 3x3, ?2 = 1 ; laplace: 3x3 center -8

% just gives black image

I = imread('cameraman.tif'); %load

Iblur = imgaussfilt(I, 1, 'FilterSize',[3 3]); %gauss filter

figure('Name','Matlab gauss_AND_laplace filter(with new function) gauss: 3x3, ?2 = 1 ; laplace: 3x3 center -8');

nrow=3; ncol=2;

subplot(nrow,ncol,1); imshow(I); title('Original Image'); subplot(nrow,ncol,2); imhist(I)

subplot(nrow,ncol,3); imshow(Iblur); title('Gaussian filtered image, \sigma = 1'); subplot(nrow,ncol,4); imhist(Iblur)

H = [0, 2, 0;

    2,  -8, 2;

    0, 2, 0]; %laplace kernel

gauss_AND_lapla = imfilter(Iblur,H,'replicate'); %laplace filter

subplot(nrow,ncol,5); imshow(gauss_AND_lapla); title('Gauss and laplace'); subplot(nrow,ncol,6); imhist(gauss_AND_lapla)

 

%% Fourier Spectrum and Average Value

I = rgb2gray(imread('4_41.png')); %load

D=double(I);

D_fft_abs=abs( fftshift(fft2(D)) );

% Orjinal, D fft abs shift, D log fft abs shift ---AND--- Histograms

figure('Name','Fourier Spectrum and Average Value')

nrow=3; ncol=4; index=1;

subplot(nrow,ncol,index); imshow(I); title('Orjinal'); index=index+1; ax1=subplot(nrow,ncol,index,'align'); imhist(I); index=index+1; title('Orjinal Histogram');

subplot(nrow,ncol,index); imshow(D_fft_abs,[]); title('D fft abs shift'); index=index+1; ax2=subplot(nrow,ncol,index); imhist(D_fft_abs); index=index+1; title('D fft abs shift Histogram');

subplot(nrow,ncol,index); imshow(log(1+D_fft_abs),[]); title('D log fft abs shift'); index=index+1; ax3=subplot(nrow,ncol,index); imhist(log(1+D_fft_abs)); index=index+1; title('D log fft abs shift Histogram');

linkaxes([ax1,ax2,ax3],'xy');

% Orjinal average, I average fft abs shift, I average log fft abs shift ---AND--- Histograms

average = fspecial('average',3);

I_average = imfilter(I,average);

I_average_fft_abs = abs( fftshift(fft2(I_average)) );

subplot(nrow,ncol,index); imshow(I_average,[]); title('I average fft shift phase'); index=index+1; ax11=subplot(nrow,ncol,index); imhist(I_average); index=index+1; title('I average fft shift phase Histogram');

subplot(nrow,ncol,index); imshow(I_average_fft_abs,[]); title('I average fft shift abs'); index=index+1; ax12=subplot(nrow,ncol,index); imhist(I_average_fft_abs); index=index+1;  title('I average fft shift abs Histogram');

subplot(nrow,ncol,index); imshow(log(1+I_average_fft_abs),[]); title('I average log fft abs shift'); index=index+1; ax13= subplot(nrow,ncol,index); imhist(log(1+I_average_fft_abs)); title('I average log fft abs shift Histogram');

linkaxes([ax11,ax12,ax13],'xy');

 

%% Gaussian Lowpass Filter

I = rgb2gray(imread('4_41.png')); %load

Iblur = imgaussfilt(I, sqrt(5), 'FilterSize',[3 3]); %gauss filter

figure('Name','Gaussian Lowpass Filter(with new function) gauss: 3x3, \sigma = sqrt(5)');

nrow=2; ncol=2;

subplot(nrow,ncol,1); imshow(I); title('Original Image'); ax1=subplot(nrow,ncol,2); imhist(I)

subplot(nrow,ncol,3); imshow(Iblur); title('Gaussian filtered image, \sigma = sqrt(5)'); ax2=subplot(nrow,ncol,4); imhist(Iblur)

linkaxes([ax1,ax2],'xy');

Functions:

My Contrast Stretching:

function [ A_F, A_EF ] = Stretching( image, smin, smax, minpixelval, percentagemin, maxpixelval, percentagemax, do_you_want_to_plot, do_you_want_to_save_images  )

% Contrast Stretching

% image

% smin=0; % output image minimum histogram point

% smax=255; % output image maximum histogram point

% minpixelval=5; % Minimum threshold value

% maxpixelval=250; % Maximum threshold value

% percentagemin=1; % How much percentage bigger than minpixelval

% percentagemax=1; % How much percentage bigger than maxpixelval

%

%

% REPORT OF THIS FUNCTION:

% it has not need to eliminate pixels if image is really low contrast.

% Otherwise image getting worse

% Controlling inputs

switch nargin

    case 1

    image=imread('low_contrast.png');

    smin=0;

    smax=255;

    minpixelval=5;

    percentagemin=250;

    maxpixelval=1;

    percentagemax=1;

    do_you_want_to_plot=true;

    do_you_want_to_save_images=true;

    warning('default parameters');

    case 2

    smin=0;

    smax=255;

    minpixelval=5;

    percentagemin=250;

    maxpixelval=1;

    percentagemax=1;

    do_you_want_to_plot=true;

    do_you_want_to_save_images=true;

    warning('default parameters');

    case 3

    smax=255;

    minpixelval=5;

    percentagemin=250;

    maxpixelval=1;

    percentagemax=1;

    do_you_want_to_plot=true;

    do_you_want_to_save_images=true;

    warning('default parameters');

    case 4

    minpixelval=5;

    percentagemin=250;

    maxpixelval=1;

    percentagemax=1;

    do_you_want_to_plot=true;

    do_you_want_to_save_images=true;

    warning('default parameters');

    case 5

    percentagemin=250;

    maxpixelval=1;

    percentagemax=1;

    do_you_want_to_plot=true;

    do_you_want_to_save_images=true;

    warning('default parameters');

    case 6

    maxpixelval=1;

    percentagemax=1;

    do_you_want_to_plot=true;

    do_you_want_to_save_images=true;

    warning('default parameters');

    case 7

    do_you_want_to_plot=true;

    do_you_want_to_save_images=true;

    warning('default parameters');

    case 8

    do_you_want_to_save_images=true;

    warning('default parameters');       

end

if ~isnumeric(image) || ~isnumeric(minpixelval) || ~isnumeric(percentagemin) || ~isnumeric(maxpixelval) || ~isnumeric(percentagemax)

    error('input must be numeric')

end

if all(image(:)<0) || all(minpixelval<0) || all(percentagemin<0) || all(maxpixelval<0) || all(percentagemax<0)

    error('input must be bigger than zero')

end

if all(minpixelval<minpixelval)

    error('minimum pixel input must be smaller than maximum pixel input')

end

% Imread and find min and max of image. uint8 operators results are wrong.

F= double( rgb2gray( image ) ); % GpuArray maybe become get faster

%F=imresize(F,[480 640]); % Images are too big for showing

% - % - % Eliminate Outliers problems

A_EF=F;

Numel_A_EF = numel(A_EF(:));

MaskMin_PixelVal = A_EF<=minpixelval; % Find smaller than minpixelval. smallest pixel value

MaskMax_PixelVal = A_EF>=maxpixelval; % Find smaller than maxpixelval. biggest pixel value

masked_Min_A_EF=A_EF(MaskMin_PixelVal);

masked_Max_A_EF=A_EF(MaskMax_PixelVal);

if 100*numel( masked_Min_A_EF)/Numel_A_EF >=percentagemin % Is MaskMin_PixelVal percentage bigger than percentagemin

    A_EF(MaskMin_PixelVal) = minpixelval;

end

if 100*numel( masked_Max_A_EF)/Numel_A_EF >=percentagemax % Is MaskMax_PixelVal percentage bigger than percentagemax

    A_EF(MaskMax_PixelVal) = maxpixelval;

end

% - % - %Appling formula

A_F=F;

A_Fmin = min(A_F(:));  % Minimum pixel of image X

A_Fmax = max(A_F(:)); % Maximum pixel of image X

A_EFmin = min(A_EF(:));  % Minimum pixel of image X

A_EFmax = max(A_EF(:)); % Maximum pixel of image X

A_F= (A_F - A_Fmin).*abs( (smax-smin)/(A_Fmax-A_Fmin) )+smin; % Just using the formula

A_EF= (A_EF - A_EFmin).*abs( (smax-smin)/(A_EFmax-A_EFmin) )+smin; % Just using the formula

F=uint8(F);

A_F=uint8(A_F);

A_EF=uint8(A_EF);

% Plots

if do_you_want_to_plot

figure('Name','My Contrast Stretching');

subplot(3,2,1); imshow(F); title('Original');

ax1=subplot(3,2,2); imhist(F); title('Original Histogram');

subplot(3,2,3); imshow(A_F); title('My Contrast Stretching');

ax2=subplot(3,2,4); imhist(A_F); title('My Contrast Stretching Histogram');

subplot(3,2,5); imshow(A_EF); title('Eliminated My Contrast Stretching');

ax3=subplot(3,2,6); imhist(A_EF); title('Eliminated My Contrast Stretching Histogram');

linkaxes([ax1 ax2 ax3],'xy');

end

if do_you_want_to_save_images

imwrite(A_F,'Stretched.png');

imwrite(A_EF,'Stretched_Eliminated.png');

end

% - % - % Comparing eliminated and non eliminated Contrast Stretching

% If writes 1 Contrast Stretching is don't work without elimination otherwise works

fprintf('F == A_F is converged is %d.\n\n', all(F(:)==A_F(:)) );

fprintf('A_F == A_EF is converged is %d.\n\n', all(A_F(:) == A_EF(:)) );

fprintf('F == A_EF is converged is %d.\n\n', all(F(:) == A_EF(:)) );

% Report of experiment

disp('it has not need to eliminate pixels if image is really low contrast.');

disp('Otherwise image getting worse');

end

My Laplacian Filter:

function [g, H]=mylaplacian_filter(image,C,D,center)

% Report working but result difeerent from matlab function.

% image is 2ND matrix

% C kernel rows

% D kernel cols

% center is in the centered value of kernel.

% nr of rows/columns of kernel h . 3 by 3 filter or whatever you want.

image=double(image);

M=size(image,1); N=size(image,2); % nr of rows/columns of image f

P=M+C-1; Q=N+D-1; % nr of rows/columns after padding

imagep=zeros(P,Q); % zero padding: start with zeroes

imagep(1:M,1:N)=image; % insert image into image imagep

hp=zeros(P,Q); % Construct filter matrix hp, same size as imagep.

hp(1,1)=-center; hp(2,1)=1; hp(1,2)=1; % Center is at (center,center)

hp(P,1)=1; hp(1,Q)=1; % Indices modulo P or Q

imagep=fft2(double(imagep), P, Q); % FFT of image imagep

Hp=fft2(double(hp), P, Q); % FFT of kernel hp

H = fftshift(Hp);% Swaps the left and right halves of Hp

image1 = abs(H); % Get the magnitude

image1 = log(image1+1); % Use log, for perceptual scaling, and +1 since log(0) is undefined

image1 = mat2gray(image1,[0 1]);% Converts the matrix image1 to the intensity image image1. I can use gpu

Gp=imagep .* Hp; % Product of FFTs

gp=ifft2(Gp); % Inverse FFT

gp=real(gp); % Take real part

g=gp(1:M, 1:N);

gnorm = g;

gshar = double(image) - gnorm;

%Plots

figure('Name','mylaplacian_filter')

subplot(4,2,1); imshow(image,[]); title('Orjinal'); ax1=subplot(4,2,2); imhist(uint8(image));

subplot(4,2,3); imshow(image1); title('Kernel in Frequency domain');  ax2=subplot(4,2,4); imhist(uint8(image1));

subplot(4,2,5); imshow(255-g,[]); title('Filtered image Real part');  ax3=subplot(4,2,6); imhist(uint8(g));

subplot(4,2,7);imshow(uint8(gshar)); title('What is difference');   ax4=subplot(4,2,8); imhist(uint8(gshar));

linkaxes([ax1,ax2,ax3,ax4],'xy');

end

My Histogram Equation:

function [out] = hist_eq(im, L)

% Formula:

% g_i = floor( (L-1) * sum_{n=0}^{i} p_i )

% function gives better contrast range

% output default dynamic range is 0-L

%

% REPORT OF THIS FUNCTION:

% It has not need to histogram equalization pixels if image is really low contrast.

% Otherwise image getting worse

% Reference: R. C. Gonzalez and R. E.Woods, Digital Image Processing, Third Edition,

% 2008.

% Controlling Inputs

if nargin<1

    L=256;

    im=rgb2gray(imread('low_contrast.png'));

elseif nargin==1

    L=256;

end

if L>256

    warning('L>256 and image goint to corrupted');

end

h = imhist(im) / numel(im);

cdf = cumsum(h);

out = (L-1)*cdf(double(im)+1); % im is uint8

out = uint8(out); % out is double

end