\section{Matrix Calculus Transformations and Some Mathematical Functions}

\subsection{3D Transformations}

In the following functions:
\begin{itemize}
\item the argument \argu{L} is either a 3D point, a polyhedron, a list of 3D points (facet), or a list of lists of 3D points (facet list),
\item a line \argu{d} is a list of two 3D points $\{A,u\}$: a point on the line ($A$) and a direction vector ($u$),
\item a plane \argu{P} is a list of two 3D points $\{A,n\}$: a point on the plane ($A$) and a normal vector to the plane ($n$).
\end{itemize}
The returned result is of the same type as \argu{L}.

\subsubsection{Apply a transformation function: ftransform3d}

The function \cmd{ld.ftransform3d(L, f)} returns the image of \argu{L} by the function \argu{f}; this must be a function from $\mathbf R^3$ to $\mathbf R^3$ (to a 3D point it associates a 3D point).

\subsubsection{Projections: proj3d, proj3dO, dproj3d}

\begin{itemize}
\item The function \cmd{ld.proj3d(L, P)} returns the image of \argu{L} by the orthogonal projection onto the plane \argu{P}.
\item The function \cmd{ld.proj3dO(L, P, v)} returns the image of \argu{L} by the projection onto the plane \argu{P} parallel to the direction of the vector \argu{v} (non-zero 3D point).
\item The function \cmd{ld.dproj3d(L, d)} returns the image of \argu{L} by the projection onto the line \argu{d}.
\end{itemize}

\subsubsection{Projections onto axes or planes related to axes}

\begin{itemize}
\item The function \cmd{ld.pxy(L \fac{, z0})} returns the image of \argu{L} by the orthogonal projection onto the $z=$\argu{z0} plane (by default $z0=0$).
\item The function \cmd{ld.pyz(L \fac{, x0})} returns the image of \argu{L} by the orthogonal projection onto the $x=$\argu{x0} plane (by default $x0=0$).
\item The function \cmd{ld.pxz(L \fac{, y0})} returns the image of \argu{L} by the orthogonal projection onto the $y=$\argu{y0} plane (by default $y0=0$).
\item The function \cmd{ld.px(L)} returns the image of \argu{L} by the orthogonal projection onto the $x$-axis.
\item The function \cmd{ld.py(L)} returns the image of \argu{L} by the orthogonal projection onto the $y$-axis.
\item The function \cmd{ld.pz(L)} returns the image of \argu{L} by the orthogonal projection onto the $z$-axis.
\end{itemize}

\subsubsection{Symmetries: sym3d, sym3dO, dsym3d, psym3d}

\begin{itemize}
\item The function \cmd{ld.sym3d(L, P)} returns the image of \argu{L} by the orthogonal symmetry about the \argu{P} plane.
\item The function \cmd{ld.sym3dO(L, P, v)} returns the image of \argu{L} by the symmetry about the \argu{P} plane and parallel to the direction of the \argu{v} vector (non-zero 3D point).
\item The function \cmd{ld.dsym3d(L, d)} returns the image of \argu{L} by the orthogonal symmetry with respect to the line \argu{d}.
\item The function \cmd{ld.psym3d(L, point)} returns the image of \argu{L} by the symmetry with respect to \argu{point} (3d point).
\end{itemize}

\subsubsection{Rotation: rotate3d, rotateaxe3d}

\begin{itemize}
\item The function \cmd{ld.rotate3d(L, angle, d)} returns the image of \argu{L} rotated along axis \argu{d} (oriented by the direction vector, which is $d[2]$), and by \emph{angle} degrees.

\item The function \cmd{ld.rotateaxe3d(L, v1, v2 \fac{, center})} returns the image of \argu{L} rotated along axis passing through the 3D point \argu{center}, which transforms the vector \argu{v1} into the vector \argu{v2}. These vectors are normalized by the function. The argument \argu{center} is optional and defaults to the point \varglob{pt3d.Origin}.
\end{itemize}

