function [helixAngleMap,colored ] = helixAnglesToHelixMap(helixAngles,X,mask,sectorMask)
    if nargin < 4
        sectorMask = mask;
    end
    helixAngleMap = nan*ones(size(mask));
    repeatCount = zeros(size(mask));
    depth = size(helixAngles,2);
    assert(size(helixAngles,1) == size(X,1));
    heartNormalField = computeNormalFieldSingleFit(X,imfill(mask,'holes'));
    fprintf('Generated heartNormal field\n');
    threshold = 100;
    for pt = 1:length(X)  
        startPoint = X(pt,1:2);
        dr = squeeze(heartNormalField(startPoint(2),startPoint(1),:));
        dr = dr';
        %for k = 1:depth
        %ugly code repetion.. need to clean this up
        currentPoint = startPoint;
        previous = currentPoint;
        count = 0;
        iter = 0;
        k = 1;
         while ((count < threshold) && (iter < depth))
             currentPoint = currentPoint + dr;
             gridPoint = round(currentPoint);
             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( sectorMask(gridPoint(2),gridPoint(1)) < 1)
                continue;
             end
            old_value = helixAngleMap(gridPoint(2),gridPoint(1));
            new_value = helixAngles(pt,k);
            k = k + 1;
            if(isfinite(old_value) && isfinite(new_value))
                N = repeatCount(gridPoint(2),gridPoint(1));
                helixAngleMap(gridPoint(2),gridPoint(1)) = ((old_value*N)/(N+1)) + (new_value/(N+1));
                repeatCount(gridPoint(2),gridPoint(1)) = N+1;
            else
                if(isfinite(new_value))
                    helixAngleMap(gridPoint(2),gridPoint(1)) = new_value;
                    repeatCount(gridPoint(2),gridPoint(1)) = 1;
                end
            end
        end
    end
    % clean up stray nan pixels
    [I,J] = find(sectorMask>0 & isnan(helixAngleMap));
    for idx = length(I)%1:size(mask,1)
        i = I(idx);
        j = J(idx);
        assert(i < size(sectorMask,1));
        assert(j < size(sectorMask,2));
        searchRadius = 1;
        while(searchRadius < 7)
            rmin = max(1,i-searchRadius);
            rmax = min(i+searchRadius,size(mask,1));
            cmin = max(1,j-searchRadius);
            cmax = min(j+searchRadius,size(mask,2));
            region = reshape(helixAngleMap(rmin:rmax,cmin:cmax),1,[]);
            if(any(isfinite(region),'all'))
                helixAngleMap(i,j) = median(region,'omitnan');
                break;
            end
            searchRadius = searchRadius + 1;
        end
    end
    
    colored = abs(helixAngleMap);
    colored(~isfinite(colored)) = 0;
    colored = colored/(pi/2);
    colored = ind2rgb(uint8(colored*255),parula(256));
    for c = 1:3
        colored(:,:,c) = colored(:,:,c).*double(isfinite(helixAngleMap));
    end

end