Chapter 1

Introduction 1.1 Context

Our research is conducted as part of the project ANR ModItère. The global objective of this project is to develop a new type of geometric modellers for computer-aided design systems (CAD), based on iterative methods following the principle of fractal geometry.

In CAD systems, objects are often modelled based on the classic matching process, that is, by combining different simple parts. Defining shapes by iteration allows us to generate new structures with specific properties, such as roughness or lacunarity, which cannot be achieved with classic modelling. However these types of objects are conceivable and may be of interest: porous structures can be used for lighter objects while maintaining satisfactory mechanical properties, rough surfaces can be used for acoustic absorption (see figure 1.1).

We propose to elaborate a set of modelling tools to adapt the fractal model, keeping the facilities of existing CAD systems, while extending their capabilities and application areas. Thus we provide algorithms for CAD systems to design new kind of shapes based on the concept of fractal geometry in order to model and fabricate these shapes.

Problems

Fractals are self-similar shapes composed of several parts, each of which is similar to the entire shape. In mathematics, fractals are sets of points in a metric space with a non-integer dimension (in the sense of Minkowski or Hausdorff metric), that is, the metric dimension differs from the topological one.

Most of existing CAD modellers are based on classic geometry. They permit modelling analytic curves and surfaces, which are usually polynomial or rational and smooth. The shapes modelled by iterative methods are generally not smooth because of their fractal structure. The topology can be controlled in order to build curves and surfaces represented by parametrized functions, but these functions are usually nowhere differentiable.

Generally, it is difficult or impossible to represent these shapes using classic geometric models, although this iterative construction principle is exploited to evaluate certain objects, such as subdivision surfaces [START_REF] Warren | Subdivision methods for geometric design : a constructive approach[END_REF][START_REF] Lanquetin | Étude des surfaces de subdivision: intersection, precision et profondeur de subdivision[END_REF]. Specific iterative models were developed to generalize the answer to this problem: Lsystem [START_REF] Prusinkiewicz | The algorithmic beauty of plants[END_REF], IFS [START_REF] Barnsley | Fractals everywhere[END_REF][START_REF] Gentil | Les fractales en synthése d'images: le modèle IFS[END_REF], LRIFS/CIFS [START_REF] Prusinkiewicz | Language-restricted iterated function systems, koch constructions, and l-systems[END_REF], BCIFS [TBSG + 06]. These models are mainly used to generate fractal structures in computer graphics, but also to compress data [START_REF] Barnsley | Methods and apparatus for image compression by iterated function system[END_REF][START_REF] Fisher | Fractal image compression, theory and application[END_REF]. In this thesis, we use the term fractal to denote a shape represented by these iterative models.

An Iterated Function System (IFS) is determined by a set of transformations. The attractor of an IFS may be evaluated by iterating these functions. A Controlled Iterated Function System (CIFS) is a more general system that allows us to control certain parts of the IFS attractor. A CIFS denotes an IFS with restrictions on the transformation sequences imposed by a control graph and defines objects whose geometry can be complex. The attractor of a CIFS can be evaluated by an automaton [START_REF] Mauldin | Hausdorff dimension in graph directed constructions[END_REF] defined on the control graph. Each state of the automaton corresponds to different parts or sub-parts of the modelled object. Transitions between states indicate that one sub-part is contained in another one. It is then possible to control the attractor more precisely. However the evaluation, analysis and characterization (location, size) of the resulting shape can be complex.

From IFS and CIFS, various shapes can be modelled. However, it is difficult to control their topological properties. Shapes are determined by a set of geometry operators. Modifying these operators leads to both global and local changes in the shape, this affects not only geometry but also topology.

To separate global and local controls, attractors can be constructed in barycentric spaces and then projected to a set of control points in the modelling space [START_REF] Schaefer | Subdivision schemes and attractors[END_REF][START_REF] Tosan | Surfaces fractales définies par leurs bords[END_REF][START_REF] Gouaty | Modélisation géométrique itérative sous contraintes[END_REF][START_REF] Zair | Fractal modeling using free form techniques[END_REF]. The global shape can therefore be modified by moving the control points, while the local geometry (geometric texture) is defined by subdivision operators, i.e. the CIFS transformations.

In order to control the topological structure of the modelled shapes, a CIFS is enriched by integrating a topological model to obtain a Boundary Controlled Iterated Function System (BCIFS) [TBSG + 06]. It is thus possible to specify incidence and adjacency constraints [START_REF] Tosan | Surfaces fractales définies par leurs bords[END_REF][START_REF] Gouaty | Modélisation géométrique itérative sous contraintes[END_REF] for the subdivision and thereby to control the resulting topology: classic topology (curve, surface,. . . ) or fractal topology. The topological model in this case is more general than the classic one [START_REF] Gouaty | Modélisation géométrique itérative sous contraintes[END_REF]. A topological cell may be a fractal object. For example, a face can be a Sierpinski triangle or an edge can be a Cantor set, but the topological structure remains consistent.

Fractal shapes can thus be modelled using BCIFS, but shape modelling in CAD systems is complex, that is not confined to simply representing the shape. To assist the user, geometric modellers are enriched by constructing operations and algorithms. These operations are usually based on geometric properties of modelled objects. Since for fractal shapes these properties are more complex and more subtle [Ben09,BGN08] and sometimes even not defined, existing CAD systems are not suitable to manipulate fractal structures. On the other hand, Zaïr and Tosan [START_REF] Zair | Fractal modeling using free form techniques[END_REF] showed that the result of certain operations can be obtained formally. For example, a tensor product of two fractal curves can be obtained by forming the tensor product of BCIFS operator matrices.

Today, there is no CAD software to design or to model fractals. This is justified by the fact that this problem was not really of interest. Indeed, there was no manufacturing process able to produce such objects. However, the emergence of new techniques, such as stereo-lithography or laser fusion, opens up new possibilities yet unused and even unexplored.

The passage of the theoretical model to the physical object requires a number of treatments. For example, curves and surfaces must be thickened (dilated with respect to Hausdorff distance) to become 3-dimensional. This must be done without loss of rendering accuracy in order to preserve the aesthetics and geometric properties such as topology (lacunarity) and differentiability (roughness). The generated shapes must thus have a good translation into the physical world.

The goal of our project is to develop a set of tools and algorithms that permits one to evaluate, to characterize and to analyse different geometric properties (localisation, convex hull, volume, fractal dimension) of fractal shapes. It is thus necessary to be able to evaluate a set operations (intersection, union, . . . ) between fractal objects and classic CAD objects (splines, NURBS, cylinders, . . . ). Finally, in order to fabricate these objects, it is also important to take into account the manufacturing constraints. The idea is to integrate these constraints as early as possible, i.e. during design. This leads to the following advantages:

• the sequence conception-simulation-fabrication is more direct;

• the manufacturing time is shortened;

• the creativity of users (designers, architects, . . . ) is improved by opening a wide range of modelling possibilities, that were not free in existing CAD modellers;

• this creativity is better exploited, i.e. designed shapes are closer to the constructed shapes, there are fewer losses and degradation of information.

Organization of thesis

This thesis consists of six parts and is organised as follows.

In chapter 2 we describe popular existing tools for geometric modelling, and mention the various constructing operations used to assist the user in computeraided shape modelling. There exist many techniques to generate solid objects, such as sweeping 2D surfaces, solid modelling, free-form surfaces. We describe basic ideas of these techniques as well as discuss approaches in the iterative geometric modelling. After listing the constructing operations presented in most of CAD systems, we identify the basic properties required for their realisation. Preliminary analysis shows that CAD operations are based on four properties: affine invariance, topological structure, parametrization and differential properties. The result is given in the form of a dependency graph.

Chapter 3 defines IFS, CIFS and BCIFS. We discuss mathematical aspects that allow us to define irregular curves and surfaces. We also identify a system of topological constraints and the notion of control polygons.

In chapter 4 we classify CAD operations adapting to BCIFS and we also introduce our algorithms to approximate the operation. The approximation algorithms are generic, i.e. defined for an arbitrary operation satisfying some constraints that we explicit. These algorithms can thus be applied once we implement the required set of operators (interface). The proof of convergence of these algorithms is based on the results of domain theory.

In chapter 5 we define a self-similarity property of the operation and introduce a generic algorithm to compute a formal representation of this result, i.e. to compute a specific CIFS with the operation image as the attractor. Most of CAD operations are semi-decidable, that is the result can be computed only approximately. The presenting algorithm computes a CIFS automaton for the required accuracy. The attractor of this CIFS is thus close enough to the required operation image. In special cases, when the automaton degenerates into an single state, the constructed model represents a specific IFS and has a linear complexity on the number of iterations.

Chapter 6 presents various examples and special cases and summarizes all the results. After a general conclusion we announce prospects of our work and discuss some open questions.

Appendix A describes the methods we used to compute a bounding ball (not necessarily a tight one) of the IFS attractor. For CIFS attractors we thus compute a set of bounding balls. We use these approaches in the initialization steps of presented approximation algorithms.

Chapter 2

CAD operations

Shape modelling in CAD systems is enriched by constructing operations and algorithms. These operations are usually based on geometric properties of modelled objects. However, for fractal shapes, these properties are more complex and sometime not defined. Since most of existing CAD modellers are based on classic geometry, the modellers are not suitable to manipulate fractal structures. We study how standard CAD operators can be defined and evaluated on fractals.

Our goal is to identify the different cases that arise when we want to implement standard CAD operations on fractals. For this, we list the main standard CAD operations presented in most of CAD systems and identify the properties on which these operations are based. We then identify the basic properties required for their realisation. Preliminary analysis shows that operations are based on four properties: affine invariance, topological structure, parametrization and differential properties. The result is given in the form of a dependency graph.

Finally, we discuss existing iterative models, such as subdivision surfaces and L-Systems, that define shapes by iteration. There are various existing algorithms to manipulate these models, that will be described, before concluding.

Geometric modelling

In this section we identify general tools widely used for geometric modelling. This survey does not pretend to be exhaustive; we considered only the methods related to our work. There exist many techniques to generate solid objects, such as sweeping 2D surfaces, solid modelling, free-form surfaces. We describe basic ideas while details are available in many resources [Req77, Chi88, Hof89, HR91, Gal00, Rao04].

Sweep or Extrusion

There are number of useful techniques for creating objects, such as extrude, revolve, sweep, loft. These 3D geometric constructing methods extended from 2D, normally can be combined into two classes:

• linear extrusion or translational sweep,

• rotational sweep.

In general, for extrusion sweeping a 2D surface along a curve (see figure 2.1) generates a 3D object. In linear extrusion, sweeping is applied along a straight line. There exist further variations in sweep for generating more complex geometry. For example, it is possible to sweep in a linear direction with a taper along the direction.

There are some cases where a sweep may fail to generate a valid solid. For example, a line when swept transforms to a surface, but not a solid. Sweeping must therefore forbid such cases by imposing some restrictions on the topology and the geometry chosen for linear extrusion in order to produce a valid solid. Another type of construction is the rotational sweep (sometimes called "revolve" or "surface of revolution"), which can be used for axi-symmetric shapes, such as bottles, glasses, vases and so on. It is also possible to add twist to sweeping in the third dimension. Similarly, the rotational sweep can be enhanced by adding axial or radial offset while sweeping to get helical or spiral objects.

Solid modelling

A well-known method for 3D solid constructing is the solid modelling technique, often called "constructive solid geometry" (CSG). The idea is to create complex objects by adding or subtracting predefined solid primitives, as illustrated in figure 2. [START_REF]Journées du Groupe de Travail en Modélisation Géométrique[END_REF].

For combining primitives one usually use basic set operators, called Boolean operators. There are following classic boolean operations:

• Union of the sets A and B, denoted A ∪ B, is the set of all objects that are a member of A, or B, or both.

• Intersection of the sets A and B, denoted A ∩ B, is the set of all objects that are members of both A and B.

• Set difference of A and B, denoted A \ B, is the set of all members of A that are not members of B. • Symmetric difference of sets A and B, denoted A B, is the set of all objects that are member of exactly one of A or B (elements which are in one of the sets, but not in both): (A ∪ B) \ (A ∩ B)

• Cartesian product of A and B, denoted A×B, is the set whose members are all possible ordered pairs (a, b) where a is a member of A and b is a member of B.

CSG is a powerful technique for representing fairly complex objects relatively easily. The boolean operators always guarantee that the object formed by these rules is physically constructable. However, note that these operations can affect both the geometry and the topology. To produce valid solid objects, one should use regularized boolean operations [START_REF] Gomes | Design and Implementation of 3D Graphics Systems[END_REF].

Free-form surfaces

Free-form surfaces are generally created by manipulating the control points, as illustrated in figure 2.3. Surfaces are not stored or defined in CAD software in terms of polynomial equations, but by their control points and a surface degree. The degree of a surface determines its geometric properties.

There are many types of free-form surfaces that are widely used in CAD systems, such as ruled, Bézier or B-spline surfaces. Some operations (chamfering, rounding/filleting) are applied directly to the model faces, edges and vertices to create a desired modification. A lofted surface is a surface constructed by transitioning between two or more curves by using a smooth blend between each section of the surface. Surface fillet is a tool to automatically generate the fillet radius between two surfaces, which could be uniform or vary linearly.

Transforms

All the modelling techniques have a number of simple and natural geometric transformations. There is an important class of geometric transformations, called affine transformations. In geometry, an affine transformation is a transformation that preserves straight lines (i.e. images of all points lying on a line, also lie on a line) and ratios of distances between points lying on a straight line (that is, the midpoint of a line segment remains the midpoint after transformation). Other examples of CAD operations simplifying geometric modelling are "copy-paste" and "array".

Operations graph

To exploit fractals in CAD systems we must be able to evaluate most standard constructing operations for fractal objects. In this section we study the constructing operations presented in the CAD system Rhino. The choice of this software was suggested by our industrial partners.

Some of these operations are composites, that is, we have to calculate several elementary operations and then combine results. For example, to determine a G 1 connection between two curves, we must be able to calculate the tangents at each curve end. But we also must be able to construct a curve with these imposed tangents.

For preliminary analysis we structure operations by determining dependencies between them. After listing the CAD operations, we construct a diagram illustrating these relations in the form of a graph (see figure 2.4). This graph illustrates the dependencies between the CAD operations and object properties. The basic properties and methods are marked by diamonds. Note that this graph does not pretend to be exhaustive neither for operations nor dependencies. However, it is constructed by considering at least one solution for each operation.

All CAD operations are divided into two categories depending on the type of the result:

• where a formal solution is provided (depicted by rectangles)

• where an approximate solution is provided (depicted by ellipses)

We now explain these categories.

Operations with a formal solution

In this category there are operations that produce a formal representation of the required result. That is, after an appropriate transform T , the NURBS surface S, for example, remains a NURBS surface, and moreover, the result represents the exact mathematical image T (S).

For instance, we consider the affine transformations. Let C be a NURBS curve, defined by a set of control points {P i | i = 1, . . . , k} as follows:

C(u) = k i=1 R i,n (u)P i ,
where R i,n are the rational basis functions:

R i,n (u) = N i,n (u)w i k j=1 N j,n (u)w j ,
with N i,n the B-spline basis functions, n the degree of the B-spline basis functions and w i the corresponding weights. The denominator is a normalizing factor. Every affine transformation is the composite of two functions: a translation and a linear transformation. The linear part is represented as a multiplication by a matrix L and the translation as the addition of a vector b, an affine transformation T acting on the curve C can thus be represented as

C (u) = T (C(u)) = b + LC(u) = b + k i=1 R i,n (u)LP i
Since all the R i,n sums to 1:

C (u) = k i=1 R i,n (u)b+ k i=1 R i,n (u)LP i = k i=1 R i,n (u)(b+LP i ) = k i=1 R i,n (u)T (P i ).
To apply an affine transformation to a NURBS curve, we can thus simply transform the control points of the curve. Moreover, the new curve obtained by transforming the control points is also a NURBS curve. We have thus the formal representation of the image under affine transformation.

Another example is determining curve endpoints. This can be done in a natural way from the parametrisation of a curve by passing the bounding values of a parameter.

Operations with an approximate solution

Sometimes it is difficult or even impossible to preserve the representation of an object after applying a CAD operation. This category consists of operations for which the result is computed approximately. A simple example is intersecting between subdivision surfaces [START_REF] Lanquetin | A graph based algorithm for intersection of subdivision surfaces[END_REF][START_REF] Biermann | Approximate boolean operations on free-form solids[END_REF]. There is no method to represent formally the result. CAD systems therefore modify the initial representation of objects after computing this operation.

Iterative modelling and manipulations

There exists another approach to represent geometric solids, called iterative modelling, where we define the shape implicitly by an algorithm that generates it. These geometric models are usually defined by a set of rules which are applied to generate the ensuing iteration in a recursive manner, in contrast to shapes represented by boundary surfaces and volumetric models (CSG and Spatial-partitioning). However, these distinctions are often blurred: for instance, geometric shapes such as circles can be defined by implicit mathematical equations, and a fractal model reduces to a constructive model when its recursive definition is truncated to a finite depth.

Iteration allows us to generate structures with complex properties, such as roughness or lacunarity, which cannot be achieved with classic modelling. In this section we discuss existing iterative models, that define shapes by iteration. There are various existing algorithms to manipulate these models, which will also be described.

Subdivision surfaces

One powerful iterative method for representing shape is called subdivision. Subdivision surfaces allow smooth free-form surface modelling without topological constraints. Subdivision have become a fundamental representation for smooth geometry, particularly in the animated movies.

Subdivision defines smooth curve or surface as the limit of a sequence of successive refinements of a mesh. The geometric domain is represented by piecewise linear objects, usually polygons or polyhedra. The iteration consists of two distinct phases:

• Splitting: Each edge or face is split into edges or faces.

• Averaging: Each new vertex introduced by splitting is positioned at a fixed affine combination of its neighbour's positions. The initial points could also be repositioned.

Subdivision thus produces a hierarchy of polyhedra that approximate the final limit shape. This process is used mostly to produce highly smooth surfaces from polyhedra of arbitrary topology. The main benefit of subdivision is that it is simple to implement. In order to implement a subdivision scheme, each vertex of a shape is tagged by specifying whether the descendants of the vertex lie on a vertex, edge or face of the final limit shape. During subdivision, the appropriate averaging mask can be chosen based on this tag. By truncating to a finite depth the curved limit shape can be approximated to any desired tolerance.

There exist a large number of different subdivision schemes: schemes that operate on quadrilateral [START_REF] Catmull | Seminal graphics. chapter Recursively generated B-spline surfaces on arbitrary topological meshes[END_REF][START_REF] Doo | A subdivision algorithm for smoothing down irregularly shaped polyhedrons[END_REF] or triangular meshes [START_REF] Loop | Smooth Subdivision Surfaces Based on Triangles[END_REF], schemes whose limit surface either interpolates a control mesh [START_REF] Dyn | A butterfly subdivision scheme for surface interpolation with tension control[END_REF][START_REF] Kobbelt | Interpolatory subdivision on open quadrilateral nets with arbitrary topology[END_REF] or approximates it [START_REF] Catmull | Seminal graphics. chapter Recursively generated B-spline surfaces on arbitrary topological meshes[END_REF][START_REF] Doo | A subdivision algorithm for smoothing down irregularly shaped polyhedrons[END_REF], as well as schemes that incorporate a rotation of the grid directions [PR97, Kob00, LG00, VZ01]. There are also schemes that use the flexibility of subdivision to mix refinement patterns in the same mesh [START_REF] Stam | Loop. Quad/triangle subdivision[END_REF][START_REF] Schaefer | On c 2 triangle/quad subdivision[END_REF][START_REF] Peters | Combining 4-and 3-direction subdivision[END_REF]. Comprehensive classifications of refinement patterns are provided by Cashman [START_REF] Cashman | Beyond catmull-clark? a survey of advances in subdivision surface methods[END_REF], Han [START_REF] Han | Classification and construction of bivariate subdivision schemes[END_REF] and Ivrissimtzis et al. [START_REF] Ivrissimtzis | A generative classification of mesh refinement rules with lattice transformations[END_REF].

Manipulation with subdivision surfaces

Given a subdivision scheme, the fundamental difficulties are to determine the properties of the limit shape, to convert other representations, such as NURBS or CSG, into a subdivision representation, and to perform geometric operations, such as intersection, lofting or fairing.

T. Cashman [START_REF] Cashman | Beyond catmull-clark? a survey of advances in subdivision surface methods[END_REF] surveys the high-performance evaluating techniques for subdivision surfaces on both CPU and GPU which makes it possible to consider subdivision surfaces in real-time applications such as computer games [START_REF] Kovacs | Real-time creased approximate subdivision surfaces[END_REF]. The survey also includes various techniques to analyse the first and the second-order smoothness of subdivision surfaces.

To trim subdivision surfaces it is usually required to compute the parametric pre-image of a trim curve in order to exclude for evaluation the relevant parts of the surface. Litke et al. [START_REF] Litke | Trimming for subdivision surfaces[END_REF] pointed out that with Levin's [Lev99] combined subdivision scheme, which can satisfy boundary constraints, it is possible to meet a desired trim curve exactly without computing an exact pre-image. The compromise they make is to modify (within a specified tolerance) the surface near the trim curve.

Much work has been devoted to interference detection and surface intersection for analytical surfaces [START_REF] Rappoport | Interactive boolean operations for conceptual design of 3-d solids[END_REF] and NURBS surfaces [KKM97, KKM99, KM97], however few previous works have attempted to solve this problem for subdivision surfaces. Subdivision surface intersection often involves truncating the evaluation tree to a finite depth and then applying a mesh intersection algorithm without considering properties of subdivision [START_REF] Biermann | Approximate boolean operations on free-form solids[END_REF], which can be costly to compute. The problem of interference detection was approximately solved for a large class of subdivision surfaces [DKT98, GS01, WP04].

Lanquetin et al. [START_REF] Lanquetin | A graph based algorithm for intersection of subdivision surfaces[END_REF] proposed an algorithm to intersect the subdivision surfaces generated by the Loop scheme. The algorithm was later improved by Severn and Samavati [START_REF] Severn | Fast intersections for subdivision surfaces[END_REF] by using the strong convex hull property in order to produce a fast intersection for subdivision surfaces with arbitrary precision. This approach can therefore be used with any subdivision scheme that has this strong convex hull property. Biermann et al. [START_REF] Biermann | Approximate boolean operations on free-form solids[END_REF] presented a method for approximating boolean operations for subdivision surfaces. There were attempts to improve robustness by using voxelization [LC07] or by operating on the limit mesh [START_REF] Jiang | Numerical validation in current hardware architectures[END_REF].

J. Stam [START_REF] Stam | Exact evaluation of catmull-clark subdivision surfaces at arbitrary parameter values[END_REF] showed that Catmull-Clark subdivision surface and all its derivatives can be evaluated in terms of a set of eigenbasis functions which depend only on the subdivision scheme and he derives analytical expressions for these basis functions. He implemented a technique to compute high quality curvature plots of subdivision surfaces. The cost of the evaluation scheme is comparable to that of a bi-cubic spline. Therefore, this method allows many algorithms developed for parametric surfaces to be applied to Catmull-Clark subdivision surfaces.

Subdivision schemes generate self-similar curves and surfaces. S. Schaefer et al. [START_REF] Schaefer | Subdivision schemes and attractors[END_REF] demonstrated that subdivision curves and surfaces are in fact special cases of the IFS attractors, that is, they can be generated by iterating contractive transformations. The authors derived the associated IFS for many different subdivision curves and surfaces without extraordinary vertices, including B-splines, piecewise Bezier, interpolatory four-point subdivision, bi-cubic subdivision. Moreover they showed how to build subdivision schemes to gener-ate traditional fractals such as the Sierpinski triangle and the Koch curve.

L-systems

