Chaotiﬁcation of 1D Maps by Multiple Remainder Operator Additions—Application to B-Spline Curve Encryption

: In this work, a chaotiﬁcation technique is proposed for increasing the complexity of chaotic maps. The technique consists of adding the remainder of multiple scalings of the map’s value for the next iteration, so that the most- and least-signiﬁcant digits are combined. By appropriate parameter tuning, the resulting map can achieve a higher Lyapunov exponent value, a result that was ﬁrst proven theoretically and then showcased through numerical simulations for a collection of chaotic maps. As a proposed application of the transformed maps, the encryption of B-spline curves and patches was considered. The symmetric encryption consisted of two steps: a shufﬂing of the control point coordinates and an additive modulation. A transformed chaotic map was utilised to perform both steps. The resulting ciphertext curves and patches were visually unrecognisable compared to the plaintext ones and performed well on several statistical tests. The proposed work gives an insight into the potential of the remainder operator for chaotiﬁcation, as well as the chaos-based encryption of curves and computer graphics.


Applications of Chaos and Chaotification
The applications of chaos in the digital world have been reported in multiple works [1][2][3][4] and span the areas of computing, encryption, surveillance, communications, fingerprinting, and many more. Within the broad range of these applications, both continuous and discrete chaotic systems are considered, depending on the application at hand.
Continuous chaotic systems have many merits, such as more degrees of freedom and a higher state space. However, this comes at the cost of computational time, which may pose problems in applications where execution time is of the essence. They are also heavily dependent on the numerical method used to compute their solution, which creates issues of reproducibility across different computers and digital devices. This is why it is more suitable to consider discrete-time chaotic systems, or chaotic maps as they are also termed. Discrete maps can achieve chaotic behaviour even in a single dimension, have a much lighter computational cost, and are solved iteratively, so they are not dependent on numerical integration methods.
Still, for both continuous and discrete systems, there are several features that need to be present, in order to make a chaotic system suitable for use in security-related applications, such as encryption, watermarking, surveillance, and others. Desirable features include the absence of parameter ranges that lead to periodic behaviour and an increased sensitivity to the initial conditions, indicated by high Lyapunov exponent values, a high key space, and an unrecognisable phase diagram with no visual structure. These are all attributes of maps with increased complexity. The pursuit of finding such maps led to the field of study known as chaotification, which in general refers to the problem of constructing novel chaotic maps or devising methods to modify known maps in order to satisfy the above properties. Since there is already an abundance of chaotic maps in the literature, chaotification usually does not involve the generation of novel maps, but techniques that can be applied to any known map to improve its performance in cryptographic and similar tasks.
In [18], the authors proposed a chaotification technique that consisted of adding a remainder term to any 1D map, in order to mix the most-and least-significant digits in each iteration and boost the chaotic behaviour of the seed map. Indeed, it was shown that, by the appropriate parameter choice, the Lyapunov exponent of the modified map can be larger than that of the seed map. Moreover, the modified map showed wide areas of uninterrupted chaotic behaviour in its parameter space. The present work is a straightforward extension of the chaotification technique in [18] by using a summation of remainder operators. The transformed maps are then considered for the encryption of B-spline curves and patches.

