<pre>

Path: news.mathworks.com!not-for-mail
From: &lt;HIDDEN&gt;
Newsgroups: comp.soft-sys.matlab
Subject: Re: Detecting spots on a butterfly
Date: Sat, 4 Apr 2009 21:36:01 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 77
Message-ID: &lt;gr8js1$4p7$1@fred.mathworks.com&gt;
References: &lt;gr4gga$d84$1@fred.mathworks.com&gt; &lt;futBl.41$TD1.22@newsfe18.iad&gt; &lt;gr5ou1$qna$1@fred.mathworks.com&gt;
Reply-To: &lt;HIDDEN&gt;
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset=&quot;ISO-8859-1&quot;
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1238880961 4903 172.30.248.38 (4 Apr 2009 21:36:01 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sat, 4 Apr 2009 21:36:01 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1343420
Xref: news.mathworks.com comp.soft-sys.matlab:530352


To Husam Aldahiyat, the butterfly scientist:
It looks like the red channel gives the most contrast.  So I thresholded on that.  This code will find the spots for you in the one image that you posted.  It may require some tweaks to make it more robust for lots of other images.  Plus, of course, you may have to fix broken line wraps if you copy and paste this code.
Good luck,
ImageAnalyst

% Demo macro to find spots on a butterfly wing.
% by ImageAnalyst
clc;
close all;
thresholdValue = 8; % Pick something
% Get a standard MATLAB demo image.
originalImage = imread('C:\Documents and Settings\My Documents\Temporary stuff\200903153435.jpg');
subplot(2, 3, 1);
imshow(originalImage, []);
title('Original image');

redband = originalImage(:,:,1);
subplot(2, 3, 2);
imshow(redband, []);
title('Red Band');

thresholdValue = 80;
binaryImage =  redband &lt; thresholdValue;
binaryImage = imfill(binaryImage, 'holes');
subplot(2, 3, 3);
imshow(binaryImage, []);
title('binary image');
labeledImage = bwlabel(binaryImage, 8);     % Label each blob so can do calc on it
coloredLabels = label2rgb (labeledImage, 'hsv', 'k', 'shuffle'); % pseudo random color labels

subplot(2,3,4); imagesc(coloredLabels); title('Pseudo colored labels');

blobMeasurements = regionprops(labeledImage, 'all');   % Get all the blob properties.
% Get a list of the areas.
area_values = [blobMeasurements.Area];
% Filter the blobs by area.  Just keep blobs in a certain area range.
% Find out which blobs have an area between 1000 and 5000 pixels.
% Note: you could also do similar for any of the measurements, such as solidity, etc.
% You can even combine filters over different measurements.
% Use whatever the know valid range for spots area is.
indexes = find((1000 &lt;= area_values) &amp; (area_values &lt;= 5000));
spotsOnlyImage = ismember(labeledImage, indexes);
subplot(2,3,5);
imshow(spotsOnlyImage, []);
title('Spots-only image');

% Get the subset of blobs that describes only the wing-spot blobs.
blobMeasurements = blobMeasurements(indexes);   % Get the spot properties only.

% 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(2,3,6); imagesc(originalImage); title('Wing spots outlined over original image');
hold on;
boundaries = bwboundaries(spotsOnlyImage);	
numberOfWingSpots = size(blobMeasurements, 1);
for k = 1 : numberOfWingSpots
	thisBoundary = boundaries{k};
	plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 2);
end
hold off;

% List the various parameters to the command window.
fprintf(1,'Wing Spot #      Mean Intensity  Area     Perimeter  Centroid\n');
for k = 1 : numberOfWingSpots           % 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(originalImage(thisBlobsPixels));             % Find mean intensity (in original image!)
	blobArea = blobMeasurements(k).Area;		% Get area.
	blobPerimeter = blobMeasurements(k).Perimeter;		% Get perimeter.
	blobCentroid = blobMeasurements(k).Centroid;		% Get centroid.
    fprintf(1,'#%d %18.1f %11.1f %8.1f %8.1f %8.1f\n', k, meanGL, blobArea, blobPerimeter, blobCentroid);
end
message = sprintf('Done processing this image.\nThere are %d wing spots.\nMaximize and check out the figure window.\nThen check out the command window for the results.', numberOfWingSpots);
msgbox(message);
<pre>