diff --git a/ReadMe.txt b/ReadMe.txt new file mode 100755 index 0000000..a9c5ff7 --- /dev/null +++ b/ReadMe.txt @@ -0,0 +1,19 @@ +This folder contains following information: + +“sampleData” contains an example image for segmentation. + +“functions” folder contains all functions I use for segmentation. + +“00_vertices” folder contains all preselected nodes of the mesh that are going to be used for automatic landmark detection. + +“loadModel.m” load the model and all information that is necessary for segmentation. + +“autoSegment.m” is the function that does the segmentation. You should provide the loaded image with some input options (information on the header of functions) and as output it will provide the surface mesh as a struct with the fields of .vertices (coordinates of the nodes) and .faces (connectivity information). + +“doMeasurements.m” is the function that gets the segmented bone as input and as output provides the measurements in a struct. + +To see how the model and images should be prepared, you can find a MATLAB script in the same folder “howToPrepareData.m” that represents how you can call different functions. + +The files with the name of “tmp” are saved during process of segmentation and are replaced for every call to “autoSegment” function. + +"Protocol Segmentation Scapula.docx" is a protocol (format word) which explains how to segmented a scapula from the beginning to the end. It is situated in the folder: "shoulder\methods\documentations\protocols" \ No newline at end of file diff --git a/autoSegment.m b/autoSegment.m new file mode 100755 index 0000000..0c2b3d8 --- /dev/null +++ b/autoSegment.m @@ -0,0 +1,537 @@ +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; + if(mod(n, 1.e5) == 0) + display('I am waiting for Video_and_Image_Block_Set Toolbox license!'); + 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) 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 \ No newline at end of file diff --git a/doMeasurements.m b/doMeasurements.m new file mode 100755 index 0000000..4d634f1 --- /dev/null +++ b/doMeasurements.m @@ -0,0 +1,374 @@ +function out = doMeasurements(fv, options) +% it works only for the mesh that was segmented using our segmentation +% method that provides correspondences to our statistical shape model. +% fv.faces: contains the connectivity information of the surface mesh (25000 x 3) +% fv.vertices: information on the vertices of the mesh (12502 x 3) +% options stores some important information regarding the bone: +% .side Can be 'left' or 'right' +% .curvature_smoothing Curvature smoothing parameter. Default value is 10 +% .verb step-by-step reporting. by default it is set to false +% out is a struct contains different measurements + +%% check the input's validity +if(~isfield(fv, 'faces')) + error('The provided mesh is not in correct format: the field .faces is not provided'); +end +if(~isfield(fv, 'vertices')) + error('The provided mesh is not in correct format: the field .vertices is not provided'); +end + +sz_F = size(fv.faces); +sz_V = size(fv.vertices); + +if (numel(sz_F)~=2 || sum(sz_F == [25000, 3])~=2 || numel(sz_V)~=2 || sum(sz_V == [12502, 3])~=2) + error('The provided mesh is not in correct format.'); +end + + +if (nargin == 1) + options.side = 'right'; + options.curvature_smoothing = 10; + options.verb = 0;% 0: do not print progress bar; 1: print progress bar in the consol +else + if(~isfield(options, 'side')) + options.side = 'right'; + end + if(~isfield(options, 'curvature_smoothing')) + options.curvature_smoothing = 10; + end + if(~isfield(options, 'verb')) + options.verb = 0; + end +end + + + + +if(strcmp(options.side, 'right')) + reflect = true; +else + reflect = false; +end +%% The important vertices needed for measurements +[selectedVertices.ScapulaToExclude, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_ScapulaToExclude.ini']); +[selectedVertices.ScapulaGroove, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_Groove.ini']); +[selectedVertices.ScapulaPillar, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_Pillar.ini']); +[selectedVertices.GlenoidToExclude, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_GlenoidToExclude.ini']); +[selectedVertices.Spinoglenoidnotch, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_Spinoglenoidnotch.ini']); +[selectedVertices.GlenoidSurface, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_GlenoidSurface.ini']); + selectedVertices.ScapulaExtremities = load(['.' filesep '00_vertices' filesep 'indiceVertices_ScapulaExtremities.txt']); + +Sref = simplePlyRead(['.' filesep '00_vertices' filesep 'ReferenceShape.ply']); + + + %% Prepare sample for measurements +S.FV = fv; +% [~, S.FV.vertices] = procrustes(Sref.vertices, fv.vertices, 'scaling', false, 'reflection', reflect); +[~, S.FV.vertices, T] = procrustes(Sref.vertices, fv.vertices, 'scaling', false, 'reflection', reflect); +%try to compute the curvature, if it does not work, we slightly smooth +%the surface and retry +temp = 1; +while(temp) + try + temp = 0; + [S.Curvature.DirectionMinCurvature,... + S.Curvature.DirectionMaxCurvature,... + S.Curvature.MinCurvature,... + S.Curvature.MaxCurvature,... + S.Curvature.MeanCurvature,... + S.Curvature.GaussianCurvature,... + S.Curvature.Normal] = ... + compute_curvature(S.FV.vertices, S.FV.faces, options); + catch %...if an error occurs, we smooth the face and recompute the curvature + S.FV = smoothpatch(S.FV, 1, 1, 0.0001); + if(options.verb) + display('\tThe registered shape is smoothed to compute the curvature\n'); + end + temp = 1; + end +end + +%% find Landmarks based on Curvatures +temp_FV = S.FV; +temp_Curvature = S.Curvature; +temp_Landmarks.Extremities = temp_FV.vertices(selectedVertices.ScapulaExtremities,:); +% ----------- find the spinoglenoid notch ----------------- +%select the vertice which has the smallest gaussian curvature +[~,temp] = min(temp_Curvature.GaussianCurvature(selectedVertices.Spinoglenoidnotch)); +temp = selectedVertices.Spinoglenoidnotch(temp); +temp_Landmarks.Extremities(6,:) = temp_FV.vertices(temp,:); + + +% ------ find the scapula groove landmarks ------ +temp_dist = [12 44 28 20 36]; +skipInd = []; +for iL=1:5 + %get all the candidate vertices + temp_indices = selectedVertices.ScapulaGroove; + temp_vertices = temp_FV.vertices(temp_indices,:); + + %get the vector front-scapula->back-scapula + [temp_n,~] = fitLineToPoints(temp_vertices); + temp_n = temp_n'; + temp_Lfront = temp_Landmarks.Extremities(6,:); + %from the front-scapula, go X mm in direction of the back, and remove + %all candidate vertices behind + temp_origine = temp_Lfront + (temp_dist(iL)-2)*temp_n; + for i=length(temp_indices):-1:1 + temp_p = temp_vertices(i,:)-temp_origine; + if(dot(temp_p,temp_n)<0) + temp_indices(i)=[]; + temp_vertices(i,:)=[]; + end + end + + %from the front-scapula, go Y mm in direction of the back, and remove + %all candidate vertices ahead + temp_origine = temp_Lfront + (temp_dist(iL)+2)*temp_n; + for i=length(temp_indices):-1:1 + temp_p = temp_vertices(i,:)-temp_origine; + if(dot(temp_p,temp_n)>0) + temp_indices(i)=[]; + temp_vertices(i,:)=[]; + end + end + + %get the best vertice based on the min curvature + if ~isempty(temp_indices) + [temp, i] = min(temp_Curvature.MinCurvature(temp_indices)); + temp_Landmarks.Groove(iL,:) = temp_vertices(i,:); + else + skipInd = [skipInd iL]; + end +end +if(numel(skipInd) == 5) + display('The shape is distorted! Could not find enough landmarks on the Groove'); + return; +else + temp_Landmarks.Groove(skipInd, :) = []; +end + +% --------- find the scapula pillar landmarks ------ +temp_dist = [30 58 44 37 51]; +for iL=1:5 + %get all the candidate vertices + temp_indices = selectedVertices.ScapulaPillar; + temp_vertices = temp_FV.vertices(temp_indices,:); + + %get the vector up->down + [temp_n,~] = fitLineToPoints(temp_vertices); + temp_n = temp_n'; + + %check the vector point in the right direction + temp_up = temp_Landmarks.Extremities(6,:); + if(dot(temp_n, temp_FV.vertices(selectedVertices.ScapulaPillar(1,:))-temp_up)<0) + temp_n =-temp_n; + end + %from the front-scapula, go X mm in direction of the back, and remove + %all candidate vertices behind + temp_origine = temp_up + (temp_dist(iL)-2)*temp_n; + for i=length(temp_indices):-1:1 + temp_p = temp_vertices(i,:)-temp_origine; + if(dot(temp_p,temp_n)<0) + temp_indices(i)=[]; + temp_vertices(i,:)=[]; + end + end + + %from the front-scapula, go Y mm in direction of the back, and remove + %all candidate vertices ahead + temp_origine = temp_up + (temp_dist(iL)+2)*temp_n; + for i=length(temp_indices):-1:1 + temp_p = temp_vertices(i,:)-temp_origine; + if(dot(temp_p,temp_n)>0) + temp_indices(i)=[]; + temp_vertices(i,:)=[]; + end + end + + %get the best vertice based on the max curvature + [temp, i] = max(temp_Curvature.MaxCurvature(temp_indices)); + temp_Landmarks.Pillar(iL,:) = temp_vertices(i,:); +end +% ------------- save everyting (and glenoid surface) in the structure ---------------- +temp_v = selectedVertices.GlenoidSurface(temp_Curvature.MaxCurvature(selectedVertices.GlenoidSurface)<1); +S.indiceVertices.GlenoidSurface = temp_v; +S.Landmarks = temp_Landmarks; + +%% compute the glenoid orientations + +S.GlenoidOrientation = ... + computeGlenoidOrientation(... + S.FV,... + S.indiceVertices.GlenoidSurface,... + S.Landmarks); + +out = rmfield(S, 'FV'); +out = rmfield(out, 'indiceVertices'); + +out.Landmarks.Pillar = bsxfun(@minus, out.Landmarks.Pillar, T.c(1, :)) * T.T'; +out.Landmarks.Extremities= bsxfun(@minus, out.Landmarks.Extremities, T.c(1, :)) * T.T'; +out.Landmarks.Groove= bsxfun(@minus, out.Landmarks.Groove, T.c(1, :)) * T.T'; + +end + + + +%% function to load vertices numbers +function [iVertices, iFaces] = loadOpenFlipperFV(filepath) + +%% read the file +fileID = fopen(filepath,'r'); +temp = textscan(fileID,'%s','Delimiter','\n'); +fclose(fileID); + + +%% extract vertices + +%select the line which list all the vertices of interest +temp_Vertices = temp{1}{7}; +temp_Vertices(1:16) = []; %remove the text 'vertexSelection=' +if(isempty(temp_Vertices)==0) + temp_Vertices(end) = []; %remove the last character + + %split the string + temp_Vertices = strsplit(temp_Vertices,';')'; + + %transform the string in nummer + iVertices = zeros(length(temp_Vertices),1); + for i=1:length(temp_Vertices) + iVertices(i,1) = str2num(temp_Vertices{i}) + 1; %add 1 because the vertices start at 0 in OpenFlipper and 1 in Matlab + end +else + iVertices = []; +end + + +%% extract faces +%select the line which list all the vertices of interest +temp_Faces = temp{1}{3}; +temp_Faces(1:14) = []; %remove the text 'faceSelection=' +if(isempty(temp_Faces)==0) + temp_Faces(end) = []; %remove the last character + + %split the string + temp_Faces = strsplit(temp_Faces,';')'; + + %transform the string in nummer + iFaces = zeros(length(temp_Faces),1); + for i=1:length(temp_Faces) + iFaces(i,1) = str2num(temp_Faces{i}) + 1; %add 1 because the vertices start at 0 in OpenFlipper and 1 in Matlab + end + +else + iFaces = []; +end + +end + + + + +function Results = computeGlenoidOrientation(FV, indiceVertices_GlenoidSurface, Landmarks) + +% *********************************************************** +% Notes: landmarks.extremities +% (1) Bottom of the scapula +% (2) Back of the scapula +% (3) Tip of acromion +% (4) Top coracoid +% (5) Back coracoid +% (6) Front scapula +% ************************************************************ + +%load Data + +% coordinate of vertices on the Glenoid Surface +GlenoidSurface.vertices = FV.vertices(indiceVertices_GlenoidSurface,:); + +% ----------------------- compute the coordinate system of the scapula ----------------------- + +% Find the transverse scapular plane (fit plane to landmarks.ScapulaPillar and landmarks.ScapulaGroove) +[ScapulaPlane.Normal, ScapulaPlane.Center] = fitPlaneToPoints([Landmarks.Groove; Landmarks.Pillar]); + +%define the X coord syst +CSYS.X = ScapulaPlane.Normal; + +%check the X axis goes in the right direction +temp = Landmarks.Extremities(3,:) - Landmarks.Extremities(4,:); +if(sign(dot(CSYS.X,temp))<0) + CSYS.X = -CSYS.X; +end + +%define the Z coord syst +temp = projectPointsToPlane(Landmarks.Groove,ScapulaPlane.Normal',[0 0 0]); % project landmarks.ScapulaGroove points on scapular plane +[CSYS.Z, ~] = fitLineToPoints(temp); % fit a line on the projected landmarks.ScapulaGroove +CSYS.Z = CSYS.Z/norm(CSYS.Z); %normalize the vector + + +%check the Z axis goes in the right direction +temp = Landmarks.Extremities(6,:) - mean(Landmarks.Groove); +if(sign(dot(CSYS.Z,temp))<0) + CSYS.Z = -CSYS.Z; +end + +%define the Y coord syst +CSYS.Y = cross(CSYS.Z,CSYS.X); % Y + + +% ----------------------- compute the Glenoid center ----------------------- + +% temp = zeros(size(GlenoidSurface.Vertices)); +% temp_distance3D = zeros(size(GlenoidSurface.Vertices,1),1); + +GlenoidSurface.Center = mean(GlenoidSurface.vertices); + +%returns the indices of the closest points in X for each point in Y +[temp_iSelectedVertice, ~]= dsearchn(GlenoidSurface.vertices,GlenoidSurface.Center); +GlenoidSurface.Center = GlenoidSurface.vertices(temp_iSelectedVertice,:); + + + +% ----------------------- Fit Sphere on Glenoid Surface ----------------------- + +[GlenoidSphere.Center,GlenoidSphere.Radius,temp_GlenoidResiduals] = fitSphereToPoints(GlenoidSurface.vertices); +GlenoidSphere.RMSE = norm(temp_GlenoidResiduals)/sqrt(length(temp_GlenoidResiduals)); + + + +% Glenoid version 3D (spherical coordinates) + +%centerline +Glenoid.Centerline = GlenoidSphere.Center'-GlenoidSurface.Center; + +%get the centerline vector in the scapula coord system +tempX = Glenoid.Centerline * CSYS.X; %project the centerline on the X axis of the scapula +tempY = Glenoid.Centerline * CSYS.Y; %project the centerline on the Y axis of the scapula +tempZ = Glenoid.Centerline * CSYS.Z; +temp_r = norm(Glenoid.Centerline); %length of the centerline + +%version +Glenoid.Version_inDeg = acosd(tempZ/temp_r); + +%version orientation +Glenoid.VersionOrientation_inDeg = 180-rad2deg(atan2(tempY,tempX)); +if(Glenoid.VersionOrientation_inDeg > 180) + Glenoid.VersionOrientation_inDeg = 360 - Glenoid.VersionOrientation_inDeg; +end + + +% ----------------------- Output ----------------------- + +% Glenoid +Results.Glenoid_Center = GlenoidSurface.Center; +Results.Glenoid_SphereCenter = GlenoidSphere.Center; +Results.Glenoid_SphereRadius = GlenoidSphere.Radius; +Results.Glenoid_Centerline = Glenoid.Centerline; +Results.Glenoid_VersionInDegree = Glenoid.Version_inDeg; +Results.Glenoid_VersionOrientationInDegree = Glenoid.VersionOrientation_inDeg; + +Results.CoordSys.X = CSYS.X; +Results.CoordSys.Y = CSYS.Y; +Results.CoordSys.Z = CSYS.Z; + +end diff --git a/howToPrepareData.m b/howToPrepareData.m new file mode 100755 index 0000000..2c2b456 --- /dev/null +++ b/howToPrepareData.m @@ -0,0 +1,145 @@ +%% How to call autoSegment function +%% init +clear all +close all +clc + +%% parameters of the segmentation. These parameters are used and you can use them as they are: +% opt.isLeft = true; +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! +opt.searchWindowSize = [65, 65, 45]; + +% The folder that contains images to segment +% folderIn = './sampleData/samplesLeft/'; +folderIn = './sampleData/samplesRight/'; +folderOut = folderIn; +modelPath = '../Model/'; + +% This line adds sift feature extraction toolbox to the path +run('../Libraries/SIFT3D-master/build/lib/wrappers/matlab/setupSift3D.m') + +% all nenessary functions are stored in the folder fonctions. We add these +% functions to the path +addpath(genpath('./functions')); + +%% load model and extract the sift features on the meanShape +model = loadModel(modelPath, opt.isLeft); + + +%% load samples one by one and perform segmentation and then save the result +% This example assumes that all images are listed in a single folder. or in +% folders inside one folder specified in folderIn. For example, let's say +% folderIn is structured as follows: +% folderIn +% Img1.am +% Img2 +% Img2.mhd +% Img2.raw +% Img3 +% Img3_0001.dcm +% Img3_0002.dcm +% ... +% Img3_0250.dcm +% Img4.am +% ... + +% If the structure of your folders are different, then you should adapt +% this folder + + +% list of all images in the folder +bones_in = dir([folderIn '*']); + +%Exclude "." and ".." from the list. Put dicom images in a folder! +bones_in = bones_in(3:end); +for b = 1 : size(bones_in, 1) + tic + boneName = [folderIn bones_in(b).name]; + + %if the folder contains DICOMDIR or each image is in one directory! + if (isdir(boneName)) + imageFile = dir(boneName); + imageFile = imageFile(3:end); + for file_counter = 1:size(imageFile, 1) + if(strcmp(imageFile(file_counter).name(end-2:end), 'mhd') || strcmp(imageFile(file_counter).name(end-2:end), 'dcm')) + boneName = [boneName, filesep imageFile(3).name]; + break; + end + end + end + + fname = boneName; + ind = strfind(fname, '.'); + if(isempty(ind)) + continue; + end + ind = ind(end); + fileFormat = fname(ind+1:end); + % right now I have code for ".am", ".mhd" and ".dcm" file formats. You + % can add more cases here to read data. + switch fileFormat + case 'am' + try + [imgInfo, Img] = LoadData_Amira(boneName); + catch + display([bones_in(b).name ': Problem in reading the file. Check the file and try again later!']); + end + case 'mhd' + try + [imgMHD, info]=read_mhd(boneName); + catch + display([bones_in(b).name ': Problem in reading the file. Check the file and try again later!']); + end + eval(['Img = ' class(imgMHD.data) '(zeros(size(imgMHD.data, 2), size(imgMHD.data, 1), size(imgMHD.data, 3)));']); + for ii = 1:size(imgMHD.data, 3) + Img(:, :, ii) = imgMHD.data(:, :, ii)'; + end + imgInfo.BoundingBox.xMin = imgMHD.origin(1); + imgInfo.BoundingBox.yMin = imgMHD.origin(2); + imgInfo.BoundingBox.zMin = imgMHD.origin(3); + imgInfo.PixelSpacing(1) = imgMHD.spacing(1); + imgInfo.PixelSpacing(2) = imgMHD.spacing(2); + imgInfo.SliceThickness = imgMHD.spacing(3); + imgInfo.xSize = imgMHD.size(1); + imgInfo.ySize = imgMHD.size(2); + imgInfo.zSize = imgMHD.size(3); + imgInfo.BoundingBox.xMax = imgMHD.origin(1) + (imgMHD.size(1)-1) * imgInfo.PixelSpacing(1); + imgInfo.BoundingBox.yMax = imgMHD.origin(2) + (imgMHD.size(2)-1) * imgInfo.PixelSpacing(2); + imgInfo.BoundingBox.zMax = imgMHD.origin(3) + (imgMHD.size(3)-1) * imgInfo.SliceThickness; + + if imgInfo.BoundingBox.zMax - imgInfo.BoundingBox.zMin > 700 + display ([ bones_in(b).name ': skipped processing, due to the large size of scan, please crop this image and then try again (' num2str(size(imgMHD.data)) ')']); + continue; + end + clear imgMHD; + case 'dcm' + dicomDir = strsplit(boneName, filesep); + dicomDir = strjoin(dicomDir(1, 1:end-1), filesep); + [Img, imgInfo] = readDicomStack(dicomDir); + if (isempty(Img)) + display ([ bones_in(b).name ': Could not read the dicom images.']); + continue; + end + otherwise + continue; + end + fname = bones_in(b).name; + ind = strfind(fname, '.'); + if(isempty(ind)) + ind = length(fname)+1; + else + ind = ind(end); + end + + finalSurface = autoSegment(Img, imgInfo, model, opt); + simplePlyWrite(finalSurface, [folderOut fname(1:ind-1) '_auto.ply']); + measurements = doMeasurements(finalSurface); + save( [folderOut fname(1:ind-1) '.am'],'measurements'); \ No newline at end of file diff --git a/howToPrepareData_AT.m b/howToPrepareData_AT.m new file mode 100755 index 0000000..94e1f6b --- /dev/null +++ b/howToPrepareData_AT.m @@ -0,0 +1,180 @@ +%% How to call autoSegment function +%% init + clear all + close all +% clc + +%% parameters of the segmentation. These parameters are used and you can use them as they are: +% opt.isLeft = true; +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! +opt.searchWindowSize = [65, 65, 45]; + +% The folder that contains images to segment +% folderIn = './sampleData/samplesLeft/'; +% folderIn = './sampleData/samplesRight/'; +folderIn = '../../../../data/N/0/3/N34-85688/CT-N34-85688-1/'; % should check it exist +folderOut = [folderIn 'matlab']; % set output in a matlab folder within folderIn +if ~exist(folderOut, 'dir') % create folder if not exist + mkdir(folderOut); +end + +modelPath = '../Model/'; + +% This line adds sift feature extraction toolbox to the path + run('../Libraries/SIFT3D-master/build/lib/wrappers/matlab/setupSift3D.m') + +% all nenessary functions are stored in the folder fonctions. We add these +% functions to the path +addpath(genpath('./functions')); + +%% load model and extract the sift features on the meanShape +tic; +display('Load SS model...') +model = loadModel(modelPath, opt.isLeft); +toc; + +if(folderIn(end)~=filesep) + folderIn = [folderIn filesep]; +end + +if(folderOut(end)~=filesep) + folderOut = [folderOut filesep]; +end +%% load samples one by one and perform segmentation and then save the result +% This example assumes that all images are listed in a single folder. or in +% folders inside one folder specified in folderIn. For example, let's say +% folderIn is structured as follows: +% folderIn +% Img1.am +% Img2 +% Img2.mhd +% Img2.raw +% Img3 +% Img3_0001.dcm +% Img3_0002.dcm +% ... +% Img3_0250.dcm +% Img4.am +% ... + +% If the structure of your folders are different, then you should adapt +% this folder + + +% list of all images in the folder +bones_in = dir([folderIn '*']); + +%Exclude "." and ".." from the list. Put dicom images in a folder! +bones_in = bones_in(3:end); + +% Restrict to folder named "dicom", for use in /home/shoulder/data/ +bones_in_idx = find(strcmp({bones_in.name},'dicom') == 1); +bones_in = bones_in(bones_in_idx); +% we should check it exists + +disp(['folderIn: ' folderIn]); +disp(['folderOut: ' folderOut]); +disp(bones_in); + + + +for b = 1 : size(bones_in, 1) +disp('Load images...'); +tic; + boneName = [folderIn bones_in(b).name]; + + %if the folder contains DICOMDIR or each image is in one directory! + if (isdir(boneName)) + imageFile = dir(boneName); + imageFile = imageFile(3:end); + for file_counter = 1:size(imageFile, 1) + if(strcmp(imageFile(file_counter).name(end-2:end), 'mhd') || strcmp(imageFile(file_counter).name(end-2:end), 'dcm')) + boneName = [boneName, filesep imageFile(3).name]; + break; + end + end + end + + fname = boneName; + ind = strfind(fname, '.'); + if(isempty(ind)) + continue; + end + ind = ind(end); + fileFormat = fname(ind+1:end); + % right now I have code for ".am", ".mhd" and ".dcm" file formats. You + % can add more cases here to read data. + switch fileFormat + case 'am' + try + [imgInfo, Img] = LoadData_Amira(boneName); + catch + display([bones_in(b).name ': Problem in reading the file. Check the file and try again later!']); + end + case 'mhd' + try + [imgMHD, info]=read_mhd(boneName); + catch + display([bones_in(b).name ': Problem in reading the file. Check the file and try again later!']); + end + eval(['Img = ' class(imgMHD.data) '(zeros(size(imgMHD.data, 2), size(imgMHD.data, 1), size(imgMHD.data, 3)));']); + for ii = 1:size(imgMHD.data, 3) + Img(:, :, ii) = imgMHD.data(:, :, ii)'; + end + imgInfo.BoundingBox.xMin = imgMHD.origin(1); + imgInfo.BoundingBox.yMin = imgMHD.origin(2); + imgInfo.BoundingBox.zMin = imgMHD.origin(3); + imgInfo.PixelSpacing(1) = imgMHD.spacing(1); + imgInfo.PixelSpacing(2) = imgMHD.spacing(2); + imgInfo.SliceThickness = imgMHD.spacing(3); + imgInfo.xSize = imgMHD.size(1); + imgInfo.ySize = imgMHD.size(2); + imgInfo.zSize = imgMHD.size(3); + imgInfo.BoundingBox.xMax = imgMHD.origin(1) + (imgMHD.size(1)-1) * imgInfo.PixelSpacing(1); + imgInfo.BoundingBox.yMax = imgMHD.origin(2) + (imgMHD.size(2)-1) * imgInfo.PixelSpacing(2); + imgInfo.BoundingBox.zMax = imgMHD.origin(3) + (imgMHD.size(3)-1) * imgInfo.SliceThickness; + + if imgInfo.BoundingBox.zMax - imgInfo.BoundingBox.zMin > 700 + display ([ bones_in(b).name ': skipped processing, due to the large size of scan, please crop this image and then try again (' num2str(size(imgMHD.data)) ')']); + continue; + end + clear imgMHD; + case 'dcm' + dicomDir = strsplit(boneName, filesep); + dicomDir = strjoin(dicomDir(1, 1:end-1), filesep); + [Img, imgInfo] = readDicomStack(dicomDir); + if (isempty(Img)) + display ([ bones_in(b).name ': Could not read the dicom images.']); + continue; + end + otherwise + continue; + end + fname = bones_in(b).name; + ind = strfind(fname, '.'); + if(isempty(ind)) + ind = length(fname)+1; + else + ind = ind(end); + end +toc; +disp('Segment CT...'); + + finalSurface = autoSegment(Img, imgInfo, model, opt); +disp('Write surface file...'); +tic; + simplePlyWrite(finalSurface, [folderOut fname(1:ind-1) '_auto.ply']); +toc; +% measurements = doMeasurements(finalSurface); +% save( [folderOut fname(1:ind-1) '.mat'],'measurements'); + + display(['Processing of sample ' fname(1:ind-1) ' is finished.']) +end