\subsubsection{Scaling: scale3d}

The function \cmd{ld.scale3d(L , k \fac{, center})} returns the image of \argu{L} by the scaling with center at the 3D point \argu{center}, and ratio \argu{k}. The argument \argu{center} is optional and defaults to the point \varglob{pt3d.Origin}.

\subsubsection{Inversion: inv3d}

The function \cmd{ld.inv3d(L, radius \fac{, center})} returns the image of \argu{L} by the inversion with respect to the sphere with center \argu{center}, and radius \argu{radius}. The argument \argu{center} is optional and defaults to the point \varglob{pt3d.Origin}.

\subsubsection{Stereography: projstereo and inv\_projstereo}

Function \cmd{ld.projstereo(L, S, N, h)}: the argument \argu{L} denotes a 3D point or a list of 3D points or a list of lists of 3D points, all belonging to the sphere \argu{S}, where \argu{S}$=\{C,r\}$ ($C$ is the center of the sphere, and $r$ the radius). The argument \argu{N} denotes a point on the sphere that will be the pole of the projection. The argument \argu{h} is a real number that defines the projection plane. This plane is perpendicular to the axis $(CN)$, and passes through the point $I=C+h \frac{\vec{CN}}{CN}$ (with $h=0$ it is the equatorial plane, with $h=-r$ it is the plane tangent to the sphere at the opposite pole). The function returns the image of \argu{L} by the stereographic projection with respect to the sphere \argu{S} with \argu{N} as pole, and on the plane $\{I,N-C\}$.

Inverse function \cmd{ld.inv\_projstereo(L, S, N)}: \argu{S}$=\{C,r\}$ is the sphere with center $C$ and radius $r$, \argu{N} is a point on the sphere \argu{S} (pole), and \argu{L} is a 3D point or a list of 3D points or a list of lists of 3D points all belonging to the same plane orthogonal to the $(CN)$ axis. The function returns the image of \argu{L} by the inverse of the stereographic projection with respect to \argu{S} and with pole \argu{N}.

\subsubsection{Translation: shift3d}

The function \cmd{ld.shift3d(L, v)} returns the image of \argu{L} by the translation of vector \argu{v} (3D point).


\subsection{Matrix Calculus}

If $f$ is an affine mapping of the space $\mathbf R^3$, we will call the list (table) of $f$ a matrix:
\codeln{ \{f(pt3d.Origin),Lf(pt3d.vecI),Lf(pt3d.vecJ),Lf(pt3d.vecK)\} }
where $Lf$ denotes the linear part of $f$ (we have \code{Lf(pt3d.vecI)=f(pt3d.vecI)-f(pt3d.Origin)}, etc.). The identity matrix is ​​denoted \varglob{ld.ID3d} in \luadrawenv; it simply corresponds to the list: \codeln{\{pt3d.Origin,pt3d.vecI,pt3d.vecJ,pt3d.vecK\}}

\subsubsection{applymatrix3d ​​and applyLmatrix3d}

\begin{itemize}
\item The function \cmd{ld.applymatrix3d(A, M)} applies the matrix \argu{M} to the 3D point \argu{A} and returns the result (which is equivalent to calculating $f(A)$ if $M$ is the matrix of $f$). If \argu{A} is not a 3D point, the function returns \argu{A}.

\item The function \cmd{ld.applyLmatrix3d(A, M)} applies the linear part of the matrix \argu{M} to the 3D point \argu{A} and returns the result (which is equivalent to calculating $Lf(A)$ if $M$ is the matrix of $f$). If \argu{A} is not a 3D point, the function returns \argu{A}.
\end{itemize}

\subsubsection{composematrix3d}
The function \cmd{ld.composematrix3d(M1, M2)} performs the matrix product \argu{M1}$\times$\argu{M2} and returns the result.

\subsubsection{invmatrix3d}
The function \cmd{ld.invmatrix3d(M)} calculates and returns the inverse of the matrix \argu{M} when possible.

