Threshold

Threshold Image

Threshold: 

Threshold follows the same concept as in basic electronics, here it is used to convert the grayscale image to black and white(binary image, consisting of 0's and 1's as pixel values.)
Threshold for an image can be set in the following ways:

• Brute force : .5 times the maximum pixel intensity, which will be 127 in most of the images.
• Otsu's method : Based on a paper published in 1979 (IEEE Link), this method makes use of the histogram shape to find the threshold for the image. This gives a good enough performance but there are a few drawbacks when we run on entire image. There are a few regions that miss out as it will not consider the local neighborhood of the image. So, I prefer to use block based processing as this splits the image into smaller blocks and this can help us in getting the good idea in neighborhood for calculating threshold.

The Matlab code is given below:

clc;clear all;close all;

%Local Function Definitions:
fun = @(block_struct) thresher(block_struct.data);

I1= imread('rice.png');

f1 = figure(1);
imshow(I1);
title('Actual image');

brute = im2bw(I1,0.5);
brute = bwareaopen(brute,50);
imshow(brute);
title('Brute force Threshold image');
imwrite(brute,'brute.jpg','jpeg');
pause();

Brute Force

thresh = graythresh(I1);
otsu_full_image = im2bw(I1,thresh);
otsu_full_image = bwareaopen(otsu_full_image,30);
imshow(otsu_full_image);
title('Otsu`s Threshold on Full Image');
imwrite(otsu_full_image,'otsu_full_image.jpg','jpeg');
pause();

Otsu's Threshold (Full Image)

 
 block_otsu = blockproc(I1,[32,32],fun);
block_otsu = bwareaopen(block_otsu,30);
imshow(block_otsu);
title('Otsu`s Threshold on block Image');
imwrite(block_otsu,'block_otsu.jpg','jpeg');
pause();

Otsu's Threshold (Block Processing)

thresher.m (Function used during block processing)

function [ out ] = thresher(in)

thresh = graythresh(in);
out = im2bw(in,thresh);
imshow(out);

end

Note that the edges near lower corner have been jagged in the first and second method. This is not seen in block processed image. In the brute force threshold image we can observe the grain missed out in the top right part of the image.

Now let's see what is block processing using the 2 lines that define it.

fun = @(block_struct) thresher(block_struct.data);
block_otsu = blockproc(I1,[32,32],fun);

 
Here we define that we should split the image I1 into smaller blocks of 32X32 image and then run the function 'fun' on it. This 'fun' was defined in the second line of the program that it has to send the data under process to the 'thresher' function. The thresher function calculates the Otsu's threshold for that 32X32 image and returns a 32X32 binary image calculated with threshold value calculated for that block. Thus we have used the ability of 32X32 neighborhood processing.Note: If block size is smaller it can take more time to run and too small block size may not always produce good results. So, its a trade-off.

In blockproc 'UseParallel',true can help sometimes if you have Parallel Computing Toolbox.

Basics of Digital Image

The blog had been quiet for a long time, I am restarting it again with a different perspective. This time around I would be learning from basics and will try to explain the stuff as I understand from the basics. If I have misunderstood any concept, I will be pleased to learn from anyone.

I will be using Gonzalez and Woods 3rd Ed. for this. I will try to give minimal explanation on Matlab code but stress more on the concepts around as Matlab Help is extensively available at www.mathworks.com.

Without wasting much time lets begin....

Few terms to begin with

Bare Basics:

Grayscale Image: A frame of 2 D distribution of intensities.
Coordinates - x,y
Any image is a product component of intensity-i(x,y) and reflectance-r(x,y)
Intensity: 0 to inf (Yes, infinity and that's how you write in Matlab).
Reflectance: 0 to 1
(Think, what would be reflectance of a 1. black screen? 2.white screen? 3. transparent plastic sheet?)

How to see in Matlab:

in = imread('rice.png');
whos in;
size(in)

The second line here whos gives variable name, size of image read into, bytes in storage, datatype per element. If you just give whos, you end up getting the details about all the variables in workspace.

This means in black and white movies you have 3 dimensions (x, y, time).

Color Image: A triplet frame of 2 D distribution of intensities. (I will explore more in color models later.)
Coordinates - x, y, z
How to see in Matlab:


in = imread('color.jpg'); %make sure there is a color image of %same name in the folder
whos in;
size(in)

This means in color movies in a basic sense you have 4 dimensions (x, y, z, time).

Image Sampling and Quantization:

Analog image to Digital image conversion uses this.
Discretization in time domain is image sampling and digitization of amplitudes after discretization is quantization. End result "A digital image".

Sampling intervals in x-axis gives your length resolution and y-axis gives you the height resolution of a picture.
More the resolution, better defined is the image for visual perception, but more is the size for storage..

Vector images that are used in Photoshop are handled differently (as to what I know the images there are defined by mathematical equations, so can be extrapolated/expanded without loss of resolution/chequer-boarding). A good read on this is in Wikipedia here...

Hope I have started off well, please post questions so that I can get different perspectives...

Point detection, Laplacian of Gaussian and High Boost Filtering

As with other posts, remove the commenting part in the below code to see the code working.
Pre-requisite You know imread, imshow and other functions in MATLAB.

clc;

%% Read Image
    %Reads the Image into an array for further processing
X = imread('build.tif');
    %X = rgb2gray(X);
   
%% Point detection
Y = imread('turbine.tif');
%Y1= bwareaopen(Y, 20);
% fplap = [1 1 1; 1 -8 1; 1 1 1]      %Laplacian at a point
fplap = [0 1 0; 1 -4 1; 0 1 0]      %Laplacian at a point
filtim = imfilter(Y,fplap,'symmetric', 'conv');

subplot(2,2,1)
imshow(Y);
title('Original');

subplot(2,2,2)
imshow(filtim);
title('Laplacian Point Filtered');

subplot(2,2,3)
imshow(Y-filtim);
title('Difference Image');

%% Gaussian Filtering of image
% fgauss = fspecial('gaussian',[25,25],4.0)
%         % 25X25 Gaussian filter with SD =25 is created.
% filtim = imfilter(X,fgauss,'symmetric', 'conv');
%         % Filter the image by convolution with the above designed filter.
%         % filtim now contains the gaussian filtered image.
%        
% subplot(2,2,1)
% imshow(X);
% title('Original');
%
% subplot(2,2,2)
% imshow(filtim);
% title('Gaussian Filtered');
%
% subplot(2,2,3)
% imshow(X-filtim);
% title('Difference Image');
%
% % flap1 = fspecial('laplacian', .2);
% % filtim1 = imfilter(filtim,flap1,'symmetric', 'conv');
%
% % fsobel = [-1 0 1; -2 0 2; -1 0 1]      %Grad in y
% fsobel = [-1 -2 -1; 0 0 0; 1 2 1]      %Grad in x
% % fsobel = [0 -1 -2; 1 0 -1; 2 1 0]        %45 degree
% filtim1 = imfilter(filtim,fsobel,'symmetric', 'conv');
%
% subplot(2,2,4)
% imshow(filtim1);
% title('LoG Image');

%% Spatial Filtering Smoothing <Weighted> &High Boost Filtering
%
% favg = [1/16 2/16 1/16; 2/16 4/16 2/16; 1/16 2/16 1/16]
% filtim = imfilter(X,favg,'symmetric', 'conv');
%
% subplot(2,2,1)
% imshow(X);
% title('Original');
%
% subplot(2,2,2)
% imshow(filtim);
% title('Weighted Filtering');
%
% subplot(2,2,3)
% imshow(X-filtim);
% title('Difference Image');
%
% subplot(2,2,4)
% imshow(X+0.2.*(X-filtim));
% title('High Boost Sharpened');

%%

Spatial Filtering etc

Assistance for Chapter 3 initial parts from Gonzalez and Woods,

Usage Help: Un-comment(Ctrl+T in Matlab) the required part of filtering and Run(F5) to understand the fun.
<Please keep a copy of a picture named as (1).jpg in the Matlab folder.

clc;

%% Read Image

X = imread('(1).jpg');

X = rgb2gray(X);

%% Histogram Equalisation of Image

% 

% figure;

% imhist(X);

% figure;

% histeq(X);

% figure;

% 

% imhist(histeq(X));

%% Gaussian Filtering of image

% fgauss = fspecial('gaussian',[3,3],1.0)

% filtim = imfilter(X,fgauss,'symmetric', 'conv');

% 

% subplot(1,2,1)

% imshow(X);

% subplot(1,2,2)

% imshow(filtim);

% figure

% imshow(X-filtim);

%% Spatial Filtering Smoothing <Weighted> &High Boost Filtering

% 

% favg = [1/16 2/16 1/16; 2/16 4/16 2/16; 1/16 2/16 1/16]

% filtim = imfilter(X,favg,'symmetric', 'conv');

% 

% subplot(2,2,1)

% imshow(X);

% title('Original');

% 

% subplot(2,2,2)

% imshow(filtim);

% title('Weighted Filtering');

% 

% subplot(2,2,3)

% imshow(X-filtim);

% title('Difference Image');

% 

% subplot(2,2,4)

% imshow(X+0.2.*(X-filtim));

% title('High Boost Sharpened');

% 

%% Median Filtering <Doubt, Anybody who knows please assist in this, right now the code is substituting %Median of first block for all the blocks>

% 

% filtim = imfilter(X,(double(median(median(X)))),'symmetric', 'conv');   %takes a single value for all the %matrices, how to make it floating

% 

% subplot(2,2,1)

% imshow(X);

% title('Original');

% 

% subplot(2,2,2)

% imshow(filtim);

% title('Median Filtering');

% 

% subplot(2,2,3)

% imshow(X-filtim);

% title('Difference Image');

%% Laplacian

% flap1 = fspecial('laplacian', .2);

% filtim1 = imfilter(X,flap1,'symmetric', 'conv');

% flap2 = -flap1;

% 

% subplot(2,2,1)

% imshow(X);

% title('Original');

% 

% subplot(2,2,2)

% imshow(filtim1);

% title('Laplacian Filtering - 1');

% 

% filtim2 = imfilter(X,flap2,'symmetric', 'conv');

% subplot(2,2,3)

% imshow(filtim2);

% title('Laplacian Filtering - 2');

% 

% subplot(2,2,4)

% imshow(X+(filtim1));

% title('Reconstructed after Laplacian');

%%