From 1a7a8e6286ceff00c6e1edc22c41d1599c9a6bdb Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Sat, 16 May 2015 16:27:22 +0200 Subject: [PATCH] =?utf8?q?Application=20de=20la=20s=C3=A9gmentation=20prop?= =?utf8?q?os=C3=A9e=20par=20l'article=20+=20plein=20d'autres=20changements?= =?utf8?q?.?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit --- src/DetectionOfParasites.m | 31 +++---- src/DetectionOfSchizonts.m | 6 +- src/DetectionOfWhiteCells.m | 4 +- src/IdentificationOfInfectedRedCells.m | 4 +- src/Main.m | 49 ++++++----- src/SegmentationOfRedCells.m | 115 ++++++++++++++++++------- src/THue.m | 4 +- src/loadImg.m | 21 +---- 8 files changed, 141 insertions(+), 93 deletions(-) diff --git a/src/DetectionOfParasites.m b/src/DetectionOfParasites.m index ab7ffc7..c24de73 100644 --- a/src/DetectionOfParasites.m +++ b/src/DetectionOfParasites.m @@ -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 diff --git a/src/DetectionOfSchizonts.m b/src/DetectionOfSchizonts.m index 8c59ace..660463f 100644 --- a/src/DetectionOfSchizonts.m +++ b/src/DetectionOfSchizonts.m @@ -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 diff --git a/src/DetectionOfWhiteCells.m b/src/DetectionOfWhiteCells.m index 5c711b2..3d4055f 100644 --- a/src/DetectionOfWhiteCells.m +++ b/src/DetectionOfWhiteCells.m @@ -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 diff --git a/src/IdentificationOfInfectedRedCells.m b/src/IdentificationOfInfectedRedCells.m index 7a0f69f..294d691 100644 --- a/src/IdentificationOfInfectedRedCells.m +++ b/src/IdentificationOfInfectedRedCells.m @@ -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 diff --git a/src/Main.m b/src/Main.m index 0456261..d4807dc 100644 --- a/src/Main.m +++ b/src/Main.m @@ -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 diff --git a/src/SegmentationOfRedCells.m b/src/SegmentationOfRedCells.m index ef12a9d..8da0e90 100644 --- a/src/SegmentationOfRedCells.m +++ b/src/SegmentationOfRedCells.m @@ -1,46 +1,103 @@ -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 diff --git a/src/THue.m b/src/THue.m index 996f816..df4f05f 100644 --- a/src/THue.m +++ b/src/THue.m @@ -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 diff --git a/src/loadImg.m b/src/loadImg.m index c985597..14e60b7 100644 --- a/src/loadImg.m +++ b/src/loadImg.m @@ -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 -- 2.45.2