extracting mean pixel value from an image

I would like to extract the mean pixel value on the image.

I would like to do it automatic by locating the brightess point next to the green rectangular point coordinates, cropping the new coordinates and finding its mean pixel vaules

Answer 1

[DGM]

This might be a start, but bear in mind how fragile this will be if the colors change or the swatch becomes rotated or something.

A = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1004350/Picture1.png');

Ahsv = rgb2hsv(A);
[H S V] = imsplit(Ahsv);
mask = S>0.9 & H>0.15; % create mask
mask = bwareaopen(mask,100); % despeckle

% use a large strel to bridge dotted lines
% this assumes that lines are roughly grid-aligned
stlen = 200;
mask = imopen(~mask,ones(stlen));

% erode to avoid including any green bits on edges
mask = imerode(mask,ones(10));

% this assumes that the largest blob is the ROI
mask = bwareafilt(mask,1);

imshow(mask)

% build an Mx3 color table from the pixels in the ROI
roipix = zeros(nnz(mask),3,class(A));
for c = 1:3 % assumes A is RGB
    thischan = A(:,:,c);
    roipix(:,c) = thischan(mask);
end

% find the mean color
meancolor = mean(roipix,1)

meancolor = 1×3

220.5722 176.5078 21.8325

% for visualization of the nonsquare ROI
% fill non-ROI with black
rgbpict = A;
rgbpict(repmat(~mask,[1 1 3])) = 0;

% pad the modified image with a color swatch
sz = size(A);
meanswatch = repmat(permute(meancolor,[1 3 2]),[50 sz(2) 1]);
meanswatch = uint8(meanswatch);
outpict = [rgbpict; meanswatch];

imshow(outpict)

Note that's just the simple arithmetic mean. It might not be what you want, considering how it will be influenced by the light bleed and the vignetting. Alternatively, you might try using median() instead of mean.

See this answer for other options. Based on that answer,

% get some different image stats from the ROI
cc = imstats(ctflop(roipix),'mean','median','mode','modecolor','modefuzzy','moderange','nmost',10);

% use those stats to construct a swatch chart for single-output stats
labels = {'mean','median','mode','modecolor','modefuzzy'};
sz = imsize(rgbpict,2);
ntiles = (numel(labels));
tilesz = [round(sz(1)/ntiles) 100];
block1 = zeros([tilesz 3 ntiles],'uint8');
for k = 1:ntiles
    thistile = colorpict([tilesz 3],cc(k,:),'uint8'); % colored swatch
    thislabel = im2uint8(1-textim(labels{k},'ibm-iso-16x9')); % text label image
    thistile = im2uint8(imstacker({thislabel thistile},'padding',1)); % match geometry
    block1(:,:,:,k) = mergedown(thistile,1,'linearburn'); % blend label and swatch
end
block1 = imtile(block1,[ntiles 1]); % vertically arrange tiles
block1 = imresize(block1,[sz(1) tilesz(2)]); % make sure it's the right size

% create another chart for moderange's multiple outputs
ntiles = (size(cc,1)-ntiles);
tilesz = [round(sz(1)/ntiles) 100];
block2 = zeros([tilesz 3 ntiles],'uint8');
for k = 1:ntiles
    thistile = colorpict([tilesz 3],cc(k+4,:),'uint8'); % colored swatch
    thislabel = im2uint8(1-textim(num2str(k),'ibm-iso-16x9')); % text label image
    thistile = im2uint8(imstacker({thislabel thistile},'padding',1)); % match geometry
    block2(:,:,:,k) = mergedown(thistile,1,'linearburn'); % blend label and swatch
end
block2 = imtile(block2,[ntiles 1]); % vertically arrange tiles
block2 = imresize(block2,[sz(1) tilesz(2)]); % make sure it's the right size

% show the combined images
imshow([rgbpict block1 block2])

Make of that what you will. See the linked answer for more details on these options.

[DGM]

Maybe this is more like the ROI that you're after

A = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1004350/Picture1.png');

Ahsv = rgb2hsv(A);
[H S V] = imsplit(Ahsv);
roimask = S>0.9 & H>0.15; % create mask
roimask = bwareaopen(roimask,100); % despeckle

% use a large strel to bridge dotted lines
% this assumes that lines are roughly grid-aligned
stlen = 200;
roimask = imopen(~roimask,ones(stlen));

% erode to avoid including any green bits on edges
roimask = imerode(roimask,ones(10));

% this assumes that the largest blob is the ROI
roimask = bwareafilt(roimask,1);

imshow(roimask)

% get mask of all light peaks
peakmask = mean(A,3)>242; % HSI intensity
peakmask = peakmask & roimask;

imshow(peakmask)

% get centroids
S = regionprops(peakmask,'centroid');
Cpeaks = vertcat(S.Centroid);
S = regionprops(roimask,'centroid');
Croi = vertcat(S.Centroid);

% find four peaks most distant from ROI center
D = sqrt(sum((Cpeaks-Croi).^2,2));
[~,idx] = maxk(D,4);
Cpeaks = Cpeaks(idx,:);

% sort by angle
th = atan2d(Cpeaks(:,2)-Croi(1,2),Cpeaks(:,1)-Croi(1,1));
[~,idx] = sort(th,'ascend');
Cpeaks = Cpeaks(idx,:);

% create polygonal ROI object
ROI = images.roi.Polygon(gca);
ROI.Position = Cpeaks;

% create mask from ROI object
mask = createMask(ROI);

imshow(mask)

% build an Mx3 color table from the pixels in the ROI
% this is all the pixels in the ROI, so you can process it however you want
roipix = zeros(nnz(mask),3,class(A));
for c = 1:3 % assumes A is RGB
    thischan = A(:,:,c);
    roipix(:,c) = thischan(mask);
end

% find the mean color
meancolor = mean(roipix,1)

meancolor = 1×3

221.2379 175.8644 20.7190

% find median color
mediancolor = median(roipix,1)

mediancolor = 1×3

222 170 3

% for visualization of the nonsquare ROI
% fill non-ROI with black
rgbpict = A;
rgbpict(repmat(~mask,[1 1 3])) = 0;

% to visualize, pad the original image with a color swatch
sz = size(A);
meanswatch = repmat(permute(meancolor,[1 3 2]),[50 sz(2) 1]);
meanswatch = uint8(meanswatch);
outpict = [rgbpict; meanswatch];

imshow(outpict)

I'm not really sure if you're trying to do perspective correction on the polygonal ROI so as to make it into a new rectangular image. If all you need to do is color analysis, it shouldn't be necessary.

[OLUWAFEMI AKINMOLAYAN]

Also, when you mention mean color, are you referring to the rgb value?

[DGM]

Using picture4.jpg, I changed the threshold value on this line:

peakmask = mean(A,3)>210; % HSI intensity

and got this composite:

... so it should work.

None of the images in these demos have lines drawn on them. The one with blue lines is showing the temporary ROI object before it's used to create the polygon mask. It's not part of the image.

The variable called roipix contains all of the RGB values in the polygonal ROI, formatted as an Mx3 matrix. If you want to export it to a spreadsheet, you should be able to use something like writematrix(). As a lot of the functionality for reading/writing to XLS/XLSX files doesn't work in my environment, I'm not the person to ask for a direct example.

When I mention mean color, yes. I'm talking about the mean RGB color. If you're after a grayscale representation of the ROI, you'll have to specify how exactly you want to reduce RGB to grayscale.

% just a smaller image for sake of comparison
smallpict = imresize(rgbpict,0.25);

% these are all different ways of reducing RGB to gray
Y601 = rgb2gray(smallpict);
Y709 = imapplymatrix([0.2126 0.7152 0.0722],double(smallpict),'uint8');
Ls = im2uint8(rgb2lightness(smallpict)/100);
V = max(smallpict,[],3);
I = mean(smallpict,3);
L = uint8((double(max(smallpict,[],3)) + double(min(smallpict,[],3)))/2);

% concatenate and display
allgrays = [Y601 Y709; Ls V; I L];
imshow(allgrays)

[OLUWAFEMI AKINMOLAYAN]

Wow! Thanks alot. I have lot of questions. I will be re-using the code and probably be modifying it for multiple use. So, I will need to understand all the synthax and why they are used in the code.

