function ptsi = NE_inversion_all(P,Q,circle,pts)
% NE_inversion_all - invert a set of points through a circle
%                  2 of the pts should be on the circle
% On input:
%     P (1x2 vector): point on circle (used when straight)
%     Q (1x2 vector): point on circle
%     circle (1x3 vector): center, radius
%     pts (nx2 array): x,y locations in the Poincare disk
% On output:
%     ptsi (nx2 array): inversion points
% Call:
%     ptsi = NE_inversion_all(cir1,pts);
% Author:
%     T. Henderson
%     UU
%     Summer 2023
%

[num_pts,~] = size(pts);
ptsi = pts;
C = circle(1:2)';
radius = circle(3);

for p = 1:num_pts
    X = pts(p,:)';
    CX = norm(C-X);
    if isinf(radius) % reflect cross straight line
        v = Q - P;
        v = v/norm(v);
        n = [-v(2),v(1)];
        line1 = [P,0;Q,0];
        X1 = X' + n;
        line2 = [X',0;X1,0];
        [ip,status] = cv_int_line_line(line1,line2,0.001);
        d = norm(ip(1:2)-X');
        w = ip(1:2) - X';
        if norm(w)>0
            w = w/norm(w);
            ptsi(p,:) = X' + 2*d*w;
        end
    elseif CX~=radius  % invert through curve
        % |CX|*|CP| = radius^2
        % ==> |CP| = (radius^2)/|CX|
        % so P is |CP| away from C in the X-C direction
        P = C + ((X-C)/CX)*((radius^2)/CX);
        ptsi(p,:) = P';
    end
end
