find contiguous region with highest values

I'm trying to get the indexes of the region covering 10% of an m x n graph that has the highest values.

For example:
0 0 0 0 0 3
2 2 2 0 0 0
0 0 0 0 0 1
0 0 0 0 1 2
0 0 0 0 0 0

Should give me rows 2,2,2 and columns 1,2,3 (the output format doesn't matter for my purposes).

The region needs to be contiguous, and not just the region with the single highest value. Could someone point me in the right direction?

[amen45]

Essentially im looking for a reliable version of this puzzler, but that searches for largest group of values rather than equal values : http://blogs.mathworks.com/videos/2008/08/18/puzzler-find-largest-connected-island/?s_tid=srchtitle that I can adjust to pick out the 10% region

[Walter Roberson]

The maximum pixel of the 1 2 1 "L" is the same as the maximum pixel of the 2 2 2 line. Shouldn't that be a tie? Or is it the total over the region that is important?

[amen45]

I'm sorry about the confusion! It's the total over the region. I think perhaps it's better if I explain the actual problem I'm trying to solve.

I have a series of images that are 20 x 200 with floating point values generated by a previous analysis (not set by me). I've identified a 30x 3 "pixel" region (ie all values at indices (7:9,60:90) that I want to evaluate across all my figures). Within that region, I want to find the "cluster" pixels of size ~5 that has the highest values. They need to be next to each other, because in my data a regional relationship is required for them to be part of the same phenomenon (so pulling out the k highest values won't work). Again, the regional value is the important part so finding the location of the maximum values won't work either.

I have some preliminary code that I'll post in a separate comment. Please let me know if I can clarify anything else?

[Walter Roberson]

This does not answer to the question of what happens for a tie, and does not answer to the question of how you determine which cluster has the "highest values" across clusters of different sizes and where simple maximum is not the key. Is the key the mean over the cluster? Is the key the sum over the cluster?

Answer 1

[Chaya N]

There may be a nice (but slightly complicated) way of getting the linear indices here. The idea is something like this: Get connected components ( blobs shall we say?) --> Get pixel index list --> For each set of pixels, compute their pixel sum --> Find the index of the maximum pixelsum --> Get linear indices for the blob with the maximum pixelsum.

The code would be something like this:

 x = [0 0 0 0 0 3      % your example here
      2 2 2 0 0 0
      0 0 0 0 0 1
      0 0 0 0 1 2
      0 0 0 0 0 0]; 
 CC = bwconncomp(x); % Blobs and their data
 SumPixels = cell2mat(cellfun(@(y) sum(x(y)),CC.PixelIdxList,'UniformOutput',0)); % Get Pixelsum
 [~,idx] = max(SumPixels);
 LinIdx = CC.PixelIdxList{idx};

And Voila! I get:

LinIdx =

       2
       7
      12

You could find the row and column indices from the above, but this is where I leave you.

PS: For pedagogical reasons, I suggest you look up the bwconncomp and cellfun functions. To learn more about linear indices, please look at the "More About" section here .

[amen45]

I think this doesn't quite do what i need- it's not a binary array, and i actually don't even always have ints. any thoughts?

[Chaya N]

It would depend on how you separate your regions of interest. How do you define that certain pixels or a certain contiguous region belongs to your object of interest and not the background?

Also, you could always adapt/modify your non-integer-valued array long enough to get your program running. You could even change all your background values to zeros with a logical statement. My point is, the entire logic of the code depends on where the 'edges' of your regions lie, and a handful of pixels more or less makes a big difference to the region that you would obtain.

Answer 2

[Image Analyst]

Try ismember() and you can do it in 2 lines. See this well commented demo:

% Define sample data consisting of integers
% (Slightly different than your example)
m = [0 0 0 0 0 3
  2 2 2 0 0 0
  0 0 0 0 0 1
  0 3 3 0 1 2
  0 0 0 0 0 0]
% 3 is the "highest value" and is what is sought.
% 3 occurs in two regions (places).
% Find those two regions of highest value
% and put locations into a logical array.
locations = ismember(m, max(m(:)))
% Get a new matrix withteh actual values.
highestValueRegions = m .* locations

Answer 3

[Image Analyst]

If you have floating point values, I think what you're trying to say is that you want a matrix with only elements that are in the highest 10% of all values in the array. You can do this like:

% Define sample data consisting of floating point numbers
m = 100 * rand(6, 5)
% Sort the array
sortedm = sort(m(:), 'descend')
% Find the 10% index
index10 = ceil(0.1 * numel(m))
% Determine the value that 10% of the values lie above
value10 = sortedm(index10)
% Get binary "map" of where these highest 10% of values live:
binaryImage = m >= value10
% Extract a matrix with only the highest 10% of floating point values:
output = m .* binaryImage

When I ran it, this is what I got. Tell me if this is what you were thinking of:

m =
       9.1113       64.762       23.623       77.029       25.644
       57.621       67.902        11.94       35.022       61.346
       68.336       63.579        60.73       66.201       58.225
       54.659       94.517       45.014       41.616       54.074
       42.573       20.893       45.873       84.193       86.994
       64.444       70.928       66.194       83.292       26.478
sortedm =
       94.517
       86.994
       84.193
       83.292
       77.029
       70.928
       68.336
       67.902
       66.201
       66.194
       64.762
       64.444
       63.579
       61.346
        60.73
       58.225
       57.621
       54.659
       54.074
       45.873
       45.014
       42.573
       41.616
       35.022
       26.478
       25.644
       23.623
       20.893
        11.94
       9.1113
index10 =
     3
value10 =
       84.193
binaryImage =
  6×5 logical array
   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0
   0   1   0   0   0
   0   0   0   1   1
   0   0   0   0   0
output =
            0            0            0            0            0
            0            0            0            0            0
            0            0            0            0            0
            0       94.517            0            0            0
            0            0            0       84.193       86.994
            0            0            0            0            0

[Chaya N]

I think the OP wants a code that would produce the indices for a contiguous region in an array that occupies at least 10% of the array (by size) and whose sum of array values is a maximum among such regions as those would satisfy the size condition.

[Image Analyst]

Yes, it's ambiguous. What you say would require that regions are defined first, so that their area can be measured with regionprops. This would require that the regions be labeled. Originally he provided an integer image which you could then assume might be labels, or at least could be thresholded and labeled to produce a labeled image. However he then said that the "image" could be anything, not just integers, so it could be any arbitrary floating point array. There are no regions in such an array. You'd first have to define what constitutes a region, for example by thresholding at some level. My Image Segmentation Tutorial goes over the steps in well commented detail, but, in short, it would go like this:

binaryImage = floatingPointImage > someThreshold;
labeledImage = bwlabel(binaryImage);
props = regionprops(labeledImage, 'Area');
allAreas = [props.Area];
area10Percent = 0.1 * numel(binaryImage ); % The 10% area.
bigAreaIndexes = find(allAreas > area10Percent); % Which blobs are bigger than 10% of image.
bigAreaImage = ismember(labeledImage, bigAreaIndexes); % Extract only the big blobs.

Now bigAreaImage is a binary image showing every blob that has an area at least as big as 10% of the whole image.

[Chaya N]

I agree about the definition of a 'region' as it relates to this question. Your code snippet above renders the original array as a binary one and might as well be treated as an image, which is sort of the point I was trying to get across earlier in my answer.

PS : Your tutorial seemingly links to the Enigma Machine on Wikipedia! Perhaps that was an oversight?

[Image Analyst]

