Thread Subject:

Find the region with the smallest mean value within a bigger matrix

Hi,
I would appreciate some hints regarding the following problem:
I am using a matrix of 512x76 values in which I am want to find a 40x40 sub-matrix that shows the lowest average value of all possible sub-matrices.
How can I achieve this?
Thanks a lot in advance,

Philipp

"Philipp Steffen" <philipp.steffen@imba.oeaw.ac.at> wrote in message <i6r7pv$h46$1@fred.mathworks.com>...
> Hi,
> I would appreciate some hints regarding the following problem:
> I am using a matrix of 512x76 values in which I am want to find a 40x40 sub-matrix that shows the lowest average value of all possible sub-matrices.
> How can I achieve this?
> Thanks a lot in advance,
>
> Philipp

If you have the Image Processing Toolbox:
help stdfilt

On 10-09-15 02:46 PM, Philipp Steffen wrote:
> Hi,
> I would appreciate some hints regarding the following problem:
> I am using a matrix of 512x76 values in which I am want to find a 40x40
> sub-matrix that shows the lowest average value of all possible
> sub-matrices.
> How can I achieve this?


"all possible submatrices" of a 40 x 40 submatrix is 2^1600 - 1 submatrices
which is about 10^333 of them. Are you sure you want to calculate the average
values of all of those ??

If so, then I offer a hint: over all possible subsets of a 40 x 40 submatrix,
the subsets that have the lowest average values are exactly those subsets that
consist of nothing except copies of the lowest value in the 40 x 40 matrix.
The searching for the lowest possible average value over all of the subsets
can thus be reduced to finding the minimum value over the 40 x 40 matrix.


If the "submatrices" was intended to refer to 40 x 40 solid blocks rather than
submatrices within the 40x40, then you need to tell us whether the submatrices
to be considered are non-overlapping or not. If they _are_ overlapping, there
are things you can do to reduce the calculations by using rolling sums --
subtract off the sum of the old first column add in the sum of the new last
column, if the total is greater than the previous total then this new
submatrix is _not_ the lowest average. If your values are not especially
delicately balanced, then you can also speed up the code by using cumsum()

Philipp Steffen:
I think Sean mean conv2() not stdfilt(), as in this demo I created:
Pay particular attention to the conv2() line, and the very last line
in the demo. BE SURE TO JOIN ANY LINES SPLIT INTO TWO BY THE
NEWSREADER.
-ImageAnalyst