An L-system or Lindenmayer [START_REF] Lindenmayer | Mathematical models for cellular interaction in development: Parts i and II[END_REF] system is a type of formal grammar. A systematic study and formalisation was undertaken by P. Prusinkiewicz [START_REF] Prusinkiewicz | L-systems in geometric modeling[END_REF][START_REF] Prusinkiewicz | The algorithmic beauty of plants[END_REF]. An L-system consists of an alphabet of symbols, a collection of production rules that expand each symbol into some larger string of symbols, an initial string from which to begin construction, and an algorithm to translate the generated strings into geometric structures. The central concept of L-systems is the notion of rewriting. Rewriting is a technique to build complex objects by replacing parts of a simple initial object using local rewriting rules. L-systems describe growing structures in terms of automatically maintained topological (neighbourhood) relations between components. The context for applying the rules is therefore always available.

Typically, L-systems are characterized by repetitive applying simple rules to discrete structures with a changing number of components. More recently, interpretations of L-systems based on affine geometry have been presented [START_REF] Derose | Theory and practice of geometric modeling. chapter A coordinate-free approach to geometric programming[END_REF][START_REF] Goldman | On the algebraic and geometric foundations of computer graphics[END_REF]. These interpretations provide a brief description of subdivision curves. In fact, the generating algorithm is similar to the one for generating fractal curves and also resembles subdivision algorithms. Their modification rules can thus be described in local terms and formalized as L-systems production rules. P. Prusinkiewicz et al. [START_REF] Prusinkiewicz | L-system description of subdivision curves[END_REF] introduced parametric context-sensitive Lsystems with affine geometry interpretations as an alternative technique for specifying and generating subdivision curves. Subdivision algorithms can thus be formalized by L-systems in an intuitive manner, which leads furthermore to their computer implementation.

Such fundamental algorithms as the Lane-Riesenfeld algorithm for generating uniform B-splines of arbitrary degree, the de Casteljau algorithm for generating Bezier curves, and their extensions to rational curves, can also be expressed using L-systems with affine geometry interpretations [START_REF] Prusinkiewicz | L-systems in geometric modeling[END_REF].

I. Zammouri et al. [START_REF] Zammouri | A theoretical framework mapping Grammar based systems and Fractal description[END_REF] showed the exact mathematical relation between L-systems and iterated function systems (IFS). They proposed a new method to describe fractal shapes using parametric L-systems.

Conclusion

In this chapter we identified the general tools widely used for geometric modelling and we described the various constructing operations used to assist users in computer-aided shape modelling. There exist many techniques to generate solid objects, such as sweeping 2D surfaces, solid modelling, free-form surfaces. We described the basic ideas of these techniques.

After enumerating all the operations presented in the CAD system Rhino, we constructed a diagram that illustrates dependencies between these operations in the form of a graph. Operations were then classified into two categories depending on type of the result: a formal representation of the result or an approximate result.

The constructing operations are usually based on geometric properties of modelled objects. A preliminary analysis shows that CAD operations are based on four basic properties: affine invariance, topological structure, parametrization and differential properties. For example, in order to align objects we must be able to compute the bounding box and to translate the objects. If the objects are not affine invariant, the computed bounding box will no longer be correct after translating. Extracting surfaces of a solid or duplicating their borders can be easily computed using the boundary representation. These basic properties are therefore essential and we must find their equivalents in fractal geometry. In fact, for fractal shapes, these properties are more complex and sometime even not defined. That is why existing CAD systems are often not suitable to manipulate fractal structures.

In fact, most of existing CAD modellers are based on classic geometry. These modellers permit modelling analytic curves and surfaces, which are usually polynomial or rational and smooth. The shapes modelled by iterative methods generally are not smooth because of their fractal structure. Their topology can be controlled in order to build curves and surfaces represented by parametrized functions, but these functions are usually nowhere differentiable.

Iterative models, such as L-systems, give a recursive definition of shapes allowing us to generate structures with specific properties, such as roughness or lacunarity. Shapes with these properties cannot be achieved with classic modelling. To exploit these shapes in CAD systems, we must be able to evaluate most standard constructing operations for fractal objects.

Being a parametric and an implicit model, a fractal reduces to a constructive geometry when its recursive definition is truncated to a finite depth. The CAD operations can therefore be applied directly to an approximation of a fractal. However, the attractor approximation can be time-consuming, because the number of objects grows exponentially with an increasing number of iterations. In the next chapter we adapt CAD operations to fractals, we formalize the notion of approximation by identifying situations for which it is possible to control this approximation. We identify the constraints on operations in order to apply general or optimised algorithms. These algorithms can thus be applied once we implement the required set of operators (interface).

Together, the principles of iterative, geometric and solid modelling form the foundation of computer-aided design and in general support the creation, visualization, animation and annotation of digital models for physical objects.

Chapter 3

Iterated Function Systems and generalisations

In this chapter we define the iterative systems that will be used to represent fractal objects. We discuss mathematical aspects that allow us to define irregular curves and surfaces. We also identify a system of topological constraints and the notion of control polygons.

IFS

Iterated Function Systems (IFS) were introduced by Hutchinson [Hut81] as a strictly mathematical model. Barnsley further developed and popularized this model for applications in computer graphics [START_REF] Barnsley | Fractals everywhere[END_REF]. Much research has followed on from this idea [START_REF] Gentil | Les fractales en synthése d'images: le modèle IFS[END_REF][START_REF] Massopust | Fractal functions and their applications[END_REF][START_REF] Barnsley | V-variable fractals: Fractals with partial self similarity[END_REF][START_REF] Tosan | Surfaces fractales définies par leurs bords[END_REF].

IFS are based on the self-similarity property. A modelled object is made up of the union of several copies of itself; each copy is transformed by a function. These functions are usually contractive, that is, bringing points closer together and making shapes smaller. Hence, the modelled object, called attractor, is made up of several possibly-overlapping smaller copies of itself, each copy also made up of copies of itself, ad infinitum.

Definition of IFS

In this section we formally define IFS and we thus establish the notation used in this thesis.

Definition. Given a complete metric space (X, d) with the associated metric d, an IFS is defined by a finite set of continuous transformations T = {T i } N -1 i=0 in the space X. Let Σ = {0, . . . , N -1} be the set of IFS transformation indices, thus |Σ| = N . The IFS is then denoted by {X; T i | i ∈ Σ}.

A simple example of an IFS can be constructed for a representation of real numbers in [0, 1]. Let

T i : x → x + i 3 : [0, 1] → [0, 1],
where i ∈ Σ = {0, 1, 2}. The IFS {[0, 1]; T 0 , T 1 , T 2 } can thus express the representation of any real number in [0, 1]. We are substantially interested in so-called hyperbolic IFS, whose transformations T i are all contractive.

Definition. A transformation T : X → X is called contractive if and only if there exists a real s, 0 s < 1 such that d(T (x), T (y)) < s • d(x, y) for all x, y ∈ X. The minimal coefficient s which satisfies the inequality is called the contraction coefficient of the transformation T . To distinguish contraction coefficients of different transformations we also use the function notation s(T ) to denote the contraction coefficient of the transformation T .

Definition. For the set of non-empty compact subsets of X, denoted H(X), we define the Hausdorff distance d X induced by the metric d: Since the metric space (X, d) is complete, the set (H(X), d X ) is also complete [START_REF] Barnsley | Fractals everywhere[END_REF].

