function circleplot(data,mask,img)
    % Form figure
    figure(1); 
        
    % Buffer for plot points
    points_plot = [];
    
    % Do image subplot first
    subplot(1,2,1);
    imshow(im2uint8(cat(3,img,img,img)),[]);  
    hold('on');
    handle_data = imshow(data,[]);
        
    m = size(get(gcf,'colormap'),1);
    n = ceil(m/4);
    u = [(1:1:n)/n ones(1,n-1) (n:-1:1)/n]';
    g = ceil(n/2) - (mod(m,4)==1) + (1:length(u))';
    r = g + n;
    b = g - n;
    g(g>m) = [];
    r(r>m) = [];
    b(b<1) = [];
    J = zeros(m,3);
    J(r,1) = u(1:length(r));
    J(g,2) = u(1:length(g));
    J(b,3) = u(end-length(b)+1:end);
    set(gcf, 'Colormap', J); 
    
    alphamap = mask.*0.7;
    set(handle_data,'AlphaData', alphamap);
    caxis(gca,[prctile(data(mask),5) prctile(data(mask),95)]);
    
    % Set ellipse
    handle_ellipse = imellipse(gca);          
    ellipse_api = iptgetapi(handle_ellipse);
    pos_ellipse = getPosition(handle_ellipse);
    setPosition(handle_ellipse,[pos_ellipse(1) pos_ellipse(2) max(pos_ellipse(3:4)) max(pos_ellipse(3:4))]);
    ellipse_api.setFixedAspectRatioMode(1);
    
    addNewPositionCallback(handle_ellipse,@(pos)callback_ellipse(pos));
    
    % Set start point
    pos_ellipse = getPosition(handle_ellipse);
    
    handle_point_start = impoint(gca,[pos_ellipse(1)+pos_ellipse(3) pos_ellipse(2)+pos_ellipse(3)/2]);
    setColor(handle_point_start,'g');
    point_api_start = iptgetapi(handle_point_start);
    point_api_start.setPositionConstraintFcn(@constraint_point);
    
    handle_point_end = impoint(gca,[pos_ellipse(1) pos_ellipse(2)+pos_ellipse(3)/2]);
    setColor(handle_point_end,'r');
    point_api_end = iptgetapi(handle_point_end);
    point_api_end.setPositionConstraintFcn(@constraint_point);
    
    % Call it to initialize
    callback_ellipse(getPosition(handle_ellipse));
    
    function callback_ellipse(pos)  
        % Update positions of points
        setPosition(handle_point_start,constraint_point(getPosition(handle_point_start)));
        setPosition(handle_point_end,constraint_point(getPosition(handle_point_end)));
        
        % Update plot
        update_plot();
    end
    
    function pos_closest = constraint_point(pos) 
        % Get position of ellipse
        pos_ellipse = getPosition(handle_ellipse);
        radius_ellipse = pos_ellipse(3)/2;
        
        % Find closest point to ellipse
        pos_center = [pos_ellipse(1)+radius_ellipse pos_ellipse(2)+radius_ellipse];
        
        % Find vector from center to the location of the point
        disp = pos-pos_center;
        
        % Scale vector to be length radius then add center 
        pos_closest = disp./(sqrt(disp(1)^2+disp(2)^2))*radius_ellipse + pos_center;             
        
        % Update plot
        update_plot();
    end

    function update_plot()
        pos_ellipse = getPosition(handle_ellipse);
        radius_ellipse = pos_ellipse(3)/2;
        pos_center = [pos_ellipse(1)+radius_ellipse pos_ellipse(2)+radius_ellipse];
        
        vec_start = getPosition(handle_point_start) - pos_center;
        vec_end = getPosition(handle_point_end) - pos_center;
        
        % Find angle between vec_start and vec_end
        angle_interior = acos((vec_start(1)*vec_end(1) + vec_start(2)*vec_end(2))/(sqrt(vec_start(1)^2+vec_start(2)^2)*sqrt(vec_end(1)^2+vec_end(2)^2)));
          
        % If sign of cross product between vec_start and vec_end is
        % positive, then angle is 2pi - angle interior
        if (sign(vec_start(1)*vec_end(2)-vec_start(2)*vec_end(1)) == 1)
            angle = 2*pi - angle_interior;
        else
            angle = angle_interior;
        end
        
        N = 100;        
        if (~isempty(points_plot))
            delete(points_plot);
            points_plot = [];
        end
              
        vec_interp = zeros(0,2);
        for i = 0:100
            angle_increment = angle*i/N;
            
            % Rotate start vector by angle_increment
            vec_interp(i+1,:) = vec_start*[cos(angle_increment) -sin(angle_increment); sin(angle_increment) cos(angle_increment)];    
        end                 
        % Get handle so they can be deleted
        subplot(1,2,1);
        hold on;        
        points_plot = plot(vec_interp(:,1)*1.1+pos_center(1),vec_interp(:,2)*1.1+pos_center(2),'yo','MarkerSize',4);        
        hold off;
                
        % Add center to vec_interp
        vec_interp(:,1) = vec_interp(:,1)+pos_center(1);
        vec_interp(:,2) = vec_interp(:,2)+pos_center(2);
        
        val_interp = interp2(data,vec_interp(:,1),vec_interp(:,2));

        subplot(1,2,2);
        plot(0:N,val_interp,'ro');
    end
end

