Homec4science
Diffusion segmentation (master) Browse

segmentation/autoSegment.m
8222343c77d4master

Last Change
Raw File

autoSegment.m
View Options

function finalSurface = autoSegment(Img, imgInfo, model, opt)
% Img a 3D image
% imgInfo a struct storing the information of the image:
% imgInfo.PixelSpacing(1)
% imgInfo.imgInfo.PixelSpacing(2)
% imgInfo.imgInfo.SliceThickness
% imgInfo.BoundingBox.xMin
% imgInfo.BoundingBox.yMin
% imgInfo.BoundingBox.zMin
% imgInfo.BoundingBox.xMax
% imgInfo.BoundingBox.yMax
% imgInfo.BoundingBox.zMax
% model Information regarding the statistical shape model:
% model.file % path to the ssmFile;
% model.bParams_train % scores of training samples
% model.basisVect % basis vectors of SSM
% model.meanShape % meanShape
% model.pcaSTD % eigen values of SSM
% model.cells % elements of shapes in the SSM
% model.confidence % confidence at of SSM fitting error for each node
% model.gVertices % nodes related to the glenoid joint
% model.trainBones % training samples;
% model.samplesPath % path to the folder containing samples
% model.LM_Extra_Fix % extra landmorks at the extremes of the scapula
% model.descSampleImg % descriptors of extracted sift features on the model
% model.indicesSampleImg % index of sift features extracted on the labels of model's mean shape
% Example:
% % segmentation parameters are set here:
% opt.isLeft = false;
% opt.doSmoothing = true;
% opt.ParallelProcess = true;
% opt.numTemplates = 5;
% opt.threshold = 190;
% opt.volThr = 35000;
% opt.ransacMaxIter = 5.0e6;
% opt.ransacDist = 6; %mm
% opt.ransacNo = 4; % number of samples for ransac!
% model = loadModel('../Model/', false); % reads a model for right
% [imgInfo, Img] = LoadData_Amira('path to amira Image');
% segmentationResult = autoSegment(Img, imgInfo, model, opt);
%
%tic
if(~isfield(imgInfo, 'xSize'))
imgInfo.xSize = size(Img, 2);
end
if(~isfield(imgInfo, 'ySize'))
imgInfo.ySize = size(Img, 1);
end
if(~isfield(imgInfo, 'xzSize'))
imgInfo.zSize = size(Img, 3);
end
if(~isfield(opt, 'isLeft'))
opt.isLeft = false;
end
if(~isfield(opt, 'doSmoothing'))
opt.doSmoothing = false;
end
if(~isfield(opt, 'ParallelProcess'))
opt.ParallelProcess = true;
end
if(~isfield(opt, 'numTemplates'))
opt.numTemplates = 3;
end
if(~isfield(opt, 'threshold'))
opt.threshold = 190;
end
if(~isfield(opt, 'volThr'))
opt.volThr = 25000;
end
if(~isfield(opt, 'ransacMaxIter'))
opt.ransacMaxIter = 8.0e6;
end
if(~isfield(opt, 'ransacDist'))
opt.ransacDist = 6;
end
if(~isfield(opt, 'ransacNo'))
opt.ransacNo = 4;
end
opt.searchWindowSize = [65, 65, 45];
Parallel = opt.ParallelProcess;
if Parallel
try
p = gcp('nocreate'); % If no pool, do not create new one.
if isempty(p)
parpool();
end
catch
Parallel = false;
end
end
if(opt.doSmoothing)
eval(['Img = ' class(Img) '(smooth3(Img, ''gaussian'', [3, 3, 3], 0.8));']);
end
if opt.isLeft % if we have to do calculations for a mirrored model! sift feature is sensitive to the flipping
T_initial = [-1 0 0; 0 1 0; 0 0 1]; % transformation for reflection
else % if not!
T_initial = eye(3);
end
Img(Img<-2000) = -2000;
Img(Img>1000) = 1000;
imgLabel = false(size(Img));
imgLabel(Img>=opt.threshold) = 1;
voxVol = imgInfo.PixelSpacing(1) * imgInfo.PixelSpacing(2) * imgInfo.SliceThickness;
imgLabel = bwareaopen(imgLabel, ceil(opt.volThr/voxVol));
if Parallel
parfor ii = 1:size(imgLabel, 2)
imgLabel(:, ii, :) = imfill(squeeze(imgLabel(:, ii, :)), 8, 'holes');
end
parfor ii = 1:size(imgLabel, 1)
imgLabel(ii, :, :) = imfill(squeeze(imgLabel(ii, :, :)), 8, 'holes');
end
else
for ii = 1:size(imgLabel, 2)
imgLabel(:, ii, :) = imfill(squeeze(imgLabel(:, ii, :)), 8, 'holes');
end
for ii = 1:size(imgLabel, 2)
imgLabel(ii, :, :) = imfill(squeeze(imgLabel(ii, :, :)), 8, 'holes');
end
end
finalSurface = [];
x = imgInfo.BoundingBox.xMin : imgInfo.PixelSpacing(1) : imgInfo.BoundingBox.xMax;
y = imgInfo.BoundingBox.yMin : imgInfo.PixelSpacing(2) : imgInfo.BoundingBox.yMax;
z = imgInfo.BoundingBox.zMin : imgInfo.SliceThickness : imgInfo.BoundingBox.zMax;
[x, y, z] = meshgrid(x, y, z);
[fv.faces, fv.vertices] = isosurface(x, y, z, imgLabel, 0);
[fv.faces, fv.vertices] = reducepatch(fv, 0.05);
InitialSegmentation = fv;
descSampleImg = model.descSampleImg;
indicesSampleImg = model.indicesSampleImg;
keys = detectSift3D(imgLabel(1:2:end, 1:2:end, 1:2:end), 'units', 2*[imgInfo.PixelSpacing(1) imgInfo.PixelSpacing(2) imgInfo.SliceThickness], 'peakThresh', 0.25, 'cornerThresh', 0.15, 'numKpLevels', 3, 'sigmaN', 1.15, 'sigma0', 1.6);
[descImg, indicesImg] = extractSift3D(keys);%, imgLabel);
clear mex
% check MATLAB license:
%C:\Program Files\MATLAB\R2016a\etc\win64>lmutil.exe lmstat -a -c 27000@id-lizenzserver-matlab.unibe.ch
MaxRatio = 0.8;
MatchThreshold = 80;
n = 1;
while(1)
try
[indexPairs, ~] = matchFeatures(descSampleImg, descImg, 'MatchThreshold', MatchThreshold, 'MaxRatio', MaxRatio, 'Unique', true, 'Method', 'Exhaustive', 'Metric', 'SSD'); % 'Approximate');
break;
catch
n = n + 1;
% disp(['autoSegment/matchFeatures/n' num2str(n)]);
if(mod(n, 1.e5) == 0)
% display('I am waiting for Video_and_Image_Block_Set Toolbox license!'); % Commented by AT
% disp('lost in a while loop, line 163 of autoSegment');
end
% 4 lines below added by AT
if n > 2.e5
fprintf('Error in matchFeatures in autoSegment');
return
end
end
end
coordsSampleImg = bsxfun(@plus, bsxfun(@times, indicesSampleImg, ...
[model.sampleImageInfo.PixelSpacing(2), model.sampleImageInfo.PixelSpacing(1), model.sampleImageInfo.SliceThickness]), ...
[model.sampleImageInfo.BoundingBox.yMin, model.sampleImageInfo.BoundingBox.xMin, model.sampleImageInfo.BoundingBox.zMin]);
coordsImg = bsxfun(@plus, bsxfun(@times, indicesImg, ...
2*[imgInfo.PixelSpacing(2), imgInfo.PixelSpacing(1), imgInfo.SliceThickness]), ...
[imgInfo.BoundingBox.yMin, imgInfo.BoundingBox.xMin, imgInfo.BoundingBox.zMin]);
movLM = [coordsImg(indexPairs(:, 2), 2), coordsImg(indexPairs(:, 2), 1), coordsImg(indexPairs(:, 2), 3)];
fixLM = [coordsSampleImg(indexPairs(:, 1), 2), coordsSampleImg(indexPairs(:, 1), 1), coordsSampleImg(indexPairs(:, 1), 3)];
indicesSampleImg = indicesSampleImg(indexPairs(:, 1), :);
numMatches = size(fixLM, 1);
options.thr = opt.ransacDist;
options.d = opt.ransacNo;
options.reflection = false;
if Parallel
possibleNumPermutations = nchoosek(numMatches, options.d);
while possibleNumPermutations>opt.ransacMaxIter
MaxRatio = MaxRatio - 0.05;
MatchThreshold = MatchThreshold - 5;
descSampleImg = descSampleImg(indexPairs(:, 1), :);
descImg = descImg(indexPairs(:, 2), :);
[indexPairs, ~] = matchFeatures(descSampleImg, descImg, 'MatchThreshold', MatchThreshold, 'MaxRatio', MaxRatio, 'Unique', true, 'Method', 'Exhaustive', 'Metric', 'SSD'); % 'Approximate');
movLM = movLM(indexPairs(:, 2), :);
fixLM = fixLM(indexPairs(:, 1), :);
indicesSampleImg = indicesSampleImg(indexPairs(:, 1), :);
numMatches = size(fixLM, 1);
possibleNumPermutations = nchoosek(numMatches, options.d);
end
options.numIterations = possibleNumPermutations;
options.type = 'exhaustive';
else
options.numIterations = 5.e5;
options.type = 'random';
end
options.parallel = Parallel;
if (size(fixLM, 1)<options.d)
message(['The shape looks very different from mean shape!' 10 ...
'Could not find enough number of landmarks for segmentation!' 10 ...
'Please check the threshold value and the initial segmentation for this bone!']);
return;
else
[~, inliers] = estimateGeometricTransform3D(fixLM, movLM, options);
end
fixLM = fixLM(inliers, :);
movLM = movLM(inliers, :);
if (size(fixLM, 1)<3)
message(['The shape looks very different from mean shape!' 10 ...
'Could not find enough number of landmarks for segmentation!' 10 ...
'Please check the threshold value and the initial segmentation for this bone!']);
return;
end
indicesSampleImg = indicesSampleImg(inliers, :)+1;
onEdgeInd = AvoidCrowdedLandmarks(fixLM, indicesSampleImg, model.sampleImage, 3, Parallel);
fixLM = fixLM * T_initial;
[~, ~, rigidTransform] = procrustes(fixLM, movLM, 'scaling', false, 'reflection', opt.isLeft);
tform = affine3d([[rigidTransform.T; rigidTransform.c(1, :)] [zeros(3, 1); 1]]);
RA = imref3d(size(Img), [imgInfo.BoundingBox.xMin imgInfo.BoundingBox.xMax], ...
[imgInfo.BoundingBox.yMin imgInfo.BoundingBox.yMax], ...
[imgInfo.BoundingBox.zMin imgInfo.BoundingBox.zMax]);
try
[ImgRigid, RB] = imwarp(Img, RA, tform, 'linear', 'FillValues', -2000);
catch
display (['Could not transform the image, most probably due to the size (' num2str(size(Img)) ')']);
return;
end
imgInfoRigid.BoundingBox.xMin = RB.XWorldLimits(1);
imgInfoRigid.BoundingBox.yMin = RB.YWorldLimits(1);
imgInfoRigid.BoundingBox.zMin = RB.ZWorldLimits(1);
imgInfoRigid.BoundingBox.xMax = RB.XWorldLimits(2);
imgInfoRigid.BoundingBox.yMax = RB.YWorldLimits(2);
imgInfoRigid.BoundingBox.zMax = RB.ZWorldLimits(2);
imgInfoRigid.PixelSpacing(1) = (RB.XWorldLimits(2) - RB.XWorldLimits(1))/(size(ImgRigid, 2)-1);
imgInfoRigid.PixelSpacing(2) = (RB.YWorldLimits(2) - RB.YWorldLimits(1))/(size(ImgRigid, 1)-1);
imgInfoRigid.SliceThickness = (RB.ZWorldLimits(2) - RB.ZWorldLimits(1))/(size(ImgRigid, 3)-1);
% resample image if necessary:
xyzResample = [imgInfoRigid.PixelSpacing(1), imgInfoRigid.PixelSpacing(2), imgInfoRigid.SliceThickness] ./ [model.sampleImageInfo.PixelSpacing(1), model.sampleImageInfo.PixelSpacing(2), model.sampleImageInfo.SliceThickness];
if sum(xyzResample < 0.90 * [1, 1, 1] | xyzResample > 1.1 * [1, 1, 1]) ~= 0
newSize = round(size(ImgRigid) .* xyzResample);
if prod(newSize)>5.e9
display (['Skipped resampling the image, due to the size of image (' num2str(size(ImgRigid)) '), crop the image and then try again!']);
return;
end
try
ImgRigid = int16(resize(double(ImgRigid), newSize, 'linear'));
catch
display ('Could not resample the image, most probably due to the size');
return;
end
imgInfoRigid.PixelSpacing(1) = model.sampleImageInfo.PixelSpacing(1);
imgInfoRigid.PixelSpacing(2) = model.sampleImageInfo.PixelSpacing(2);
imgInfoRigid.SliceThickness = model.sampleImageInfo.SliceThickness;
end
imgInfoRigid.xSize = size(ImgRigid, 2);
imgInfoRigid.ySize = size(ImgRigid, 1);
imgInfoRigid.zSize = size(ImgRigid, 3);
%%%%%%%%%%%%%%%%%%%%% find 3 essential landmarks %%%%%%%%%%%%%%%%%%%%%%%%
tmpLM_Move = round(bsxfun(@times, bsxfun(@minus, model.LM_Extra_Fix, [imgInfoRigid.BoundingBox.xMin, imgInfoRigid.BoundingBox.yMin, imgInfoRigid.BoundingBox.zMin]), 1./[imgInfoRigid.PixelSpacing(2) imgInfoRigid.PixelSpacing(1) imgInfoRigid.SliceThickness]));
tmpLM_Move = tmpLM_Move(:, [2, 1, 3]);
LM_Extra_Mov = zeros(size(model.LM_Extra_Fix));
skipInd = [];
sw = false;
LMWindowSz = size(model.AtlasPatch_LM);
LMWindowSz = LMWindowSz(2:end);
for ii = 1:size(model.LM_Extra_Fix, 1)
searchRangeX = [max(tmpLM_Move(ii, 1) - 60, floor(LMWindowSz(2)/2)), min(tmpLM_Move(ii, 1) + 60, size(ImgRigid, 1) - floor(LMWindowSz(2)/2))];
searchRangeY = [max(tmpLM_Move(ii, 2) - 60, floor(LMWindowSz(1)/2)), min(tmpLM_Move(ii, 2) + 60, size(ImgRigid, 2) - floor(LMWindowSz(1)/2))];
searchRangeZ = [max(tmpLM_Move(ii, 3) - 60, floor(LMWindowSz(3)/2)), min(tmpLM_Move(ii, 3) + 60, size(ImgRigid, 3) - floor(LMWindowSz(3)/2))];
if(searchRangeX(2)<=searchRangeX(1) || searchRangeY(2)<=searchRangeY(1) || searchRangeZ(2)<=searchRangeZ(1))
LM_Extra_Mov(ii, :) = tmpLM_Move(ii, :);
sw = true;
continue;
end
Patch = ImgRigid(searchRangeX(1):searchRangeX(2), searchRangeY(1):searchRangeY(2), searchRangeZ(1):searchRangeZ(2));
if(size(Patch, 1) <= size(model.AtlasPatch_LM, 2) || size(Patch, 2) <= size(model.AtlasPatch_LM, 3) || size(Patch, 3) <= size(model.AtlasPatch_LM, 4))
LM_Extra_Mov(ii, :) = tmpLM_Move(ii, :);
sw = true;
continue;
end
CC = normxcorr3(double(squeeze(model.AtlasPatch_LM(ii, :, :, :))), double(Patch), 'same');
[v, ind] = max(CC(:));
if(v<0.5)
LM_Extra_Mov(ii, :) = tmpLM_Move(ii, :);
sw = true;
continue;
end
if(length(ind)~=1)
h = fspecial3('gaussian', size(CC), 2);
[~, ind] = max(CC(:).*h(:));
ind = ind(1);
end
[sub(1), sub(2), sub(3)] = ind2sub(size(CC), ind);
sub = sub + [searchRangeX(1), searchRangeY(1), searchRangeZ(1)];
LM_Extra_Mov (ii, :) = sub;
end
LM_Extra_Mov = bsxfun(@plus, bsxfun(@times, LM_Extra_Mov, ...
[imgInfoRigid.PixelSpacing(2), imgInfoRigid.PixelSpacing(1), imgInfoRigid.SliceThickness]), ...
[imgInfoRigid.BoundingBox.yMin, imgInfoRigid.BoundingBox.xMin, imgInfoRigid.BoundingBox.zMin]);
LM_Extra_Mov = LM_Extra_Mov(:, [2, 1, 3]);
LM_Extra_Mov = bsxfun(@minus, LM_Extra_Mov, rigidTransform.c(1, :))*rigidTransform.T';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% The extra lanmarks were not used for aligning the image so I won't
% include it for aligning initial segmentation either!
ind = 1:size(model.LM_Extra_Fix, 1);
ind(skipInd) = [];
fixLM = [fixLM; model.LM_Extra_Fix(ind, :)];
movLM = [movLM; LM_Extra_Mov(ind, :)];
movLM = bsxfun(@plus, movLM * rigidTransform.T, rigidTransform.c(1, :));
fv = InitialSegmentation;
fv.vertices = bsxfun(@plus, fv.vertices * rigidTransform.T, rigidTransform.c(1, :));
options.mirror = opt.isLeft;
options.initialColor = 'none';
options.fitColor = 'none';
options.replace = false;
options.InitialTransform.T = eye(3);%T_initial;
options.InitialTransform.c = [0 0 0];
options.Extras = numel(ind);
options.fitLM = [onEdgeInd; true(options.Extras, 1)];
options.Name = '';
if (sw)
options.fitRegularization = 1;
confidence = model.confidence * 2.0;
else
options.fitRegularization = 0.0005;
confidence = model.confidence;
end
[fvFit, bParams] = fitModel(model.file, fv, fixLM, movLM, options);
similarity = zeros(1, size(model.bParams_train, 1));
weight = model.pcaSTD./sum(model.pcaSTD);
for ii = 1:size(model.bParams_train, 1)
similarity(ii) = sum(abs(bParams' - model.bParams_train(ii, :)) .* weight);
end
[~, idx] = sort(similarity, 'ascend');
numVertices = size(model.meanShape, 1);
gVertices = model.gVertices;
Par = load([model.samplesPath 'Patches_' model.trainBones{idx(1)}]);
AtlasPatches = zeros(numVertices, size(Par.Patches, 2), size(Par.Patches, 3), size(Par.Patches, 4), opt.numTemplates);
AtlasGPatches = zeros(length(gVertices)-1, size(Par.GPatches, 2), size(Par.GPatches, 3), size(Par.GPatches, 4), opt.numTemplates); %%
AtlasPatches(:, :, :, :, 1) = Par.Patches;
AtlasGPatches(:, :, :, :, 1) = Par.GPatches;
for ii = 2:opt.numTemplates
Par = load([model.samplesPath 'Patches_' model.trainBones{idx(ii)}]);%squeeze(Patches(closestSample, :, :, :, :));
AtlasPatches(:, :, :, :, ii) = Par.Patches;
AtlasGPatches(:, :, :, :, ii) = Par.GPatches;
end
clear Par;
eps = 0.00001;
patchSz = [size(AtlasPatches, 2), size(AtlasPatches, 3), size(AtlasPatches, 4)];
GpatchSz = [size(AtlasGPatches, 2), size(AtlasGPatches, 3), size(AtlasGPatches, 4)]; %%
imgInfo = imgInfoRigid;
Img = ImgRigid;
voxCenter = [imgInfo.PixelSpacing(1)/2, imgInfo.PixelSpacing(2)/2, imgInfo.SliceThickness/2] + eps;
% zero padding
tmp = zeros(size(Img)+opt.searchWindowSize) - 2000;
tmp ((opt.searchWindowSize(1)-1)/2+1:end-(opt.searchWindowSize(1)-1)/2-1, ...
( opt.searchWindowSize(2)-1)/2+1:end-(opt.searchWindowSize(2)-1)/2-1, ...
( opt.searchWindowSize(3)-1)/2+1:end-(opt.searchWindowSize(3)-1)/2-1) = Img;
Img = tmp;
tmp = [];
imgInfo.BoundingBox.xMin = imgInfo.BoundingBox.xMin - (opt.searchWindowSize(1)-1)/2 * imgInfo.PixelSpacing(1);
imgInfo.BoundingBox.yMin = imgInfo.BoundingBox.yMin - (opt.searchWindowSize(2)-1)/2 * imgInfo.PixelSpacing(2);
imgInfo.BoundingBox.zMin = imgInfo.BoundingBox.zMin - (opt.searchWindowSize(3)-1)/2 * imgInfo.SliceThickness;
imgInfo.xSize = imgInfo.xSize + opt.searchWindowSize(1);
imgInfo.ySize = imgInfo.ySize + opt.searchWindowSize(2);
imgInfo.zSize = imgInfo.zSize + opt.searchWindowSize(3);
% we will use these vectors to quickly find the closest voxel to the
% vertices
xImg = imgInfo.BoundingBox.xMin + (0:imgInfo.xSize-1) * imgInfo.PixelSpacing(1);
yImg = imgInfo.BoundingBox.yMin + (0:imgInfo.ySize-1) * imgInfo.PixelSpacing(2);
zImg = imgInfo.BoundingBox.zMin + (0:imgInfo.zSize-1) * imgInfo.SliceThickness;
disp = zeros(numVertices, 3);
g = 1; %%
if Parallel %&& (ispc() || isunix())
P{numVertices} = zeros(opt.searchWindowSize);
Atlas{numVertices} = 0;
for v = 1:numVertices
indX = find(abs(xImg - fvFit.vertices(v, 1)) <= voxCenter(1));
indY = find(abs(yImg - fvFit.vertices(v, 2)) <= voxCenter(2));
indZ = find(abs(zImg - fvFit.vertices(v, 3)) <= voxCenter(3));
if (isempty(indX) || isempty(indY) || isempty(indZ)) % partial image!
P{v} = [];
Atlas{v} = [];
continue;
end
indX = indX(1); indY = indY(1); indZ = indZ(1);
if(gVertices(g) == v) %%
searchWindowSize = GpatchSz + 2*ceil(0.5 + confidence(v)./[imgInfo.PixelSpacing(1) imgInfo.PixelSpacing(2) imgInfo.SliceThickness]);
if(indX < (searchWindowSize(1)-1)/2 || indY < (searchWindowSize(2)-1)/2 || indZ < (searchWindowSize(3)-1)/2 || ...
indX > size(Img, 2) - (searchWindowSize(1)-1)/2 || indY > size(Img, 1) - (searchWindowSize(2)-1)/2 || indZ > size(Img, 3) - (searchWindowSize(3)-1)/2)
P{v} = [];
Atlas{v} = [];
continue;
end
P{v} = Img(indY-(searchWindowSize(2)-1)/2+1 : indY+(searchWindowSize(2)-1)/2+1, ... %%
indX-(searchWindowSize(1)-1)/2+1 : indX+(searchWindowSize(1)-1)/2+1, ... %%
indZ-(searchWindowSize(3)-1)/2+1 : indZ+(searchWindowSize(3)-1)/2+1); %%
Atlas{v} = squeeze(AtlasGPatches(g, :, :, :, :));
g = g+1; %%
else %%
searchWindowSize = patchSz + 2*ceil(0.5 + confidence(v)./[imgInfo.PixelSpacing(1) imgInfo.PixelSpacing(2) imgInfo.SliceThickness]);
if(indX < (searchWindowSize(1)-1)/2 || indY < (searchWindowSize(2)-1)/2 || indZ < (searchWindowSize(3)-1)/2 || ...
indX >= size(Img, 2) - (searchWindowSize(1)-1)/2 || indY >= size(Img, 1) - (searchWindowSize(2)-1)/2 || indZ >= size(Img, 3) - (searchWindowSize(3)-1)/2)
P{v} = [];
Atlas{v} = [];
continue;
end
P{v} = Img(indY-(searchWindowSize(2)-1)/2+1 : indY+(searchWindowSize(2)-1)/2+1, ...
indX-(searchWindowSize(1)-1)/2+1 : indX+(searchWindowSize(1)-1)/2+1, ...
indZ-(searchWindowSize(3)-1)/2+1 : indZ+(searchWindowSize(3)-1)/2+1);
Atlas{v} = squeeze(AtlasPatches(v, :, :, :, :));
end %%
disp(v, :) = [xImg(indX), yImg(indY), zImg(indZ)] - fvFit.vertices(v, :);
end
disp = disp + correction(P, Atlas, imgInfo, Parallel);
clear Atlas;
clear P;
clear Img;
clear ImgRigid;
else
sub = [0, 0, 0];
% handle = waitbar(0, 'Please wait ...', 'Name','Correction in progress...');
for v = 1:numVertices
indX = find(abs(xImg - fvFit.vertices(v, 1)) <= voxCenter(1));
indY = find(abs(yImg - fvFit.vertices(v, 2)) <= voxCenter(2));
indZ = find(abs(zImg - fvFit.vertices(v, 3)) <= voxCenter(3));
if (isempty(indX) || isempty(indY) || isempty(indZ)) % partial image!
disp(v, :) = [0, 0, 0];
continue;
end
indX = indX(1); indY = indY(1); indZ = indZ(1);
if(gVertices(g) == v) %%
searchWindowSize = GpatchSz + 2*ceil(0.5 + confidence(v)./[imgInfo.PixelSpacing(1) imgInfo.PixelSpacing(2) imgInfo.SliceThickness]);
if(indX < (searchWindowSize(1)-1)/2 || indY < (searchWindowSize(2)-1)/2 || indZ < (searchWindowSize(3)-1)/2 || ...
indX >= size(Img, 2) - (searchWindowSize(1)-1)/2 || indY >= size(Img, 1) - (searchWindowSize(2)-1)/2 || indZ >= size(Img, 3) - (searchWindowSize(3)-1)/2)
continue;
end
P = Img(indY-(searchWindowSize(2)-1)/2+1 : indY+(searchWindowSize(2)-1)/2+1, ... %%
indX-(searchWindowSize(1)-1)/2+1 : indX+(searchWindowSize(1)-1)/2+1, ... %%
indZ-(searchWindowSize(3)-1)/2+1 : indZ+(searchWindowSize(3)-1)/2+1); %%
ncc3 = normxcorr3_multi(squeeze(AtlasGPatches(g, :, :, :, :)), P, 'valid');
ncc3(ncc3<0.5) = 0;
ncc = squeeze(ncc3(:, :, :, 1) + ncc3(:, :, :, 2) + ncc3(:, :, :, 3));
g = g+1; %%
else %%
searchWindowSize = patchSz + 2*ceil(0.5 + confidence(v)./[imgInfo.PixelSpacing(1) imgInfo.PixelSpacing(2) imgInfo.SliceThickness]);
if(indX < (searchWindowSize(1)-1)/2 || indY < (searchWindowSize(2)-1)/2 || indZ < (searchWindowSize(3)-1)/2 || ...
indX >= size(Img, 2) - (searchWindowSize(1)-1)/2 || indY >= size(Img, 1) - (searchWindowSize(2)-1)/2 || indZ >= size(Img, 3) - (searchWindowSize(3)-1)/2)
continue;
end
P = Img(indY-(searchWindowSize(2)-1)/2+1 : indY+(searchWindowSize(2)-1)/2+1, ...
indX-(searchWindowSize(1)-1)/2+1 : indX+(searchWindowSize(1)-1)/2+1, ...
indZ-(searchWindowSize(3)-1)/2+1 : indZ+(searchWindowSize(3)-1)/2+1);
ncc3 = normxcorr3_multi(squeeze(AtlasPatches(v, :, :, :, :)), P, 'valid');
ncc3(ncc3<0.5) = 0;
ncc = squeeze(ncc3(:, :, :, 1) + ncc3(:, :, :, 2) + ncc3(:, :, :, 3));
end %%
h = fspecial3('gaussian',size(ncc), 2);
% ncc(ncc<0.5) = 0;
[~, ind] = max(ncc(:).* h(:));
if(ncc(ind(1)) == 0 || length(ind)~=1)
disp(v, :) = [0, 0, 0];
else
[sub(2), sub(1), sub(3)] = ind2sub(size(ncc), ind(1));
disp(v, :) = [xImg(indX), yImg(indY), zImg(indZ)] - fvFit.vertices(v, :) ... % displacement that comes from maping the nodes to the closest image voxel
+ (sub - (size(ncc))/2).* [imgInfo.PixelSpacing(1), imgInfo.PixelSpacing(2), imgInfo.SliceThickness];
end
% waitbar(v/numVertices, handle);
end
% close(handle);
clear P;
end
% in case the nodes are out of image boundary, we set the displacement to
% zero, so that we have only prediction from statistical shape model.
indOut = fvFit.vertices(:, 1) < imgInfoRigid.BoundingBox.xMin | fvFit.vertices(:, 1)>imgInfoRigid.BoundingBox.xMax | ...
fvFit.vertices(:, 2) < imgInfoRigid.BoundingBox.yMin | fvFit.vertices(:, 2)>imgInfoRigid.BoundingBox.yMax | ...
fvFit.vertices(:, 3) < imgInfoRigid.BoundingBox.zMin | fvFit.vertices(:, 3)>imgInfoRigid.BoundingBox.zMax;
disp(indOut, 1) = 0;
disp(indOut, 2) = 0;
disp(indOut, 3) = 0;
% option.niter_averaging = 2;
% disp = perform_mesh_smoothing(fvFit.faces, fvFit.vertices, disp, option);
option.niter_averaging = 1;
disp = perform_mesh_smoothing(fvFit.faces, fvFit.vertices, disp, option);
newDisp = smoothBasedonSuperNodes(fvFit, disp, model.superNodes);
finalSurface.vertices = fvFit.vertices + newDisp;
finalSurface.vertices = bsxfun(@minus, finalSurface.vertices, rigidTransform.c(1, :)) * rigidTransform.T';
finalSurface.faces = fvFit.faces;
finalSurface.vertices = taubinsmooth( finalSurface.faces, finalSurface.vertices, 15);
% finalSurface.vertices = fvFit.vertices + disp;
% finalSurface.vertices = bsxfun(@minus, finalSurface.vertices, rigidTransform.c(1, :)) * rigidTransform.T';
% finalSurface.faces = fvFit.faces;
%t = toc;
%display (['The sample was segmented in ' num2str(t)])
end
c4science · Help