d X (A, B) = max{sup
In the early 80's, Hutchinson [START_REF] Hutchinson | Fractals and self-similarity[END_REF] used the Banach fixed point theorem to deduce the existence and the uniqueness of an attractor for a hyperbolic IFS, i.e. the fixed point of the associated contractive map. It is thus possible to define an operator T : H(X) → H(X), called Hutchinson operator [START_REF] Barnsley | Fractals everywhere[END_REF], as the union of the IFS transformations T i :

T(K) = N -1 i=0 T i (K).
The contraction coefficient of the Hutchinson operator is thus defined by:

s max = max i=0,••• ,N -1 s(T i ).
Definition. If s max < 1 then T is also contractive in the complete metric space (H(X), d X ). According to the Banach fixed point theorem [START_REF] Barnsley | Fractals everywhere[END_REF], T has a unique fixed point A. This fixed point is named the attractor of the IFS:

A = T(A) = N -1 i=0 T i (A). (3.1)

Examples of an IFS

In the following examples we use the Euclidean plane R 2 with the associated Euclidean metric.

The first example is the Sierpinski triangle illustrated in figure 3.1. We define an IFS with the three contractive transformations T 0 , T 1 and T 2 , defined by the homotheties with centres p i and ratios 1 2 as follows: Here H p s denotes the homothety centred at point p with ratio s. The Sierpinski triangle is the attractor of this IFS, that is a non-empty compact set satisfying the equation (3.1):

T 0 = H p0 1 2 T 1 = H p1 1 2 T 2 = H p2
A = H p0 1 2 (A) ∪ H p1 1 2 (A) ∪ H p2 1 2 (A).
The shape of the IFS attractor depends entirely on the transformations T i . The Sierpinski triangle has a number of particularities, which can generally be found in other IFS attractors. One of the main characteristics, common to all IFS attractors, is self-similarity with changing scale. Inside the attractor there are subsets similar to the whole set, and this pattern occurs at any scale.

Figure 3.2 illustrates the same triangle with a slightly modified geometry. Here, IFS transformations have the same ratio, but we modify the centre of the homothety T 2 .

The following example demonstrates an IFS composed of three homotheties Figure 3.3: Attractor of an IFS composed of two homotheties with positive ratios 1 2 and one homothety with negative ratio -1 2 . The centres of these homotheties are p i for i = 0, . . . , 2.

whose centres form an equilateral triangle, but we change their ratios as follows:

T 0 = H p0 1 2 T 1 = H p1 1 2 T 2 = H p2 -1 2 ,
The attractor of this IFS is shown in figure 3.3. These examples show that a simple IFS system can produce complex selfsimilar structures, but it is quite difficult to control the geometry and topology of these objects, even though the geometry depends continuously on the transformations [START_REF] Barnsley | Fractals everywhere[END_REF][START_REF] Gentil | Les fractales en synthése d'images: le modèle IFS[END_REF].

IFS construction algorithms

The Banach fixed point theorem [START_REF] Barnsley | Fractals everywhere[END_REF] describes the existence and uniqueness of the attractor for a hyperbolic IFS {X; T i | i ∈ Σ}. In this section we describe methods to approximate and to visualize the attractor.

There are two canonical ways to evaluate the attractor of a dynamical system:

• By constructing a decreasing sequence of non-empty compact sets which shrink down to the attractor. The smaller the compact set in the sequence, the better the attractor approximation with respect to the Hausdorff distance. Here, the initial compact set has to cover the attractor. This method provides an approximation of the fractal exterior.

• By constructing an increasing sequence of non-empty compact sets whose union has the attractor as its closure. Here, we have to start with a compact set included in the attractor. In this case, the bigger the compact set in the sequence, the better the attractor approximation.

Other construction and visualisation methods exist [START_REF] Prusinkiewcz | Koch curves as attractors and repellers[END_REF][START_REF] Hepting | Rendering methods for iterated function systems[END_REF], such as the dynamic chaos game or non-deterministic algorithms [START_REF] Barnsley | Fractals everywhere[END_REF][START_REF] Barnsley | Iterated function systems and the global construction of fractals[END_REF]. However, we are substantially interested in the deterministic method allowing us to approximate most operations. In fact, in both these cases, we use a similar approximation algorithm, i.e. we recursively apply transformations to the initial non-empty compact subset.

The attractor of an IFS may therefore be evaluated recursively. That is, the attractor can be approximated by a sequence {K n } n∈N converging to A. The initial element in the sequence is defined by means of a primitive K ∈ H(X).

The following elements are defined recursively:

K 0 = K K n+1 = i∈Σ T i (K n ).
The elements K n are images of composite functions applied to K.

Each element in the sequence represents an approximation of the IFS attractor. Each term K n is composed of N n images of K by a composite of n functions. For example, a sequence of the attractor approximations for an IFS {X; T 0 , T 1 } is presented here:

K 0 = K, K 1 = T(K 0 ) = T 0 (K) ∪ T 1 (K), K 2 = T(K 1 ) = T 0 T 0 (K) ∪ T 0 T 1 (K) ∪ T 1 T 0 (K) ∪ T 1 T 1 (K), K 3 = T(K 2 ) = T 0 T 0 T 0 (K) ∪ T 0 T 0 T 1 (K) ∪ T 0 T 1 T 0 (K) ∪ T 0 T 1 T 1 (K) ∪ T 1 T 0 T 0 (K) ∪ T 1 T 0 T 1 (K) ∪ T 1 T 1 T 0 (K) ∪ T 1 T 1 T 1 (K), . . . . . . . . . K n = T(K 2 ) = αi∈Σ T α1 T αn (K n-1 ).
In this iterative algorithm, a set of transformed primitives K is constructed recursively and calculations can be represented by an evaluation tree. Each node on the i-th level of the tree corresponds to the image of a composite of i IFS transformations. This tree is traversed up to a given depth n, where we display the image of K by the composite function associated with the current node, as shown in figure 3.4.

Note that these composite functions are calculated from left to right. The primitive K is finally transformed by a constructed composite function. In practice, the IFS transformations T i are affine operators and can therefore be represented by matrices. A composite affine transformation can thus be represented by a product of transformation matrices.

In chapter 4 we use the notion of the evaluation tree to decompose an operation image and to optimize the approximation.

Controlled IFS

In IFS all the transformations are applied at each iteration. It is possible to enrich this model by adding rules to control the iterations. This is the principle of a CIFS (Controlled IFS).

CIFS are more general systems allowing us to control certain parts of the IFS attractor. A CIFS denotes an IFS with restrictions on transformation sequences imposed by a control graph. This system is similar to "Restricted IFS" (RIFS) [START_REF] Barnsley | Fractals everywhere[END_REF], and is also described [START_REF] Prusinkiewicz | Language-restricted iterated function systems, koch constructions, and l-systems[END_REF][START_REF] Thollot | Construction of fractales using formal languages and matrices of attractors[END_REF] by means of formal languages, called LRIFS (Language-Restricted Iterated Function System). CIFS allow us to handle fractals with subsets of different dimensions. CIFS define objects whose geometry can be complex. However CIFS attractors are more convenient and controllable than IFS attractors for manufacturability purposes.

The attractor of a CIFS can be evaluated by an automaton [MW88] defined on the control graph. Each state of the automaton corresponds to different parts of the modelled object. States are associated with construction spaces. Transitions between states indicate that one sub-part is contained in another one. It is then possible to control the attractor more precisely. However, it is usually NP-complex to evaluate, to analyse and to characterize (location, size) the resulting shape.

Definition of CIFS

Definition. A CIFS is given by an automaton, where each state q is associated with an attractor A q ∈ X q , and each transition from q to w is associated with an operator X w → X q . Here is a list of parameters describing the CIFS:

• An automaton (Σ, Q, δ), where Σ is an alphabet, Q is a set of states and δ is a transition function δ : Q × Σ → Q;

• A set of complete metric spaces associated with the automaton states {X q } q∈Q ;

• An operator associated to each transition T q i : X δ(q,i) → X q ; • A compact set K q ∈ H(X q ), called primitive, associated with each state q ∈ Q. Primitives are not used to define the attractor, but only to approximate it;

• Finally, the automaton is provided by an initial state, noted by .

In the following, we denote by Σ q the restriction of Σ by outgoing transitions from the state q, i.e.:

Σ q = {i ∈ Σ, δ(q, i) ∈ Q}

Attractor of a CIFS

The attractor of an IFS can be evaluated recursively.

Definition. A CIFS defines a family of attractors associated with the states: {A q } q∈Q , where A q ∈ H(X q ). The attractors A q are mutually defined recursively:

A q = i∈Σ q T q i (A δ(q,i) ) (3.2)
A CIFS attractor, denoted A, is the attractor associated with the initial state , i.e.: A = A .

As for IFS, each CIFS attractor can be approximated by a sequence {K q n } n∈N converging to A q . Each state q ∈ Q is associated with a primitive K q ∈ H(X q ), which defines the initial element in the sequence. The following elements are mutually defined recursively:

K q 0 = K q K q n+1 = i∈Σ q T q i (K δ(q,i) n )
In this iterative algorithm a set of transformed primitives K q is recursively constructed and the calculations can also be represented by an evaluation tree. Each node on the i-th level of the tree corresponds to the image of a composition of i CIFS transformations, i.e. a path of length i in the automaton. This tree is traversed up to a given depth n, where we display the image of K q by the composite function associated with the current node, as shown in figure 3.5.

Consider an example of the CIFS attractor, illustrated in figure 3.6. The system is described by an automaton with two states: a and b. There are two subdividing operators from the state a and three from b. Each transition δ(q, i) = w of the automaton is associated with an operator T q i : X w → X q . Assume that X a and X b are both the same Euclidean affine plane. Let the operator between these spaces be the identity map Id and the operators defined from the space into itself be the following homotheties:

T a 0 = H p3 0.42 T b 0 = H p0 0.5 T a 1 = Id T b 1 = H p1 0.5 T b 2 = H p2 0.5 ,
where H p s denotes the homothety centred at point p with ratio s. Thus, the attractors A a and A b satisfy the following equations:

A a = i∈Σ a T a i (A δ(a,i) ) = T a 0 (A a ) ∪ T a 1 (A b ) = H p3 0.42 (A a ) ∪ A b A b = i∈Σ b T b i (A δ(b,i) ) = T b 0 (A b ) ∪ T b 1 (A b ) ∪ T b 2 (A b ) = = H p0 0.5 (A b ) ∪ H p1 0.5 (A b ) ∪ H p2 0.5 (A b )
As can be deduced from defining A b (see equation 3.2), this attractor is the Sierpinski triangle with vertices p 0 , p 1 and p 2 .

In figure 3.6 we can see that the attractor A a is an infinite set of contracted Sierpinski triangles.

Remark Attractors of a CIFS are uniquely defined if operators associated with each cycle in the control graph are contractive. In the previous example, the transition between the states does not need to be contractive.

Projection onto the modelling space

From IFS and CIFS, various shapes can be modelled. However, it is difficult to control their topological properties. The shapes are determined by a set of geometry operators. Modifying these operators leads to both global and local changes in the shape and affects not only geometry but also topology.

To separate global and local control, we construct the attractor in a barycentric space and then project this attractor according to a set of control points in the modelling space. The global shape can therefore be modified by moving control points, while the local geometry (geometric texture) is defined by subdividing operators, i.e. CIFS transformations [START_REF] Zair | Fractal modeling using free form techniques[END_REF][START_REF] Schaefer | Subdivision schemes and attractors[END_REF].

Figure 3.8: The automaton of the CIFS considered in the example. There are three subdividing operators from state S to itself and one projecting operator from S to .

Consider an example of the Sierpinski triangle as an attractor of a projecting CIFS. The system is described by an automaton with two states: (for the modelling space) and S (to construct the triangle). There are three subdividing operators from state S to itself and one projecting operator from S to . Figure 3.8 shows the automaton of this CIFS. Transition operators are the following:

T 0 = P T S 0 = H (1,0,0) 0.5 T S 1 = H (0,1,0) 0.5
T S 2 = H (0,0,1) 0.5 , Figure 3.9: Attractor of a CIFS described by an automaton with three states , C and S. where P is the projecting operator and H p s denotes the homothety centred at point p with ratio s.

In this case, the Sierpinski triangle is constructed in a 3D barycentric space. We choose the centres of homotheties as the unit vectors to simplify the following projection onto the modelling space. The projecting operator P has the following matrix representation:

P =      p 0 [0] p 1 [0] p 2 [0] p 0 [1] p 1 [1] p 2 [1] . . . . . . . . . p 0 [d] p 1 [d] p 2 [d]     
, where d is the dimension of the modelling space and {p 0 , p 1 , p 2 } is a set of control points, called the control polygon. Here p i [j] denotes the j-th coordinate of point p i . The projecting matrix is thus composed of the control point coordinates in columns.

According to definition (3.2), the attractor A of this CIFS satisfies the following equation:

A = P (A S ).
Note that applying affine transformations becomes simple. In order to apply a transformation R to the attractor, we modify the projecting operator as follows:

P ← R • P,
where P is the initial projecting operator and P is the new one. Indeed, the attractor of the modified CIFS, denoted by A , can be determined as follows:

R(A ) = RP (A S ) = P (A S ) = A .

Consider another example of a CIFS, illustrated in figure 3.9. The system is described by an automaton with three states: , C and S. State C is to construct a square, and is associated with the 4D barycentric space X C , illustrated in figure 3.10. The square is then projected to the modelling space X . State S is to construct the Sierpinski triangle in the 3D barycentric space X S , which is then projected onto X C . Figure 3.11: The automaton of the CIFS considered in the example. There are four subdividing operators from C to construct a square, but the Sierpinski triangle is inserted as one quarter. There are three subdividing operators from S to construct the Sierpinski triangle. The result is then projected onto the modelling space X .

There are four subdividing operators from state C and three from S.

δ(C, 0) = C δ(S, 0) = S δ(C, 1) = C δ(S, 1) = S δ(C, 2) = C δ(S, 2) = S δ(C, 3) = S
Each transition δ(q, i) = w of the automaton is associated with an operator T q i : X w → X q . The subdividing operators are defined as follows:

T C 0 = H
(1,0,0,0) 0.5

T S 0 = H (1,0,0) 0.5 T C 1 = H (0,1,0,0) 0.5 T S 1 = H (0,1,0) 0.5 T C 2 = H (0,0,1,0) 0.5 T S 2 = H (0,0,1) 0.5
, where H p s denotes the homothety centred at point p with ratio s. The projecting map T C 3 of the Sierpinski triangle onto the space X C is defined as follows:

T C 3 =     0 0 0 0 0.5 0 0.5 0 0 0.5 0.5 1     .
This matrix projects the 3D barycentric space X S , in which the triangle is constructed, onto the 4D barycentric space X C by conforming the unit vectors of X S to the control points (0, 0, 0.5, 0.5), (0, 0.5, 0, 0.5) and (0, 0, 0, 1), as illustrated in figure 3.10. Finally, the projecting map P is defined as follows:

P =      p 0 [0] p 1 [0] p 2 [0] p 3 [0] p 0 [1] p 1 [1] p 2 [1] p 3 [1] . . . . . . . . . . . . p 0 [d] p 1 [d] p 2 [d] p 3 [d]     
, where d is the dimension of the modelling space and {p i | i = 0, . . . , 3} is the control polygon. According to definition (3.2), the attractors of this CIFS satisfy the following equations:

A = P (A C ), A C = T C 0 (A C ) ∪ T C 1 (A C ) ∪ T C 2 (A C ) ∪ T C 3 (A S ), A S = T S 0 (A S ) ∪ T S 1 (A S ) ∪ T S 2 (A S ).
Figure 3.12 illustrates the same attractor projected onto a different control polygon.

BCIFS

From IFS and CIFS, various shapes can be modelled. However, it is difficult to control their topological properties. The shapes are determined by a set of geometry operators. Modifying these operators leads to both global and local changes in the shape and affects not only geometry but also topology.

In order to control the topological structure of the modelled shape, one can enrich CIFS by integrating the topological model to obtain BCIFS (Boundary Controlled Iterated Function System) [TBSG + 06]. It is thus possible to specify incidence and adjacency constraints [START_REF] Tosan | Surfaces fractales définies par leurs bords[END_REF][START_REF] Gouaty | Modélisation géométrique itérative sous contraintes[END_REF] for subdivision and thereby to control the resulting topology: classic topology (curve, surface,. . . ) or fractal topology. The topological model used is more general than the classic one. A topological cell may be fractal. For example, a face may be the Sierpinski triangle or an edge may be the Cantor set, but the topological structure remains consistent.

Definition of BCIFS

The topological model is commonly used in CAD systems to represent geometric objects. This formalism was introduced to separate the topology and the geometric properties of a considered object. The topological structure is encoded by a set of topological cells (faces, edges, vertices) interconnected by a set of incidence and adjacency relations. The incidence relations are based on nesting of cells: each face is bounded by a set of edges, and each edge is bounded by two vertices. The adjacency relations are based on sharing cells: two adjacent faces share a common edge, and two adjacent edges are bounded by a common vertex. BCIFS is thus an extension of a CIFS enriched by the topological model.

Classic shape description includes a finite number of components. Each shape is represented by a finite set of cells with incidence and adjacency relations using polynomial or rational (NURBS) functions. Topological structure of these shapes is thus encoded by a finite graph. An iterative description of fractal shapes thus consists of an infinite number of components generated by subdivision. At each step, the algorithm generates an approximate mesh with an increasing (potentially infinite) number of components. The topological struc- ture is encoded by a finite set of rules, according to which the current level of approximation is obtained from the previous one. All topological rules on the n-th subdivision level can thus be derived from the rules on the first level.

Each topological cell corresponds to an attractor in a certain space. A face corresponds to the attractor

A f of the IFS (X f ; T f 0 , T f 1 , T f 2 , T f 3 ) (see figure 3.13
). An edge corresponds to the attractor A e of the IFS (X e ; T e 0 , T e 1 ) (see figure 3.14). A vertex corresponds to the attractor A v of the IFS (X v ; T v 0 ). Figure 3.15 illustrates how a square can be decomposed into the topological cells. Incidence operators are used to generalize the concept of edge. Each included cell is defined by an embedding operator. Each face of the square is bounded by four edges and each edge is bounded by two vertices. Embedding of edges A e in X f is defined by four embedding operators f i : X e → X f for i = 0, . . . , 3:

f i A e ⊂ A f .
Embedding of vertices A v in X e is defined by two embedding operators e i : X v → X e for i = 0, 1:

e i A v ⊂ A e .
There are two types of transition in the BCIFS automaton:

• transitions subdividing a topological cell;

• transitions embedding a topological cell in another one.

The alphabet Σ is also divided into:

• symbols of subdivision Σ ÷ = {÷ i | i ∈ N}; • symbols of incidence Σ ∂ = {∂ i | i ∈ N}.
Each subdividing transition δ(q, ÷ i ) = w is associated with a subdividing operator T q i : X q → X w , where q, w ∈ Q and ÷ i ∈ Σ ÷ . Similarly, each embedding transition δ(q, ∂ i ) = w is associated with an embedding operator q i : X q → X w , where q, w ∈ Q and ∂ i ∈ Σ ∂ . Usually edges share certain vertices. An example of shared vertices is illustrated in figure 3.15. These relations can be translated into adjacency constraints on the embedding operators f i , i = 0, . . . , 3:

f 0 e 0 = f 3 e 1 f 0 e 1 = f 1 e 1 f 1 e 0 = f 2 e 1 f 2 e 0 = f 3 e 0 ,
There is a correspondence between subdivided borders of a face and the borders of subdivided faces, which can be expressed by the following incidence constraints:

T f 0 f 0 = f 0 T e 0 T f 1 f 0 = f 0 T e 1 T f 1 f 1 = f 1 T e 1 T f 3 f 1 = f 1 T e 0 T f 3 f 2 = f 2 T e 1 T f 2 f 2 = f 2 T e 0 T f 2 f 3 = f 3 T e 0 T f 0 f 3 = f 3 T e
1 , One can prove that these constraints imply coherence between the borders on all levels of subdivision.

To ensure the continuity of the shape, subdivided faces can also share edges. These relations can be translated into adjacency constraints between borders of the subdivided faces:

T f 0 f 1 = T f 1 f 3 T f 1 f 2 = T f 3 f 0 T f 3 f 3 = T f 2 f 1 T f 2 f 0 = T f 0 f 2 ,
All the incidence and adjacency constraints are thus illustrated in figure 3.16. In this case, there are twelve constraints that guarantee the absence of holes in the surface. These constraints define a system of constraints on subdividing and embedding operators and establish an equivalence of paths in the automaton for the first level of subdivision. The equivalence on the next levels of subdivision can thus be derived from this system of equations.

Objects modelled using BCIFS do not necessarily have specific fractal properties; shapes with a classic topology can be modelled as well. The topological constraints of BCIFS guarantee the topological structure of the limit shape. However, by choosing appropriate initial compact sets, these constraints guarantee the coherence and integrity of the topological structure for all the approximation levels as well. Moreover, the topological constraints are easy to solve regardless of the nature of the embedding operators [START_REF] Gouaty | Modélisation géométrique itérative sous contraintes[END_REF].

For simplification purposes, in this thesis, we primarily describe the ideas for an IFS, keeping in mind that generalization to a CIFS is immediate.

Manipulations with iterative models

Despite the difficulty of controlling the final shape defined by an iterative system there exist construction algorithms and methods to manipulate their attractors. In this section we describe some existing algorithms.

Model reconstruction

Much work has been done on the so-called "inverse problem" for IFS and its generalisations. This problem involves developing methods to find a set of contractive mappings whose attractor is close to a given shape in the sense of a predefined error measure.

Barnsley [START_REF] Barnsley | Fractals everywhere[END_REF] proposed a method to determine the transformations associated with a given attractor, called the Collage Theorem. The theorem states that to find an IFS whose attractor is close to a given set, one must determine a set of contraction mappings on a suitable space, within which the given set lies. The union of the images of the given set under these transformations must be close to the given set with respect to the Hausdorff metric. More precisely, the theorem implies that finding an IFS {X; T i | i ∈ Σ} whose attractor is close to a given image K is equivalent to the minimization of the distance

d X (K, i∈Σ T i (K)).
Roll, Lutton et al. [RLG + 95] developed a method to solve the "general" inverse problem for IFS within a reasonable computation time by using variable-sized structures in the genetic programming algorithm, which seems to perform a more efficient search in a large space.

Bounding the attractor

Methods to approximate the convex hull of an IFS attractor have already been developed.

Strichartz, Wang, Kenyon et al. [START_REF] Strichartz | Geometry of self-affine tiles i. Information science research[END_REF][START_REF] Kenyon | Geometry of self-affine tiles ii[END_REF] studied the boundaries and the convex hulls of self-affine tiles that can be considered as attractors of very special affine IFS, for which all the transformations have the same linear part.

Lawlor and Hart [START_REF] Lawlor | Bounding recursive procedural models using convex optimization[END_REF] presented an algorithm to construct a tight bounding polyhedron of the IFS attractor. An algorithm expresses the IFS-bounding as a set of linear constraints on a linear objective function, which can then be solved via standard techniques for linear convex optimization. This method works for a predetermined number of convex hull faces and shows the interactive rate only when this number is small.

More recently, Duda [START_REF] Duda | Analysis of the convex hull of the attractor of an ifs[END_REF] and Martyn [START_REF] Martyn | The attractor-wrapping approach to approximating convex hulls of 2d affine ifs attractors[END_REF] have presented methods to calculate the approximate convex hull of an affine IFS attractor. The two methods are similar. The first is based on so-called "width function" that maps a point to the nearest bounding half-space in a given direction. The second is based on constructing a sequence of balls that bound corresponding parts of the attractor to approximate the convex hull. This latter approach is presented in 2D only.

Another well studied problem involves determining a bounding ball for the attractor of an IFS.

Gentil [START_REF] Gentil | Les fractales en synthése d'images: le modèle IFS[END_REF] described an approach based on the dichotomous search for the minimal radius of a ball bounding the IFS attractor. The approach can be applied in a multi-dimensional space and the result is calculated for a given precision.

Hart and DeFanti [START_REF] Hart | Efficient anti-aliased rendering of 3d linear fractals[END_REF] introduced a method which starts with a unit ball centred at the origin. The algorithm iteratively produces a sequence of balls converging to the limit ball that bounds the attractor. Rice [START_REF] Rice | Spatial bounding of self-affine iterated function system attractor sets[END_REF] improved the approach of Hart and DeFanti by optimizing the radius of the bounding ball using the generic optimization package. He also showed that the centre of the limit ball can be determined analytically by solving a system of linear equations. Martyn [START_REF] Martyn | The attractor-wrapping approach to approximating convex hulls of 2d affine ifs attractors[END_REF] showed that the solution of this system is the centroid of the attractor with particular weights. To obtain a better approximation, he presented an heuristic iterative method, called "balancing the attractor". The algorithm is not limited by the dimension of the space in which the attractor lies.

More recently, Martyn [START_REF] Martyn | The smallest enclosing disc of an affine ifs fractal[END_REF] presented a novel approach to approximate the smallest disc to enclose an affine IFS attractor at any accuracy. The method is based on a concept of spanning points he introduced to describe the extent of the IFS attractor.

We studied the convex hull approximation for a given affine IFS attractor [START_REF] Mishkinis | Approximate convex hull of affine iterated function system attractors[END_REF]. Our work can be considered as a generalization and an optimization of Martyn's method. The algorithm constructs a sequence of convex hull approximations using the self-similarity property of the attractor in order to reduce the number of necessary operations. In addition, we have introduced a method to simplify the approximate convex hull without losing accuracy.

Combination approaches

Thollot [START_REF] Thollot | Extension du modèle IFS pour une géométrie fractale constructive[END_REF] proposed a model generalising IFS based on Languages Theory. This model permits us to compose two attractors by means of specific language operations, such as union, condensation, mixing and intersection. The author studied the relation between attractors defined by languages L and L and also the attractor L L defined by combining these languages. Certain operations were introduced and studied. The attractor defined by condensation of these languages can be interpreted as the attractor associated with the first language fabricated by iterating the attractor of the second one. Mixing is the operation that "fills the space" between two initial attractors. The language intersection of two IFS is equivalent to the intersection of their sets of transformations. Note that the language intersection is not equivalent to a classic set intersection. In fact two different languages may have the same attractor and the result thus depends on the choice of languages. In addition, this study presented an algorithm to construct the result of these operations in the case of rational languages.

Zaïr and Tosan [ZT96], Thollot et al. [START_REF] Thollot | Modeling fractal shapes using generalizations of IFS techniques[END_REF] presented an approach to fractal modelling which is based on IFS theory. The authors used free-form techniques to provide a practical and efficient way to build controlled fractal attractors. They also showed that it is possible to manipulate fractals with methods used for classic geometric design, such as control point editing, deformations or tensor products.

Martyn [START_REF] Martyn | A new approach to morphing 2d affine ifs fractals[END_REF] studied the affine stable morphing of 2D affine IFS fractals. The affine stability of metamorphosis can be described by the two following properties:

• intermediate shapes should preserve the topology under affine deformations of the key fractals;

• intermediate shapes should be at least as similar to the key fractals as they are to each other.

The approach presented in the paper satisfies these conditions. The author also considered morphing fractals specified by IFS with different numbers of transformations and gave a number of solutions. Recently, Podkorytov et al. [START_REF] Podkorytov | Geometry control of the junction between two fractal curves[END_REF] developed tools based on BCIFS to construct and to control a junction between two fractal curves. The authors also deduced conditions that guarantee continuity of the intermediate curve by studying the eigenvalues of the subdividing operators. Moreover, the authors suggested a method to control the differential behaviour at the connection points between the initial curves and the intermediate one.

Gentil and Neveu [START_REF] Gentil | Mixed-aspect fractal surfaces[END_REF] developed a formalism based on IFS that builds parametric shapes (curves, surfaces, . . . ) with a non-uniform local aspect: to every point is assigned a "geometric texture" that evolves continuously from a smooth to a rough aspect. The principle is to blend shapes with uniform aspects to define the shape with a variable aspect. A blending function controls the influence of each initial shape.

Chapter 4

Calculation of CAD operations on fractals

Our goal is to identify the different cases that arise when we want to evaluate CAD operations on fractals and to study the approximation of these operations. In section 2.2, we listed the main standard CAD operations presented in most CAD systems and identified the properties on which the operations are based. A preliminary analysis shows that the operations are based on four basic properties: affine invariance, topological structure, parametrization and differential properties. These properties are therefore essential and we must determine their equivalents in fractal geometry.

In this chapter, we classify these operations by adapting them to fractals in three categories. The first category consists of operators for which the result can be formally represented. The second is composed of operators that require extending or adapting to fractals because of their specific properties. And the third category corresponds to cases for which we can approximate the result.

We then discuss in detail the last category, that in which we can approximate the operation. In studying and identifying the constraints on operations to evaluate this result, we introduce several algorithms. The approximation algorithms are generic, i.e. defined for an arbitrary operation satisfying certain constraints that we state. These algorithms can be applied once we implement the required set of operators (interface).

The proof of convergence of these algorithms is based on Domain theory and computability of functions, as presented by Edalat in solid modelling in [START_REF] Edalat | Foundation of a computable solid modelling[END_REF]. In studying measure theory and IFS fractals [Eda95], A. Edalat described a computational model for dynamical systems using the upper space and the probabilistic power domains. More particularly, he developed a constructive framework to study IFS attractors and he described when the attractor of a dynamical system is computable.

Classification in the context of fractals

Fractals have more complex and more subtle properties than classic shapes, and existing CAD systems are not suitable to manipulate fractal structures. In this section, we study to what extent the standard CAD operators can be defined and evaluated on fractals. We are primarily interested in shapes represented by BCIFS described in section 3.3.

By analysing the graph of operations presented in section 2.2, we identify the basic properties on which the CAD operations are based. These properties are thus required for their realisation.

A preliminary analysis shows that CAD operations are based on four basic properties: affine invariance, topology, parametrization and differential properties. For shapes described by a BCIFS, the affine invariance property is immediately verified, because objects are defined using control points. In order to apply an affine transformation we can simply modify the control points in the same way as for NURBS surfaces. The topology is formalized by the model and extends concepts of classic topology to fractal topology [Gou09,GTWS10,GHT11]. Parametrization is obtained in terms of a BCIFS address function [START_REF] Barnsley | Fractals everywhere[END_REF]. Finally, differential properties have been rarely studied [Ben09, BGN08] and, for example, the notion of Frenet basis has still not been adapted for fractal curves.

These basic properties are thus partially defined and it is then possible to adapt CAD operations and algorithms to fractals. Here, we distinguish three cases:

1. The operation result has a formal expression; 2. The operation requires a generalization because of the specific properties of fractals;

3. The operation result cannot be represented formally at the moment, but it can be approximated.

The first category is the same as in the classification described in section 2.2. An operation Op preserves the BCIFS of an object S, and the output BCIFS has T (S) as its attractor. The first example is consists of applying an affine transformation. As for NURBS surfaces, we simply need to apply the transformation to control points, hence the BCIFS attractor is projected according to the set of transformed control points (see section 3.

2.3).

There are certain other examples mentioned in section 3.4, such as tensor products of fractal curves [START_REF] Zair | Fractal modeling using free form techniques[END_REF], specific language operations [START_REF] Thollot | Extension du modèle IFS pour une géométrie fractale constructive[END_REF] (union, condensation, mixing, language intersection), some studies on blending between two fractal curves [START_REF] Podkorytov | Geometry control of the junction between two fractal curves[END_REF].

The second category corresponds to operations that have not yet been defined for fractals. For example, "pulling curve onto surface" cannot be used in the original form. This operation projects a given curve onto a given surface along the normal to the surface. Since fractal surfaces are nowhere differentiable, this projection would be discontinuous. To obtain a more intuitive result, it is therefore necessary to redefine or to generalise the operation.

In fact, most of operations in this category depend on the differential properties of shapes. Because of the lack of continuity, applying this kind of operation becomes inconvenient. Examples are "network", determining tangents, or the extrusion by the normal direction. Nevertheless, some of these operations may be generalised to a more convenient analogue. Inversely, in some cases the operations can be restricted to a set of suitable cases, in which they may be applied. For example, we can require a surface to have a tangent plane at any point of an edge to apply a "patch". Furthermore, some operations, such as the curvature graph or calculating curve length, are based on properties that are not defined for fractals. For instance, consider the attractor of the von Koch curve. Starting with a single line segment, the von Koch curve can be constructed by recursively altering each line segment as follows (see figure 4.1):

1. divide the line segment into three segments of equal length; 2. draw an equilateral triangle that has the middle segment obtained from step 1 as its base and points outward;

3. remove the line segment that is the base of the triangle from step 2.

The number of sides N n on the n-th iteration is therefore:

N n = N n-1 • 4 = 4 n ,
where N 0 = 1. After each iteration, the number of sides increases four times, because we add two segments in the middle of each segment. Let l be the length of the initial segment. The length of each segment L n on the n-th iteration is thus:

L n = L n-1 3 = l 3 n
After each iteration, the length of each segment is three times smaller. The length of the von Koch curve after n iterations can therefore be determined by the following equation:

P n = N n • L n = l • 4 3 n ----→ n→∞ ∞.
The von Koch curve itself is of infinite length. Operations such as the chamfer operation or dividing a curve by its length cannot be applied directly. These operations must either be restricted to fractal curves with a finite length, or redefined in a more convenient way so that they can be applied to curves with an infinite length (the same for surface and volume).

Certain new operations that are natural for fractals can also be included in this category, such as fractal dimensions, for example. After all, it is possible to apply operations from this category to an attractor approximation as a classic object. Approximating the attractor is described in detail in chapter 3.3. The time complexity, in this case, depends exponentially on the number of iterations.

For the last category, we develop generic algorithms to approximate the operator result. These algorithms are defined for any operator satisfying certain explicit properties. We show that with some additional constraints these algorithms can be optimized by exploiting the self-similarity property.

Approximation of a CAD operator image

In this section, we formalise application of an arbitrary CAD operation and we establish the notation used in this thesis. To simplify explanations we describe the idea on IFS; a generalization to CIFS is thus immediate.

To denote a CAD operation we use the symbol Op and the function notation. Note that the operator Op may be any algorithm that takes a non-empty compact subset B in a complete metric space X as input, and generates the resulting set Op(B). We then describe the algorithm to approximate the image Op(A), where A is the attractor of our iterative system, and we show how to optimize the algorithm when the operator verifies certain specific properties.

Definition. Given a hyperbolic IFS {T

i } N -1
i=0 defined in a complete metric space (X, d X ) with the attractor A. Let Σ = {0, . . . , N -1} be the set of IFS transformation indices. An arbitrary operator Op could be any algorithm that takes a non-empty compact subset in a complete metric space X as input, i.e.:

Op : H(X) → Y,
where Y denotes the target space of the operator. This space may be a metric space, as for most CAD operations. However, it may also be a non-metric space, as for the membership predicate or any other classifying operation.

In studying measure theory and IFS [START_REF] Edalat | Dynamical systems, measures, and fractals via domain theory[END_REF], A. Edalat developed a constructive framework to study IFS attractors. He determined when the attractor of a dynamical system is computable. As mentioned in section 3.1.3, there are two canonical ways to evaluate the attractor of a dynamical system:

• By constructing a decreasing sequence of non-empty compact sets which shrinks down to the attractor.

• By constructing an increasing sequence of non-empty compact sets whose union has the attractor as its closure.

In fact, in both cases we use similar approximation algorithms, i.e. we recursively apply transformations to an initial non-empty compact subset.

Starting with an arbitrary compact subset B ∈ H(X), we construct a sequence of compact sets converging to the attractor A:

(B n ) n∈N = T n (B),
where T is the Hutchinson operator corresponding to a given IFS. Statement. We denote r = d X (A, B). The Hutchinson operator is contractive, we therefore have:

d X (A, B n ) s(T) n • r ----→ n→∞ 0,
where s(T) < 1 is the contraction coefficient of the Hutchinson operator and n is the number of iterations.

Therefore, ∀ε > 0 ∃m ∈ N : ∀n m d X (A, B n ) < ε, or equivalently,

B n ----→ n→∞ A.
The Hutchinson operator is thus continuous at the attractor A.

The idea of our generic algorithm is to apply the operation to attractor approximations, i.e. to calculate Op(B n ), in order to obtain a sequence of approximate results. If the operator Op is continuous at the point A, then this sequence converges to the required result Op(A).

For instance, we consider the operator determining the distance from a given point to the attractor. Let p ∈ X be a point; the operator Op is then defined as follows:

Op p (A) : A → min In this case, the space Y is the set of real numbers (R). One can show that this operator is continuous at point A and, moreover, that it is a Lipschitz function.

We can therefore construct the attractor approximations B n and compute the images Op(B n ).

Statement. For a given accuracy ε, the number of iterations is determined as follows:

n = log( ε Lr ) log(s(T)) , (4.1)
where r = d X (A, B), s(T) < 1 is the contraction coefficient of the Hutchinson operator and L is the Lipschitz constant of the operation (L = 1 in this example).

Continuity of the generic algorithm

Definition. An operator Op is called continuous at point A iff for any sequence (B n ) n∈N of non-empty compact subsets of X converging to attractor A, the corresponding sequence (Op(B n )) n∈N converges to Op(A).

In general, most CAD operations are not continuous. However, in our case we construct only particular sequences defined by applying the Hutchinson operator T and convergence for any sequence is unnecessary. This issue is closely related to computability in solid modelling [START_REF] Edalat | Foundation of a computable solid modelling[END_REF] and more particularly to the question of domain choice.

For example, such basic predicates and operations as membership, subset inclusion and intersection are not continuous and, therefore, are not computable. However, in some cases we can restrict the domain of the operator Op (that is, the choice of the sequences (B n ) n∈N ) to guarantee the continuity at the attractor A.

Directed complete and continuous posets

In this section we introduce the basic notation of domain theory and we explain how continuous domains can be effectively presented to guarantee the convergence of our algorithms.

Definition. A set P with a binary relation is called a partially ordered set (poset), if the following conditions are satisfied ∀x, y, z ∈ P : There is a natural topology on dcpos given by Scott. Definition. An operator F : P → P from a dcpo P to another dcpo P is monotone if for any directed set D ⊆ P and for any pair of elements x, y ∈ D,

• x x (Reflexivity) • x y & y z ⇒ x z (Transitivity) • x y & y x ⇒ x = y (Antisymmetry)
x P y ⇒ F (x) P F (y).

Definition. An operator F is continuous iff F is monotone and it preserves a least upper bound of directed sets, i.e. The key idea of transitioning to domain theory is that functions may be discontinuous with respect to the Hausdorff metric, but continuous in the less restrictive Scott topology.

Consider the following example. Let f : R → R be a real-valued function defined by:

f (x) = x + 1 if x > a, x otherwise.
The plot of this function is shown in figure 4.5. In the standard topology this function is not continuous. As illustrated in figure 4.6, in the standard topology, a neighbourhood of the point a is bounded on both sides. By considering the pre-image of this neighbourhood, we do not obtain an open set. Thus, there exists a sequence x n → a whose tail does not fall entirely in the neighbourhood.

However, this function is lower semi-continuous. Note that a function F : X → R for a topological space X is lower semi-continuous iff the set {x ∈ X | F (x) > r} is open for all r ∈ R. This function is thus continuous with respect to the Scott topology on R. Indeed, a neighbourhood of the point a is only bounded on the bottom. As illustrated in figure 4.7, the pre-image of this neighbourhood is always open. Thus, we remove from consideration the sequences converging on the right side.

This simple example demonstrates that convergence in the Scott topology is less restrictive than in the classic topology. In general, most CAD operations are not continuous. Continuity in the Scott topology requires that a smaller set of sequences converge. Moreover, when applying the operation to fractals we construct only a single sequence defined by applying an iterative algorithm and convergence for all sequences is unnecessary.

We will now define the so-called order of approximation, which is also more suggestively called the way-below relation.

Definition. An element x ∈ P is way-below y ∈ P , or equivalently x approximates y, denoted by x y, iff whenever y D for a directed set D with a lub, there exists z ∈ D such that x z. This condition implies that x y, since the singleton set {y} is directed.

Consider the interval domain with the subset inclusion order:

[a, b] [c, d] ⇔ [a, b] ⊇ [c, d],
where However, being way-below of an element is a relative notion and does not reveal much about the element alone. For example, one would like to characterize finite sets in an order-theoretic way, but even infinite sets can be way-below any other set. The special property of these finite elements is that they are way below themselves, i.e. x

x. An element with this property is also called compact. Yet, such elements need not be finite nor compact in any other mathematical usage of the terms. Compact elements of a domain have the important special property that they cannot be obtained as a limit of a directed set in which these elements did not already occur.

The previous thoughts raise another question: is it possible to guarantee that any element of a domain can be obtained as a limit of much simpler elements? This is quite relevant in practice, since we cannot compute infinite objects but we may still hope to approximate these objects arbitrarily closely. More generally, we would like to restrict to a certain subset of elements as being sufficient to obtain all other elements as least upper bounds.

Definition. A basis of a poset P is a subset B ⊆ P , such that, for each x ∈ P , the set of elements in B that are way-below x contains a directed set with lub x. We say P is continuous if it has a basis. The poset is ω-continuous if it has a countable basis.

A dcpo is pointed if it has a least element. We can always add a bottom element to a domain to make it pointed.

Solid domain

In this section we introduce a domain chosen to represent the attractor of iterative systems and to extend the operator in a continuous and computable way.

Definition. The solid domain of a topological space X is the set of ordered pairs (A, B) of disjoint open subsets of X endowed with the information order:

(A 1 , B 1 ) (A 2 , B 2 ) ⇔ A 1 ⊆ A 2 and B 1 ⊆ B 2 .
An element (A, B) is called a partial solid: A and B are intended to capture, respectively, the interior and the exterior of a solid, possibly, at some finite stage of computation (see figure 4.8). The border of the solid is thus somewhere in between. The solid domain is a mathematical model to represent rigid solids. A. Edalat and A. Lieutier have shown [START_REF] Edalat | Foundation of a computable solid modelling[END_REF] that the basic predicates and operations are continuous and computable in the solid domain (SX, ). We can thus represent the attractor in the solid domain and then apply continuous operations to it.

A. Edalat presented a domain-theoretic model for iterated function systems [START_REF] Edalat | Domains for computation in mathematics, physics and exact real arithmetic[END_REF]. He constructed computational models for classic spaces using continuous dcpo and provided applications in IFS theory. A domain-theoretic model for IFS gives a unifying framework to study various properties and provides a set of new results in the theory and applications of these systems. In particular, this domain-theoretic approach allows us to prove convergence of the developing algorithms and to approximate operations on fractals.

In the study of computable partial solids, A. Lieutier [Lie99] introduced conditions on the domain to be an ω-continuous dcpo and showed when a continuous function f : D → E on the ω-continuous domains D and E is computable.

Definitions.

• A function is called recursive if the program computing this function ends within a finite number of steps for each input.

• A function is called partial recursive iff for some values given as inputs, a program computing this function may run forever.

• Computable functions are partial functions that can be computed by the Turing machine.

Definitions. We use the equivalent notion of partial recursive functions as a map that can be computed by a general purpose computer or a (finite length) program written in a general purpose language. For a countable set X:

• a subset O ⊆ X is called recursive, if
, n ∈ N, m f (n) ⇒ d X (O, O m ) 2 -n .
Note that, in defining the computability of O ⊆ X, the existence alone of a computable sequence of approximations converging to O is not sufficient; the convergence has to be effective.

If X is a second countable locally compact Hausdorff space, the solid domain is thus ω-continuous bounded complete dcpo's. An ω-continuous domain can be effectively presented with respect to the enumeration of a basis by requiring that the way-below relation restricted to the basis elements be recursively enumerable. This requirement allows us to compute the number of iterations n needed to achieve accuracy ε.

Attractor in the solid domain

Generally speaking, most CAD operations are not continuous and therefore not computable. Figure 4.9 illustrates the overall idea of our method. We construct a representation of the attractor in the solid domain and then apply a continuous extended operation to it. We use the symbol SOp to denote the continuous extension of the operator Op to the solid domain:

SOp : SX → SY.
For any subset K of a topological space, ∂K, K, K o and K c denote, respectively, the border, the closure, the interior and the complement of K.

To represent the attractor A in the solid domain we consider the partial solid (A o , A c ). However, constructing the attractor interior is quite complex, and we therefore often use the empty set as an approximate interior, i.e. (∅, A c ), as illustrated in figure 4.10. This partial solid is not the maximal element and, generally, is not a classic solid. However, to apply for example the boolean operations or the Minkowski sum, it is sufficient to consider the exterior only. Note that for such basic operations as membership predicate, subset inclusion or difference between two fractals, one must also find a non-empty approximation of the attractor interior. (∅,A c ). We use the empty set as an interior approximation.

In fact, when the interior of a partial solid is empty, this means that we do not know the attractor interior A o . The operation is therefore computed using only the exterior approximations. The membership predicate, for example, is semi-decidable in this case, because we can say if a point is outside the attractor, but cannot determine if the point is inside.

Statement. When avoiding the approximation of the attractor interior, we do not distinguish the objects (A 1 , B) and (A 2 , B) if the following equality holds:

( Ā1 ) o = ( Ā2 ) o .
The objects illustrated in figure 4.11 will not be distinguished without approximating the interior. However, this is never the case for IFS (or CIFS) attractors. These models are defined on a set of non-empty compact subsets of a complete metric space X. The attractors are therefore compact sets and the attractors with the same exterior cannot have different interiors.

Figure 4.11: Two partial solids which differ on a set of a zero measure. These partial solids will not be distinguished without approximating the interior.

Algorithm convergence

In this section we discuss convergence of constructed partial solids to the attractor of a given iterative system. We identify necessary conditions for the initial partial solid to constitute a directed set with the attractor as the lub. Here we assume the operator SOp to be continuous in the Scott topology.

Theorem. Let A n be a sequence of open sets, B n be a sequence of non-empty compact sets and SA = (A o , A c ) be the attractor of our iterative system in the Scott topology. Let SOp be a continuous operator (for example, a boolean operator or the Minkowski sum [START_REF] Edalat | Foundation of a computable solid modelling[END_REF]). If the partial solids SB n = (A n , B c n ) constitute a directed set with the lub SA then we can write the following equation:

n SOp((A n , B c n )) = SOp( n (A n , B c n )) = SOp((A o , A c )).
The proof is based on defining continuous functions in the solid domain.

The simplest and most commonly used example is when this sequence SB n constitutes an increasing chain, where we thus improve the information about the attractor with each iteration. The sequence SOp(SB n ) thus converges to the image of the attractor SA and this sequence can be computed with respect to the effective structure on SX, as described in section 4.3.2. Note that the image under the Hutchinson operator T(SB) gives (T(A), T(B) c ). We can then apply the IFS transformation to partial solids in the same manner.

In fact, constructing SOp((A n , B c n )) is equivalent to constructing Op(A n ) and Op(B n ). In this case, we use very particular sequences of approximations and we therefore converge to Op(A).

To exploit the solid domain, we should firstly verify that the sequence of the partial solids (A n , B c n ) constitutes a directed set and, more precisely, an increasing chain. We need to verify the following inequality:

(A 0 , B c 0 ) (A 1 , B c 1 ) . . . (A n , B c n ),
or equivalently:

A 0 ⊆ A 1 ⊆ . . . ⊆ A n B c 0 ⊆ B c 1 ⊆ . . . ⊆ B c n ,
It is easy to show that by choosing an initial compact subset B such that T i (B) ⊆ B ∀i ∈ Σ (see appendix A) the following inclusions hold:

T n (B) ⊆ . . . ⊆ T 2 (B) ⊆ T(B) ⊆ B,
or equivalently:

B c = B c 0 ⊆ B c 1 ⊆ B c 2 ⊆ . . . ⊆ B c n = (T n (B c )).
The interior part A = A 0 of the initial partial solid must thus be chosen to satisfy the following constraint:

A ⊆ T(A).
If the interior A is chosen properly, the sequence of partial solids (A n , B c n ) constitutes a directed set and the sequence of Op(B n ) converges to the required image Op(A). Note that if the initial interior part A 0 is chosen as the empty set, the sequence of partial solids (A n , B c n ) also constitutes a directed set. However, in this case the lub of this set would be (∅, A c ).

Let B ∈ H(X) be a non-empty compact subset and A ⊂ B be a non-empty open set. Let the sequence SB n = (A n , B c n ) be a directed set constructed by applying the Hutchinson operator to (A, B c ). In order to identify the lub of this sequence, consider the partial solid (Q i , Q e ) such that:

(A n , B c n ) (Q i , Q e ) ∀n.
By the definition of partial order in the solid domain, this system of inequalities is equivalent to the following:

A n ⊆ Q i and B c n ⊆ Q e , or A n ⊆ Q i and Q c e ⊆ B n . (4.2)
Now we consider the relation between the attractor A and the set Q e . Let us assume that A ⊂ Q c e . Let a be a point of Q c e \ A. This means that a does not belong to the attractor and there exists a non-zero distance between them:

d(a, A) = > 0.
Since the sequence of compact subsets B n converges to the attractor A with respect to the Hausdorff distance, there exists a number m ∈ N such that:

∀n m d X (A, B n ) <
This implies that the point a does not belong to any of these sets B n for ∀n m:

a / ∈ B n but a ∈ Q c e \ A ⊂ Q c e ⊆ B n .
The contradiction means that the set Q c e \ A = ∅ or equivalently:

Q c e ⊂ A ⇔ A c ⊂ Q e .
Both sets Q e and A c are open and the set Q e can thus be represented as the union of A c with an open set C ⊆ X, i.e.:

Q e = A c ∪ C. (4.3)
Since the open sets Q i and Q e are disjoint and by inclusion (4.2), we can therefore write the following:

A n ⊆ Q i ⊂ Q c e = A ∩ C c ⊆ A.
Thus:

Q i ⊆ A o . (4.4)
Hence the set A n is open and non-empty, there exists a compact subset A ⊂ A 0 and we can write the following:

T n (A ) ⊂ A n ⊂ A ∩ C c ⊆ A.
One can thus derive the following inequality for the distances:

d X (T n (A ), A ∩ C c ) d X (T n (A ), A) ----→ n→∞ 0.
By uniqueness of the limit, the sets A ∩ C c and A are equal, i.e.:

A ⊆ C c ⇔ C ⊆ A c .
We can thus simplify the definition (4.3) as follows:

Q e = A c ∪ C = A c . (4.5)
One can similarly show that:

Qi = A.
Thus we identified the following restrictions on the lub of the directed set SB n :

Q i ⊆ A o ; (4.6) Qi = A;
(4.7)

Q e = A c . (4.8)
That is, the lub (Q i , Q e ) exists only when the exterior part Q e is equal to A c and the closure of the interior part equals A.

Note that if the attractor interior is empty, that is A o = ∅, one can write: On the n-th iteration we obtain the set B n , which is close enough to A c , and the set A n , for which there is always some empty space near the middle of the square.

A n ⊆ Q i ⊆ A o = ∅,
Objects described by BCIFS are usually so-called just-touching IFS (CIFS). This means that the subdivided attractor parts intersect only by the boundary. The interiors of these parts do not intersect with each other:

∀i, j ∈ Σ, i = j T i (A o ) ∩ T j (A o ) = ∅.
By choosing the initial interior set A ⊆ A o we always obtain a case similar to that considered in the example. For these BCIFS, the empty set is the single option for the initial interior part in order to construct a directed set:

SB n = (∅, B c n ).
However, this directed set has the lub (∅, A c ).

Nevertheless, for some operations such as the distance from a point for example (see section 6.2), the interior part of the operator image SOp(SB) can be constructed even without approximating the interior part of the attractor.

In addition, constructing a directed set is not the only way to have the convergence SOp(SB) → SOp(SA). In some cases, even if the partial solids SB n = (A n , B c n ) do not constitute a directed set the convergence can still be stated.

Statement. If the sequence SB n converges to the attractor SA with respect to the Scott topology, then by the Scott-continuity of the operator SOP we can write the following:

SOp((A n , B c n )) → SOp((A o , A c )).
Consider the following example. Let {R; T 0 ; T 1 } be an IFS with the following transformations:

T 0 : R → R : x → x/2, T 1 : R → R : x → (x + 1)/2.
The attractor of this IFS system is thus the segment [0, 1].

Consider the images of an arbitrary interval (a, b) under the IFS transformations:

T 0 ((a, b)) = (a/2, b/2), T 1 ((a, b)) = ((a + 1)/2, (b + 1)/2).
We note that these images do not intersect when ba > 1.

Let (A 0 , B c 0 ) be the initial partial solid defined as follows:

(A 0 , B c 0 ) = ((1, 3), (-∞, 0) ∪ (3, ∞)).
Several first elements of the sequence SB n are illustrated in figure 4.16. One can prove that the sequence SB n = (1/2 n , 1 + 1/2 n ) converges to the attractor (A o , A c ) according to the Scott topology. Indeed, for any neighbourhood V of the attractor SA in the Scott topology there exists an open set U = (U i , U e ) SA. By the definition of open sets, the interior part U i of U has no common borders with the attractor A = [0, 1]:

Ūi ⊂ A.
There is a non-zero distance between these sets, which implies that ∃m ∈ N from which U i is included in all A n .

By analogy, there exists a number after which the exterior part U e of U is included in all B c n . This consequence implies that after a certain number of iterations we will exceed U e , that is U SB n , or by the definition of Scott open sets:

∃m ∈ N ∀n m SB n ∈ V.
The sequence SB n thus converges to the attractor A according to the Scott topology. The elements SB n are not comparable to each other and, therefore, do not constitute a directed set. Finally, constructing a directed set is not limited to constructing an increasing chain. To satisfy the definition of directed sets, we need to find a common maximum for each pair of elements:

(A n , B c n ) (C (n,m) , D c (n,m) ), (A m , B c m ) (C (n,m) , D c (n,m) ).
These inequalities are equivalent to the following:

C (n,m) ⊇ A n ∪ A m , D (n,m) ⊆ B n ∩ B m .
We consider a special case of this property, called IFS with condensation. Let A = A 0 ⊆ X be an initial open set and B = B 0 be the attractor bounding ball such that T(B) ⊆ B. Let A ⊂ B for (A, B c ) be a partial solid. All the following partial solids are constructed by applying the Hutchinson operator.

Let

C 0 = A 0 and C 1 = A 0 ∪ A 1 = A 0 ∪ T(A 0 ) = T (C 0 )
, where T is the Hutchinson operator associated with the IFS: {X; F, T i | i ∈ Σ}. Thus we add a transformation F : K → A 0 , for K ⊆ X to the set of transformations.

Introduce the two following sequences:

C n = T (C n-1 ) = T n (C 0 ) = A 0 ∪ A 1 ∪ . . . ∪ A n , D n = B n-1 ∩ B n = B n .
The sequence (C n , D c n ) thus constitutes a directed set. In the beginning of this section, we identified the restrictions on the lub of the directed set SB n . The new directed set (C n , D c n ) thus has the same restrictions and does not converge to the attractor SA according to the Scott topology (see figure 4.17).

Finally, one can imagine a completely different method to construct a directed set converging to the attractor. For example, after applying the Hutchinson operator we can merge certain subdivided parts based on the adjacency constraints of the BCIFS. We leave this question open for future studies. In conclusion, we have described a generic algorithm calculating an approximate operation image for a fractal shape represented by an IFS or a CIFS. We use the solid domain to ensure the calculability of most standard CAD operations, and we can therefore state that whenever the initial approximation SB is properly chosen, the sequence of operator images SOp(SB n ) converges to the required result SOp(SA). However, approximating the attractor itself can be time consuming. Applying the operator directly to SB n is thus not always the optimal way to compute the result. In the following sections we consider certain properties of the operator that allow us to optimize the SOp(SB n ) construction.

Image decomposition

In the previous section we established the generic algorithm and proved convergence to the required result. In the following sections we show how this algorithm can be optimized by exploiting certain specific properties of fractals.

Approximating the attractor can be time-consuming, because the number of objects (in T n (SB)) grows exponentially as the number of iterations increases. First of all, each IFS transformation is applied to the initial partial solid SB. On the next iteration, we repeat the same procedure. Since T (A ∪ B) = T (A) ∪ T (B) ∀A, B ∈ SX, the calculation of SB n can be illustrated by a tree whose nodes are images of the initial compact subset SB obtained by applying sequences of the IFS transformations (see figure 4.2) [START_REF] Edalat | Dynamical systems, measures, and fractals via domain theory[END_REF]. The union of compact subsets in the n-th level of the tree represents the corresponding approximation SB n .

If SOp satisfies certain conditions, the construction of SOp(SB n ) can be optimized by exploiting the self-similarity property of our iterative system. In the following explanations we assume that the required number of iterations n is chosen according to the effective structure of the domain used, as described in section 4.3.2 (in simple cases, this method reduces to the expression 4.1).

In the algorithm, we calculate the approximation SB n and then apply the operator SOp to this approximation. The complexity of the algorithm is thus O (f SOp (N n )), where f SOp (x) is the complexity of the operator SOp evaluated on the union of x objects and N is the number of IFS transformations. Generally, applying SOp to the union of an exponential number of objects is expensive. The idea of this optimization is to evaluate SOp(SB n ) by applying the operator to more simple objects, that is, to each leaf of the tree (see figure 4.18). Statement. The property of the image decomposition is not always verified. For example, calculating the volume V ol cannot be decomposed by Θ:

∀Θ : SY × SY → SY ∃A, B ∈ SX : V ol(A ∪ B) = V ol(A) Θ V ol(B).
Indeed, the volume operator is translation invariant and on the right side we therefore do not have the information about the relative positions of A and B or their intersection. That is, we cannot determine the common part of these sets to identify V ol(A ∪ B). Nevertheless, in some cases, the volume of the IFS attractor can be computed as a sum of the volumes of its parts:

∀i, j ∈ Σ V ol(T i (A) ∩ T j (A)) = 0 ⇔ V ol(A) = i∈Σ V ol(T i (A)).
Defining the operator Θ allows us to redefine the algorithm in many ways, depending on how we traverse the evaluation tree. This choice primarily affects the order in which we combine the results of the tree leaves. There exist two canonical ways to traverse the tree defined up to the ordering of branches: the iterative breadth-first and the recursive depth-first algorithms. Our generic algorithm can accordingly be redefined in these two ways. We describe the two algorithms in the following sections.

Breadth-first algorithm

In this section we describe the breadth-first algorithm for approximating the operator SOp image. We choose the required number of iterations n according to a given accuracy ε, as described in section 4.3.2. The idea of the breadth-first algorithm is to first calculate all the operator results at the n-th level of the tree, and then compose these results by means of the operation Θ.

The following function illustrates an implementation of this algorithm:

SOp(T (SB))

End

Note that in the breadth-first algorithm we apply the operator SOp to a transformed initial partial solid in each node of the n-th tree level, i.e. N n times. The operation Θ is evaluated only once to compose all the operator results. The complexity of this algorithm is therefore:

O (N n • f SOp (1) + f Θ (N n )) = O (N n + f Θ (N n )) ,
where f SOp is the complexity of the operator SOp evaluation, f Θ is the complexity of the operation Θ as a function of the number of parameters and N is the number of IFS transformations. Note also that f SOp (1) does not depend on the number of iterations.

Depth-first algorithm

In this section we describe the depth-first algorithm to approximate the image of the operator SOp. Here also, we choose the required number of iteration n according to a given accuracy ε, as described in section 4.3.2. A depth-first algorithm recursively calculates results for each branch of the tree and then composes the results by using the operation Θ.

The following recursive function illustrates an implementation of this algorithm:

Function DFAlgo Global variables: SB is the initial partial solid SOp is the operator required to evaluate {T i } i∈Σ is the set of IFS transformations Input: T is a composite IFS transformation n is the number of levels left to pass Output: the approximate operator result SOp(T n (SB)) Body:

if n = 0 // End of the recursion return SOp(T (SB)) endif return Θ i∈Σ DFAlgo(T • T i , n -1)

End

In the depth-first algorithm we also apply the operator SOp to a transformed initial partial solid in each node of the n-th tree level, i.e. N n times. However, in this case, the operation Θ is evaluated in each recursive call, i.e. 1+N +N 2 + . . . + N n-1 = N n -1 N -1 times to compose N results. The algorithm complexity is therefore:

O N n • f SOp (1) + N n-1 • f Θ (N ) = O (N n ) ,
where f SOp is the complexity of the operator SOp evaluation, f Θ is the complexity of the operation Θ as a function of the number of parameters and N is the number of IFS transformations.

Statement. The described algorithms have different complexities and the optimal choice depends on the operation Θ. If it is better to compose the N n results once, we choose the breadth-first algorithm, and if it is better to combine N operator results N n -1 N -1 times, we use the depth-first algorithm. One can prove that if evaluating the operator Θ is more complex than linear, i.e.:

lim x→∞ f Θ (x) x = ∞, (4.10)
then the depth-first algorithm is obviously the better choice.

Adaptive tree cutting

In this section we describe how to optimise the algorithms, described in section 4.4, by collecting information about the solution in intermediate levels of the evaluation tree. Both algorithms have an exponential complexity according to the number of iterations n, since we explore the entire tree to approximate the operator result. Due to the self-similarity of IFS attractors, in some cases we can predict the operator results in the leaves of the tree by calculating SOp in intermediate levels of the tree. If these results do not affect the overall composition by Θ, then it is not necessary to compute the corresponding branch of the tree. We call this method adaptive tree cutting optimization.

To describe this optimization more precisely, we introduce the following notation. Firstly, we define the subdomain S b Y of bounded partial solids:

S b Y = {(A, B) ∈ SY | B c is compact} ∪ {(∅, ∅)}, ordered by inclusion. Re- call that SY is the target space of the operator SOp. The continuous predi- cate [EL02] ⊆: S b Y × SY → {+, -} ⊥ is then defined by: (A, B) ⊆ (C, D) =      + if B ∪ C = Y, -if A ∩ D = ∅, ⊥ otherwise.
Definition. The information about where the solution can be found, contained in each node of the evaluation tree, is called solution set. We denote the solution sets by SOp(T (SB)), where T is a composite IFS transformation. The solution set thus corresponds to an approximate part of SOp(SA). All these sets are subsequently composed to approximate the resulting image SOp(SA). Since this image is composed iteratively, to optimize we can define the approximate result as a global variable and use it in the following iterations.

Statement. Consider a solution set of some node SOp(T (SB)), where T is a composite IFS transformation. By composing this set with the global solution set S, we can deduce whether it is suitable, i.e. whether it changes the overall composition. More precisely, if the following conditions hold:

SOp(T (SB)) SOp(T • T(SB)), (4.11)
SOp(T (SB)) ⊆ S, (4.12)

we can stop iterations for this branch of the tree without changing the resulting solution set S. Relation here is the partial order of the solid domain. Condition (4.11) implies that by applying the Hutchinson operator we increase the precision of the solution set and that if this current solution does not improve the result (see condition (4.12)), we can then avoid unnecessary computations. Note that condition (4.11) is always verified when the partial solids SB n constitute a directed set, due to the continuity of SOp and T . Otherwise, this condition has to be verified additionally. Hereafter, we use the term stopping criterion to refer to condition (4.12) and assume that the initial partial solid SB is properly chosen (see section 4.3).

Composing the solution sets by Θ at each iteration, we thus restrict the set of possible solutions, or, in terms of domain theory, we approximate the resulting partial solid SOp(SA).

The algorithms described in section 4.4 can both be optimized by constructing a global solution set S, i.e. an approximation of the solution, and by verifying the stopping criterion (4.12). The following recursive function illustrates this optimization applied to the breadth-first algorithm:

Function BFAlgoOpt Global variables: SB = (∅, B c )

SOp(T (SB))

End

The price we must pay for this optimization is applying the operator N n+1 -1 N -1 times, instead of N n , that is roughly N N -1 times more, and applying Θ roughly 2 times more, as illustrated in figure 4.19. More precisely, the complexity of the optimized breadth-first algorithm is:

O N n • f SOp (1) + f Θ (N n ) + N n-1 • f Θ (2) = O (N n + f Θ (N n )) ,
where f SOp is the complexity of the operator SOp evaluation, f Θ is the complexity of the operation Θ as a function of the number of parameters and N is the number of IFS transformations. Note that the best-case running time for this optimized algorithm is O(1). This is the case in which SOp(SB) returns the empty partial solid (∅, Y).

One can also, by analogy, adapt the depth-first algorithm. The following function thus shows the complexity of the optimized depth-first algorithm:

O N n • f SOp (1) + N n-1 • f Θ (N ) + N n • f Θ (2) = O (N n ) .
In optimistic non-trivial scenarios the depth-first algorithm shows better running-times. However, on average, the breadth-first algorithm cuts the tree much earlier than the depth-first algorithm and the breadth-first algorithm therefore has a better average running-time (see section 7). The evaluation tree of the optimized depth-first algorithm is presented in figure 4.20. Figure 4.20: Modifying the evaluation tree that allows us to optimize the depthfirst generic algorithm. We apply the operator SOp on intermediate levels of the tree to prevent unnecessary calculations as early as possible. During the calculation we construct the solution set S and verify the stopping criterion 4.12.

Image with an unknown interior

Let T be a composite transformation considered in an arbitrary step of the algorithm. To test the inclusion SOp(T (SB)) ⊆ S we must be able to compute the interior part of these partial solids.

We have already described the constraints required for the interior part of the initial partial solid SB (see section 4.3.4). It is sometimes difficult to construct the interior of SOp(SB), and accordingly the interior of all SOp(T (SB)).

We therefore use the empty set ∅ to approximate the interior; that is, only the exterior is used to approximate the solution. In this case, the stopping criterion (4.12) is verified only for the empty partial solid (∅, Y). To implement, we can then simply verify if SOp(T (SB)) is equal to (∅, Y) and avoid calculating the solution set S, as illustrated in figure 4.21. This leads to slightly different algorithms with the same complexities. For instance, we consider an operator determining whether a given point belongs to the attractor. Let p ∈ X be a point; the extended operator SOp is then defined as follows:

SOp((A, B)) =      + if p ∈ A, -if p ∈ B, ⊥ otherwise.
In this case, the space Y is a set {+, -} ⊥ without a metric. The operator SOp satisfies all conditions and we can write the following:

p ∈ i T i (B) ⇔ OR i (p ∈ T i (B)) .
The composing operator Θ is therefore the boolean operator OR.

Note that, since we do not construct the interior of the initial partial solid SB, the operator SOp never returns +. We therefore cannot determine whether or not the point belongs to the attractor. The problem is therefore semidecidable. However, we can bound the distance from point p to the attractor, and after n iterations, we can determine whether it is small enough, i.e.:

SOp(SB n ) = ⊥ ⇒ p / ∈ B c n ⇒ p ∈ B n ⇒ d X (p, A) < ε.
In this example, the stopping criterion (4.12) is only verified for the element -. That is, we need to determine if SOp(SB) = -. If this is the case, we stop the algorithm and deduce that SOp(T n (SB)) = -, that is, the operator result is false. Otherwise, we determine SOp(T i (SB)), ∀i ∈ Σ. If all results are equal towe deduce that SOp(T(SB)) =and that the point does not belong to T n (B). If there exists i 0 such that SOp(T i0 (SB) = -, we then continue explorations by determining SOp(T i0 T j (SB), ∀j ∈ Σ. For the images that give we stop the exploration, for the others we continue up to level n.

Conclusion

In this chapter we studied the possibility of applying an arbitrary CAD operation to fractals represented by a BCIFS (see section 3.3). We classify most standard constructing operations by adapting them to BCIFS. Thus we encountered three categories depending on the type of result. The first category corresponds to cases, in which we have a formal expression of the operation result. Applying operations from the second category is not immediate. The properties of fractals are more specific than those of conventional geometric objects (concepts of roughness or lacunarity, for example) and we have to study the underlying concepts of these operations in order to generalize or to adapt. The last category consists of operations for which we can construct an approximate result without specifying a formal representation.

In studying and identifying constraints on the operations to evaluate an approximate result, we introduced seven algorithms, which are the following:

• direct application to the approximate attractor;

• the breadth-first and the depth-first algorithms using image decomposition (see section 4.4);

• two optimized algorithms using the additional stopping criterion (see section 4.5);

• two simplified algorithms without constructing the interior part (see section 4.5.1).

These algorithms are generic, i.e. they can be applied to any CAD operation satisfying certain constraints once we implement the required interface. The first step is to verify the continuity of the operator. If it is possible to construct a continuous extension of the operator on the solid domain, we can thus construct the approximate image by using the effective structure on this domain. Moreover, this allows us to optimize the approximations. Otherwise, we cannot say with certainty that the sequence of approximate images converges to the required result.

If there exists an associative operation Θ described in detail in section 4.4, the approximate image can be obtained by composing images of each leaf, i.e. we can decompose the approximate result on separate parts. This leads to the breadth-first and depth-first algorithms.

Image decomposition allows us to subdivide the calculations in order to abandon the unnecessary parts, therefore noticeably optimizing the calculations. By using the solid domain we formalize solution set notation. Thus, we identify a set of possible solutions which are suitable for the current approximation level. This set is therefore restricted while iterating the algorithm. The presented algorithms can thus be optimized by verifying the stopping criterion (4.12). In certain cases, this criterion can also be simplified by simply comparing the result with (∅, Y), as described in section 4.5.1.

At the moment, we cannot provide the full automation of this process. This means that to apply a CAD operator we have to perform some preliminary work. Firstly, we need to define the extended operator, which must be continuous on the solid domain. We then have to design and code the image decomposition operator Θ, which may be any algorithm satisfying the definition (4.9). Finally, we must be able to verify the stopping criterion (4.12), by calculating the inclusion. Once we implement this interface we can apply the presented generic algorithms.

In the next chapter, we define a self-similarity property of the operation and introduce a new generic algorithm to formally represent an approximate operation result, i.e. to compute a specific CIFS with the approximate image as the attractor. In special cases, when the automaton degenerates into a single state, the constructed model is a specific IFS and has a linear complexity on the number of iterations.

In chapter 6, we consider some examples and special cases to evaluate different CAD operations. Based on the operation properties, we choose the appropriate algorithms presented in this section. We then analyse results and compare the performance of our algorithms.

Chapter 5

Formal representation of the operator image

In chapter 4 we presented several algorithms which could be used to approximate an operation. Once we implement the required interface, the operation image can be computed for the required accuracy ε. In this chapter, we identify some additional properties of the operation used to formally represent the approximate result. The idea is to build a specific iterated system (CIFS or IFS) for which the operation image would be the attractor. To this end we have to formalize the self-similarity property of the constructing set.

If the operation image satisfies the self-similarity property, we can then construct a CIFS with the generalised Hutchinson operator generating the attractor close to the operation image. First, we define this property and then introduce a new generic algorithm to construct this CIFS. In order to keep the process as simple as possible, we describe the idea on IFS. Since CIFS can contain transitions between spaces of different dimensions, generalizing to CIFS involves additional manipulations that we discuss in section 5.1.1.

Because of the semi-decidability of the real numbers comparison, the result can only be computed approximately. The final construction algorithm thus computes a CIFS automaton for a given accuracy γ > 0. In certain cases, when the automaton degenerates into a single state, the constructed model represents an IFS and has a linear complexity on the number of iterations.

Self-similarity property

In all the algorithms described in chapter 4, we apply the extended operator SOp to a transformed initial partial solid B in each node of the n-th tree level, i.e. N n times (see figure 4.18), or even more in the optimized versions. These calculations can be expensive if the operator SOp has a relatively high complexity or if the number of iterations n is large, sometimes these calculations can also be optimized.

Let us remind you that we approximate the operation by constructing a sequence of partial solids:

SOp(SB n ) = SOp((A n , B c n )) ----→ n→∞ SOp((A o , A c )).
Here, (B n ) n∈N = T n (B), (A n ) n∈N = T n (A) and T is the Hutchinson operator corresponding to a given iterative system. The initial partial solid SB = (A, B c ) is chosen as described in section 4.3.4.

In this section we define the self-similarity property of the images SOp(SB n ). The constructing algorithm uses the image decomposition, described in section 4.4, and we also assume that the operator SOp is continuous to guarantee the convergence of the sequence SOp(SB n ) to SOp(SA) (see section 4.3.4).

The principle of the self-similarity property is that instead of applying the operator to transformed partial solids, we factor out the transformation, and thus always apply the operator to the same partial solid SB before transforming it, as illustrated in figure 5.1. However, the operator has to be modified, in this case, preserving the continuity and image decomposition, described in chapter 4.

Definition. The self-similarity property of an operator SOp image is defined for all transition operators T of the model as follows:

∀T ∃ T : SY → SY such that SOp(T (SB)) = T SOp T (SB), (5.1)
where SOp T is the operator SOp modified using the corresponding transformation T .

Figure 5.1: Self-similarity property of the operator SOp. The principle is that instead of applying the operator to transformed partial solids, we factor out the transformation, and apply the operator to the initial partial solid SB before transforming it. However, the operator has to be modified preserving continuity and image decomposition.

If the operator SOp is, for example, a binary operation applied to a fractal A and an object K:

SOp : A → A K,
we can write the following:

SOp(T (SB)) = T (SB) K.
In simple cases, when transformation T has distributivity under binary operation , i.e.:

T (A B) = T (A) T (B),
and the transformation is right-invertible, by factoring it out we obtain

SOp(T (SB)) = T (SB) K = T (SB T -1 right (K)) = T SOp T (SB),
where T -1 right is the right-inverse operator. However, generally, transformation T does not have this distributivity property, and we therefore have to perform some additional actions to maintain the equality, that is:

SOp(T (SB)) = T (SB) K = T (SB T -1 right (K)) = T SOp T (SB).
Here T -1 right is the right-inverse operator, transformation T can differ from T . If transformation T is not invertible, by factoring it out we obtain

SOp(T (SB)) = T (SB) K = T (SB K ) = T SOp T (SB),
where

T (K ) = K.
We should thus find all points that are transformed into K under T .

Examples and special cases

As a first example, consider the operator intersecting an IFS {X;

T i | i ∈ Σ} with a line L: SOp : SA → SA ∩ L : SX → SX.
Assume that all the IFS transformations are invertible. One can prove that for any non-empty partial solid SB ∈ SX and an invertible transformation T : SX → SX the following property holds:

T (SB) ∩ L = T (SB ∩ T -1 (L)),
or in terms of the operator SOp:

SOp(T (SB)) = T (SOp T (SB)),
where SOp T is the same intersecting operator applied to another object T -1 (L), that is, SOp T (SB) = SB ∩ T -1 (L). This operator thus has the same properties (continuity, image decomposition) as the initial operator SOp and the image is self-similar [START_REF] Hart | Efficient anti-aliased rendering of 3d linear fractals[END_REF][START_REF] Gentil | Les fractales en synthése d'images: le modèle IFS[END_REF].

Consider another example introduced in section 3.2.3, which corresponds to a case of non-invertible transformation. In a CIFS each state of the automaton is associated with a space. A transformation between two spaces with different dimensions is generally non-invertible. In this example we demonstrate how to deal with such cases. Figure 5.2 illustrates the attractor of a given CIFS and a line to intersect with it. The CIFS is described by an automaton with three states: , C and S. State C corresponds to a square, and is associated with the 4D barycentric space X C , illustrated in figure 3.10. The result is then projected onto the 2D modelling space X = R 2 . State S is used to construct the Sierpinski triangle in the 3D barycentric space X S , which is then projected onto X C as one of the quarters. Let L be a line in the modelling space defined parametrically for t ∈ R as follows:

L = 0 1 + t 1 -1 = t 1 -t .
To simplify notation, we here use the equality sign, meaning that each point of set L has this structure for a given value of parameter t ∈ R, that is, in fact:

L = {p = t 1 -t | t ∈ R}.
Let SOp : SX → SX be extended to the solid domain operator calculating the intersection with a line L.

SOp : SA → SA ∩ SL,
where SL is a representation of the line L in the solid domain. We approximate the operation by constructing a directed set {SOp(SB n ) | n ∈ N}:

SOp(SB n ) ----→ n→∞ SOp(SA) = SOp((∅, A c )),
where SB n is an approximation of the attractor A = A in the solid domain, such that (see section 4.3.4):

SB 0 SB 1 . . . SB n .
Each transition δ(q, i) = w of the automaton is associated with an operator T q i : X w → X q . Looped subdividing operators are defined as follows:

T C 0 = H (1,0,0,0) 0.5 T S 0 = H (1,0,0) 0.5 T C 1 = H (0,1,0,0) 0.5 T S 1 = H (0,1,0) 0.5 T C 2 = H (0,0,1,0) 0.5 T S 2 = H (0,0,1) 0.5
, where H p s denotes the homothety centred at point p with ratio s. These operators are invertible, and we thus have the self-similarity property (5.1):

T q i = T q i , for i = 0, 1, 2 and q ∈ C, S, since the operators are continuous and therefore have distributivity under the intersection.

Let P be a projection to the 2D modelling space defined as follows:

P = 0 1 0 1 1 1 0 0 .
Projection P is not invertible and for the self-similarity property we therefore have to find a set of points, denoted L , such that:

P (L ) = L.
Let L be defined parametrically in a barycentric space for a, b, c ∈ R as follows:

L =     a b c 1 -a -b -c     .
We therefore have

P (L ) = 1 -a -c a + b = t 1 -t = L, which implies 1 -a -c = t a + b = 1 -t.
Here we have two equations for four variables. Set L thus defines a plane. Let t and c be parameters, we can write the following:

L =     1 -c -t c c t -c     .
This plane L is projected onto line L by P . This operator thus verifies the self-similarity property (5.1) and the operator SOp P is defined as follows: Since P projects the 4D subset L to the 2D subset L, the solution L of the equation P (L ) = L can always be found. The last transformation corresponds to another case, where we have to project a 3D subset to the 4D barycentric space. The projection T C 3 of the Sierpinski triangle to the space X C is defined as follows:

SOp P : SB → SB ∩ L ,
T C 3 =     0 0 0 0 0.5 0 0.5 0 0 0.5 0.5 1     .
Also we have to find a set of points, denoted L , that gives L after transformation by T C 3 :

T C 3 (L ) = L .
This equation generally has no solution, because we have to match the 3D subset L with the 4D subset L and the system would consist of four equations for three variables. However, the transformation T C 3 places the attractor A S following the control polygon ((0, 0, 0.5, 0.5), (0, 0.5, 0, 0.5) and (0, 0, 0, 1)). Thus we know that the set T C 3 (A S ) lies in a plane H formed by these control points. We have to find a set L that gives L ∩ H after transforming by T C 3 : T C 3 (L ) = L ∩ H. This is equivalent to finding a set L such that P T C 3 (L ) = L. We will use the second equation, as it is easier to implement.

One can derive the following:

P T C 3 = 0.5 1 1 0 0.5 0 .
Let L be defined parametrically as follows:

L =   e f 1 -e -f   ,
where e, f ∈ R are real variables.

We therefore have

P T C 3 (L ) = 1 -0.5e 0.5f = t 1 -t = L, which implies 1 -e 2 = t f 2 = 1 -t Therefore L =   e e 1 -2e   .

By construction, this line projects by T C

3 onto the intersection L ∩ H and then projects onto the line L. Instead of intersecting the transformed attractor T C 3 (A S ) with the plane L in the 4D barycentric space X C , we have now to intersect the attractor A S with the line L in the 3D barycentric space X S .

The self-similarity property (5.1) is satisfied for all the transition operators and the algorithm can therefore be applied to construct a CIFS automaton to approximate the operator SOp image. In the next section we introduce this construction algorithm. The intersection between the attractor and the considered line is illustrated in figure 5.4.

Approximate iterative model

In this section we introduce a generic algorithm to formally represent the approximate operation result. For simplicity purposes, we first describe the idea on IFS, a generalization to CIFS requires additional manipulations to verify the self-similarity property (see section 5.1.1).

Preliminaries

Let {X; T i | i ∈ Σ} be a hyperbolic IFS defined in a complete metric space (X, d X ) with the attractor A. An extension of an arbitrary operator Op : H(X) → Y to the solid domain SOp : SX → SY could be any algorithm that takes a nonempty partial solid as input (see section 4.3.4). We approximate the operation by constructing a directed set {SOp(SB n ) | n ∈ N}:

SOp(SB n ) ----→ n→∞ SOp(SA) = SOp((A o , A c )),
where SB n = T n (SB) is an approximation of the attractor A in the solid domain, and SB ∈ SX is an initial partial solid chosen as described in section 4.3.4.

Assume also that the operator SOp satisfies the image decomposition (see section 4.4) and is continuous to guarantee the convergence of the sequence SOp(SB n ) to SOp(SA). The self-similarity property of the operator SOp image is defined for each IFS transformation T i , i ∈ Σ as follows:

∃ Ti : SY → SY such that SOp(T i (SB)) = Ti SOp Ti (SB).

(5.2)

Here SOp Ti is the operator SOp modified somehow using the corresponding transformation T i (see section 5.1).

Construction algorithm

In this section we describe an algorithm to construct a specific CIFS with an approximate operator image as the attractor. The attractors of this CIFS are mutually recursively defined for each state q ∈ Q by the generalised Hutchinson operator:

SA q = Θ i∈Σ q
T q i (SA δ(q,i) ), (5.3)

We describe in detail the steps for constructing the CIFS automaton. In general, the algorithm produces an automaton with an infinite number of states, and we therefore have to bound the approximating accuracy to be able to compute the result in finite time.

The key idea of the algorithm is to apply a classic evaluation and verify the similarities between constructed sub-attractors, as illustrated in figure 5.5. Thus to the composition T α1 • . . . • T α k . If there exists a state q associated to the composition T β1 • . . . • T β l for which the following equality holds:

SOp Tα 1 •...•Tα k (SB) = SOp T β 1 •...•T β l (SB),
then the corresponding states are associated with the same sub-attractors. It is thus not necessary to add a new state, and we simply close the transition to the state q.

Figure 5.6 illustrates an example of a CIFS automaton produced by the presented algorithm. In general, the operator image SOp(SA) is an aperiodic set, and the algorithm thus produces an automaton with an infinite number of states. Sub-attractors equality is semi-decidable and is related to the issue of comparing real numbers, where we can only say if the numbers are close enough, but we cannot determine if the numbers are equal. By setting the accuracy at γ > 0 we instead verify if the sub-attractors are close enough. This allows us to determine whether the sub-attractors are relatively similar and, moreover, to bound the required number of iterations. Thus we verify the following inequality:

d SY SOp Tα 1 •...•Tα k (SB), SOp T β 1 •...•T β l (SB) < γ,
To simplify the notation, we use a relation γ to denote this inequality, i.e.:

SOp Tα 1 •...•Tα k (SB) γ SOp T β 1 •...•T β l (SB).
(5.4)

Note that in the evaluation algorithm, the constructing sub-attractors will be placed later on the output attractor by a composite of contractive transformations (see figure 5.5). The accuracy for this inequality, denoted , could be greater than the initial γ. More precisely, the accuracy should satisfy the following inequality:

γ • min{s(T α1 • . . . • T α k ), s(T β1 • . . . • T β l )}.
(5.5)

The inequality (5.4) holds when the sub-attractors are relatively close to each other. In this case, it is not necessary to overly complicate the automaton. By definition, the attractor A is a compact set, it therefore has a finite covering by sub-attractors which are not close to each other. Our algorithm thus always stops.

Consider the previous example of intersecting with a line L ∈ SX:

SOp : SA → SA ∩ L : SX → SX.
As described in section 5.1, for an invertible transformation T : SX → SX with a distributivity under the boolean intersection, i.e.

T (SB) ∩ L = T (SB ∩ T -1 (L)),
the self-similarity property (5.2) is satisfied, as we can rewrite this in terms of the operator SOp:

SOp(T (SB)) = T (SOp T (SB)),
where SOp T is the operator calculating the intersection with a line T -1 (L), that is:

SOp T (SB) = SB ∩ T -1 (L).
This operator has the same properties (continuity, image decomposition, selfsimilarity) as the initial operator SOp.

To determine similar sub-attractors we verify the following inequality:

T -1 α k • . . . • T -1 α1 (L) T -1 β l • . . . • T -1 β1 (L), (5.6)
where the tolerance can be chosen by the inequality (5.5). It implies the inequality (5.4) and does not contain applications of SOp, here we only verify whether two transformed lines are close enough. The constructing algorithm thus goes through the IFS evaluation tree and checks the emptiness of sub-attractors and the inequality (5.4). In the first case, we stop the iterations, in the second one we close the transition to the corresponding state.

S ← [] for ∀(t, s) ∈ S for ∀T i i ∈ Σ t ← t • T i if SOp t (SB) = (∅, Y)
go to next iteration endif // Find similar sub-attractors if ∃s : (tt, s ) ∈ Q and SOp t SOp tt , where = γ • min{s(t ), s(tt)} F.push(a new transition s --→

T i s ) else s ← a new state Q.push((t , s )) S .push((t , s )) F.push(a new transition s --→ T i s ) endif endfor endfor S ← S endwhile return (Q[1], F ) End
Here Q is a list of pairs (t, s), where s is a CIFS state and t is the corresponding composite transformation. In the return statement we denote by

Q[1] a set of states in Q, that is: Q[1] = {s | ∃(t, s) ∈ Q}.
The algorithm runs as follows. We start with an initial state . The corresponding composite transformation is obviously the identity map. We then loop through all the IFS transformations T i for i ∈ Σ, and check which of the alternatives the corresponding sub-attractor hits:

• the sub-attractor is empty. We do not add any special state to construct the empty set, so we continue the iteration.

• the sub-attractor is similar to one that was already parsed, denoted s . In this case, we close the transition at the state s .

• otherwise, we should continue the iteration recursively.

The algorithm returns a set of states and transitions composing a CIFS automaton. By construction, all the sub-attractors are close enough to corresponding parts of SOp(SA).

Statement. For the attractor SB of the output CIFS the following inequality holds:

SB γ SOp(SA) ⇔ d SY (SB, SOp(SA)) < γ.
As usual, we can approximate the CIFS attractor SB by an increasing chain of partial solids.

Let {X; T i | i ∈ Σ} be a hyperbolic IFS defined in a complete metric space (X, d X ) with the attractor A. Let SOp : SX → SY be an extended operation that satisfied all the required conditions to apply the presented algorithm, and γ > 0 be a calculation accuracy. The algorithm returns a CIFS automaton (Q, F ), where Q is a set of states and F is a set of CIFS transitions. Denote the attractor of this system by SB. By construction, one can state the following:

SB γ SOp(SA).
Statement. The number of states |Q| is equal to the number of iterations, for which we did not find any similar (for the accuracy γ) sub-attractor or for which the sub-attractor was empty. In other words, it is equal to a number of sub-attractors that make up the attractor, which are not similar to each other. Since the attractor is compact and the distances between sub-attractors are greater than γ, this number of states |Q| is finite.

The maximal number of states can roughly be calculated as follows. The idea is that we cover the attractor SB by sub-attractors that are not similar to each other. the similarity is defined with respect to the Hausdorff distance between the sub-attractors which cannot thus be less than γ. Assuming that all sub-attractors are points, we cover the attractor B by balls of radius γ. In the evaluation algorithm we can use a bounding ball as the initial compact subset. The radius r of this ball can be mutually recursively determined (see appendix A). By applying the evaluation algorithm we obtain a coverage of the attractor by balls of radius at most γ. The maximal number of iterations k can be computed as follows:

k = log(γ/r) log(s) ,
where s is the maximal contraction coefficient of the IFS transformations. In the worst case, we thus need N k balls to cover the entire attractor B. This means that the maximal number of iterations required to obtain a resulting CIFS for the accuracy γ, is equal to N k+1 .

|Q| |F | N k+1 ,
where N is the number of the IFS transformations.

In the worst case, the number of iterations exponentially depends on k, and k logarithmically depends on the accuracy γ, the algorithm complexity is thus:

O(N k ) = O γ r 1/ log N (s)
.

Exact iterative model for solution

In the previous section we introduced a generic algorithm to formally represent the operation result. In this section we discuss how to construct an exact iterative model of the CAD operation image, that is the CIFS with the attractor SOp(SA), defined by the generalised Hutchinson operator. Recall that the attractors of this CIFS are mutually recursively defined for each state q ∈ Q by the generalised Hutchinson operator:

SA q = Θ i∈Σ q T q i (SA δ(q,i) ), (5.7)
In fact, the presented algorithm expresses a finite list of well-defined instructions for formally representing the operation image. The only instruction that cannot be calculated exactly is the one used to verify sub-attractors selfsimilarity (5.4), because it involves real numbers comparisons, which are semidecidable. However, all these instructions can be carried out manually, that is if we can prove that two sub-attractors are self-similar, then the closed transition represents a part of the automaton with the exact solution as the attractor. The problem is that, in general, there is an infinite number of such verifications. Nevertheless, there exist cases where the solution SOp(SA) itself can be represented by an automaton with a finite number of states. In these cases, performing the algorithm instructions manually allows us to determine an exact model of the solution. Consider an example of intersecting the Sierpinski triangle A with a line L in the 3D barycentric space X, illustrated in figure 5.7. We define an IFS with the three contractive transformations T 0 , T 1 and T 2 , defined by homotheties centered at the unit vectors with the ratios 1 2 as follows:

T 0 = H (1,0,0) 1 2 T 1 = H (0,1,0) 1 2 T 2 = H (0,0,1) 1 2
Let the line L be defined as follows:

L =   t t 1 -2t   ,
where t ∈ R is a real variable. Let SL be a representation of the line L in the solid domain SX. Let SB ∈ SX be an initial partial solid, which represents a bounding ball centered in the point

c =   1/3 1/3 1/3   ,
with radius 2/3 and with an empty interior. In order to verify self-similarity on the first iteration (see figure 5.8) we must compare the line SL with T -1 i (SL) for i ∈ Σ. To avoid overly complicated explanations, we present here the calculations for the line L, instead of SL. For the first IFS transformation T 0 we thus have the following:

T -1 0 (L) =   1 -1 -1 0 2 0 0 0 2     t t 1 -2t   =   2t -1 2t 2 -4t   =   x x 1 -2x   ,
where x ∈ R. We therefore have the system of equations:

   2t -1 = x 2t = x 2 -4t = 1 -2x
, which obviously has no solution. This means that the lines are not the same, and neither are the sub-attractors. For the second IFS transformation T 1 we have the same, by symmetry.

Finally, for the last IFS transformation T 2 we can write the following:

T -1 2 (L) =   2 0 0 0 2 0 -1 -1 1     t t 1 -2t   =   2t 2t 1 -4t   =   x x 1 -2x   ,
where x ∈ R. We therefore have the following system of equations:

   2t = x 2t = x 1 -4t = 1 -2x
, which has the solution x = 2t. This means that the line T -1 2 (L) is the same as L, and the corresponding sub-attractors are also the same. We can thus close the current transition to the initial state.

On the next iteration, we compare the line

T -1 0 (L) with T -1 i T -1 0 (L) for i ∈ Σ. For i = 0 we have: T -2 0 (L) =   1 -3 -3 0 4 0 0 0 4     t t 1 -2t   =   4t -3 4t 4 -8t   .
This line does not intersect with the ball SB as well as the line T -1 2 T -1 0 (L). For i = 1 we can write the following:

T -1 1 T -1 0 (L) =   2 -2 -2 -1 3 -1 0 0 4     t t 1 -2t   =   4t -2 4t -1 4 -8t   .
Comparing this line with the line T -1 0 (L) we obtain the following system of equations:

   4t -2 = 2x -1 4t -1 = 2x 4 -8t = 2 -4x
, where x ∈ R. This system has a solution x = 2t -1 2 , which implies that the line T -1

1 T -1 0 (L) equals to the line T -1 0 (L), and the corresponding subattractors are therefore the same. We can thus close the current transition to the corresponding state. By analogy, we can deduce that T -1 Since all other lines do not intersect with the ball SB the algorithm stops and we thus obtain an exact CIFS with the attractor SOp(SA) = SA ∩ SL. The automaton is illustrated in figure 5.9. In special cases, when the automaton degenerates into a single state, the constructed model represents an IFS and has a linear complexity on the number of iterations. We discuss conditions for constructing an IFS in the following section 5.3.1.

1 (L) = T -1 0 T -1 1 (L).

Iterated function system

In section 5.2 we described the algorithm constructing a specific CIFS with the approximate operator image as the attractor. To determine similar subattractors we verify the inequality (5.4). In section 5.3 we showed that if the operation image is self-similar, it is possible to manually construct the exact CIFS of the image. In certain cases, when the inequality (5.4) holds for all the transformations, the output CIFS degenerates to an IFS. In this section, we identify the condition for constructing iterative model to be an IFS.

To avoid overly complicated explanations, we describe the idea in IFS, a generalization to CIFS is thus immediate with minor additional manipulations (see section 5.1.1).

Given a hyperbolic IFS {X; T i | i ∈ Σ} defined in a complete metric space (X, d X ) by a finite set of contractive transformations T i . Let SB = (A, B c ) ∈ SX be an initial partial solid such that T(B) ⊆ B and A ⊆ A o (see section 4.3.3). Suppose that the extended operator SOp is continuous and verifies the selfsimilarity property (5.2), i.e. ∀i ∈ Σ:

∃ Ti : SY → SY such that SOp(T i (SB)) = Ti SOp Ti (SB).
Statement. If all the modified operators are equal, i.e.:

SOp Ti (SB) = SOp Tj (SB) for ∀i, j ∈ Σ, (5.8) then the presented algorithm (see section 5.2) produce an IFS with the attractor SOp(SA) as shown in figure 5.11. The produced IFS has the same set of transformations as the initial one and it is possible to redefine a generalised Hutchinson operator T : SY → SY, by composing the transformations Ti :

T(SY ) = Θ i∈Σ Ti (SY ),
where SY ∈ SY. Θ here is the decomposition operator chosen as described in section 4.4.

Figure 5.11: Graph for operators satisfying the condition (5.8). The automaton is degenerated to a single state associated with the attractor SOp(SA). All the transitions are closed to this state. The result is an IFS.

The sequence defined by applying this operator thus converges to the operator image SOp(SA). Consider a sequence of partial solids Tn (SOp(SB)). For the first iteration (n = 1) we have:

T(SOp(SB)) = Θ i∈Σ Ti (SOp(SB)) = Θ i∈Σ SOp(T i (SB)) = SOp(T(SB))
and on the n-th iteration:

Tn (SOp(SB)) = SOp(T n (SB)) = SOp(SB n ) ----→ n→∞ SOp(SA).
To approximate the operator image SOp(SA) we thus apply the generalised Hutchinson operator n times, where the number of iterations n is chosen according to the effective structure of the used domain, as described in section 4.3.2 (in simple cases, by expression 4.1).

To verify the condition (5.8) for the constructing iterative model to be an IFS one must perform a theoretical analysis of the operation to verify the selfsimilarity property (5.2). Once the transformations Ti are determined and the condition (5.8) holds, we can define the generalised Hutchinson operator and thus construct an IFS with the attractor SOp(SA). Note that we thus obtain a precise formal model of the operation image.

In this case, the complexity of constructing SOp(SA) equals to the complexity of evaluating the IFS, that is it linearly depends on the number of iterations n:

O (n • f Θ (N )) ,
where f Θ is the complexity of the operation Θ as a function of parameters count and N is the number of IFS transformations. For instance, consider the convex hull operator Op = conv and an extended operator SOp defined by:

SOp : (A, B) → (conv(A), (conv(B c )) c ).
There exists an associative binary operation Θ merging two convex hulls, i.e.:

(A, B) Θ (C, D) = (merge(A, B), merge(C c , D c ) c ).
Suppose that the IFS transformations are affine, we can therefore write the following equality:

conv(T i (A)) = T i (conv Ti (A)) = T i (conv(SA)) ∀i ∈ Σ,
the self-similarity property (5.2) of the operator is therefore verified.

Moreover, the equality of all sub-attractors holds:

SOp Ti (SA) = SOp(SA) ∀i ∈ Σ.
All the conditions for constructing the IFS with the attractor SOp(SA) are satisfied. Let SB ∈ SX be an initial non-empty partial solid. It is thus possible to define the generalised Hutchinson operator:

T(SB) = merge i∈Σ T i (SB).
And the operator image can be determined by constructing the corresponding attractor:

SOp(SA) = T(SOp(SA)) = lim n→∞ ( Tn (SB)).
Since merging N convex hulls requires O(hN 2 ), where h is the number of convex hull vertices and N = |Σ| is the number of IFS transformations, the total time complexity for constructing the convex hull is O(nhN 2 ). Note that it linearly depends on the number of iterations n.

Constructing an iterative function system significantly reduces the number of operations performed to compute the operator image, and this therefore provides a qualitative gain in execution time. This is the fundamental advantage of this method, since there exist many tasks in which it is necessary to calculate the operator result in real time. Another important advantage is that this method does not depend on computational accuracy. Indeed, the operation image is represented formally.

Conclusion

In this chapter we introduced an algorithm constructing a CIFS with the generalised Hutchinson operator (5.7), whose attractor is close enough to the operation result with respect to the Hausdorff metric. The output automaton consists of states for sub-attractors of the operation image and transitions to embed these sub-attractors together. Each transition is associated with one of the initial transformations, the algorithm does thus not produce any new transformation.

On each iteration we check, which of the alternatives the corresponding subattractor hits:

• the sub-attractor is empty.

• the sub-attractor is similar to one that was already parsed.

• otherwise, we should continue the iteration recursively.

In the case of the empty sub-attractor, we do not add any special state to construct the empty set, so we continue the iteration. To preserve the selfsimilarity of the output set we verify the similarity of sub-attractors that are already constructed. Thus we defined a condition (5.1) to determine the similar sub-attractors. The self-similarity property is defined for every transition operator in the input model.

Because of the semi-decidability of real number comparisons, we cannot verify this equation automatically. In order to construct an exact model of the required solution, if it exists, we must prove all the necessary similarities manually.

To evaluate the algorithm in finite time one can set the computation accuracy γ. Thus we identified the self-similarity property in terms of approximate calculations and proved that the algorithm requires a finite number of iterations in this case. For a given accuracy γ the algorithm constructs a CIFS with the attractor, defined by the generalised Hutchinson operator, that is close enough to the operation image. Using the evaluation algorithm, described in section 3, we can approximate the attractor by the union of the transformed primitives. However, in the case of the set calculation tolerance γ, the distance between the attractor and the operation image cannot be less that γ.

In some cases, as demonstrated in section 5.3, the output iterative model degenerates to an IFS. We specify a condition when this is the case, however to verify this condition one must perform a theoretical analysis of the operation. Once the condition (5.8) holds, we can thus construct an IFS with the operation image as attractor. Note that we thus obtain an exact formal model of the operation image. Evaluating the operation is thus reduced to evaluating this IFS which has a linear complexity on the number of iterations.

Chapter 6

Examples and special cases

In this chapter we present various examples and special cases. We consider several examples of CAD operators applied to fractal shapes. We use the algorithms presented in this thesis based on operator properties. After discussing the examples we provide experimental results by comparing the running-time of the approximation algorithms and we then summarize all the results.

For purposes of simplicity, we describe most examples for IFS, a generalization to the CIFS is thus immediate (see section 5.1.1).

Prerequisites

In this section we summarize the results already presented in chapters 4 and 5. More precisely, we describe methods which need to be implemented manually in order to apply the presented generic algorithms.

Given a CIFS defined by an automaton (Σ, Q, δ), where Σ is an alphabet, Q is a set of states and δ is a transition function δ : Q × Σ → Q. An initial partial solid SB q ∈ SX q is associated to each state q ∈ Q. Thus we have to implement a function GetInitial : q → SB q that gives the initial partial solid for a given state q ∈ Q.

After that, the operator SOp : SX → SY has to be defined. For this we should decide how the solution will be represented, that is to say that we define the nature of solids in SY. In the presented algorithms, the following methods are applied to partial solids in SY:

• The inclusion predicate ⊆;

• GetAccuracy -returns an accuracy value for a given partial solid (C,

D) ∈ SY, i.e. d Y (C, D c ).
The first method is required to optimize calculations, and the second is used to stop the algorithm when the required accuracy is achieved. In cases where the solution interior is empty we use the simplified stopping criterion (see section 4.5.1) and do not use the GetAccuracy method. Note that we do not need to be able to represent any partial solid in SY. We only deal with the transformed initial bounding balls SB q and their compositions by Θ. For example, when Θ is the union operator, the partial solid can be represented as a list of pairs (T, SB q ), T is a composite CIFS transformation and SB q is the initial partial solid for a given state q ∈ Q. After defining these methods the operator SOp itself should be implemented.

If the operator image can be decomposed, as described in section 4.4, the corresponding operator Θ : SY × SY → SY has to be implemented as well. In the breadth-first algorithm this operator could also compose a whole list of intermediate results, i.e. Θ : SY k → SY, where k is the number of input results.

Figure 6.1 illustrates a class diagram for our algorithms. We also implemented GetDFA, GetBFA, GetDirectAlgo methods returning the corresponding algorithm to apply to a CIFS. The SetX and SetY classes represent partial solids respectively of SX and SY. The SetY class is thus abstract, because its implementation depends on the operation.

Constructing the formal representation, as it is described in chapter 5, requires additional theoretical analysis of the operator. First we should verify the self-similarity property (5.1) and define all the corresponding transformations T : SY → SY. Second, we have to be able to verify the relation γ (see the inequality (5.4)) between the automaton's states. Here we can store the approximate attractor for each state, however in some cases the inequality can be transformed to an even simpler comparison (see the inequality (5.6)).

During the study of the self-similarity property (5.1) one can also verify the condition (5.8). If this condition is verified, applying the operator SOp is thus reduced to applying the generalised Hutchinson operator. Otherwise the approximate solution can be generated. Once we identify the transformations T and implement the method to verify self-similarities, the algorithm, presented in section 5.2.2, can be applied to produce a CIFS with an attractor that is close to the operator image SOp(SA).

Finally, after computing the approximate partial solid SOp(SA) one should be able to obtain the required approximate solution Op(A) by applying the BuildSolution method.

These are all the methods which have to be implemented to apply the operator to fractals.

Distance from a point

In this section we consider the operator for computing the distance from a given point. We also present the implementation of prerequisites, described in section 6.1.

Given an IFS defined by a finite set of contractive transformations {T i } i∈Σ and the operator Op : X → R + defined as follows:

p ∈ X Op : A → d(p, A) = min a∈A d(p, a).
It is easy to show, that for any sequence (B n ) n∈N of non-empty compact subsets of X converging to A, the corresponding sequence of the distances (d(p, B n )) n∈N converges to d(p, A), that is, the distance operator Op is continuous at A. We can thus apply the generic algorithm approximating the distance to the attractor.

Since the following equality holds: there exists an associative operation Θ = min{. . .} decomposing the image Op(B n ). The complexity of the operator Θ as a function of a count of parameters is therefore:

d(p, A ∪ B) = min{d(p, A), d(p, B)},
f Θ (x) = O(x).
The initial compact set B can be chosen as a point, that is to say, each time we apply the operator to a single transformed point. The complexity of the operator Op is constant:

f Op (1) = O(1).

As described in section 4.4, we defined two approximation algorithms. There is no difference between applying Θ one time to many elements and applying Θ many times to N elements (see 4.4.2). Both the breadth-first and depth-first algorithms have the same complexity O (N n ), where N is the number of the IFS transformations and n is the number of iterations chosen by equation (4.1). However, both algorithms have the same complexity as that of simply computing Op(B n ).

The optimization described in section 4.5 requires the notion of interior and exterior bounds of the approximate solution. In order to optimize the computation we should thus extend the operator Op to the solid domain SR + .

As described in section 4.3, we choose the initial compact subset B as an attractor bounding ball, such that T i (B) ⊆ B ∀i ∈ Σ. Thus we consider the operator defined by: On the next iterations we will precise these bounds, that is to say that we will restrict the set of possible solutions up to the required precision ε. The complexity of the extended operator SOp evaluated for a single transformed ball is the same.

SOp : (∅, B c ) → ((max
In the optimized algorithm we always compute the operator on a single transformed initial partial solid. By choosing the initial compact subset B as an attractor bounding ball B(c, r) as described in section 4.3, the lower and upper bounds can be simply computed as follows:

min b∈B d(p, b) = d(p, c) -r, max b∈B d(p, b) = d(p, c) + r
The operator Θ, in this case, is the continuous union operation:

(A, B) Θ (C, D) = (A ∪ B, C ∩ D).
The complexity of the operator Θ as a function of the number of parameters is therefore f Θ (x) = O(x), due to the fact that each composition requires two comparisons of real numbers.

The algorithm steps are the following. At the beginning, we identify the initial solution set, as illustrated in figure 6.2: To be able to restrict the set of possible solutions, we should not compare solution sets between different tree levels. Thus we reset S to be able to obtain a more accurate result on the next iteration. We then start a loop for all the IFS transformations T i , ∀i ∈ Σ.

S = SOp(SB) = ((d + r, +∞), (0, d -r)),
Let us consider, for example, the first transformation T 0 . We compute the solution set: and we continue the iterations for this branch of the tree by adding T 0 to the list (ts in the implementation BFAlgoOpt) of transformations for the next iteration.

For the next transformation T 1 we execute the same steps. We compute the corresponding solution set:

SOp(T 1 (SB)) = ((d 1 + r 1 , +∞), (0, d 1 -r 1 )),
where d 1 = d(p, T 1 (c)), r 1 = r • s 1 and s 1 is the contraction coefficient of T 1 . We verify the stopping criterion (4.12) for this solution set:

SOp(T 1 (SB)) ⊆ S =      + if d 1 -r 1 d 0 + r 0 , -if d 1 + r 1 d 0 -r 0 , ⊥ otherwise.
If it returns +, we can deduce that for all solutions in this branch of the tree there exists a better one in the previous branch associated with T 0 . We can therefore avoid unnecessary computations and stop iterations for this branch T 1 . Otherwise, we also add T 1 to the list (ts in the implementation BFAlgoOpt) of transformations for the next iteration.

After n iterations of the algorithm we obtain a solution set of the following form:

((b, +∞), (0, a)),

where a, b ∈ R + are correspondingly the lower and upper solution bounds. The exact solution is therefore somewhere in the segment [a, b]. We also know that the distance between the exact solution and a (the same for b) is inferior to the given precision ε. That is to say, any point between a and b is an appropriate result.

In the pessimistic case where the stopping criterion is never verified, we have to go over the entire tree. The worst-case running time is therefore exponential O(N n ), where N is the number of the IFS transformations and n is the required number of iterations. However, the best-case running time of the closest ball search is equal to the depth of the tree, that is to say, O(log N n ) = O(n).

Implementation

In this section we present the implementation of prerequisites for the operator used in computing the distance from a given point.

As described in section 6.1, the first function to be implemented determines the initial partial solid for a given state q ∈ Q, i.e. q → SB q . We compute the bounding balls B q for each CIFS attractor A q using the algorithm, presented in appendix A. The initial partial solids SB q are thus the pairs (∅, B q c ).

As described in section 6. In simpler terms, for a transformation T we approximate the transformed ball T (B(c, r)) by the ball B(T (c), s(T ) • r). Thus we denote the center of the transformed input ball SB by SB.center and its radius (or the semi-major axis for ellipses) by SB.radius.

The decomposing operator Θ is simply defined for two partial solids in SY by the following function: In the breadth-first algorithm this operator could also compose a whole list of intermediate results, i.e. Θ : SY k → SY, where k is the number of input results. For this we simply apply the operator Θ to the first two elements in the list, and we then recursively compose the result with the following elements.

Finally, after computing the approximate partial solid SOp(SA) of the type DistanceSetY, one should be able to obtain the real distance from the given point p to the attractor A. We therefore define a mapping SY → R as follows:

Function BuildSolution Input: A is an object of the type DistanceSetY approximating SOp(SA). Output: An approximate operator image Op(A) Body: return A.in+A.ex 2

End

These are all the methods we need to implement in order to use the presented algorithm.

Convex hull

In this section we consider the convex hull computing operator for a given affine IFS.

Consider an IFS defined by a finite set of affine contractive operators {T i } i∈Σ and the convex hull operator Op = conv. There exists an associative binary operation Θ merging two convex hulls, i.e.: conv(A ∪ B) = merge(conv(A), conv(B)).

One can prove that the operator conv is continuous at point A. Merging convex hull runs faster that recalculating the overall convex hull, because in merging we use information according to which the initial sets are convex [START_REF] Mishkinis | Approximate convex hull of affine iterated function system attractors[END_REF].

However, in order to apply the optimized algorithm, we should extend the operator to the solid domain. Let us consider the convex hull of the fixed points f ix(T i ) for ∀i ∈ Σ. We can only find fixed points approximately. Let δ > 0 be the accuracy used in calculating the fixed point. We then compute the offset on a distance -δ to ensure that the resulting set will be inside the exact convex hull conv({f ix(T i ) | i ∈ Σ}). This set is an interior part of the initial partial solid SB. The exterior can be chosen simply as a bounding ball of the attractor.

Let us consider an extended operator SOp defined by:

SOp : (A, B) → (conv(A), conv(B c ) c ).
One can see that SOp(SB) = SB and the interior part of the image can thus be constructed. The interior of partial solids SOp(T (SB)) for an arbitrary composite transformation T can be constructed as well. We therefore have both the interior and the exterior parts of the solution. Moreover, the operation images constitute an increasing chain, that is we can apply the optimization described in section 4.5. The operator Θ is merging:

(A, B) Θ (C, D) = (merge(A, B), merge(C c , D c ) c ).
Merging two convex hulls is linear on the number of points in the convex hulls. So, merging x convex hulls requires O(h • x 2 ) time, where h is the number of convex hull vertices. The complexity of the operator Θ as a function of the parameters count is therefore:

f Θ (x) = O(x 2 ).
Note that the equation (4.10) is satisfied, that is to say that, the depth-first algorithm is better than the breadth-first one.

The algorithm performs the following steps. We start by identifying the initial partial solid SB as described in section 4.3.4 and get down to the n-th level of the IFS tree. We then compute the solution set SOp(T n 0 (SB)) and therefore assign the global solution set to S = SOp(T n 0 (SB)). similarly, on the next N iterations we compute the solution sets SOp(T n-1 0 T i (SB)) for ∀i ∈ Σ and merge the corresponding convex sets by Θ.

On the next sub-branch of the tree we have to compute the solution set SOp(T n-2 0 T 1 (SB)) to verify the stopping criterion (4.12). If this solution set is inside the global solution set S, we can deduce that all the tree leaves will also lie inside S, that is to say, all this sub-tree has already been counted in the solution. We therefore stop iterations for this sub-tree. Otherwise we should continue to explore the tree, because each leaf can potentially extend the final result. Given an ε iapproximation of the convex hull, we calculate its images by applying all the IFS transformations to its vertices. We then merge all these images to obtain the ε i+1 -approximation of the convex hull.

T 2 (E i ) T 0 (E i ) T 1 (E i ) E i+1 E i
The process is repeated iteratively to achieve an approximate convex hull of the required precision. Figure 6.3 illustrates these iterations.

At a given time, we compute the operator Θ for N convex hulls maximum. This takes O(hN 2 ) time, where h is the number of convex hull vertices, and N is the number of the IFS transformations. The complexity of evaluating Θ is therefore constant on the number of iterations n. On the other side, the breadthfirst algorithm compute one final operation Θ on at most N n elements. This takes O(hN 2n ) time, which is exponentially worse than the chosen algorithm.

Note that if the output convex hull has an infinite number of vertices, these algorithms quickly go beyond the space and time limits, because merging depends on the number of convex hull vertices and is evaluated at most N n times to achieve the n-th level of the tree.

As described in section 5.3.1, it is possible to formally represent the convex hull. Since the IFS transformations are affine, the self-similarity property of the operator is verified:

∀i ∈ Σ conv(T i (SB)) = T i (conv Ti (SB)) = T i (conv(SB)).
Moreover, we also have the equality of all the sub-attractors (see 5.8):

conv Ti (SB) = conv(SB) ∀i ∈ Σ.
In other words, the algorithm starts with the initial state . Since all the sub-attractors conv Ti (SB) are equal, for each transformation T i we add a cyclic transition from the state to itself and then the algorithm stops. The produced generalised IFS constructs the attractor conv(SA) by iteratively applying the generalised Hutchinson operator T(SB) = merge i∈Σ (T i (conv(SB))), as described in section 5.3.1. Evaluating the convex hull is thus reduced to evaluating this IFS.

Note that in this method we do not construct the IFS evaluation tree with its exponential number of objects. Instead, we compose all the results on each iteration, which is, in fact, a more natural way of evaluating the IFS. In this case, both the space and time complexities are reduced, because our composing operator simplifies the resulting approximate object in contrast to the union operator.

Since merging N convex hulls requires O(hN 2 ), where h is the number of convex hull vertices and N = |Σ| is the number of IFS transformations, the total time complexity necessary for constructing the convex hull is O(nhN 2 ). Note that it linearly depends on the number of iterations n.

Simplification of the convex hull approximation

The presented results confirm our previous analysis for the convex hull of affine IFS [START_REF] Mishkinis | Approximate convex hull of affine iterated function system attractors[END_REF]. However, sometimes the algorithms could be optimized by using operator specific properties. For instance, the approximate convex hull can be simplified without any loss of accuracy. We simplify the approximate convex hull at each iteration by replacing vertices that are close enough to one another with one single vertex. To identify such vertices, we simply need to check the neighbouring points on the approximate convex hull.

r new ≤ ε i c new c ′ c r c ′
If two neighbouring points c and c satisfy:

d(c, c ) + r c + r c 2ε i (6.1)
we can then replace these points with the middle point without any loss of accuracy, as shown in figure 6.4. Thus, the number of points considered in constructing intermediate convex hulls is reduced and intermediate convex hulls are simplified without losing any accuracy. To avoid a constant rebuilding of data structures the simplification can be performed during the merge. This method of simplification can therefore easily be integrated into our approach by redefining the composing operator Θ in the following way:

(A, B) Θ (C, D) = (mergeAndSimplif y(A, B), mergeAndSimplif y(C c , D c ) c ),
where the method mergeAndSimplif y takes into account distances between neighbouring points, as described above.

Results and discussion

We examined the presented algorithm for various IFS attractors and different precision levels in 2D and 3D spaces. It should be pointed out that all tests were performed on an ordinary PC (Intel R Core TM 2 Duo T7500 2.2GHz 3GB RAM) without any additional computational facilities.

Several results for various IFS attractors are presented in figures 6.5 and 6.6. For each attractor, the program generates an approximation of the convex hull at the specified accuracy. The execution time of our algorithm is presented in Table 1 Here N is the number of the IFS transformations, s max is their maximal contraction coefficient and dim is the dimension of the space in which the attractor lies.

An approximation of the convex hull is computed in an interactive rate for all of the 2D attractors and also for the attractors whose convex hulls have a finite number of vertices. We compare our method for computing an approximate convex hull with Martyn's approach [START_REF] Martyn | The attractor-wrapping approach to approximating convex hulls of 2d affine ifs attractors[END_REF]. 2. Execution time of the algorithms for several IFS attractors with a given precision ε = 10 -3 . Here N is the number of the IFS transformations and s max is their maximal contraction coefficient.

One can see that Martyn's approach strongly depends on the maximal contraction coefficient s max whereas our algorithm depends more on the number of the approximate convex hull vertices h. The stopping criterion (4.12) can significantly reduce the number of operations, and our optimization therefore provides a qualitative gain in execution time.

Simplifying the convex hull is another important process. To improve collision tests performance it is necessary to verify whether an object intersects with the convex hull of our attractor. This process requires testing all the faces of the convex hull faces. The simplified convex hull involves a smaller number of vertices and faces, verifying the intersections will therefore take less time. The number of approximate convex hull vertices for ε = 10 On average, the efficiency of the simplification is about 46%, in 2D, but it decreases, when we increse the level of precision. In 3D, this process is much more efficient -we eliminate about 90% of vertices. However, such efficiency can be achieved only for the attractors whose convex hulls have an infinite number of vertices. Another important advantage is that we can improve the precision of the obtained approximation by further iterating our algorithm. That is to say, to improve the precision of the process it is not necessary to recalculate all k iterations, we start with the approximation first obtained and continue to iterate the algorithm. This can be useful advantage when it is necessary to compute a sequence of the approximate convex hulls (to obtain different levels of detail for example).

Boolean intersection

In this section we present an example of the operator used in calculating the intersection between a given IFS and a given set of X. We also present the implementation of prerequisites, as it was described in section 6.1.

Given an IFS defined by a finite set of contractive transformations {T i } i∈Σ and the operator Op : X → X defined as follows:

L ⊆ X Op : A → A ∩ L,
where L is an object intersect a given IFS.

The boolean intersection is not continuous for the Hausdorff metric. To apply the generic algorithms and also to optimize the process as described in section 4.5, we should extend the operator Op to the solid domain SX. More precisely, we consider the operator defined by:

SOp : (∅, B c ) → (∅, (B ∩ L) c ).
Here, the initial compact set SB is chosen as (∅, B(c, r) c ), where B(c, r) is an attractor bounding ball. In this case, we do not construct the interior part, because it is difficult and not always possible. The interior is therefore approximated by the empty set ∅, that is to say, only the exterior is used to approximate the solution, as described in section 4.5.1. We also consider some modified algorithms without constructing the global solution set S.

It can easily be shown, that the operator Θ is the union operator:

Op(T i (SB) ∪ T j (SB)) = (∅, ((T i (B) ∪ T j (B)) ∩ L) c ) = = (∅, ((T i (B) ∩ L) ∪ (T j (B) ∩ L)) c ) = = (∅, (T i (B) ∩ L) c ∩ (T j (B) ∩ L) c ) = = Op(T i (SB)) ∪ Op(T j (SB)),
where T i and T j are the IFS transformations, L is the object intersecting them, and B is the initial compact subset. The complexity of the operator Θ as a function of the parameters count is therefore linear:

f Θ (x) = O(x).
While evaluating the algorithm, we apply the operator SOp only to transformed balls. Intersecting between classic shapes does not depend on the number of iterations. The complexity of the operator SOp is therefore:

f SOp (1) = O(1).
As described in section 4.4, we defined two approximation algorithms. As composing by Θ is linear on the number of parameters, there is no difference between applying Θ one time to many elements or applying Θ many times to N elements (see 4.4.2). Both the breadth-first and depth-first algorithms have the same complexity O (N n ), where N is the number of IFS transformations and n is the chosen number of iterations as described in section 4.3.2.

We do not construct the interior and the stopping criterion (4.12) therefore reduces to the following:

SOp(T (SB)) = (∅, Y).
This simplification reduces the running-time, because we do not construct the global solution set S. Nevertheless, the asymptotic complexity remains the same.

Moreover, the self-similarity property is also verified, which allows us to factor out the transformation, and to always apply the operator to the same partial solid SB before transforming, as described in section 5.1. One can thus find the corresponding transformations Ti (see section 5.1.1) in order to construct the approximate iterative model (see section 5.2).

The modified operator SOp T , for a corresponding composite transformation T , differs from the original operator SOp by intersecting the initial partial solid SB with the modified object L. In the case of an invertible transformation T we have:

SOp T : (∅, B c ) → (∅, (B ∩ T -1 (L)) c ).
How to deal with the non-invertible transformation was described in section 5.1.1. The generic algorithm, described in section 5.2, produces a CIFS with the attractor that is close enough to the operator image SOp(SA).

In some cases, when the solution is exactly self-similar, the algorithm stops for any approximation accuracy. This is the case when only a finite number of semi-decidable equations have to be verified. If all these equations are proved manually the resulting CIFS formally represents the exact solution, that is to say, its attractor SOp(SA) is defined by applying the generalised Hutchinson operator. One can also compute a set difference and a symmetric difference by expressing these operators in terms of union, intersection and complement.

where B(0, r) is a ball centered at the origin with the radius r. The operator thus performs offsetting at a negative distance r.

Consider Op(A ∪ B) calculating, where A and B are two intersecting sets, as illustrated in figure 6.9. There is no relatively simple method to compose the results Op(A) and Op(B) in order to obtain Op(A ∪ B), because Op(A) and Op(B) do not contain any information about the intersection A ∩ B. The operation Op image cannot thus be decomposed and to offset on a negative distance we have to apply the operation directly to the approximate attractor, as described in section 4.2.

(A 1 , B 1 ) ⊕ (A 2 , B 2 ) = (A 1 ⊕ A 2 , (B c 1 ⊕ B c 2 ) c ).
This map is thus well defined and continuous [START_REF] Edalat | Foundation of a computable solid modelling[END_REF].

For this example, the initial compact set SB is chosen as (∅, B(c, r) c ), where B(c, r) is an attractor bounding ball, as for the distance example (see 6.2). We do not construct the interior part, because it is difficult and not always possible. The interior is therefore approximated by the empty set ∅, that is to say, only the exterior is used to approximate the solution, as described in section 4.5.1. The operator extended to the solid domain is thus:

SOp : SB → SB ⊕ SR,
where SR = (B(0, r) o , B(0, r) c ).

It can be shown, that the operator Θ is the union operator:

SOp(T i (SB) ∪ T j (SB)) = {x + b | x ∈ T i (SB) ∪ T j (SB), b ∈ B(0, r)} = = {x + b | x ∈ T i (SB), b ∈ B(0, r)}∪ ∪ {x + b | x ∈ T j (SB), b ∈ B(0, r)} = = SOp(T i (SB)) ∪ SOp(T j (SB)),
where T i and T j are the IFS transformations, SB is the initial partial solid. The complexity of the operator Θ as a function of the parameters count is therefore linear:

f Θ (x) = O(x).
Evaluating the algorithm only involves applying the operator SOp to transformed balls. If the IFS transformations are affine, the initial bounding ball is transformed into ellipses and offsetting these ellipses can be computed in constant time.

As described in section 4.4, we defined two approximation algorithms. Composing by Θ is linear on the number of parameters, there is thus no difference between applying Θ one time to many elements or applying Θ many times to N elements (see 4.4.2). Both the breadth-first and depth-first algorithms have the same complexity O (N n ), where N is the number of the IFS transformations and n is the number of iterations chosen as described in section 4.3.2.

If we do not construct the interior, the stopping criterion (4.12) is simplified to the following:

SOp(T (SB)) = (∅, Y),
This is never the case because:

(T (B) ⊕ B(0, r)) c = Y ⇔ T (B) ⊕ B(0, r) = ∅.
This simplification does not reduce the running-time, because we always reach all the leaves of the evaluation tree and there are N n leaves. Now we consider the self-similarity of the Minkowski addition. Assume that all the IFS transformations are affine and can therefore be represented as follows:

T

i (x) = L i (x) + v i , for i ∈ Σ, x ∈ X,
where L i is the linear part of the transformation and v i is a vector used for translation. We can thus write the following for any composite IFS transformation

T (x) = L(x) + v, x ∈ X: T (B) ⊕ B(0, r) = {T (x) + b | x ∈ B, b ∈ B(0, r)} = = {L(x) + v + b | x ∈ B, b ∈ B(0, r)} = = {L(x) + b -v + 2v | x ∈ B, b ∈ B(0, r)} = = {L(x + L -1 (b -v)) + 2v | x ∈ B, b ∈ B(0, r)} = = {T (x + T -1 (b)) + v | x ∈ B, b ∈ B(0, r)} = = { T (x + T -1 (b)) | x ∈ B, b ∈ B(0, r)} = T (B ⊕ T -1 (B(0, r))),
where T (x) = T (x) + v is the corresponding transformation providing the selfsimilarity. The operator SOp thus satisfies the self-similarity property and the corresponding modified operators are the following:

SOp T : SB → SB ⊕ T -1 (SR).
As described in section 5.1.1, in the case of non-invertible transformation T we have to find a set B such that:

T (B ) = SR,
where SR = (B(0, r) o , B(0, r) c ). Once we determine this set we can write the following:

T We could thus continue the CIFS automaton construction as described in section 5.2.

Figure 6.10: Worst-case running time of the optimized algorithms and the direct Op(B n ) calculating for the Sierpinski triangle. The operation is the distance from a point, defined as described in section 6.2.

Results and discussion

In this section we show the gain derived from using the optimized algorithms and we compare the presented algorithms on different examples. All the presented algorithms were designed and coded in C# language and tested on an ordinary laptop (Intel Core i7-2720QM CPU @2.2GHz).

Figure 6.10 illustrates the worst-case running-time for the presented algorithms. All the algorithms have an exponential running-time. However, the tree-cutting optimisation of the breadth-first algorithm is triggered earlier compared with the depth-first algorithm. As illustrated on the chart, the running time of the breadth-first algorithm is better than that of the depth-first algorithm. But the space complexity of the breadth-first algorithm is exponential, since on the i-th iteration we have to keep all the N i tree nodes in accessible memory to compute the approximate result. The algorithm therefore runs out of memory very quickly. This does not happen with the depth-first algorithm, since on each iteration we need to keep only one branch of the evaluation tree in accessible memory. Figure 6.11: Worst-case running time of the optimized algorithms without constructing the solution interior part, and the direct Op(B n ) calculating for the Sierpinski triangle. The operation computes the intersection with a line, as described in section 6.4. Figure 6.11 illustrates the worst-case running-time for the presented simplified algorithms. As described in section 4.5.1, we do not construct the interior part of the approximate result. Only the exterior is therefore used to approximate the solution. In this case, the stopping criterion (4.12) is only verified for the empty partial solid allowing us to simplify the verifications. The algorithms also have an exponential running-time. As in the previous example, compared to the depth-first algorithm, the breadth-first algorithm is triggered earlier and has a better running-time. However, the algorithm also runs out of memory very quickly, which does not happen with the depth-first algorithm. Compared to the direct algorithm both optimized algorithms still show a considerable gain in running time.

The best-case running time of both the optimized and simplified algorithms is constant, representing the gain over the direct algorithm which exponentially depends on the number of iterations.

The following table illustrates the complexities of the algorithms described in this thesis:

Algorithm Time Space Direct O (f SOp (N n )) O (N n ) Breadth-first (BFA) O (N n + f Θ (N n )) O (N n ) Optimized BFA O (N n + f Θ (N n )) O (N n ) Simplified BFA O (N n + f Θ (N n )) O (N n ) Depth-first (DFA) O (N n ) O (N ) Optimized DFA O (N n ) O (N ) Simplified DFA O (N n ) O (N ) CIFSAlgo O (N n ) O (N n ) Generalised IFS O (n • f Θ (N )) O (N )
We recall that f SOp (x) is the complexity of the operator SOp evaluated on the union of x objects, f Θ is the complexity of the operation Θ as a function of parameters count, N is the number of the IFS transformations and n is the required number of iterations. Optimized algorithms do not depend on the operator SOp complexity because the operation is evaluated only for a single object. The function f SOp (x) is therefore constant on the number of iterations n.

Chapter 7

Conclusion and prospects

Our research is conducted as part of the ANR ModItère project . The global objective of this project is to develop a new type of geometric modellers for computer-aided design systems (CAD), based on iterative methods following the principles of fractal geometry. Most existing CAD modellers are based on classic geometry. These modellers allow a user to model analytic curves and surfaces, which are usually polynomial or rational and smooth. Fractal shapes have more complex and specific geometric properties. Fractal shapes can generally be rough, lacunar, porous, they are not smooth and usually not differentiable in any way. It is therefore difficult or impossible to represent these shapes by using classic geometric models. To generalise the answer to this problem specific iterative models such as IFS, CIFS and BCIFS were developed, in which we are substantially interested. Defining shapes by iteration allows us to generate fractal structures with all the diversity of their specific properties, which cannot be achieved with classic modelling.

Shape modelling in CAD systems is complex, in that it is not confined to the simple shape representation. To assist the user, geometric modellers are improved by construction operations and algorithms. These operations are usually based on the geometric properties of the modelled objects. Existing CAD systems are not suitable to manipulate fractal structures because of their specific properties. After listing the constructing operations presented in most of CAD systems, we identify the basic properties required for their realization. A preliminary analysis shows that the operations are based on four properties: affine invariance, topological structure, parametrization and differential properties. The result is given in the form of a dependency graph.

From this graph, we classified CAD operations by adapting to BCIFS : certain operations require generalization because of the specific properties of fractals, for certain operations the result can be approximated. For approximating CAD operations, we thus introduced generic algorithms defined for arbitrary operations satisfying some constraints that we explicit.

The proof of convergence of these algorithms is based on domain theory. Most CAD operations are not continuous. We consider continuous extensions of CAD operations to the solid domain. Thus we represent the attractor in this domain, apply the continuous operation to the approximate attractor and then return to the modelling space (see figure 4.9).

The presented approximation algorithms are iterative. Domain theory allows us to formalize the notion of solution sets that represent the information on the solution at the current approximation level. On each iteration we compute a local solution set representing a part of the global solution. In the optimized algorithms we use a criterion determining whether a local solution set could potentially improve the global one. The global solution set is restricted up to the required accuracy while iterating the algorithm. We then identify cases in which this criterion can also be simplified to a comparison with a constant. Given a CAD operation which satisfies the image decomposition, it is possible to construct a CIFS with a generalised Hutchinson operator, whose attractor is close enough to the operation result with respect to the Hausdorff metric. Thus we introduce a generic algorithm to compute such CIFS for a given accuracy. We define the self-similarity property of the operation. This property defines a set of transformations, which are used in the output iterative system. Starting with the initial state, the algorithm iteratively traverses the evaluation tree of the input BCIFS. In order to provide a self-similar nature of the output set the algorithm verifies the similarity of the sub-attractors that were already parsed.

In order to construct an exact CIFS (with a generalized Hutchinson operator) of the operation, if it exists, we must prove all the necessary similarities manually. We also explicit the condition of the operation to be represented by an IFS with a generalised Hutchinson operator. In this case, only this condition should be proved manually.

We studied the complexity of all the presented algorithms and performed some experiments: a convex hull calculation, boolean set operations (intersection, union, difference) between fractal objects and classic CAD objects (lines, planes, splines, cylinders). A primary goal for our future works could be the integration of these approaches to our geometric modeller.

Several questions remain unanswered. As we described in section 4.3.4, there are restrictions on the choice of the initial partial solid. Constructing a suitable initial interior part is generally challenging, and we often have no choice except to work without approximating the interior. This makes it impossible to apply important operations such as the membership predicate and the subset inclusion.

Another question is related to representing the operation result. As described in chapter 5, when the self-similarity property (5.1), is satisfied, it is possible to construct a CIFS whose the attractor is close enough to the operation image. However, it is also interesting to derive the topological model of the output iterative system, that is to say to construct a BCIFS representing the operator image. For example, in applying a boolean intersection, some constraints of incidence and adjacency can be induced from the initial BCIFS constraints. However, the additional constraints for the new edges and faces somehow have to be identified.

Finally, a whole class of operations working with multiple attractors remains unexplored. Of course, for approximate computations we could replace the attractors with their approximations, and calculate the result using available techniques. However, the effectiveness of this approach may to be questioned. The key difficulty here is that we have to increase the level of approximation in all attractors in parallel, that is to say, to find a balance between the rates of their evolution. efficiency, this approach allows us to achieve better convergence as it computes the minimum bounding ball for a given accuracy.

A.2 Bounding ball for the attractor of a CIFS

Here we describe a method for computing a set of bounding balls for the attractor of a CIFS. For each CIFS attractor we thus find a bounding ball (not necessarily a tight one) centered in a point of the attractor.

A path in the CIFS graph is a sequence of the automaton states {q k } l+1 k=1 which satisfies the following conditions:

• ∀k = 1, • • • , l ∃i : q k+1 = δ(q k , i)

• q i = q j ∀i, j, i = j

Here l ∈ N is called length of the path and denotes a number of transformations applied to reach the state q l+1 . A simple cycle in the CIFS graph is a path which is closed on its first state:

• q 1 = q l+1

An operator of the simple cycle is a composite of the transformations T q k i between the cycle states. Let us denote the simple cycle from the state a (i.e. q 1 = q l+1 = a) by ξ a .

Finally, let θ a be a set of all simple cycles from the state a.

The CIFS system defines a family of attractors associated with the states. Since the attractors are mutually recursively defined, we construct the balls by depth-first searching in the CIFS control graph. Let us consider an example of a CIFS attractor, illustrated in figure A.1. The system is described by an automaton with the initial state a and several simple cycles. where T ξ a is an operator of the simple cycle ξ a and s ξ a is its contraction coefficient. A ball center c a is computed as a fixed point of T ξ a , where ξ a is any cycle in θ a . It is convenient to compute the radius r a for each simple cycle ξ a ∈ θ a separately and then take the maximum one. Let us consider another example of a CIFS attractor, illustrated in figure A.2. There is a path {q 1 = , q 2 , . . . , q l , q l+1 = a} to the several simple cycles ξ a i ∈ θ a . As described above, we firstly compute a ball B a using the simple cycle transformations. Secondly, we compute the balls B qi recursively verifying the following inclusions: B(c qi , r qi ) = B qi ⊇ T qi (B qi+1 ), ∀i = l, . . . , 1.

That is, we compute the radius r qi as follows: r qi = max(r qi , d(T qi (c qi+1 ), c qi ) + r qi+1 • s(T qi )), If the ball center c qi has not defined yet, we compute it as the image of c qi+1 .

The attractor of any CIFS is formed by schemes shown in figure A.2. Since each simple cycle is considered separately, we can then find all schemes by depthfirst searching in the CIFS control graph. For each simple cycle ξ a found we firstly compute the ball B a , and secondly we compute the balls through the path to this cycle. There is a path {q 1 = , q 2 , . . . , q l , q l+1 = a} to the several simple cycles ξ a i ∈ θ a .

To perform a depth-first search in the control graph, we use a recursive function. The algorithm starts off by calling this function for the initial state . When the whole graph is gone through, all the required radii and centres will be defined.

Given a CIFS with an automaton (Σ, Q, δ) and an initial state . Here is the described initialization method implemented on pseudocode:

Résumé :

La définition de formes par ces procédés itératifs génère des structures avec des propriétés spécifiques intéressantes : rugosité, lacunarité. . . . Cependant, les modèles géométriques classiques ne sont pas adaptés à la description de ces formes. Dans le but de développer un modeleur itératif pour concevoir des objets fractals décrits à l'aide du BCIFS, nous avons développé un ensemble d'outils et d'algorithmes génériques qui nous permettent d'évaluer, de caractériser et d'analyser les différentes propriétés géométriques (la localisation, le calcul de l'enveloppe convexe, de la distance à partir d'un point, etc) de fractals. Nous avons identifié les propriétés des opérations standards (intersection, union, offset, . . . ) permettant de calculer une approximation d'image des fractales et de plus d'optimiser ces algorithmes d'approximation. Dans certains cas, il est possible de construire un CIFS avec l'opérateur de HUTCHINSON généralisé dont l'attracteur est suffisamment proche du résultat de l'opération par rapport à la métrique de Hausdorff. Nous avons développé un algorithme générique pour calculer ces CIFS pour une précision donnée. Nous avons défini la propriété d'auto-similarité de l'opération, qui définie un ensemble de transformations utilisé dans un système itératif résultant. Pour construire un CIFS exact de l'image, si il existe, il faut prouver tous les similitudes nécessaires manuellement. Nous explicitons également la condition de l'opération, quand le résultat peut être représenté par un IFS avec un opérateur de HUTCHINSON généralisé. Dans ce cas, il n'est que cette condition à prouver manuellement.

Mots-clés : Fractal, modélisation géométrique, conception assistée par ordinateur, géométrie algorithmique, informatique graphique

Abstract:

Defining shapes by iteration allows us to generate new structures with specific properties (roughness, lacunarity), which cannot be achieved with classic modelling.

For developing an iterative modeller to design fractals described by a BCIFS, we developed a set of tools and algorithms that permits one to evaluate, to characterize and to analyse different geometric properties (localisation, convex hull, volume, fractal dimension) of fractals. We identified properties of standard CAD operations (intersection, union, offset, . . . ) allowing us to approximate them for fractals and also to optimize these approximation algorithms.

In some cases, it is possible to construct a CIFS with generalised HUTCHINSON operator, whose attractor is close enough to the operation result with respect to the HAUSDORFF metric. We introduced a generic algorithm to compute such CIFS for a given accuracy. We defined the self-similarity property of the operation defining a set of transformations, which are used in the output iterative system. In order to construct an exact CIFS of the image, if it exists, we must prove all the necessary similarities manually. We explicit also the condition of the operation to be represented by an IFS with a generalised HUTCHINSON operator. In this case, only this condition should be proved manually.

Keywords: fractal, geometric modelling, computer-aided design, computational geometry, computer graphics

Figure 1

 1 Figure 1.1: Several applications of fractal structures. (a) Example of application in architecture, EPFL 2010 (b) The study of acoustical properties of fractal structures (Emmanuel Redon LVA LYON I) (c) Examples of fractal shapes fabrication (Project ANR ModItre).

Figure 2 .

 2 Figure 2.1: A 2D surface swept along a curve

Figure 2 . 2 :

 22 Figure 2.2: Complex objects are created by adding or subtracting predefined solid primitives. For combining these primitives one usually uses basic set operators, called Boolean operators.

Figure 2 .

 2 Figure 2.3: Free-form surfaces are generally created by manipulating the surface control points

Figure 2 . 4 :

 24 Figure 2.4: Graph of the relations between the CAD operations and the object properties.

  a∈A inf b∈B d(a, b), sup b∈B inf a∈A d(a, b)}.

Figure 3

 3 Figure 3.1: Sierpinski triangle

Figure 3

 3 Figure 3.4: The IFS evaluation tree calculated to the third level. Internal nodes correspond to the calculation of a composite function. Leaves correspond to subsets of K 3 to construct or to visualise.

  Figure 3.5: CIFS evaluation tree calculated to the third level. Internal nodes correspond to the calculation of a composite function. Leaves correspond to subsets of K 3 to construct or to visualise.

Figure 3 . 6 :

 36 Figure 3.6: Attractor of a CIFS described by an automaton with two states: a and b.

Figure 3 . 7 :

 37 Figure 3.7: The automaton of the CIFS considered in the example. There are two subdividing operators from the state a to construct an array of the Sierpinski triangles and three operators from b to construct each triangle.

Figure 3 .

 3 Figure 3.10: 4D barycentric space X C for constructing a square. Three subdividing operators transform the space into itself by scaling with a ratio 0.5. The fourth operator projects the Sierpinski triangle constructed in the space X S .

  Figure 3.11 shows the automaton of this CIFS. Transition functions are the following:

Figure 3 .

 3 Figure 3.12: Attractor of the CIFS considered in the previous example with a different projecting map P .

Figure 3 .

 3 Figure 3.13: Subdivision of a face in the considered example of a BCIFS.

Figure 3 .

 3 Figure 3.14: Subdivision of an edge in the considered example of a BCIFS.

Figure 3 .

 3 Figure 3.15: Decomposition of the square into the topological cells (face, edge, vertex).

Figure 3 .

 3 Figure 3.16: All the incidence and adjacency constraints for a square subdivided into four similar squares.

Figure 4

 4 Figure 4.1: Construction steps for the Koch curve (from left to right, and from top to bottom). The curve is constructed by starting with a single line segment, and by repeating the following three steps recursively: (1) divide the line segment into three segments of equal length, (2) draw an equilateral triangle that has the middle segment obtained from step 1 as its base and points outward, (3) remove the line segment that is the base of the triangle from step 2.

  Figure 4.2 illustrates this process.

Figure 4 . 2 :

 42 Figure 4.2: Approximation of the attractor illustrated by a tree whose leaves are the images of an initial compact subset B by applying sequences of IFS transformations.

Figure 4 .

 4 Figure 4.3 illustrates two simple examples of partially ordered sets.

Figure 4

 4 Figure 4.3: Two simple partially ordered sets. The pointed domain of booleans {+, -} ⊥ is illustrated on the left, and an increasing chain on the right.

Figure 4

 4 Figure 4.4: In the theory of information, a directed set corresponds to a consistent set of inputs and outputs of a given program.

  Definition. Let (P, ) be a dcpo. A set O ⊆ P is open iff: 1. x ∈ O and x y ⇒ y ∈ O; 2. whenever D ⊂ P is directed set and D ∈ O then D ∩ O = ∅. The first condition thus introduces the nature of open sets in the Scott topology and the second condition means that the border of the set O is not included.

x∈DF

  (x) = F ( x∈D x),for all directed subsets D ⊆ P .

Figure 4 . 5 :

 45 Figure 4.5: Plot of the function f considered in the example.

Figure 4

 4 Figure 4.6: In the standard topology this function is not continuous, but lower semi-continuous.

  [a, b] and [c, d] are intervals. One can show that [a, b] [c, d] ⇔ (a, b) ⊇ [c, d]. Indeed, if the intervals share the left border a, there exists a directed set D = {[a -1 n , a] | n ∈ N} with the lub a. There is no element z ∈ D such that

Figure 4

 4 Figure 4.7: Lower semi-continuous functions are continuous in the Scott topology.

Figure 4 . 8 :

 48 Figure 4.8: Partial solid (A, B). The sets A and B are intended to capture, respectively, the interior and the exterior of a solid, possibly, at some finite stage of computation.

Figure 4

 4 Figure 4.9: We consider a continuous extension of the operator SOp on an ωcontinuous dcpo (D, ). Thus we construct a representation of the attractor in this domain, apply the continuous operation to it and then return to the modelling space X by approximating the attractor with the accuracy ε.

Figure 4 .

 4 Figure 4.10: Attractor A represented in the solid domain by a pair of open sets(∅, A c). We use the empty set as an interior approximation.

  which implies that the set A 0 must also be the empty set.Thus we identified the lub of the sequence SB n to be the partial solid (D, A c ) SA, where D = A. The only possibility for SB n to have the lub SA is to choose the initial open set A 0 = A o , which corresponds to the case D = A o . In fact, another set D prevents the convergence SB n → SA with respect to the Scott topology.Consider the following example. Let {X; T i | i = 0, . . . , 3} be an IFS subdividing square in four parts, as illustrated in figure4.12.

Figure 4 .

 4 Figure 4.12: Subdivision of a square by four transformations. Let A = A 0 ⊆ A o be an open set and B = B 0 be a bounding ball of this square such that A ⊆ T(B) ⊆ B, as illustrated in figure 4.13. On the following iteration we apply the Hutchinson operator to this partial solid (see figure 4.14). The resulting set is no longer connected, since the images of a border appear in the middle of the square.

Figure 4 .

 4 Figure 4.13: Initial partial solid (A 0 , B c 0 ) considered in the example.

Figure 4 .

 4 Figure 4.14: The first iteration partial solid (A 1 , B c 1 ) considered in the example.

Figure 4 .

 4 Figure 4.15: The n-th iteration partial solid (A n , B c n ) considered in the example.

Figure 4 .

 4 Figure 4.16: Several elements of the sequence SB n in the considered example. The interior and the exterior parts are in blue and red, respectively.

Figure 4 .

 4 Figure 4.17: Example of the directed set constructed by applying an IFS with condensation. Condensation set is the circle in the center. The interior and the exterior parts are in blue and red, respectively.

Statement.

  If there exists an associative operation Θ defined on SY as follows: SOp(A ∪ B) = SOp(A) Θ SOp(B), (4.9) the image of SB n = T n (SB) can be obtained by composing the images of each leaf, i.e.: SOp( i T i (SB)) = Θ i SOp(T i (SB)).

Figure 4 .

 4 Figure 4.18: Calculation of the approximate operation result by applying the operation to each leaf of the tree. We then compose the results by means of an associative binary operation Θ.

Figure 4 .

 4 Figure 4.19: Modifying the evaluation tree allows us to optimize the generic algorithms. We apply the operator SOp on intermediate levels of the tree to prevent unnecessary calculations as early as possible.

Figure 4 .

 4 Figure 4.21: Simplification of the stopping criterion. We also apply the operator SOp on the intermediate levels of the tree, but we compare the result with the empty partial solid (∅, Y) only. Thus we avoid constructing the global solution set S.

Figure 5 . 2 :

 52 Figure 5.2: Example of an intersection between a CIFS and a line L. The system is described by an automaton with three states: , C and S. The image illustrates attractor A in the 2D modelling space.

  which implies: P (SOp P (SB)) = P (SB ∩ L ) = P (SB) ∩ L = SOp(P (SB)).

Figure 5

 5 Figure 5.3: Instead of intersecting attractor A with line L in the 2D modelling space, now we have to intersect the attractor A C with the plane L in the 4D barycentric space. Plane L is then projected onto L by the operator P .

Figure 5 . 4 :

 54 Figure 5.4: Intersection between a CIFS and a line L. The output CIFS automaton consists of three states, we denote these states by the corresponding calculating sets A ∩ L, A C ∩ L and A S ∩ L . The intersection between the line and the Sierpinski triangle is thus computed in the state A S ∩ L , it is then placed in the 4D barycentric space according to the control points, and finally projected onto the modelling space.

Figure 5

 5 Figure 5.6: CIFS automaton constructed by the presented algorithm. Each node corresponds to a state of the automaton with the associated CIFS attractor. Attractors are defined by the generalised Hutchinson operator. Directed edges correspond to the transitions and they are marked by the associated transformations. The initial state of the automaton is associated with the attractor SOp(SA).

  Figure 5.6 illustrates an example of a CIFS automaton produced by this algorithm. The following function illustrates the algorithm implementation: Function CIFSAlgo Input: SOp is the operator required to evaluate γ is the calculation tolerance Output: a CIFS with the attractor close to SOp(SA) Body: Q ← [(Id, )] // an output set of CIFS states and composite transformations F ← [] // an output set of CIFS transitions S ← [(Id, )] // a temporary list of states and composite transformations while (S is not empty)

Figure 5 . 7 :

 57 Figure 5.7: Intersecting the Sierpinski triangle with a line L in the 3D barycentric space.

Figure 5 . 8 :

 58 Figure 5.8: Intersecting the Sierpinski triangle with the lines T -1 i (L) for i ∈ Σ.

Figure 5 . 9 :

 59 Figure 5.9: Automaton of the CIFS considered in the example. The attractor of this CIFS is exactly the set SL ∩ SA.

Figure 5 .

 5 10 illustrates another example of intersection between the Sierpinski triangle and different lines.

Figure 5 .

 5 Figure 5.10: Example of intersection between the Sierpinski triangle and different lines.

Figure 6 . 1 :

 61 Figure 6.1: Class diagram for the presented solution.

  b∈B d(p, b), +∞), (0, min b∈B d(p, b))). The values min b∈B d(p, b) and max b∈B d(p, b) correspondingly indicate the lower and upper bounds for the operator result SOp(SA).

Figure 6 . 2 :

 62 Figure 6.2: Initial solution set. The lower and upper bounds of the solution are computed as the minimal and maximal distances to the attractor bounding ball. The exact result SOp(SA) is thus situated somewhere in between.

  SOp(T 0 (SB)) = ((d 0 + r 0 , +∞), (0, d 0r 0 )), where d 0 = d(p, T 0 (c)), r 0 = r • s 0 and s 0 is the contraction coefficient of T 0 . The stopping criterion obviously does not return +. We therefore specify the global solution set S by composing with SOp(T 0 (SB)):S ← S Θ SOp(T 0 (SB)) = = ((min{d + r, d 0 + r 0 }, +∞), (0, min{dr, d 0r 0 }))= ((d 0 + r 0 , +∞), (0, dr)),

2 ,

 2 we represent the solution sets in SY by pairs of the following form: ((a, +∞), (0, b)), where a, b ∈ R are two real values. The following pseudo-code illustrates how we define the solids of SY: Object DistanceSetY implements SetY Body: in ∈ R is the interior part of the partial solid ex ∈ R is the exterior part of the partial solid Function GetAccuracy Output: An accuracy of a given partial solid. Body: return |in -ex| End Function Inclusion predicate ⊆ Input: A, B are objects of type DistanceSetY. Output: A ⊆ B. B.in A.ex, -if A.in B.ex, ⊥ otherwise. End End Thus we defined the nature of solids in SY, and the operator SOp itself can now be implemented. The following function illustrates the implementation of the operator SOp: Function OP Input: SB is the input partial solid in SX representing a transformed ball Output: The distance from the point p ∈ X to the partial solid SB Body: dist ← |p -SB.center| return DistanceSetY with in = dist + SB.radius ex = max{0, dist -SB.radius} End

  Function Θ Input: A, B are objects of the type DistanceSetY. Output: A Θ B Body: return DistanceSetY with in = min{A.in, B.in} ex = min{A.ex, B.ex} End

Figure 6

 6 Figure 6.3: Iterative construction of approximations.Given an ε iapproximation of the convex hull, we calculate its images by applying all the IFS transformations to its vertices. We then merge all these images to obtain the ε i+1 -approximation of the convex hull.

Figure 6

 6 Figure 6.4: Illustration of the simplification. If two neighboring points c and c are close, i.e. d(c, c ) + r c + r c 2ε i we can then replace these points with the middle point without any loss of accuracy.

Figure 6 . 7 :

 67 Figure 6.7: Subdivision of the intersection between a plane and the Menger sponge.

Figure 6 .

 6 Figure 6.7 shows an example of intersection between a plane and the Menger

Figure 6 . 9 :

 69 Figure 6.9: Counterexample where the offsetting on a negative distance cannot be decomposed by Θ operator.

  B) ⊕ B(0, r) = {T (x) + b | x ∈ B, b ∈ B(0, r)} = = {T (x) + T (b ) | x ∈ B, b ∈ B } = = {T (x + b ) + v | x ∈ B, b ∈ B } = = { T (x + b ) | x ∈ B, b ∈ B } = = T (B ⊕ B ).

