Application de la ségmentation proposée par l'article + plein d'autres changements.
authorGreg Burri <greg.burri@gmail.com>
Sat, 16 May 2015 14:27:22 +0000 (16:27 +0200)
committerGreg Burri <greg.burri@gmail.com>
Sat, 16 May 2015 14:27:22 +0000 (16:27 +0200)
src/DetectionOfParasites.m
src/DetectionOfSchizonts.m
src/DetectionOfWhiteCells.m
src/IdentificationOfInfectedRedCells.m
src/Main.m
src/SegmentationOfRedCells.m
src/THue.m
src/loadImg.m

index ab7ffc7..c24de73 100644 (file)
@@ -15,30 +15,28 @@ function [THS, c, nucleiRadius] = DetectionOfParasites(imgRGB)
     imgH = uint8(mod(255 - imgH + hueShiftValue, 256));
     imwrite(imgH, '../output/imgH.png')
 
-    imgS = uint8(imgHSV(:, :, 2) * 255);
+    imgS = uint8(255 - 255 * imgHSV(:, :, 2));
     imwrite(imgS, '../output/imgS.png')
     
-    imgV = uint8(255 - imgHSV(:, :, 3) * 255);
+    imgV = uint8(255 - 255 * imgHSV(:, :, 3));
     imwrite(imgV, '../output/imgV.png')
 
     % p. 136: median filter to smooth the noise and 
     % area closing to enhance the bright objects and make flatter, darker and cleaner the image background.
     medianFilterWindow = [5, 5];
-    imgFiltered{1} = mmareaclose(medfilt2(imgH, medianFilterWindow), 800); % Hue.
-    % We use an opening fot the saturation because the parasites are darker than the rest of the image.
-    imgFiltered{2} = mmareaopen(medfilt2(imgS, medianFilterWindow), 800); % Saturation.
-    imgFiltered{3} = mmareaclose(medfilt2(imgV, medianFilterWindow), 800); % Value. This component isn't used.
+    imgFiltered{1} = mmareaclose(medfilt2(imgH, medianFilterWindow), 200); % Hue.
+    imgFiltered{2} = mmareaclose(medfilt2(imgS, medianFilterWindow), 200); % Saturation.
+    imgFiltered{3} = mmareaclose(medfilt2(imgV, medianFilterWindow), 200); % Value. This component isn't used.
     
     imwrite(imgFiltered{1}, '../output/imgHFiltered.png');
     imwrite(imgFiltered{2}, '../output/imgSFiltered.png');
     imwrite(imgFiltered{3}, '../output/imgVFiltered.png');
     
-
     %% Granulometry
     % We use the saturation component to find the red cells mean size.
     saturationComponent = imgFiltered{2};
 
-    redCellMaxSize = round(size(imgRGB, 2) / 35); % Radius [px], depending of the input size.
+    redCellMaxSize = uint32(size(imgRGB, 2) / 35); % Radius [px], depending of the input size.
     funVolume = @(m) sum(sum(m));
 
     volImg = funVolume(saturationComponent);
@@ -46,13 +44,13 @@ function [THS, c, nucleiRadius] = DetectionOfParasites(imgRGB)
     sizeDistribution = zeros(redCellMaxSize, 1);
     parfor k = 1:redCellMaxSize 
         SE = mmsedisk(k, '2D', 'OCTAGON'); % 'EUCLIDEAN' is more precise.
-        imgOpened = mmopen(saturationComponent, SE);
-        A = funVolume(imgOpened);
+        imgClosed = mmclose(saturationComponent, SE); % Here we use a closing because the cells are darker than the background.
+        A = funVolume(imgClosed);
         N = 1 - A / volImg;
         sizeDistribution(k) = N;
     end
 
-    % TODO : voir diff
+    % TODO : voir diff.
     patternSpectrum = zeros(redCellMaxSize - 1, 1);
     for k = 1:redCellMaxSize - 1
         patternSpectrum(k) = abs(sizeDistribution(k + 1) - sizeDistribution(k));
@@ -60,7 +58,7 @@ function [THS, c, nucleiRadius] = DetectionOfParasites(imgRGB)
 
     % The paper chooses the biggest red cell size among the possible red cell sizes. (FIXME)
     [~, c] = max(patternSpectrum);
-    nucleiRadius = c / 8; % We admit the size of a parasite is eight time smaller than a red blood cell. Find a way to extract this information from the pattern spectrum histogram. (FIXME)
+    nucleiRadius = c / 5; % We admit the size of a parasite is five time smaller than a red blood cell. Find a way to extract this information from the pattern spectrum histogram. (FIXME)
 
     bar(patternSpectrum);
 
