function [X] = generateBoundaryPoints(boundary,HeartC)
    % boundary : 3d/2d volume mask which is 1 if boundary pixel else 0
    % HeartC : 1x2 vector center of LV.
    % returns N x 4  matrix X where N = number of boundary points
    % X(n,1:3) = x,y,z coordinate, X(n,4) is the angle w.r.t heart center.
    % The angle computation assumes aligned x/y axes
    linIdx = find(boundary);
    [i,j,k] = ind2sub(size(boundary),linIdx);
    X = [j,i,k];
    %sort by angle..
    for i = 1:length(X)
        xx = X(i,1) - HeartC(1);
        yy = X(i,2) - HeartC(2);
        theta = atan2(yy,xx);
        if theta < 0
            theta = 2*pi + theta;
        end
        X(i,4) = theta;
    end
    X = sortrows(X,4,'descend');
end