clc; % Clear command window.
clear all; % Delete all variables.
close all; % Close all figure windows except those created by imtool.
imtool close all; % Close all figure windows created by imtool.
workspace; % Make sure the workspace panel is showing.
fontSize = 20;
% Change the current folder to the folder of this m-file.
if(~isdeployed)
cd(fileparts(which(mfilename)));
end
% Read in a standard MATLAB gray scale demo image.
folder = 'C:\Program Files\MATLAB\R2010a\toolbox\images\imdemos';
baseFileName = 'cameraman.tif';
fullFileName = fullfile(folder, baseFileName);
grayImage = imread(fullFileName);
% Get the dimensions of the image. numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
% Display the original gray scale image.
subplot(1, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
set(gcf, 'Position', get(0,'Screensize')); % Enlarge figure to full
screen.
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off');

% Find the local averages over all locations that a
% sliding 41 by 41 window can have.
% First find the convolution kernel.
windowSize = 41;
kernel = ones(windowSize, windowSize) / windowSize^2;
localMinImage = conv2(grayImage, kernel, 'same');
% Show the local average image.
subplot(1, 2, 2);
imshow(localMinImage, []);
title('Local Mean Image', 'FontSize', fontSize);

halfWindow = ceil(windowSize / 2)
% Find the valid subimage - which will have no edge effects
% and the window completely fits inside the original image.
subImage = localMinImage(halfWindow:rows-halfWindow+1,
halfWindow:columns-halfWindow+1);
% Find the lowest local average
minGL = min(min(subImage));
% Find row and column where this occurs
% and print it to the command window.
[minRows, minCols] = find(subImage == minGL)

Unfortunately I dont have the Image Processing Toolbox. But you are right, the purpose is indeed to find a region in an image (represented as a matrix) that has the lowest intensity values. Perhaps I was not quite clear about it. I dont really have to calculate the average of all possible sub-matrices. I would probably be fine with using a simpler approach that helps me to find the region with the lowest values. I will have a look at the rolling sum approach but I dont quite see how cumsum() should help me?
Thanks so far for your help. I might come back later with further questions.
Philipp

Then just write your own convolution/filtering routine - it's not that
hard, just 4 nested for loops in it's most intuitive yet brainless
form.

"Philipp Steffen" <philipp.steffen@imba.oeaw.ac.at> wrote in message <i6reqs$1eg$1@fred.mathworks.com>...
> Unfortunately I dont have the Image Processing Toolbox. But you are right, the purpose is indeed to find a region in an image (represented as a matrix) that has the lowest intensity values. Perhaps I was not quite clear about it. I dont really have to calculate the average of all possible sub-matrices. I would probably be fine with using a simpler approach that helps me to find the region with the lowest values.
=============

Is this what you want?



myMatrix=reshape(1:(572*76), 572,76); %Example 572x76 data
 

k=ones(1,40)/40; %1D averager

Averages=conv2(k,k,myMatrix,'valid'); %40x40 block averages

[M,N]=size(Averages);
Map=inf(size(myMatrix));
Map(1:M,1:N)=Averages;

[~,minBlock]=min(Map(:)), %minBlock is coordinate of upperleft corner of minimum average block

Matt:
From the help description, I don't see how that does it in blocks
(jumps) of 40 pixels (like blockproc would). It looks just like it's
doing a regular convolution, although it's forcing the kernel to be
separable. But it still goes along one pixel at a time I believe.

Anyway, I also mentioned conv2() and he said he didn't want it because
he doesn't have the image processing toolbox. However he could also
do this with just the standard conv() that comes in base MATLAB
because the 40x40 kernel is separable. He can just convolve in one
direction and then convolve with the same kernel rotated 90 degrees.
He'll need to convert his image to 1D row or column vectors before
calling conv() with the reshape command, but that's easy.

In case, Phillipp comes back, here is the version that doesn't
require the Image Processing Toolbox and uses the conv() function with
separable kernels:
(Looks bigger than it really is because of all the demo & tutorial
stuff to help explain it better).

clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
imtool close all; % Close all imtool figures.
clear; % Erase all existing variables.
workspace; % Make sure the workspace panel is showing.
fontSize = 14;

% Change the current folder to the folder of this m-file.
% (The line of code below is from Brett Shoelson of The Mathworks.)
if(~isdeployed)
cd(fileparts(which(mfilename)));
end
clc; % Clear command window.
clear all; % Delete all variables.
close all; % Close all figure windows except those created by
imtool.
imtool close all; % Close all figure windows created by imtool.
workspace; % Make sure the workspace panel is showing.
fontSize = 20;
% Change the current folder to the folder of this m-file.
if(~isdeployed)
        cd(fileparts(which(mfilename)));
end
% Read in a standard MATLAB gray scale demo image.
folder = 'C:\Program Files\MATLAB\R2010a\toolbox\images\imdemos';
baseFileName = 'cameraman.tif';
fullFileName = fullfile(folder, baseFileName);
grayImage = imread(fullFileName);
% Get the dimensions of the image. numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
% Display the original gray scale image.
subplot(1, 2, 1);
image(grayImage);
axis square;
title('Original Grayscale Image', 'FontSize', fontSize);
set(gcf, 'Position', get(0,'Screensize')); % Enlarge figure to full
screen.
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off');

% Find the local averages over all locations that a
% sliding 41 by 41 window can have.
% First find the convolution kernel.
windowSize = 41;
halfWindow = ceil(windowSize / 2)
kernel = ones(windowSize, 1) / windowSize;
% Turn into a column vector.
grayImage1D = single(reshape(grayImage, [rows*columns 1]));
localMeanImage = conv(grayImage1D, kernel);
localMeanImage = localMeanImage(halfWindow:end-halfWindow+1);
% Reshape back to a 2D matrix.
localMeanImage = reshape(localMeanImage, [rows columns]);
% Turn into a row vector.
grayImage1D = reshape(localMeanImage', [rows*columns 1]);
localMeanImage = conv(grayImage1D, kernel);
localMeanImage = localMeanImage(halfWindow:end-halfWindow+1);
% Reshape back to a 2D matrix.
localMeanImage = reshape(localMeanImage, [rows columns]);
% Transpose
localMeanImage = localMeanImage';
% Show the local average image.
subplot(1, 2, 2);
localMeanImage8 = uint8(localMeanImage); % Convert to uint8 for
display only.
image(localMeanImage8);
axis square;
cmap = linspace(0, 255, 256)/256;
grayMap = [cmap;cmap;cmap]';
colormap(grayMap);
title('Local Mean Image', 'FontSize', fontSize);

% Find the valid subimage - which will have no edge effects
% and the window completely fits inside the original image.
subImage = localMeanImage(halfWindow:rows-halfWindow+1,
halfWindow:columns-halfWindow+1);
% Find the lowest local average
minGL = min(min(subImage));
% Find row and column where this occurs
% and print it to the command window.
[minRows, minCols] = find(subImage == minGL)

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <8cbcc692-80e0-43f8-97bf-0c50d1543716@c16g2000vbp.googlegroups.com>...
> Matt:
> From the help description, I don't see how that does it in blocks
> (jumps) of 40 pixels (like blockproc would). It looks just like it's
> doing a regular convolution, although it's forcing the kernel to be
> separable. But it still goes along one pixel at a time I believe.
===========

Yes, but I see nothing in Philip's posts suggesting that he is looking for a blockproc-style search. I tlooks to me like he wants all 40x40 blocks searched. Regardless, my approach is easily modified to accomodate that.


> Anyway, I also mentioned conv2() and he said he didn't want it because
> he doesn't have the image processing toolbox.
============

conv2 doesn't require the Image Processing Toolbox, only regular MATLAB.

On Sep 15, 9:02 pm, "Matt J " <mattjacREM...@THISieee.spam> wrote:
> conv2 doesn't require the Image Processing Toolbox, only regular MATLAB.
--------------------------------------
You're right. (I thought it used to but maybe I was wrong). In that
case, he can use the conv2(). My apologies that I thought you were
suggesting a block-by-block approach. In this forum I just always
associate block with the blockproc way of doing things, and "kernel,"
"window," or "structuring element" with the pixel-by-pixel approach.

I am sorry about this, but I think I did not really make clear what my problem is.
Given the following matrix:
Matrix=[1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 1 2 1 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 0 0 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 0 0 0 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 1 0 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
        1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;]

I would need a method to find the coordinates of a 3x3 submatrix that has the lowest average or sum of values. In my example I would expect to find the following submatrix:

1 0 0
1 0 0
1 1 0

I would need for example the coordinates of the upper left corner. In my example this would therefore be Column #3, Row #6

So far I could not figure out which of your suggestions would be the most appropriate on?
Nevertheless I am really thankful for all the effort you put into helping me with this!

"Philipp Steffen" <philipp.steffen@imba.oeaw.ac.at> wrote in message <i6sk65$g6b$1@fred.mathworks.com>...
> I am sorry about this, but I think I did not really make clear what my problem is.

Oh it is clear, I think rather you don't pay enough attention to the hints ImageAnalyst and Matt have provided:

Matrix=[1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
1 1 1 1 2 1 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
1 1 1 0 0 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
1 1 1 0 0 0 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
1 1 1 1 0 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;
1 1 1 2 2 2 3 3 3 4 4 4 3 3 3 1 1 1 2 2 2;]

szblk=3;
meanMat=conv2(Matrix,ones(szblk),'valid');
[i j]=find(meanMat==min(meanMat(:)),1)
minblock=Matrix(i(1)+(0:szblk-1),j(1)+(0:szblk-1))

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <i6sl72$kuv$1@fred.mathworks.com>...
> "Philipp Steffen" <philipp.steffen@imba.oeaw.ac.at> wrote in message <i6sk65$g6b$1@fred.mathworks.com>...
> > I am sorry about this, but I think I did not really make clear what my problem is.
>
> Oh it is clear, I think rather you don't pay enough attention to the hints ImageAnalyst and Matt have provided:

I did not want to say that ImageAnalysts and Matts replies were not clear. Since I am a beginner I had troubles figuring out what the different approaches were and I could not clearly see which one would be best for my problem. Therefore I tried to clarify what I really aim for. But thanks a lot. Brunos code is exactly doing what I need.
Thank you all for your help!