Figure A. 1 :

 1 Figure A.1: An example of a CIFS. The system is described by an automaton with the initial state a and several simple cyclesThe attractor A a is defined as follows:

Figure A. 2 :

 2 Figure A.2: An automaton of the CIFS.There is a path {q 1 = , q 2 , . . . , q l , q l+1 = a} to the several simple cycles ξ a i ∈ θ a .

  is an initial partial solid such that T(B) ⊆ B SOp is the operator required to evaluate S is the solution set Input: n is a number of levels to pass Output: an approximate operator result SOp(T n (SB))

	Body:
	ts ← [Id] // A set of composite IFS transformations
	for k = 1 to n do

ts ← [] // Build a new set of composite IFS transformations S ← (∅, Y) // Reset S to avoid the comparison with the precedent tree level foreach T ∈ ts such that SOp(T (SB)) S // Stopping criterion (4.12) S ← S Θ SOp(T (SB)) ts .push(T T i ), ∀i ∈ Σ endfor ts ← ts endfor return Θ T ∈ts

Table 1 .

 1 . Execution time of our algorithms for several IFS attractors with the precision ε.

	Attractor	N	s max	dim	ε	time(ms)
	Sierpinski triangle von Koch curve Dragon curve Tree Sierpinski tetrahe-	3 4 2 5 4	0.5 0.25 0.70711 0.66762 0.5	2 2 2 2 3	10 -9 10 -9 10 -9 10 -9 10 -9	16 27 62 94 109
	dron
	IFS dragon Barnsley fern FIF 3D Barnsley fern	2 4 4 4	0.91944 0.85094 0.86603 0.85586	2 2 3 3	10 -9 10 -9 10 -3 10 -3	265 484 3541 10951

  Table 2 demonstrates the execution time of our algorithm compared to Martyn's for 2D IFS attractors.

	Attractor	N	s max	Martyn's method	Our algorithm
				time, ms	time, ms
	von Koch curve	4	0.25	125	16
	Sierpinski triangle 3	0.5	647	16
	Dragon curve	2	0.70711	756	16
	Tree	5	0.66762	1131	31
	Barnsley fern	4	0.85094	45630	250
	IFS dragon	2	0.91944	more than 60 sec	78
	Table
				Fractal interpolation function
	3D Barnsley fern		3D Barnsley fern
	Figure 6.6: Results for 3D IFS attractors at precision ε = 10 -2
			97