I will be glad if i can get detail explaination of each line in the code.

However, if it is too much, i can ask the specifics.

1. Why do we need to convert to HSV before we mask?

2. How do you determine the mask value s> 0.9 & H>0.15?

3. Can us explain what this does "peakmask = mean(A,3)>242;"

4. Why is 242 chosen for the image?

5. Can you explain a little more about this part.

"% find four peaks most distant from ROI center

D = sqrt(sum((Cpeaks-Croi).^2,2));

[~,idx] = maxk(D,4);

Cpeaks = Cpeaks(idx,:);

% sort by angle

th = atan2d(Cpeaks(:,2)-Croi(1,2),Cpeaks(:,1)-Croi(1,1));

[~,idx] = sort(th,'ascend');

Cpeaks = Cpeaks(idx,:);"

I understood what was done, but dont understand most of the lines and synthax

6. I dont really undertsand the use of the synthax in the other parts of the code.

7. Also which of the grayscale conversion converts to 8bits?

8. how can i include the pixel coordinates associated with the rgb values in the roipix variable?

9. And lastly, a general question. Is pixel value equal to the average of the rgb values?

[DGM]

I'm going to try to walk through the prior example and show how things progress. It'll be a lot of images, but maybe it'll clear up a few uncertain things to see the intermediate steps.

% the input image
A = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1004350/Picture1.png');

% convert to HSV so that we can use color information
% to discriminate between green marks and surroundings
Ahsv = rgb2hsv(A);
[H S V] = imsplit(Ahsv);

It should be easy to see some advantage to being able to use color information to isolate the green marks from the predominantly yellow background. If this were done strictly in grayscale, it would be very difficult to reliably find them.

% show both H and S for visualization:
imshow(H)

imshow(S)

Notice that the green marks have a hue which is only slightly higher than their surroundings. Also notice that colors near black and white may have any number of hues. By excluding low-saturation regions, we can avoid picking up spurious hues in bright and dark regions -- all while using relatively open threshold constraints.

roimask = S>0.9 & H>0.15; % create mask
roimask = bwareaopen(roimask,100); % despeckle
imshow(roimask)

% use a large strel to bridge dotted lines
% this assumes that lines are roughly grid-aligned
stlen = 200;
roimask = imopen(~roimask,ones(stlen));

% erode to avoid including any green bits on edges
roimask = imerode(roimask,ones(10));
imshow(roimask)

% this assumes that the largest blob is the ROI
roimask = bwareafilt(roimask,1);
imshow(roimask)

Now we have a rough mask that should contain all the bright spots of interest. Need to reduce to a binary mask. The value 242 (or later revised to 210) was simply chosen by using a datatip and looking at the mean image. If the goal here is to find the peaks based on their brightness, why not use V as calculated earlier in the HSV conversion? Note from the prior example what the V image looks like. It has very poor contrast because value (V) only represents the maximum component of RGB. Out of all the available "brightness" metrics, V would be one of the worst for this.

In the absence of knowledge that suggests that perceived brightness is of importance, I chose to consider a uniform weighting of the RGB components (the simple arithmetic mean, aka "I" in HSI), as opposed to a nonuniform weighting as with luma (Y).

It's worth noting that while I couldn't use V from prior calculations used for finding H and S, I could have reused prior calculations if I had used HSI instead of HSV. The reason I didn't is simply that MATLAB doesn't have HSI tools. It would complicate the example to include them.

% get mask of all light peaks
peakmask = mean(A,3)>210; % HSI intensity
peakmask = peakmask & roimask;
imshow(peakmask) % the forum editor displays this down below with the ROI overlay

% get centroids
% Cpeaks is a list of [x y] positions of the centers 
% of each white blob in peakmask
% Croi is the [x y] position of the center of the blob in roimask
S = regionprops(peakmask,'centroid');
Cpeaks = vertcat(S.Centroid);
S = regionprops(roimask,'centroid');
Croi = vertcat(S.Centroid);

% find four peaks most distant from ROI center
% calculate the euclidean distance between Croi and each peak
D = sqrt(sum((Cpeaks-Croi).^2,2)); 
[~,idx] = maxk(D,4); % find the indices of the 4 largest distances
Cpeaks = Cpeaks(idx,:); % keep only those center positions

At this point, we have the coordinates of the four corner blobs in peakmask, but they're not necessarily sorted in any useful way. Because the ROI polygon needs to have ordered vertices in order to not self-intersect, the list of coordinates needs to be sorted. The difference Cpeaks(:,2)-Croi(1,2) is the y-offset between the center of roimask and each peak blob; similarly, Cpeaks(:,1)-Croi(1,1) is the x-offset. Using atan2d(), we have a list of 4 angles between -180 and 180. We don't need the sorted angles for anything, but we do need to know the order of the indices associated with the sorted angles. Like the above example with maxk(), we can discard the first output and just keep the second output (the indices).

The ROI tools are generally intended to be used interactively with a mouse (e.g. drawPolygon(gca)), but they can also be used programmatically. In this case, an empty object is created and then the centroids are used as a vertex list to create a polygonal region that can then be converted into a new binary mask.

% sort by angle
th = atan2d(Cpeaks(:,2)-Croi(1,2),Cpeaks(:,1)-Croi(1,1));
[~,idx] = sort(th,'ascend'); % get indices
Cpeaks = Cpeaks(idx,:); % use indices to reorder centroid (vertex) list

% create polygonal ROI object
ROI = images.roi.Polygon(gca);
ROI.Position = Cpeaks;

% create mask from ROI object
mask = createMask(ROI);
imshow(mask)

So now we don't need peakmask or roimask anymore. We just have a clean polygonal mask to extract a region from the image. Step through each channel (a 2D array) and extract the elements associated with the new mask.

% build an Mx3 color table from the pixels in the ROI
% this is all the pixels in the ROI, so you can process it however you want
roipix = zeros(nnz(mask),3,class(A));
for c = 1:3 % assumes A is RGB
    thischan = A(:,:,c);
    roipix(:,c) = thischan(mask);
end

% just dump the contents of roipix
roipix % [R G B]

roipix = 323598×3

255 255 255 255 255 254 255 255 245 253 255 218 254 254 187 254 254 157 254 255 133 254 254 118 255 254 105 255 252 93

You bring up a good point. What about the coordinates of the elements in roipix? We can get those like so:

% get a vector of linear indices for all nonzero elements of the mask
idx = find(mask); 
% convert from linear indices to row,column subscripts
[row col] = ind2sub(size(mask),idx);

Since we're going to write this to a spreadsheet, they'll probably be concatenated with the RGB data:

outdata = [row col roipix] % [row column R G B]

outdata = 323598×5

82 200 255 255 255 83 200 255 255 254 84 200 255 255 245 85 200 253 255 218 86 200 254 254 187 87 200 254 254 157 88 200 254 255 133 89 200 254 254 118 90 200 255 254 105 91 200 255 252 93

[DGM]

To wrap up a few questions:

6. I dont really undertsand the use of the synthax in the other parts of the code.

That's a bit too nonspecific to answer directly. I chose to be extra-verbose in the prior comment to hopefully address some points of confusion. If there are other issues, you'll have to be more specific.

7. Also which of the grayscale conversion converts to 8bits?

In the examples I gave, all the outputs are 8 bit (uint8). Bear in mind that both A and rgbpict are RGB uint8 images. Some of the given functions will preserve the class and scaling of their input, so the outputs are implicitly also uint8. Others don't, or might require intermediate conversions due to the limitations of working with integers:

% this is a uint8 image (RGB)
smallpict = imresize(A,0.25);

% the output of rgb2gray() will inherit the class of the input (convenient!)
Y601 = rgb2gray(smallpict);

% imapplymatrix()'s class handling is a confusing disaster
% it would take a paragraph to explain, but the output is uint8 in this case
Y709 = imapplymatrix([0.2126 0.7152 0.0722],double(smallpict),'uint8');

% the output of rgb2lightness() is 'double', scaled [0 100]
% so it needs to be rescaled and cast to uint8 [0 255]
Ls = im2uint8(rgb2lightness(smallpict)/100);

% these are simple and inherit the input class
V = max(smallpict,[],3);
I = mean(smallpict,3);