Encryption
The principles of chaos-based encryption are applicable to any type of data and, so far, have indeed been applied to the encryption of images [11,12,24,25], text [25], and continuous signals [26]. Image encryption seems to be the most-popular application, mainly due to its visual aspect, which makes it easier to depict and compare plaintext and ciphertext data and interpret the results of the various statistical tests performed on them.
In addition to the above, recently, there has been an increasing interest in encrypting 3D objects. Three-dimensional graphics can represent any type of data, from mechanical and biological objects, to video game graphics, [27][28][29]. The format used to represent such graphics is the stereolithography (STL) format. In this work, in order to complement the work on computer graphics encryption, B-spline curves and patches will be chosen as the plaintext data.
B-spline curves, and their simpler form, Bézier curves, are a well-established family of parametric curves, which are extensively used for modelling 2D and 3D graphics [30][31][32][33][34][35][36]. Their use originated in the automotive industry, but currently, their usage has been expanded to a plethora of applications, most notably graphic design, robotics, and navigation, as they can be used to effectively generate motion trajectories for a robotic manipulator, a UVG, or a UAV.
There are several properties that make B-spline curves so useful. First is the ease of their implementation, as they are simply described by a closed-form formula. That formula requires only the specification of a set of control points characterising the shape of the curve. These points also define the convex hull of the curve, so it is easy to specify its boundaries, making it easier to avoid intersections with obstacles. Another property is that they are locally controllable, meaning that changing a single control point of the B-spline will affect only a part of the curve. Finally, their smoothness can also be controlled by choosing the appropriate curve order and control point placement.
In this work, a symmetric encryption technique of B-spline curves and patches was considered, which consists of two separate steps. The first step involves the shuffling of the control points, by the application of a chaotic pseudorandom rule. The second consists of combining the control points with the values of a transformed chaotic map, to mask their value. The resulting ciphertext curves were visually unrecognisable and also performed well on several statistical tests. This technique complements previous works that considered the encryption of 3D objects and proposes a setup for encrypting parametric curves.
The main contributions of this work is the construction and analysis of mappings using multiple remainder operators and the presentation of a new B-spline encryption algorithm. The use of multiple remainder operators can efficiently mix the least-and most-significant digits between iterative mappings and can be further explored and modified in the future. The encryption of B-spline curves using a two-step process can be a starting point for further experimentation into encrypting geometrical data structures in general.
The rest of the work is structured as follows. In Section 2, the proposed chaotification technique is introduced. In Section 3, the technique is applied to several known maps. In Section 4, the effect of the proposed chaotification technique in expanding the key space of the resulting encryption protocol is discussed. In Section 5, the computational cost of applying the technique is presented. In Sections 6 and 7, the symmetric encryption of B-spline curves and surfaces is considered. The simulation results of encryption are shown in Section 8. Finally, Section 9 concludes the present work and discusses topics for future research.

The Proposed Chaotification Technique
Consider a one-dimensional discrete time map described by the iterative formula: The proposed chaotification technique is an extension of [18] and consists of using the remainder operator to add multiple scalings of the map's value x i−1 into the next iteration x i as follows: where a j > 0 are the control parameters, N j ∈ R + * map the remainder result on the interval (−N j , N j ), and m ∈ N. The technique is also illustrated in Figure 1.

Theorem 1. Let
Then, the following hold: where D m is the domain of F .

Proof.
We begin by considering the derivative of the remainder operation, which is for almost all x. Using (6), the derivative of H m (x i ) as in (3) for almost all Observe that ∑ m j=1 a j is independent of x i , and thus, if we denote it as which is exactly the case studied in [18], from which the conclusion is that, selecting results in a system with a greater Lyapunov exponent. Hence, in our case, the terms a j must be chosen so that m ∑ j=1 a j > 2 sup x i ∈D m Ḟ (x i ) . We can extend the proof in that work by observing that a selection of A m < −2 sup yields the same result. Thus, combining (10) and (11) yields that, selecting m ∑ j=1 a j > 2 sup results in a system that has a higher Lyapunov exponent than that of the original system. Observe that, until this point, no specific assumption has been made for each of the coefficients a j individually, only for their sum.
Next, we are interested in determining conditions such that the sequence of Lyapunov exponents created by adding extra terms in (3) is increasing. First of all, we observe that for every m ∈ N, the domains of F and H m coincide. With that in mind, we can then write = H m (x i ) + rem(a m+1 x i , N m+1 ).
Combining (15) with the first part of the proof and the fact that the domains of F and H coincide yields that selecting a m+1 as |a m+1 | > 2 sup (16) guarantees that λ m+1 > λ m for any m, which concludes the proof.