we start with the initial state and then loop through all the IFS transformations checking which of the alternatives the corresponding sub-attractor hits:

• the sub-attractor is empty. We do not add any special state to construct the empty set, so we continue the iteration.

• the sub-attractor is similar to one that was already parsed, denotes s . In this case, we close the transition to the state s .

• otherwise, we should continue the iteration recursively.

The algorithm thus finds the similar parts in the output set and constructs a CIFS automaton with the corresponding attractor.

Figure 5.5: The key idea of the algorithm is to apply a classic evaluation and verify the similarities between constructed sub-attractors. If the sub-attractor is similar to one that was already parsed, we process these sub-attractors as the same set and unite the corresponding states.

If the operator SOp verifies the self-similarity property (see equality 5.2) we can construct a specific CIFS with the attractor SOp(SA) as follows. We start with the initial state associated to the attractor SOp(SA) construction. We then add the transitions T i for all i ∈ Σ from the initial state to the new states associated to the corresponding attractors SOp Ti (SA) and continue this process recursively.

If a partial solid SOp Ti (SB) is empty, that is equal to (∅, Y), then it is not necessary to add a new state, because such state is associated with the empty attractor.

Consider the k-th iteration of this recursion. Let consider two composite transformations given by two sets of indices:

Here k and l, l k represent the lengths of these composites.

