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:

newPositionsContactCandidates = zeros(3,length(contactCandidate));

% 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');
    
    % 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))
        centroidPositionContact = ccProps(componentIdxsContacts(j)).Centroid;
        newContactPosition = coeff * centroidPositionContact' + centroidPositionHead';
   
        disp(['new Contact Position at ' num2str(j) ' is: ' num2str(newContactPosition')]);
        
        newPositionsContactCandidates(:,j) = newContactPosition;
    end
    
    % compute the centroid of all contact points together to later get the
    % direction in which we define our estimated electrode
    centroidOfAllContacts = mean(newPositionsContactCandidates, 2);

    % new origin is centroid point of head and new (local) coordinate
    % system is coeff from the pca
    
    % compute vector between the two points
    vectorBetweenHeadCentroidCentroidOfAll = centroidOfAllContacts - centroidPositionHead';
    
    % compute both "end points" on the bolt head, were the x-axis of our
    % new coordinate system is going through
    % then we can use them to determine, to which end the mean is nearer
    % and therefore we now on which side of the bolt head the contact
    % candidates lie
    endOneOfBoltHead = coeff(:,1) * latent(1,:) + centroidPositionHead';
    endTwoOfBoltHead = coeff(:,1) * latent(1,:) * (-1.0) + centroidPositionHead';
    
    % compute the lengths of the vectors between both endpoints and the
    % centroid point of the contact candidates
    distMeanToEndOne = norm(centroidOfAllContacts - endOneOfBoltHead);
    distMeanToEndTwo = norm(centroidOfAllContacts - endTwoOfBoltHead);
    
    % we set the origin of our future line to the end which is furthest
    % away from the centroid point and go with it in the direction of the
    % other endpoint, because now we know, there are all the contact
    % candidates
    if(distMeanToEndOne < distMeanToEndTwo) 
        originOfLine = endTwoOfBoltHead;
        secondPointOnLine = endOneOfBoltHead;
        
    else
        originOfLine = endOneOfBoltHead;
        secondPointOnLine = endTwoOfBoltHead;
    end
    
    % can this be 3-dimensional
    % baseline = polyfit(originOfLine, secondPointOnLine, 1);
    
    % vector pointing from origin point of line to second one
    originToSecondPoint = secondPointOnLine - originOfLine;
    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;
    % 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);
    
    % iterate through all contact candidates to find all of them, that have
    % a smaller distance to the baseline than some defined threshold
    % TODO: sort out possible problems with voxels/mm between line and
    % point
    for k = 1:(length(contactCandidate))
        contactPoint = newPositionsContactCandidates(:,k);
        xValueOfContactPoint = contactPoint(1,:);
        % wanted to compute distance between point and corresponding point
        % on the baseline, but there may be conflicts with position of
        % point and baseline measurement in mm
        dist = norm(componentIdxsContacts(k) - baseline(k));
        
        if(dist < threshold)
            %add point to matrix
        end
    end
    
    
    figure
    plotv(baseline);
    
    figure
    scatterMatrix3(newPositionsContactCandidates);
    hold on
    scatterMatrix3(centroidOfAllContacts');
    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


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