% but this approach required a wider numeric class to hold intermediate 
% calculations (the sum).  This could be simplified, but I don't
% think that's the focus here.  The output is still uint8.
L = uint8((double(max(smallpict,[],3)) + double(min(smallpict,[],3)))/2);

% the reason that I made sure that they were all uint8
% was so that I could concatenate them all into one image
allgrays = [Y601 Y709; Ls V; I L]; 

9. And lastly, a general question. Is pixel value equal to the average of the rgb values?

It depends what you mean by "pixel" and "value". Maybe I'm being pedantic, but I understand how they get used loosely and that confusion results when context gets lost or misunderstood.

Let's consider an image called mypurple:

mypurple(1,1,:) % look at one pixel

1×1×3 uint8 array

ans(:,:,1) =

100

ans(:,:,2) =

20

ans(:,:,3) =

255

So the result is a 1x1x3 vector representing the color content at a given spatial location. When I say "pixel", I usually mean something like this. A "pixel" in a grayscale image would be 1x1 (scalar). A pixel in RGB would be 1x1x3. in RGBA, it would be 1x1x4. It's the smallest spatial subdivision of the image.

In other contexts, people might use the word "pixel" to also describe the elements of the 1x1x3 vector that I called a "pixel". It might depend on the context, as images aren't always represented as MxNxC rectangular arrays. If an image is vectorized for transport or storage in a file, the spatial relationships aren't really intuitive anymore.

Long story short, I'm usually talking about "whole" pixels when I use the word. If I want to talk about the elements of mypurple(1,1,:), I might talk about the components or elements of a pixel.

But those are just numbers representing a physical quantity. You could also say that each element of that vector contains numeric values. That's a perfectly valid description. The confusion might arise when talking about something like those grayscale images. In that case, "value" might also refer to the brightness metric used by the HSV color model. Similarly, "intensity" and "lightness" have multiple specific technical meanings which aren't interchangeable. HSI "intensity" isn't equivalent to actual luminous intensity (or other measures of physical intensity). Similarly, HSL lightness isn't equivalent to L* lightness as used in CIELAB/CIELUV. I think there just aren't enough synonyms for "brightness" in the English language.

So back to the question as worded: Is pixel value equal to the average of the rgb values?

Considering an RGB image, and my interpretations of these words, the value (V in HSV) of a "whole" pixel is the maximum of the components of the pixel. For the example with mypurple, the HSV value of [100 20 255] is 255. As mentioned before, the simple average would be what is called I (intensity) in HSI. The HSI intensity would be mean([100 20 255]), which is 125.

I honestly can't tell if this will make things more or less confusing. If there's a specific thing I said that's ambiguous, feel free to point it out and I can clarify what I mean in-context.

[OLUWAFEMI AKINMOLAYAN]

Thanks a lot. I now understood some of it now. I might have to digest the comments more than once to fully understand all. Wish i can pay you for the great work you have done.

I have few last clarifications/questions.

1. You mentioned something about it. However, if i am to use another image, what is the metric for choosing the threshold for the peakmask. You chose different values (210, 242) for the different image. In order word, How can i chose the value myself for diifferent images. Is it by guessing or i have to do some intuitive calculation?
2. If i am to populate the code with multiple images, how can i do that. For example, if i have like 50 images and i want the code to extract the roi the same way and store the coordinates and rgb values of all the images' roi, how can i do that?
3. And last can i show the roipix only by using imshow(roipix)?

[OLUWAFEMI AKINMOLAYAN]

One other thing. The outpict that showed the visualisation of the roi. Is there a way i can crop out the ROI to form a new standalone image. If i can crop out the ROI extracted from the non-ROI, I could just easily extract the rgb values of the new cropped image without going through the roipix code

[OLUWAFEMI AKINMOLAYAN]

The reason i asked the question "Is pixel value equal to the average of the rgb values?" is that when onw want to find the histogram of an image, the x-axis show values from 0 - 256. Are those number pixel values or average of the rgb values?

[DGM]

For picking the threshold for peakmask, you don't need to be too terribly specific. All it needs to be is low enough that all the peaks get picked up, and high enough that the peaks are well-defined and isolated from each other. For the two example images given, 210 will work for both. Consider what would happen if I went down to 150:

So you might need to look at some of your images and tweak the value as necessary to find a threshold that works for all of them. As far as how I picked them, I just did this:

imshow(mean(double(A),3),[0 255])

and then used a tooltip to inspect the corner peaks

If you need to process multiple files (multple image files being read and multiple spreadsheets being written), you'll end up needing to set up some sort of loop to handle reading each image, processing it, and then writing each spreadsheet.

https://www.mathworks.com/matlabcentral/answers/57446-faq-how-can-i-process-a-sequence-of-files

https://www.mathworks.com/help/images/process-image-sequences.html

Note that some examples may show reading all the images at once. If you can avoid it, it's probably better in this case to read them and process them one at a time. That way you don't have a bunch of stuff occupying memory.

It may be worthwhile to encapsulate the core of the main code in one or more functions, though if you're more comfortable working with a long script, that can also work.

You may be able to get the coordinates as described and write them to an XLSX file like so:

writematrix(outdata,'myfilename.xlsx')

or considering that you'll be doing this in a loop:

for k = 1:number_of_files
    % read file
    % ...

    % process file
    % ...
    
    % write file
    thisfname = sprintf('myfilename_%04d.xlsx',k);
    writematrix(outdata,thisfname)
end

... though you may need to come up with some convention for output file naming.

You can't practically view roipix directly with imshow(), since it's essentially one long stripe. It would show up as a nearly-invisible line. If you wanted to visualize what exactly the calculated mask is selecting, you can do this

% for visualization of the nonsquare ROI
% fill non-ROI with black
rgbpict = A;
rgbpict(repmat(~mask,[1 1 3])) = 0;

imshow(rgbpict)

... though bear in mind that this is really just for visualization. The problem with trying to use this directly is that it's full of black pixels that are needed for padding (since the ROI is not rectangular). The nice thing about roipix is that it doesn't matter what the ROI is shaped like. It's a list of all pixels in the mask region, taken columnwise. It can be of any length -- including lengths which are prime numbers (meaning that the result can't be reshaped into a rectangular image).

If you wanted to crop the above visualization, you could, but the fact that the ROI is nonrectangular means that there will necessarily be padding. Arrays must be rectangular.

Since the polygon is a quadrilateral, you could transform it into a rectangle if you wanted to. At that point, it could fit in a rectangular array without padding.

% these are the coordinates of the box corners
boxm = Cpeaks;

% assert that this is where they're supposed to be
boxf = [Cpeaks([1 2 2 1],1) Cpeaks([1 1 3 3],2)];

% transform the image
TF = fitgeotrans(boxm,boxf,'projective');
outview = imref2d(size(A));
B = imwarp(A,TF,'fillvalues',255,'outputview',outview);

imshow(B)

% crop out the transformed region based on output coordinates
croprect = [min(boxf,[],1) range(boxf,1)];
B = imcrop(B,croprect);

imshow(B)

That works, but it adds a lot of extra complexity. It might look nice, but it doesn't really change the color content. If what you need is to do is to store it in a tabular form, this doesn't really help either. If you don't really care about storing it in tabular form (a spreadsheet), then maybe this would be desirable if you wanted to store the output as an image instead?

[OLUWAFEMI AKINMOLAYAN]

Thanks once a again for the detailed explanation.

From you last paragraph, do you mean one can not store the histogram or pixel information of the ''B'' image in a spreadsheet?

Also, in one of the image labelled above, you have [x, y], index, and [rgb]. the index was equal to 245.3 which show the intensities of light in that spot. I always think that value is the pixel/greyscale value measure the level of dark to bright using 0 - 255 index. I just need you to clarify my remaining doubts on this.

[OLUWAFEMI AKINMOLAYAN]

If i could get the the numbers you tagged as index from the ROI with their coordinate, then i will be glad. This is because, I believe the index shows some kind of brightness level from the image which is integral to my analysis.

I initially thought average of rgb values will give me the index, but after going through your explanations, i realised i was wrong

[DGM]

You can certainly do processing on the transformed/cropped image. I was just saying that if you want the pixel data as a Mx3 list in a spreadsheet, then you'd still have to reshape it. It also doesn't really matter what shape the region is if your goal is to do some sort of color analysis, so transforming it into a perfect rectangle may be unnecessary toward that goal.