@@ -72,11 +70,7 @@ function [THS, c, nucleiRadius] = DetectionOfParasites(imgRGB)
 
     % The functions 'mmregmax' and 'mmregmin' are very time consuming!
     parfor i = 1:3
-        if i == 2
-            M{i} = mmregmin(imgFiltered{i}, regionalExtremumSE); % Special case for the saturation component.
-        else
-            M{i} = mmregmax(imgFiltered{i}, regionalExtremumSE);
-        end
+        M{i} = mmregmax(imgFiltered{i}, regionalExtremumSE);
     end
 
     nucleiSE = mmsedisk(nucleiRadius, '2D', 'OCTAGON');
@@ -87,7 +81,7 @@ function [THS, c, nucleiRadius] = DetectionOfParasites(imgRGB)
     muS = funVolume(uint8(MHS) .* imgFiltered{2}) / volumeMHS;
 
     TH = im2bw(imgFiltered{1}, double(muH) / 255);
-    TS = ~im2bw(imgFiltered{2}, double(muS) / 255);    
+    TS = im2bw(imgFiltered{2}, double(muS) / 255);    
     
     THS = TH & TS;
 
@@ -95,6 +89,5 @@ function [THS, c, nucleiRadius] = DetectionOfParasites(imgRGB)
     imwrite(M{2}, '../output/MS.png')
     imwrite(MHS, '../output/MHS.png')
     imwrite(THS, '../output/THS.png')
-
 end
 
index 8c59ace..660463f 100644 (file)
@@ -3,7 +3,7 @@
 %  c: (TODO: should be the average red cell size)
 function [schizontsSmoothed] = DetectionOfSchizonts(THS, c)    
     
-    schizontsConnectivityMin = 3; % Minimum cluster size when detecting schizonts.
+    schizontsConnectivityMin = 6; % Minimum cluster size when detecting schizonts.
     redCellsSE = mmsedisk(c, '2D', 'OCTAGON');    
     schizonts = zeros(size(THS), 'uint8');    
     
@@ -41,7 +41,7 @@ function [schizontsSmoothed] = DetectionOfSchizonts(THS, c)
         end
     end
     
-    %% Fill 'schizonts' by testing the connectivity of each cluster set.
+    %% Fill the variable 'schizonts' by testing the connectivity of each cluster set.
     [~, C] = graphconncomp(adjacencyMatrix, 'Directed', false, 'Weak', true);
     n = 1;    
     while 1
@@ -65,5 +65,5 @@ function [schizontsSmoothed] = DetectionOfSchizonts(THS, c)
     % Dilation + Morphological smoothing (Hole filling).
     schizontsSmoothed = logical(mmclohole(mmdil(schizonts, mmsedisk(4))));
     
-    imwrite(schizontsSmoothed, '../output/schizonts.png')
+    WriteImageGradient(schizontsSmoothed, '../output/schizonts.png', [255, 30, 30]);
 end
\ No newline at end of file
index 5c711b2..3d4055f 100644 (file)
@@ -2,7 +2,7 @@
 %  c: The typical red cell radius [px].
 function [WBCSmoothed] = DetectionOfWhiteCells(THS, c)    
     % We use a radius near the smallest red blood cell size (80% of the radius)
-    redCellsSE = mmsedisk(c * 0.8, '2D', 'OCTAGON'); 
+    redCellsSE = mmsedisk(c * 0.5, '2D', 'OCTAGON'); 
 
     % The white cells may not be correctly marked in the case they have a very special shape (oblong).
     WBCMarker = mmero(THS, redCellsSE); % Erosion to keep only fragments of white cells.
@@ -11,5 +11,5 @@ function [WBCSmoothed] = DetectionOfWhiteCells(THS, c)
     % Morphological smoothing (Hole filling).
     WBCSmoothed = mmclohole(WBC);
     
-    imwrite(WBCSmoothed, '../output/whiteCells.png')
+    WriteImageGradient(WBCSmoothed, '../output/whiteCells.png', [150, 255, 150]);
 end
\ No newline at end of file
index 7a0f69f..294d691 100644 (file)
@@ -1,5 +1,7 @@
 function [infectedRedCells] = IdentificationOfInfectedRedCells(redCells, parasites)
 
     infectedRedCells = mminfrec(parasites, redCells);    
-    imwrite(infectedRedCells, '../output/infected red cells.png')
+    
+    
+    WriteImageGradient(infectedRedCells, '../output/infected red cells.png', [50, 50, 255]);
 end