Remark 1.
Note that the choice of the remainder operator is made over the more common modulo operator, as the remainder operator does not change the sign of its input. This results in a more uniform distribution of the map's values, compared to the use of the modulo. Figure 2 illustrates this difference using a modified version of the Sine Map. In addition, it is worth noting that the remainder operator results in a symmetric distribution in this case. In contrast, the map distribution with the modulo operator does not have such symmetry.

Chaotification Examples for Different Chaotic Maps
In this section, we apply the above chaotification technique to several well-known maps as shown below.

Sine Map
Consider the well-known Sine Map [37,38]: Using the Sine Map (17) as the base function F (x), two cases of mappings of the form (2) are further analysed, for m = 2 and m = 3: The first map (18), apart from the base function (17), consists of two additional components: rem(10x i−1 , 1) and rem(100x i−1 , 1). In turn, the second map (19) modifies the first with an additional component: rem(1000x i−1 , 1). One of the tools for analysing dynamical systems is to present the so-called phase diagrams that show the system's relationship (x i−1 , x i ). The phase diagrams of all three cases (Sine Map (17) and its two modifications (18) and (19)) for fixed parameter values are presented in Figure 3a. These graphs show that adding multiple rem(·) operators causes a significant increase in the complexity of the base mapping. In turn, Figure 3b presents bifurcation diagrams of the analysed mappings. This figure depicts the existence of the so-called periodic windows for the Sine Map (17) and their respective absence for the modified mappings (18) and (19). This means that, in the case of mappings (18) and (19), there is no worry about choosing the proper mapping parameters, as any choice will always result in chaotic behaviour.
The increased complexity of the resulting maps is also confirmed by the three-dimensional phase diagrams, which show the relation of points (x i , x i−1 , x i−2 ). These charts are presented in Figure 4. In the case of the Sine Map (17), the graph in Figure 4a shows that points fall on a curve, while for the modified maps (19) (Figure 4b), they fill up a 3D region, forming a cloud of points. Therefore, the mapping with the rem(·) operators is much more complex than the base mapping.  Chaotic systems are characterised by the system's sensitivity to a change in initial conditions, which is measured using the Lyapunov exponent λ given by the formula: The plot of the Lyapunov exponent is shown in Figure 5a. This figure clearly shows the presence of stable/periodic regions for the Sine Map (17) and their absence for mappings with the rem(·) operator. This figure also shows that the sensitivity to the change of initial conditions for the mappings (18) and (19) is roughly constant for the entire range of the parameter k. Thus, after the same number of iterations, the system will lose information about the initial conditions. Moreover, to better illustrate the chaotic nature of the analysed dynamical systems, the Lyapunov exponent value λ's dependence on (k, N) is calculated. The results are shown in Figure 5b. The obtained values show that λ > 0 for the entire analysed parameter range for the mapping (19). Moreover, the λ value is roughly constant for the whole range of parameter values. These results are consistent with the value of the Lyapunov exponent from Figure 5a. They also confirm that the mapping (19) is characterised by interrupted chaotic behaviour for large ranges of parameter changes and, thus, has a much more complex behaviour than the base function Sine Map (17). for the map (19). The colour bar corresponds to the value of the exponent.

