function [helixAngles,penetrationMask,X] = fieldToHelixAngle2(field, mask)
    % 2d field plane (MxNx3)
    % 3d mask (MxN)
    maxLength = 900;
    lvcenter = computeLVCenter(mask);
    mask = imfill(mask,'holes');
    boundary = mask - imerode(mask,strel('disk',1));
    X = generateBoundaryPoints(boundary,lvcenter);
    fprintf('Generated boundary points\n');
    %heartNormalField = computeNormalField(X,imfill(mask,'holes'),250);
    heartNormalField = computeNormalFieldSingleFit(X,imfill(mask,'holes'));
    fprintf('Generated heartNormal field\n');
    helixAngles = zeros(length(X),maxLength);
    threshold = 100;
    penetrationMask = zeros(size(mask));
    penetrationMask(lvcenter) = 1;
    for pt = 1:length(X)  
        startPoint = X(pt,1:2);
        dr = squeeze(heartNormalField(startPoint(2),startPoint(1),:));
        dr = dr';
        count = 0;
        iter = 0;
        %line = startPoint;
        previous = startPoint;
        current = startPoint;
        normal = [dr,0];
        tangent = [-normal(2),normal(1),0];
        axial = cross(normal, tangent);
        %create frame where each row is a unit vector
        frame = [normal;tangent;axial];
        orientation = squeeze(field(startPoint(2),startPoint(1),:));
        helix_angle_row = nan*zeros(1,maxLength);
        k = 1;
        helix_angle_row(k) = computeAngle(frame,orientation); 
        %helix_angle = [];
        while ((count < threshold) && (iter < maxLength))
             current = current + dr;
             gridPoint = round(current);
             iter = iter + 1;
             if(gridPoint(1) == previous(1) && gridPoint(2) == previous(2))
                 count = count+1;
                 iter = iter - 1;
                 continue;
             end
             if (any(~isfinite(gridPoint)))
                break;
             end
             if (gridPoint(1) < 1 || gridPoint(2) < 1  || ...
                 gridPoint(1) >= size(mask,2) || gridPoint(2) >= size(mask,1) )
                break;
             end
             previous = gridPoint;
             if( mask(gridPoint(2),gridPoint(1)) < 1)
                 count = count + 1;
             else
                 count = 0;
             end
             orientation = squeeze(field(gridPoint(2),gridPoint(1),:));
             computed_angle = computeAngle(frame,orientation);
             if ((computed_angle > -1.57)  && (computed_angle < 1.57))
                helix_angle_row(k) = computed_angle;
             else
                 helix_angle_row(k) = nan;
             end

             k = k + 1;
             penetrationMask(gridPoint(2),gridPoint(1)) = pt;
        end
        helixAngles(pt,:) = helix_angle_row;
    end
end