\ No newline at end of file
index 0456261..d4807dc 100644 (file)
@@ -13,56 +13,65 @@ delete('../output/*');
 % imageFolder = '1405022890'; % 4.
 % imageFolder = '1409191647'; % 9.
 imageFolder = '1412151257'; % 10. (Gamma)
+% imageFolder = '13.05.2015';
 
 %% Gamma
-% imageNumber = 'Gamma-0.8-0002';
-% imageNumber = 'Gamma-0.8-0005';
-% imageNumber = 'Gamma-0.8-0008';
+% imageName = '1412151257-Gamma-0.8-0002.png';
+% imageName = '1412151257-Gamma-0.8-0005.png';
+% imageName = '1412151257-Gamma-0.8-0008.png';
 
-% imageNumber = 'Gamma-1-0001';
-% imageNumber = 'Gamma-1-0004';
-imageNumber = 'Gamma-1-0007';
+% imageName = '1412151257-Gamma-1-0001.png';
+% imageName = '1412151257-Gamma-1-0004.png';
+imageName = '1412151257-Gamma-1-0007.png';
  
-% imageNumber = 'Gamma-1.2-0003';
-% imageNumber = 'Gamma-1.2-0006';
-% imageNumber = 'Gamma-1.2-0009';
+% imageName = 'Gamma-1.2-0003';
+% imageName = 'Gamma-1.2-0006';
+% imageName = 'Gamma-1.2-0009';
+
+%% 13 mai 2015
+
+% imageName = '1412151257-100x-Teinte30-Saturation3-0013.tif';
+% imageName = '1412151257-100x-Teinte30-Saturation3-0014.tif'; % Good result.
+% imageName = '1412151257-100x-Teinte0-Saturation0-0015.tif'; % Bad result.
+
+segmentationMethod = SegmentationMethod.WatershedByMorphologicalGradient;
 
 %% Main
-fprintf('Folder: %s, Number: %s\n', imageFolder, imageNumber);
+fprintf('Folder: %s, Image name: %s\n', imageFolder, imageName);
 
-scaleFactor = 1;
+scaleFactor = 0.6;
 
 % Load the image and its ground truth.
-[imgRGB, gt] = loadImg(imageFolder, imageNumber);
+imgRGB = loadImg(imageFolder, imageName);
 
 % Resample the image and its ground truth.
 imgRGBResampled = imresize(imgRGB, scaleFactor);
-gtResampled = imresize(gt, scaleFactor);
 
 imwrite(imgRGBResampled, '../output/imgRGB.png')
-imwrite(gtResampled, '../output/gt.png')
 
 disp('1) Detection of parasites ...')
 [THS, c, nucleiRadius] = DetectionOfParasites(imgRGBResampled);
 
 disp('2) Detection of white cells ...')
 [WBC] = DetectionOfWhiteCells(THS, c);
+THSWithoutWBC = THS & ~WBC; % We remove the white cells from THS.
 
 disp('3) Detection of schizonts ...')
-[schizonts] = DetectionOfSchizonts(THS, c);
-[imgRGBRedCells] = RemoveWBCAndSchizonts(imgRGBResampled, WBC, schizonts);
+[schizonts] = DetectionOfSchizonts(THSWithoutWBC, c);
+THS(schizonts) = 0; % We remove the white cells from THS.
+THSWithoutWBCandSchizonts = THSWithoutWBC & ~schizonts; % We remove the shizonts from THS.
 
 disp('4) Segmentation of red cells ...')
-[redCells] = SegmentationOfRedCells(imgRGBRedCells, c);
+[redCells] = SegmentationOfRedCells(imgRGBResampled, c, WBC, schizonts, segmentationMethod);
 
 disp('5) Finding infected red cells ...')
-parasites = THS & ~WBC & ~schizonts;
-[infectedRedCells] = IdentificationOfInfectedRedCells(redCells, parasites);
+[infectedRedCells] = IdentificationOfInfectedRedCells(redCells, THSWithoutWBCandSchizonts);
 
+disp('6) Counting red cells and infected red cells ...')
 nbRedCells = mmstats(mmlabel(redCells), 'max');
 nbInfectedRedCells = mmstats(mmlabel(infectedRedCells), 'max');
 infectionPercentage = 100 * nbInfectedRedCells / nbRedCells;
-fprintf('Percentage of infected cell: %.1f%%\n', infectionPercentage)
+fprintf('Percentage of infected cell: %.1f%% (%d/%d)\n', infectionPercentage, nbInfectedRedCells, nbRedCells)
 
 disp('Finished')
 toc