Sine-Sine Map
In [6] a sine chaotification was proposed. As an example of this method, consider the Sine-Sine Map: As in the case of the Sine Map (17), the Sine-Sine Map (21) was also used as the base function F (x), and two new maps are defined, which will be analysed further in the paper: Their design is very similar to the previous case, i.e., the first one, apart from the base function having two additional components rem(10x i−1 , 1) and rem(100x i−1 , 1). In turn, the second one has one more component equal to rem(1000x i−1 , 1).
As previously, the phase diagrams of all three cases (Sine-Sine Map (21) and its two modified versions (22) and (23)) for various fixed parameter values are presented in Figure 6a. These graphs confirm that adding multiple rem(·) operators causes a significant increase in the complexity of the base mapping. In addition, Figure 6b presents bifurcation diagrams of the analysed mappings. This figure also demonstrates that, although there are visible periodic windows for the Sine-Sine Map (21), they are absent for the modified mappings (22) and (23). Similar to the Sine Map, this means that any choice of mapping parameters in (22) or (23) will always result in chaotic behaviour. The three-dimensional phase diagram confirms the increase of complexity also in this case. These relevant plots are presented in Figure 7. In the case of the Sine-Sine Map (21), the graph in Figure 7a consists of points forming a curve, while for the map (23) (Figure 7b), points fill up a 3D region indicating that the mapping with the rem(·) operators is much more complex than the base mapping.
Similar to the sine map analysis, the Lyapunov exponent was also calculated, showing the significant increase of the system's sensitivity to initial conditions after the chaotification is applied. The plot of the Lyapunov exponent is shown in Figure 8a. The results are very similar to the Sine Map case. This figure clearly shows the presence of stable regions for the Sine-Sine Map (21) and their absence for mappings with the rem(·) operator. This figure also shows that the sensitivity to the change of initial conditions for the mappings (22) and (23) is roughly constant for the entire range of the parameter k. Thus, after the same number of iterations, the system will lose information about the initial conditions. Moreover, for these mappings, the Lyapunov exponent value λ was also determined as a function of parameters (k, N). The results are shown in Figure 8b. The obtained values show that λ > 0 for the entire analysed parameter range for the mapping (23). Moreover, λ is roughly constant for the whole range of parameter values. These results are consistent with the value of the Lyapunov exponent from Figure 8a. The obtained results confirmed that the mapping (23) is characterised by interrupted chaotic behaviour for large ranges of parameter changes and, therefore, has a much more complex behaviour than the base function Sine-Sine Map (21).

Cosine-Logistic Map
In [5] a cosine chaotification was proposed. As an example of this method, consider the Cosine-Logistic Map: As in the case of the analysed previous maps, the Cosine-Logistic Map (24) was also used as the base function F (x), and two new maps are defined, which will be analysed further in the paper: Their design is very similar to the previous cases, i.e., the first one, apart from the base function, has two additional components rem(10x i−1 , 1) and rem(100x i−1 , 1). In turn, the second one has one more component equal to rem(1000x i−1 , 1).
As before, phase diagrams of all three cases (the base Cosine-Logistic Map (24) and its two modifications (25) and (26))) for fixed parameter values are presented in Figure 9a. These graphs confirm that adding multiple rem(·) operators causes a significant increase in the complexity of the base mapping. Figure 9b presents bifurcation diagrams of the analysed mappings. This figure shows that the behaviour of the Cosine-Logistic Map (24) is more complex than the behaviour of the previous base functions, i.e., after some value of the k parameter, there are no visible periodic windows. However, the mappings (25) and (26) Figure 10a,b also depict the difference between the three-dimensional phase diagrams of the Cosine-Logistic Map (24) and the map (26). The original map has a curve diagram, but the modified map fills a wide 3D region.
The respective plots of the Lyapunov exponent are shown in Figure 11a. The results are very similar to the Sine-Sine Map. This figure clearly shows the presence of stable regions for the Cosine-Logistic Map (24) and their absence for the respective modified mappings with the rem(·) operator. This figure also shows that the sensitivity to the change of initial conditions for the mappings (25) and (26) is roughly constant for the entire range of the parameter k. Thus, after the same number of iterations, the system will lose information about the initial conditions. Moreover, for these mappings, the Lyapunov exponent value λ was determined as a function of (k, N). The results are shown in Figure 11b. The obtained values show that λ > 0 for the entire analysed parameter range for the mapping (26). Moreover, the λ value is roughly constant for the whole range of parameter values. These results are consistent with the value of the Lyapunov exponent from Figure 11a. Again, the obtained results confirm that the mapping (26) is characterised by interrupted chaotic behaviour for large ranges of parameter changes and, thus, has a much more complex behaviour than the base function Cosine-Logistic Map (24).
The design is very similar to the previous cases, i.e., the first one, apart from the base function, has two additional components rem(10x i−1 , 1) and rem(100x i−1 , 1). In turn, the second one has one more component equal to rem(1000x i−1 , 1).
As in the earlier analysis, the phase diagrams of all three cases (Renyi Map (27) and its two modifications (28) and (29)) for fixed parameter values are presented in Figure 12a. These graphs confirm that adding multiple rem(·) operators causes a significant complication in the base mapping. In turn, Figure 12b presents bifurcation diagrams of the analysed mappings. This figure shows that the behaviour of the Renyi Map (27) is also very complicated, with no visible periodic windows (the stable states occur in this map for k ∈ N). However, the mappings (28) and (29) also do not have visible periodic windows. As in the previous cases, this situation means that choosing the mapping parameters in (28) or (29) will always result in a chaotic behaviour of the system.
The three-dimensional phase diagrams are shown in Figure 13 for the Renyi Map (27) and the map (29). The Renyi Map diagram consists of piecewise linear parts, while the modified map has a phase diagram that is more dense and covers a 3D area. Similar to the analysis of the previous maps, plots of the Lyapunov exponent were also determined in this case, showing the system's sensitivity to the change of the initial conditions. The plot of the Lyapunov exponent is shown in Figure 14a. The obtained results are very similar to the previous case. This figure clearly shows the presence of stable regions for the Renyi Map (27) and their absence for mappings with the rem(·) operator. This figure also shows that the sensitivity to the change of initial conditions for the mappings (28) and (29) is roughly constant for the entire range of the parameter k. Thus, after the same number of iterations, the system will lose information about the initial conditions. Again, the Lyapunov exponent value λ was also determined with respect to (k, N). The results are shown in Figure 14b. The obtained values show that λ > 0 for the entire analysed parameter range for the mapping (29). Moreover, λ is roughly constant for the whole range of parameter values. These results are consistent with the value of the Lyapunov exponent from Figure 14a. The obtained results confirm that the mapping (29) is characterised by interrupted chaotic behaviour for large ranges of parameter changes and, thus, has a much more complex behaviour than the base function Renyi Map (27).

