Blockproc not working as I expected

Hello, I am having problems with blockproc.

I have an black image with white spots that I want to basic threshold locally rather than globally - that is to perform a basic threshold on tiles.

IM=double(getimage(handles.axes4)); [rows,cols]=size(IM);     %Get image from axes component
IM=IM-min(IM(:));                                             %Background subtract
% Number of tiles or chunks you want it divide up into.
xtiles =4;  ytiles= 4; 

%Numbers pixels per tile to maintain integer number of tiles
nx=floor(cols/xtiles);    ny=floor(rows/ytiles);

Call my "LocallythresholdImage" function to define a level for binarising the image where

level = average + 2* standard deviations of the image

centroids=LocallythresholdImage(IM,nx,ny,f)            %f=2 i.e. number of standard deviatins above the average

function centroids=LocallythresholdImage(IM,m,n,f)
%Uses the thresh function below and then creates a binary image using local 
%threshold of size mxn. I
IM=double(IM);                %make sure its a double
 
%Use blockproc uses my "thresh" function
myfun = @(block_struct) thresh3(block_struct.data,f);      % myfun = @(block_struct) thresh2(block_struct.data,f);
%I2 = blockproc(Orig,[200 200],myfun,'BorderSize', [3 3], 'TrimBorder', true);
centroids = blockproc(IM,[m n],myfun,'BorderSize', [3 3], 'TrimBorder', true);

function [centroids]=thresh3(IM,f)

%n=str2num(get(handles.editThreshn,'String'));    
%Thresholds image  based on mean and S.D.    
    image=double(IM)-min(IM(:));    
    threshold=mean(image(:)) + f*std(image(:));
    lessThanThreshold = image < threshold;  
    imageThresholded = image; 
    imageThresholded(lessThanThreshold) = 0; 
    regmax = imregionalmax(imageThresholded);
    s = regionprops(regmax,double(image), 'WeightedCentroid');
    format bank
    centroids = cat(1, s.WeightedCentroid);

for each tile, the centroids are calculated correctly, but in my top level function "centroid=LocallythresholdImage(IM,nx,ny,f) ", centroids is empty

[Jason]

Maybe if I re-phrase my problem, someone may kindly comment.

When using blockproc, can the function fun that is applied to each block return an array of coordinates.

The aim is to plot the coordinates from each block on the original image.

Answer 1

[Guillaume]

Change your 'TrimBorder', true to 'TrimBorder', false. As it is you're removing a border of size 3 on each side of your output of thresh3.

