smiuser/cytocube
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

malaria detection code in matlab

Browse Source
This commit is contained in:
moonanjum26
2018-08-22 22:10:26 +05:30
parent b178a5db75
commit ad7e0b7a0c
15 changed files with 2846 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
src/gopa kumar code/CNNTRainFczRGBCall.m Normal file
View File

@@ -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

42
src/gopa kumar code/CnnSvmOnSlideCalling.m Normal file
View File

@@ -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

431
src/gopa kumar code/MalEx6ColorRewrittenTrainFcZRGB.m Normal file
View File

@@ -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');

25
src/gopa kumar code/Read Me.txt Normal file
View File

@@ -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'

48
src/gopa kumar code/SVMTrainBFcsRGB.m Normal file
View File

@@ -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

74
src/gopa kumar code/countTPFNFPTN.m Normal file
View File

@@ -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

289
src/gopa kumar code/generate32X32FcsStck_Fcsd_Feat_RGB.m Normal file
View File

@@ -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

282
src/gopa kumar code/generateDstAndWBCPatches.m Normal file
View File

@@ -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

649
src/gopa kumar code/getSegmentation.m Normal file
View File

@@ -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

37
src/gopa kumar code/initializeCharacterCNNBNNEx6ColorFczRGB.m Normal file
View File

@@ -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') ;

310
src/gopa kumar code/makeDataSetForTrainingBest.m Normal file
View File

@@ -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

30
src/gopa kumar code/reName.m Normal file
View File

@@ -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

25
src/gopa kumar code/setup.m Normal file
View File

@@ -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

13
src/gopa kumar code/testCNNBNNEx6Color.m Normal file
View File

@@ -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

582
src/gopa kumar code/testTrainedCNNSVMAutoCnt.m Normal file
View File

@@ -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