\subsubsection{matrix3dof}

The function \cmd{ld.matrix3dof(f)} calculates and returns the matrix of \argu{f} which must be an affine mapping of the space $\mathbf R^3$ (to a 3D point it associates a 3D point).

\subsubsection{mtransform3d and mLtransform3d}

\begin{itemize}
\item The function \cmd{ld.mtransform3d(L, M)} applies the matrix \argu{M} to the list \argu{L} and returns the result. \argu{L} must be a list of 3D points (a facet) or a list of lists of 3D points (a list of facets).
\item The function \cmd{ld.mLtransform3d(L, M)} applies the linear part of the matrix \argu{M} to the list \argu{L} and returns the result. \argu{L} must be a list of 3D points (a facet) or a list of lists of 3D points (a list of facets).
\end{itemize}

\subsection{Matrix associated with the 3D graph}

When creating a graph in the \luadrawenv environment, for example:
\begin{Luacode}
local ld = luadraw
local g = ld.graph3d:new{size={10,10}}
\end{Luacode}

The created \emph{g} object has a 3D transformation matrix that is initially the identity. All graphics methods automatically apply the graph's 3D transformation matrix. To manipulate this matrix, the following methods are available.

\subsubsection{g:Composematrix3d()}

The \cmd{g:Composematrix3d(M)} method multiplies the 3D matrix of the graph \emph g by the matrix \argu{M} (with \argu{M} on the right), and the result is assigned to the graph's 3D matrix. The argument \argu{M} must therefore be a 3D matrix.

\subsubsection{g:Det3d()}

The \cmd{g:Det3d()} method returns $1$ when the 3D transformation matrix has a positive determinant, and $-1$ otherwise. This information is useful when we need to know whether the orientation of space has been changed or not.

\subsubsection{g:IDmatrix3d()}

The \cmd{g:IDmatrix3d()} method reassigns the identity to the 3D matrix of the graph \emph g.

\subsubsection{g:Mtransform3d()}

The \textbf{g:Mtransform3d(L)} method applies the 3D graph matrix of \emph g to \emph{L} and returns the result. The argument \emph L must be a list of 3D points (a facet) or a list of lists of 3D points (a list of facets).

\subsubsection{g:MLtransform3d()}

The \cmd{g:MLtransform3d(L)} method applies the linear part of the 3D matrix of the graph \emph g to \argu{L} and returns the result. The argument \argu{L} must be a list of 3D points (a facet) or a list of lists of 3D points (a list of facets).

\subsubsection{g:Rotate3d()}

The method \cmd{g:Rotate3d(angle, axis)} modifies the 3D transformation matrix of the graph \emph g by composing it with the rotation matrix of angle \argu{angle} (in degrees) and axis \argu{axis}.

\subsubsection{g:Scale3d()}

The method \cmd{g:Scale3d(factor \fac{, center})} modifies the 3D transformation matrix of the graph \emph g by composing it with the homothety matrix of ratio \argu{factor} and center \argu{center}. The argument \argu{center} is a 3D point that defaults to \varglob{pt3d.Origin}.

\subsubsection{g:Setmatrix3d()}

The \cmd{g:Setmatrix3d(M)} method assigns the matrix \argu M to the 3D transformation matrix of the graph \emph g.

\subsubsection{g:Shift3d()}

The \cmd{g:Shift3d(v)} method modifies the 3D transformation matrix of the graph \emph g by composing it with the translation matrix of vector \argu{v}, which must be a 3D point.


\subsection{Additional Mathematical Functions}

\subsubsection{clippolyline3d()}

The function \cmd{ld.clippolyline3d(L, poly \fac{, exterior, close})} clips the 3D polygonal line \argu{L} to the \textbf{convex} polyhedron \argu{poly}. If the optional argument \argu{exterior} is \true, then the part outside the polyhedron is returned (\false by default). If the optional argument \argu{close} is \true, then the polygonal line is closed (\false by default). \argu{L} is a list of 3D points or a list of lists of 3D points.