Note that this thresh3 function is not well written. image should not be used as a variable name (it's a matlab function), format bank has no place in there (it only affects command line display) and the cat(1, ...) is completely wrong. I'm not sure what you're trying to do with that. There shouldn't be more than one centroid, but if there were, concatenating them all would result in an error in blockproc since the blocks would return arrays of different sizes.

[Walter Roberson]

The outputs of the function you pass into blockproc must be numeric or logical, but not necessarily the same class as the input.

You will get an error if you return different sizes of output for different blocks, with one exception: if processing the first block returns the empty matrix, then the final result will be the empty matrix even if processing other blocks would result in something non empty.

[Jason]

Hi Walter , thanks for your reply.

so each block (=one tile on my image) can have a different nmber of spots (the centroids ,im trying to find) so the number of centroids returned will be different. My aim is to plot on the composite image the centre of all the spots, so is this not possible using blockproc?

[Jason]

Hi Guillaume, thanks for your response.

1 - I will test your trim border suggestion

2 - I will not use image as variable name

3 - format bank was just for debugging so I can output troids in the command window.

4 - As each tile has several spots there will be more than 1 centroid hence i cat them.cen

5 - As I am returning an array of different size for each block is this my problem, is there no workaround?

thanks

[Walter Roberson]

Why not threshold, then bwareafilt, then do a regionprops asking for centroids, and the scatter() the centroids ? It is not obvious to me that you need to divide into blocks.

[Guillaume]

When you specify 'BorderSize', matlab trims the border from the input and output unless you specify 'TrimBorder' as false. Since you were removing 6 rows from an output with only one row, you were always going to get an empty output.

Yes, you cannot use blockproc to return arrays of different sizes. Considering the discussion in IA's answer, what I'd recommend is that instead your tresh3 function return the image you were passing to regionprops. You can then do the regionprops afterward on all the blocks at once. If you forget about the local background subtraction:

function filteredblock = thresh3(IM, f)
    threshold = mean(IM(:)) + f*std(IM(:));
    IM(IM < threshold) = 0;
    filteredblock = imregionalmax(IM);
end

function centroids = LocallythresholdImage(IM,m,n,f)
     filteredimage = blockproc(double(IM), [m n], @(block_struct) thresh3(block_struct.data, f), 'BorderSize', [3 3], 'TrimBorder', false);
     props = regionprops(filteredimage, IM, 'WeightedCentroid');
     centroids = vertcat(props.WeightedCentroid);
end

If you do need that local background subtraction, then return a 3D output (first page: the regional max, 2nd page, the local background subtraction image:

function filteredblock = thresh3(IM, f)
    IM = IM - min(IM(:));
    threshold = mean(IM(:)) + f*std(IM(:));
    thresholdimage = IM;
    thresholdimage(IM < threshold) = 0;
    filteredblock = cat(3, imregionalmax(thresholdimage), IM);
end

function centroids = LocallythresholdImage(IM, m, n, f)
    filteredimage = blockproc(double(IM), [m n], @(bs) thresh3(bs.data, f), 'BorderSize', [3 3], 'TrimBorder', false);
    props = regionpros(filteredimage(:, :, 1), filteredimage(:, :, 2), 'WeightedCentroid');
    centroids = vertcat(props.WeightedCentroid);
end

[Jason]

Thankyou, thats a really nice ay to get round the different size array issue - just return the binary image itself and then perform the threshold on that!

Jason

Answer 2

[Image Analyst]

"The aim is to plot the coordinates from each block on the original image."

Then don't complicate it. You can use a global threshold.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
% clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 18;

folder = pwd; % Current folder.
baseFileName = 'index.png';
% % Have user browse for a file, from a specified "starting folder."
% % For convenience in browsing, set a starting folder from which to browse.
% startingFolder = pwd;  % or 'C:\wherever';
% if ~exist(startingFolder, 'dir')
% 	% If that folder doesn't exist, just start in the current folder.
% 	startingFolder = pwd;
% end
% % Get the name of the file that the user wants to use.
% defaultFileName = fullfile(startingFolder, '01.png');
% [baseFileName, folder] = uigetfile(defaultFileName, 'Select a file');
% if baseFileName == 0
% 	% User clicked the Cancel button.
% 	return;
% end
fullFileName = fullfile(folder, baseFileName)

%===============================================================================
% Check if file exists.
if ~exist(fullFileName, 'file')
	% The file doesn't exist -- didn't find it there in that folder.
	% Check the entire search path (other folders) for the file by stripping off the folder.
	fullFileNameOnSearchPath = baseFileName; % No path this time.
	if ~exist(fullFileNameOnSearchPath, 'file')
		% Still didn't find it.  Alert user.
		errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
		uiwait(warndlg(errorMessage));
		return;
	end
end
rgbImage = imread(fullFileName);

% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(rgbImage)
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	% 	grayImage = rgb2gray(rgbImage);
	% ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
	% which in a typical snapshot will be the least noisy channel.
	grayImage = rgbImage(:, :, 3); % Take blue channel.
else
	grayImage = rgbImage; % It's already gray scale.
end
% Now it's gray scale with range of 0 to 255.

% Display the image.
subplot(1, 2, 1);
imshow(grayImage, []);
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
hp = impixelinfo();

%------------------------------------------------------------------------------
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
drawnow;

% Threshold.
mask = grayImage  == 255;
% Take the largest blobs only
mask = bwareaopen(mask, 100);
% Display the mask.
subplot(1, 2, 2);
imshow(mask);
hp = impixelinfo();
axis('on', 'image'); % Make sure image is not artificially stretched because of screen's aspect ratio.
title('Mask', 'FontSize', fontSize);

% Find the centroids
props = regionprops(mask, 'Centroid', 'Area');
allCentroids = vertcat(props.Centroid)
hold on;
plot(allCentroids(:, 1), allCentroids(:, 2), 'r+', 'MarkerSize', 30, 'LineWidth', 2);

You could do a little better if you worked on the image without the number annotation on it so that you didn't get the 6 and 16.

[Jason]

Walter / Image Analyst, thank you very much for you comments.

1: I love the bwareafilt function that Walter suggested, I had actually created my own area filtering function but using the built in function is more elegant. (its shown in the image below)

2: You both mention keep it simple and use Global thresholding. Well, I am building a GUI that allows for global thresholding & local thresholding (as well as other methods such as normalised cross correlation). The image I uploaded was a simple one where global thresholding works. However, I also need to analyse more complicated images as below which have spots much more closer together and some images have a non-uniform background, so I absolutely need a localised thresholding option.

(Sorry, there doesn't seem to be an option to shrink these images)

As you can see global thresholding is failing. Notice the bwareafilt together with imshowpair is a very nice way to show the thresholded objects that have been excluded via area classification (green blobs are excluded as too large)

I.A, in your example you nicely show bwareaopen, but often I need to go the other way and exclude the larger objects, I don't see a bwareaclose function to remove large objects from a binary image.

Also to add, the grid lines and number annotations are added after the processing and so don't affect the thresholding / binarisation.

So it still comes downto my problem. I want to break up my image into smaller tiles and do the thresholding and then reconstruct the centroids from each tile onto the main image. Im taking that this just isn't possible with blockproc which is a shame.

Thanks

[Image Analyst]

Jason, my code worked with your original example because you had well separated blobs that were all saturated (255), so global thresholding was sufficient. If your blobs are not saturated, then it wouldn't work.

There is no bwareaclose(), unless you wrote one. You can use bwareafilt() to explude larger blobs:

smallBlobs = bwareafilt(binaryImage, [1, maxAllowableBobSize]);

I'm attaching a program that cleans up a document by doing a local Otsu thresholding. It may be just what you want (at least the thresholding part of the program).

It looks like you could benefit by a higher resolution camera since some of those blobs are pretty blocky and only 3 pixels across.

There is an option to shrink your images. When you browse to it, turn the black twisty that says "display" downwards and select 25% or 50%, or type in something like 33%.