//WIP changed some passages, added new functions to compute centroids

sEEG_Segmentation_Script.m

+% 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, '-');