So far I haven't said anything about histogramming anything (I don't think I did), but if you want to get the histogram, you can get a histogram from either roipix or from the transformed image above.

The bit with the datatip can be confusing. Bear in mind that imagesc() and imshow() display single-channel images using a colormap (regardless of numeric class), unlike RGB images which are shown in truecolor and don't use a colormap. The "Index" value it's showing is the actual value of the pixel. The RGB tuple shown below is the associated color at that position in the colormap (which is unit-scale floating point by convention).

So there's two things to note. First is the fact that the 245.3 isn't an integer when you might expect it to be. That's just because I took the mean() of the image, so the values aren't rounded. For sake of inspection, it doesn't really matter. Second is that the color tuple isn't generally useful information. Consider if I had done this:

imshow(uint8(mean(A,3)),[0 255])
colormap(parula)

Note that the value is an integer in the range [0 255], and the RGB tuple is no longer a neutral gray, even though the image obviously is gray. So yes, the "index" field in the tooltip refers to the image values when using a grayscale image. Datatips behave differently for RGB images. The scale of the values depends on the image class. Consider what it would look like if I did this:

imshow(A)

You can also use impixelinfo() to inspect the image values. I can't really show a screenshot because it's always interacting with the mouse, but if you do

imshow(A)
impixelinfo

you can just mouse over the image and the coordinates and pixel values will be displayed at the bottom of the figure.

[OLUWAFEMI AKINMOLAYAN]

Thanks for the clarification. I still have some confusion. in the last image, does this mean the average of the rgb value is equal to the index (pixel value)?

if its not the average, then how do I store the index (pixel value) ranging from 0 - 255 for the output image "B"?

You used variable "B" for the cropped ROI and showed how i can store the rgb value and the coordinate. What if i just want the index (pixel value).

[DGM]

The values in the three datatips aren't really comparable since I didn't pick the exact same spot each time. The takeaway is that the things that are displayed in the datatip depend on how many channels the image has. If it has 3 channels, it just gives the RGB values at that location. If it has 1 channel, it gives the pixel value at that location (called "index") and the associated colormap tuple. The fact that it's labeled as an index is misleading, since it's not necessarily an indexed image, it's not necessarily an index into the colormap that's in use, and it's obviously not the linear index of the pixel (the position in the array). It's just whatever the pixel value is.

In the one example I gave, the pixel value happens to be the mean because i'm explicitly calculating the mean and displaying it. I was using the mean because it was sufficient for the task of finding the peaks without complications. Sometimes it's beneficial to use color information for identifying parts of an image; other times color information is just extraneous information. Color helped isolating the green marks, but it would only complicate finding prominent features like the bright peaks.

I'm not sure what you mean by index. The output image B is class 'uint8'. Its an integer image with values in the range of [0 255]. It's not an indexed image. Storing the "RGB values" is the same thing as storing the "pixel values" because each pixel contains three R,G, and B components.

[OLUWAFEMI AKINMOLAYAN]

Still confused somehow, So what is the x- axis values in a image histogram showing? if its not the average of rgb values, what values are they?

[OLUWAFEMI AKINMOLAYAN]

I read through your comments again. So i guess its a pixel in a grayscale image? Can you put me through how i can output the pixel values after converting the image to a grayscale image

[OLUWAFEMI AKINMOLAYAN]

While i was working around the codes, the number of pixels for roipix (324569) is not equal to the cropped image "B" (315262). Do you know why this happened?

Also, while working with the cropped image "B", I could not retain the orignally pixel coordinates from my image "A". The pixel coordinates are refreshed. How do i retain the coordinates?

[OLUWAFEMI AKINMOLAYAN]

ARE YOU THERE?

[DGM]

Sorry about that. Your last comments got buried in my notifications. Sometimes that happens, depending on the order that things show up in my feed.

It depends how you calculated the histogram. Consider the image:

A = imread('10A_22.jpg');
imhist(A)

That's the histogram of all elements of A, regardless of which channel they belong to.

subplot(3,1,1); imhist(A(:,:,1)); % red
subplot(3,1,2); imhist(A(:,:,2)); % green
subplot(3,1,3); imhist(A(:,:,3)); % blue

The difference between roipix and the cropped/reshaped subimage B is fundamental. They are not the same image anymore. They do not contain the same number of pixels and each pixel in roipix has been interpolated. Its relative location within the image has been altered due to the projective transformation and removal of a large offset due to cropping.

If you want the coordinates of the transformed image, it depends whether you want the coordinates of the pixels of B with respect to the extent of B, or whether you want the coordinates of the pixels of B with respect to the original untransformed image.

If you want the coordinates of pixels of B with respect to the extent of B:

sz = size(B);
[xx yy] = meshgrid(1:sz(2),1:sz(1)); % these are the coordinates

If you want the coordinates of pixels of B with respect to the original image:

sz = size(A);
[xx yy] = meshgrid(1:sz(2),1:sz(1));

% these are the coordinates of the box corners
boxm = Cpeaks;

% assert that this is where they're supposed to be
boxf = [Cpeaks([1 2 2 1],1) Cpeaks([1 1 3 3],2)];

% transform the image
TF = fitgeotrans(boxm,boxf,'projective');
outview = imref2d(size(A));
xx = imwarp(xx,TF,'fillvalues',0,'outputview',outview);
yy = imwarp(yy,TF,'fillvalues',0,'outputview',outview);

% crop out the transformed region based on output coordinates
croprect = [min(boxf,[],1) range(boxf,1)];
xx = imcrop(xx,croprect);
yy = imcrop(yy,croprect);

As before, you can format these rectangular arrays into columnar format:

% reshape B and desired coordinates into a columnar format
outdata = [xxtr(:) yytr(:) double(reshape(B,[],3))];

Again, if all you need is color information, the transformation is simply adding needless complexity.

As for the new images, some adjustments need to be made. 40A_12.jpg seems to process just fine with the prior settings. The other two images require a larger strel and a different masking approach, since there isn't sufficient chroma information to use the difference between green/yellow. They do have very distinct V differences, so that can be used instead.

clf

% the input image
%A = imread('Picture1.png'); % worked fine with higher Ithresh
%A = imread('picture4.jpg'); % needed Ithresh reduced to 210
%A = imread('40A_12.jpg'); % works fine with stlen=200, Ithresh=210
%A = imread('25A_13.jpg'); % use V for mask, stlen=600, Ithresh=210
A = imread('10A_22.jpg'); % use V for mask, stlen=600, Ithresh=210

% convert to HSV so that we can use color information
% to discriminate between green marks and surroundings
Ahsv = rgb2hsv(A);
[H S V] = imsplit(Ahsv);
%roimask = S>0.9 & H>0.15; % create mask
roimask = V<0.3; % mask for 25A_13.jpg, 10A_22.jpg
roimask = bwareaopen(roimask,100); % despeckle

% use a large strel to bridge dotted lines
% this assumes that lines are roughly grid-aligned
stlen = 600; % this may need adjusted
roimask = imopen(~roimask,ones(stlen));

% erode to avoid including any green bits on edges
roimask = imerode(roimask,ones(10));

% this assumes that the largest blob is the ROI
roimask = bwareafilt(roimask,1);

imshow(roimask)

% get mask of all light peaks
peakmask = mean(A,3)>210; % HSI intensity
peakmask = peakmask & roimask;

imshow(peakmask)

This is probably turning into a good example of how it's often beneficial to spend time designing the physical part of the photographing/staging process in such a way that the software part of the process (the image segmentation) is more consistent and reliable. If there are still a bunch of photos that need to be taken, that might be something to think about.

[OLUWAFEMI AKINMOLAYAN]

Thanks alot for the response. You are a genius.

I have the following questions.

I am still not clear on the H, S and V threshold usage. perhaps. Is there a reference i could read to understand it more?

What is the effect of the strel value?

And lastly, is there a way i can come the codes to work for all case. Perhaps using conditional statements will work?

[DGM]

I don't really have a good reference for applying basic HSV thresholding. It's a combination of familiarity with the color model itself and the usage of basic tools such as datatips as described. You could also try to use the color thresholder app that comes with Image Processing toolbox, though I find it cumbersome.

