%% tikz-nfold.sty
%% Copyright 2023 Jonathan Schulz
%
% This work may be distributed and/or modified under the
% conditions of the LaTeX Project Public License, either version 1.3c
% of this license or (at your option) any later version.
% The latest version of this license is in
% http://www.latex-project.org/lppl.txt
% and version 1.3c or later is part of all distributions of LaTeX
% version 2008-05-04 or later.
%
% This work has the LPPL maintenance status 'maintained'.
%
% The Current Maintainer of this work is Jonathan Schulz.
%
% This work consists of the files pgflibrarybezieroffset.code.tex,
% tikzlibrarynfold.code.tex, tikz-nfold-doc.tex, and tikz-nfold-doc.pdf.


% stores the current \pgf@x and \pgf@y in #1
\def\pgfstorepoint#1{\edef#1{\noexpand\pgfpoint{\the\pgf@x}{\the\pgf@y}}}
% global version in case we need it:
\def\pgfglobalstorepoint#1{\xdef#1{\noexpand\pgfpoint{\the\pgf@x}{\the\pgf@y}}}

% Split a Bezier curve (de Casteljau's algorithm)
% #1 = time (between 0 and 1)
% #2-#5: control points
% Outputs the first part into \pgf@splitbezier@i@i, ... , \pgf@splitbezier@i@iv,
% and the second part into \pgf@splitbezier@ii@i, ... , \pgf@splitbezier@ii@iv.

% TODO this code can be simplified to splitting at 0.5 only;
% this is partially implemented in some pgf file, possibly decorations or basic paths.
% But maybe I will need the general case in the future, maybe with some advanced fully simple detection.
% Leave it in for now
\newcommand{\pgf@splitbezier}[5]{
  % based on pgfcorepoints.code.tex, \pgfpointcurveattime
  \pgfmathparse{#1}%
  \let\pgf@time@s=\pgfmathresult%
  \global\pgf@x=\pgfmathresult pt%
  \global\pgf@x=-\pgf@x%
  \advance\pgf@x by 1pt%
  \edef\pgf@time@t{\pgf@sys@tonumber{\pgf@x}}%
  % P^0_3
  \pgf@process{#5}%
  \pgf@xc=\pgf@x%
  \pgf@yc=\pgf@y%
  \pgfstorepoint{\pgf@splitbezier@ii@iv}
  % P^0_2
  \pgf@process{#4}%
  \pgf@xb=\pgf@x%
  \pgf@yb=\pgf@y%
  % P^0_1
  \pgf@process{#3}%
  \pgf@xa=\pgf@x%
  \pgf@ya=\pgf@y%
  % P^0_0
  \pgf@process{#2}%
  \pgfstorepoint{\pgf@splitbezier@i@i}
  % First iteration:
  % P^1_0
  \global\pgf@x=\pgf@time@t\pgf@x\global\advance\pgf@x by\pgf@time@s\pgf@xa%
  \global\pgf@y=\pgf@time@t\pgf@y\global\advance\pgf@y by\pgf@time@s\pgf@ya%
  \pgfstorepoint{\pgf@splitbezier@i@ii}
  % P^1_1
  \pgf@xa=\pgf@time@t\pgf@xa\advance\pgf@xa by\pgf@time@s\pgf@xb%
  \pgf@ya=\pgf@time@t\pgf@ya\advance\pgf@ya by\pgf@time@s\pgf@yb%
  % P^1_2
  \pgf@xb=\pgf@time@t\pgf@xb\advance\pgf@xb by\pgf@time@s\pgf@xc%
  \pgf@yb=\pgf@time@t\pgf@yb\advance\pgf@yb by\pgf@time@s\pgf@yc%
  \edef\pgf@splitbezier@ii@iii{\noexpand\pgfpoint{\the\pgf@xb}{\the\pgf@yb}}
  % P^2_0
  \global\pgf@x=\pgf@time@t\pgf@x\global\advance\pgf@x by\pgf@time@s\pgf@xa%
  \global\pgf@y=\pgf@time@t\pgf@y\global\advance\pgf@y by\pgf@time@s\pgf@ya%
  \pgfstorepoint{\pgf@splitbezier@i@iii}
  % P^2_1
  \pgf@xa=\pgf@time@t\pgf@xa\advance\pgf@xa by\pgf@time@s\pgf@xb%
  \pgf@ya=\pgf@time@t\pgf@ya\advance\pgf@ya by\pgf@time@s\pgf@yb%
  \edef\pgf@splitbezier@ii@ii{\noexpand\pgfpoint{\the\pgf@xa}{\the\pgf@ya}}
  % P^3_0
  \global\pgf@x=\pgf@time@t\pgf@x\global\advance\pgf@x by\pgf@time@s\pgf@xa%
  \global\pgf@y=\pgf@time@t\pgf@y\global\advance\pgf@y by\pgf@time@s\pgf@ya%
  \pgfstorepoint{\pgf@splitbezier@i@iv}
  \pgfstorepoint{\pgf@splitbezier@ii@i}
}


% computes the cross product and puts it into \pgfmathresult
\newcommand{\pgfcrossproduct}[2]{
  \pgf@process{#1}%
  \pgf@xa=\pgf@x%
  \pgf@ya=\pgf@y%
  \pgf@process{#2}%
  \pgfmathparse{\pgf@xa*\pgf@y-\pgf@ya*\pgf@x}%
}

\newcommand{\pgfdotproduct}[2]{
  \pgf@process{#1}%
  \pgf@xa=\pgf@x%
  \pgf@ya=\pgf@y%
  \pgf@process{#2}%
  \pgfmathparse{\pgf@xa*\pgf@x+\pgf@ya*\pgf@y}%
}

\newcommand{\pgfcrossdot}[2]{
  \pgf@process{#1}%
  \pgf@xa=\pgf@x%
  \pgf@ya=\pgf@y%
  \pgf@process{#2}%
  \pgfmathsetlengthmacro{\pgf@tmp@dot}{\pgf@xa*\pgf@x+\pgf@ya*\pgf@y}%
  \pgfmathsetlengthmacro{\pgf@tmp@cross}{\pgf@xa*\pgf@y-\pgf@ya*\pgf@x}%
}


% Computes the normalised tangents of a given Bezier curve and stores them in \pgf@tmp@tang@i and \pgf@tmp@tang@ii.
% All degenerate cases are covered. For a triple degenerate curve (all points equal), the vector (1,0) is returned.
\newcommand{\pgf@offset@compute@tangents}[4]{
  \pgfpointdiff{#1}{#2}  % unintuitively, this is PTii - PTi
  \pgfmathparse{abs(\pgf@x) + abs(\pgf@y)}
  \ifdim\pgfmathresult pt<0.1pt\relax
    % edge case: first point and first control point are equal
    \pgfpointdiff{#1}{#3}
    \pgfmathparse{abs(\pgf@x) + abs(\pgf@y)}
    \ifdim\pgfmathresult pt<0.1pt\relax
      % edge case: first three points are equal
      \pgfpointdiff{#1}{#4}
    \fi
  \fi
  \pgfpointnormalised{}
  \pgfstorepoint\pgf@tmp@tang@i
  \pgfpointdiff{#3}{#4}
  \pgfmathparse{abs(\pgf@x) + abs(\pgf@y)}
  \ifdim\pgfmathresult pt<0.1pt\relax
    \pgfpointdiff{#2}{#4}
    \pgfmathparse{abs(\pgf@x) + abs(\pgf@y)}
    \ifdim\pgfmathresult pt<0.1pt\relax
      \pgfpointdiff{#1}{#4}
    \fi
  \fi
  \pgfpointnormalised{}
  \pgfstorepoint\pgf@tmp@tang@ii
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Offsetting a simple section %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\pgf@offset@bezier@segment}[5]{%
  % TODO would it make sense to use \pgf@process here?
  % normalise tangents and normals; this avoids overflow issues later, and we need
  % the normal vector to be of length 1 anyway
  \pgf@offset@compute@tangents{#1}{#2}{#3}{#4}
  % offset A1
  % compute the normal
  \pgf@tmp@tang@i
  \pgf@xa=\pgf@x
  \pgf@x=-\pgf@y
  \pgf@y=\pgf@xa
  \pgfstorepoint\pgf@tmp@normal@i
  % Leaving this parameter empty amounts to working directly on the register
  \pgfpointadd{\pgfpointscale{#5}{}}{#1}
  \pgfstorepoint{\pgf@bezier@offset@i}
  % offset A4
  \pgf@tmp@tang@ii
  \pgf@xa=\pgf@x
  \pgf@x=-\pgf@y
  \pgf@y=\pgf@xa
  \pgfstorepoint\pgf@tmp@normal@ii
  \pgfpointadd{\pgfpointscale{#5}{}}{#4}
  \pgfstorepoint\pgf@bezier@offset@iv
  % now compute A'_2 and A'_3
  \pgfpointdiff{#1}{#4}
  \pgfmathsetmacro{\pgf@tmp@secantlen}{veclen(\pgf@x,\pgf@y)}
  \ifdim\pgf@tmp@secantlen pt<0.1pt\relax
    % Edge case: Either the curve is degenerate to a point or it is not simple.
    % Either way we offset A1 and A4, and preserve the vectors A1A2 and A3A4.
    \pgfwarning{pgf-offset: first and last point are too close, expect glitches}
    \pgfpointadd{\pgf@bezier@offset@i}{\pgfpointdiff{#1}{#2}}
    \pgfstorepoint\pgf@bezier@offset@ii
    \pgfpointadd{\pgf@bezier@offset@iv}{\pgfpointdiff{#4}{#3}}
    \pgfstorepoint\pgf@bezier@offset@iii
  \else
    \pgfpointnormalised{}
    \pgfstorepoint\pgf@tmp@secant
    \pgfcrossdot{}{\pgf@tmp@tang@ii}
     \ifdim\pgf@tmp@dot<.5pt\relax%
       % this can only happen in non-simple curves
       \pgfwarning{pgf-offset: cosine of \pgf@tmp@dot\space clamped to 0.5 in non-simple segment}%
       \def\pgf@tmp@dot{.5pt}%
    \fi%
    \pgfmathsetmacro{\pgf@tmp@tanbeta}{\pgf@tmp@cross/\pgf@tmp@dot}%
    \pgfcrossdot{\pgf@tmp@secant}{\pgfpointnormalised{\pgfpointdiff{#1}{#2}}}
    % There are cases where we want #5/secantlen to be quite large, so we should not clamp the value here
    \pgfmathparse{1 + #5/\pgf@tmp@secantlen*(\pgf@tmp@cross - \pgf@tmp@dot*\pgf@tmp@tanbeta)}
    \pgfpointadd%
      {\pgf@bezier@offset@i}%
      {\pgfpointscale{\pgfmathresult pt}{\pgfpointdiff{#1}{#2}}}%
    \pgfstorepoint\pgf@bezier@offset@ii
    % third control point
    \pgfcrossdot{\pgf@tmp@secant}{\pgf@tmp@tang@i}
    \ifdim\pgf@tmp@dot<.5pt\relax%
      \pgfwarning{pgf-offset: cosine of \pgf@tmp@dot\space clamped to 0.5 in non-simple segment}%
      \def\pgf@tmp@dot{.5pt}%
    \fi%
    \pgfmathsetmacro{\pgf@tmp@tanbeta}{\pgf@tmp@cross/\pgf@tmp@dot}%
    \pgfcrossdot{\pgf@tmp@secant}{\pgfpointnormalised{\pgfpointdiff{#4}{#3}}}%
    \pgfmathparse{1 + #5/\pgf@tmp@secantlen*(\pgf@tmp@cross - \pgf@tmp@dot*\pgf@tmp@tanbeta)}%
    \pgfpointadd%
      {\pgf@bezier@offset@iv}%
      {\pgfpointscale{\pgfmathresult pt}{\pgfpointdiff{#4}{#3}}}%
    \pgfstorepoint\pgf@bezier@offset@iii%
  \fi
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subdividing and offsetting %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Maximum level of recursion. The theoretical limit to the number of subdivisions in the final curve
% is given by 2^\pgf@offset@max@recursion
\newcount\pgf@offset@max@recursion
\pgf@offset@max@recursion=4

%
% Subdivides a Bezier curve into "simple" segments (according to the definition below),
% offsets the segments, and draws them.  Because offsetting also involves relocating
% the starting points, these macros come in two variants: with and without a \pgfmoveto{}
% to the new starting point.
%
% Interface:
% #1-#4: control points of the whole Bezier curve
% #5: offset

\newcommand{\pgfoffsetcurve}[5]{%
  \pgf@subdivideandoffsetcurve{#1}{#2}{#3}{#4}{#5}{\pgf@offset@max@recursion}{0}{\pgf@nfold@callback@move}%
}
\newcommand{\pgfoffsetcurvenomove}[5]{%
  \pgf@subdivideandoffsetcurve{#1}{#2}{#3}{#4}{#5}{\pgf@offset@max@recursion}{0}{\pgf@nfold@callback@nomove}%
}

% Arguments:
% #1-#4: control points of the segment
% #5: =0 if this is the first segment of the curve, =1 otherwise
%   (checking for #5=0 allows us to draw the curve without interruptions)
\newcommand{\pgf@nfold@callback@move}[5]{%
  \ifnum#5=0\relax\pgfpathmoveto{#1}\fi%
  \pgfpathcurveto{#2}{#3}{#4}%
}
% this version never does a moveto at the start. Useful for drawing a path consisting of
% multiple Bezier curves.
\newcommand{\pgf@nfold@callback@nomove}[5]{\pgfpathcurveto{#2}{#3}{#4}}

% Like the previous macro, but with a custom callback macro for each segment instead of
% executing \drawsegment as defined above. See \drawsegment for the arguments.
\newcommand{\pgfoffsetcurvecallback}[6]{%
  \pgf@subdivideandoffsetcurve{#1}{#2}{#3}{#4}{#5}{\pgf@offset@max@recursion}{0}{#6}%
}


% 6 parameters:
% #1-#4: control points
% #5: offset
% #6: recursion limit (decreases on every recursive call)
% #7: =0 if this is the start of the curve, =1 otherwise;
% #8: callback for output (see above)
\newif\ifpgf@offset@subdivide
\newcommand{\pgf@subdivideandoffsetcurve}[8]{%
  % we need a group to avoid overwriting variables in recursive calls
  \begingroup%
  \pgf@offset@subdividefalse%
  \c@pgf@counta=#6%
  \advance\c@pgf@counta by-1%
  \pgf@process{#1}\pgfstorepoint{\pgf@ctrl@i}%
  \pgf@process{#2}\pgfstorepoint{\pgf@ctrl@ii}%
  \pgf@process{#3}\pgfstorepoint{\pgf@ctrl@iii}%
  \pgf@process{#4}\pgfstorepoint{\pgf@ctrl@iv}%
  % Use the non-degenerate tangents for the simplicity check
  \pgf@offset@compute@tangents{\pgf@ctrl@i}{\pgf@ctrl@ii}{\pgf@ctrl@iii}{\pgf@ctrl@iv}%
  \pgfpointdiff{\pgf@ctrl@i}{\pgf@ctrl@iv}\pgfstorepoint{\pgf@itoiv}%
  \pgfcrossproduct{\pgf@itoiv}{\pgf@tmp@tang@i}%
  \edef\firstcross{\pgfmathresult}%
  \pgfcrossproduct{\pgf@itoiv}{\pgf@tmp@tang@ii}%
  % First simplicity check: Are A2 and A3 on the same side of the A1-A4 line?
  % -> compute the sign of the cross products, use the sign function to avoid overflows
  \pgfmathparse{sign(\firstcross)*sign(\pgfmathresult)}%
  \ifdim\pgfmathresult pt>0pt\relax%
    \pgf@offset@subdividetrue%
  \else%
    % Second simplicity check: How large is the angle between the tangents in A1 and A4?
    \pgfdotproduct{\pgf@tmp@tang@i}{\pgf@tmp@tang@ii}%
    \ifdim\pgfmathresult pt<.5pt\relax%
      \pgf@offset@subdividetrue%
    \else
      % Third simplicity check: Put a limit on the lengths of the i-ii and iii-iv vectors combined
      \pgf@itoiv
      \pgfmathsetmacro{\pgf@tmp@len@i@iv}{veclen(\pgf@x,\pgf@y)}
      \pgfpointdiff{\pgf@ctrl@i}{\pgf@ctrl@ii}
      \pgf@xb=\pgf@x\pgf@yb=\pgf@y
      \pgfpointdiff{\pgf@ctrl@iii}{\pgf@ctrl@iv}
      \pgfmathparse{\pgf@tmp@len@i@iv < veclen(\pgf@xb,\pgf@yb) + veclen(\pgf@x,\pgf@y)}
      \ifnum\pgfmathresult=1\relax
        \pgf@offset@subdividetrue%
      \fi
    \fi%
  \fi%
  \ifpgf@offset@subdivide%
    \ifnum\c@pgf@counta<0%
      % We hit the recursion limit but the segment is not simple
      \pgfwarning{pgf-offset: Recursion limit reached, glitches may occur. %
        Consider increasing \string\pgf@offset@max@recursion}%
      % Try to offset the curve anyway. The result will not be precise,
      % but the code is sufficiently robust to not crash
      \pgf@offset@bezier@segment{\pgf@ctrl@i}{\pgf@ctrl@ii}{\pgf@ctrl@iii}{\pgf@ctrl@iv}{#5}%
      #8{\pgf@bezier@offset@i}{\pgf@bezier@offset@ii}{\pgf@bezier@offset@iii}{\pgf@bezier@offset@iv}{#7}%
    \else
      % split the non-simple segment and execute recursive calls
      \pgf@splitbezier{.5}{\pgf@ctrl@i}{\pgf@ctrl@ii}{\pgf@ctrl@iii}{\pgf@ctrl@iv}%
      % pass on the "start of the curve flag" only to the first term
      \pgf@subdivideandoffsetcurve{\pgf@splitbezier@i@i}{\pgf@splitbezier@i@ii}{\pgf@splitbezier@i@iii}{\pgf@splitbezier@i@iv}{#5}{\c@pgf@counta}{#7}{#8}%
      \pgf@subdivideandoffsetcurve{\pgf@splitbezier@ii@i}{\pgf@splitbezier@ii@ii}{\pgf@splitbezier@ii@iii}{\pgf@splitbezier@ii@iv}{#5}{\c@pgf@counta}{1}{#8}%
    \fi%
  \else%
    % curve is simple
    \pgf@offset@bezier@segment{\pgf@ctrl@i}{\pgf@ctrl@ii}{\pgf@ctrl@iii}{\pgf@ctrl@iv}{#5}%
    #8{\pgf@bezier@offset@i}{\pgf@bezier@offset@ii}{\pgf@bezier@offset@iii}{\pgf@bezier@offset@iv}{#7}%
  \fi%
  \endgroup%
}


%
% Offsetting straight lines
% -------------------------
%
% For convenience we also provide macros that offset straight lines. These also come in two variants
% similar to the macros for curves.
%
\newcommand{\pgfoffsetline}[3]{
  \pgfpointscale{#3}{\pgfpointnormalised{\pgfpointdiff{#1}{#2}}}
  \pgf@xc=-\pgf@y
  \pgf@yc=\pgf@x
  \pgfpathmoveto{\pgfpointadd{#1}{\pgfqpoint{\pgf@xc}{\pgf@yc}}}
  \pgfpathlineto{\pgfpointadd{#2}{\pgfqpoint{\pgf@xc}{\pgf@yc}}}
}
\newcommand{\pgfoffsetlinenomove}[3]{
  \pgfpointscale{#3}{\pgfpointnormalised{\pgfpointdiff{#1}{#2}}}
  \pgf@xc=-\pgf@y
  \pgf@yc=\pgf@x
  \pgfpathlineto{\pgfpointadd{#2}{\pgfqpoint{\pgf@xc}{\pgf@yc}}}
}