Implementation

In this section we present an implementation of the prerequisites for the operator to compute the intersection with a given line L. Let the line L be defined by a point L.ref and a vector L.dir representing the reference point on the line and the direction correspondingly. As for the distance example in section 6.2.1, the first function which has to be implemented determines the initial partial solid as a bounding ball B q of the attractor A q for a given state q ∈ Q using the algorithm, presented in appendix A. The initial partial solids SB q are thus the pairs (∅, B q c ).

Solution sets in SY are represented by a set of segments on the line L as the exterior part and the empty interior part. The following pseudo-code illustrates how to define partial solids in SY:

Object IntersectionSetY implements SetY Body: segments is the list of segments Function GetAccuracy Output: an accuracy of a given partial solid. Body: // We do not use this stopping criteria return ∞ End // We use the simplified stopping criterion (see section 4.5.1).

Function Inclusion predicate ⊆ Not applicable

End

Thus we defined the nature of partial solids in SY, and the operator SOp itself can now be implemented. The following function illustrates the operator SOp implementing:

The decomposing operator Θ is simply defined for two partial solids in SY by the following function: (c,d) intersect with each other we union them, i.e. (min{a, c}, max{b, d}). The segments, which do not intersect, may be left unchanged.

In the breadth first algorithm the operator Θ could also compose a whole list of the intermediate results, i.e. Θ : SY k → SY, where k is the number of input results. We simply apply the defined operator Θ to the first two elements in the list, and then we recursively compose the result with each the following element.