Ultimately it's a matter of finding some particular way to describe the colors in an image such that the features you want to isolate have a clearly unique description. Sometimes you can get away with separating things by brightness alone; other times it's better to use hue or chroma/saturation. Often times there aren't any simple ways to do segmentation via global color-based thresholding alone.

You can think of the strel (structuring element) like a filter kernel. It's a small 2D logical array which defines a neighborhood around a given pixel during a filtering process. In this case, it's just a simple square defined by ones(stlen). Consider the case for 10A_22.jpg. After thresholding, we have this image:

If we try to filter it using a 200x200 strel, the dots can't be connected, since they're farther apart than 200px.

Try a 400x400 strel:

Try a 600x600 strel:

There are probably plenty of more complicated ways to get the dots connected, but in this method, the strel needs to be big enough that it bridges the dots consistently. As the required strel becomes larger, the process becomes less robust against image rotation.

Coming up with code that works for all cases requires that you know what all cases are -- or at least that the manner in which the cases vary has some known constraints. As I mentioned, if you have the opportunity to design the process such that the photos you're taking can be easily and consistently processed, that would be ideal. If the spacing and apparent color of important features like the green dots can vary so widely, you are going to have a difficult time coming up with a single process that can account for all cases.

Some thoughts:

If the marker dots were lines (even if there are breaks due to the texture), it would help with segmentation. Using a straightedge might help.

Using a stencil (or maybe tape) instead of marker might be an alternative way to physically mark out the ROI when taking the picture.

You might be able to put in some tests in the segmentation process to check that the intermediate results are sane. At the very least, this could alert you to make some manual adjustment. In a more elaborate setup, it might allow you to automatically try a different method to perform the segmentation.

Toward that last point, consider the example of the case where roimask was created with a 400x400 strel. In that scenario, this would be what peakmask looks like:

If you expect this to be 81 dots in a 9x9 grid, you might check the number of elements in the struct S returned by regionprops. In the malformed case shown, you have 88 dots instead of 81. You might use that information as some form of process indicator.

Are there more elaborate ways of finding the pattern using machine learning? There probably are, but I am completely unfamiliar with the subject.

[OLUWAFEMI AKINMOLAYAN]

Hello DGM,

I was working on modifying the code for multiple image.

What if we assume that the grid line are not there.

Can we use the method used in locating the four point cordinate of the ROI rectangle?

Can we do that without using the green lines to mask? Just thinking of other options too.

Thanks

[DGM]

... maybe. It would depend on how you consistent the framing/scale is between photos. Let's consider the image above. If you assume that the ROI is defined by a 9x9 grid of dots in the center ...

I have an idea. It'll probably be flaky, but give me a bit to throw something together...

[DGM]

This might not actually be too bad. It works for all the given images so far without any consideration of the marker dots. You can tweak the angle tolerance quite a bit. I modified one of the images to demonstrate that it still works even if there's quite a bit of distortion.

% the input image
A = imread('10A_22_skew.jpg');

% get mask of all light peaks
peakmask = mean(A,3)>210; % HSI intensity
peakmask = imclearborder(peakmask);
peakmask = bwareaopen(peakmask,100);

imshow(peakmask); hold on

% get centroids
S = regionprops(peakmask,'centroid');
Cpeaks = vertcat(S.Centroid);
Cimg = fliplr(size(peakmask))/2;

% find the peak closest to the image center
D = sqrt(sum((Cpeaks-Cimg).^2,2));
[~,idx] = min(D);
Ccenter = Cpeaks(idx,:);

% mark the center peak
plot(Ccenter(1),Ccenter(2),'r*','linewidth',2)

% find the 8 peaks closest to this peak
D = sqrt(sum((Cpeaks-Ccenter).^2,2));
[~,idx] = mink(D,9);
Cnh = Cpeaks(idx(2:end),:); % ignore self-distance

% find points which are oriented roughly 45 
thtol = 16; % allowable angle deviation from 45
thc = atan2d(Cnh(:,2)-Cimg(1,2),Cnh(:,1)-Cimg(1,1));
th = mod(thc,90);
idxc = th>(45-thtol) & th<(45+thtol); % find angles which are roughly multiples of 45
Ccorners = Cnh(idxc,:); % there should only be 4 rows here
thc = thc(idxc);

% walk outwards along those diagonal trajectories
sqsize = 9; % the square is 9x9 points
for kann = 1:floor(sqsize/2)-1
    for k = 1:4 % assuming Ccorners has 4 rows
        % find the 8 peaks closest to this peak
        thiscorner = Ccorners(k,:);
        D = sqrt(sum((Cpeaks-thiscorner).^2,2));
        [~,idx] = mink(D,9);
        Cnh = Cpeaks(idx(2:end),:); % ignore self-distance
        
        % find point which is oriented roughly in the same direction
        % as the prior point was oriented WRT to the image center
        th = atan2d(Cnh(:,2)-thiscorner(1,2),Cnh(:,1)-thiscorner(1,1));
        [~,idx] = min(abs(th-thc(k)));
        Ccorners(k,:) = Cnh(idx,:); % there should only be one row here
        thc(k) = th(idx);
    end
end

% show where the estimated corners are
plot(Ccorners(:,1),Ccorners(:,2),'mo','linewidth',2,'markersize',20)
	
% sort by angle
[~,idx] = sort(thc,'ascend');
Ccorners = Ccorners(idx,:);

% create polygonal ROI object
ROI = images.roi.Polygon(gca);
ROI.Position = Ccorners;

At that point, you could use createMask() as before to get a logical mask from ROI, or whatever else is needed.

Answer 2

[OLUWAFEMI AKINMOLAYAN]

Great work DGM. Thanks alot. This works for most of the images. I just needed to change the peakmask for some.

However, very few of the images (one of the images and outcome is attached) have their grid lines shifted a bit to one side. This somehow affected the outcome. I tried to tweak the code but it gave me errors. I think one would have to use the grid line to mask the ROI in this case. I might be wrong. You will best know how to tweak the code to get the desire result if its possible.

Also, I have a question. In an attempt to crop out the ROI to give the rectangular image "B", do you know what information might be lost?

I hope i am not troubling you with too much of my questions? lols

For some reason i couldn't add comment to previous answer. I had to comment under more answers. I hope you see this.

Thanks once again

Answer 3

[DGM]

I was assuming that the green marks weren't going to be the ROI delimiters anymore, but if they still are, you can't use image geometry to estimate the ROI center. You still have to deal with finding the marks -- but maybe you can get away without needing to do as much masking.

Still, finding the region center is questionable if there's significant perspective distortion. I slapped it in a loop and made it at least try to fix itself if it runs into the edge of the image, but I don't doubt that this can still break easily. It works for all the images so far at least. It doesn't seem too picky about the masking yet.

% the input image
%A = imread('Picture1.png');
%A = imread('picture4.jpg');
%A = imread('40A_12.jpg');
%A = imread('25A_13.jpg');
%A = imread('10A_22.jpg');
%A = imread('10A_22_skew.jpg');
A = imread('40A_23.jpg');

sqsize = 9; % the square is 9x9 points
thtol = 15; % allowable angle deviation from 45

% try to find dark-ish areas that might be marks
V = max(A,[],3);
roimask = V<100;
roimask = bwareaopen(roimask,100); % despeckle

% find the point that's furthest away from any dark spots
D = bwdist(roimask);
[mx,idx] = max(D(:));
[r c] = ind2sub(size(D),idx);
Cimg = [c r];

% get mask of all light peaks
peakmask = mean(A,3)>200; % HSI intensity
peakmask = imclearborder(peakmask);
peakmask = bwareaopen(peakmask,100);

imshow(peakmask); hold on

% get centroids
S = regionprops(peakmask,'centroid');
Cpeaks = vertcat(S.Centroid);

% distance from every object to every other object
D = sqrt((Cpeaks(:,1)-Cpeaks(:,1).').^2 + (Cpeaks(:,2)-Cpeaks(:,2).').^2);
D(D<1E-6) = NaN; % remove self-distances
Dmean = mean(min(D,[],2)); % average grid spacing