Sorry. Yes (copy apparently didn't work after an earlier copy for an earlier post). Fixed. Thanks for the catch.

[Image Analyst]

Let's hope amen45 comes back and clears up what "region with highest values" means exactly. Like does it mean "region with the biggest area" or "region with the greatest numerical values" or whatever. And also hopefully he'll give a floating point (non-integer) example and tell us where the "regions" in it are, and what data he wants out.

[amen45]

i mean region of size k with the greatest numerical values. For example, if k =4 and the array is as below:

0.115754943948525  0.254835774407905  0.418797235529874
0.133596512865232  0.267704112017176  0.427967388090852
0.148849791641271  0.274078404255313  0.425934632215066
0.166527249110295  0.977300000000000  0.998700000000000
0.196980352387009  0.923140000000000  0.987890000000000
0.224172402663761  0.323278947773164  0.431610265045853
0.283912660670702  0.380237101362344  0.466341084111908

What I'm trying to pull out are the indexes (in any format) and mean values of the 4 adjacent pixels with the highest numerical values. In the example above, those pixels are at (4:5,2:3)- that is the "region" I would like to identify.

Does this answer your question? Thank you all for all your help!

[amen45]

A rough example:

sets = [2,2;4,1];

for round = 1:size(sets,1)
    mrow = sets(round,1);
    mcol = sets(round,2);
    for cindex = 1:(size(temp_holder,2)-(mcol-1))
        for rindex = 1:(size(temp_holder,1)-(mrow-1))
            mean_holder(rindex,cindex,round) = mean2(temp_holder(rindex:rindex+(mrow-1),cindex:cindex+(mcol-1)));
        end
    end
    %     end
end

max_mean_value = (max(max(max(mean_holder))));

Here I'm presetting the different array shape combinations to go through the image and store mean values of that shape array at each possible location. However, ideally this would be automatically generated for all possible array shapes for array element number k that fit in the region we're looking in.

[Image Analyst]

The greatest value is 0.998700000000000. It occurs at row 4, column 3. It has only three 4-connected neighbors, and five 8-connected neighbors. Or four and six if you also include the value itself. But if you're not at an edge, you will have 8 8-connected pixels (9 total) and 4 four-connected pixels (5 total with the center one). So I'm still not clear what region you want. Perhaps if you described why you wanted it, or what you would do as the next step if you had this information.

By the way, round is a built in function and you should not use it as a variable name. And you can get the local mean image with a call to conv() or imfilter() instead of that double for loop.

[amen45]

The value in each pixel represents the relationship between a wave of frequency x_index and a wave of frequency y_index. I'm trying to find the range of wave pairs that have the strongest relationship, since wave pairs that are close together are likely produced by the same phenomenon. I'll then use this information to compare the distribution of the wave-pair bands with the strongest relationship across all my subjects.

I don't have any image processing experience-- after looking into what you described, I think part of the reason that this is confusing is because I don't have a predetermined shape/orientation that I'm looking for. The region with the highest values could be 4x1 or 2x2. I just want the number of wave-wave pairs (the pixels) in the region I pull out to be consistent. Does that make sense?

[amen45]

Essentially, my ideal function will take 2d array and some integer k:

#1 Find all "scrolling" array sizes possible using k that don't exceed the size of the 2d array (for example if array has dimensions 8x10 and k is 12, "scrolling arrays" could be 2x6,6x2,3x4 and 4x3)

#2 find the region of size k in whatever dimensions that has the highest values, and pull out indices and mean

[Chaya N]

Are you perhaps trying to track a particular region and its changes through a set of images? A quick and dirty way might be to simply define a square mask (>k) such that your elements of interest are given a higher 'weight'. This is only so you could define a clear threshold or shape for the connected region you are interested in through a little numerical manipulation.

Consider this:

x = [0.115754943948525  0.254835774407905  0.418797235529874
0.133596512865232  0.267704112017176  0.427967388090852
0.148849791641271  0.274078404255313  0.425934632215066
0.166527249110295  0.977300000000000  0.998700000000000
0.196980352387009  0.923140000000000  0.987890000000000
0.224172402663761  0.323278947773164  0.431610265045853
0.283912660670702  0.380237101362344  0.466341084111908];
h = [1 1 1   % A 3x3 mask that weights the center value greater than its surroundings
    1 100 1
    1 1 1];
y = conv2(x,h,'same'); % Filter your original array 
z = y.*x; % This is only to exaggerate the differences between the array values

From z, you would just pick out the k highest valued elements. This is just a quick illustrative example that I have not tried on any other data (it may mostly work or not).

[Chaya N]

To add to the above, I would suggest a square mask (or scrolling array, if you'd like) with odd number of elements so you would have something of a 'center' value. This is a common tactic used in filtering operations. Also, the mask size and shape would influence the data, so choose as small a square size as possible so you would not end up manipulating your data over the necessary amount. In my example, I have chosen a 3x3 mask because it is the smallest odd square that I could use that is greater than k.

[Image Analyst]

To get the mean of the interactions (pixels) in a horizontal 1-by-2 sliding window, do this:

[rows, columns] = size(im)
for col = 1 : columns
  windowWidth = col;
  kernel = ones(1, windowWidth) / windowWidth;
  meanImage = conv2(im, kernel, 'same');
  % Now do something with meanImage.....
end

where m is your interaction matrix.

To do the same but in a vertical 2-by-1 sliding window, do this:

[rows, columns] = size(im)
for row = 1 : rows
  windowHeight = row;
  kernel = ones(windowHeight, 1) / windowHeight;
  meanImage = conv2(im, kernel, 'same');
  % Now do something with meanImage.....
end

If you want, you could use a double for loop over both rows and columns to create all possible kernel shapes:

[rows, columns] = size(im)
for col = 1 : columns
  for row = 1 : rows
    kernel = ones(row, col)/(row * col);
    meanImage = conv2(im, kernel, 'same');
    % Now do something with meanImage.....
  end
end

Not sure what you want to do with the mean image once you have it though.