From ad7e0b7a0c5428532a9b0be6d8c9d0ad0ecffd4c Mon Sep 17 00:00:00 2001 From: moonanjum26 Date: Wed, 22 Aug 2018 22:10:26 +0530 Subject: [PATCH] malaria detection code in matlab --- src/gopa kumar code/CNNTRainFczRGBCall.m | 9 + src/gopa kumar code/CnnSvmOnSlideCalling.m | 42 ++ .../MalEx6ColorRewrittenTrainFcZRGB.m | 431 ++++++++++++ src/gopa kumar code/Read Me.txt | 25 + src/gopa kumar code/SVMTrainBFcsRGB.m | 48 ++ src/gopa kumar code/countTPFNFPTN.m | 74 ++ .../generate32X32FcsStck_Fcsd_Feat_RGB.m | 289 ++++++++ .../generateDstAndWBCPatches.m | 282 ++++++++ src/gopa kumar code/getSegmentation.m | 649 ++++++++++++++++++ .../initializeCharacterCNNBNNEx6ColorFczRGB.m | 37 + .../makeDataSetForTrainingBest.m | 310 +++++++++ src/gopa kumar code/reName.m | 30 + src/gopa kumar code/setup.m | 25 + src/gopa kumar code/testCNNBNNEx6Color.m | 13 + .../testTrainedCNNSVMAutoCnt.m | 582 ++++++++++++++++ 15 files changed, 2846 insertions(+) create mode 100644 src/gopa kumar code/CNNTRainFczRGBCall.m create mode 100644 src/gopa kumar code/CnnSvmOnSlideCalling.m create mode 100644 src/gopa kumar code/MalEx6ColorRewrittenTrainFcZRGB.m create mode 100644 src/gopa kumar code/Read Me.txt create mode 100644 src/gopa kumar code/SVMTrainBFcsRGB.m create mode 100644 src/gopa kumar code/countTPFNFPTN.m create mode 100644 src/gopa kumar code/generate32X32FcsStck_Fcsd_Feat_RGB.m create mode 100644 src/gopa kumar code/generateDstAndWBCPatches.m create mode 100644 src/gopa kumar code/getSegmentation.m create mode 100644 src/gopa kumar code/initializeCharacterCNNBNNEx6ColorFczRGB.m create mode 100644 src/gopa kumar code/makeDataSetForTrainingBest.m create mode 100644 src/gopa kumar code/reName.m create mode 100644 src/gopa kumar code/setup.m create mode 100644 src/gopa kumar code/testCNNBNNEx6Color.m create mode 100644 src/gopa kumar code/testTrainedCNNSVMAutoCnt.m diff --git a/src/gopa kumar code/CNNTRainFczRGBCall.m b/src/gopa kumar code/CNNTRainFczRGBCall.m new file mode 100644 index 0000000..87029de --- /dev/null +++ b/src/gopa kumar code/CNNTRainFczRGBCall.m @@ -0,0 +1,9 @@ +function CNNTRainFczRGBCall() + DataSetLoc = 'WithDstWBC\'; + DataSet = 'FczRGB_MalImdb'; MalEx6ColorRewrittenTrainFcZRGB(DataSetLoc, DataSet); + DataSetLoc = 'WithOutDstWBC\'; + DataSet = 'FczRGB_MalImdb'; MalEx6ColorRewrittenTrainFcZRGB(DataSetLoc, DataSet); + DataSetLoc = 'WithWBCNoDst\'; + DataSet = 'FczRGB_MalImdb'; MalEx6ColorRewrittenTrainFcZRGB(DataSetLoc, DataSet); +% CnnSvmOnSlideCalling(); +end \ No newline at end of file diff --git a/src/gopa kumar code/CnnSvmOnSlideCalling.m b/src/gopa kumar code/CnnSvmOnSlideCalling.m new file mode 100644 index 0000000..03a8208 --- /dev/null +++ b/src/gopa kumar code/CnnSvmOnSlideCalling.m @@ -0,0 +1,42 @@ +function CnnSvmOnSlideCalling() +close all; clear all; clc; + %If 3rd param is true (SVMModel), second arg has no significance + CNet (:, :, 1) = [20 20;38 26;30 34]; + CNet (:, :, 2) = [36 40;38 32;33 37]; + params.alreadySegmented = false; + for NetExpt = 1:1 + NetIds = CNet (:, :, NetExpt); + for Mdls = 1:1 + if (Mdls == 1) + params.UseMdlFrm = 'WithDstWBC'; + elseif (Mdls == 2) + params.UseMdlFrm = 'WithOutDstWBC'; + else + params.UseMdlFrm = 'WithWBCNoDst'; + end + + for CNNSVM = 2:2 + if (CNNSVM == 1) + params.isSVMModel = false; + for subModel = 1:2 + params.NetId = NetIds(Mdls, subModel); + if (subModel == 1) + params.dataSet = 'BFczdRGB'; + else + params.dataSet = 'FczRGB'; + end + display (['Mdls = ' params.UseMdlFrm '; isSVM = ' num2str(params.isSVMModel) '; SubMOdel = ' params.dataSet]); + testTrainedCNNSVMAutoCnt(params); + end + else + %No need for different NetIds %continue; + params.loadFeat = false; + params.isSVMModel = true; + params.dataSet = 'BFczdFeatRGB'; + display (['Mdls = ' params.UseMdlFrm '; isSVM = ' num2str(params.isSVMModel) '; SubMOdel = ' params.dataSet]); + testTrainedCNNSVMAutoCnt(params); + end + end + end + end +end \ No newline at end of file diff --git a/src/gopa kumar code/MalEx6ColorRewrittenTrainFcZRGB.m b/src/gopa kumar code/MalEx6ColorRewrittenTrainFcZRGB.m new file mode 100644 index 0000000..0d7b404 --- /dev/null +++ b/src/gopa kumar code/MalEx6ColorRewrittenTrainFcZRGB.m @@ -0,0 +1,431 @@ +function MalEx6ColorRewrittenTrainFcZRGB(DataSetLoc, bkUpName, varargin) + tileGenFlag = true; testFlag = false; ConfnMatxs = []; + clc; close all; + %initialize + setup ; + %Create the network + net = initializeCharacterCNNBNNEx6ColorFczRGB() ; + + maxEpoch = 100; + + %Train Options + trainOpts.batchSize = 100 ; + trainOpts.numEpochs = maxEpoch; + trainOpts.continue = true ; + trainOpts.useGpu = false ; + trainOpts.learningRate = 0.001 ; + folderOfIntst = ['myExp6\' DataSetLoc bkUpName]; + mkdir(folderOfIntst); + trainOpts.expDir = folderOfIntst ; + delete([folderOfIntst '\*mat']); + trainOpts.errorType = 'multiclass'; %'binary' ; + trainOpts = vl_argparse(trainOpts, varargin); + + % Take the average image out + imdb = getPreLoadedDataset(DataSetLoc, bkUpName);%PercPstvTrainVaidnTest, PercNgtvTrainVaidnTest);%populateANewDataSetAB(); % + imageMean = mean(imdb.images.data(:)) ; + save([folderOfIntst '\imageMean'], 'imageMean'); + imdb.images.data = imdb.images.data - imageMean ; + + % Call training function in MatConvNet + [net, info] = cnn_train(net, imdb, @getBatch, trainOpts); + + % Save the result for later use + net.layers(end) = [] ; + net.imageMean = imageMean ; + save([folderOfIntst '\latest.mat'], '-struct', 'net'); + + %on training set. As we are testing and inside testing code we are + %subtracting the mean, we add ir back + imdb.images.data = imdb.images.data + imageMean ; + trainIndxes = (imdb.images.set == 1); + trainImages = imdb.images.data(:, :, :, trainIndxes); + trainLbls = imdb.images.label(trainIndxes); + + trainImdb.images.id = 1:length(trainLbls); + trainImdb.images.label = trainLbls; + trainImdb.images.data = trainImages; + Lbls = testCNNBNNEx6Color(trainImdb, net); + [TnConfnMat, TnFPTile, TnFNTile] = getConfnMatrix(Lbls, trainImdb, tileGenFlag); + clear trainImdb; + + %Validation stat + validnIndxes = (imdb.images.set == 2); + validnImages = imdb.images.data(:, :, :, validnIndxes); + validnLbls = imdb.images.label(validnIndxes); + + validnImdb.images.id = 1:length(validnLbls); + validnImdb.images.label = validnLbls; + validnImdb.images.data = validnImages; + Lbls = testCNNBNNEx6Color(validnImdb, net); + [VnConfnMat, VnFPTile, VnFNTile] = getConfnMatrix(Lbls, validnImdb, tileGenFlag); + clear validnImdb; + + if (testFlag) + %on testing set + testIndxes = (imdb.images.set == 3); + testImages = imdb.images.data(:, :, :, testIndxes); + testLbls = imdb.images.label(testIndxes); + + testImdb.images.id = 1:length(testLbls); + testImdb.images.label = testLbls; + testImdb.images.data = testImages; + Lbls = testCNNBNNEx6Color(testImdb, net); + [TtConfnMat, TtFPTile, TtFNTile] = getConfnMatrix(Lbls, testImdb, tileGenFlag); + clear testImdb; + if (~isempty(TtFPTile)) + TtNFP = (TtFPTile - min(TtFPTile(:)))/(max(TtFPTile(:)) - min(TtFPTile(:))); + imwrite (TtNFP, [folderOfIntst '\TtNFP' num2str(maxEpoch) '.jpg']); + end + if (~isempty(TtFNTile)) + TtNFN = (TtFNTile - min(TtFNTile(:)))/(max(TtFNTile(:)) - min(TtFNTile(:))); + imwrite (TtNFN, [folderOfIntst '\TtNFN' num2str(maxEpoch) '.jpg']); + end + save ([folderOfIntst '\TtConfnMat'], 'TtConfnMat'); + ConfnMatxs = TtConfnMat; + end + clear FullImdb; + + if (~isempty(TnFPTile)) + TnNFP = (TnFPTile - min(TnFPTile(:)))/(max(TnFPTile(:)) - min(TnFPTile(:))); + imwrite (TnNFP, [folderOfIntst '\TnNFP' num2str(maxEpoch) '.jpg']); + end + if (~isempty(TnFNTile)) + TnNFN = (TnFNTile - min(TnFNTile(:)))/(max(TnFNTile(:)) - min(TnFNTile(:))); + imwrite (TnNFN, [folderOfIntst '\TnNFN' num2str(maxEpoch) '.jpg']); + end + + if (~isempty(VnFPTile)) + VnNFP = (VnFPTile - min(VnFPTile(:)))/(max(VnFPTile(:)) - min(VnFPTile(:))); + imwrite (VnNFP, [folderOfIntst '\VnNFP' num2str(maxEpoch) '.jpg']); + end + if (~isempty(VnFNTile)) + VnNFN = (VnFNTile - min(VnFNTile(:)))/(max(VnFNTile(:)) - min(VnFNTile(:))); + imwrite (VnNFN, [folderOfIntst '\VnNFN' num2str(maxEpoch) '.jpg']); + end + + + + save ([folderOfIntst '\TnConfnMat'], 'TnConfnMat'); + save ([folderOfIntst '\VnConfnMat'], 'VnConfnMat'); + + ConfnMatxs = [TnConfnMat ;VnConfnMat; ConfnMatxs]; + [Sensitivity, Specificity, FScore] = getStatistics(ConfnMatxs) +end + +function [Sensitivity, Specificity, FScore] = getStatistics(ConfnMatxs) + Sensitivity = zeros(2, 1); Specificity = zeros(2, 1); FScore = zeros(2, 1); + for i = 1:2 + CConfnMat = ConfnMatxs((i-1)*2+1:i*2, 1:2); + Sensitivity(i) = CConfnMat(1, 1)/(CConfnMat(1, 1) + CConfnMat(1, 2)); + Specificity(i) = CConfnMat(2, 2)/(CConfnMat(2, 2) + CConfnMat(2, 1)); + FScore(i) = 2*CConfnMat(1, 1)/(2*CConfnMat(1, 1)+ CConfnMat(1, 2)+CConfnMat(2, 1)); + end +end +function bkUpName = getBackUpFolderName(PercPstvTrainVaidnTest, PercNgtvTrainVaidnTest) +% PercPstvTrainVaidnTest = [1 45]; PercNgtvTrainVaidnTest = [10.3 4.5]; + strname = num2str([PercPstvTrainVaidnTest PercNgtvTrainVaidnTest]); + bkUpName = ''; flag = 0; + for i = 1:length(strname) + if (flag == 1) + if (strname(i) ~= ' ') + flag = 0; + else + continue; + end + end + if (strname(i) == '.') + flag = 0; + bkUpName(end+1) = 'P'; + elseif (strname(i) == ' ') + bkUpName(end+1) = '_'; + flag = 1; + else + flag = 0; + bkUpName(end+1) = strname(i); + end + end +end +function [ConfnMat, FPTile, FNTile] = getConfnMatrix(Lbls, imdb, genTile) + FPTile = []; FNTile = []; + ActlLbls = Lbls(2, :); ObtndLbls = Lbls(1, :); + numCls = max(ActlLbls); + ConfnMat = zeros(numCls, numCls); + for i = 1:numCls + lookngFr = (ActlLbls == i); + for j = 1:numCls + ConfnMat(i, j) = sum(ObtndLbls(lookngFr) == j); + if (i ~= j) + wrngIndx = find(lookngFr & (ObtndLbls == j)); + FlseImgStck = imdb.images.data(:, :, :, wrngIndx); + if (genTile && ~isempty(FlseImgStck)) + if (i == 1) + FNTile = generateTileFromStack(FlseImgStck); + else + FPTile = generateTileFromStack(FlseImgStck); + end + end + end + end + end +end +function imdbN = populateANewDataSetAB() + imdb = load('data/charsdb.mat') ; + %Select all As and Bs + imdbN = imdb; + numImages = length(imdb.images.id); + idN = []; dataN = single(zeros(32, 32, 1862)); labelN = []; setN = []; cnt = 1; + for i = 1:numImages + cLabel = imdb.images.label(i); + if (cLabel == 1 || cLabel == 2) + idN = [idN cnt]; + dataN(:, :, cnt) = imdb.images.data(:, :, cnt); + labelN = [labelN cLabel]; + setN = [setN imdb.images.set(i)]; + cnt = cnt+1; + end + end + imdbN.images.id = idN; + imdbN.images.data = dataN; + imdbN.images.label = labelN; + imdbN.images.set = setN; +end + +% -------------------------------------------------------------------- +function [im, labels] = getBatch(imdb, batch) +% -------------------------------------------------------------------- +im = imdb.images.data(:,:,:, batch) ; +im = 256 * reshape(im, 32, 32, 9, []) ; +labels = imdb.images.label(1,batch) ; +end +% -------------------------------------------------------------------- +function [im, labels] = getBatchWithJitter(imdb, batch) + % -------------------------------------------------------------------- + im = imdb.images.data(:,:,batch) ; + labels = imdb.images.label(1,batch) ; + + n = numel(batch) ; + train = find(imdb.images.set == 1) ; + + sel = randperm(numel(train), n) ; + im1 = imdb.images.data(:,:,sel) ; + + sel = randperm(numel(train), n) ; + im2 = imdb.images.data(:,:,sel) ; + + ctx = [im1 im2] ; + ctx(:,17:48,:) = min(ctx(:,17:48,:), im) ; + + dx = randi(11) - 6 ; + im = ctx(:,(17:48)+dx,:) ; + sx = (17:48) + dx ; + + dy = randi(5) - 2 ; + sy = max(1, min(32, (1:32) + dy)) ; + + im = ctx(sy,sx,:) ; + + % Visualize the batch: + % figure(100) ; clf ; + % vl_imarraysc(im) ; + + im = 256 * reshape(im, 32, 32, 1, []) ; +end +function MalImdb = getMyDataset(PercPstvTrainVaidnTest, PercNgtvTrainVaidnTest) + firstRun = false; writeFlag = false; + DifCultDataSet = true; + if (DifCultDataSet) + load DifSamSlctdByPgm; + load ('E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\Patch32By32\DataSet32By32Color\postvStack_'); + + + ImgStack = DifSamSlctdByPgm.ImgsDifcltStack; + TrthStack = DifSamSlctdByPgm.TrthDifcltStack; + LblsStack = DifSamSlctdByPgm.LblsDifcltStack; + TrthLbl = uint8(zeros(size(LblsStack))); + [SzM, SzN, SzO, numImgs] = size(ImgStack); + PtchSzM = 32; PtchSzN = 32; ofst = 7; + for i = 1:numImgs + TthMskPatch = TrthStack(:, :, i); + cntrPatch = TthMskPatch(PtchSzM/2-ofst:PtchSzM/2+ofst, PtchSzN/2-ofst:PtchSzN/2+ofst); + cntCentr = sum(cntrPatch(:)); + if (cntCentr == 1) + TrthLbl(i) = 1; + else + TrthLbl(i) = 2; + end + end + + numPstvImgs = sum(TrthLbl == 1); + numNgtvImgs = numImgs - numPstvImgs; + negtvStack = ImgStack(:, :, :, TrthLbl == 2); + postvStack = augmentDataSetByRotation(postvStack); + [~, ~, ~, numPstvImgs] = size(postvStack); + + %Create a full Databasclear alle + FullImdb.meta.classes = 'Malaria,Healthy'; + FullImdb.meta.sets = {'train', 'val', 'test'}; + FullImdb.meta.infn = 'Difficult Cases For Second CNN'; + FullImdb.images.id = 1:(numPstvImgs+numNgtvImgs); + Fdata = postvStack; + clear postvStack; + Fdata(:, :, :, numPstvImgs+1:(numPstvImgs+numNgtvImgs)) = negtvStack; + clear negtvStack; + FullImdb.images.data = Fdata; + clear Fdata; + FullImdb.images.label = [ones(1, numPstvImgs) 2*ones(1, numNgtvImgs)]; + + pstvSplit = divideDataInRatio(numPstvImgs, PercPstvTrainVaidnTest); + ngtvSplit = divideDataInRatio(numNgtvImgs, PercNgtvTrainVaidnTest); + FullImdb.images.set = [pstvSplit ngtvSplit]; + save('FullImdb', 'FullImdb'); + %Select dataset for training + FllData = FullImdb.images.data; + FllLbl = FullImdb.images.label; + lbelCatgry = FullImdb.images.set; + slectdFrTraining = (lbelCatgry == 1) | (lbelCatgry == 2); + + MalImdb.meta.classes = 'Malaria,Healthy'; + MalImdb.meta.sets = {'train', 'val'}; + + MalImdb.images.id = 1:sum(slectdFrTraining); + MalImdb.images.data = FllData(:, :, :, slectdFrTraining); + MalImdb.images.label = FllLbl(slectdFrTraining); + MalImdb.images.set = lbelCatgry(slectdFrTraining); + else + if (firstRun) + load ('E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\MalResize40By45\negtvStackSngleColor'); + load ('E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\MalResize40By45\postvStackSngleColor'); + postvStack = augmentDataSetByRotation(postvStack); + [~, ~, ~, numPstvImgs] = size(postvStack); + [~, ~, ~, numNgtvImgs] = size(negtvStack); + + + %Create a full Databasclear alle + FullImdb.meta.classes = 'Malaria,Healthy'; + FullImdb.meta.sets = {'train', 'val', 'test'}; + FullImdb.meta.infn = 'Cell As Whole Placed on Bgnd Intnsty 200'; + FullImdb.images.id = 1:(numPstvImgs+numNgtvImgs); + Fdata = postvStack; + clear postvStack; + Fdata(:, :, :, numPstvImgs+1:(numPstvImgs+numNgtvImgs)) = negtvStack; + clear negtvStack; + FullImdb.images.data = Fdata; + clear Fdata; + FullImdb.images.label = [ones(1, numPstvImgs) 2*ones(1, numNgtvImgs)]; + save('E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\MalResize40By45\FullImdb', 'FullImdb'); + % save('E:\Gopakumar\GopakumarIISTDrive\Dataset\4th IIST Visit\Mal01\Image Stack\FullImdb', 'FullImdb'); + end + + load ('E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\MalResize40By45\FullImdb'); + numPstvImgs = sum(FullImdb.images.label == 1); + numNgtvImgs = sum(FullImdb.images.label == 2); + + pstvSplit = divideDataInRatio(numPstvImgs, PercPstvTrainVaidnTest); + ngtvSplit = divideDataInRatio(numNgtvImgs, PercNgtvTrainVaidnTest); + FullImdb.images.set = [pstvSplit ngtvSplit]; + save('E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\MalResize40By45\FullImdb', 'FullImdb'); + + %Select dataset for training + FllData = FullImdb.images.data; + FllLbl = FullImdb.images.label; + lbelCatgry = FullImdb.images.set; + slectdFrTraining = (lbelCatgry == 1) | (lbelCatgry == 2); + + MalImdb.meta.classes = 'Malaria,Healthy'; + MalImdb.meta.sets = {'train', 'val'}; + + MalImdb.images.id = 1:sum(slectdFrTraining); + MalImdb.images.data = FllData(:, :, :, slectdFrTraining); + MalImdb.images.label = FllLbl(slectdFrTraining); + MalImdb.images.set = lbelCatgry(slectdFrTraining); + if (writeFlag) + % load imdbS; + % postvStack = imdb.images.data(:, :, 1:1440); + % selNegtiveStack = imdb.images.data(:, :, 1441:end); + + pstvTiledImg = generateTileFromStack(postvStack); + ngtvTiledImg = generateTileFromStack(selNegtiveStack); + figure; imshow(pstvTiledImg, []); + figure; imshow(ngtvTiledImg, []); + end + end +end + +function augStack = augmentDataSetByRotation(postvStack) + [SzM, SzN, SzO, numImgs] = size(postvStack); + augStack = single(zeros(SzM, SzN, SzO, 4*numImgs)); + augCnt = 0; + for i = 1:numImgs + currIm = postvStack(:, :, :, i); + augStack(:, :, :, augCnt+1) = currIm; %0 degree + augStack(:, :, :, augCnt+2) = rot90(currIm, 1); %90 + augStack(:, :, :, augCnt+3) = rot90(currIm, 2); %180 + augStack(:, :, :, augCnt+4) = rot90(currIm, 3); %270 + augCnt = augCnt + 4; + end +end +function augStack = augmentDataSetByRotationColor(postvStack) + [SzM, SzN, SzO, numImgs] = size(postvStack); + augStack = single(zeros(SzM, SzN, SzO, 4*numImgs)); + augCnt = 0; + for i = 1:numImgs + currIm = postvStack(:, :, :, i); + augStack(:, :, :, augCnt+1) = currIm; %0 degree + augStack(:, :, :, augCnt+2) = rot90(currIm, 1); %90 + augStack(:, :, :, augCnt+3) = rot90(currIm, 2); %180 + augStack(:, :, :, augCnt+4) = rot90(currIm, 3); %270 + augCnt = augCnt + 4; + end +end +function TnTtVn = divideDataInRatio(MaxLt, Ratio) + blck = randperm(MaxLt); strt = 1; + cmSum = cumsum(Ratio); + TnTtVn = zeros(1, MaxLt); + for i = 1:length(Ratio) + intstd = round(cmSum(i)/100*MaxLt); + TnTtVn((blck >= strt) & (blck <= intstd)) = i; + strt = intstd+1; + end +end +function currStackTile = generateTileFromStack(cellStack) + [Szm, Szn, SzO, numCells] = size(cellStack); + numCellsPrRow = round(sqrt(numCells)); + numCellsPrCol = ceil(numCells/numCellsPrRow); + rwGap = 1; clGap = 1; %pixels + TileCl = (numCellsPrRow - 1)*clGap + numCellsPrRow*Szn; + TileRw = (numCellsPrCol - 1)*rwGap + numCellsPrCol*Szm; + currStackTile = zeros(TileRw, TileCl); + startRw = 1; cntCells = 1; + for i = 1:numCellsPrCol+1 + startCl = 1; + endRw = startRw + Szm - 1; + for j = 1:numCellsPrRow + if (cntCells > numCells) + return; + end + endCl = startCl + Szn - 1; + if (SzO == 9) + bndRed = 3; + currCell = cellStack(:, :, 4:6, cntCells); + elseif (SzO == 3) + bndRed = 3; + currCell = cellStack(:, :, :, cntCells); + else + bndRed = 1; + currCell = cellStack(:, :, cntCells); + end + cntCells = cntCells+1; + currStackTile(startRw:endRw, startCl:endCl, 1:bndRed) = currCell; + startCl = endCl + clGap + 1; + end + startRw = endRw + rwGap+1; + end +end +function imdb = getPreLoadedDataset(DataSetLoc, bkUpName) + load (['CNNPatches\' DataSetLoc bkUpName]); + imdb = MalImdb; +end +% load ('E:\Gopakumar\GopakumarIISTDrive\Dataset\4th IIST Visit\Mal01\Image Stack\FullImdb'); +% SplitLbels = FullImdb.images.set; +% save ('E:\Gopakumar\Dataset IIST RSDAY\IIST RSDAY MY WORK\Work\Code\matconvnet-1.0-beta11\practical-cnn-2015a\practical-cnn-2015a\myExp6\lastStat\SplitLbels', 'SplitLbels'); \ No newline at end of file diff --git a/src/gopa kumar code/Read Me.txt b/src/gopa kumar code/Read Me.txt new file mode 100644 index 0000000..e325566 --- /dev/null +++ b/src/gopa kumar code/Read Me.txt @@ -0,0 +1,25 @@ +1. Run SVMTrainBFcsRGB to train SVM + +2. Run CnnSvmOnSlideCalling to test the classifier on each slide image. + +3. Output will be stored in ./PrintPath\WithDstWBC\IndVid\BFczdFeatRGB + +4. Confusion matrices and Total statistics will be stored in ./Programs\myExp6\WithDstWBC\BFczdFeatRGB_MalImdb + + + +%%%%%%%%%%%%%%%%% + + +If base path is 'E:\Gopakumar\GopakumarIISTDrive\Dataset\JCMIG_Dataset_4th IISC\Mal01\Images\GndTruth\Separate\NewEasy\' + +This program access +the best focussed images from 'base path/ImgesMinAcrsStck' + +the focusstack from 'E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\Patch32By32\DataSet32By32Color\Manual_But_Sel_Training\AllFcsStack\' + +DstLocByPGM from 'base path/DstLocnsPgm' + +GndTruths are accessed from 'base path/GndTrth' + + diff --git a/src/gopa kumar code/SVMTrainBFcsRGB.m b/src/gopa kumar code/SVMTrainBFcsRGB.m new file mode 100644 index 0000000..3e64071 --- /dev/null +++ b/src/gopa kumar code/SVMTrainBFcsRGB.m @@ -0,0 +1,48 @@ +function SVMTrainBFcsRGB() + DataSetLoc = 'WithDstWBC\'; saveTrainedModel (DataSetLoc); + if (false) + DataSetLoc = 'WithOutDstWBC\'; saveTrainedModel (DataSetLoc); + DataSetLoc = 'WithWBCNoDst\'; saveTrainedModel (DataSetLoc); + end +end +function saveTrainedModel (DataSetLoc) + imdb = getPreLoadedDataset(DataSetLoc); + setCategory = imdb.images.set; + tnIndx = (setCategory == 1); + vnIndx = (setCategory == 2); + %By default it normalizes + trainData = imdb.images.data(tnIndx, :); + vnData = imdb.images.data(vnIndx, :); + trainLabl = imdb.images.label(tnIndx); + vnnLabl = imdb.images.label(vnIndx); + + sigMa = 0.6; + SvMModel.Meta.Kf = 'rbf'; + if (strcmp(SvMModel.Meta.Kf, 'rbf')) + SvMModel.Meta.Sigma = sigMa; + svmModel = svmtrain(trainData, trainLabl, 'kernel_function', SvMModel.Meta.Kf, 'rbf_sigma', SvMModel.Meta.Sigma); + else + svmModel = svmtrain(trainData, trainLabl); + end + clfdLabl = svmclassify(svmModel, vnData); + VnConfnMat = getConfnMat(clfdLabl, vnnLabl); + VnConfnMat + %0.2813; + SvMModel.model = svmModel; + save (['myExp6\' DataSetLoc 'BFczdFeatRGB_MalImdb\VnConfnMat'], 'VnConfnMat'); +% save (['myExp6\' DataSetLoc 'BFczdFeatRGB_MalImdb\SvMModel_' SvMModel.Meta.Kf], 'SvMModel'); + save (['myExp6\' DataSetLoc 'BFczdFeatRGB_MalImdb\SvMModel'], 'SvMModel'); +end +function confnMat = getConfnMat(clfdLabl, vnnLabl) + confnMat = zeros(2, 2); + for i = 1:2 + curIntst = (vnnLabl == i); + for j = 1:2 + confnMat(i, j) = sum(clfdLabl(curIntst) == j); + end + end +end +function imdb = getPreLoadedDataset(DataSetLoc) + load (['CNNPatches\' DataSetLoc 'BFczdFeatRGB_MalImdb']); + imdb = MalImdb; +end \ No newline at end of file diff --git a/src/gopa kumar code/countTPFNFPTN.m b/src/gopa kumar code/countTPFNFPTN.m new file mode 100644 index 0000000..e8e1ca9 --- /dev/null +++ b/src/gopa kumar code/countTPFNFPTN.m @@ -0,0 +1,74 @@ +function [TP, FN, FP, TN] = countTPFNFPTN(segmentedImg, GndTruth, Classfd) +% imRGB = imread(['E:\Gopakumar\GopakumarIISTDrive\Dataset\4th IIST Visit\Mal01\Images\GndTruth\Separate\NewEasy\ImgesMinAcrsStck\' '0_5_B.jpg']); +% currImgUint8 = imRGB; +% load Classfd; load GndTruth; + ofstSeg = 7; + [X, Y ] = meshgrid(-ofstSeg:ofstSeg, -ofstSeg:ofstSeg); distM = sqrt(X.^2+Y.^2); +% [segmentedImg, segWithBndry, WBCs] = getSegmentation(currImgUint8); + GndTrthClLbl = segmentedImg(GndTruth); + excludeGThCnt = 0; + %sepecial case and is due to inaccurate segmentation + if (sum(GndTrthClLbl == 0) > 0) + %select the closest label from 15x15 neighbourhood. + %for each such points + [r, c] = find(GndTruth); + for i = 1:length(r) + if (segmentedImg(r(i), c(i)) == 0) + candLbls = segmentedImg(r(i)-ofstSeg:r(i)+ofstSeg, c(i)-ofstSeg:c(i)+ofstSeg); + candLblDst = distM .* (candLbls > 0); + candLblDst(candLblDst == 0) = Inf; + candLblDst = (candLblDst == min(candLblDst(:))); + closestLbl = max(candLbls(candLblDst)); %believing that there is a cell in 15x15 loc + if (closestLbl > 0) + GndTrthClLbl = [GndTrthClLbl; closestLbl]; + else + %No cell in the visinity + excludeGThCnt = excludeGThCnt + 1; + end + end + end + end + GndTrthClLbl = unique(GndTrthClLbl); + if (~isempty(GndTrthClLbl) && GndTrthClLbl(1) == 0); GndTrthClLbl(1) = []; end; + ClassfdClLbl = segmentedImg(Classfd); + excludeClsfdCnt = 0; + %sepecial case and is due to inaccurate segmentation + if (sum(ClassfdClLbl == 0) > 0) + %select the closest label from 15x15 neighbourhood. + %for each such points + [r, c] = find(Classfd); + for i = 1:length(r) + if (segmentedImg(r(i), c(i)) == 0) + candLbls = segmentedImg(r(i)-ofstSeg:r(i)+ofstSeg, c(i)-ofstSeg:c(i)+ofstSeg); + candLblDst = distM .* (candLbls > 0); + candLblDst(candLblDst == 0) = Inf; + candLblDst = (candLblDst == min(candLblDst(:))); + closestLbl = max(candLbls(candLblDst)); %believing that there is a cell in 15x15 loc + if (closestLbl > 0) + ClassfdClLbl = [ClassfdClLbl; closestLbl]; + else + %No cell in the visinity + excludeClsfdCnt = excludeClsfdCnt + 1; + end + end + end + end + ClassfdClLbl = unique(ClassfdClLbl); + if (~isempty(ClassfdClLbl) && ClassfdClLbl(1) == 0); ClassfdClLbl(1) = []; end; + + + %count TP and FN + TP = 0; FN = 0; + for i = 1:length(GndTrthClLbl) + curCell = GndTrthClLbl(i); + if (sum(curCell == ClassfdClLbl) > 0) + TP = TP + 1; + else + FN = FN + 1; + end + end + + FP = length(ClassfdClLbl) - TP; + TN = double(max(segmentedImg(:)) - (TP + FN + FP)); +% imshow(segWithBndry); +end diff --git a/src/gopa kumar code/generate32X32FcsStck_Fcsd_Feat_RGB.m b/src/gopa kumar code/generate32X32FcsStck_Fcsd_Feat_RGB.m new file mode 100644 index 0000000..5fecb04 --- /dev/null +++ b/src/gopa kumar code/generate32X32FcsStck_Fcsd_Feat_RGB.m @@ -0,0 +1,289 @@ +function generate32X32FcsStck_Fcsd_Feat_RGB() +alreadySegmented = true; + getTrainingSet32X32FcsStk_BFcsd_Feat_RGB('Type1', false, alreadySegmented); %This function get the fcsRGB patches, BFczdPathes and Features on RGB and stores it + getTrainingSet32X32FcsStk_BFcsd_Feat_RGB('Type2', false, alreadySegmented); + getTrainingSet32X32FcsStk_BFcsd_Feat_RGB('Type1', true, alreadySegmented); %This function get the fcsRGB patches, BFczdPathes and Features on RGB and stores it + getTrainingSet32X32FcsStk_BFcsd_Feat_RGB('Type2', true, alreadySegmented); +end + +function [FnlFeatSet, CellId] = getMyFeatures(cellsImg, imUint8RGB, slideName) + %Lets do the processing on Green channel for the time being + avMsk = fspecial('average', 3); + FnlFeatSet = single(zeros(1, 14*3)); %3 Channels + for ch = 1:3 + imGreen = imUint8RGB(:, :, ch); + imGreenDble = im2double(imGreen); + [grdMg, ~] = imgradient(imGreenDble); + meanLclMsk = imfilter(imGreenDble, avMsk); + varnLclMsk = imfilter(imGreenDble.^2, avMsk) - meanLclMsk.^2; + [SzM, SzN] = size(imGreen); + props = regionprops(cellsImg, {'PixelIdxList', 'Centroid'}); + [numObjs, ~] = size(props); + myBgndImgI = uint8(200* ones(SzM, SzN)); + FeatSet = single(zeros(numObjs, 14)); + CellId = []; + for i = 1:numObjs + pxlIdxLst = props(i).PixelIdxList; + [SR, SC] = ind2sub(size(imGreen), pxlIdxLst); + tmpI = myBgndImgI; + tmpI(pxlIdxLst) = imGreen(pxlIdxLst); + %Get the patch + minSR = min(SR); maxSR = max(SR); + minSC = min(SC); maxSC = max(SC); + imgPatch = tmpI(minSR:maxSR, minSC:maxSC); + mskPatch = cellsImg(minSR:maxSR, minSC:maxSC); + FeatSet(i, 1:4) = getGLCMFeatPrPatch(imgPatch, mskPatch); + + imGrnDblPchPxls = imGreenDble(pxlIdxLst); + meanGlblPxls = mean(imGrnDblPchPxls); + varGlblPxls = var(imGrnDblPchPxls); + minPxls = min(imGrnDblPchPxls); + maxPxls = max(imGrnDblPchPxls); + minGrdMag = min(grdMg(pxlIdxLst)); + maxGrdMag = max(grdMg(pxlIdxLst)); + minLclMean = min(meanLclMsk(pxlIdxLst)); + maxLclMean = max(meanLclMsk(pxlIdxLst)); + minLclVarn = min(varnLclMsk(pxlIdxLst)); + maxLclVarn = max(varnLclMsk(pxlIdxLst)); + FeatSet(i,5:14) = [meanGlblPxls varGlblPxls minPxls maxPxls ... + minGrdMag maxGrdMag minLclMean maxLclMean ... + minLclVarn maxLclVarn]; + cellCentr = round(props(i).Centroid); + CellId(i).name = [slideName '_' num2str(cellCentr(2)) '_' num2str(cellCentr(1))]; + end + FnlFeatSet(1, (ch-1)*14+1:ch*14) = FeatSet; + end +end + +function GLCMFeat = getGLCMFeatPrPatch(imUint8Gray, Msk) + im = imUint8Gray; + %im has to be Single band [0 255] + [SzR, SzC, SzO] = size(im); + MskX = [Msk(:, 2:end) Msk(:, end)]; + MskY = [Msk(2:end, :); Msk(end, :)]; + MskD = Msk & MskX & MskY; %MskD = bwmorph(MskD, 'erode', 3); + %Compute the GLCM for each band for the reg and return the props. + Lvl = 32; LvlDiv = 256/Lvl; + regIntst = MskD; + GLCMMat = zeros(Lvl, Lvl); + GLCMFeat = zeros(SzO, 4); + for band = 1:SzO + imBand = double(im); + imgL = ceil((imBand+1)/LvlDiv); + imgL(~regIntst) = -1; + for i = 1:SzR + for j = 1:SzC-1 + valLeft = imgL(i, j); + valRight = imgL(i, j+1); + if (valLeft ~= -1 && valRight ~= -1) + GLCMMat(valLeft, valRight) = GLCMMat(valLeft, valRight)+1; + end + end + end + GFeat = graycoprops(GLCMMat); + GLCMFeat(band, :) = [GFeat.Contrast GFeat.Correlation GFeat.Energy GFeat.Homogeneity]; + end +end + +function getTrainingSet32X32FcsStk_BFcsd_Feat_RGB(Folder, dustExclFrTrain, alreadySegmented) +%This function saves the 32x32 focus stack patches of 5th, 21st, and 37th +%images of the stck FczdRGB, stote the 21st alone as BFczdRGB, and the +%features of the BFcsdRGB as BFcsdFeatRGB. +% + + if (dustExclFrTrain) + svPath = 'CNNPatches\DntUseDst4NgPatc\'; + else + svPath = 'CNNPatches\DoUseDst4NgPatc\'; + end + + if (alreadySegmented) + load ('CNNPatches\Segment\alreadySegImg'); + load ('CNNPatches\Segment\alreadyWBCImg'); + load ('CNNPatches\Segment\segIndxNames'); + end + + btchSz = 12000; + BasePath = ['E:\Gopakumar\GopakumarIISTDrive\Dataset\' ... + '4th IIST Visit\Mal01\Images\GndTruth\Separate\NewEasy\']; + fnames = dir([BasePath '\ImgsMinAcrsStckTwoCls\' Folder '\*B.jpg']); + FczStkPath = 'E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\Patch32By32\DataSet32By32Color\Manual_But_Sel_Training\AllFcsStack\'; + + numfids = length(fnames); M = 32; N = 32; R = 1; Rad = 0.75*M; + dskMskPstv = strel('disk', Rad); dskMskDust = strel('disk', Rad); + posCnt = 0; negCnt = 0; negChnk = 1; prevVidIndx = ''; + for K = 1:numfids + K + AbsFNme = [BasePath '\ImgsMinAcrsStckTwoCls\' Folder '\' fnames(K).name]; + %Get details needed to acces the Ground Truth file + [~, FileName, ~] = fileparts(AbsFNme); + for i = 1:length(FileName); if (FileName(i) == '_'); vid = i-1; break; end; end; + for j = vid+2:length(FileName); if (FileName(j) == '_'); stck = j-1; break; end; end; + vidIndx = FileName(1:vid); stckIndx = FileName(vid+2:stck); + slideName = [vidIndx '_' stckIndx]; + %Read the min RGB image + img = imread(AbsFNme); [SizeR, SizeC, ~] = size(img); + if(isTherWBC([vidIndx '_' stckIndx])) + if (~alreadySegmented) + [segmentedImg, segWithBndry, WBCs] = getSegmentation(img); + else + IndxPosSlctd = getSegmentedIndx([vidIndx '_' stckIndx], segIndxNames); + segmentedImg = alreadySegImg(:, :, IndxPosSlctd); + WBCs = alreadyWBCImg(:, :, IndxPosSlctd); + end + else + WBCs = false(SizeR, SizeC); + end + imF1 = im2single(imread([FczStkPath vidIndx '_' stckIndx '_5.jpg'])); + imF = imread([FczStkPath vidIndx '_' stckIndx '_21.jpg']); imF2 = im2single(imF); + imF3 = im2single(imread([FczStkPath vidIndx '_' stckIndx '_37.jpg'])); + imG(:, :, 1:3) = imF1; imG(:, :, 4:6) = imF2; imG(:, :, 7:9) = imF3; +% imG = 0; +% imG = img; %(:, :, 2); +% imG = im2single(imG); + [SzM, SzN, SzO] = size(imG); + %Get the positive patches + load ([BasePath 'GndTrth\' vidIndx '_' stckIndx]); + if (~strcmp(prevVidIndx, vidIndx)); + prevVidIndx = vidIndx; + DstMsk = getDstLocnsFor(vidIndx); + end + + load ([BasePath 'GndTrth\Dst_' vidIndx '_' stckIndx]); + [PlocR, PlocC] = find (Msk); + falseMsk = false(size(Msk)); + for i = 1:length(PlocR) + + %For the time being, if there is a patch of the required size + %around the point, then only we are considering it. + minRw = PlocR(i) - M/2; minCl = PlocC(i) - N/2; + maxRw = PlocR(i) + M/2-1; maxCl = PlocC(i) + N/2-1; + if (minRw > 0 && minCl > 0 && maxRw <= SzM && maxCl <= SzN) + currPatch = imG(minRw:maxRw, minCl:maxCl, :); + posCnt = posCnt + 1; + postvStackFcsRGB(1:M, 1:N, 1:SzO, posCnt) = currPatch; + postvBestFczdRGB(1:M, 1:N, 1:3, posCnt) = imG(minRw:maxRw, minCl:maxCl, 4:6); + falseMsk (minRw:maxRw, minCl:maxCl) = true; + [FeatSet, ~] = getMyFeatures(falseMsk, imF, slideName); + falseMsk (minRw:maxRw, minCl:maxCl) = false; + postvBstFczdFeat(posCnt, :) = FeatSet; + end + end + + %Get the negative patches by Masking out positive regions + %and Masking out Dust Regions + mskdPstv = imdilate(Msk, dskMskPstv); + mskdDust = imdilate(DstMsk, dskMskDust); + %Find regional minima + regMin = getMySuspectedRegion(img); + %Get negative Positions + if (dustExclFrTrain) + mskNgtv = regMin & ~mskdPstv & ~mskdDust & ~WBCs; + else + mskNgtv = regMin & ~mskdPstv & ~WBCs; + end + %the regmin can be a group of pixels. So select pixels with respect + %to the centroids +% mySusp = false(size(mskNgtv)); + Cntrids = regionprops(mskNgtv, 'centroid'); + [numNgtvs, ~] = size(Cntrids); + for i = 1:numNgtvs + + %For the time being, if there is a patch of the required size + %around the point, then only we are considering it. + CCntrids = round(Cntrids(i).Centroid); + minRw = CCntrids(2) - M/2; minCl = CCntrids(1) - N/2; + maxRw = CCntrids(2) + M/2-1; maxCl = CCntrids(1) + N/2-1; + if (minRw > 0 && minCl > 0 && maxRw <= SzM && maxCl <= SzN) + %Because of the size constraint lets sample one out of five + decideTTake = rand(); outOf = 1; + if (decideTTake >= (1 - 1/outOf)) + currPatch = imG(minRw:maxRw, minCl:maxCl, :); + negCnt = negCnt + 1; + negtvStackFcsRGB(1:M, 1:N, 1:SzO, negCnt) = currPatch; + negtvBestFczdRGB(1:M, 1:N, 1:3, negCnt) = imG(minRw:maxRw, minCl:maxCl, 4:6); + falseMsk (minRw:maxRw, minCl:maxCl) = true; + [FeatSet, ~] = getMyFeatures(falseMsk, imF, slideName); + falseMsk (minRw:maxRw, minCl:maxCl) = false; + negtvBstFczdFeat(negCnt, :) = FeatSet; + + + if (negCnt == btchSz) + save ([svPath 'FczRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvStackFcsRGB'); + save ([svPath 'BFczdRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvBestFczdRGB'); + save ([svPath 'BFczdFeatRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvBstFczdFeat'); + negChnk = negChnk+1; + negCnt = 0; + end + end +% mySusp(CCntrids(2), CCntrids(1)) = true; + end + end +% dilMySusp = bwmorph(mySusp, 'dilate', 5); +% r = img(:, :, 1); g = img(:, :, 2); b = img(:, :, 3); +% r(dilMySusp) = 0; g(dilMySusp) = 0; b(dilMySusp) = 255; +% mySuspIm(:, :, 1) = r; mySuspIm(:, :, 2) = g; mySuspIm(:, :, 3) = b; +% imshow(img); title('original'); figure; +% imshow(mySuspIm); title('original'); +% close all; + end + save ([svPath 'FczRGBPsStk_' Folder], 'postvStackFcsRGB'); + save ([svPath 'BFczdRGBPsStk_' Folder], 'postvBestFczdRGB'); + save ([svPath 'BFczdFeatRGBPsStk_' Folder], 'postvBstFczdFeat'); + + negtvStackFcsRGB = negtvStackFcsRGB(1:M, 1:N, 1:SzO, 1:negCnt); %+1:end) = []; + negtvBestFczdRGB = negtvBestFczdRGB(1:M, 1:N, :, 1:negCnt); + negtvBstFczdFeat = negtvBstFczdFeat(1:negCnt, :); + + save ([svPath 'FczRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvStackFcsRGB'); + save ([svPath 'BFczdRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvBestFczdRGB'); + save ([svPath 'BFczdFeatRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvBstFczdFeat'); +end + +function IndxPosSlctd = getSegmentedIndx(indxName, segIndxNames) + [~, numNames] = size(segIndxNames); + for i = 1:numNames + if (strcmp(indxName, segIndxNames(i).name)) + IndxPosSlctd = i; + break; + end + end +end + +function stat = isTherWBC(fileName) + stat = false; + indx = {'1_9', '2_19', '3_3', '3_4', '4_8', '5_1', 'A_6', 'A_14', 'BB_3', 'BB_13', ... + 'BB_14', 'C_15', 'C_16', 'D_13', 'E_2', 'EE_4', 'F_3', 'F_5', 'F_9', 'FF_6', ... + 'GG_4', 'G_16', 'G_17', 'H_5', 'H_16', 'I_1', 'J_18', 'KK_9', 'L_6', ... + 'O_8', 'P_12', 'P_13', 'S_16', 'T_12', 'V_8', 'V_9', 'W_8', 'W_18', 'X_11', ... + 'X_13', 'Y_2', 'Z_11'}; + [~, num] = size(indx); + for i = 1:num + if (strcmp(fileName, indx(i))) + stat = true; + break; + end + end +end +function cDstLoc = getDstLocnsFor(vidIndx) + Path = 'E:\Gopakumar\GopakumarIISTDrive\Dataset\4th IIST Visit\Mal01\Images\GndTruth\Separate\NewEasy\DstLocnsPgm\'; + load ([Path 'DstLocByPGM']); cnt = 0; + vidIndxs = DstLocByPGM.vidIndx; + [~, numVidIndxs] = size(vidIndxs); + for i = 1:numVidIndxs + if (strcmp(vidIndx, vidIndxs(i).name)) + cDstLoc = DstLocByPGM.DstLocByPgm(:, :, i); + break; + end + end +end +function susReg = getMySuspectedRegion(img) + MnSz = 900; + hsv = rgb2hsv(img); + vlue = hsv(:, :, 3); + stDsk = strel('disk', 7); + openimg = imopen(vlue, stDsk); + mask = imregionalmin(openimg); + %First Filtering Exclude all the Bgnd + susReg = (bwareaopen((vlue < graythresh(vlue)), MnSz)) & mask; +end \ No newline at end of file diff --git a/src/gopa kumar code/generateDstAndWBCPatches.m b/src/gopa kumar code/generateDstAndWBCPatches.m new file mode 100644 index 0000000..e36ae9f --- /dev/null +++ b/src/gopa kumar code/generateDstAndWBCPatches.m @@ -0,0 +1,282 @@ +function generateDstAndWBCPatches(alreadySegmented) + getTrainingSet32X32FcsStk_BFcsd_Feat_RGB('Type1', true, alreadySegmented); %This function get the fcsRGB patches, BFczdPathes and Features on RGB and stores it + getTrainingSet32X32FcsStk_BFcsd_Feat_RGB('Type2', true, alreadySegmented); +end + +function [FnlFeatSet, CellId] = getMyFeatures(cellsImg, imUint8RGB, slideName) + %Lets do the processing on Green channel for the time being + avMsk = fspecial('average', 3); + FnlFeatSet = single(zeros(1, 14*3)); %3 Channels + for ch = 1:3 + imGreen = imUint8RGB(:, :, ch); + imGreenDble = im2double(imGreen); + [grdMg, ~] = imgradient(imGreenDble); + meanLclMsk = imfilter(imGreenDble, avMsk); + varnLclMsk = imfilter(imGreenDble.^2, avMsk) - meanLclMsk.^2; + [SzM, SzN] = size(imGreen); + props = regionprops(cellsImg, {'PixelIdxList', 'Centroid'}); + [numObjs, ~] = size(props); + myBgndImgI = uint8(200* ones(SzM, SzN)); + FeatSet = single(zeros(numObjs, 14)); + CellId = []; + for i = 1:numObjs + pxlIdxLst = props(i).PixelIdxList; + [SR, SC] = ind2sub(size(imGreen), pxlIdxLst); + tmpI = myBgndImgI; + tmpI(pxlIdxLst) = imGreen(pxlIdxLst); + %Get the patch + minSR = min(SR); maxSR = max(SR); + minSC = min(SC); maxSC = max(SC); + imgPatch = tmpI(minSR:maxSR, minSC:maxSC); + mskPatch = cellsImg(minSR:maxSR, minSC:maxSC); + FeatSet(i, 1:4) = getGLCMFeatPrPatch(imgPatch, mskPatch); + + imGrnDblPchPxls = imGreenDble(pxlIdxLst); + meanGlblPxls = mean(imGrnDblPchPxls); + varGlblPxls = var(imGrnDblPchPxls); + minPxls = min(imGrnDblPchPxls); + maxPxls = max(imGrnDblPchPxls); + minGrdMag = min(grdMg(pxlIdxLst)); + maxGrdMag = max(grdMg(pxlIdxLst)); + minLclMean = min(meanLclMsk(pxlIdxLst)); + maxLclMean = max(meanLclMsk(pxlIdxLst)); + minLclVarn = min(varnLclMsk(pxlIdxLst)); + maxLclVarn = max(varnLclMsk(pxlIdxLst)); + FeatSet(i,5:14) = [meanGlblPxls varGlblPxls minPxls maxPxls ... + minGrdMag maxGrdMag minLclMean maxLclMean ... + minLclVarn maxLclVarn]; + cellCentr = round(props(i).Centroid); + CellId(i).name = [slideName '_' num2str(cellCentr(2)) '_' num2str(cellCentr(1))]; + end + FnlFeatSet(1, (ch-1)*14+1:ch*14) = FeatSet; + end +end + +function GLCMFeat = getGLCMFeatPrPatch(imUint8Gray, Msk) + im = imUint8Gray; + %im has to be Single band [0 255] + [SzR, SzC, SzO] = size(im); + MskX = [Msk(:, 2:end) Msk(:, end)]; + MskY = [Msk(2:end, :); Msk(end, :)]; + MskD = Msk & MskX & MskY; %MskD = bwmorph(MskD, 'erode', 3); + %Compute the GLCM for each band for the reg and return the props. + Lvl = 32; LvlDiv = 256/Lvl; + regIntst = MskD; + GLCMMat = zeros(Lvl, Lvl); + GLCMFeat = zeros(SzO, 4); + for band = 1:SzO + imBand = double(im); + imgL = ceil((imBand+1)/LvlDiv); + imgL(~regIntst) = -1; + for i = 1:SzR + for j = 1:SzC-1 + valLeft = imgL(i, j); + valRight = imgL(i, j+1); + if (valLeft ~= -1 && valRight ~= -1) + GLCMMat(valLeft, valRight) = GLCMMat(valLeft, valRight)+1; + end + end + end + GFeat = graycoprops(GLCMMat); + GLCMFeat(band, :) = [GFeat.Contrast GFeat.Correlation GFeat.Energy GFeat.Homogeneity]; + end +end + +function getTrainingSet32X32FcsStk_BFcsd_Feat_RGB(Folder, dustExclFrTrain, alreadySegmented) +%This function saves the 32x32 focus stack patches of 5th, 21st, and 37th +%images of the stck FczdRGB, stote the 21st alone as BFczdRGB, and the +%features of the BFcsdRGB as BFcsdFeatRGB. +% + getWBCDataSet = true; getDstDataSet = true; + if (dustExclFrTrain) + svPath = 'CNNPatches\DntUseDst4NgPatc\'; + else + svPath = 'CNNPatches\DoUseDst4NgPatc\'; + end + if (alreadySegmented) + load ('CNNPatches\Segment\alreadySegImg'); + load ('CNNPatches\Segment\alreadyWBCImg'); + load ('CNNPatches\Segment\segIndxNames'); + end + btchSz = 12000; + BasePath = ['E:\Gopakumar\GopakumarIISTDrive\Dataset\' ... + '4th IIST Visit\Mal01\Images\GndTruth\Separate\NewEasy\']; + fnames = dir([BasePath '\ImgsMinAcrsStckTwoCls\' Folder '\*B.jpg']); + FczStkPath = 'E:\Gopakumar\GopakumarIISTDrive\Dataset\Malaria\Patch32By32\DataSet32By32Color\Manual_But_Sel_Training\AllFcsStack\'; + + numfids = length(fnames); M = 32; N = 32; R = 1; Rad = 0.75*M; + dskMskPstv = strel('disk', Rad); dskMskDust = strel('disk', Rad); + posCnt = 0; negCnt = 0; negChnk = 1; prevVidIndx = ''; + for K = 1:numfids + K + AbsFNme = [BasePath '\ImgsMinAcrsStckTwoCls\' Folder '\' fnames(K).name]; + %Get details needed to acces the Ground Truth file + [~, FileName, ~] = fileparts(AbsFNme); + for i = 1:length(FileName); if (FileName(i) == '_'); vid = i-1; break; end; end; + for j = vid+2:length(FileName); if (FileName(j) == '_'); stck = j-1; break; end; end; + vidIndx = FileName(1:vid); stckIndx = FileName(vid+2:stck); + slideName = [vidIndx '_' stckIndx]; + %Read the min RGB image + img = imread(AbsFNme); [SizeR, SizeC, ~] = size(img); + if(isTherWBC([vidIndx '_' stckIndx])) + if (~alreadySegmented) + [segmentedImg, segWithBndry, WBCs] = getSegmentation(img); + else + IndxPosSlctd = getSegmentedIndx([vidIndx '_' stckIndx], segIndxNames); + segmentedImg = alreadySegImg(:, :, IndxPosSlctd); + WBCs = alreadyWBCImg(:, :, IndxPosSlctd); + if (sum(WBCs(:) > 0)) + thereIsWBC = true; + end + end + else + WBCs = false(SizeR, SizeC); + thereIsWBC = false; + end + imF1 = im2single(imread([FczStkPath vidIndx '_' stckIndx '_5.jpg'])); + imF = imread([FczStkPath vidIndx '_' stckIndx '_21.jpg']); imF2 = im2single(imF); + imF3 = im2single(imread([FczStkPath vidIndx '_' stckIndx '_37.jpg'])); + imG(:, :, 1:3) = imF1; imG(:, :, 4:6) = imF2; imG(:, :, 7:9) = imF3; +% imG = 0; +% imG = img; %(:, :, 2); +% imG = im2single(imG); + [SzM, SzN, SzO] = size(imG); + %Get the positive patches + load ([BasePath 'GndTrth\' vidIndx '_' stckIndx]); + if (~strcmp(prevVidIndx, vidIndx)); + prevVidIndx = vidIndx; + DstMsk = getDstLocnsFor(vidIndx); + end + + load ([BasePath 'GndTrth\Dst_' vidIndx '_' stckIndx]); + [PlocR, PlocC] = find (Msk); + falseMsk = false(size(Msk)); + mskdPstv = imdilate(Msk, dskMskPstv); + mskdDust = imdilate(DstMsk, dskMskDust); + + regmin = getMySuspectedRegion(img); + if (thereIsWBC && getWBCDataSet) + WBCMsks = WBCs & regmin; + Cntrids = regionprops(WBCMsks, 'centroid'); + [numWBCPts, ~] = size(Cntrids); + %Add any missed location + newIntstRegMn = ~WBCs & regmin; + mYR = imF(:, :, 1); mYG = imF(:, :, 2); mYB = imF(:, :, 3); + mYR(newIntstRegMn) = 255; mYG(newIntstRegMn) = 0; mYB(newIntstRegMn) = 0; + showIm(:, :, 1) = mYR; showIm(:, :, 2) = mYG; showIm(:, :, 3) = mYB; + imshow(showIm); [NC, NR, Btn] = ginput(1); + while(Btn ~= 3) %until right button is pressed + numWBCPts = numWBCPts + 1; + Cntrids(numWBCPts).Centroid = [NC NR]; + [NC, NR, Btn] = ginput(1); + end + + for i = 1:numWBCPts + CCntrids = round(Cntrids(i).Centroid); + %For the time being, if there is a patch of the required size + %around the point, then only we are considering it. + minRw = CCntrids(2) - M/2; minCl = CCntrids(1) - N/2; + maxRw = CCntrids(2) + M/2-1; maxCl = CCntrids(1) + N/2-1; + if (minRw > 0 && minCl > 0 && maxRw <= SzM && maxCl <= SzN) + currPatch = imG(minRw:maxRw, minCl:maxCl, :); + posCnt = posCnt + 1; + WBCStackFcsRGB(1:M, 1:N, 1:SzO, posCnt) = currPatch; + WBCBestFczdRGB(1:M, 1:N, 1:3, posCnt) = imG(minRw:maxRw, minCl:maxCl, 4:6); + falseMsk (minRw:maxRw, minCl:maxCl) = true; + [FeatSet, ~] = getMyFeatures(falseMsk, imF, slideName); + falseMsk (minRw:maxRw, minCl:maxCl) = false; + WBCBstFczdFeat(posCnt, :) = FeatSet; + end + end + thereIsWBC = false; + end + + + if (getDstDataSet) + Cntrids = regionprops(DstMsk, 'centroid'); + [numDstPts, ~] = size(Cntrids); + + for i = 1:numDstPts + CCntrids = round(Cntrids(i).Centroid); + %For the time being, if there is a patch of the required size + %around the point, then only we are considering it. + minRw = CCntrids(2) - M/2; minCl = CCntrids(1) - N/2; + maxRw = CCntrids(2) + M/2-1; maxCl = CCntrids(1) + N/2-1; + if (minRw > 0 && minCl > 0 && maxRw <= SzM && maxCl <= SzN) + currPatch = imG(minRw:maxRw, minCl:maxCl, :); + posCnt = posCnt + 1; + DustStackFcsRGB(1:M, 1:N, 1:SzO, posCnt) = currPatch; + DustBestFczdRGB(1:M, 1:N, 1:3, posCnt) = imG(minRw:maxRw, minCl:maxCl, 4:6); + falseMsk (minRw:maxRw, minCl:maxCl) = true; + [FeatSet, ~] = getMyFeatures(falseMsk, imF, slideName); + falseMsk (minRw:maxRw, minCl:maxCl) = false; + DustBstFczdFeat(posCnt, :) = FeatSet; + end + end + end + end + if (getDstDataSet) + save ([svPath 'FczRGBDsStk_' Folder], 'DustStackFcsRGB'); + save ([svPath 'BFczdRGBDsStk_' Folder], 'DustBestFczdRGB'); + save ([svPath 'BFczdFeatRGBDsStk_' Folder], 'DustBstFczdFeat'); + end + if (getWBCDataSet) + save ([svPath 'FczRGBWcStk_' Folder], 'WBCStackFcsRGB'); + save ([svPath 'BFczdRGBWcStk_' Folder], 'WBCBestFczdRGB'); + save ([svPath 'BFczdFeatRGBWcStk_' Folder], 'WBCBstFczdFeat'); + end + +% negtvStackFcsRGB = negtvStackFcsRGB(1:M, 1:N, 1:SzO, 1:negCnt); %+1:end) = []; +% negtvBestFczdRGB = negtvBestFczdRGB(1:M, 1:N, :, 1:negCnt); +% negtvBstFczdFeat = negtvBstFczdFeat(1:negCnt, :); +% +% save ([svPath 'FczRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvStackFcsRGB'); +% save ([svPath 'BFczdRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvBestFczdRGB'); +% save ([svPath 'BFczdFeatRGBNgStk_' Folder '_' num2str(negChnk)], 'negtvBstFczdFeat'); +end +function IndxPosSlctd = getSegmentedIndx(indxName, segIndxNames) + [~, numNames] = size(segIndxNames); + for i = 1:numNames + if (strcmp(indxName, segIndxNames(i).name)) + IndxPosSlctd = i; + break; + end + end +end +function stat = isTherWBC(fileName) + stat = false; + indx = {'1_9', '2_19', '3_3', '3_4', '4_8', '5_1', 'A_6', 'A_14', 'BB_3', 'BB_13', ... + 'BB_14', 'C_15', 'C_16', 'D_13', 'E_2', 'EE_4', 'F_3', 'F_5', 'F_9', 'FF_6', ... + 'GG_4', 'G_16', 'G_17', 'H_5', 'H_16', 'I_1', 'J_18', 'KK_9', 'L_6', ... + 'O_8', 'P_12', 'P_13', 'S_16', 'T_12', 'V_8', 'V_9', 'W_8', 'W_18', 'X_11', ... + 'X_13', 'Y_2', 'Z_11'}; + [~, num] = size(indx); + for i = 1:num + if (strcmp(fileName, indx(i))) + stat = true; + break; + end + end +end +function cDstLoc = getDstLocnsFor(vidIndx) + Path = 'E:\Gopakumar\GopakumarIISTDrive\Dataset\4th IIST Visit\Mal01\Images\GndTruth\Separate\NewEasy\DstLocnsPgm\'; + load ([Path 'DstLocByPGM']); cnt = 0; + vidIndxs = DstLocByPGM.vidIndx; + [~, numVidIndxs] = size(vidIndxs); + for i = 1:numVidIndxs + if (strcmp(vidIndx, vidIndxs(i).name)) + cDstLoc = DstLocByPGM.DstLocByPgm(:, :, i); + break; + end + end +end +function susReg = getMySuspectedRegion(img) + MnSz = 900; + hsv = rgb2hsv(img); + vlue = hsv(:, :, 3); + stDsk = strel('disk', 7); + openimg = imopen(vlue, stDsk); + mask = imregionalmin(openimg); + %First Filtering Exclude all the Bgnd + susReg = (bwareaopen((vlue < graythresh(vlue)), MnSz)) & mask; +end \ No newline at end of file diff --git a/src/gopa kumar code/getSegmentation.m b/src/gopa kumar code/getSegmentation.m new file mode 100644 index 0000000..3a8e418 --- /dev/null +++ b/src/gopa kumar code/getSegmentation.m @@ -0,0 +1,649 @@ +function [segmentedImg, segWithBndry, WBCs, remSetGTh] = getSegmentation(currImgUint8) +% load 'imgUint8'; currImgUint8 = imgUint8; +%2.19 micrometer 15 pixels + %note that the overlapping ellipses overwrite the labels. Fr the + %time being OK. need to improve later + currImg = currImgUint8; + [clearSetLTh, clearSetGTh, remSetGTh] = getBinaryImage(im2double(currImg)); + clearSetLTh = bwfill(clearSetLTh, 'holes'); clearSetGTh = bwfill(clearSetGTh, 'holes'); + + clsLThP = bwperim(clearSetLTh); clsGThP = bwperim(clearSetGTh); remGThP = bwperim(remSetGTh); + + fnalSeg = clearSetLTh | clearSetGTh| remSetGTh; + fnalSeg(clsLThP | clsGThP | remGThP) = false; + %Move WBCs to clear set. + [WBCs, pbRBC, pbWBC, pbBGD] = getWBCs(im2double(currImg), fnalSeg); + WBCs = bwfill(WBCs, 'holes'); + wbcP = bwperim(WBCs); + clearSet = clearSetLTh | clearSetGTh | WBCs; + clearSet(clsLThP | clsGThP | wbcP) = false; + remSetGTh(WBCs) = false; + remSetGTh = bwareaopen(remSetGTh, 500); +% [clustImg, elpseLbld, fitImg] = kMeansByGeoDesicAutomatic(remSetGTh); + + [segLbld, elpseLbld, fitImg] = kMeansByGeoDesicCentrWatershed(currImgUint8, clearSet, remSetGTh); + %Now assign a lable starting next to maxLbl assigned to ellipses + maxLb = max(segLbld(:)); + cmp = bwconncomp(clearSet); + currLb = maxLb; segmentedImg = segLbld; + for lblOfst = 1:cmp.NumObjects + currLb = currLb+1; + segmentedImg(cmp.PixelIdxList{lblOfst}) = currLb; + end + segmentedImg = uint8(segmentedImg); + remP = bwperim(segLbld); + %For Display Purpose + fnalPerm = wbcP | remP | clsLThP | clsGThP; + R = currImgUint8(:, :, 1); G = currImgUint8(:, :, 2); B = currImgUint8(:, :, 3); + R(fnalPerm) = 255; G(fnalPerm) = 0; B(fnalPerm) = 0; + segWithBndry (:, :, 1) = R; segWithBndry (:, :, 2) = G; segWithBndry (:, :, 3) = B; + %Believing that there is max 255 cells +end + +function remSetGTh = closelyFilter(remSetGTh, im) +close all; + intstLoc = im .* remSetGTh; + im = im2double(im); + figure; imshow(remSetGTh); + loG = fspecial('log'); + GF = imfilter(im, loG); + avgK = fspecial('average', 3); + meanIm = imfilter(GF, avgK); + figure; imshow(GF.*remSetGTh); +end + +function [clearSetLTh, clearSetGTh, remSetGTh] = getBinaryImage(im) + imG = im2double(im(:, :, 2)); Sz = 21; + msk = fspecial('gaussian', Sz, Sz/6); hlfSz = round(Sz/2); + avG1 = imfilter(imG, msk); + msk(1:hlfSz, 1:hlfSz) = rot90(msk(1:hlfSz, 1:hlfSz), 2); + msk(1:hlfSz, hlfSz+1:end) = rot90(msk(1:hlfSz, hlfSz+1:end), 2); + msk(hlfSz+1:end, 1:hlfSz) = rot90(msk(hlfSz+1:end, 1:hlfSz), 2); + msk(hlfSz+1:end, hlfSz+1:end) = rot90(msk(hlfSz+1:end, hlfSz+1:end), 2); + +% msk(:, 1:8) = fliplr(msk(:, 1:8)); msk(:, 9:15) = fliplr(msk(:, 9:15)); + %Lets find weighted Mean + avG2 = imfilter(imG, msk); + avGB(:, :, 1) = avG1; avGB(:, :, 2) = avG2; + avG = min(avGB, [], 3); + %if the avg Belongs to White region? Exclude : In effect (set at least 0.1 + + %0.01) thresh + avG(avG > 0.7) = 0.6; + fildLThImg = bwfill(bwareaopen(imG < (avG - 0.005), 200), 'holes'); %note that simple avg with 0.015 thresh gave good results 0.01 + fildGThImg = lOtThresh(im); + [clearSetLTh, ~, remClsSet] = getClearSet(bwfill(fildLThImg, 'holes')); + remSetGTh = bwareaopen(bwmorph(~(clearSetLTh | remClsSet) & fildGThImg, 'open', 3), 500); + [clearSetGTh, remSetGTh, clsGth] = getClearSet(bwfill(remSetGTh, 'holes')); + remSetGTh = remSetGTh | bwmorph(clsGth, 'erode', 1) | bwmorph(remClsSet, 'erode', 1); + clearSet = clearSetLTh | clearSetGTh; clearSet(bwperim(bwmorph(clearSetGTh, 'dilate', 1))) = false; +end + +function [WBCsActual, pbRBC, pbWBC, pbBGD] = getWBCs(imRGBDble, segImg) + form = makecform ('srgb2lab'); + LAB = applycform(imRGBDble, form); + A = LAB(:, :, 2); + B = LAB(:, :, 3); + pbRBC = (exp(-0.5*((A - 14.6585).^2/10.1368))) .* (exp(-0.5*((B - -2.8583).^2/8.6485))); %((1/sqrt(2*pi*8.6485)) * + pbWBC = (exp(-0.5*((A - 28.4699).^2/106.8322))) .* (exp(-0.5*((B - -31.1815).^2/207.2493))); %((1/sqrt(2*pi*97.4714)) * + pbBGD = (exp(-0.5*((A - -6.6678).^2/3.3932))) .* (exp(-0.5*((B - 1.5055).^2/3.6731))); + RBCs = (pbRBC > pbWBC) & (pbRBC > pbBGD) & (pbRBC > 0.1); + WBCs = (pbWBC > pbRBC) & (pbWBC > pbBGD) & (pbWBC > 0.01); + %Take better RBCs + WBCsErde = bwareaopen(bwmorph(WBCs, 'erode', 3), 500); + [~, WBCsActual] = maskAllSharingObjects(segImg, WBCsErde); WBCsActual = bwfill(bwmorph(WBCsActual, 'dilate', 1), 'holes'); + WBCsCand = bwfill(bwmorph(bwareaopen(bwmorph(WBCsActual & ~RBCs, 'erode', 3),200), 'dilate', 3), 'holes'); + if (sum(WBCsCand(:))) + WBCsActual = ensureWBCs(imRGBDble, WBCsCand); + end + %If the solididty of this one identified is too small, probably it is + %an RBC infected by parasite + + %Note that the WBC from pbWBC is much better than the WBC just found. + %except that it contain some RBC parts. This can be solved by taking out + %RBCs and masking out. Rather than masking in WBCs. +end + +function WBCsActual = ensureWBCs(imRGBDble, WBCsCand) + regProps = regionprops(WBCsCand, {'PixelIdxList', 'Solidity'}); + [numObj, ~] = size(regProps); + RC = imRGBDble(:, :, 1); GC = imRGBDble(:, :, 2); BC = imRGBDble(:, :, 3); + WBCsActual = WBCsCand; thresh = 100/255; + for i = 1:numObj + currObj = regProps(i).PixelIdxList; + %If heavily infected by parasite, there will be darker pixels + numPixBlack = sum((RC(currObj) < thresh) & (GC(currObj) < thresh) & (BC(currObj) < thresh)); + if(numPixBlack > 10) + %Infected RBC + WBCsActual(currObj) = false; + elseif (regProps(i).Solidity < 0.75) + WBCsActual(currObj) = false; + end + end +end +function [clearSet, remSet, remClseSet] = getClearSet(bwImg) + remSet = bwImg; clearSet = false(size(remSet)); remClseSet = clearSet; + cmps = regionprops(bwImg, {'PixelIdxList', 'Solidity', 'ConvexImage', 'BoundingBox', 'Eccentricity'}); + [numObj, ~] = size(cmps); + lowThresh = 750; highThresh = 1750; maxThresh = 2200; + for i = 1:numObj + currObj = cmps(i).PixelIdxList; + objArea = length(currObj); + if (objArea >= highThresh && objArea < maxThresh) + Tlr = 0.1; + else + Tlr = 0; + end + if ((objArea > lowThresh && objArea < maxThresh)) + if (cmps(i).Solidity > (0.85 + Tlr)) + clearSet(currObj) = true; + remSet(currObj) = false; + elseif (cmps(i).Solidity > (0.85 + Tlr)) + remClseSet (currObj) = true; + remSet(currObj) = false; + end + end + end +end + +function ThImg = lOtThresh(im) + globTh = graythresh(im); + div = 2; + [SzM, SzN, SzO] = size(im); + if (SzO == 3) + im = rgb2gray(im); %im(:, :, 2); + end + im = im2double(im); + rOfst = round(SzM/div); cOfst = round(SzN/div); + ThImg = false(SzM, SzN); + for i = 1:div + startR = (i-1)*rOfst+1; + if (i == div) + endR = SzM; + else + endR = i*rOfst; + end + for j = 1:div + startC = (j-1)*cOfst+1; + if (j == div) + endC = SzN; + else + endC = j*cOfst; + end + divImg = im(startR:endR, startC:endC); + locTh = graythresh(divImg); + if (locTh > 1.1*globTh) + ThImg(startR:endR, startC:endC) = false(size(divImg)); %probably there is no cell here + else + ThImg(startR:endR, startC:endC) = divImg < (locTh + 0.01); + end + end + end + ThImg = bwareaopen(ThImg, 500); %bwfill(, 'holes'); +% else +% +% end +end +function ClustImg = waterShedBasedOnGeoDesicDist(img) +ClustImg = []; +close all; + imb = img; %img(22:129, 393:525); %img; % +% imb = imb(20:35, 70:90); + imshow(imb); +% imb = true(7, 8); imb(1, 1:5) = false; imb(2, 3:4) = false; imb(3, 4) = false; + [SzR, SzC] = size(imb); + imagesc(imb); title('Click On rough Center Posns : To return, click right button'); C = []; R = []; + while (true) + [Cc, Rr, Btn] = ginput(1); + if (Btn ~= 1) + break; + end + C = [C Cc]; R = [R Rr]; + end +% C = [38 66 96]; R = [71 46 28]; + C = round(C); R = round(R); numCells = length(C); + imbc = imb; imbc(R(1), C(1)) = 0; imbc(R(2), C(2)) = 0; imbc(R(3), C(3)) = 0; + imagesc(imbc); + + [X, Y] = meshgrid(1:SzC, 1:SzR); + FlImg = true(size(imb)); + figure; + +% for i = 1:numCells +% CGD = bwdistgeodesic(imb, C(i), R(i), 'quasi-euclidean'); GD(:, :, i) = CGD; +% TED = bwdistgeodesic(FlImg, C(i), R(i), 'quasi-euclidean'); ED(:, :, i) = TED; +% CCD = CGD + abs(TED - CGD); CD(:, :, i) = CCD; %CGD+ +% imagesc(CGD > TED); figure; +% end +% +% Cm = min(CD, [], 3); + Cm = bwdistgeodesic(imb, C, R, 'quasi-euclidean'); + Cm(isnan(Cm)) = Inf; + lb = watershed(Cm); + figure; imshow(lb .* uint8(imb), []); + +% DistU = D; DistU(isnan(D)) = Inf; +% Lbl = watershed(DistU); +% ClustImg = uint8(imb) .* uint8(Lbl); +% figure; imshow(ClustImg, []); +end + +function [ClustImg, elpseLbld, fitImg] = kMeansByGeoDesic(img) +% close all; +% imb = img(22:129, 393:525); %img; %img(182:228, 127:198); % + imb = img; %im2bw(imread('C:\Users\GGG\Desktop\b.jpg')); % + [SzR, SzC] = size(imb); + figure; + imshow(imb); title('Click On rough Center Posns : To return, click right button'); C = []; R = []; + while (true) + [Cc, Rr, Btn] = ginput(1); + if (Btn ~= 1) + break; + end + C = [C Cc]; R = [R Rr]; + end + C = round(C); R = round(R); numCells = length(C); + [X, Y] = meshgrid(1:SzC, 1:SzR); + MaxItrn = 4; FlImg = true(size(imb)); + for itrn = 1:MaxItrn + for i = 1:numCells + CGD = bwdistgeodesic(imb, C(i), R(i), 'quasi-euclidean'); + FGD = bwdistgeodesic(FlImg, C(i), R(i), 'quasi-euclidean'); + GD(:, :, i) = CGD + abs(CGD - FGD).*(0.75*CGD); + end + %Assign to cluster + [~, lbl] = min(GD, [], 3); + Lbls = lbl.*imb; +% figure; imshow(Lbls, []) + %Update Centroid + for i = 1:numCells + cLbl = (Lbls == i); + C(i) = mean(X(cLbl)); R(i) = mean(Y(cLbl)); + end + C = round(C); R = round(R); + end + ClustImg = Lbls; + + [elpseLbld, fitImg] = improveSheds(Lbls, imb); + +end + +function [ClustImg, elpseLbld, fitImg] = kMeansByGeoDesicAutomatic(bwImg) + %Clean my binary image. I.e., remove spurious projections + dsk = strel('disk', 5); bwImg = imdilate(imerode(bwImg, dsk), dsk); + ClustImg = false(size(bwImg)); elpseLbld = ClustImg; fitImg = ClustImg; + imb = bwImg; %im2bw(imread('C:\Users\GGG\Desktop\b.jpg')); % + [SzR, SzC] = size(imb); +% figure; +% imshow(imb); title('Click On rough Center Posns : To return, click right button'); C = []; R = []; + %Initialis the cluster centroids + [R, C] = initClutsreCentres(bwImg, 'distT'); + load ('flDet'); flDet.remSet.R = R; flDet.remSet.C = C; save('flDet', 'flDet'); + numCells = length(R); + if (numCells > 0) + [X, Y] = meshgrid(1:SzC, 1:SzR); + MaxItrn = 1; FlImg = true(size(imb)); + for itrn = 1:MaxItrn + showTime = bwImg; showTime(sub2ind(size(bwImg), R,C)) = false; %figure; imshow(showTime); title('hi'); + % itrn + imb = doTheSmallestClosing(imb, C, R, 'line'); + for i = 1:numCells + % itrn + % [C R] + CGD = bwdistgeodesic(imb, C(i), R(i), 'quasi-euclidean'); + % [C R] + FGD = bwdistgeodesic(FlImg, C(i), R(i), 'quasi-euclidean'); + GD(:, :, i) = CGD + abs(CGD - FGD).*(0.75*CGD); + end + %Assign to cluster + [~, lbl] = min(GD, [], 3); + Lbls = lbl.*imb; + % figure; imshow(Lbls, []) + %Update Centroid + for i = 1:numCells + cLbl = (Lbls == i); + C(i) = mean(X(cLbl)); R(i) = mean(Y(cLbl)); + end + C = round(C); R = round(R); + % if (sum(isnan(C)) > 0) + % hld = 1; + % end + end + imwrite(showTime, 'showTime.jpg'); + + ClustImg = Lbls; + [elpseLbld, fitImg] = improveSheds(Lbls, imb); + end +end + +function [SegLbld, elpseLbld, fitImg] = kMeansByGeoDesicCentrWatershed(imgUint8, clearSet, remSet) + %Clean my binary image. I.e., remove spurious projections + dsk = strel('disk', 5); remSet = imdilate(imerode(remSet, dsk), dsk); + [R, C] = initClutsreCentres(remSet, 'distT'); + centRemSetRC = [R C]; + SegLbld = markrBasedWshed(imgUint8, clearSet, remSet, centRemSetRC); + %repeat once again for those missed + SegLbldN = doASecondPass(imgUint8, SegLbld, clearSet, remSet); + %Make a single seg image + SegLbld = SegLbld + SegLbldN; + [elpseLbld, fitImg] = improveSheds(SegLbld, remSet); +end + +function SegLbldN = doASecondPass(imgUint8, SegLbld, clearSet, remSet) + cmpLbl = 0; stDsk = strel('disk', 5); + firstPassSet = (SegLbld > 0); + clearSet = clearSet | firstPassSet; + remSet = remSet & ~firstPassSet; + %remove spurious projections + remSet = bwareaopen(imdilate(imerode(remSet, stDsk), stDsk), ceil(0.5*(41*41))); + [R, C] = initClutsreCentres(remSet, 'distT'); + centRemSetRC = [R C]; + SegLbldN = markrBasedWshed(imgUint8, clearSet, remSet, centRemSetRC); + %update the labels with respect to SegLbld + SegLbldN(SegLbldN > 0) = SegLbldN(SegLbldN > 0) + max(SegLbld(:)); +end +function SegLbld = markrBasedWshed(img, clearSet, remSet, centRemSetRC) + %%Make our image where watershed is to be performed. + img = rgb2gray(img); remSetImg = img; + dilremSet = bwmorph(remSet, 'dilate'); + remSetImg(~dilremSet) = 0; + remSetImg(clearSet) = 0; + %% + %%Select the watershed function. We use gradient + hy = fspecial('sobel'); hx = hy'; + Iy = imfilter(double(remSetImg), hy, 'replicate'); + Ix = imfilter(double(remSetImg), hx, 'replicate'); + gradmag = sqrt(Ix.^2 + Iy.^2); + + %%get the foreground markers + stdsk = strel('disk', 3); + centDil = false(size(remSet)); centDil(sub2ind(size(remSet), centRemSetRC(:, 1), centRemSetRC(:, 2))) = true; + fgm = imdilate(centDil, stdsk); + %% + %%get background markers + bw = remSet | clearSet; + D = bwdist(bw); + DL = watershed(D); + bgm = DL == 0; + %% + %%apply watershed + gradmagMdfd = imimposemin(gradmag, bgm | fgm); + L = watershed(gradmagMdfd); +% remSetImg(L == 0) = 0; + %% + %%Post processing, remove small extensions, merge adjacent small cells + SegLbld = postprocessing(L); +end + +function SegLbld = postprocessing(L) + %Remove all components which are bigger & smaller compared to 41*41 + cmpLbl = 0; stDsk = strel('disk', 3); + SegLbld = uint8(zeros(size(L))); + typArea = 41*41; maxArea = 1.5*typArea; minArea = 0.5*typArea; + numCmp = max(L(:)); + for i = 1:numCmp + currCmp = (L == i); areaCmp = sum(currCmp(:)); + if (areaCmp > minArea && areaCmp < maxArea) + cmpLbl = cmpLbl + 1; + %remove small extensions +% currCmp = bwmorph(currCmp, 'open', 5); + currCmp = bwareaopen(imdilate(imerode(currCmp, stDsk), stDsk), 100); + SegLbld(currCmp) = cmpLbl; + end + end +end + +function imb = doTheSmallestClosing(imb, C, R, by) + [SzM, SzN] = size(imb); + len = length(C); baseIm = false(size(imb)); + stDsk = strel('disk', 5); + for i = 1:len + if (imb(R(i), C(i))) + continue; + else + if (strcmp(by, 'line')) + len = 30; strt = floor(len/2); + lin000SL = -Inf; lin000EL = Inf; + lin045SL = -Inf; lin045EL = Inf; + lin090SL = -Inf; lin090EL = Inf; + lin135SL = -Inf; lin135EL = Inf; + + cR = R(i); cC = C(i); + %for zero degree line + flag = false; for j = 0:-1:-1*strt; if ((cC+j < 1) || (cC + j > SzN)); flag = false; break; end; if (imb(cR, cC+j) == false); indx = j; else; flag = true; break; end; end; if (flag); lin000SL = indx; end; + flag = false; for j = 0:1:strt; if ((cC+j < 1) || (cC + j > SzN)); flag = false; break; end; if (imb(cR, cC+j) == false); indx = j; else; flag = true; break; end; end; if (flag); lin000EL = indx; end; + %for 90 degree line + flag = false; for j = 0:-1:-1*strt; if ((cR+j < 1) || (cR+j > SzM)); flag = false; break; end; if (imb(cR+j, cC) == false); indx = j; else; flag = true; break; end; end; if (flag); lin090SL = indx; end; + flag = false; for j = 0:1:strt; if ((cR+j < 1) || (cR+j > SzM)); flag = false; break; end; if (imb(cR+j, cC) == false); indx = j; else; flag = true; break; end; end; if (flag); lin090EL = indx; end; + %for 45 degree line + flag = false; for j = 0:-1:-1*strt; if ((cC-j < 1) || (cC - j > SzN) || (cR+j < 1) || (cR+j > SzM)); flag = false; break; end; if (imb(cR+j, cC-j) == false); indx = j; else; flag = true; break; end; end; if (flag); lin045SL = indx; end; + flag = false; for j = 0:1:strt; if ((cC-j < 1) || (cC - j > SzN) || (cR+j < 1) || (cR+j > SzM)); flag = false; break; end; if (imb(cR+j, cC-j) == false); indx = j; else; flag = true; break; end; end; if (flag); lin045EL = indx; end; + %for 135 degree line + flag = false; for j = 0:-1:-1*strt; if ((cC+j < 1) || (cC + j > SzN) || (cR+j < 1) || (cR+j > SzM)); flag = false; break; end; if (imb(cR+j, cC+j) == false); indx = j; else; flag = true; break; end; end; if (flag); lin135SL = indx; end; + flag = false; for j = 0:1:strt; if ((cC+j < 1) || (cC + j > SzN) || (cR+j < 1) || (cR+j > SzM)); flag = false; break; end; if (imb(cR+j, cC+j) == false); indx = j; else; flag = true; break; end; end; if (flag); lin135EL = indx; end; + lngths = [lin000EL - lin000SL + 1; lin045EL - lin045SL + 1; lin090EL - lin090SL + 1; lin135EL - lin135SL + 1]; + [minVl, pos] = min(lngths); + if (~isfinite(minVl)) + imb(cR, cC) = true; + else + if (pos == 1) + %0 degree qualified + imb(cR, cC+lin000SL:cC+lin000EL) = true; + elseif (pos == 2) + %45 degree qualified +% imb(cR+lin045SL:cR+lin045EL, cC+lin045SL:cC+lin045EL) = true; + for k = lin045SL:lin045EL + imb(cR+k, cC - k) = true; + end + elseif (pos == 3) + %90 degree qualified + imb(cR+lin090SL:cR+lin090EL, cC) = true; + else + %135 dgree qualified +% imb(cR+lin135SL:cR+lin135EL, cC+lin135SL:cC+lin135EL) = true; + for k = lin135SL:lin135EL + imb(cR+k, cC + k) = true; + end + end + imb = bwfill(imb, 'holes'); + end + + elseif (strcmp(by, 'disk')) + tmpIm = baseIm; tmpIm(R(i), C(i)) = true; + while(true) + tmpIm = imdilate(tmpIm,stDsk); + stat = (tmpIm & imb); + if (sum(stat(:)) > 0) + imb = imb | tmpIm; + break; + end + end + end + imb = bwfill(imb, 'holes'); + end + end + +end +function [R, C] = initClutsreCentres(bwImg, Method) + R = []; C = []; flseIm = false(size(bwImg)); + if (strcmp(Method, 'manual')) + while (true) + [Cc, Rr, Btn] = ginput(1); + if (Btn ~= 1) + break; + end + C = [C; Cc]; R = [R; Rr]; + end + C = round(C); R = round(R); numCells = length(C); + elseif (strcmp(Method, 'random')) + lbls = bwlabel(bwImg); + cellArea = 37*37; + numObj = max(lbls(:)); + for i = 1:numObj + [Rs, Cs] = find(lbls == i); + objArea = length(Rs); + numCells = max(round(objArea/cellArea), 1); %We set 0.2 (500/(40*40)) +% numCells = 9; + slctdIndx = randperm(objArea, numCells); + R = [R; Rs(slctdIndx)]; + C = [C; Cs(slctdIndx)]; + end + elseif (strcmp(Method, 'distT')) + distTImg = bwdist(~bwImg); Radius = 12; + [KeyPoints, Rows, Cols] = doNonMaximaSuppression(distTImg, Radius); + KeyPoints = KeyPoints & bwImg; + Centroids = getCentroids(KeyPoints); + initKPImg = flseIm; indxKP = sub2ind(size(bwImg), Centroids(2, :), Centroids(1, :)); + initKPImg(indxKP) = true; + fnalKPImg = confirmCentroids(bwImg, initKPImg, Radius); + [R, C] = find(fnalKPImg); + end +end +function Centroids = getCentroids(KeyPoints) + Centroids = []; + propsCent = regionprops(KeyPoints, 'centroid'); + [numCmp, ~] = size(propsCent); + dp = [Centroids; propsCent.Centroid]; + Centroids = round(reshape(dp, 2, numCmp)); +end +function fnalKPImg = confirmCentroids(bwImg, initKPImg, Radius) + fnalKPImg = initKPImg; + typcalArea = 41*41; Radius = 0.75*Radius; + diskSt = strel('disk', Radius); + falseIm = false(size(bwImg)); + cmps = bwconncomp(bwImg); + numSeg = cmps.NumObjects; + for i = 1:numSeg + pxls = cmps.PixelIdxList{i}; + cArea = length(pxls); + cSeg = falseIm; cSeg(pxls) = true; + cCentr = cSeg & initKPImg; + numClsExpctd = max(round(cArea/typcalArea), 1); + numClsIdntfd = sum(cCentr(:)); + if (numClsExpctd > 1) + TBIdntfd = numClsExpctd - numClsIdntfd; + if (TBIdntfd > 0) + cCentrD = imdilate(cCentr, diskSt); + distImage = bwdist(~(~cCentrD & cSeg)); + [KeyPoints, Rows, Cols] = doNonMaximaSuppression(distImage, Radius); + Centroids = getCentroids(KeyPoints); + indxKP = sub2ind(size(bwImg), Centroids(2, :), Centroids(1, :)); + [distnCs, posns] = sort(distImage(indxKP), 'descend'); + TBIdntfd = min(TBIdntfd, length(posns)); + indxUpdtd = indxKP(posns(1:TBIdntfd)); + fnalKPImg(indxUpdtd) = true; + end + end + end +end +function applyMarkrBasedWatershed(img, ROIntst, foreGnd, backGnd) + +end +function ClustImg = waterShedBasedOnGeoDesicDistNew(img) +ClustImg = []; +close all; + imb = img; %img(22:129, 393:525); %img; % + imshow(imb); + [SzR, SzC] = size(imb); + imagesc(imb); title('Click On rough Center Posns : To return, click right button'); C = []; R = []; + while (true) + [Cc, Rr, Btn] = ginput(1); + if (Btn ~= 1) + break; + end + C = [C Cc]; R = [R Rr]; + end + C = round(C); R = round(R); numCells = length(C); + imbc = imb; imbc(R(1), C(1)) = 0; imbc(R(2), C(2)) = 0; imbc(R(3), C(3)) = 0; + imagesc(imbc); + + [X, Y] = meshgrid(1:SzC, 1:SzR); + FlImg = true(size(imb)); + figure; + Cm = bwdistgeodesic(imb, C, R, 'quasi-euclidean'); + Cm(isnan(Cm)) = Inf; + lb = watershed(Cm); + figure; imshow(lb .* uint8(imb), []); + figure; imshow(lb == 0); + fnlLbl = improveSheds(lb, imb) +end + +function [elpseLbld, fitImg] = improveSheds(Lbls, img) + numCells = max(Lbls(:)); + [SzR, SzC] = size(img); + [X, Y] = meshgrid(1:SzC, 1:SzR); + FnalImg = img; + for i = 1:numCells + cBin = (Lbls == i); +% figure; imshow(cBin); + FnalImg(bwperim(cBin)) = false; + %Compute Major and Minor Axis of this + end + props = regionprops(FnalImg, {'MajorAxisLength', ... + 'MinorAxisLength', 'Orientation', 'Centroid'}); + [elpseLbld, fitImg] = getTheEllipses(props, img); +end + +function [ellipses, fitImg] = getTheEllipses(props, img) +close all; +% imshow(img); +[MxR, MxC] = size(img); + [numEllpse, ~] = size(props); + falsEImg = false(size(img)); + ellipses = zeros(size(img)); + for i = 1:numEllpse + CCC = props(i).Centroid(1); CCR = props(i).Centroid(2); + CMJ = round(props(i).MajorAxisLength); + CMN = round(props(i).MinorAxisLength); + CON = props(i).Orientation; + %[X, Y] = meshgrid(1:CMN, 1:CMJ); + + if (abs(CON) >= 45) + X = repmat(-CMN/2:CMN/(CMN-1):CMN/2, CMJ, 1); + Y = repmat((-CMJ/2:CMJ/(CMJ-1):CMJ/2)', 1, CMN); + a = CMJ/2; b = CMN/2; + %vertical eelipse + el = (Y/a).^2 + (X/b).^2 < 1; + else + Y = repmat(-CMN/2:CMN/(CMN-1):CMN/2, CMJ, 1); + X = repmat((-CMJ/2:CMJ/(CMJ-1):CMJ/2)', 1, CMN); + a = CMJ/2; b = CMN/2; + el = (X/a).^2 + (Y/b).^2 < 1; + end +% imshow(el); + rotImg = imrotate(el, CON+90); +% imshow(rotImg); + %placeAtCentroid + [NSr, NSc] = size(rotImg); + %This is the minimum sampling. Mke it higher + CX = repmat(-NSc/2:NSc/(NSc-1):NSc/2, NSr, 1); + CY = repmat((-NSr/2:NSr/(NSr-1):NSr/2)', 1, NSc); + + CIX = CX(rotImg); CIY = CY(rotImg); + CIXs = round(CIX + CCC); CIYs = round(CIY + CCR); + outSub = ((CIYs > MxR) | (CIXs > MxC) | (CIYs < 1) | (CIXs < 1)); + CIYs(outSub) = []; CIXs(outSub) = []; + indX = sub2ind(size(falsEImg), CIYs, CIXs); + falsEImg(indX) = true; + ellipses(indX) = i; + end + ellipses = medfilt2(ellipses); + fitImg = ellipses > 0; +end +function [remMsk, Mskd] = maskAllSharingObjects(BaseImg, ShareObjs) + ToBMaskd = BaseImg & ShareObjs; + LblsFrmBase = bwlabel(BaseImg); + Mskd = false(size(BaseImg)); + lbls2BMaskd = LblsFrmBase(ToBMaskd); + unqLbls = unique(lbls2BMaskd(:)); + for lbl = 1:length(unqLbls) + Mskd(LblsFrmBase == (unqLbls(lbl))) = true; + end + remMsk = BaseImg & ~Mskd; +end +function [KeyPoints, Rows, Cols] = doNonMaximaSuppression(distImage, Radius) + Threshold = Radius; + SizeOfMask = 2*Radius+1; % Size of mask. + Max = ordfilt2(distImage,SizeOfMask^2,ones(SizeOfMask)); % Grey-scale dilate. + KeyPoints = (distImage == Max)&(distImage > Threshold); % Find maxima. + [Rows, Cols] = find(KeyPoints); +end \ No newline at end of file diff --git a/src/gopa kumar code/initializeCharacterCNNBNNEx6ColorFczRGB.m b/src/gopa kumar code/initializeCharacterCNNBNNEx6ColorFczRGB.m new file mode 100644 index 0000000..1f52ef2 --- /dev/null +++ b/src/gopa kumar code/initializeCharacterCNNBNNEx6ColorFczRGB.m @@ -0,0 +1,37 @@ +function net = initializeCharacterCNNBNNEx6ColorFczRGB() + +f=1/100 ; +net.layers = {} ; +net.layers{end+1} = struct('type', 'conv', ... + 'filters', f*randn(5,5,9,20, 'single'), ... + 'biases', zeros(1, 20, 'single'), ... + 'stride', 1, ... + 'pad', 0) ; +net.layers{end+1} = struct('type', 'relu') ; +net.layers{end+1} = struct('type', 'pool', ... + 'method', 'max', ... + 'pool', [2 2], ... + 'stride', 2, ... + 'pad', 0) ; +net.layers{end+1} = struct('type', 'conv', ... + 'filters', f*randn(5,5,20,50, 'single'),... + 'biases', zeros(1,50,'single'), ... + 'stride', 1, ... + 'pad', 0) ; +net.layers{end+1} = struct('type', 'pool', ... + 'method', 'max', ... + 'pool', [2 2], ... + 'stride', 2, ... + 'pad', 0) ; +net.layers{end+1} = struct('type', 'conv', ... + 'filters', f*randn(4,4,50,500, 'single'),... + 'biases', zeros(1,500,'single'), ... + 'stride', 1, ... + 'pad', 0) ; +net.layers{end+1} = struct('type', 'relu') ; +net.layers{end+1} = struct('type', 'conv', ... + 'filters', f*randn(2,2,500,2, 'single'),... + 'biases', zeros(1,2,'single'), ... + 'stride', 1, ... + 'pad', 0) ; +net.layers{end+1} = struct('type', 'softmaxloss') ; diff --git a/src/gopa kumar code/makeDataSetForTrainingBest.m b/src/gopa kumar code/makeDataSetForTrainingBest.m new file mode 100644 index 0000000..40330e9 --- /dev/null +++ b/src/gopa kumar code/makeDataSetForTrainingBest.m @@ -0,0 +1,310 @@ +function makeDataSetForTrainingBest() + clear all; clc; close all; + params.PstvTnVnTtRaio = [60 20 20]; + params.DatasetLoc = 'CNNPatches\AllDataSet\'; + params.typeStrt = 1; params.typeEnd = 2; + params.wbc = true; + params.dst = false; + params.ngtvProp = 0.8; + params.wbcProp = 0.2; + params.dstProp = 0; + + if (params.wbc) + if (params.dst) + params.DatasetSve = 'CNNPatches\WithDstWBC\'; + else + params.DatasetSve = 'CNNPatches\WithWBCNoDst\'; + end + else + if (params.dst) + params.DatasetSve = 'CNNPatches\WithDstNoWBC\'; + else + params.DatasetSve = 'CNNPatches\WithOutDstWBC\'; + end + end + params.isForCNN = true; params.cpy = false; + IdString = 'FczRGB'; makeDataSetCNNTrain(IdString, params); + display('First Pass'); + params.isForCNN = true; params.cpy = true; params.cpyString = 'FczRGB'; + IdString = 'BFczdRGB'; makeDataSetCNNTrain(IdString, params); + display('Second Pass'); + params.isForCNN = false; params.cpy = true; params.cpyString = 'FczRGB'; + IdString = 'BFczdFeatRGB'; makeDataSetCNNTrain(IdString, params); + display('Third Pass'); +end +function makeDataSetCNNTrain(IdString, params) + fnames = dir([params.DatasetLoc IdString 'Ps*Type1*.mat']); + numPstvStk(1) = length(fnames); + fnames = dir([params.DatasetLoc IdString 'Ps*Type2*.mat']); + numPstvStk(2) = length(fnames); + + fnames = dir([params.DatasetLoc IdString 'Ng*Type1*.mat']); + numNgtvStk(1) = length(fnames); + fnames = dir([params.DatasetLoc IdString 'Ng*Type2*.mat']); + numNgtvStk(2) = length(fnames); + + fnames = dir([params.DatasetLoc IdString 'Wc*Type1*.mat']); + numWbcStk(1) = length(fnames); + fnames = dir([params.DatasetLoc IdString 'Wc*Type2*.mat']); + numWbcStk(2) = length(fnames); + + fnames = dir([params.DatasetLoc IdString 'Ds*Type1*.mat']); + numDstStk(1) = length(fnames); + fnames = dir([params.DatasetLoc IdString 'Ds*Type2*.mat']); + numDstStk(2) = length(fnames); + + %Process positive stack + numPstvImgs = 0; + for i = params.typeStrt:params.typeEnd + load ([params.DatasetLoc IdString 'PsStk_' 'Type' num2str(i)]); + if (params.isForCNN) + ptvStack = augmentDataSetByRotation(postvStack); + [~, ~, ~, cPstvCnt] = size(ptvStack); + POSTVSTACK(:, :, :, numPstvImgs+1: numPstvImgs+cPstvCnt) = ptvStack; + else + ptvStack = augmentDataSetSimply(postvStack); + [cPstvCnt, ~] = size(ptvStack); + POSTVSTACK(numPstvImgs+1: numPstvImgs+cPstvCnt, :) = ptvStack; + end + numPstvImgs = numPstvImgs+cPstvCnt; + end + clear ptvStack; clear postvStack; + + PstvTrain = params.PstvTnVnTtRaio(1); + PstvValdn = params.PstvTnVnTtRaio(2); + PstvTest = params.PstvTnVnTtRaio(3); + numSamPTn = round(numPstvImgs*PstvTrain/100); + numSamPVn = round(numPstvImgs*PstvValdn/100); + numSamPTt = numPstvImgs - (numSamPTn + numSamPVn); + PtvTnVnTt = [PstvTrain PstvValdn PstvTest]; + + %negative can be many and hence read in batches + numNgtvImgs = 0; + for i = params.typeStrt:params.typeEnd + for j = 1:numNgtvStk(i) + load ([params.DatasetLoc IdString 'NgStk_' 'Type' num2str(i) '_' num2str(j)]); + if (params.isForCNN) + [~, ~, ~, currNcnt] = size(negtvStack); + else + [currNcnt, ~] = size(negtvStack); + end + numNgtvImgs = numNgtvImgs + currNcnt; + end + end + + numSamNTn = round(numSamPTn * params.ngtvProp); + numSamNVn = round(numSamPVn * params.ngtvProp); + numSamNTt = numNgtvImgs - (numSamNTn + numSamNVn); + NgtvTrain = (numSamNTn/numNgtvImgs)*100; + NgtvVlidn = (numSamNVn/numNgtvImgs)*100; + NgtvTest = 100 - (NgtvTrain + NgtvVlidn); + NtvTnVnTt = [NgtvTrain NgtvVlidn NgtvTest]; + + + if (params.wbc) + numWbcImgs = 0; + for i = params.typeStrt:params.typeEnd + load ([params.DatasetLoc IdString 'WcStk_' 'Type' num2str(i)]); + if (params.isForCNN) + [~, ~, ~, cWbcCnt] = size(wbcStack); + WBCSTACK(:, :, :, numWbcImgs+1: numWbcImgs+cWbcCnt) = wbcStack; + else + [cWbcCnt, ~] = size(wbcStack); + WBCSTACK(numWbcImgs+1: numWbcImgs+cWbcCnt, :) = wbcStack; + end + numWbcImgs = numWbcImgs+cWbcCnt; + end + clear wbcStack; + + numSamWTn = round(numSamPTn * params.wbcProp); + numSamWVn = round(numSamPVn * params.wbcProp); + numSamWTt = numNgtvImgs - (numSamWTn + numSamWVn); + WbcTrain = (numSamWTn/numWbcImgs)*100; + WbcVlidn = (numSamWVn/numWbcImgs)*100; + WbcTest = 100 - (WbcTrain + WbcVlidn); + WbcTnVnTt = [WbcTrain WbcVlidn WbcTest]; + else + WBCSTACK = []; numWbcImgs = 0; + numSamWTn = 0; numSamWVn = 0; numSamWTt = 0; + WbcTnVnTt = [0 0 0]; + end + + if (params.dst) + numDstImgs = 0; + for i = params.typeStrt:params.typeEnd + load ([params.DatasetLoc IdString 'DsStk_' 'Type' num2str(i)]); + if (params.isForCNN) + [~, ~, ~, cDstCnt] = size(dustStack); + DUSTSTACK(:, :, :, numDstImgs+1: numDstImgs + cDstCnt) = dustStack; + else + [cDstCnt, ~] = size(dustStack); + DUSTSTACK(numDstImgs+1: numDstImgs + cDstCnt, :) = dustStack; + end + numDstImgs = numDstImgs+cDstCnt; + end + clear dustStack; + + numSamDTn = round(numSamPTn * params.dstProp); + numSamDVn = round(numSamPVn * params.dstProp); + numSamDTt = numDstImgs - (numSamDTn + numSamDVn); + DstTrain = (numSamDTn/numDstImgs)*100; + DstVlidn = (numSamDVn/numDstImgs)*100; + DstTest = 100 - (DstTrain + DstVlidn); + DstTnVnTt = [DstTrain DstVlidn DstTest]; + else + DUSTSTACK = []; numDstImgs = 0; + numSamDTn = 0; numSamDVn = 0; numSamDTt = 0; + DstTnVnTt = [0 0 0]; + end + + %Create a full Databas clear all + if (~(params.cpy)) + FullImdb.meta.classes = 'Malaria,Healthy'; + FullImdb.meta.sets = {'train', 'val', 'test'}; + FullImdb.meta.infn = ['Patches Cells : ' IdString]; + FullImdb.images.id = 1:(numPstvImgs+numNgtvImgs+numWbcImgs+numDstImgs); + %Not Storing the data as it is huge. May be difficult to load + FullImdb.images.label = [ones(1, numPstvImgs) 2*ones(1, numNgtvImgs+numWbcImgs+numDstImgs)]; + %Identify and store the train, validn and Test indexes + + pstvSplit = divideDataInRatio(numPstvImgs, PtvTnVnTt); + ngtvSplit = divideDataInRatio(numNgtvImgs, NtvTnVnTt); + wbcSplit = divideDataInRatio(numWbcImgs, WbcTnVnTt); + dstSplit = divideDataInRatio(numDstImgs, DstTnVnTt); + + FullImdb.images.set = [pstvSplit ngtvSplit wbcSplit dstSplit]; + save([params.DatasetSve IdString '_FullImdb'], 'FullImdb'); + else + load ([params.DatasetSve params.cpyString '_FullImdb'], 'FullImdb'); + end + + %Select dataset for training + FllLbl = FullImdb.images.label; + lbelCatgry = FullImdb.images.set; + slectdFrTraining = (lbelCatgry == 1) | (lbelCatgry == 2); + + MalImdb.meta.classes = 'Malaria,Healthy'; + MalImdb.meta.sets = {'train', 'val'}; + MalImdb.images.id = 1:sum(slectdFrTraining); + MalImdb.images.label = FllLbl(slectdFrTraining); + MalImdb.images.set = lbelCatgry(slectdFrTraining); + + pstvTrainRegion = slectdFrTraining(1:numPstvImgs); + if (params.isForCNN) + pstVDataSlctd = POSTVSTACK(:, :, :, pstvTrainRegion); %we did not clear the positive set + [~, ~, ~, cntPstvTnSlctd] = size(pstVDataSlctd); + else + pstVDataSlctd = POSTVSTACK(pstvTrainRegion, :); %we did not clear the positive set + [cntPstvTnSlctd, ~] = size(pstVDataSlctd); + end + + nextStart = 1+numPstvImgs; + NgStart = 1; + for i = params.typeStrt:params.typeEnd + for j = 1:numNgtvStk(i) + load ([params.DatasetLoc IdString 'NgStk_' 'Type' num2str(i) '_' num2str(j)]); + if (params.isForCNN) + [~, ~, ~, currNcnt] = size(negtvStack); + else + [currNcnt, ~] = size(negtvStack); + end + + ngtvTrainRegion = slectdFrTraining(nextStart:nextStart-1+currNcnt); + nextStart = nextStart+currNcnt; + currSegmntNgtvs = sum(ngtvTrainRegion); + NgEnd = NgStart-1+currSegmntNgtvs; + if (params.isForCNN) + ngtVDataSlctd(:, :, :, NgStart:NgEnd) = negtvStack(:, :, :, ngtvTrainRegion); + else + ngtVDataSlctd(NgStart:NgEnd, :) = negtvStack(ngtvTrainRegion, :); + end + NgStart = NgEnd+1; + end + end + if (params.isForCNN) + [~, ~, ~, cntNgtvTnSlctd] = size(ngtVDataSlctd); + else + [cntNgtvTnSlctd, ~] = size(ngtVDataSlctd); + end + + %Now select WBC data + nextStart = 1+numPstvImgs+numNgtvImgs; + wbcTrainRegion = slectdFrTraining(nextStart:nextStart-1+numWbcImgs); + if (params.isForCNN) + wbcDataSlctd = WBCSTACK(:, :, :, wbcTrainRegion); + [~, ~, ~, cntWbcTnSlctd] = size(wbcDataSlctd); + else + wbcDataSlctd = WBCSTACK(wbcTrainRegion, :); + [cntWbcTnSlctd, ~] = size(wbcDataSlctd); + end + + %Now select Dst data + nextStart = 1+numPstvImgs+numNgtvImgs+numWbcImgs; + dstTrainRegion = slectdFrTraining(nextStart:nextStart-1+numDstImgs); + if (params.isForCNN) + dstDataSlctd = DUSTSTACK(:, :, :, dstTrainRegion); + [~, ~, ~, cntDstTnSlctd] = size(dstDataSlctd); + else + dstDataSlctd = DUSTSTACK(dstTrainRegion, :); + [cntDstTnSlctd, ~] = size(dstDataSlctd); + end + + if (params.isForCNN) + strt = 1; endLc = cntPstvTnSlctd; + FllData(:, :, :,strt:endLc) = pstVDataSlctd; + strt = endLc+1; endLc = strt-1+cntNgtvTnSlctd; + FllData(:, :, :,strt:endLc) = ngtVDataSlctd; + strt = endLc+1; endLc = strt-1+cntWbcTnSlctd; + FllData(:, :, :,strt:endLc) = wbcDataSlctd; + strt = endLc+1; endLc = strt-1+cntDstTnSlctd; + FllData(:, :, :,strt:endLc) = dstDataSlctd; + else + strt = 1; endLc = cntPstvTnSlctd; + FllData(strt:endLc, :) = pstVDataSlctd; + strt = endLc+1; endLc = strt-1+cntNgtvTnSlctd; + FllData(strt:endLc, :) = ngtVDataSlctd; + strt = endLc+1; endLc = strt-1+cntWbcTnSlctd; + FllData(strt:endLc, :) = wbcDataSlctd; + strt = endLc+1; endLc = strt-1+cntDstTnSlctd; + FllData(strt:endLc, :) = dstDataSlctd; + end + + MalImdb.images.data = FllData; + save ([params.DatasetSve IdString '_MalImdb'], 'MalImdb'); +end + + +function TnTtVn = divideDataInRatio(MaxLt, Ratio) + blck = randperm(MaxLt); strt = 1; + cmSum = cumsum(Ratio); + TnTtVn = zeros(1, MaxLt); + for i = 1:length(Ratio) + intstd = round(cmSum(i)/100*MaxLt); + TnTtVn((blck >= strt) & (blck <= intstd)) = i; + strt = intstd+1; + end +end +function augStack = augmentDataSetByRotation(postvStack) + [SzM, SzN, SzO, numImgs] = size(postvStack); + augStack = single(zeros(SzM, SzN, SzO, 4*numImgs)); + augCnt = 0; + for i = 1:numImgs + currIm = postvStack(:, :, :, i); + augStack(:, :, :, augCnt+1) = currIm; %0 degree + augStack(:, :, :, augCnt+2) = rot90(currIm, 1); %90 + augStack(:, :, :, augCnt+3) = rot90(currIm, 2); %180 + augStack(:, :, :, augCnt+4) = rot90(currIm, 3); %270 + augCnt = augCnt + 4; + end +end + +function augStack = augmentDataSetSimply(postvStack) + [numImgs, FeatD] = size(postvStack); + augStack = single(zeros(4*numImgs, FeatD)); + augCnt = 0; + for i = 1:numImgs + currIm = postvStack(i, :); + augStack(augCnt+1:augCnt+4, :) = [currIm; currIm;currIm;currIm]; %0 degree + augCnt = augCnt + 4; + end +end \ No newline at end of file diff --git a/src/gopa kumar code/reName.m b/src/gopa kumar code/reName.m new file mode 100644 index 0000000..f3dc731 --- /dev/null +++ b/src/gopa kumar code/reName.m @@ -0,0 +1,30 @@ +function reName() +clear all; clc; close all; + IdString = 'FczRGB'; + Categry = 'Ds'; % 'Wc'; % 'Ng' % 'Ps'; % % % + DatasetLoc = 'CNNPatches\AllDataSet\'; + fnames = dir([DatasetLoc IdString '*' Categry '*.mat']); + [cnt, ~] = size(fnames); + for i = 1:cnt + load ([DatasetLoc fnames(i).name]); + if (strcmp(Categry, 'Ps')) + postvStack = postvStackFcsRGB; + save ([DatasetLoc fnames(i).name], 'postvStack'); + clear postvStack; + elseif (strcmp(Categry, 'Ng')) + negtvStack = negtvStackFcsRGB; + save ([DatasetLoc fnames(i).name], 'negtvStack'); + clear negtvStack; + elseif (strcmp(Categry, 'Wc')) + wbcStack = WBCStackFcsRGB ; + save ([DatasetLoc fnames(i).name], 'wbcStack'); + clear wbcStack; + elseif (strcmp(Categry, 'Ds')) + dustStack = DustStackFcsRGB; + save ([DatasetLoc fnames(i).name], 'dustStack'); + clear dustStack; + else + df = 1; + end + end +end \ No newline at end of file diff --git a/src/gopa kumar code/setup.m b/src/gopa kumar code/setup.m new file mode 100644 index 0000000..af88bcc --- /dev/null +++ b/src/gopa kumar code/setup.m @@ -0,0 +1,25 @@ +function setup(varargin) +FeetBase = 'E:\Gopakumar\Dataset IIST RSDAY\IIST RSDAY MY WORK\Work\Code\matconvnet-1.0-beta11\practical-cnn-2015a\practical-cnn-2015a'; +run ([FeetBase '\vlfeat\toolbox\vl_setup']); +run ([FeetBase '\matconvnet\matlab\vl_setupnn']); +addpath ([FeetBase '\matconvnet\examples']); + +opts.useGpu = false ; +opts.verbose = false ; +opts = vl_argparse(opts, varargin) ; + +try + vl_nnconv(single(1),single(1),[]) ; +catch + warning('VL_NNCONV() does not seem to be compiled. Trying to compile it now.') ; + vl_compilenn('enableGpu', opts.useGpu, 'verbose', opts.verbose) ; +end + +if opts.useGpu + try + vl_nnconv(gpuArray(single(1)),gpuArray(single(1)),[]) ; + catch + vl_compilenn('enableGpu', opts.useGpu, 'verbose', opts.verbose) ; + warning('GPU support does not seem to be compiled in MatConvNet. Trying to compile it now') ; + end +end diff --git a/src/gopa kumar code/testCNNBNNEx6Color.m b/src/gopa kumar code/testCNNBNNEx6Color.m new file mode 100644 index 0000000..ecbbeb5 --- /dev/null +++ b/src/gopa kumar code/testCNNBNNEx6Color.m @@ -0,0 +1,13 @@ +function Lbls = testCNNBNNEx6Color(testImdb, net) + numImages = length(testImdb.images.id); + clasfdLb = zeros(1, numImages); + for i = 1:numImages + im = testImdb.images.data(:, :, :, i); + im = 256 * (im - net.imageMean) ; + % Apply the CNN to the larger image + res = vl_simplenn(net, im); + [~, clsfd] = max(res(end).x); + clasfdLb(i) = clsfd; + end + Lbls = [clasfdLb; testImdb.images.label]; +end \ No newline at end of file diff --git a/src/gopa kumar code/testTrainedCNNSVMAutoCnt.m b/src/gopa kumar code/testTrainedCNNSVMAutoCnt.m new file mode 100644 index 0000000..dc38206 --- /dev/null +++ b/src/gopa kumar code/testTrainedCNNSVMAutoCnt.m @@ -0,0 +1,582 @@ +function testTrainedCNNSVMAutoCnt(params) %bkUpName, NetId, isModelSvm) +lstAccsd = 0; featCnt = 0; +alreadySegmented = params.alreadySegmented; +useMyMdlFrm = params.UseMdlFrm; +if (alreadySegmented) + load('CNNPatches\Segment\alreadySegImg'); + load('CNNPatches\Segment\alreadyWBCImg'); +end + %%This function reads each image, Identify patches 32x32 where there is + %%a chance of parasites (Regional Minima) and then test the location + %%for the possible parasite and mark it on the slide if it is there. + TPSTV = 0; FPSTV = 0; FNGTV = 0; uWntToCnt = true; diFCnt = 0; uwntToPrint = true; + wrngCnt = 0; WrngImExclDiffclt =single([]); WrngTrtExclDiffclt = single([]); LblsWrng = uint8([]); + TrthDifcltStack = logical([]); ImgsDifcltStack = single([]); LblsDifcltStack = uint8([]); + filtrOutDifclt = false; + if (uWntToCnt) + TotStat.Tp = 0; TotStat.Fn = 0; + TotStat.Fp = 0; TotStat.Tn = 0; + IndStat.stat = []; + end + + bkUpName = params.dataSet; + trainedPath = ['myExp6\' useMyMdlFrm '\' bkUpName '_MalImdb\']; +% ['E:\Gopakumar\Dataset IIST RSDAY\IIST RSDAY MY WORK\Work\' ... +% 'Code\matconvnet-1.0-beta11\practical-cnn-2015a\practical-cnn-2015a\myExp6\60_20_20_60_20_20']; +% = 60; + if (params.isSVMModel) + load ([trainedPath 'SvMModel.mat']); + SVMModel = SvMModel.model; + NetId = 0; + if (params.loadFeat) + load ('CNNPatches\Features\saveFeat'); + end + else + %OK. Run the set up for MATCONVNET ENVIRONMENT + setup ; + NetId = params.NetId; + %%%%load the trained net and make it ready for testing + load ([trainedPath 'net-epoch-' num2str(NetId)]); + load ([trainedPath 'imageMean']); + net.layers(end) = [] ; net.imageMean = imageMean ; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + end + %Get the minima R,G,B across stack + BasePath = ''; + fnames = dir([BasePath 'ImgesMinAcrsStck\*B.jpg']); + fczPath = 'AllFcsStack\'; + + numfids = length(fnames); PtchSzM = 32; PtchSzN = 32; + prevVidIndx = ''; + fh = figure; + for K = 1:numfids +% tic + K + numPstv = 0; +% fnames(K).name + AbsFNme = [BasePath 'ImgesMinAcrsStck\' fnames(K).name]; + %Get details needed to acces the Ground Truth file + [~, FileName, ~] = fileparts(AbsFNme); + for i = 1:length(FileName); if (FileName(i) == '_'); vid = i-1; break; end; end; + for j = vid+2:length(FileName); if (FileName(j) == '_'); stck = j-1; break; end; end; + vidIndx = FileName(1:vid); stckIndx = FileName(vid+2:stck); + load ([BasePath 'GndTrth\' vidIndx '_' stckIndx]); + if (~strcmp(prevVidIndx, vidIndx)); + prevVidIndx = vidIndx; + cDstLoc = getDstLocnsFor(vidIndx); + end + %Read the image and Make it Single precision + imgUint8 = imread([fczPath vidIndx '_' stckIndx '_21.jpg']); imSingle = im2single(imgUint8); + impUint8 = imread([fczPath vidIndx '_' stckIndx '_5.jpg']); impSingle = im2single(impUint8); + imnUint8 = imread([fczPath vidIndx '_' stckIndx '_37.jpg']); imnSingle = im2single(imnUint8); + if (strcmp(bkUpName, 'BFczdRGB')) + imGCombinedFcz(:, :, 1:3) = imSingle; + else + imGCombinedFcz(:, :, 1:3) = impSingle; imGCombinedFcz(:, :, 4:6) = imSingle; imGCombinedFcz(:, :, 7:9) = imnSingle; + end + + clear impUint8; clear impSingle; clear imnUint8; clear imnSingle; +% tic; + if (~alreadySegmented) + [segmentedImg, segWithBndry, WBCs] = getSegmentation(imgUint8); + alreadySegImg(:, :, K) = segmentedImg; alreadyWBCImg(:, :, K) = WBCs; + if (mod(K, 50) == 0) + save('CNNPatches\Segment\alreadySegImg', 'alreadySegImg'); + save('CNNPatches\Segment\alreadyWBCImg', 'alreadyWBCImg'); + end + else + segmentedImg = alreadySegImg(:, :, K); + WBCs = alreadyWBCImg(:, :, K); + end +% toc; + RC = imgUint8(:, :, 1); GC = imgUint8(:, :, 2); BC = imgUint8(:, :, 3); + [SzM, SzN, ~] = size(imgUint8); LblMsk = false(SzM, SzN); + %Compute regional Minima + [regMin, bgndMsk] = getMyRegionalMinima(imgUint8); + %Note that regMin is already filtered wrt bgnd + %This is for SVM + cellsImg = false(size(Msk)); + slideName = [vidIndx '_' stckIndx]; + + %Filter out regional minima at the location of WBC s + regMin = regMin & ~WBCs; + + %If postv to be checkd closer + chkClsr = cDstLoc & regMin; + + %Now test Each suspected Locns + Cntrids = regionprops(regMin, 'centroid'); + [numPosns, ~] = size(Cntrids); + vlidPosCnt = 0; + for posCnt = 1:numPosns + %For the time being, if there is a patch of the required size + %around the point, then only we are considering it. + CCntrids = round(Cntrids(posCnt).Centroid); + CR = CCntrids(2); CC = CCntrids(1); + minRw = CR - PtchSzM/2; minCl = CC - PtchSzN/2; + maxRw = CR + PtchSzM/2-1; maxCl = CC + PtchSzN/2-1; + if (minRw > 0 && minCl > 0 && maxRw <= SzM && maxCl <= SzN) + vlidPosCnt = vlidPosCnt + 1; + %Get curr patch + currPatch = imSingle(minRw:maxRw, minCl:maxCl, :); + iptPatch = imGCombinedFcz(minRw:maxRw, minCl:maxCl, :); + if (filtrOutDifclt) + tmpDel = currPatch; + TthMskPatch = Msk(minRw:maxRw, minCl:maxCl, :); + end + %Test it and Get the label + if (params.isSVMModel) + if(params.loadFeat) + [FeatSet, lstAccsd] = loadFeat(saveFeat, [vidIndx '_' stckIndx], CR, CC, lstAccsd); + else + cellsImg (minRw:maxRw, minCl:maxCl) = true; + [FeatSet, ~] = getMyFeatures(cellsImg, imgUint8, slideName); + cellsImg (minRw:maxRw, minCl:maxCl) = false; + featCnt = featCnt+1; + saveFeat(featCnt).FeatSet = FeatSet; + saveFeat(featCnt).FeatName = [vidIndx '_' stckIndx]; + saveFeat(featCnt).FeatR = CR; + saveFeat(featCnt).FeatC = CC; + end + Lbl = svmclassify(SVMModel, FeatSet); + else + iptPatch = 256*(iptPatch - net.imageMean); + res = vl_simplenn(net, iptPatch); + [~, Lbl] = max(res(end).x); + end + + + if (filtrOutDifclt) + rsp = res(end).x; + vld = abs(rsp(:, :, 1) - rsp(:, :, 2)); + if (vld < 3.0) + diFCnt = diFCnt + 1; + TrthDifcltStack (:, :, diFCnt) = TthMskPatch; + ImgsDifcltStack (:, :, :, diFCnt) = tmpDel; + LblsDifcltStack (:, diFCnt) = Lbl; + else + %Irrespective of Dst Loc + ofst = 5; + cntrPatch = TthMskPatch(PtchSzM/2-ofst:PtchSzM/2+ofst, PtchSzN/2-ofst:PtchSzN/2+ofst); + fpWrng = (Lbl == 1) && (sum(cntrPatch(:)) == 0); %flase pstv wrong + fnWrng = (Lbl == 2) && (sum(cntrPatch(:)) ~= 0); %flase ngtv wrong + if (fpWrng || fnWrng) + wrngCnt = wrngCnt + 1; + WrngTrtExclDiffclt(:, :, wrngCnt) = TthMskPatch; + WrngImExclDiffclt (:, :, :, wrngCnt) = tmpDel; + LblsWrng(:,wrngCnt) = Lbl; + end + end + end + + %If the label is postive, closely examine if it is dust? + if (Lbl == 1) + LblMsk(CR, CC) = true; + end + end + end + %Now Examine the Msks with the Mask of Dst Pos we have + LblMskDil = bwmorph(LblMsk, 'dilate', 3); + ToBExmnd = (LblMskDil & chkClsr); + [~, ToBExmnd] = maskAllSharingObjects(LblMskDil, ToBExmnd); + [~, ToBExmnd] = maskAllSharingObjects(LblMsk, ToBExmnd); + %If the average pixel intensity in the region 7x7 is not below 128 + %Exclude it from Positive + [clsR, clsC] = find(ToBExmnd); + for cntCls = 1:length(clsR) + cClsR = clsR(cntCls); cClsC = clsC(cntCls); + patchClse = GC(cClsR-3:cClsR+3, cClsC-3:cClsC+3); + meanPatch = mean(patchClse(:)); + if (meanPatch > 100) + LblMsk(cClsR , cClsC) = false; + end + end + %Show + LD = bwmorph(LblMsk, 'dilate', 4); + RC(LD) = 255; GC(LD) = 0; BC(LD) = 0; + RGBCmbnd1(:, :, 1) = RC; RGBCmbnd1(:, :, 2) = GC; RGBCmbnd1(:, :, 3) = BC; + subplot(1, 2, 1); imshow(RGBCmbnd1); title('Detected'); + + %Mask out the area that we did not consider ; boundry and bgnd + exclRfst = ceil(PtchSzM/2); exclCfst = ceil(PtchSzN/2); + MskExclded = false(size(Msk)); MskExclded(exclRfst:SzM-exclRfst, exclCfst:SzN-exclCfst) = true; + remFrmTruth = Msk & MskExclded & bwmorph(bwfill(bgndMsk, 'holes'), 'dilate', 3); %lets consider a little more towards the wall + TrthDidNtConsder = ~remFrmTruth & Msk; + + MD = bwmorph(remFrmTruth, 'dilate', 4); RD = bwmorph(TrthDidNtConsder, 'dilate', 4); + + RC = imgUint8(:, :, 1); GC = imgUint8(:, :, 2); BC = imgUint8(:, :, 3); + RC(MD) = 0; GC(MD) = 255; BC(MD) = 0; RC(RD) = 0; GC(RD) = 0; BC(RD) = 255; + RGBCmbnd2(:, :, 1) = RC; RGBCmbnd2(:, :, 2) = GC; RGBCmbnd2(:, :, 3) = BC; + subplot(1, 2, 2); imshow(RGBCmbnd2); + + if (uwntToPrint) + PrintPath = ['PrintPath\' useMyMdlFrm '\']; + print (fh, [PrintPath bkUpName '\' vidIndx '_' stckIndx], '-djpeg'); + imwrite (RGBCmbnd1, [PrintPath 'IndVid\' bkUpName '\' vidIndx '_' stckIndx '_Det.jpg']); + imwrite (RGBCmbnd2, [PrintPath 'IndVid\' bkUpName '\' vidIndx '_' stckIndx '_Trt.jpg']); + end + + %Lets count here + if (uWntToCnt) + Classfd = LblMsk; GndTruth = remFrmTruth; + [TP, FN, FP, TN] = countTPFNFPTN(segmentedImg, GndTruth, Classfd); + TotStat.Tp = TotStat.Tp+TP; TotStat.Fn = TotStat.Fn+FN; + TotStat.Fp = TotStat.Fp+FP; TotStat.Tn = TotStat.Tn+TN; + IndStat.name{K} = fnames(K).name; + IndStat.stat = [IndStat.stat; TP, FN, FP, TN]; +% IndStat + end + if (params.isSVMModel && ~params.loadFeat && (mod(K, 300) == 0)) + save('CNNPatches\Features\saveFeat', 'saveFeat'); + end + end +% alreadySegImg + if (~alreadySegmented) + save('CNNPatches\Segment\alreadySegImg', 'alreadySegImg'); + save('CNNPatches\Segment\alreadyWBCImg', 'alreadyWBCImg'); + end + if (params.isSVMModel && ~params.loadFeat) + save('CNNPatches\Features\saveFeat', 'saveFeat'); + end + if (uWntToCnt) + save (['myExp6\' useMyMdlFrm '\' bkUpName '_MalImdb\TotStat_' num2str(NetId)], 'TotStat'); + save (['myExp6\' useMyMdlFrm '\' bkUpName '_MalImdb\IndStat_' num2str(NetId)], 'IndStat'); + end + if (filtrOutDifclt) + DifSamSlctdByPgm.TrthDifcltStack = TrthDifcltStack; + DifSamSlctdByPgm.ImgsDifcltStack = ImgsDifcltStack; + DifSamSlctdByPgm.LblsDifcltStack = LblsDifcltStack; + WnglyClasfdExclDiff.WrngTrtExclDiffclt = WrngTrtExclDiffclt; + WnglyClasfdExclDiff.WrngImExclDiffclt = WrngImExclDiffclt; + WnglyClasfdExclDiff.LblsWrng = LblsWrng; + save ([bkUpName 'DifSamSlctdByPgm'], 'DifSamSlctdByPgm'); + save ([bkUpName 'WnglyClasfdExclDiff'], 'WnglyClasfdExclDiff'); + end +end +function [FeatSet, lstAccsd] = loadFeat(saveFeat, slideName, R, C, lstAccsd) + [~, totCnt] = size(saveFeat); + for i = lstAccsd+1:totCnt + if (strcmp(slideName, saveFeat(i).FeatName) && R == saveFeat(i).FeatR && C == saveFeat(i).FeatC) + FeatSet = saveFeat(i).FeatSet; + lstAccsd = i; + return; + end + end +end +function [TPstv, FNgtv, FPstv] = getCMat(ClsfdD, GnD, PtchSzM, imgUint8) +% [clearSetLTh, clearSetGTh, remSetGTh] = getBinaryImage(imUint8RGB); +% fnh = figure; + [SzM, SzN] = size(GnD); + HdeMask = false (SzM, SzN); + HdeMask(PtchSzM/2+1:SzM - PtchSzM/2, PtchSzM/2+1:SzN - PtchSzM/2) = true; + Trth = GnD & HdeMask; MClsfdD = ClsfdD & HdeMask; + %Check whether manual intervention is needed or not. + dilMsk = strel('disk', 17); + ManualTrth = imdilate(Trth, dilMsk); + NmManTr = bwconncomp(ManualTrth); NmManTrth = NmManTr.NumObjects; + ManualClfd = imdilate(MClsfdD, dilMsk); + NmManCd = bwconncomp(ManualClfd); NmManClfd = NmManCd.NumObjects; + + + TrthNumObj = round(sum(Trth(:))/81); + ClsfdNmObj = round(sum(MClsfdD(:))/81); + if ((NmManTrth ~= TrthNumObj) || (NmManClfd ~= ClsfdNmObj)) + TPstv = input('TPSTV = '); + FNgtv = input('FNgtv = '); + FPstv = input('FPstv = '); + else + [~, Mskd] = maskAllSharingObjects(MClsfdD, Trth); + TPstv = round(sum(Mskd(:))/81); %81 pixel for one marking + FNgtv = TrthNumObj - TPstv; + FPstv = ClsfdNmObj - TPstv; + end +end + +function [susReg, bgndMsk] = getMyRegionalMinima(imgUint8) + MnSz = 900; + hsv = rgb2hsv(imgUint8); + vlue = hsv(:, :, 3); + stDsk = strel('disk', 11); + openimg = imopen(vlue, stDsk); + mask = imregionalmin(openimg); + %First Filtering Exclude all the Bgnd + bgndMsk = bwareaopen((vlue < graythresh(vlue)), MnSz); + susReg = bgndMsk & mask; +end +function cDstLoc = getDstLocnsFor(vidIndx) + BasePath = 'DstLocnsPgm\'; + load ([BasePath 'DstLocByPGM']); cnt = 0; + vidIndxs = DstLocByPGM.vidIndx; + [~, numVidIndxs] = size(vidIndxs); + for i = 1:numVidIndxs + if (strcmp(vidIndx, vidIndxs(i).name)) + cDstLoc = DstLocByPGM.DstLocByPgm(:, :, i); + break; + end + end +end +function DstLoc = getDustLocns() +close all; + BasePath = ['E:\Gopakumar\GopakumarIISTDrive\Dataset\' ... + '4th IIST Visit\Mal01\Images\GndTruth\Separate\EasySeg\']; + fnames = dir([BasePath '\ImgesMinAcrsStck\*B.jpg']); + + numfids = length(fnames); M = 32; N = 32; R = 1; Rad = 0.75*M; + prevVidIndx = ''; DstCnt = 0; + for K = 388:numfids + K + AbsFNme = [BasePath 'ImgesMinAcrsStck\' fnames(K).name]; + %Get details needed to acces the Ground Truth file + [~, FileName, ~] = fileparts(AbsFNme); + for i = 1:length(FileName); if (FileName(i) == '_'); vid = i-1; break; end; end; + for j = vid+2:length(FileName); if (FileName(j) == '_'); stck = j-1; break; end; end; + vidIndx = FileName(1:vid); stckIndx = FileName(vid+2:stck); + if (strcmp(prevVidIndx, vidIndx)); continue; end + prevVidIndx = vidIndx; + figure; imshow(imread(AbsFNme)); title(num2str(K)); + ImgStack = getImageStackFromVideoIndx(vidIndx); + DstCnt = DstCnt+1; + DstLocByPgm = computeDustLocFrmStk(ImgStack); + save (['DstLocnsPgm\DstLocByPGM_' vidIndx], 'DstLocByPgm'); + if (strcmp(vidIndx, 'I') || strcmp(vidIndx, 'J')) + stp =1; + end + end +% getDstLocByMajority(); +end +function DstLocByPGM = getDstLocByMajority() + BasePath = ''; %['E:\Gopakumar\GopakumarIISTDrive\Dataset\' ... + %'E:\Gopakumar\GopakumarIISTDrive\Dataset\4th IIST Visit\Mal01\Images\GndTruth\Separate\EasySeg\']; + fnames = dir([BasePath 'DstLocnsPgm\DstLocByPGM_*.mat']); + + numfids = length(fnames); M = 32; N = 32; R = 1; Rad = 0.75*M; + VoteSecnd = uint8(zeros(480, 720)); + VoteFirst = uint8(zeros(480, 720)); + KTop = [25 27 29 31 33]; + close all; + for K = 1:numfids + AbsFNme = [BasePath 'DstLocnsPgm\' fnames(K).name]; + [~, FileName, ~] = fileparts(AbsFNme); + load (AbsFNme); + vidIndx = FileName(13:end); + [SzM, SzN] = size(DstLocByPgm); + if (K == 26 || K == 28 || K == 30 || K == 32 || K == 34) +% vidIndx + VoteSecnd(DstLocByPgm) = VoteSecnd(DstLocByPgm)+1; +% figure; imshow(VoteSecnd, []); + else + VoteFirst(DstLocByPgm) = VoteFirst(DstLocByPgm)+1; +% figure; imshow(VoteFirst, []); + end + display([num2str(K) ' ' vidIndx ' ' num2str(SzM) ' ' num2str(SzN)]); + end + DstLocByPgmFirst = VoteFirst > 15; + DstLocByPgmSecond = VoteSecnd > 3; +% DstLocByPgm = false(480, 720); + for K = 1:numfids + AbsFNme = [BasePath 'DstLocnsPgm\' fnames(K).name]; + [~, FileName, ~] = fileparts(AbsFNme); + vidIndx = FileName(13:end); + DstLocByPGM.vidIndx(K).name = vidIndx; + if (K == 26 || K == 28 || K == 30 || K == 32 || K == 34) + vidIndx + DstLocByPGM.DstLocByPgm(:, :, K) = DstLocByPgmSecond; + else + DstLocByPGM.DstLocByPgm(:, :, K) = DstLocByPgmFirst; + end + end + DstLocByPGM.meta = 'Genrtd By PGM testTrainedCNNNetOnSlide in MatCNetMalClass Fldr'; + save('DstLocnsPgm\DstLocByPGM', 'DstLocByPGM'); +end +function DstLocByPgm = computeDustLocFrmStk(ImgStackUint8) +close all; DstLocByPgm = 0; + [SzM, SzN, SzO, numCells] = size(ImgStackUint8); + numAssns = zeros(SzM, SzN); + Addns = zeros( SzM, SzN); + currImD = 0; Diffr = zeros( SzM, SzN); + objsIntst = false(SzM, SzN, numCells); + Vote = uint8(zeros(SzM, SzN)); + close all; + for i = 1:numCells + imUint8 = ImgStackUint8(:, :, :, i); + [clearSetLTh, clearSetGTh, remSetGTh] = getBinaryImage(imUint8); + segImg = clearSetLTh | clearSetGTh | remSetGTh; + bgndImg = ~segImg; +% figure; imshow(bgndImg); + bgndErde = bwmorph(bgndImg, 'erode', 25); + G = imUint8(:, :, 2); + avgbgnd = mean(G(bgndErde)); + objsCand = G < (avgbgnd - 0.05*avgbgnd); + %remove all bigger ones + intstd = objsCand & ~bwareaopen(objsCand, 200); + objsIntst(:, :, i) = intstd; + Vote(intstd) = Vote(intstd)+1; +% figure; imshow(imUint8); title(num2str(i)); +% figure; imshow(intstd); + end + %Select those having 25% Support + mnSprt = 0.25; %support + DstLocByPgm = Vote > max(round(mnSprt*numCells), 4); +% figure; imshow(imUint8); +% figure; imshow(fnalDstLocs); +end +function [clearSetLTh, clearSetGTh, remSetGTh] = getBinaryImage(imUint8RGB) + imG = im2double(imUint8RGB(:, :, 2)); + msk = fspecial('average', 15); + avG = imfilter(imG, msk); + fildLThImg = bwfill(bwareaopen(imG < (avG - 0.01), 200), 'holes'); + fildGThImg = lOtThresh(imUint8RGB); + [clearSetLTh, ~] = getClearSet(fildLThImg); + remSetGTh = bwareaopen(bwmorph(~clearSetLTh & fildGThImg, 'open', 3), 200); + [clearSetGTh, remSetGTh] = getClearSet(remSetGTh); +end + +function [clearSet, remSet] = getClearSet(bwImg) + remSet = bwImg; clearSet = false(size(remSet)); + cmps = regionprops(bwImg, {'PixelIdxList', 'Solidity', 'ConvexImage', 'BoundingBox'}); + [numObj, ~] = size(cmps); + lowThresh = 750; highThresh = 2000; + for i = 1:numObj + currObj = cmps(i).PixelIdxList; + objArea = length(currObj); + if (cmps(i).Solidity > 0.9 && objArea > lowThresh && objArea < highThresh) + clearSet(currObj) = true; + remSet(currObj) = false; + end + end +end + +function ThImg = lOtThresh(im) + div = 2; + [SzM, SzN, SzO] = size(im); + if (SzO > 3) + im = rgb2gray(im); %im(:, :, 2); + end + im = im2double(im); + rOfst = round(SzM/div); cOfst = round(SzN/div); + ThImg = false(SzM, SzN); + for i = 1:div + startR = (i-1)*rOfst+1; + if (i == div) + endR = SzM; + else + endR = i*rOfst; + end + for j = 1:div + startC = (j-1)*cOfst+1; + if (j == div) + endC = SzN; + else + endC = j*cOfst; + end + divImg = im(startR:endR, startC:endC); + ThImg(startR:endR, startC:endC) = divImg < (graythresh(divImg)+0.01); + end + end + ThImg = bwfill(bwareaopen(ThImg, 200), 'holes'); +end + +function [remMsk, Mskd] = maskAllSharingObjects(BaseImg, ShareObjs) + ToBMaskd = BaseImg & ShareObjs; + LblsFrmBase = bwlabel(BaseImg); + Mskd = false(size(BaseImg)); + lbls2BMaskd = LblsFrmBase(ToBMaskd); + unqLbls = unique(lbls2BMaskd(:)); + for lbl = 1:length(unqLbls) + Mskd(LblsFrmBase == (unqLbls(lbl))) = true; + end + remMsk = BaseImg & ~Mskd; +end + +function ImgStack = getImageStackFromVideoIndx(vidIndx) + BasePath = ['E:\Gopakumar\GopakumarIISTDrive\Dataset\' ... + '4th IIST Visit\Mal01\Images\GndTruth\Separate\EasySeg\']; + fnames = dir([BasePath '\ImgesMinAcrsStck\' vidIndx '*B.jpg']); + + numfids = length(fnames); M = 32; N = 32; R = 1; Rad = 0.75*M; + for K = 1:numfids + AbsFNme = [BasePath 'ImgesMinAcrsStck\' fnames(K).name]; + im = imread(AbsFNme); + ImgStack(:, :, :, K) = im; + end +end + +function [FnlFeatSet, CellId] = getMyFeatures(cellsImg, imUint8RGB, slideName) + %Lets do the processing on Green channel for the time being + avMsk = fspecial('average', 3); + FnlFeatSet = single(zeros(1, 14*3)); %3 Channels + for ch = 1:3 + imGreen = imUint8RGB(:, :, ch); + imGreenDble = im2double(imGreen); + [grdMg, ~] = imgradient(imGreenDble); + meanLclMsk = imfilter(imGreenDble, avMsk); + varnLclMsk = imfilter(imGreenDble.^2, avMsk) - meanLclMsk.^2; + [SzM, SzN] = size(imGreen); + props = regionprops(cellsImg, {'PixelIdxList', 'Centroid'}); + [numObjs, ~] = size(props); + myBgndImgI = uint8(200* ones(SzM, SzN)); + FeatSet = single(zeros(numObjs, 14)); + CellId = []; + for i = 1:numObjs + pxlIdxLst = props(i).PixelIdxList; + [SR, SC] = ind2sub(size(imGreen), pxlIdxLst); + tmpI = myBgndImgI; + tmpI(pxlIdxLst) = imGreen(pxlIdxLst); + %Get the patch + minSR = min(SR); maxSR = max(SR); + minSC = min(SC); maxSC = max(SC); + imgPatch = tmpI(minSR:maxSR, minSC:maxSC); + mskPatch = cellsImg(minSR:maxSR, minSC:maxSC); + FeatSet(i, 1:4) = getGLCMFeatPrPatch(imgPatch, mskPatch); + + imGrnDblPchPxls = imGreenDble(pxlIdxLst); + meanGlblPxls = mean(imGrnDblPchPxls); + varGlblPxls = var(imGrnDblPchPxls); + minPxls = min(imGrnDblPchPxls); + maxPxls = max(imGrnDblPchPxls); + minGrdMag = min(grdMg(pxlIdxLst)); + maxGrdMag = max(grdMg(pxlIdxLst)); + minLclMean = min(meanLclMsk(pxlIdxLst)); + maxLclMean = max(meanLclMsk(pxlIdxLst)); + minLclVarn = min(varnLclMsk(pxlIdxLst)); + maxLclVarn = max(varnLclMsk(pxlIdxLst)); + FeatSet(i,5:14) = [meanGlblPxls varGlblPxls minPxls maxPxls ... + minGrdMag maxGrdMag minLclMean maxLclMean ... + minLclVarn maxLclVarn]; + cellCentr = round(props(i).Centroid); + CellId(i).name = [slideName '_' num2str(cellCentr(2)) '_' num2str(cellCentr(1))]; + end + FnlFeatSet(1, (ch-1)*14+1:ch*14) = FeatSet; + end +end + +function GLCMFeat = getGLCMFeatPrPatch(imUint8Gray, Msk) + im = imUint8Gray; + %im has to be Single band [0 255] + [SzR, SzC, SzO] = size(im); + MskX = [Msk(:, 2:end) Msk(:, end)]; + MskY = [Msk(2:end, :); Msk(end, :)]; + MskD = Msk & MskX & MskY; %MskD = bwmorph(MskD, 'erode', 3); + %Compute the GLCM for each band for the reg and return the props. + Lvl = 32; LvlDiv = 256/Lvl; + regIntst = MskD; + GLCMMat = zeros(Lvl, Lvl); + GLCMFeat = zeros(SzO, 4); + for band = 1:SzO + imBand = double(im); + imgL = ceil((imBand+1)/LvlDiv); + imgL(~regIntst) = -1; + for i = 1:SzR + for j = 1:SzC-1 + valLeft = imgL(i, j); + valRight = imgL(i, j+1); + if (valLeft ~= -1 && valRight ~= -1) + GLCMMat(valLeft, valRight) = GLCMMat(valLeft, valRight)+1; + end + end + end + GFeat = graycoprops(GLCMMat); + GLCMFeat(band, :) = [GFeat.Contrast GFeat.Correlation GFeat.Energy GFeat.Homogeneity]; + end +end \ No newline at end of file