Effect on Key Space
Let K denote the key space of a seed map, that is the sum of all possible key parameter combinations for that map, which can result in unique solution trajectories. The proposed chaotification technique (2) introduces 2m additional parameters, a j , N j , j = 1, . . . , m. Assuming a D digit accuracy, this results in a multiplicative increase of the key space by 10 2mD . If we let D = 16, the key space increase is 10 32m ≈ (10 3 ) 10.6m ≈ (2 10 ) 10.6m = 2 106m . In the proposed examples, there are three added terms, so m = 3, and the increase is approximately 2 318 . Table 1 lists the key spaces for all the considered maps. Table 1. Key space of the considered maps.

Operations per Iteration
Here, the number of operations that the technique introduces in a single iteration is counted. The results are summarised in Table 2. For m remainder terms, the chaotification technique requires 3 m additional operations, so the computational complexity of the chaotification technique increases linearly with the number of remainder operators.
Note that, in Table 2, operations that can be precomputed only once are not counted. For example, the term πk in the Sine-Sine Map (21) only needs to be computed once and not repeated in every iteration.

B-Spline Basics
A B-spline curve of degree n is characterised by its control points. Given p + 1 control points P i ∈ R d in the d-dimensional space, the curve is described by where N n i (t) are the B-spline basis polynomials of degree n. The curve is defined over a knot sequence T in the interval t ∈ [0, 1], which we assumed here as uniform, described by The repetition of the first and last knots is performed so that the curve interpolates the first and last control points and is also tangent to the control polygon.

