extracting mean pixel value from an image
Afficher commentaires plus anciens
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
Réponse acceptée
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)
% 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.
% 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.
25 commentaires
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)
% find median color
mediancolor = median(roipix,1)
% 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
le 20 Mai 2022
DGM
le 20 Mai 2022
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
le 20 Mai 2022
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]
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]
DGM
le 21 Mai 2022
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
le 21 Mai 2022
OLUWAFEMI AKINMOLAYAN
le 21 Mai 2022
OLUWAFEMI AKINMOLAYAN
le 21 Mai 2022
DGM
le 21 Mai 2022
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.
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
le 21 Mai 2022
OLUWAFEMI AKINMOLAYAN
le 21 Mai 2022
DGM
le 22 Mai 2022
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
le 22 Mai 2022
DGM
le 22 Mai 2022
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
le 23 Mai 2022
OLUWAFEMI AKINMOLAYAN
le 23 Mai 2022
OLUWAFEMI AKINMOLAYAN
le 23 Mai 2022
OLUWAFEMI AKINMOLAYAN
le 27 Mai 2022
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
le 28 Mai 2022
DGM
le 28 Mai 2022
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
le 2 Juin 2022
DGM
le 2 Juin 2022
... 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...
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.
Plus de réponses (3)
OLUWAFEMI AKINMOLAYAN
le 3 Juin 2022
0 votes
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.
3 commentaires
OLUWAFEMI AKINMOLAYAN
le 4 Juin 2022
Image Analyst
le 4 Juin 2022
Try
mask = ROI.createMask;
or else
mask = poly2mask(Ccorners(:, 1), Ccorners(:, 2), imageRows, imageColumns);
OLUWAFEMI AKINMOLAYAN
le 4 Juin 2022
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
13 commentaires
OLUWAFEMI AKINMOLAYAN
le 4 Juin 2022
Image Analyst
le 5 Juin 2022
@OLUWAFEMI AKINMOLAYAN you can download an interactive/visual thresholding app from my File Exchange:
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
le 5 Juin 2022
... 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
le 5 Juin 2022
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
le 7 Juin 2022
OLUWAFEMI AKINMOLAYAN
le 7 Juin 2022
OLUWAFEMI AKINMOLAYAN
le 9 Juin 2022
DGM
le 9 Juin 2022
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
le 16 Juin 2022
OLUWAFEMI AKINMOLAYAN
le 30 Juin 2022
Modifié(e) : OLUWAFEMI AKINMOLAYAN
le 30 Juin 2022
DGM
le 1 Juil 2022
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?
Catégories
En savoir plus sur Image Arithmetic dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
