diff --git a/sEEG_Segmentation_Script.m b/sEEG_Segmentation_Script.m
index c9a3d09f696d77392ec77f4d1781149fdef7ee43..fbe7b3a05d17d1eb9c51a55596cd4c1cdde967a5 100644
--- a/sEEG_Segmentation_Script.m
+++ b/sEEG_Segmentation_Script.m
@@ -1,18 +1,9 @@
-% PCA
+%% load the data (CT and BrainMask)
%CT_post = NiftiMod('Seeg_chinese_1_before_final.nii');
%CT_post = NiftiMod('seeg_seg_bul_inside_binary.nii');
CT_post = NiftiMod('seeg_electrode_additional_tryout_high_threshold_2.nii.gz');
-Brain_mask = NiftiMod('brain_mask.nii.gz');
-Brain_mask_img = Brain_mask.img;
-
-bm = bwconncomp(Brain_mask_img,26);
-
-bmProps = regionprops(bm, Brain_mask.img, 'Area', 'PixelIdxList', 'PixelList', 'PixelValues', 'BoundingBox', 'Centroid'); % TODO not really needed as that info is already in cc
-[areas, idxs] = sort([bmProps.Area], 'descend'); % sort by size
-
-bm_centroid = bmProps(idxs(:,1)).Centroid
-
-
+Brain_mask = NiftiMod('brain_mask_chinese_2.nii.gz');
+brainMaskImg = Brain_mask.img;
%% CONSTANTS
% for the contact region (from the paper)
@@ -27,19 +18,18 @@ MAXVOXELNUMBER_CONTACTS = 150;
MINVOXELNUMBER_HEAD = 200;
MAXVOXELNUMBER_HEAD = 2000;
-%% TODO: rename brainmask files in all folders to not get them confused!!!
-%
-% CT_post = NiftiMod('seeg_electrode_additional_tryout_high_threshold_2.nii.gz');
-% Brain_mask = NiftiMod('brain_mask.nii.gz');
-% Brain_mask_img = Brain_mask.img;
-%
-% bm = bwconncomp(Brain_mask_img,26);
-%
-% bmProps = regionprops(bm, Brain_mask.img, 'Area', 'PixelIdxList', 'PixelList', 'PixelValues', 'BoundingBox', 'Centroid');
-% [areas, idxs] = sort([bmProps.Area], 'descend'); % sort by size
-%
-% bm_centroid = bmProps(idxs(:,1)).Centroid
+% define areas of contact regions and heads
+% TODO: define areas at beginning of code
+contactCandidate = areas(areas >= MINVOXELNUMBER_CONTACTS & areas <= MAXVOXELNUMBER_CONTACTS);
+headCandidate = areas(areas >= MINVOXELNUMBER_HEAD & areas <= MAXVOXELNUMBER_HEAD);
+% get all possible candidates
+%componentIdxs = idxs((areas >= MINVOXELNUMBER_CONTACTS) &( areas <= MAXVOXELNUMBER_CONTACTS) | (areas >= MINVOXELNUMBER_HEAD) & (areas <= MAXVOXELNUMBER_HEAD));
+
+componentIdxsContacts = idxs(areas >= MINVOXELNUMBER_CONTACTS & areas <= MAXVOXELNUMBER_CONTACTS);
+componentIdxsHeads = idxs(areas >= MINVOXELNUMBER_HEAD & areas <= MAXVOXELNUMBER_HEAD);
+
+%% TODO: rename brainmask files in all folders to not get them confused!!!
%% TODO: find some solution for electrodes, where the contacts regions are "melted" together and no seperate
%% regions and head are defined
@@ -47,18 +37,19 @@ MAXVOXELNUMBER_HEAD = 2000;
%% get nifti img to binary img, so we can apply some functions
-maskedImg = CT_post.img;
-
-% structEle = strel('sphere', ceil(3 / max(CT_post.voxsize)) ); % make sure brain mask contains no skull
-% brainMask = imerode(brainMask,structEle);
-% maskedImg(~(brainMask)) = NaN;
+binCTImg = CT_post.img;
-% change accordingly
-croppedImg = expandBrainMask(maskedImg);
+% Expand the brain mask to get head bolts of electrodes and get rid of
+% artefacts around brain (pre-selection of candidate)
+% Step 1
+[croppedImg, expBM] = expandBrainMask(binCTImg, brainMaskImg);
-cc = bwconncomp(maskedImg,26); % 26 for 3D images
+% Compute the connected components of the new expanded brain mask
+% cc = bwconncomp(binCTImg,26); % 26 for 3D images (plotting the artefacts)
+cc = bwconncomp(croppedImg,26);
disp([num2str(cc.NumObjects) ' potential metal components detected within brain.']);
+% Sort out the area and idex
ccProps = regionprops(cc, CT_post.img, 'Area', 'PixelIdxList', 'PixelList', 'PixelValues', 'BoundingBox', 'Centroid'); % TODO not really needed as that info is already in cc
[areas, idxs] = sort([ccProps.Area], 'descend'); % sort by size
@@ -66,16 +57,6 @@ ccProps = regionprops(cc, CT_post.img, 'Area', 'PixelIdxList', 'PixelList', 'Pix
elecIdxs = [];
-% define areas of contact regions and heads
-% TODO: define areas at beginning of code
-contactCandidate = areas(areas >= MINVOXELNUMBER_CONTACTS & areas <= MAXVOXELNUMBER_CONTACTS);
-headCandidate = areas(areas >= MINVOXELNUMBER_HEAD & areas <= MAXVOXELNUMBER_HEAD);
-
-% get all possible candidates
-%componentIdxs = idxs((areas >= MINVOXELNUMBER_CONTACTS) &( areas <= MAXVOXELNUMBER_CONTACTS) | (areas >= MINVOXELNUMBER_HEAD) & (areas <= MAXVOXELNUMBER_HEAD));
-
-componentIdxsContacts = idxs(areas >= MINVOXELNUMBER_CONTACTS & areas <= MAXVOXELNUMBER_CONTACTS);
-componentIdxsHeads = idxs(areas >= MINVOXELNUMBER_HEAD & areas <= MAXVOXELNUMBER_HEAD);
%% PCA to get biggest eigenvalue from latent and principle direction from coeff matrix
% TODO:
@@ -83,8 +64,7 @@ componentIdxsHeads = idxs(areas >= MINVOXELNUMBER_HEAD & areas <= MAXVOXELNUMBER
newPositionsContactCandidates = zeros(3,length(contactCandidate));
newHeadPositions = zeros(3, length(headCandidate));
-% this should be the same as length of componentIdxs?
-numCandidates = length(contactCandidate)+length(headCandidate);
+% Compute the PCA for head candidate
for i = 1:(length(headCandidate))
[coeff,~,latent] = pca(ccProps(componentIdxsHeads(i)).PixelList .* repmat(CT_post.voxsize', length(ccProps(componentIdxsHeads(i)).PixelList) ,1));
@@ -98,8 +78,6 @@ for i = 1:(length(headCandidate))
disp(['centroid at ' num2str(i) ' is: ' num2str(centroidPositionHead)]);
- % shiftcoord = global2localcoord([0;0;0], 'rr', centroidPositionHead');
-
% get the centroids of the contactCandidates and shift them by the
% vector of the centroid of the current head candidate
for j = 1:(length(contactCandidate))
@@ -152,12 +130,12 @@ for i = 1:(length(headCandidate))
lengthBetweenPoints = norm(originToSecondPoint);
% get unit vector along the line
originToSecondPoint = originToSecondPoint / lengthBetweenPoints;
- % compute y = a * x + b to get the line, where a is 0 < a < lengthBetweenPoints
- % baseline = originOfLine + (lengthBetweenPoints/2) * originToSecondPoint;
+ % compute y = a * x + b to get the line
% description of our baseline going out from one end of the bolt through the other end into the head space
% with the contact candidates
baseline = @(t)(originOfLine + originToSecondPoint * t);
-
+
+ % matrix to store all candidate inside the tube
% concreteCandidatesForElectrode = zeros(3,:);
% iterate through all contact candidates to find all of them, that have a smaller distance to the baseline
@@ -196,37 +174,21 @@ for i = 1:(length(headCandidate))
%
% figure
% plotv(baseline);
-
- figure
-
-
+
figure
axis([0 500 0 500 0 500])
axis equal
scatterMatrix3(newPositionsContactCandidates);
+ isosurface(binCTImg, 'blue');
+ isosurface(expBM, 'yellow');
+ isosurface(brainMaskImg, 'red');
+ alpha(0.3);
hold on
scatterMatrix3(centroidOfAllContacts');
scatterMatrix3(thresholdElecMaxValue', {}, '*');
scatterMatrix3(endOneOfBoltHead');
scatterMatrix3(endTwoOfBoltHead');
- % hold off
- %
- % de=det(coeff);
- %
- % newcoeff = [coeff(:,2), -coeff(:,3), coeff(:,1)];
- %
- % figure;
- % hold on
- % x=[0,coeff(:,1)];
- % y=[0,coeff(:,2)];
- % z=[0,coeff(:,3)];
- % quiver3(0,0,0, coeff(:,1), coeff(:,2), coeff(:,3));
- % plot3(x);
- % plot3(coeff(:,2));
- % plot3(coeff(:,3));
- % plot3(coeff(:,2), coeff(:,3), coeff(:,1), '*-', 'LineWidth', 10);
- %
- % te=det(newcoeff);
+
end
scatterMatrix3(newHeadPositions');