Encryption of a B-Spline Curve
As described in the previous section, the information required to completely describe a B-Spline curve is its control points P 0 , . . . , P p arranged in order, the polynomial degree n, and a knot sequence. Thus, the encryption of a curve should aim at masking these characteristics. As changing the polynomial order would effect the proper definition of the knot sequence, here we assumed that the encryption process will not affect the degree of the curve. The same holds for the knot sequence. The proposed encryption will not affect it, as it is fairly common to assume a uniform knot sequence. If required though, an optional step can be added, encrypting the knot sequence as well.
The procedure is described in the steps below. There are two main steps involved: a shuffling of the plaintext values and a modulation through addition. The purpose of the shuffling step is to disperse the information regarding the sequence of control points, for all their coordinates. The second step hides this information, by combining it with the values of a chaotic map. For both steps, a chaotic map is used, chosen from the ones presented in the previous section. Here, the modified Sine Map (19) is considered as an example.
Here, it is worth noting that a common, and very effective substitution action is the exclusive or (XOR) operation, which is performed on the binary level of the plaintext. Here, since the plaintext consists of rational numbers, the application of the XOR operation is possible only after transforming to binary format, for example using the IEEE-754 standard. However, as this is deemed time-consuming, the more computationally efficient option of additive modulation was chosen:

Input.
A plaintext B-spline curve described by its control points P 0 , . . . , P p ∈ R d arranged in a matrix Π = (P 0 , . . . , P p ) ∈ R d×(p+1) , the polynomial degree n, and a knot sequence U. It is assumed that all the control point coordinates lie in the interval A chaotic map of the form (19), with key parameters k, a 1 , a 2 , a 3 , N 1 , N 2 , N 3 , to be computed below.

Key setup.
To make the design resistant to plaintext attacks, the map keys should be plaintext-dependent. This achieves the confusion task in encryption. They are initialised as a 3 = 1000 (37) where || · || denotes the Euclidean norm.
Step 1. The first step is to shuffle the sequence of the control points. This is also part of the usual confusion step in encryption, which disperses the plaintext information.
For the shuffling operation, the following pseudo-random number generator (PRNG) is utilised: which generates integers in the interval [1, d · (p + 1)], until d · (p + 1) discrete values are generated. If an integer is generated twice, it is discarded the second time. The result is a list of d · (p + 1) non-repeating integers, denoted as S = (s 0 , . . . , s d·(p+1)−1 ). These integers represent the new order of the elements of the vector Π row . Therefore, the s 0 element of Π row is moved to the first position, the s 1 element of Π row is moved to the second position, and so on. The resulting matrix, which contains the re-arranged elements of Π row , is denoted as Σ row : Step 2. The next step is the diffusion of the shuffled plaintext. The values of the vector Σ row are masked using the values of the chaotic map used in Step 1. Continuing the iteration of the map from Step 2 (without making use of the PRNG), the elements of Σ row are masked as follows: where the i index in x i is used to denote that the map is being iterated, continuing from the previous Step 1. The resulting values are arranged in a vector: Step 3. After the substitution step is finished, the elements of Ψ row are reshaped into an d × (p + 1) matrix Ψ = (E 0 , . . . , E p ), where each column represents a ciphertext control point, similar to the matrix Π.

Output.
A ciphertext B-spline curve described by its control points E 0 , . . . , E p , the polynomial degree n, and a knot sequence U.

Decryption Process
The decryption process follows the exact same steps, in reverse order. These are described below:

Input.
A ciphertext B-spline curve described by its control points E 0 , . . . , E p , the polynomial degree n, and a knot sequence U.
After this process is finished, the map is iterated again d × (p + 1) times, to generate the values x i , . . . , x i+d·(p+1)−1 . The i index here denotes that the map is iterated consecutively from the previous process of generating S.
Step 3. The confusion is reversed as The resulting elements are combined in a vector as Eqauiton (42): Step 4. The next step is to reverse the shuffling of the elements of Σ row .
The elements of S from Step 2 indicate the shuffling that was originally performed in the encryption process. The procedure is reversed, by moving the first element into the s 0 position, the second element into the s 1 position, and so on, until the last element is moved into the s d·(p+1)−1 position. The new matrix is denoted as Π row .

Output.
A plaintext B-spline curve described by its control points P 0 , . . . , P p ∈ R d arranged in a matrix Π = (P 0 , . . . , P p ) ∈ R d×(p+1) , the polynomial degree n, and a knot sequence U.

