One simple way would be a threshold image using some known value after converting the image to grayscale. The problem with this approach is that we apply the global threshold, and therefore some of the paper at the bottom of the image will be lost if you make the threshold too high. If you set the threshold too low, you will definitely get paper, but you will also turn on a lot of background pixels, and it will probably be difficult to remove these pixels with further processing.
One thing I can offer is to use the adaptive threshold algorithm. The algorithm that worked for me in the past is the Bradley-Roth adaptive threshold algorithm . You can read about it here in a post that I commented a while ago:
Bradley Adaptive Thresholding - Confused (questions)
However, if you want its essence, first make an integral image in black and white version of the image. The integrated image is important because it allows you to calculate the sum of the pixels in the window in complexity O(1) . However, the calculation of the integral image is usually O(n^2) , but you only need to do this once. With an embedded image, you look at the neighborhood of sxs pixels and you check to see if the average intensity is less than t% actual average in this sxs window, then this is a pixel classified as background. If it is larger, then it is classified as part of the foreground. This is adaptive because the threshold value is executed using local pixels and not using the global threshold.
I have encoded for you the implementation of the Bradley-Roth algorithm. The default parameters for algorithm s are 1/8 of the image width, and t is 15%. Therefore, you can simply call it in such a way as to call the default parameters:
out = adaptiveThreshold(im);
im is the input image and out is a binary image indicating that it belongs to the foreground ( logical true ) or background ( logical false ). You can play with the second and third input parameters: s is the size of the threshold window and t percent, which we talked about above, and we can call the function as follows:
out = adaptiveThreshold(im, s, t);
Therefore, the code for the algorithm is as follows:
function [out] = adaptiveThreshold(im, s, t) %// Error checking of the input %// Default value for s is 1/8th the width of the image %// Must make sure that this is a whole number if nargin <= 1, s = round(size(im,2) / 8); end %// Default value for t is 15 %// t is used to determine whether the current pixel is t% lower than the %// average in the particular neighbourhood if nargin <= 2, t = 15; end %// Too few or too many arguments? if nargin == 0, error('Too few arguments'); end if nargin >= 4, error('Too many arguments'); end %// Convert to grayscale if necessary then cast to double to ensure no %// saturation if size(im, 3) == 3 im = double(rgb2gray(im)); elseif size(im, 3) == 1 im = double(im); else error('Incompatible image: Must be a colour or grayscale image'); end %// Compute integral image intImage = cumsum(cumsum(im, 2), 1); %// Define grid of points [rows, cols] = size(im); [X,Y] = meshgrid(1:cols, 1:rows); %// Ensure s is even so that we are able to index the image properly s = s + mod(s,2); %// Access the four corners of each neighbourhood x1 = X - s/2; x2 = X + s/2; y1 = Y - s/2; y2 = Y + s/2; %// Ensure no co-ordinates are out of bounds x1(x1 < 1) = 1; x2(x2 > cols) = cols; y1(y1 < 1) = 1; y2(y2 > rows) = rows; %// Count how many pixels there are in each neighbourhood count = (x2 - x1) .* (y2 - y1); %// Compute row and column co-ordinates to access each corner of the %// neighbourhood for the integral image f1_x = x2; f1_y = y2; f2_x = x2; f2_y = y1 - 1; f2_y(f2_y < 1) = 1; f3_x = x1 - 1; f3_x(f3_x < 1) = 1; f3_y = y2; f4_x = f3_x; f4_y = f2_y; %// Compute 1D linear indices for each of the corners ind_f1 = sub2ind([rows cols], f1_y, f1_x); ind_f2 = sub2ind([rows cols], f2_y, f2_x); ind_f3 = sub2ind([rows cols], f3_y, f3_x); ind_f4 = sub2ind([rows cols], f4_y, f4_x); %// Calculate the areas for each of the neighbourhoods sums = intImage(ind_f1) - intImage(ind_f2) - intImage(ind_f3) + ... intImage(ind_f4); %// Determine whether the summed area surpasses a threshold %// Set this output to 0 if it doesn't locs = (im .* count) <= (sums * (100 - t) / 100); out = true(size(im)); out(locs) = false; end
When I use your image and I set s = 500 and t = 5 , here is the code and this is the image I get:
im = imread('http://i.stack.imgur.com/MEcaz.jpg'); out = adaptiveThreshold(im, 500, 5); imshow(out);
You can see that there are some false white pixels in the lower white image, and there are some holes that we need to fill inside the paper. So, let’s use some morphology and declare a structuring element that has a 15 x 15 square, make an opening to remove the noisy pixels, and then fill the holes when we finish:
se = strel('square', 15); out = imopen(out, se); out = imfill(out, 'holes'); imshow(out);
This is what I get after all of this:
Not bad, huh? Now, if you really want to see how the image looks with segmented paper, we can use this mask and multiply it by the original image. Thus, any pixels that belong to the paper are saved, and those that belong to the background go away:
out_colour = bsxfun(@times, im, uint8(out)); imshow(out_colour);
We get the following:
You will have to play with the parameters until they work for you, but the above parameters were the ones I used to make it work on the page that you showed us. Image processing is all about cork and error, as well as processing the processing steps in the correct sequence until you get something good enough for your purposes.
Happy image filtering!