1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
| tic
clear all;close all;clc;
%% Sélection de l'image XXXXX
% dossier_img = ('XXX');
% nom_img = ('XXX');
%% Création chemin et lecture de l'image
chemin_img = fullfile(dossier_img,nom_img);
I = imread(chemin_img);
I = rgb2gray(I);
captionFontSize = 14;
%% XXXXXXXXXXXX
subplot(3, 3, 1);
imshow(I);
% Maximize the figure window.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Force it to display RIGHT NOW (otherwise it might not display until it's all done, unless you've stopped at a breakpoint.)
drawnow;
caption = sprintf('Original');
title(caption, 'FontSize', captionFontSize);
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
% Just for fun, let's get its histogram and display it.
[pixelCount, grayLevels] = imhist(I);
subplot(3, 3, 2);
bar(pixelCount);
title('Histogram of original image', 'FontSize', captionFontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.
grid on;
thresholdValue = 95; (c'est fixé un peu arbitrairement pour l'instant)
binaryImage = I > thresholdValue; % Bright objects will be chosen if you use >.
% ========== IMPORTANT OPTION ============================================================
% Use < if you want to find dark objects instead of bright objects.
% binaryImage = originalImage < thresholdValue; % Dark objects will be chosen if you use <.
% Do a "hole fill" to get rid of any background pixels or "holes" inside the blobs.
binaryImage = imfill(binaryImage, 'holes');
% Show the threshold as a vertical red bar on the histogram.
hold on;
maxYValue = ylim;
line([thresholdValue, thresholdValue], maxYValue, 'Color', 'r');
% Place a text label on the bar chart showing the threshold.
annotationText = sprintf('Thresholded at %d gray levels', thresholdValue);
% For text(), the x and y need to be of the data class "double" so let's cast both to double.
text(double(thresholdValue + 5), double(0.5 * maxYValue(2)), annotationText, 'FontSize', 10, 'Color', [0 .5 0]);
text(double(thresholdValue - 70), double(0.94 * maxYValue(2)), 'Background', 'FontSize', 10, 'Color', [0 0 .5]);
text(double(thresholdValue + 50), double(0.94 * maxYValue(2)), 'Foreground', 'FontSize', 10, 'Color', [0 0 .5]);
% Display the binary image.
subplot(3, 3, 3);
imshow(binaryImage);
title('Binary Image, obtained by thresholding', 'FontSize', captionFontSize);
% Identify individual blobs by seeing which pixels are connected to each other.
% Each group of connected pixels will be given a label, a number, to identify it and distinguish it from the other blobs.
% Do connected components labeling with either bwlabel() or bwconncomp().
labeledImage = bwlabel(binaryImage, 8); % Label each blob so we can make measurements of it (j'aurais vraiment aimé mettre plus, genre fixer une valeur très haute pour avoir uniquement de grandes zones mais dans l'aide c'est écrit 4 ou 8...)
% labeledImage is an integer-valued image where all pixels in the blobs have values of 1, or 2, or 3, or ... etc.
subplot(3, 3, 4);
imshow(labeledImage, []); % Show the gray scale image.
title('Labeled Image, from bwlabel()', 'FontSize', captionFontSize);
% Let's assign each blob a different color to visually show the user the distinct blobs.
coloredLabels = label2rgb (labeledImage, 'hsv', 'k', 'shuffle'); % pseudo random color labels
% coloredLabels is an RGB image. We could have applied a colormap instead (but only with R2014b and later)
subplot(3, 3, 5);
imshow(coloredLabels);
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
caption = sprintf('Labels colorés\n Numérotation de haut en bas, de gauche à droite.');
title(caption, 'FontSize', captionFontSize);
% Get all the blob properties. Can only pass in originalImage in version R2008a and later.
blobMeasurements = regionprops(labeledImage, I, 'all');
numberOfBlobs = size(blobMeasurements, 1);
% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the coins on the original grayscale image using the coordinates returned by bwboundaries.
subplot(3, 3, 6);
imshow(I);
title('Outlines, from bwboundaries()', 'FontSize', captionFontSize);
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
hold on;
boundaries = bwboundaries(binaryImage);
numberOfBoundaries = size(boundaries, 1);
for k = 1 : numberOfBoundaries
thisBoundary = boundaries{k};
plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 2);
end
hold off;
textFontSize = 14; % Used to control size of "blob number" labels put atop the image.
labelShiftX = -7; % Used to align the labels in the centers of the coins.
blobECD = zeros(1, numberOfBlobs);
% Print header line in the command window.
fprintf(1,'Blob # Mean Intensity Area Perimeter Centroid Diameter\n');
% Loop over all blobs printing their measurements to the command window.
for k = 1 : numberOfBlobs % Loop through all blobs.
% Find the mean of each blob. (R2008a has a better way where you can pass the original image
% directly into regionprops. The way below works for all versions including earlier versions.)
thisBlobsPixels = blobMeasurements(k).PixelIdxList; % Get list of pixels in current blob.
meanGL = mean(I(thisBlobsPixels)); % Find mean intensity (in original image!)
meanGL2008a = blobMeasurements(k).MeanIntensity; % Mean again, but only for version >= R2008a
blobArea = blobMeasurements(k).Area; % Get area.
blobPerimeter = blobMeasurements(k).Perimeter; % Get perimeter.
blobCentroid = blobMeasurements(k).Centroid; % Get centroid one at a time
blobECD(k) = sqrt(4 * blobArea / pi); % Compute ECD - Equivalent Circular Diameter.
fprintf(1,'#%2d %17.1f %11.1f %8.1f %8.1f %8.1f % 8.1f\n', k, meanGL, blobArea, blobPerimeter, blobCentroid, blobECD(k));
% Put the "blob number" labels on the "boundaries" grayscale image.
text(blobCentroid(1) + labelShiftX, blobCentroid(2), num2str(k), 'FontSize', textFontSize, 'FontWeight', 'Bold');
end
% Now, I'll show you another way to get centroids.
% We can get the centroids of ALL the blobs into 2 arrays,
% one for the centroid x values and one for the centroid y values.
allBlobCentroids = [blobMeasurements.Centroid];
centroidsX = allBlobCentroids(1:2:end-1);
centroidsY = allBlobCentroids(2:2:end);
% Put the labels on the rgb labeled image also.
subplot(3, 3, 5);
for k = 1 : numberOfBlobs % Loop through all blobs.
text(centroidsX(k) + labelShiftX, centroidsY(k), num2str(k), 'FontSize', textFontSize, 'FontWeight', 'Bold');
end
% Now I'll demonstrate how to select certain blobs based using the ismember() function.
% Let's say that we wanted to find only those blobs
% with an intensity between 150 and 220 and an area less than 2000 pixels.
% This would give us the three brightest dimes (the smaller coin type).
allBlobIntensities = [blobMeasurements.MeanIntensity];
allBlobAreas = [blobMeasurements.Area];
% Get a list of the blobs that meet our criteria and we need to keep.
% These will be logical indices - lists of true or false depending on whether the feature meets the criteria or not.
% for example [1, 0, 0, 1, 1, 0, 1, .....]. Elements 1, 4, 5, 7, ... are true, others are false.
allowableIntensityIndexes = (allBlobIntensities > 20) & (allBlobIntensities < 220); (fixé là aussi un peu arbitrairement...)
allowableAreaIndexes = allBlobAreas > 400000; % Take the small objects. (j'ai fixé 400 000 pixels pour mon image de 1MPixels, mais sur 256MPixels, il faudrait que je fasse un essai. Bon le ratio sera le même donc une règle de trois et c'est bon mais j'aurais aimé éviter de calculer tout ce qu'il y a avant pour rien....)
% Now let's get actual indexes, rather than logical indexes, of the features that meet the criteria.
% for example [1, 4, 5, 7, .....] to continue using the example from above.
keeperIndexes = find(allowableIntensityIndexes & allowableAreaIndexes);
% Extract only those blobs that meet our criteria, and
% eliminate those blobs that don't meet our criteria.
% Note how we use ismember() to do this. Result will be an image - the same as labeledImage but with only the blobs listed in keeperIndexes in it.
keeperBlobsImage = ismember(labeledImage, keeperIndexes);
% Re-label with only the keeper blobs kept.
labeledDimeImage = bwlabel(keeperBlobsImage, 8); % Label each blob so we can make measurements of it
% Now we're done. We have a labeled image of blobs that meet our specified criteria.
subplot(3, 3, 7);
imshow(labeledDimeImage, []);
axis image;
title('"Keeper" blobs (3 brightest dimes in a re-labeled image)', 'FontSize', captionFontSize);
% Plot the centroids in the original image in the upper left.
% Dimes will have a red cross, nickels will have a blue X.
message = sprintf('Now I will plot the centroids over the original image in the upper left.\nPlease look at the upper left image.');
reply = questdlg(message, 'Plot Centroids?', 'OK', 'Cancel', 'Cancel');
% Note: reply will = '' for Upper right X, 'OK' for OK, and 'Cancel' for Cancel.
if strcmpi(reply, 'Cancel')
return;
end
subplot(3, 3, 1);
hold on; % Don't blow away image.
for k = 1 : numberOfBlobs % Loop through all keeper blobs.
% Identify if blob #k is a dime or nickel.
itsADime = allBlobAreas(k) < 2200; % Dimes are small.
if itsADime
% Plot dimes with a red +.
plot(centroidsX(k), centroidsY(k), 'r+', 'MarkerSize', 10, 'LineWidth', 2);
else
% Plot dimes with a blue x.
plot(centroidsX(k), centroidsY(k), 'bx', 'MarkerSize', 10, 'LineWidth', 2);
end
end
% Now use the keeper blobs as a mask on the original image.
% This will let us display the original image in the regions of the keeper blobs.
maskedImageDime = I; % Simply a copy at first.
maskedImageDime(~keeperBlobsImage) = 0; % Set all non-keeper pixels to zero.
subplot(3, 3, 8);
imshow(maskedImageDime);
axis image;
title('Only the 3 brightest dimes from the original image', 'FontSize', captionFontSize);
% Now let's get the nickels (the larger coin type).
keeperIndexes = find(allBlobAreas > 2000); % Take the larger objects.
% Note how we use ismember to select the blobs that meet our criteria.
nickelBinaryImage = ismember(labeledImage, keeperIndexes);
% Let's get the nickels from the original grayscale image, with the other non-nickel pixels blackened.
% In other words, we will create a "masked" image.
maskedImageNickel = I; % Simply a copy at first.
maskedImageNickel(~nickelBinaryImage) = 0; % Set all non-nickel pixels to zero.
subplot(3, 3, 9);
imshow(maskedImageNickel, []);
axis image;
title('Only the nickels from the original image', 'FontSize', captionFontSize);
elapsedTime = toc;
% Alert user that the demo is done and give them the option to save an image.
message = sprintf('Done making measurements of the features.\n\nElapsed time = %.2f seconds.', elapsedTime);
message = sprintf('%s\n\nCheck out the figure window for the images.\nCheck out the command window for the numerical results.', message);
message = sprintf('%s\n\nDo you want to save the pseudo-colored image?', message);
reply = questdlg(message, 'Save image?', 'Yes', 'No', 'No');
% Note: reply will = '' for Upper right X, 'Yes' for Yes, and 'No' for No.
if strcmpi(reply, 'Yes')
% Ask user for a filename.
FilterSpec = {'*.PNG', 'PNG Images (*.png)'; '*.tif', 'TIFF images (*.tif)'; '*.*', 'All Files (*.*)'};
DialogTitle = 'Save image file name';
% Get the default filename. Make sure it's in the folder where this m-file lives.
% (If they run this file but the cd is another folder then pwd will show that folder, not this one.
thisFile = mfilename('fullpath');
[thisFolder, baseFileName, ext] = fileparts(thisFile);
DefaultName = sprintf('%s/%s.tif', thisFolder, baseFileName);
[fileName, specifiedFolder] = uiputfile(FilterSpec, DialogTitle, DefaultName);
if fileName ~= 0
% Parse what they actually specified.
[folder, baseFileName, ext] = fileparts(fileName);
% Create the full filename, making sure it has a tif filename.
fullImageFileName = fullfile(specifiedFolder, [baseFileName '.tif']);
% Save the labeled image as a tif image.
imwrite(uint8(coloredLabels), fullImageFileName);
% Just for fun, read image back into the imtool utility to demonstrate that tool.
tifimage = imread(fullImageFileName);
imtool(tifimage, []);
end
end
message = sprintf('Would you like to crop out each coin to individual images?');
reply = questdlg(message, 'Extract Individual Images?', 'Yes', 'No', 'Yes');
% Note: reply will = '' for Upper right X, 'Yes' for Yes, and 'No' for No.
if strcmpi(reply, 'Yes')
figure; % Create a new figure window.
% Maximize the figure window.
set(gcf, 'Units','Normalized','OuterPosition',[0 0 1 1]);
for k = 1 : numberOfBlobs % Loop through all blobs.
% Find the bounding box of each blob.
thisBlobsBoundingBox = blobMeasurements(k).BoundingBox; % Get list of pixels in current blob.
% Extract out this coin into it's own image.
subImage = imcrop(originalImage, thisBlobsBoundingBox);
% Determine if it's a dime (small) or a nickel (large coin).
if blobMeasurements(k).Area > 2200
coinType = 'nickel';
else
coinType = 'dime';
end
% Display the image with informative caption.
subplot(3, 4, k);
imshow(subImage);
caption = sprintf('Coin #%d is a %s.\nDiameter = %.1f pixels\nArea = %d pixels', ...
k, coinType, blobECD(k), blobMeasurements(k).Area);
title(caption, 'FontSize', textFontSize);
end
% Display the MATLAB "peaks" logo.
logoFig = subplot(3, 4, 11:12);
caption = sprintf('A MATLAB Tutorial by ImageAnalyst');
text(0.5,1.15, caption, 'Color','r', 'FontSize', 18, 'FontWeight','b', 'HorizontalAlignment', 'Center') ;
positionOfLowerRightPlot = get(logoFig, 'position');
L = 40*membrane(1,25);
logoax = axes('CameraPosition', [-193.4013 -265.1546 220.4819],...
'CameraTarget',[26 26 10], ...
'CameraUpVector',[0 0 1], ...
'CameraViewAngle',9.5, ...
'DataAspectRatio', [1 1 .9],...
'Position', positionOfLowerRightPlot, ...
'Visible','off', ...
'XLim',[1 51], ...
'YLim',[1 51], ...
'ZLim',[-13 40], ...
'parent',gcf);
s = surface(L, ...
'EdgeColor','none', ...
'FaceColor',[0.9 0.2 0.2], ...
'FaceLighting','phong', ...
'AmbientStrength',0.3, ...
'DiffuseStrength',0.6, ...
'Clipping','off',...
'BackFaceLighting','lit', ...
'SpecularStrength',1, ...
'SpecularColorReflectance',1, ...
'SpecularExponent',7, ...
'Tag','TheMathWorksLogo', ...
'parent',logoax);
l1 = light('Position',[40 100 20], ...
'Style','local', ...
'Color',[0 0.8 0.8], ...
'parent',logoax);
l2 = light('Position',[.5 -1 .4], ...
'Color',[0.8 0.8 0], ...
'parent',logoax);
end
toc |
Traitement d'images avec la ToolBox
.
Répondre avec citation
Envoyé par RomainDF42
Partager