Homec4science
Diffusion segmentation (master) Browse

segmentation/howToPrepareData_AT.m
8222343c77d4master

Last Change
Raw File

howToPrepareData_AT.m
View Options

%
% This script can be executed from shell with the command
% matlab -nodisplay -r "howToPrepareData_AT;quit"
%
%% 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
% folderIn = '../../../../data/P/2/2/P221-657235/CT-P221-657235-1/';
% folderIn = '../../../../data/N/2/0/N209-149341/CT-N209-149341-1/';
%folderIn = '../../../../data/P/0/0/P421-430069/CT-P421-430069-1/';
folderIn = '../../../../data/P/0/0/P1-490313/CT-P1-490313-1/';
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
display('Load SS models (right & left)...')
tic;
% model = loadModel(modelPath, opt.isLeft);
modelLeft = loadModel(modelPath, true); % isLeft = true;
modelRight = loadModel(modelPath, false); % isLeft = false;
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(file_counter).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;
% Identify scapula side
disp('Identify scapula side');
tic;
[Left, Right] = IdentifyRightAndLeftScapula(Img, imgInfo, modelLeft, modelRight); % Img and imgInfo is similar to autoSegment function.
toc;
% Segment scapula
disp('Segment CT...');
tic;
if(isempty(Right) && ~isempty(Left))
disp('Image contains only the left scapula.');
opt.isLeft = true;
finalSurface = autoSegment(Img, imgInfo, modelLeft, opt);
end
if(isempty(Left) && ~isempty(Right))
disp('Image contains only the right scapula.');
opt.isLeft = false;
finalSurface = autoSegment(Img, imgInfo, modelRight, opt);
end
if (~isempty(Right) && ~isempty(Left))
disp('Image contains both left and right scapulae.');
opt.isLeft = false;
finalSurface_Right = autoSegment(Right.img, Right.imgInfo, modelRight, opt);
opt.isLeft = true;
finalSurface_Left = autoSegment(Left.img, Left.imgInfo, modelLeft, opt);
end
if (isempty(Right) && isempty(Left))
disp('Could not recognize, which scapulae it is.');
end
toc;
% Write surface file
disp('Write surface file...');
tic;
simplePlyWrite(finalSurface, [folderOut fname(1:ind-1) '_auto.ply']);
toc;
disp('Measure anatomy...');
tic;
measurements = doMeasurements(finalSurface);
% save( [folderOut fname(1:ind-1) '.mat'],'measurements');
save( [folderOut 'anatData.mat'],'measurements');
toc;
disp(['Processing of sample ' fname(1:ind-1) ' is finished.'])
end
c4science · Help