% try to find corners
% if we run into the image edges, adjust center and try again
% adjust at most twice before continuing/failing
Ccorrection = 0;
for attempt = 1:3
    [Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol);
    if Ccorrection == 0; break; end
    
    if attempt == 1
        % distance from every object to every other object
        D = sqrt((Cpeaks(:,1)-Cpeaks(:,1).').^2 + (Cpeaks(:,2)-Cpeaks(:,2).').^2);
        D(D<1E-6) = NaN; % remove self-distances
        Dmean = mean(min(D,[],2)); % average grid spacing
    end
    
    % adjust estimate of image center
    Cimg(2) = Cimg(2) - Ccorrection*Dmean;
end

% show where the estimated corners are
plot(Ccorners(:,1),Ccorners(:,2),'mo','linewidth',2,'markersize',20)

% show ROI on original image
figure
imshow(A); hold on

% sort by angle
[~,idx] = sort(thc,'ascend');
Ccorners = Ccorners(idx,:);

% create polygonal ROI object
ROI = images.roi.Polygon(gca);
ROI.Position = Ccorners;

function [Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol)
    % find the peak closest to the image center
    D = sqrt(sum((Cpeaks-Cimg).^2,2));
    [~,idx] = min(D);
    Ccenter = Cpeaks(idx,:);
    
    % mark the center peak
    plot(Ccenter(1),Ccenter(2),'c*')
    
    % find the 8 peaks closest to this peak
    D = sqrt(sum((Cpeaks-Ccenter).^2,2));
    [~,idx] = mink(D,9);
    Cnh = Cpeaks(idx(2:end),:); % ignore self-distance
    
    % find points which are oriented roughly 45 
    thc = atan2d(Cnh(:,2)-Cimg(1,2),Cnh(:,1)-Cimg(1,1));
    th = mod(thc,90);
    idxc = th>(45-thtol) & th<(45+thtol); % find angles which are roughly multiples of 45
    Ccorners = Cnh(idxc,:); % there should only be 4 rows here
    thc = thc(idxc);
    
    % walk outwards along those diagonal trajectories
    for kann = 1:floor(sqsize/2)-1
        for k = 1:4 % assuming Ccorners has 4 rows
            % find the 8 peaks closest to this peak
            thiscorner = Ccorners(k,:);
            D = sqrt(sum((Cpeaks-thiscorner).^2,2));
            [~,idx] = mink(D,9);
            Cnh = Cpeaks(idx(2:end),:); % ignore self-distance
            
            % find point which is oriented roughly in the same direction
            % as the prior point was oriented WRT to the image center
            th = atan2d(Cnh(:,2)-thiscorner(1,2),Cnh(:,1)-thiscorner(1,1));
            [minth idx] = min(abs(th-thc(k)));
            
            % if there are no neighboring points in this direction, 
            % that's probably because we just ran into the image edge
            if minth > thtol
                if thc(k)<0 
                    % if walking downward, try shifting initial center upwards
                    Ccorrection = -1;
                else
                    % if walking upward, try shifting initial center downward
                    Ccorrection = 1;
                end
                return;
            else
                Ccorrection = 0;
            end
            
            % update outputs
            Ccorners(k,:) = Cnh(idx,:); % there should only be one row here
            thc(k) = th(idx);
        end
    end
end

If you're cropping out B (and transforming it, etc), you obviously lose information outside of the ROI. You're also losing the exact original values within the ROI, since it's being interpolated. Whether that interpolation is meaningful for your analysis, I don't know.

[OLUWAFEMI AKINMOLAYAN]

Thanks DGM.

The current code works. However, I could not create mask with it. It was giving me the error below.

Error using images.roi.internal.mixin.CreateMask/validateInputs

createMask expects a current figure containing an image.

Error in images.roi.internal.mixin.CreateMask/createMask

Error in (line 90)

mask = createMask(ROI);

[Image Analyst]

Try

mask = ROI.createMask;

or else

mask = poly2mask(Ccorners(:, 1), Ccorners(:, 2), imageRows, imageColumns);

[OLUWAFEMI AKINMOLAYAN]

@Image Analyst, None of them works.

Answer 4

[DGM]

It seems to work for me. The figure in which the ROI object is created needs to be present and needs to contain only one image.

% the input image
%A = imread('Picture1.png');
%A = imread('picture4.jpg');
%A = imread('40A_12.jpg');
%A = imread('25A_13.jpg');
%A = imread('10A_22.jpg');
%A = imread('10A_22_skew.jpg');
A = imread('40A_23.jpg');

sqsize = 9; % the square is 9x9 points
thtol = 15; % allowable angle deviation from 45

% try to find dark-ish areas that might be marks
V = max(A,[],3);
roimask = V<100;
roimask = bwareaopen(roimask,100); % despeckle

% find the point that's furthest away from any dark spots
D = bwdist(roimask);
[mx,idx] = max(D(:));
[r c] = ind2sub(size(D),idx);
Cimg = [c r];

% get mask of all light peaks
peakmask = mean(A,3)>200; % HSI intensity
peakmask = imclearborder(peakmask);
peakmask = bwareaopen(peakmask,100);
imshow(peakmask); hold on

% get centroids
S = regionprops(peakmask,'centroid');
Cpeaks = vertcat(S.Centroid);

% distance from every object to every other object
D = sqrt((Cpeaks(:,1)-Cpeaks(:,1).').^2 + (Cpeaks(:,2)-Cpeaks(:,2).').^2);
D(D<1E-6) = NaN; % remove self-distances
Dmean = mean(min(D,[],2)); % average grid spacing

% try to find corners
% if we run into the image edges, adjust center and try again
% adjust at most twice before continuing/failing
Ccorrection = 0;
for attempt = 1:3
    [Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol);
    if Ccorrection == 0; break; end
    
    if attempt == 1
        % distance from every object to every other object
        D = sqrt((Cpeaks(:,1)-Cpeaks(:,1).').^2 + (Cpeaks(:,2)-Cpeaks(:,2).').^2);
        D(D<1E-6) = NaN; % remove self-distances
        Dmean = mean(min(D,[],2)); % average grid spacing
    end
    
    % adjust estimate of image center
    Cimg(2) = Cimg(2) - Ccorrection*Dmean;
end

% show where the estimated corners are
plot(Ccorners(:,1),Ccorners(:,2),'mo','linewidth',2,'markersize',20)

% show ROI on original image
figure
imshow(A); hold on

% sort by angle
[~,idx] = sort(thc,'ascend');
Ccorners = Ccorners(idx,:);

% create polygonal ROI object
ROI = images.roi.Polygon(gca);
ROI.Position = Ccorners;

% create mask from ROI object
mask = createMask(ROI);

% show mask
figure
imshow(mask)

function [Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol)
    % find the peak closest to the image center
    D = sqrt(sum((Cpeaks-Cimg).^2,2));
    [~,idx] = min(D);
    Ccenter = Cpeaks(idx,:);
    
    % mark the center peak
    plot(Ccenter(1),Ccenter(2),'c*')
    
    % find the 8 peaks closest to this peak
    D = sqrt(sum((Cpeaks-Ccenter).^2,2));
    [~,idx] = mink(D,9);
    Cnh = Cpeaks(idx(2:end),:); % ignore self-distance
    
    % find points which are oriented roughly 45 
    thc = atan2d(Cnh(:,2)-Cimg(1,2),Cnh(:,1)-Cimg(1,1));
    th = mod(thc,90);
    idxc = th>(45-thtol) & th<(45+thtol); % find angles which are roughly multiples of 45
    Ccorners = Cnh(idxc,:); % there should only be 4 rows here
    thc = thc(idxc);
    
    % walk outwards along those diagonal trajectories
    for kann = 1:floor(sqsize/2)-1
        for k = 1:4 % assuming Ccorners has 4 rows
            % find the 8 peaks closest to this peak
            thiscorner = Ccorners(k,:);
            D = sqrt(sum((Cpeaks-thiscorner).^2,2));
            [~,idx] = mink(D,9);
            Cnh = Cpeaks(idx(2:end),:); % ignore self-distance
            
            % find point which is oriented roughly in the same direction
            % as the prior point was oriented WRT to the image center
            th = atan2d(Cnh(:,2)-thiscorner(1,2),Cnh(:,1)-thiscorner(1,1));
            [minth idx] = min(abs(th-thc(k)));
            
            % if there are no neighboring points in this direction, 
            % that's probably because we just ran into the image edge
            if minth > thtol
                if thc(k)<0 
                    % if walking downward, try shifting initial center upwards
                    Ccorrection = -1;
                else
                    % if walking upward, try shifting initial center downward
                    Ccorrection = 1;
                end
                return;
            else
                Ccorrection = 0;
            end
            
            % update outputs
            Ccorners(k,:) = Cnh(idx,:); % there should only be one row here
            thc(k) = th(idx);
        end
    end