Finally, after computing the approximate partial solid SOp(SA) in the form of the IntersectionSetY object, the output result can be obtained by using its segments:

Function BuildSolution Input:

A is the objects of the type IntersectionSetY approximating SOp(SA). Output: an approximate operator image Op(A) Body: return A.segments End As described in section 5.2.2, to determine similar sub-attractors we verify the inequality (5.6). It implies the inequality (5.4) and does not contain any applications of the operator SOp, that is to say that we verify it only if two transformed lines are close enough.

These are all the methods which have to be implemented to apply the operator using the presented algorithm. Moreover, for this example we can construct a CIFS with the attractor that is close to the solution SOp(SA).

Offset

The offsetting map produces a shape driven out from an original shape in equal distance and direction, normally or perpendicularly to the original shape. One usually distinguishes between the offsetting on a positive distance and on a negative distance. In this section we consider both offsetting maps and describe whether the maps could be computed by our approach.

Consider first the operator Op : X → X defined through the Minkowski sum as follows:

Appendices

Appendix A

Computing the attractor bounding ball

A.1 Bounding ball for the attractor of an IFS Methods to determine a bounding ball for the attractor of an IFS have already been developed. Gentil [START_REF] Gentil | Les fractales en synthése d'images: le modèle IFS[END_REF] described an approach based on the dichotomous search for the minimal radius of a ball that bounds the attractor of an IFS. The approach can be applied in a multi-dimensional space and calculates the result for a given precision.

