From f99cf539d1f05399f36606a4df20924ad9f62064 Mon Sep 17 00:00:00 2001
From: Marie Bofferding <mc.bofferding@web.de>
Date: Tue, 10 Sep 2019 09:44:29 +0200
Subject: [PATCH] //WIP changed some passages, added new functions to compute
centroids
---
sEEG_Segmentation_Script.m | 139 +++++++++++++++++++++++++++++++++++++
1 file changed, 139 insertions(+)
create mode 100644 sEEG_Segmentation_Script.m
diff --git a/sEEG_Segmentation_Script.m b/sEEG_Segmentation_Script.m
new file mode 100644
index 0000000..8335018
--- /dev/null
+++ b/sEEG_Segmentation_Script.m
@@ -0,0 +1,139 @@
+% PCA
+CT_post = NiftiMod('Seeg_chinese_1_before_final.nii');
+
+%% CONSTANTS
+% for the contact region (from the paper)
+% contact region areas have defined interval of voxels, because they are
+% the smallest parts in the data
+MINVOXELNUMBER_CONTACTS = 4; % at least 4 as minimum value
+MAXVOXELNUMBER_CONTACTS = 50;
+
+% for head
+% bolt head areas have defined minimum of voxels, as they are bigger than
+% contact regions
+% a maximum is also defined, as there are sometimes cable artifacts left in
+% the data, so the heads are not the biggest areas in the data
+MINVOXELNUMBER_HEAD = 300;
+MAXVOXELNUMBER_HEAD = 900;
+
+% TODO: find some solution for electrodes, where the contacts regions are
+% "melted" together and no seperate regions and head are defined
+
+%% 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;
+
+cc = bwconncomp(maskedImg,26); % 26 for 3D images
+disp([num2str(cc.NumObjects) ' potential metal components detected within brain.']);
+
+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
+
+% number of index of electrodes
+
+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:
+
+A = [];
+
+% this should be the same as length of componentIdxs?
+numCandidates = length(contactCandidate)+length(headCandidate);
+
+for i = 1:(length(headCandidate))
+ [coeff,~,latent] = pca(ccProps(componentIdxsHeads(i)).PixelList .* repmat(CT_post.voxsize', length(ccProps(componentIdxsHeads(i)).PixelList) ,1));
+ if(length(latent) < 3)
+ continue
+ end
+ latent = sqrt(latent) * 2; % to get full direction
+
+ % get centroids of the heads to shift coordinate system to and get
+ % starting point for electrode line
+ centroidPositionHead = ccProps(componentIdxsHeads(i)).Centroid;
+
+ disp(['centroid at ' num2str(i) ' is: ' num2str(centroidPositionHead)]);
+
+ shiftcoord = global2localcoord([0;0;0], 'rr', centroidPositionHead');
+
+ for j = 1:(length(contactCandidate))
+
+ centroidPositionContact = ccProps(componentIdxsContacts(j)).Centroid;
+
+ newContactPosition = coeff * centroidPositionContact' + centroidPositionHead';
+
+
+
+ end
+
+
+% 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
+
+
+% figure
+% for i = 1:(length(contactCandidate)+length(headCandidate))
+% X = ccProps(componentIdxs(i)).PixelList .* repmat(CT_post.voxsize', length(ccProps(componentIdxs(i)).PixelList) ,1);
+% scatterMatrix3(X)
+% hold on
+% end
+
+
+
+x=coeff(:,1);
+y=coeff(:,2);
+z=coeff(:,3);
+
+inter = intersect(intersect(x,y),z);
+
+plot(x);
+plot(y);
+plot(z);
+
+t=det(coeff);
+
+plot3(x,y,z);
+
+%coeff(:,1)
+
+
+M = [0,0,1;0,0,-1;0,1,0];
+
+a = M(:,1);
+b = M(:,2);
+c = M(:,3);
+
+vec_inter = intersect(intersect(a,b), c);
+
+figure;
+plotv(M, '-');
--
GitLab