\ No newline at end of file
index ef12a9d..8da0e90 100644 (file)
-function [redCells] = SegmentationOfRedCells(imgRGB, c)
-    cMin = uint32(c * 0.8); % 80%.
+function [segmentedCells] = SegmentationOfRedCells(imgRGB, c, WBC, schizonts, segmentationMethod)
+    cMin = uint32(c * 0.75); % 75%.
     cMinArea = pi * cMin ^ 2;
 
-    gray = mmareaopen(rgb2gray(imgRGB), 2 * cMinArea); % The gray image is better than the green component alone.    
+    gray = imgRGB(:,:,2);
+    gray(WBC) = 255;
+    gray(schizonts) = 255;
     imwrite(gray, '../output/red cells segmentation - gray.png')
     
-    % Thresholding using the Otsu's method .
-    threshold = ~im2bw(gray, graythresh(gray));
+    grayFlat = mmareaopen(gray, 2 * cMinArea, mmsebox(1)); % The gray image is better than the green component alone.    
+    imwrite(grayFlat, '../output/red cells segmentation - gray flat.png')
+    
+    % Thresholding using the Otsu's method.
+    % TODO : montrer l'histogramme
+    threshold = ~im2bw(grayFlat, graythresh(grayFlat));
     
     % We remove all objects with a smaller area than a red blood cell.
     threshold = mmareaopen(threshold, cMinArea);
     
     % We set the background to 0 by using the previous threshold.
-    redCellsFiltered = gray;
+    redCellsFiltered = grayFlat;
     redCellsFiltered(~threshold) = 0;
     
     imwrite(redCellsFiltered, '../output/red cells segmentation - input.png')
     
     %% Segmentation
-    opened = mmopen(redCellsFiltered, mmsedisk(cMin, '2D', 'EUCLIDEAN', 'NON-FLAT', 0));
+    
+    % To smooth the borders and to remove little elements.
+    opened = mmopen(redCellsFiltered, mmsedisk(cMin, '2D', 'EUCLIDEAN', 'NON-FLAT'));
     imwrite(opened, '../output/red cells segmentation - opened.png')
     