Patch Encryption
B-spline patches are obtained by considering a net of control points P i,j ∈ R 3 , i = 0, . . . , P, j = 0, . . . , Q, and two B-spline polynomial bases N n i (u), N m j (v) of degrees n, m, defined over knot sequences U, V as in (31). The B-spline patch is defined as or in matrix form: To encrypt a B-spline patch, the procedure followed is a direct extension of Section 6.2. The idea is again to encrypt the control points of the patch, and this is performed following the exact same steps. To do so, the only modification required is to concatenate all the control points in a single row vector. This can be performed by arranging the coordinates of all points in a row Π row similar to (32) as follows: Π row = P 0,0,1 P 0,0,2 P 0,0,3 P 1,0,1 · · · P P,Q,3 , where the third index in P i,j,1 , P i,j,2 , P i,j,3 denotes the x, y, z coordinates, respectively. With this modification and the appropriate reshaping of Ψ row , Section 6.2 can be applied to encrypt a B-spline patch. The same process can be performed when there is a collection of curves or patches to be encrypted. Since in (33), the encryption keys are plaintext-dependent, repeating the process of key generation for a collection of curves or patches would be impractical. This can be easily resolved, by arranging first all of the given control points into a single vector Π row , as in (32) or (52), and then, performing the encryption on the single vector Π row , which carries the plaintext information on all the curves. Of course, the appropriate reverse reshaping would have to be performed during the decryption stage.

Simulation Results
The encryption design is tested on a collection of several 2D and 3D cubic B-splines. Here, B = 5. Figures 15-17 show a collection of cubic B-spline curves to be used for encryption, along with their convex hulls, and the corresponding encrypted curves. The encrypted curves have no visual similarity to the corresponding plaintext ones. The curves in Figure 16 are taken from [31]. The curves were generated using the algorithms described in [30,40].   [31], along with their convex hull, the encrypted curve (middle), and both curves (right). Figure 18 shows a collection of three patches and the respective encrypted ones. The encrypted patches are again visually dissimilar from the corresponding plaintext patch. Moreover, the encryption of the Utah teapot, teacup, and teaspoon is considered, with data taken from [41]. The data available consist of a .txt file of control points and a secondary .txt file indicating the control point indexes in the first file, for each patch. The results are shown in Figure 19. The teapot consists of 32 patches; the cup is 26 patches; the spoon is 16 patches. The resulting ciphertext collection of patches clearly has no distinguishable shape. Note that, for the ease of implementation, the reshaping of the control points into a single vector Π row (52) is performed directly to the .txt file provided in [41].
In the following subsections, a series of statistical tests were performed on the plaintext and ciphertext curves. For this, the curves in Figure 16 and the shapes in Figure 19 were chosen, due to having the most control points. The "G" curve consists of 50 control points; the music key consists of 30 control points; the sinusoidal curve consists of 50 control points. For the teapot, all 32 patches have overall 306 control points; the teacup has 251 points; the teaspoon 251 control points.

Map Choice
In (34), the map parameter k is chosen based on the value where the curve lies. Depending on the map, the choice of k can affect the map's range. This is the case for the sine map (19) and the Cosine-Logistic Map (24). Having a map with a wider range will affect how wide the resulting curve will be, after the modulation step. Naturally, the map's range should be wide enough to hide the plaintext curve's range after the modulation. See, for example, Figures 15 and 16, where the encrypted curves lie on a much wider area in the plane, effectively hiding the range of the plaintext control points.
If a map's control parameter k cannot affect its mapping range, as for example in the sine-sine (21) and Renyi (27) maps, then before the additive modulation step, the chaotic values can be scaled by an appropriate factor.