Hart and DeFanti [START_REF] Hart | Efficient anti-aliased rendering of 3d linear fractals[END_REF] introduced a method which starts with the unit ball centered at the origin. The algorithm iteratively produces a sequence of balls converging to the limit ball that bounds the attractor.

Rice [START_REF] Rice | Spatial bounding of self-affine iterated function system attractor sets[END_REF] improved on Hart and DeFanti's approach by optimizing the radius of the bounding ball using a generic optimization package. He also showed that the centre of the limit ball can be determined analytically by solving a system of linear equations.

Martyn [START_REF] Martyn | Tight bounding ball for affine ifs attractor[END_REF] showed that the solution of this system is the centroid of the attractor with particular weights. To obtain a better approximation, he presented a heuristic iterative method called "balancing the attractor". The algorithm is not limited by the dimension of the space in which the attractor lies.

More recently, Martyn [START_REF] Martyn | The smallest enclosing disc of an affine ifs fractal[END_REF] presented a novel approach to approximate the smallest disc to enclose the affine IFS attractor at any accuracy. The method is based on the concept of spanning points he introduced to describe the extent of an IFS attractor.

In this thesis, we require the initial compact subset B to be not only the bounding ball (not necessarily a tight one), but also the set, for which the following inclusion holds:

T(B) ⊆ B, because with this property we have the continuous approximation of the attractor exterior. Each element of this sequence is included in the next one.

We have chosen the method suggested in [START_REF] Gentil | Les fractales en synthése d'images: le modèle IFS[END_REF] because it satisfies all the requirements and also because of it is simple to implement. In addition to // To perform a depth-first search in the automaton's // graph, we use a recursive function. Algorithm starts // off by calling this function with q = and empty lists Function CheckState Input: q is a current state q ∈ Q states is a list of traversed states t is a list of applied transformations Output: initial bounding balls B q = B(c q , r q ) ∀q ∈ Q Body: { // We are looking for a cycle index ← states.indexof(q)