-    openedGradient = mmgradm(opened);
-    imwrite(openedGradient, '../output/red cells segmentation - opened - gradient.png')
-        
-    openedGradientThreshold = im2bw(openedGradient, 0.02);
-    imwrite(openedGradientThreshold, '../output/red cells segmentation - opened - gradient - threshold.png')
-    
-    holesClosed = mmareaclose(openedGradientThreshold, 10 * cMinArea);
-    imwrite(holesClosed, '../output/red cells segmentation - holes closed.png')
-        
-    % Watershed algorithm from here: http://mmorph.com/mxmorph/html/mmdemos/mmdcells.html.
-    d = holesClosed;
-    e1 = mmdist(d, mmsebox,'EUCLIDEAN');
-    e3 = mmregmax(e1);
-    e  = mmdil(e3, mmsedisk(5, '2D', 'EUCLIDEAN'));
-    f = mmneg(e1);
-    g = mmcwatershed(f, e, mmsebox);
-    h = mmintersec(d, mmneg(g));
-    i = mmedgeoff(h);
-    j = mmareaopen(i, cMinArea / 3);
-    imwrite(j, '../output/red cells segmentation - individual.png')
-    
-    redCells = j;
+    openedThreshold = im2bw(opened, 0.02);
+    imwrite(openedThreshold, '../output/red cells segmentation - opened - threshold.png')    
+    
+    switch segmentationMethod        
+        case SegmentationMethod.WatershedByDistanceTransform            
+            %% Simplified version.
+    
+            % Remove cells which are touching the edges.
+            segmentedCells = mmedgeoff(WatershedsByDistanceTransform(openedThreshold));
+                    
+        case SegmentationMethod.WatershedByMorphologicalGradient  
+            
+            %% Don't know the meaning. (end of the right column of page 139).
+%             openedGradient = mmgradm(opened);
+%             imwrite(openedGradient, '../output/red cells segmentation - opened - gradient.png')
+
+%             openedGradientThreshold = im2bw(openedGradient, 0.02);
+%             imwrite(openedGradientThreshold, '../output/red cells segmentation - opened - gradient - threshold.png')
+
+%             holesClosed = mmareaclose(openedGradientThreshold, 10 * cMinArea);
+%             imwrite(holesClosed, '../output/red cells segmentation - holes closed.png')
+            
+            % Watershed by distance transform.
+            watershed = mmedgeoff(WatershedsByDistanceTransform(openedThreshold));
+            imwrite(watershed, '../output/red cells segmentation - watershed.png')
+            
+            watershedLabeled = mmlabel(watershed);
+            axisProperties = regionprops(watershedLabeled, 'MajorAxisLength', 'MinorAxisLength');
+            
+            % The area opening remove the little artifacts.
+            watershedSingles = mmareaopen(watershed, cMinArea / 5);
+            
+            for k = 1:length(axisProperties)
+                % We remove each composite cell from 'watershedSingle'.
+                if axisProperties(k).MajorAxisLength / axisProperties(k).MinorAxisLength >= 1.3
+                    watershedSingles(watershedLabeled == k) = 0;
+                end
+            end            
+            imwrite(watershedSingles, '../output/red cells segmentation - watershed singles.png')
+
+            % We remove the single cells from the initial image to obtain the composite cells.
+            composites = mmedgeoff(and(openedThreshold, ~watershedSingles));
+            imwrite(composites, '../output/red cells segmentation - composites.png')
+            
+            compositesOpened = grayFlat;
+            compositesOpened(~composites) = 0;
+            compositesOpened = mmopen(compositesOpened, mmsedisk(cMin, '2D', 'EUCLIDEAN', 'FLAT'));
+            
+            imwrite(compositesOpened, '../output/red cells segmentation - composites opened.png')
+            
+            compositesOpenedGradient = mmgradm(compositesOpened);
+            imwrite(compositesOpenedGradient, '../output/red cells segmentaion - composites opened gradient.png');
+            
+            compositesOpenedGradientThreshold = im2bw(compositesOpenedGradient, 0.00001);
+            imwrite(compositesOpenedGradientThreshold, '../output/red cells segmentation - composites opened gradient threshold.png')
+            
+            compositesOpenedGradientThresholdHolesFilled = ~mmareaopen(~compositesOpenedGradientThreshold, cMinArea / 2);
+            imwrite(compositesOpenedGradientThresholdHolesFilled, '../output/red cells segmentation - composites opened gradient threshold holes filled.png')
+            
+            compositesThinning = mmthin(compositesOpenedGradientThresholdHolesFilled);
+            compositesThinningFlooded = mmareaclose(compositesThinning, 5 * cMinArea);
+            compositesThinningFlooded(compositesThinning) = 0;
+            compositesThinningFlooded = mmero(compositesThinningFlooded); % We erode by a 3x3 cross to avoid the composites to touch the singles when fusionning later.
+            imwrite(compositesThinningFlooded, '../output/red cells segmentation - composites thinning.png')
+            
+            % Compute the ratio of the major axis over the minor axis and keep
+            % cells with a ration greater than or equal to 1.3.
+            segmentedCells = or(compositesThinningFlooded, watershedSingles);            
+    end
+    
+    WriteImageGradient(segmentedCells, '../output/red cells segmentation - individual.png', [255, 255, 255]);
 end
\ No newline at end of file
index 996f816..df4f05f 100644 (file)
@@ -7,7 +7,7 @@ close all
 %% Data
 % RGB = imread ('/Users/michel.kocher/michel/Data/Movable/imgs_corrMK/Gamma/1412151257-Gamma-0.8-0002.png');
 
-RGB = imread('../imgs_corrMK/Gamma/1412151257-Gamma-0.8-0002.png');
+RGB = imread('../imgs_corrMK/1412151257/1412151257-Gamma-0.8-0002.png');
 
 %% HSV
 HSV = rgb2hsv (RGB);
@@ -44,4 +44,4 @@ imshow (S, []), title ('Saturation')
 linkaxes (ax)
 
 %% Print 2 file
-print ('-f1', '-dpng', 'THue.png')
\ No newline at end of file
+print ('-f1', '-dpng', 'THue.png')
\ No newline at end of file
index c985597..14e60b7 100644 (file)
@@ -1,20 +1,7 @@
-% Load the image and the ground-truth from the given foldername and suffix.
-% If the ground-truth isn't found then it is put to white.
-function [img, gt] = loadImg(name, num)
-
+% Load the image from the given foldername and filename.
+function img = loadImg(dirname, filename)
     imagesFolder = '../imgs_corrMK';
-    groundTruthFolder = '../gt';
-    filename = [name, '-', num, '.png'];
-    filepath_img = [imagesFolder, '/', name, '/', filename];
-    filepath_gt = [groundTruthFolder, '/', name, '/', filename];
-    
-    img = imread(filepath_img);
-    
-    try
-        gt = imread(filepath_gt);
-    catch exception
-        fprintf('Can''t load the ground-truth : %s\n', exception.message)
-        gt = ones(size(img));
-    end        
+    filepath_img = [imagesFolder, '/', dirname, '/', filename];    
+    img = imread(filepath_img);  
 end