end

[OLUWAFEMI AKINMOLAYAN]

Thanks @DGM. It works now.

The peak mask threshold values varies extensively with most of the images in the new code.

When i run the code with few images. It gave the error below until i either reduce or increase the peakmask value in the code. This sometimes gives wierd output. Example image is attached.

Index in position 1 exceeds array bounds (must not exceed 2).

Error in DGMREFINED_3>findcorners (line 112)

thiscorner = Ccorners(k,:);

Error in DGMREFINED_3 (line 55)

[Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol);

[Image Analyst]

@OLUWAFEMI AKINMOLAYAN you can download an interactive/visual thresholding app from my File Exchange:

Image Analyst's File Exchange

[DGM]

I changed the marker threshold a lot higher. With the new (half) markerless method, having a tight mask on the markers isn't critical so long as the markers are found. I added a couple checks to provide clearer information about what broke when it breaks. I haven't changed the peak threshold since the last setting (200).

Mind you, this is still written in a very ad-hoc fashion. You'll probably want to incorporate some plotting routines into this so you can keep track of the intermediate masks that are being generated. That way you'll know what's going on when it breaks. If I wrote it to produce a big figure with a bunch of subplots, it would basically be unreadably tiny when I ran it on the forum.

% the input image
%A = imread('Picture1.png');
%A = imread('picture4.jpg');
%A = imread('40A_12.jpg');
%A = imread('25A_13.jpg');
%A = imread('10A_22.jpg');
%A = imread('10A_22_skew.jpg');
%A = imread('40A_23.jpg');
A = imread('40A_35.jpg');

sqsize = 9; % the square is 9x9 points
thtol = 15; % allowable angle deviation from 45
markerthr = 150; % threshold used to find markers in V
peakthr = 200; % threshold used to find bright peaks in I

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% try to find dark-ish areas that might be marks
V = max(A,[],3);
roimask = V<markerthr;
roimask = bwareaopen(roimask,100); % despeckle
imshow(roimask)

% find the point that's furthest away from any dark spots
D = bwdist(roimask);
[mx,idx] = max(D(:));
[r c] = ind2sub(size(D),idx);
Cimg = [c r];
imshow(D,[])

% get mask of all light peaks
peakmask = mean(A,3)>peakthr;
peakmask = imclearborder(peakmask);
peakmask = bwareaopen(peakmask,100);

% mark the estimated region center
imshow(peakmask); hold on
plot(Cimg(1),Cimg(2),'yd','linewidth',3,'markersize',10)

% get centroids
S = regionprops(peakmask,'centroid');
Cpeaks = vertcat(S.Centroid);

% distance from every object to every other object
D = sqrt((Cpeaks(:,1)-Cpeaks(:,1).').^2 + (Cpeaks(:,2)-Cpeaks(:,2).').^2);
D(D<1E-6) = NaN; % remove self-distances
Dmean = mean(min(D,[],2)); % average grid spacing

% try to find corners
% if we run into the image edges, adjust center and try again
% adjust at most twice before continuing/failing
Ccorrection = 0;
for attempt = 1:3
    [Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol);
    if Ccorrection == 0; break; end
    
    if attempt == 1
        % distance from every object to every other object
        D = sqrt((Cpeaks(:,1)-Cpeaks(:,1).').^2 + (Cpeaks(:,2)-Cpeaks(:,2).').^2);
        D(D<1E-6) = NaN; % remove self-distances
        Dmean = mean(min(D,[],2)); % average grid spacing
    end
    
    % adjust estimate of image center
    Cimg(2) = Cimg(2) - Ccorrection*Dmean;
end

% show where the estimated corners are
plot(Ccorners(:,1),Ccorners(:,2),'mo','linewidth',2,'markersize',20)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% show ROI on original image
figure
imshow(A); hold on

% sort by angle
[~,idx] = sort(thc,'ascend');
Ccorners = Ccorners(idx,:);

% create polygonal ROI object
ROI = images.roi.Polygon(gca);
ROI.Position = Ccorners;

% create mask from ROI object
mask = createMask(ROI);


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol)
% find the peak closest to the image center
D = sqrt(sum((Cpeaks-Cimg).^2,2));
[~,idx] = min(D);
Ccenter = Cpeaks(idx,:);

% mark the center peak
plot(Ccenter(1),Ccenter(2),'c*','linewidth',3)

% find the 8 peaks closest to this peak
D = sqrt(sum((Cpeaks-Ccenter).^2,2));
[~,idx] = mink(D,9);
Cnh = Cpeaks(idx(2:end),:); % ignore self-distance

% find points which are oriented roughly 45 
thc = atan2d(Cnh(:,2)-Cimg(1,2),Cnh(:,1)-Cimg(1,1));
th = mod(thc,90);
idxc = th>(45-thtol) & th<(45+thtol); % find angles which are roughly multiples of 45
if nnz(idxc)~=4
    sort(th)
    error('could not find four neighboring corners')
end
Ccorners = Cnh(idxc,:); % there should only be 4 rows here
thc = thc(idxc);



% walk outwards along those diagonal trajectories
for kann = 1:floor(sqsize/2)-1
    for k = 1:4 % assuming Ccorners has 4 rows
        % find the 8 peaks closest to this peak
        thiscorner = Ccorners(k,:);
        D = sqrt(sum((Cpeaks-thiscorner).^2,2));
        [~,idx] = mink(D,9);
        Cnh = Cpeaks(idx(2:end),:); % ignore self-distance
        
        % find point which is oriented roughly in the same direction
        % as the prior point was oriented WRT to the image center
        th = atan2d(Cnh(:,2)-thiscorner(1,2),Cnh(:,1)-thiscorner(1,1));
        [minth idx] = min(abs(th-thc(k)));
        
        % if there are no neighboring points in this direction, 
        % that's probably because we just ran into the image edge
        if minth > thtol
            if thc(k)<0 
                % if walking downward, try shifting initial center upwards
                warning('ran into bottom edge of image; adjusting center')
                Ccorrection = -1;
            else
                % if walking upward, try shifting initial center downward
                warning('ran into top edge of image; adjusting center')
                Ccorrection = 1;
            end
            return;
        else
            Ccorrection = 0;
        end
        
        % update outputs
        Ccorners(k,:) = Cnh(idx,:); % there should only be one row here
        thc(k) = th(idx);
    end
end
end

[DGM]

... I have no idea how I ended up starting another answer, but I guess I did. Sometimes it's hard to tell what's going on when the page is super laggy, and it's been extremely laggy the last few days.

[OLUWAFEMI AKINMOLAYAN]

Yeah. I make the mistake of starting another answer sometimes.

Thanks. If i get the error, "could not find four neighboring corners", what adjustment do i need to do. I get it for some of the image. E.g is attached

[DGM]

I modified it so that the initial angles aren't found by angle discrimination, but by distance maximization within the initial neighborhood around the central peak. It works on N_B_42 and the skewed copy I used in this example. I tried to improve the ROI centerfinding a bit.

% the input image
%A = imread('Picture1.png');
%A = imread('picture4.jpg');
%A = imread('40A_12.jpg');
%A = imread('25A_13.jpg');
%A = imread('10A_22.jpg');
%A = imread('10A_22_skew.jpg');
%A = imread('40A_23.jpg');
%A = imread('40A_35.jpg');
A = imread('40A_35_skew.jpg');
%A = imread('N_B_42.jpg');

sqsize = 9; % the square is 9x9 points
thtol = 25; % allowable difference in segment angles along a diagonal
markerthr = 150; % threshold used to find markers in V
peakthr = 200; % threshold used to find bright peaks in I

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% try to find dark-ish areas that might be marks
V = max(A,[],3);
roimask = V<markerthr;
roimask = bwareaopen(roimask,100); % despeckle
imshow(roimask)

