views:

86

answers:

3

What is a quick, elegant way of using MatLab to form a subarray around a particular element? Element are selected randomly from the data, so you can't take a subarray in the normal way (it has to be generalized for the elements that are selected). What I mean is, forming an array for example 5x5 or 7x7 or something, where the middle element is the one you want.


Thanks, but I mean a more general approach. Sorry for not being clear. What I mean is something like, let's say

a = magic(10)          
b = find(a<8) %giving linear indices     
m = size(a)      
[r,c] = ind2sub(m,b) %giving rows and columns of the elements that you later want to select

and then I thought about doing the rest using for loops but am not really sure how to generalize it and wondered, like I said, if there's an easier way to do it using MatLab.

The reason for such a general approach is I want to analyze any array for elements with certain properties and then be able to form a subarray around each element to basically zoom in on it for further data analysis if that makes sense.

+1  A: 

Here are some examples of how you can approach this problem:

mat = rand(10);  %# A 10-by-10 array of random elements

vec = -3:3;      %# To make a 7-by-7 subarray
subArray = mat(5+vec,4+vec);  %# Get subarray centered on element (5,4)

vec = -2:2;      %# To make a 5-by-5 subarray
subArray = mat(3+vec,3+vec);  %# Get subarray centered on element (3,3)

Of course, you'll have to be wary of conditions where your subscripts extend beyond the bounds of the matrix you are indexing, and then decide how you want to deal with such a situation (i.e. only return the part of the subarray that is within the larger matrix, or pad the edges of the larger matrix using PADARRAY to avoid such a situation, etc.).

EDIT:

To generalize the above for sets of points that you want to get subarrays around, the solutions given by Jonas should work well. The only drawback is that you would need to have access to the Image Processing Toolbox. Here's a solution that will work without additional toolboxes, using a, r, and c as you have defined them in the question:

blockSize = 3;                             %# Block size around each point
padSize = floor(blockSize/2);              %# Padding size for the matrix
aPad = padarray(a,[padSize padSize],NaN);  %# Pad the matrix with NaNs
rPad = r+padSize;                          %# Shift the row indices
cPad = c+padSize;                          %# Shift the column indices
blockIndices = (1:blockSize)-ceil(blockSize/2);  %# Create indices for a block
getBlockFcn = @(r,c) aPad(r+blockIndices,c+blockIndices);  %# Function to get
                                                           %#   a block
subArrays = arrayfun(getBlockFcn,rPad,cPad,...  %# Create a cell array with
                     'UniformOutput',false);    %#   one block per cell

This will give you an N-by-1 cell array subArrays, where N is the number of points you have. If, for example, you want to get the maximum value of the subarray for point 3, you can do the following:

maxValue = max(subArrays{3}(:));  %# MAX will ignore NaN values

To take a mean of the entire subarray, you either have to remove NaN values first (using the function ISNAN), or use the function NANMEAN from the Statistics Toolbox:

mat = subArray{3}(:);                %# Values from subarray
meanValue = mean(mat(~isnan(mat)));  %# Mean of non-NaN values
%# Or...
meanValue = nanmean(subArrays{3}(:));  %# Mean of non-NaN values
gnovice
+1  A: 

If you want to gather the individual (possibly overlapping) neighborhoods into an array so that you can process each neighborhood individually, you can use IM2COL. If you want to create a 'halo' around each element to easily capture all candidate regions, so that overlaps are not counted twice, you can use IMDILATE.

Example for im2col

%# create array and get linear index of interesting elements
a = magic(10);
b = find(a<8);

%# define region size
subSz = 3; %# for 3x3 array
halfSz = floor(subSz/2); %# for a 3x3 array, there window extends 1 pix to each side

%# pad a with NaNs so that sub-window never goes outside the image
ap = padarray(a,[halfSz halfSz],NaN);

%# convert ap to columns - there are 9 rows (for 3x3 window) and as many columns as there
%# are elements in a (100 in this example)
apCol = im2col(ap,[subSz,subSz]);

%# select columns of interest
%# the first column contains the first neighborhood of interest, the second column 
%# the second neighborhood etc.
interestingColumns = apCol(:,b);

%# take e.g. the average of each neighborhood
avg = nanmean(interestingColumns,1);

Example for imdilate

%# create array and get logical mask where 1 indicates an interesting element
a = magic(10);
bw = (a<8);

%# define region size
subSz = 3; %# for 3x3 array

%# grow each element, so that a single 'good' pixel is in the center of a 3x3 mask
bwDil = imdilate(bw,strel('square',subSz));

%# show the centers in white, 'halo' in grey
%# note the overlaps
figure,imshow(bw+bwDil,[])

%# take the average of all 'labeled' pixels
avg = mean(a(bwDil));
Jonas
+1: Nice solutions. The only limitation is the need for the Image Processing Toolbox to use the functions IM2COL and IMDILATE.
gnovice
With respect to the first solution, I think you forgot to adjust the indices `b` to take into account the padding added to the matrix.
gnovice
I think not. im2col with the default sliding window gives (nRowArray-nRowWindow+1)*(nColArray-nColWindow+1) columns. If you pad with half the window size rounded down, you get exactly nRowArray*nColArray sub-windows. Thus, you index directly with b.
Jonas
@Jonas: You're right. I mistakenly thought that the default block type used by IM2COL was 'distinct', when it's actually 'sliding'.
gnovice
A: 

Thanks. Both solutions are very nice, but for an array padded with NaN's how can you still perform operations on it such as mean2 or fit functions?

JJ
@JJ: You should add things like this as a comment on an answer, or edit your question to include them. You shouldn't post it as an answer.
gnovice
@JJ: Also, I added an example to my answer covering ways you can deal with NaN values.
gnovice
@JJ: My solution actually nan-pads the array, so it is designed to work, anyway. If you don't have the statistics toolbox, check gnovice's post.
Jonas