Histogram Analysis
Any statistically unbiased random time series should have a histogram that is uniform, indicating that there is no bias in the average values that occur over a specific subinterval of the state space. Here, the histograms were computed for the concatenated coordinate vectors Π row , Ψ row . The row Σ row , being a shuffling of Π row , has the exact same histogram, so it is not shown here. The histograms are displayed in Figure 20. For the curves, the histogram consists of 10 groups, while for the surfaces, there are 30 groups. The ones corresponding to the ciphertext curves/patches are clearly more uniform compared to the plaintext ones. Table 3 also displays the histogram variance for three curves and shapes. The variance is much lower for the ciphertext points, indicating a more uniform histogram. Notice, however, that this analysis was performed in a very small ciphertext and, therefore, should only be interpreted as a measure of improved confusion in the ciphertext.

Correlation
Correlation indicates whether there is any linear relationship between two datasets. For two independent datasets, their correlation should be very close to zero, indicating that their values are linearly independent. Here, the correlation was measured between consecutive control points in the plaintext and ciphertext curves. Note here that, as a plaintext curve can have any form or shape, it is possible that the consecutive control points of the plaintext curve may also showcase low correlation. However, even in this case, it is expected that the correlation of the ciphertext control points will be even lower. The results are shown in Table 4.

Noise
The process of decryption described in Section 6.3 can be performed even if the ciphertextΨ = Ψ + W has been corrupted by noise W during transmission. The noise will still be present in the decrypted set of control points though. When the channel noise is too intense, the decryption process can be independently complemented by a suitable filtering process to reduce the noise. Figure 21 shows the result of decryption from a ciphertext setΨ of control points under different noise intensities. The noise was modelled as uniform random numbers in the interval [−z, z], where z = 0.001, 0.01, 0.1 for the "G" curve and z = 0.01, 0.1, 0.5 for the teapot. The error norms between the plaintext and decrypted control points are 0.0041, 0.0380, and 0.4505, respectively, for the "G" curve and 0.1043, 1.0768, 5.1801 for the teapot.

Key Space
The key space of the encryption (Section 6.2) is the same as the key space of the source map. As shown in Table 1, all the maps have the required key space to resist brute force attacks.

Key Sensitivity
The encryption process inherits the avalanche property of the chaotic map it uses as its source of randomness. As any chaotic map showcases sensitivity to parameter changes, perturbing any key value in (33) in the decryption process will result in the failure of the decryption. Figure 22 shows two such examples, were it can be observed that, indeed, the decryption fails to recreate the original plaintext curve and patch. Figure 22. Examples of plaintext (red, left), encrypted (green, middle), and decrypted (black, right) curve/patch, where the key value x 0 is perturbed asx 0 = x 0 + 10 −9 during the decryption.

Conclusions and Future Goals
This work studied a chaotification technique, which is an extension of [18] and consists of combining the remainder of multiple scalings of the map's value for the next iteration. The resulting maps can achieve a high Lyapunov exponent value, which can be controlled by appropriate parameter tuning.
An application of the proposed maps to the encryption of B-spline curves and patches was considered. The encryption consists of two steps: the first is the shuffling of the control points, and the second is an addition modulation step.
Future studies will be directed towards two main research areas. The first is the problem of chaotification. The use of the remainder operator will be explored for the creation of chaotic maps with a controllable mapping domain and without any equilibria. The control of the mapping interval will be of use in developing maps with modifiable statistical characteristics, which is desirable in various applications. The absence of equilibria is also desirable for achieving robust chaotic behaviour, as the Period 1 behaviour completely vanishes. Moreover, maps and systems without equilibria are themselves a distinct topic of study in dynamical systems' analysis.
A second extension of this work is to generalise the encryption for rational B-spline curves and patches [30]. Here, the process should include the encryption of the control points along with the corresponding weights. Moreover, other types of geometric data structures can be considered, for example 3D objects. Finally, apart from the additive modulation, which is a linear operator, other operations can be considered for increased complexity, such as multiplication or cosine transformation. It must be noted though that highly nonlinear operations may introduce numerical errors during the decryption, so their choice must be carefully made. Any additional encryption steps can be considered either directly on the floating point plaintext data or in their binary representation. It would also be interesting to study if chaotic maps can be used for watermarking B-spline curves. Funding: This research received no external funding.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.