% find the point that's furthest away from any dark spots
% pad array to make sure that border-connected regions are less likely to get picked up
% i'm reluctant to mask or weight based on position within image
D = bwdist(padarray(roimask,[1 1],1,'both'));
D = D(2:end-1,2:end-1);
[mx,idx] = max(D(:));
[r c] = ind2sub(size(D),idx);
Cimg = [c r];
imshow(D,[])

% get mask of all light peaks
peakmask = mean(A,3)>peakthr;
peakmask = imclearborder(peakmask);
peakmask = bwareaopen(peakmask,100);

% mark the estimated region center
imshow(peakmask); hold on
plot(Cimg(1),Cimg(2),'yd','linewidth',3,'markersize',10)

% get centroids
S = regionprops(peakmask,'centroid');
Cpeaks = vertcat(S.Centroid);

% distance from every object to every other object
D = sqrt((Cpeaks(:,1)-Cpeaks(:,1).').^2 + (Cpeaks(:,2)-Cpeaks(:,2).').^2);
D(D<1E-6) = NaN; % remove self-distances
Dmean = mean(min(D,[],2)); % average grid spacing

% try to find corners
% if we run into the image edges, adjust center and try again
% adjust at most twice before continuing/failing
Ccorrection = 0;
for attempt = 1:3
    [Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol);
    if Ccorrection == 0; break; end
    
    if attempt == 1
        % distance from every object to every other object
        D = sqrt((Cpeaks(:,1)-Cpeaks(:,1).').^2 + (Cpeaks(:,2)-Cpeaks(:,2).').^2);
        D(D<1E-6) = NaN; % remove self-distances
        Dmean = mean(min(D,[],2)); % average grid spacing
    end
    
    % adjust estimate of image center
    Cimg(2) = Cimg(2) - Ccorrection*Dmean;
end

% show where the estimated corners are
plot(Ccorners(:,1),Ccorners(:,2),'mo','linewidth',2,'markersize',20)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% show ROI on original image
figure
imshow(A); hold on

% sort by angle
[~,idx] = sort(thc,'ascend');
Ccorners = Ccorners(idx,:);

% create polygonal ROI object
ROI = images.roi.Polygon(gca);
ROI.Position = Ccorners;

% create mask from ROI object
mask = createMask(ROI);


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ccorners thc Ccorrection] = findcorners(Cpeaks,Cimg,sqsize,thtol)
    % find the peak closest to the image center
    D = sqrt(sum((Cpeaks-Cimg).^2,2));
    [~,idx] = min(D);
    Ccenter = Cpeaks(idx,:);
    
    % mark the center peak
    plot(Ccenter(1),Ccenter(2),'c*','linewidth',3)
    
    % find the 8 peaks closest to this peak
    D = sqrt(sum((Cpeaks-Ccenter).^2,2));
    [~,idx] = mink(D,9);
    Dnh = D(idx(2:end)); % ignore self-distance
    Cnh = Cpeaks(idx(2:end),:); % ignore self-distance
    
    % find corners by picking the four most distant neighbors
    [~,idxc] = maxk(Dnh,4);
    Ccorners = Cnh(idxc,:);
    thc = atan2d(Ccorners(:,2)-Cimg(1,2),Ccorners(:,1)-Cimg(1,1));
    
    % walk outwards along those diagonal trajectories
    for kann = 1:floor(sqsize/2)-1
        for k = 1:4 % assuming Ccorners has 4 rows
            % find the 8 peaks closest to this peak
            thiscorner = Ccorners(k,:);
            D = sqrt(sum((Cpeaks-thiscorner).^2,2));
            [~,idx] = mink(D,9);
            Cnh = Cpeaks(idx(2:end),:); % ignore self-distance
            
            % find point which is oriented roughly in the same direction
            % as the prior point was oriented WRT to the image center
            th = atan2d(Cnh(:,2)-thiscorner(1,2),Cnh(:,1)-thiscorner(1,1));
            [minth idx] = min(abs(th-thc(k)));
            
            % if there are no neighboring points in this direction, 
            % that's probably because we just ran into the image edge
            if minth > thtol
                if thc(k)<0 
                    % if walking downward, try shifting initial center upwards
                    warning('ran into bottom edge of image; adjusting center')
                    Ccorrection = -1;
                else
                    % if walking upward, try shifting initial center downward
                    warning('ran into top edge of image; adjusting center')
                    Ccorrection = 1;
                end
                return;
            else
                Ccorrection = 0;
            end
            
            % update outputs
            Ccorners(k,:) = Cnh(idx,:); % there should only be one row here
            thc(k) = th(idx);
        end
    end
end

[OLUWAFEMI AKINMOLAYAN]

Thanks a bunch @DGM . This works fine for majority of my images except for few that are badly distorted. I am currently work of populating multiple images using the code in order to save output "B" and the pixel values and coordinate data in an excel file of all the images at once. I will coomment if i hit a dead end.

[OLUWAFEMI AKINMOLAYAN]

filePattern = fullfile(Folder, '*.jpg');
jpegFiles   = dir(filePattern);

for k = 1:length(jpegFiles)
  baseFileName = jpegFiles(k).name;
  fullFileName = fullfile(Folder, baseFileName);
  fprintf('Now reading %s\n', fullFileName);
%    subplot(1,length(jpegFiles),k)
  A = imread(fullFileName);
  
  %process images and extract data...

  FileName = sprintf('C:\\Users\\%d.jpg', k)
  imwrite(B,FileName); %B is the output images
  
  thisfname = sprintf('C:\\Users\\\\myfilename_%04d.xlsx',k);
    writematrix(xy,thisfname) % xy is the output data
    
end

In the above code

The output images are saved as 1...n. I want to save the images using same name format as the input naming in order to distinct which output is for which input. Please how do i do that.

Also, the data save to excel in different xls files. I want to save the data in one excel. But in different sheets. How do i do that. Couldnt find an specific answer for this.

Or if it would be possible to save the array of data in one sheet with the column names to distinguish the different data.

[OLUWAFEMI AKINMOLAYAN]

ANY IDEA? @DGM

[DGM]

This is perhaps a start. I'm not really sure how you want to break down the filename patterns. I just inserted '_0001' (and so on) between the old filename and the extension, so 'mypicture.jpg' becomes 'mypicture_0001.jpg'. You can change the number of digits if you want. If there are parts of the original filename that you want to omit, you'll have to describe the pattern.

% you'll have to set these directories
Folder = 'sources'; % this is just where i had some test files
outFolder = './'; % the folder to place the outputs in
outfiletype = '.jpg'; % jpg isn't very good, this makes it easy to change

filePattern = fullfile(Folder, '*.jpg');
jpegFiles   = dir(filePattern);
for k = 1:length(jpegFiles)
    baseFileName = jpegFiles(k).name;
    fullFileName = fullfile(Folder, baseFileName);
    
    fprintf('Now reading %s\n', fullFileName);
    A = imread(fullFileName);
    
    %process images and extract data...
    B = A; % placeholder
    
    % write images
    % this assumes all input images are .jpg or .png; add cases as needed
    fnamenoext = regexprep(baseFileName,'(.jpg$)|(.png$)','');
    % build filename out of old filename, k, and specified path/extension
    FileName = fullfile(outFolder,sprintf('%s_%04d%s',fnamenoext,k,outfiletype))
    imwrite(B,FileName); %B is the output images
    
    % not sure if you want to put these in the same directory
    % but if so, you could do something similar to construct the path
    %thisfname = sprintf('C:\\Users\\\\myfilename_%04d.xlsx',k);
    %writematrix(xy,thisfname) % xy is the output data
end

[OLUWAFEMI AKINMOLAYAN]

THANK YOU VERY MUCH

[OLUWAFEMI AKINMOLAYAN]

@DGM, Is there a way to make the ROI the same size, that is the pixel coordinate arrays are the same size irrespective of the image. Thank in advance.

[DGM]

I'm sure anything is possible, but at this point it's not really clear how that would change the task requirements. The last incarnation of the script finds an ROI as a irregular quadrilateral with vertices located on image features. In what manner is the ROI size to be fixed (e.g. area/height/width)? Does fixing the size also fix the shape or orientation? What parts of the fixed ROI would correspond to the image features?