\textbf{Note}: The result is not always satisfactory for the exterior part.

\paragraph{Special case}: Clipping a 3D polygonal line \argu L with the current 3D window can be done with this function as follows:
\codeln{L=ld.clippolyline3d(L,g:Box3d())}

Indeed, the \cmd{g:Box3d()} method returns the current 3D window as a parallelepiped.

\subsubsection{clipline3d()}

The function \cmd{ld.clipline3d(line, poly)} clips the line \argu{line} with the \textbf{convex} polyhedron \argu{poly}; the function returns the part of the line inside the polyhedron. The argument \argu{line} is a table of the form $\{A,u\}$ where $A$ is a point on the line and $u$ is a direction vector (two 3D points).

\paragraph{Special case}: Clipping a line \argu d with the current 3D window can be done with this function as follows:
\codeln{d=clipline3d(d,g:Box3d())}

Indeed, the \cmd{g:Box3d()} method returns the current 3D window as a parallelepiped (\argu d then becomes a segment).

\subsubsection{cutpolyline3d()}

The function \cmd{ld.cutpolyline3d(L, plane \fac{, close})} intersects the 3D polygonal line \argu{L} with the plane \argu{plane}. If the optional argument \argu{close} is \true, then the line is closed (\false by default). \argu{L} is a list of 3D points or a list of lists of 3D points, \argu{plane} is a table of the form $\{A,n\}$ where $A$ is a point in the plane and $n$ is a normal vector (two 3D points).

The function returns three things:
\begin{itemize}
\item the part of \argu{L} that is in the half-space containing the vector $n$,
\item followed by the part of \argu{L} that is in the other half-space,
\item followed by the list of intersection points.
\end{itemize}

\subsubsection{getbounds3d()}

The function \cmd{ld.getbounds3d(L)} returns the bounds \emph{xmin, xmax, ymin, ymax, zmin, zmax} (sequence) of the 3D polygonal line \argu{L} (list of 3D points or a list of lists of 3D points).

\subsubsection{interDP()}

The function \cmd{ld.interDP(d, P)} calculates and returns (if it exists) the intersection between line \argu d and plane \argu P.

\subsubsection{interPP()}

The function \cmd{ld.interPP(P1, P2)} calculates and returns (if it exists) the intersection between planes \argu{P1} and \argu{P2}.

\subsubsection{interDD()}

The function \cmd{ld.interDD(D1, D2 \fac{, epsilon})} calculates and returns (if it exists) the intersection between lines \argu{D1} and \argu{D2}. The argument \argu{epsilon} is $10^{-10}$ by default (used to test whether a given float is zero).

\subsubsection{interCS()}

The function \cmd{ld.interCS(C, S)} calculates and returns (if it exists) the intersection between the circle \argu{C}$=\{A,r,n\}$ ($A$ is the center of the circle, $r$ the radius, and $n$ a normal vector to the plane of the circle), and the sphere \argu{S}$=\{B,R\}$ ($B$ is the center of the sphere and $R$ the radius). The function returns either \nil (empty intersection), a single point, or two points (sequence).

\subsubsection{interDS()}

The function \cmd{ld.interDS(d, S)} calculates and returns (if it exists) the intersection between the line \argu{d} and the sphere \argu{S} which is a table of the form $\{C,r\}$ with $C$ as the center (3D point) and $r$ as the radius of the sphere. The function returns either \nil (empty intersection), a single point, or two points.

\subsubsection{interPS()}

The function \cmd{ld.interPS(P, S)} calculates and returns (if it exists) the intersection between the plane \argu{P} and the sphere \argu{S}$=\{C,r\}$ with $C$ as the center (3d point) and $r$ as the radius of the sphere. The function returns either \nil (empty intersection) or a sequence of the form $I,r,n$, where $I$ is a 3D point representing the center of a circle, $r$ its radius, and $n$ a normal vector to the plane of the circle. This circle is the desired intersection.

\subsubsection{interSS()}

The function \cmd{ld.interSS(S1, S2)} calculates and returns (if it exists) the intersection between the sphere \argu{S1}$=\{C1,r1\}$ and \argu{S2}$=\{C2,r2\}$. The function returns either \nil (empty intersection) or a sequence of the form $I,r,n$, where $I$ is a 3D point representing the center of a circle, $r$ its radius, and $n$ a normal vector to the plane of the circle. This circle is the desired intersection.

\subsubsection{interSSS()}

The function \cmd{ld.interSSS(S1, S2, S3)} calculates and returns (if it exists) the intersection between the spheres \argu{S1}$=\{C1,r1\}$, \argu{S2}$=\{C2,r2\}$ and \argu{S3}$=\{C3,r3\}$. The function returns either \nil (empty intersection), a single point, or two points (sequence).

\subsubsection{merge3d()}

The function \cmd{ld.merge3d(L \fac{, epsilon})} combines, if possible, the connected components of \argu{L}, which must be a list of lists of 3D points. The function returns the result. The argument \argu{epsilon} defaults to $10^{-10}$, it is used during comparisons.

\subsubsection{split\_points\_by\_visibility()}

The function \cmd{ld.split\_points\_by\_visibility(L, visible\_function)}, where \argu L is a list of 3D points, or a list of lists of 3D points, and where \argu{visible\_function} is a function such that \code{visible\_function(A)} returns \true if the 3D point $A$ is visible, \false otherwise, sorts the points of \argu L according to whether they are visible or not. The function returns a sequence of two tables: \emph{visible\_points}, \emph{hidden\_points}.

\begin{demo}{A curve on a cylinder}
\begin{luadraw}{name=curve_on_cylinder}
local ld = luadraw
local pt3d = ld.pt3d
local Origin, vecI, vecJ, vecK, M, Mc = pt3d.Origin, pt3d.vecI, pt3d.vecJ, pt3d.vecK, pt3d.M, pt3d.Mc

local g = ld.graph3d:new{adjust2d=true,bbox=false,size={10,10}, viewdir="central"}
g:Labelsize("footnotesize")
ld.Hiddenlines = true; ld.Hiddenlinestyle = "dashed"
local curve_on_cylinder = function(curve,cylinder)
-- curve is a 3D polyline on a cylinder,
-- cylinder = {A,r,V,B}
    local A,r,V,B = table.unpack(cylinder)
    if B == nil then B = V; V = B-A end
    local U = B-A
    local visible_function
    if ld.projection_mode == "central" then
        visible_function = function(N)
            local I = ld.dproj3d(N,{A,U})
            local M1, M2 = ld.interCS({I,r,U},{ (I+ld.camera)/2, pt3d.abs(I-ld.camera)/2})
            local alpha = pt3d.angle3d(M1-I,ld.camera-I)
            return pt3d.angle3d(N-I,ld.camera-I) <= alpha
        end
    else
        visible_function = function(N)
            local I = ld.dproj3d(N,{A,U})
            return (pt3d.dot(N-I,g.Normal) >= 0)
        end
    end
    return ld.split_points_by_visibility(curve,visible_function)
end
-- test
local A, r, B = -5*vecJ, 4, 5*vecJ -- cylinder
local p = function(t) return Mc(r,t,t/3) end
local Curve = ld.rotate3d( ld.parametric3d(p,-4*math.pi,4*math.pi),90,{Origin,vecI})
local Vi, Hi = curve_on_cylinder(Curve,{A,r,B})
local curve_color = "DarkGreen"
g:Dboxaxes3d({grid=true,gridcolor="gray",fillcolor="LightGray"})
g:Dcylinder(A,r,B,{color="orange"})
g:Dpolyline3d(Vi,curve_color)
g:Dpolyline3d(Hi,curve_color..","..ld.Hiddenlinestyle)
g:Show()
\end{luadraw}
\end{demo}