Ekkehard Ramm

F Gauss

P.-L M De Maupertuis

Principles of Least Action and of Least Constraint

Keywords: Least Action, Principle of Least Constraint, Method of Least Squares, C

The present contribution describes the evolution of two major extremum principles in me-chanics proposed in the 18 th and the first half of the 19 th century, namely the Principle of Least Action associated with the name Maupertuis and Gauss's Principle of Least Constraint. We will briefly mention the d'Alembert Principle strongly related to the Principle of Gauss.

Preliminary Remark

Like the other contributions in this special issue the present paper is also based on a lecture given in a session on the history of mechanics at the GAMM annual meeting in 2010. The author not being a historian apologizes for a rather short and partly superficial essay on both principles, based mainly on secondary sources. It can by far not replace any of the excellent rigorous and detailed descriptions on the subject, some of them mentioned in the list of references. This holds in particular for the description of the scientific controversy on the Principle of Least Action, probably one of the most documented affairs in the history of sciences.

Introduction

The search for the existence of extremum principles in nature and technology can be traced back to the ancient times. The scholars were not only driven by scientific observations; often their objective was based on metaphysical arguments. In Figure 1 major figures involved in the development of extremum principles in mechanics are arranged in the respective time course. In the following essay the era of the mid 18 th and the early 19 th century is selected; in particular we will concentrate on Pierre-Louis Moreau de Maupertuis related to the Principle of Least Action in 1744/46 and on Carl Friedrich Gauss who stated the Principle of Least Constraint in 1829. 2 Principle of Least Action

Maupertuis and his Principle

The discussion on the Principle of Least Action is one of the most extensive examinations of a scientific controversy. The protagonists (Figure 2) are Maupertuis and Euler on one side and Fig. 2 The Protagonists. König and Voltaire on the other side. Leibniz' work also played a role in the fight along with Frederick the Great, King of Prussia, as a leading authority at that time in Berlin.

The French mathematician Pierre-Louis Moreau de Maupertuis (Figure 3) was very much influenced by the work of Newton. Already in papers to the French Academy of Sciences in 1741 and 1744 he mentioned a principle minimizing a quantity which he called action. After he became the President of the Prussian Academy of Sciences in Berlin, on invitation of King Frederick the Great in 1746, he presented the book "Les Loix du Movement et du Repos" (The Laws of Movement and of Rest derived from a Metaphysical Principle) [START_REF] De Maupertuis | Les Loix du Mouvement et du Repos déduites d'un Principe Metaphysique[END_REF], see Figure 4, left. In the introduction he points to his previous work at the Paris Academy in 1744 Fig. 4 Maupertuis' Les Loix (1746) [START_REF] De Maupertuis | Les Loix du Mouvement et du Repos déduites d'un Principe Metaphysique[END_REF] and the Principe de la Moindre Action (Principle of Least Action). Fig. 5 Euler's Methodus Inveniendi Lineas Curvas (1744) [START_REF] Euler | Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes[END_REF] and Supplement II. and adds the remark "At the end of the same year, Professor Euler published his excellent book Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes" (A method for finding curved lines enjoying properties of maximum or minimum) [START_REF] Euler | Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes[END_REF], Figure 5. He then continues "In a Supplement to his book, this illustrious Geometer showed that, in the trajectory of a particle acted on by a central force, the velocity multiplied by the line element of the trajectory is always a minimum". In other words Maupertuis was aware of maximum and minimum properties however he reduced his further considerations to a minimum principle.

In his "Les Loix. . . " Maupertuis first critized the usual proofs of the existence of God and refers to the fundamental laws of nature. His self-confident statements are remarkable when he writes "After so many great men have worked on this subject, I almost do not dare to say that I have discovered the universal principle upon which all these laws are based. This is the principle of least action, a principle so wise and so worthy of the Supreme Being, and intrinsic to all natural phenomena; one observes it at work not only in every change, but also in every constancy that Nature exhibits".

Finally he states the general principle as follows (Figure 4, right): "When a change occurs in Nature, the Quantity of Action necessary for that change is as small as possible", and defines "The Quantity of Action is the product of the Mass of Bodies times their velocity and the distance they travel. When a Body is transported from one place to another, the Action is proportional to the Mass of the Body, to its velocity and to the distance over which it is transported". Expressing the Action A in a formula yields

A ∼ M • v • ds (1)
With the mass of a particle M , the velocity v and the transported distance ds. The relation of the Action to the kinetic energy is apparent if the distance ds is replaced by velocity times the time increment dt.

In the sequel of this definition Maupertuis adds three examples as a proof for the generality of his principle, namely on the laws of motion for inelastic as well as elastic bodies and the law of mechanical equilibrium. In the mentioned Supplement II "De motu projectorum. . . " (On a motion of particles in a non-resistant medium, determined by a Method of maxima and minima), Figure 5, right, Leonard Euler (Figure 6) says "1. Since all natural phenomena obey a certain maximum or minimum law, there is no doubt that some property must be maximized or minimized in the trajectories of particles acted upon by external forces. However, it does not seem easy to determine which property is minimized from metaphysical principles known a priori. Yet if the trajectories can be determined by a direct method, the property being minimized or maximized by these trajectories can be determined, provided that sufficient care is taken. After considering the effects of external forces and the movements they generate, it seems most consistent with experience to assert that the integrated momentum, i.e. the sum of all momenta contained in the particle's movement, is the minimized quantity".

In section 2 of the Additamentum Euler defines mass M , velocity √ v (he used √ v instead of v), infinitesimal distance ds, the momentum integrated over the distance ds as M • ds • √ v and writes: "Now I assert that the true trajectory of the moving particle is the trajectory to be described, from among all possible trajectories connecting the same end point, that minimizes M ds √ v or since M is constant ds √ v. Since the velocity √ v resulting from the external forces can be calculated a posteriori from the trajectory itself, a method of maxima and minima should suffice to determine the trajectory a priori. The minimized integral can be expressed in terms of the momentum (as above), but also in terms of the living forces (kinetic energies). For, given an infinitesimal time dt during which the element ds is traversed, we have ds = dt √ v. Hence ds √ v = vdt, i.e. the true trajectory of a moving particle minimizes the integral over time of its instantaneous living forces (kinetic energies). Thus, this minimum principle should appeal both to those who favor momentum for mechanics calculations and to those who favor living forces".

Euler realized that for the entire path the action has to be a sum along all segments ds and essentially defined as quantity (velocity defined as v as usual)

A = Mvds = Mv 2 dt (2)
which is equal to the kinetic energy up to the factor 1 2 . Thus it turned out that Maupertuis' pretension on universality of his principle was wrong, confer [START_REF] Szabo | Geschichte der mechanischen Prinzipien und ihrer wichtigsten Anwendungen[END_REF][START_REF] Knobloch | Das große Spargesetz der Natur: Zur Tragikkomödie zwischen Euler, Voltaire und Maupertuis[END_REF]; the examples were not well chosen and partially not correct. His claim that the minimization is an economy principle did not hold because some problems lead to a maximum. Chevalier P. D'Arcy reproached Maupertuis for this error; for example the refraction of a concave mirror is based on a maximum property, [START_REF] Szabo | Geschichte der mechanischen Prinzipien und ihrer wichtigsten Anwendungen[END_REF]. Consequently his interpretation of Euler's statements was not correct.

The critics

It was his fellow student and friend Johann Samuel König (Figure 7) who became Maupertuis' strongest critic. Two years after he was appointed member of the Berlin Academy, strongly supported by the president, the Swiss lawyer and later mathematician König (1712-1757) passed an essay to Maupertuis who approved it for publication without reading. In his article König pled for a principle of extremal energy from which an extremum for the action can be derived. The essay appeared in Nova Acta Eruditorum 1751 [START_REF] Koenigio | De universali principio aequilibrii et motus, in Vi viva reperto, deque nexu inter Vim vivam et Actionem, utriusque minimo, Dissertatio[END_REF]. When Maupertuis read it he became upset because not only the limitation of his principle was questioned but it also claimed that G.W. Leibniz (1646Leibniz ( -1716) ) has given a more precise formulation already in a letter to the Swiss mathematician and theologian Jakob Hermann (1678-1733) in 1707. At the end of the paper [START_REF] Koenigio | De universali principio aequilibrii et motus, in Vi viva reperto, deque nexu inter Vim vivam et Actionem, utriusque minimo, Dissertatio[END_REF] he quotes from Leibniz' letter to Hermann [START_REF] Leibniz | Appel au Public, du Ingement de L'[END_REF] (see Figure 8) "But action is in no way what you think, there the consideration of the time enters; it is the product of the mass and the time 1 ; or of the time by the living force. I have pointed out that in the modifications of movements, it is usually derived as a Maximum, or a Minimum. From this can be deduced several important propositions; it can be used to determine the curves describing the bodies that are attached to one or several centres".

According to König the "force vive" (vis viva) has been already introduced by Leibniz expressing M • v 2 , i.e. it is two times the kinetic energy. It was not his intention though to accuse Maupertuis being a plagiarist or to disparage his merits; he simply wanted to point out the more general formulation of Leibniz (minimum or maximum). Nevertheless it was the beginning of intensive mutual disparagements ending in one of the ugliest of all scientific disputes [START_REF] Harnack | Geschichte der Koniglich Preussischen Akademie der Wissenschaften, Band 1.1: Von der Grundung bis zum Tod Friedrich des Großen[END_REF].

The letter König presented was a copy, a fact which would play a key role in the subsequent vehement quarrel. All attempts to find the original letter ended only in further copies; König claimed that he got his copy from the Swiss poet and politician Samuel Henzi in Bern who collected letters of Leibniz among other things. Henzi was decapitated in 1749 because of conspiration and all his original papers were supposed to be burned. The Berlin Academy met in spring 1752 charging König of forgery. It is noteworthy that the meeting was headed by Euler who from the very beginning was a strong supporter of Maupertuis despite his own more rigorous work on the same principle in 1744 (Supplement II, Figure 5); the reason for Euler's reaction is still a matter of conjectures; see e.g. [START_REF] Szabo | Geschichte der mechanischen Prinzipien und ihrer wichtigsten Anwendungen[END_REF][START_REF] Knobloch | Das große Spargesetz der Natur: Zur Tragikkomödie zwischen Euler, Voltaire und Maupertuis[END_REF][START_REF] Ariga | Science and its Public Perception: The Principle of Least Action in Eighteenth-Century Europe[END_REF]]. König's appeal was rejected; he resigned as member of the Academy in summer 1752.

King Frederick used all his authority to support Maupertuis, the President of "his" Academy. It was the time when another figure entered the scene, namely Voltaire (Figure 9), the famous French writer. Voltaire liked satires and polemics and was known for his sharp, malicious, sometimes also witty remarks. Although he was invited by Frederick to come to Berlin he took part of König; he published an anonymous letter "A reply from an Academician of Berlin to an Academician of Paris" in fall 1752 defending König (Figure 13, right); he refers to Monsieur Moreau de Maupertuis and says: "He asserts that in all possible cases, Action is always a Minimum, which has been demonstrated false; and he says he discovered this law of Minimum, what is not less false. Mr. Koenig, as well as other Mathematicians, wrote against this strange assertion, & he cited among other things, a fragment of Leibnitz in which the great man has remarked, that in modifications of movement, the action usually becomes either a Maximum, or a Minimum". König once more referred to copies of Leibniz' letter in a further pamphlet. Frederick was provoked and defended Maupertuis in a "Letter of an Academician of Berlin to an Academician in Paris" (Figure 13, right).

Voltaire published further polemics, the most famous being the "Diatribe due Docteur Akakia, Médicin du Pape" (Figure 10, left), which was a defamatory piece of writing. The Greek word Akakia means "without guile". The Diatribe printed in Potsdam was burned by Le Président et M. Le Professeur, 1st January 1753 [START_REF] Voltaire | Histoire du Docteur Akakia et du Natif de St[END_REF] was included saying among other things: "In future we promise, not to put the Germans down and admit that the Copernicus's, the Kepler's, the Leibnitz's. . . , are something, and that we have studied under the Bernoulli's, and shall study again; and that, finally, Professor Euler, who was very anxious to serve us as a lieutenant, is a very great geometer who has supported our principle with formulae which we have been quite unable to understand, but which those who do understand have assured us they are full of genius, like the published works of the professor referred to, our lieutenant", [START_REF] Dugas | A History of Mechanics[END_REF].

Thus Voltaire also included Euler in this dispute and holds him up to ridicule. Euler was indeed a vehement advocate of the Academy President despite his own more rigorous exposition of the subject. Dugas refers in [START_REF] Dugas | A History of Mechanics[END_REF] to an essay of Euler and quotes: "This great geometer has not only established the principle more firmly than I had done but his method, more ubiquitous and penetrating than mine, has discovered consequences that I had not obtained. After so many vested interests in the principle itself, he has shown, with the same evidence, that I was the only one to whom the discovery could be attributed".

Euler discussed his view of the development of the principle in [START_REF] Euler | Sur la principe de la moindre action[END_REF], presented already in 1751 but not printed in the Histoires of the Prussian Academy until February 1753 (Figure 11). He referred to König's role and made a strong statement: "But there is no one with whom we would be less likely to have dispute than Professor Koenig, who boldly denies that there is in nature such a universal law and pushes this absurdity to the point of mocking the Principle of Conservation, which constitutes the Minimum that nature appoints. In addition, he introduces the great Leibnitz as claiming and explaining that he was far from knowing such a principle. From this we see that Mr. Koenig cannot deny our President the discovery of the principle that he himself considers to be false". Finally he said: "The principle that Mr. de Maupertuis discovered is therefore worthy of the greatest of praise; and without a doubt it is far superior to all discoveries that have been made in dynamics up until now" (Figure 11, right). This essay was included in the "Dissertatio de Principio Minimae Actionis. . . ", a bilingual edition in Latin and French [START_REF] Eulero | Dissertatio de principio minimae actionis una cum examine obiectionum CL. Prof. Koenigii contra hoc principium factarum. Berolini, Ex Officina Michaelis[END_REF] together with an examination of the objections of M. Professor Koenig made against the principle (Figure 12), see also [START_REF] Voltaire | Collection of pamphlets related to controversy on Principle of Least Actions[END_REF].

Shortly thereafter in April 1753 sixteen polemics and letters in dispute were put together in a publication initiated by Voltaire [START_REF] Voltaire | Collection of pamphlets related to controversy on Principle of Least Actions[END_REF]; it was entitled "Maupertuisiana" (Figure 13) and printed in a fictitious location Hambourg; among them are the aforementioned "Reply" and "Letter". The figure on the title page shows Don Quixote (Maupertuis) fighting the windmills with a broken lance shouting "Tremblez" (tremble); behind him Sancha Panza (Euler) is riding a donkey. On the right side a satyr is depicted saying "sic itur ad astra" (thus one reaches the stars), meaning "thus one is blamed forever". Above the picture a line of Virgil's "Aeneis" is added: "Discite Justitiam, moniti" (be warned and learn justice).

The "Traité de Paix" was also included in the Maupertuisiana; in the title page of that volume the first part of a quote from Horace's Satires was added: "ridiculum acri forties ac melius (plerumque secat res)" or "ridicule often settles matters of importance better (and with more effect than severity)".

With the publication of the Maupertuisiana the quarrel came more or less to an end. Maupertuis was ill, spent a year in Paris and St. Malo, returned to Berlin which he left in 1756. After another year in St. Malo he was accommodated by Johann II Bernoulli in Basel where he died in 1759. S. König being the moral victor in this battle died from a stroke in Holland two years earlier in 1757. Voltaire had still some active and eventful years; he left Berlin in 1753 for Paris, but was banned by Louis X. He then lived in Geneva and in Ferney across the French border. In 1778 he returned for the first time in 20 years to Paris where he died after three months. Euler as the fourth figure in this dispute stayed a couple of years in Berlin, but returned to St. Petersburg in 1766. He was extremely active and published several fundamental treatises. He passed away in 1783.

On Leibniz' Alleged Letter to Hermann

The claim for forgery was taken up again 140 years later when C.I. Gerhardt, editor of Leibniz' mathematical oevre at the Prussian Academy, reinvestigated the case. The Physicist H. von Helmholtz, who gave a speech on the history of the Principle of Least Action in the Academy in 1887 [START_REF] Helmholtz | Rede über die Entstehungsgeschichte des Princips der kleinsten Action[END_REF], refers to Gerhardt's remark, that the letter does not fit to the correspondence with Hermann. Leibniz may have postponed the publication because he planned a later application. In a detailed report [START_REF] Gerhardt | Uber die vier Briefe von Leibniz, die Samuel Konig in dem Appel au public, Leide MDCCLIII, veroffentlicht hat[END_REF] in 1898 (see minutes for a meeting of the Academy) Gerhardt claims that three of the four letters presented by König are genuine and concludes at the end of his paper: "From the above explanations it should undoubtedly follow that the letter fragment, which König published, is not made up, and the entire letter has not been foisted. The letter is written by Leibniz. Also with a probability close to certainty it has been proven that it had been directed to Varignon (Pierre de Varignon (1654-1722), French mathematicien and physicist). The coincidental letters of the correspondence with Varignon in the Royal Library of Hannover have disappeared as the letters of Leibniz in the respective correspondence with Hermann did". The assumption for the authenticity was underlined by another argument in 1913, when W. Kabitz [START_REF] Kabitz | Uber eine in Gotha aufgefundene Abschrift des von S. Konig in seinem Streite mit Maupertuis und der Akademie veroffentlichten, seinerzeit fur unecht erklarten Leibnizbriefes[END_REF] found a further copy of the letter in Gotha among a collection of other genuine letters of Leibniz. This was accepted as a "proof" for the originality of the letter at that time.

Recently Breger has re-examined the case in a remarkably thorough discussion [START_REF] Breger | Über den von Samuel König veröffentlichten Brief zum Prinzip der kleinsten Wirkung[END_REF] in which he elaborated on eight arguments against the authenticity of the letter. Among those is the argument that Leibniz never applied the notion "limites" in the sense as it is used in the letter. Leibniz also often referred to minimum and maximum properties and also mentioned the term action, however combining both these concepts would be a singular event; for further arguments see [START_REF] Breger | Über den von Samuel König veröffentlichten Brief zum Prinzip der kleinsten Wirkung[END_REF], a paper worth reading. Breger recomments that now it is time to inversely analyse which arguments are speaking in favor of the genuineness. He also mentions that in case of forgery the rather strange reaction of Euler is easier to comprehend; he understood the letter as an action against his own person. A further question namely with respect to the identity of the falsifier still remains open. Opposite to the engineering community the notion of Principle of Least (Stationary) Action is still present in modern physics. Entire chapters of treatises are devoted to the principle, e.g. [START_REF] Pars | A Treatise on Analytical Dynamics[END_REF][START_REF] Berdichevsky | Variational Principles of Continuum Mechanics. Part I Fundamentals[END_REF]; see also the famous Feynman Lectures on Physics [START_REF] Feynman | [END_REF]. Occasionally they are looked upon as generalization of Hamilton's principle rather than as its predecessor. The Principle has been applied for example in the theory of relativity or in quantum mechanics, see [START_REF] Taylor | Principles of Least Action[END_REF]. 3 Principle of Least Constraint

Method of Least Squares

Carl Friedrich Gauss (1777-1855) (Figure 14) applied his method of least squares in 1801 when he determined the elliptical orbit of the astroid Ceres; however he had developed the basis of this method already in 1795 when he was 18 years old. It was not until 1809 that he published the method in the second volume of his book on the Theory of Celestial Bodies [START_REF] Gauss | Theoria Motus Corporum Coelestium in sectionibus conicus solum ambientium[END_REF] (Figure 15). He said "Our principle which we have made use of since the year 1795 has lately been published by Legendre in the work Nouvelles methodes. . . ". It was a further priority argument in the history of mechanics; however Gauss could present his correspondence with colleagues concerning its use much earlier.

Principle of Least Constraint

Twenty years later in 1829 Gauss wrote an essay on " Über ein neues allgemeines Grundgesetz der Mechanik" (On a new Fundamental Law of Mechanics) [START_REF] Gauss | Über ein neues allgemeines Grundgesetz der Mechanik[END_REF]. The paper was published in the 4th issue of Journal für die reine und angewandte Mathematik (ed. A.L. Crelle), a journal still existing today. He describes his new law on four pages that became known as the Principle of Least Constraint (Figure 16), by the way with very little algebra. It says "The motion of a system of material points. . . takes place in every moment in maximum accordance with the free movement or under least constraint;. . . ". He continues "the measure of constraint, . . . , is considered as the sum of products of mass and the square of the deviation to the free motion".

Applying the usual mathematical notation the acceleration a free i of the free unconstrained motion is defined by force F i divided by mass m i ; r is the position vector. free motion

a free i = F i m i a i = ri (3) constrained motion a i = a free i ( 4
)
The acceleration of the constraint motion is called a i where kinematic conditions constrain the corresponding material point i. For example a kinematic constraint may be given by a prescribed trajectory. The measure Z for the "constraint" (Zwang) is proportional to the sum of squares of the differences between free and constrained accelerations. Each term in the sum is weighted by mass m i .

Z ∼ N i=1 a i -a free i 2 (5) or Z = N i=1 m i a i - F i m i 2 = MIN + kinematical conditions, e.g. prescribed displacements (6)
According to the principle Z is supposed to be a Minimum. Thus from all possible motions (accelerations) the actual motion leads under given conditions to the least constraint.

Figure 17 applies the principle to the motion of a pendulum, the circle being the constrained trajectory. As expected the principle yields the equation of motion. This example elucidates that there is a strong relation with d'Alembert's Principle (Jean d'Alembert, 1717-1783). We start from the Principle of Least Constraint and its variation with respect to its free parameters a

i constraint Z = N i=1 m i a i - F i m i fix 2 = MIN (7)
variation δZ = N i=1 2m i a i - F i m i δa i =δr i =δd i = 0 (8)
The variation or virtual displacements satisfy the essential boundary conditions. δZ = 0 directly yields the Principle of d'Alembert (1743) as a variational principle of the equation of motion (Figure 18). Gauss finishes his paper with the remark (Figure 19) "It is strange that the free movements, when they cannot withstand the necessary conditions, are modified in the same way as the analyzing mathematician, applying the method of least squares, balances experiences which are based on parameters depending on necessary interactions", see the comparison in Figure 20. He adds a final sentence to his paper which says: "This analogy could be further followed up, but this is currently not my intention". However he never picked up this matter again.

Application of Principle of Least Constraint

Gauss's Principle is not very well known although it is mentioned as a fundamental principle in many treatises, e. g. [START_REF] Szabo | Geschichte der mechanischen Prinzipien und ihrer wichtigsten Anwendungen[END_REF][START_REF] Lanczos | The Variational Principles of Mechanics[END_REF][START_REF] Evans | Statistical Mechanics of Nonequilibrium Liquids[END_REF][START_REF] Päsler | Prinzipe der Mechanik[END_REF], see also [START_REF] Linkwitz | Numerical Computation and Principles of Mechanics: C.F. Gauss, the Great Mathematician, Muses Comprehensively on Figures and Formulas[END_REF]; correspondingly it has not been applied too often. Evans and Morriss [START_REF] Evans | Statistical Mechanics of Nonequilibrium Liquids[END_REF] discuss in detail the application of the Principle for holonomic (constraints depend only on co-ordinates) and nonholonomic constraints (non-integrable constraints on velocity) and conclude "The correct application of Gauss's principle is limited to arbitrary holonomic constraints and apparently, to nonholonomic constraint functions which are homogeneous functions of the momenta". Evans et al. [START_REF] Evans | Nonequilibrium molecular dynamics via Gauss's principle of least contraint[END_REF] applied the principle in the context of statistical mechanics, see also [START_REF] Bright | New observations regarding deterministic, time reversible thermostats and Gauss's principles of least constraint[END_REF]. As another example the work of Glocker [START_REF] Glocker | The Principles of d'Alembert, Jourdain and Gauss in Nonsmooth Dynamics. Part I: Scleronomic Multibody Systems[END_REF][START_REF] Glocker | Set-Value Force Laws: Dynamics of Non-Smooth Systems[END_REF] is mentioned where accelerations in rigid multibody systems under set-valued forces are evaluated applying Gauss's principle. Recently Udwadia and co-workers extended the Principle to underdetermined systems [START_REF] Fan | Reflections on the Gauss Principle of Least Constraint[END_REF].

As a final remark the author would like to point out to the possibility of using Gauss's Principle in optimization and design of structures under dynamic excitation where the system ought to be tuned for selected accelerations.

Fig. 1

 1 Fig. 1 Evolution of Extremum Principles.

Fig. 3

 3 Fig. 3 Pierre-Louis Moreau de Maupertuis.

Fig. 6

 6 Fig. 6 Leonard Euler.

Fig. 7

 7 Fig. 7 Johann Samuel König.

Fig. 8

 8 Fig. 8 König's reference letter of Leibniz to Hermann (1708) in Nova Acta Eruditorum (1751) [5].

Fig. 9

 9 Fig. 9 Franc ¸ois-Marie Arouet, called himself Voltaire.

Fig. 10

 10 Fig.10Voltaire's Diatribe[START_REF] Voltaire | Diatribe du Docteur Akakia[END_REF] and Histoire du Docteur Akakia[START_REF] Voltaire | Histoire du Docteur Akakia et du Natif de St[END_REF] (1753).

Fig. 11

 11 Fig. 11 Euler's Sur le Principe de la Moindre Action 1751/53 [12].

Fig. 12

 12 Fig.12Bilingual edition of Euler's Dissertatio[START_REF] Eulero | Dissertatio de principio minimae actionis una cum examine obiectionum CL. Prof. Koenigii contra hoc principium factarum. Berolini, Ex Officina Michaelis[END_REF].

Fig. 13

 13 Fig. 13 Maupertuisiana (1753), List of Essays (figure from [3]).

Fig. 14

 14 Fig. 14 Carl Friedrich Gauss.

2. 4

 4 Application of Principle of Least Action As said above the Principle of Least Action has been defined in more rigorous form already by Euler in 1744 as an energy principle. It turned out to be rather a Principle of Stationary Action. Lagrange derived the equation of motion based on his newly developed calculus of variations (1760). In 1834/35 Hamilton introduced the Lagrangian function into the variational principle deriving what has been called later the Euler-Lagrange equations. Whereas Euler and Maupertuis concentrated in the Principle of Least Action on the kinetic energy (vis viva, living forces) tacitly assuming constant potential energy the difference between the kinetic and the potential energy entered the Lagrangian in Hamilton's principle. Despite this close relationship the Principle of Least Action has not found its way into the engineering community which mostly refers to the works of Lagrange, d'Alembert and Hamilton. In the 19th century the well known Principles of Minimum of Total Potential Energy (Dirichlet-Green) and Minimum of Total Complementary Energy (Menabrea-Castigliano) became basic theorems in mechanics. The Principles of Virtual Work and Complementary Work entered the scene in particular because of their generality being applicable also for problems where no potentials exist. They are nowadays together with mixed Variational Principles (Hellinger-Reissner, Fraeijs deVeubeke-Hu-Washizu, etc.) the foundations for the derivation of discretization methods.

Fig. 15 C

 15 Fig. 15 C.F. Gauss's Method of Least Square (1809), extract from [23].

Fig. 16 C

 16 Fig. 16 C.F. Gauss in 1828 and his paper on the Principle of Least Constraint (1829) [24].

Fig. 17

 17 Fig. 17Application of the Principle of Least Constraint for Pendulum.

Fig. 18

 18 Fig. 18 Principle of d'Alembert.

Fig. 19

 19 Fig. 19 Conclusions of Paper [24].

Fig. 20

 20 Fig. 20 Principle of Least Constraint and Least Square Method.

E-mail: ramm@ibb.uni-stuttgart.de Phone: +49 711 685 66124 Fax: +49 711 685 66130

The quote contains a misprint; in the first definition of the action time must be replaced by distance and speed.


Jean-Daniel Boissonnat

Jane Tournois Geometryfactory

Mariette Yvinec

Anisotropic Delaunay Meshes of Surfaces

Anisotropic simplicial meshes are triangulations with elements elongated along prescribed directions. Anisotropic meshes have been shown to be well suited for interpolation of functions or solving PDEs. They can also significantly enhance the accuracy of a surface representation. Given a surface S endowed with a metric tensor field, we propose a new approach to generate an anisotropic mesh that approximates S with elements shaped according to the metric field. The algorithm relies on the well-established concepts of restricted Delaunay triangulation and Delaunay refinement and comes with theoretical guarantees. The star of each vertex in the output mesh is Delaunay for the metric attached to this vertex. Each facet has a good aspect ratio with respect to the metric specified at any of its vertices. The algorithm is easy to implement. It can mesh various types of surfaces like implicit surfaces, polyhedra or isosurfaces in 3D images. It can handle complicated geometries and topologies, and very anisotropic metric fields.

Introduction

Anisotropic simplicial meshes are triangulations with elements elongated along prescribed directions. Anisotropic meshes have been shown to be particularly well suited for interpolation of functions and solving PDE's [START_REF] Simpson | On optimal interpolation triangle incidences[END_REF][START_REF] Richard | What is a good linear finite element? Interpolation, conditioning, anisotropy, and quality measures[END_REF][START_REF] Mirebeau | Optimal meshes for finite elements of arbitrary order[END_REF]. They can also significantly enhance the accuracy of a surface representation if the anisotropy of the mesh conforms to the curvature tensor of the surface [START_REF] Heckbert | Optimal triangulation and quadric-based surface simplification[END_REF]. In this paper, we consider the problem of sampling and meshing a surface S endowed with a metric tensor field that varies over S. The metric tensor can be chosen arbitrarily provided that it satisfies a continuity condition over the surface S. The metric tensor may, for instance, be the curvature tensor of S or the Hessian matrix of a function defined over S.

Background

Due to its practical importance, several methods have been proposed in the literature to address the problem of approximating surfaces using anisotropic surface meshes. The case of a parametric surface S can be handled by working in the 2-dimensional space U of the parameters of S, using an appropriate metric field to account for the deformation induced by the mapping U → S. The problem is then 2-dimensional [BGH + 97, BH96, LT Ü99]. Several methods have also been reported to remesh a polyhedral surface according to some specified metric tensor. Alliez et al. [ACSD + 03] first estimate the principal direction fields and trace a network of lines of curvature from which they deduce a quad-dominant anisotropic mesh. Jiao et al. [START_REF] Jiao | Anisotropic mesh adaptation for evolving triangulated surfaces[END_REF] use the quadratic metric tensor of Heckbert and Garland [START_REF] Heckbert | Optimal triangulation and quadric-based surface simplification[END_REF] and propose a set of operations to allow anisotropic adaptation of a triangular mesh to static or evolving surfaces. The case of implicit surfaces can be handled by first polygonizing the surface in a standard way and then remeshing the obtained surface using one of the anisotropic methods mentioned above. Although this approach is clearly not optimal, very few methods have been proposed to directly generate anisotropic meshes on implicit surfaces. One exception is the work of Azernikov and Fischer [START_REF] Azernikov | Anisotropic meshing of implicit surfaces[END_REF] who propose a grid-based approach to mesh implicit surfaces according to their curvature tensor. The method uses a deformed grid and produces quad-meshes. The anisotropic grid is obtained by deforming a cartesian grid while maintaining the structured topology, which is a limitation when considering complex shapes.

In order to settle more rigorous foundations and to precisely characterize anisotropic triangulations, several authors have introduced anisotropic Voronoi diagrams (AVD). An anisotropic triangulation, called anisotropic Delaunay triangulation (ADT), is then defined as the dual of the AVD. Various AVD and ADT have been proposed. Leibon and Letscher [START_REF] Leibon | Delaunay triangulations and voronoi diagrams for riemannian manifolds[END_REF] introduced Delaunay triangulations and Voronoi diagrams for Riemannian manifolds. This approach requires to compute geodesic paths and intrinsic balls, which may be quite complicated in practice. To overcome this issue, two approaches have been proposed that approximate the geodesic distance between a given site and any point by considering the metric constant and equal to the metric at the site [START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF] or at the point [START_REF] Du | Anisotropic centroidal voronoi tessellations and their applications[END_REF][START_REF] Lévy | Lp centroidal voronoi tessellation and its applications[END_REF]. Although nice results have been reported, it is recognized in [START_REF] Du | Anisotropic centroidal voronoi tessellations and their applications[END_REF] that "more precise theoretical analysis are needed in identifying conditions on the metric, geometry, and the generator distributions that guarantee the well-defined duality between the AVDs and the ADTs". The Labelle and Shewchuk approach [START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF] is easier to analyze and, in the 2-dimensional case, the authors have proposed a refinement algorithm that can provably produce anisotropic meshes. This approach has been extended by Cheng et al. [START_REF] Cheng | Anisotropic surface meshing[END_REF] to produce anisotropic meshes of surfaces embedded in R 3 . The algorithm is however very complicated and no implementation has been reported. Although some progress has been recently reported in [START_REF] Canas | Orphan-free anisotropic Voronoi diagrams[END_REF], it remains unclear under which conditions an anisotropic Voronoi diagram admits a dual embedded triangulation for 3 or higher dimensional domains.

A different approach has been proposed by Boissonnat, Wormser and Yvinec [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF]. They introduced a new variant of anisotropic Delaunay meshes, called locally uniform anisotropic meshes. Such meshes are not defined by duality from some anisotropic Voronoi diagram. Instead, it is defined as a simplicial mesh in which the star of each vertex is Delaunay for the metric attached to the vertex. It is shown that, when the metric field satisfies some continuity property, such a mesh can be obtained for bounded polyhedral domains of R 2 or R 3 with no sharp edges. The extension of the approach to surfaces has not been considered so far and no implementation has been reported for 3-dimensional domains.

Our Contributions

Our work uses the notion of locally uniform anisotropic mesh defined in [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF], combined with the notion of Delaunay triangulation restricted to a surface. As in the previous approaches, we assume that the anisotropy is prescribed by a metric field that associates to each point p of the surface, a symmetric positive definite 3 × 3 matrix M p describing the metric at point p. Given a set of (punctual) sites V and v ∈ V , it is easy to compute the 3D Delaunay triangulation Del v (V ) for the metric M v of v. Indeed, Del v (V ) is just the image of an Euclidean Delaunay triangulation under a stretching transformation. For each site v ∈ V , we consider the star T v of v in Del v (V ), and the surface star S v . The star T v is the subcomplex of Del v (V ) formed by the tetrahedra that are incident to v, while the surface star S v is the subcomplex formed by the triangular facets incident to v and whose dual Voronoi edges intersect S. The surface star S v is thus a collection of triangles incident to p.

In general, there are inconsistencies among the surface stars of the sites : a facet τ , appearing in the surface stars of some of its vertices, may not appear in the surface stars of all of them. As a result, the set of surface stars of V does not form a triangulated surface.

In [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF][START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF], we showed that it is possible to refine the set V of sites until there is no more inconsistencies among the set T (V ) of 3D stars. This leads to a locally uniform anisotropic 3D mesh in which the star of each vertex is Delaunay with respect to the metric of this vertex. In the present paper, we are interested in generating surface meshes. The algorithm will incrementally compute a set V of sites on the surface and will maintain the associated surface star set S(V ). The set V is refined until there is no inconsistencies among the facets of S(V ). In contrast to [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF], we do not consider inconsistencies among the 3D stars. When there is no more inconsistencies among the surface stars, the facets of these surface stars form a locally uniform anisotropic surface mesh that approximates S and conforms to the specified metric field. This algorithm is provably correct and much more efficient than the 3D anisotropic mesh generator of [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF].

In addition to conforming to the given metric field, this new method has several notable advantages over previous methods.

-The algorithm relies on the well-established concepts of restricted Delaunay triangulations and Delaunay refinement and comes with theoretical guarantees.

-It can handle complicated geometries and topologies, as well as very anisotropic metric fields.

-It can handle various types of surfaces, e.g. implicit surfaces, polyhedra, isosurfaces in 3D images.

-It is easy to implement and robust since it relies on the usual Delaunay predicates (applied to some stretched spaces).

-The facets of the output mesh have a good aspect ratio as measured with respect to the metric field.

The idea of maintaining independent stars for each vertex of a mesh has been first proposed by Shewchuk [START_REF] Shewchuk | Star splaying: an algorithm for repairing Delaunay triangulations and convex hulls[END_REF] for maintaining triangulations of moving points. This idea has also been used to build the dual of an anisotropic Voronoi diagram as defined by Labelle and Shewchuk [START_REF] Schoen | Robust, guaranteed-quality anisotropic mesh generation[END_REF]. The idea of stitching stars restricted to a surface can be found in the work on the tangential complex [START_REF] Boissonnat | Manifold reconstruction using Tangential Delaunay Complexes[END_REF]. The context is however different from ours since no anisotropy was considered and the main motivation was to reconstruct a k-manifold embedded in a d-dimensional space without constructing any global d-dimensional triangulation, which is not a critical issue in 3 dimensions (see however the discussion in Section 5.3). A major technical difference with our paper is that the stars considered in [START_REF] Boissonnat | Manifold reconstruction using Tangential Delaunay Complexes[END_REF] are all contained in a global Euclidean d-dimensional Delaunay triangulation, which is not the case here. Combining ideas from the present paper and [START_REF] Boissonnat | Manifold reconstruction using Tangential Delaunay Complexes[END_REF] would allow to define and generate anisotropic meshes of smooth submanifolds embedded in high dimensional spaces.

Preliminaries

Anisotropic Metric

An anisotropic metric in R 3 is defined by a symmetric positive definite quadratic form represented by a 3 × 3 matrix M . The distance between two points a and b, as measured using a metric M , is defined as

d M (a, b) = (a -b) t M (a -b).
This definition can be extended to M -lengths and M -areas.

In the following, we often use the same notation, M , for a metric and the associated matrix in a given basis. Given the symmetric positive definite matrix M , we denote by F M any matrix such that det(F M ) > 0 and

F t M F M = M . Note that d M (a, b) = (a -b) t F t M F M (a -b) = F M (a -b) (1)
where the notation . stands for the Euclidean norm.

Given some metric M , an M -sphere C M (c, r), with center c and radius r, is defined as the set of points p such that d M (c, p) = r. Likewise, an M -ball B M (c, r), is defined as the set of points p such that d M (c, p) ≤ r. Note that an M -sphere is an Euclidean ellipsoid, with its axes aligned along the eigenvectors of M .

Given a simplex τ in R 3 and a metric M , we define the M -circumsphere C M (τ ) as the smallest M -sphere that circumscribes τ . The M -circumball B M (τ ) is the M -ball bounded by C M (τ ) and the M -circumradius r M (τ ) of a simplex τ is the radius of C M (τ ).

Distortion between Metrics

The following definition is due to Labelle and Shewchuk [START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF]. Given two metrics M and N , and their square-roots F M and F N , the distortion between M and N is defined as

γ(M, N ) = max{ F M F -1 N , F N F -1 M },
where . is the matrix norm operator associated with the Euclidean norm, i.e. for a 3 × 3 square matrix A, A = sup x∈R d

Ax

x . Observe that γ ≥ 1 and γ = 1 iff M = N . A fundamental property of the distortion γ(M, N ) is to relate distances d M and d N . Specifically, for any x, y, we have

1/γ(M, N ) d M (x, y) ≤ d N (x, y) ≤ γ(M, N ) d M (x, y).
In the rest of the paper, we consider a compact closed surface S of R 3 endowed with a continuous metric field {M p , p ∈ S} defined over S. For convenience, we will simply write X p instead of X Mp when a quantity X is measured using metric M p . Hence, we will write for instance F p for F Mp and d p (a, b) for d Mp (a, b).

Given two points p and q of S, the relative distortion between p and q is defined as γ(p, q) = γ(M p , M q ). Note that Γ = sup x,y∈S γ(x, y) is finite since S is compact.

Let γ 0 be some positive constant called the distortion bound. The distortion radius at p, dr(p, γ 0 ) is the upper bound on distances such that for all q and r in S, max(d p (p, q), d p (p, r)) ≤ ⇒ γ(q, r) ≤ γ 0 .

It can be shown that the distortion radius satisfies the following continuity property [BWY11] : 1 γ(p, q) (dr(p, γ 0 ) -d p (p, q)) ≤ dr(q, γ 0 ) dr(q, γ 0 ) ≤ γ(p, q) (dr(p, γ 0 ) + d p (p, q)) (2)

Sizing Field

To control the distorsion of the metric on mesh elements, we use a sizing field that is pointwise smaller than the distorsion radius. In addition to controlling the distortion of the metric, it is important to control the topology of the output mesh. To do so, we adapt the notion of local feature size to the anisotropic setting. The (anisotropic) local feature size at a point p ∈ S, denoted by lfs(p), is defined as the M p -distance from p to the M p -medial axis of S. It can be proved that the anisotropic local feature size enjoys a continuity property analoguous to (2).

To control both distortion and topology, we use a sizing field sf(p, γ 0 ) defined as sf(p, γ 0 ) = min(lfs(p), dr(p, γ 0 )). Plainly, the sizing field sf(p, γ 0 ) satisfies a continuity property analog to (2). In fact, the user may include his own sizing requirement by using as sizing field any function provided that it is pointwise smaller than sf(p, γ 0 ) and satisfies a continuity property analog to (2).

Stars and Refinement

Anisotropic Star Set

Given a metric M and a set of punctual sites V , the Delaunay triangulation Del M (V ) for metric M is the triangulation of V such that the M -circumball of each tetrahedron is empty. By empty, we mean that the interior of the M -circumball contains no points of V . By (1), the Delaunay triangulation Del M (V ) of a finite set of points V for metric M is simply obtained by computing the Euclidean Delaunay triangulation of the stretched image F M (V ), and then stretching back the result using F -1 M . For each site v in V , we consider the Delaunay triangulation Del v (V ) of V for metric M v . We define the star T v of site v as the subcomplex of Del v (V ) formed by the tetrahedra that are incident to v. The collection of all the stars T v , v ∈ V , is called the star set of V . We denote it by T (V ).

We now adapt to the anisotropic setting the definition of a restricted Delaunay triangulation. The surface star of v, denoted by S v , is the subcomplex of T (v) formed by the facets of T (v) that are incident to v and whose dual Voronoi edge intersects S. The collection of the surface stars S v , for v ∈ V , is called the surface star set of V . We denote it by S(V ).

Let N be a metric and τ be a triangle with its three vertices p 0 , p 1 and p 2 on S. We define the N -axis of τ as the set of centers of all the N -circumballs of τ . If p 0 , p 1 and p 2 are sufficiently close to each other, and if τ is well-shaped with respect to metric N , the N -axis of τ intersects S. We call surface N -circumball of τ , and note

B N (τ ), the smallest N -circumball of τ centered on S. A facet τ ∈ T v incident to v belongs to S v iff its surface M v -circumball is empty.
The facet τ with vertices p 0 , p 1 and p 2 on S has three surface circumballs B i (τ ) associated respectively to the metrics M i of each of its vertices. We define the distortion γ(τ ) of τ as the maximal distortion between any pair of points of S which are both inside the same surface circumball B i (τ ).

Inconsistencies

Two surface stars S v and S w are said to be inconsistent if some facet incident to v and w appears in only one of the two stars S v and S w . Any facet that appears in the stars of some of its vertices, but not in the stars of all of them, is also said to be inconsistent (see Figure 1). Our algorithm incrementally inserts new sites in V to resolve inconsistencies among the surface stars.

Quasi-Cosphericity

Let γ 0 > 1 be a bound on the distortion and M be a metric. Following [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF], we now introduce the notion of (γ 0 , M )-cosphericity and show its link to inconsistent simplices.

Let τ be a facet of some surface star and p ∈ V \ τ . We call U = (τ, p) a (γ 0 , M )-cospherical configuration iff there exists two metrics N and N such that

1. γ(M, N ) ≤ γ 0 , γ(M, N ) ≤ γ 0 and γ(N, N ) ≤ γ 0 ;
2. p is included in the interior of the surface N -circumball B N (τ ) and not in the interior of the surface N -circumball B N (τ ).

If M is clear from the context, we simply say that U is a γ 0 -cospherical configuration and if both M and γ 0 are understood, we say that U is a quasi-cospherical configuration. See Figure 2 for an illustration where S is the plane.

The following lemma ensures that, when all simplices in the stars have small distorsion, inconsistencies only result from the occurence of quasi-cospherical configurations.

Lemma 3.1 (Inconsistencies and quasi-cosphericities) Let τ be an inconsistent facet that appears in star S v but not in star S w where v and w are two vertices of τ . If γ(τ ) < γ 0 , then there exists a site q ∈ S w such that the configuration (τ, q) is (γ

0 , M v )-cospherical. a c b d N' (τ) N (τ) τ Figure 2: The configuration (τ, d) is quasi-cospherical because point d is outside the N -circumball B N (τ ) but inside the N -circumball B N (τ ).
Proof Take N = M v and N = M w . Because the distortion of τ is less than γ 0 , we have γ(v, w) ≤ γ 0 , and therefore the distorsions γ(M v , N ), γ(M v , N ) and γ(N, N ) are all less than γ 0 . Assume first that all the vertices of τ appear in S w . Then if τ ∈ S w , the surface M w -circumball B w (τ ) of τ contains some sites of S w in its interior. Consider now the case where some vertex v of τ does not appear in S w . This means that v w is not an edge of Del w (V ). Hence, any surface M w -ball circumscribing edge v w, and in particular B w (τ ), contains some sites of S w . Thus, in both cases, there is some site q in S w such that q is included in the interior of B w (τ ) and not in the interior of B v (τ ). It follows from Definition 3.3 that (τ, q) is a (γ 0 , M v )-cospherical configuration.

Meshing Algorithm

Algorithm Outline

Let S be a given surface we want to mesh. The algorithm constructs the set of sites V by inserting new sites on S in a greedy way. While there remain inconsistent or bad facets in the surface star set, the algorithm selects one such facet τ and kills τ by inserting a new site on S called the refinement point of τ . Inserting a new point is performed by the following Insert procedure that maintains the star sets T (V ) and S(V ). We say that a tetrahedron σ of a star T w conflicts with a point p if p belongs to the M w -circumball of p.

Insert(p)

1. insert p in each star T w that contains a simplex in conflict with p and update the surface star S w ;

2. create the new stars T p and S p .

The refinement algorithm consists in applying the following rules with a priority order : Rule (i) is applied only if no Rule (j), j < i, can be applied. The algorithm ends when no rule applies anymore. The algorithm relies on procedure Insert to insert a new point in the data structures, and on procedure Pick valid to select the location of the new site (see the next section).

If τ is a facet in some star S v , we write B v (τ ) for the surface M v -circumball of τ , and c v (τ ) and r v (τ ) for the center and the radius of this ball. The quality of facet τ with respect to metric M v is characterized by the M v -radius-edge ratio, ρ v (τ ), which measures the ratio between the M v -circumradius of τ and the M v -length of the edge of τ of shortest M -length.

The refinement rules depend on the constants α, ρ 0 and γ o .

1. Sizing field -Distortion :

If ∃v ∈ V and τ ∈ S v such that r v (τ ) ≥ α sf(c v (τ ), γ 0 ), then Insert(c v (τ ));
2. Shape -Radius-edge ratio :

If ∃v ∈ V and τ ∈ S v such that ρ v (τ ) > ρ 0 , then Insert (Pick valid(τ, M v ));
3. Inconsistency :

If ∃v ∈ V, τ ∈ S v such that τ is inconsistent, then Insert (Pick valid(τ, M v ));
Assume that the algorithm terminates. Then, no rule applies anymore. In particular rule (3) does not apply which means that the star sets S(V ) includes no inconsistent simplex. All surface stars are consistent and can be stitched together to form a triangulated surface T with the property that the star of any vertex v in T is Delaunay for metric M v . Hence, T is a locally uniform anisotropic mesh of surface S. The next step is therefore to establish that the algorithm terminates (Section 4.3).

Picking Region and Hitting Sets

The simplest idea to kill a bad triangle τ in some surface star S v is to insert a new site at the center of the surface M v -circumball of τ . However, this simple strategy may lead to endless occurrences of quasi-cosphericities and inconsistencies. We present in this section an alternative strategy, inspired by the strategy introduced by Li and Teng [START_REF] Li | Generating well-shaped delaunay meshed in 3d[END_REF][START_REF] Li | Generating well-shaped d-dimensional delaunay meshes[END_REF] to avoid slivers in isotropic meshes and also used in [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF]. The basic idea is to relax the choice of the refinement point of a bad facet. Instead of using the center of the surface circumball, the refinement point is picked from a small region around this center, and carefully chosen so as to avoid the formation of new quasi-cosphericities.

Let δ < 1 be a constant called the picking ratio. The M v -picking region of a facet τ in star S v , denoted byP v (τ ), is the intersection of the M v -ball B v (c v (τ ), δr v (τ )) with the surface S.

In fact, it is not possible, when choosing a refinement point in the picking region P v (τ ) of a triangle τ of S v to avoid the formation of any new quasi-cospherical configurations. The Pick valid procedure will only avoid the creation of small quasi-cospherical configurations where the meaning of small depends on an additional parameter β ≥ 1 and is made precise below.

Assume we are choosing a new site p in the M v -picking region of some triangle τ in surface star S v . A subset σ of three sites of the current set V is said to hit p if the configuration U = (σ, p) is both :

1. (γ 0 , M )-cospherical for some metric M such that γ(M p , M ) ≤ γ 0 , 2. small, meaning that the smallest M -circumradius of all the triangles with vertices in U is less than βr v (τ ). A subset σ of three sites of V that hits some point p in the picking region P v (τ ) is called a hitting set of P v (τ ). We associate to each hitting set σ a so-called forbidden region that consists of the points of P v (τ ) hit by σ. A point of the picking region P v (τ ) is called a valid refinement point when it does not belong to any forbidden region. See Figure 3.

Note that the definition of a valid refinement point depends on the constants δ and β that respectively defines the size of the picking regions and bounds the minimum size of new quasicospherical configurations. It also depends on the constant γ 0 that defines quasi-cosphericity. The next lemma is fundamental for our algorithm.

Lemma 4.1 (Picking lemma) For any positive values of the parameters α ≤ 1, β ≥ 1, δ < 1 and ρ 0 , it is possible to choose a distortion parameter γ 0 close enough to 1, so that the valid refinement points form a subregion of positive area inside any picking region considered by the algorithm.

The proof, only sketched here, is given in [START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF]. It uses a volume argument based on the two following facts. 1. The number of hitting sets of a picking region is bounded. 2. For given values of parameters α, β, δ and ρ 0 , the volume of the forbidden region associated to a hitting set tends to zero when γ 0 tends to 1.

To find a valid refinement point in the M v -picking region P v (τ ) of some bad k-simplex τ , the insertion algorithm calls the following Pick valid procedure that randomly selects a point in the picking region P v (τ ) until a valid refinement point is found. This procedure depends on the constants γ o , δ and β, to be fixed later in Section 4.3.

Pick valid(τ, M v ) :

1. Pick randomly a point p in the picking region P v (τ ) 2. If there exists a subset of three sites in V that hits p, then discard p and go back to 1.

3. Return p.

Lemma 4.1 ensures that the picking process will succeed with a positive probability at each trial. Observe also that Pick valid takes care of small configurations only. Hence, checking if a point of some picking region is valid requires only to consider sites that are close to the picking region. Furthermore, the number of triplets we need to consider as possible hitting sets is bounded by a constant (see the proof of Lemma 4.1 in [START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF]).

Termination of the Algorithm and Quality of the Mesh

We now prove that, under some conditions on the parameters ρ 0 , γ 0 , δ and β, the number of inserted vertices is bounded, which implies that the algorithm terminates. Let us define the intersite distance between two sites p and q as d(p, q) = min(d p (p, q), d q (p, q)), and the insertion radius r(p) of a site p as the minimum intersite distance between p and the current set of sites, i.e. the set of sites that have been inserted before p. Recall that Γ = sup x,y∈S γ(x, y) is finite since S is compact. The proofs of the following lemmas, available in [START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF], are only sketched here.

Lemma 4.2 (Insertion radius lemma) If the two conditions

(1 -δ)ρ 0 Γγ 2 0 ≥ 2 and (1 -δ)β Γγ 2 0 (1 + δ) ≥ 2
hold, then there exists a constant C such that the insertion radius r(p) of any vertex p of the mesh is at least C sf(p).

The proof uses the fact that the new point is inserted close to the center of the surface circumball of a facet in a surface star, and therefore relatively far from the current set of sites. The precise proof is a case analysis that studies the three rules of the algorithm. From this lemma, we can bound from below the intersite distance between any two points p and q in the final mesh. We consider whether p or q has been inserted last. Lemma 4.2 and the continuity property (2) of sf then allow to prove the following lemma.

Lemma 4.3 (Smallest intersite distance) If the conditions of Lemma 4.2 are satisfied, the smallest intersite distance between a site p and all the other vertices of the final mesh is at least :

d(p, q) ≥ C (1 + C)Γ sf(p).
n ≤ γ 2 0 π 1 + 2Γ(1 + C) C 2 S d M p sf(p) 2 .
When the algorithm terminates, the surface stars are consistent and can be stitched so as to make a triangulated surface mesh Ŝ. Each triangular facet of Ŝ is Delaunay for the metric of any of its three vertices, conforms to the sizing field sf (Rule 1) and has a radius-edge ratio smaller than ρ 0 (Rule 2) when measured using the metric of any of its vertices. The radius-edge ratio being a fair quality measure for triangles [START_REF] Richard | What is a good linear finite element? Interpolation, conditioning, anisotropy, and quality measures[END_REF], the last property ensures that the facets of the output triangulated surface are well-shaped with respect to the local anisotropy.

In addition, we can ensure that the topology of Ŝ is the same as the topology of S. This can be done by choosing a sizing field that is sufficiently small with respect to the local feature size of the surface. This condition, together with the bound on the radius-edge ratio of the facets, ensures that Ŝ is homeomorphic (and even ambient isotopic) to S. See for example [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF][START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF].

Implementation

The algorithm can handle different kinds of surfaces provided that we are given an oracle able to :

(1) compute a few initial points on the surface, (2) detect and compute the intersections between a given line segment and the surface.

A small set of initial points is enough. We create the initial stars and put the bad facets of the surface stars in a global queue containing all the facets to be refined. The refinement process can then start. Our implementation uses the cgal library [CGA].

Figure 4: A unit sphere endowed with a metric field derived from the sine-shaped shock function E(x, y, z) = tanh( 1 λ (2x -sin(5y))) + x 3 + xy 2 . λ = 0.6, ρ 0 = 3.0, r 0 = 0.1, γ 0 = 1.4, β = 2.5 and δ = 0.3. The mesh has 18323 vertices and 36642 facets.

The algorithm parameters

Our anisotropic surface mesh generation algorithm depends on five parameters : α, β, δ, ρ 0 and γ 0 . A priori, the parameters values should be chosen to satisfy the conditions of Lemma 4.2 and γ 0 should be set close enough to 1 to ensure the existence of valid refinement points in any picking region (Lemma 4.1). Note that, if the conditions of Lemma 4.2 are satisfied for a given choice of α, β, δ, ρ 0 and γ 0 , they are also satisfied for the same values of α, β, δ, ρ 0 and a lower γ 0 .

Conditions in Lemmas 4.1 and 4.2 are sufficient to ensure the termination of the meshing algorithm. However they are not necessary. These conditions arise from our proof of termination and our will to keep the proof simple makes them over-restrictive. In particular, the conditions in Lemma 2 depend on a global distortion bound Γ. This makes the proof easier but a more thorough proof would involve local distortion bounds. In practice, in our experiments we do not try to evaluate a global distortion bound, nor to respect the conditions of Lemmas 4.1 and 4.2. In fact, the meshing algorithm runs smoothly and yields good meshes for a quite large range of parameter values.

Parameter α should be chosen ≤ 1. This parameter is needed to prove the picking lemma. In practice, this parameter is not critical and has been set to 1 in all our experiments. Parameters β and δ, which define the picking region and the validity of a refinement point, have little influence on the quality of the final mesh. Those parameters have to be set in such a way that procedure Pick valid. finds a valid refinement point in a reasonable amount of time. By definition, we must have δ < 1, and the second condition in Lemma 4.2 suggests to choose β > 2. In all experiments reported later in this paper, we use β = 2.5 and δ = 0.3. Further experiments showed that varying the value of β in the range [2, 5] and the value of δ in the range [0.1, 0.5] has basically no effect on the results (number of vertices, computing time, and quality of the output mesh). More information on these experiments are given in the appendix.

Parameter ρ 0 , on one hand, should be as small as possible to ensure that the simplices of the final mesh are well-shaped according to the metric field. In the isotropic case, the Delaunay refinement algorithm for smooth surface [START_REF] Chew | Guaranteed-quality mesh generation for curved surfaces[END_REF] is known to terminate for any value ρ 0 > 1. The first condition in Lemma 4.2 suggests to take ρ 0 > 2 . In all our experiments, we use ρ 0 = 3.

The size of the simplices in the final mesh mostly depends on the parameter γ 0 and on the sizing field . As defined above, the sizing field takes into account the variation of the metric (distortion radius) and the geometry of the surface (local feature size). Estimating the local feature size is a notoriously difficult problem, for which existing solutions can be used [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF][START_REF] Antani | Mesh sizing with additively weighted Voronoi diagrams[END_REF]. In this paper, the emphasis is put on the metric variation and, in our experiments, we assume the local feature size to be uniformly bounded from below by some parameter r 0 that we use to specify the desired size of the mesh simplices. To take into account the metric variation, we replace the distortion radius by a direct and simpler control on the distortion of the simplices. Specifically, Rule 1 is replaced in practice by the following rule :

1. Sizing field -Distortion :

If ∃v ∈ V and τ ∈ S v such that either r v (τ ) ≥ r 0 or γ(τ ) ≥ γ 0 , then Insert(c v (τ ));

Parameter γ 0 has to be greater than 1 and close enough to 1. The influence of the value of parameter γ 0 is complex. On one hand, choosing for γ 0 a value very close to one triggers many applications of the Sizing field -Distorsion rule and leaves only a few inconsistencies to be removed by the Inconsistency rule. Such a choice also keeps the area of the forbidden regions small, which eases picking valid refinement points. On the other hand, such a choice may lead to an overly dense mesh. Raising the value of γ 0 results in a larger number of inconsistencies to be removed by the Inconsistency rule, an increasing number of calls to procedure Pick valid, which therefore slows down the meshing process. In all our experiments, γ 0 is chosen between 1.3 and 1.8. Also, for efficiency reasons, we approximate γ(τ ) by the maximum distortion between two vertices of τ . The influence of parameter γ 0 on the resulting mesh density and algorithm efficiency is shown in Figure 7. In summary, our implementation depends on parameters β, δ, ρ 0 , γ 0 and r 0 . In all our experiments, β, δ, ρ 0 have been set to fixed values.

Computing the Metric Field

The algorithm requires to compute the metric field at any refinement point on S. It can either be prescribed by the user or approximate the curvature tensor of the surface at the refinement point. We detail below several cases to illustrate the versatility of our approach.

Three-Dimensional Scalar Fields

Three-dimensional scalar fields are common in physical or mathematical applications. If E denotes such a scalar field defined in R 3 and S a surface, we may want to mesh S so as to provide the best approximation of E on S. The triangles should then be elongated orthogonally to the gradient ∇E. We define the metric field as follows.

F M = U • diag{1/(1 + ϕ ∇E ), 1, 1} • U t where U = [ ∇E ∇E , U 1 , U 2 ],
U 1 and U 2 being two arbitrary unit vectors that form an orthogonal frame with ∇E. If the gradient is zero, the metric field is isotropic. Parameter ϕ controls how the facets are stretched with respect to the norm of the gradient.

Figure 4 presents an anisotropic Delaunay mesh of a sphere endowed with a sine-shaped shock scalar field. The triangles are elongated orthogonally to the gradient direction and the density depends both on the norm of the gradient and on the changing rate of the gradient.

f (x, y, z) = x 2 + y 2 + z 2 -ak 2 2 -b (z -k) 2 -2x 2 (z + k) 2 -2y 2 . In

Curvature Tensor on an Implicit Surface

Let S be an implicit surface defined by f (x, y, z) = 0. At each point p of S, we denote by N the normal vector -∇f / ∇f and by H the Hessian of f . The principal curvatures c max , c min are the non-zero eigenvalues of the matrix

C = P • H • P where H = H/ ∇f , P = I 3 -N • N t ,
see e.g. [START_REF] Hughes | Differential geometry of implicit surfaces in 3-space-a primer[END_REF]. Let U max , U min and N be the normalized eigenvectors of C, and let U = [U max , U min , N ]. The metric at point p is defined by M p = F t p F p , where:

F p = U • ∆ • U t ,
with ∆ = diag{e max , e min , e n },

e max = max{ , c max }, e min = max{ , c min }.
Here, is a small positive constant that prevents the eigenvalues of the metric to vanish when c max and/or c min vanish, e.g., at planar or parabolic points. e n is a global constant to be discussed in the next subsection, and ∆ = diag{e max , e min , e n } stands for a diagonal matrix whose entries are {e max , e min , e n }.

Figure 5 shows an anisotropic Delaunay mesh of the "chair" implicit surface.

Curvature Tensor on Polyhedral Surfaces

To remesh a polyhedral surface S that may be an approximation of a smooth surface S, we can estimate the curvature tensor of S from S. This is done in two steps.

1. Curvature estimation at the vertices of S We first estimate the curvature tensor at each vertex v ∈ S as follows. We compute the Euclidean k-nearest vertices of v, and then apply the jet-fitting algorithm of [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] to estimate Figure 6: Result on the polyhedral surface "fertility", a model from AIM@SHAPE repository. We took r 0 = 0.1, γ 0 = 1.5. The mesh has 12,480 vertices and 24,972 facets.

the two principal curvatures c max and c min , the corresponding unit vectors U max and U min , and the normal vector N at v. From the estimated curvatures, we compute e max and e min as above.

We choose e n as explained in subsection 5.4 and store at each vertex v a resulting estimated metric tensor M (v) = U ∆ 2 U t , with U = [U max , U min , N ] and ∆ = diag{e max , e min , e n }.

Blending metric tensors

To estimate the metric tensor at any prescribed query point on S, we smoothly blend the metric tensors computed in the previous step using the following formula :

M p = v∈S ω p (v)M (v), with ω p (v) = µ p (v) v∈S µ p (v), µ p (v) = e -p-v 2 /r 2 B
where r B is a parameter that controls the locality of the blending. By default, r B can be computed by averaging the distance to the k-nearest points of each vertex in S.

M p is a positive combination of symmetric tensors and therefore has no negative eigenvalues. This property can be proved by contradiction. Assume that M p has some negative eigenvalues. Then there exists a vector x that satisfies

x t M p x = v∈S ω p (v)x t M (v)x < 0. Since all M (v)
are symmetric tensors, x t M (v)x ≥ 0 and, since according to the definition, ω p (v) ≥ 0, we have x t M p x ≥ 0, which provides a contradiction.

Figure 6 shows an anisotropic Delaunay mesh of a polyhedral surface conforming to the estimated curvature metric field.

Computing Surface Stars

Our algorithm needs to maintain the surface star S v of each vertex v in the mesh. As described above, S v can be extracted from the star T v of v in Del v (V ). Yet, maintaining stars T v is not a fully satisfactory solution. Indeed, since the vertices of the mesh all lie on the surface, most of the tetrahedra in stars T v spread through space and their diameters do not decrease significantly when the sample density increases. As a consequence, the algorithm is not local and the computing time is not optimized. To remedy this problem, we add a few vertices close to the medial axis of the surface and maintain for each vertex v a triangulation T v that does not necessarily include the star T v but nevertheless has the same restricted star S v as T v . The vertices close to the medial axis are computed from the initial points following Amenta and Bern's poles method [START_REF] Amenta | Surface reconstruction by voronoi filtering[END_REF]. The triangulation T v is defined by the following insertion rule: a site p is inserted in the triangulation T v of a site v only if p belongs to the so-called conflict zone of S v . The conflict zone of a surface star S v , denoted by Z v , is the union of the surface M v -circumballs B v (τ ) for all the facets τ of S v . The correctness of this variant of the algorithm follows from the fact that the conflict zone of star S v can only shrink when inserting a new site p, and from the fact that the surface star S v has to be updated iff p belongs to Z v . When the sample density increases, the mesh facets become small and have small surface circumballs. At each insertion of a new site p, we use the conflict zones as filters to quickly report the vertices v whose small triangulations T v have to be updated and the vertices w that have to be inserted in the triangulation T p created for the new site.

To each conflict zone Z v , we associate its (Euclidean) axis-parallel bounding box that we simply call the bounding box of v and denote by B v . We maintain in a data structure Q the bounding boxes of the current set of vertices V and we maintain V in a kd-tree Kd(V ). When inserting a new site p, we first query the data structure Q and report all the boxes that contain p. If p belongs to some box B v , we further check whether p belongs to the conflict zone Z v and, in the affirmative, we insert p in the triangulations T v . We then build the triangulation T p of the new site. T p is initialized as the star of p in the 3D Delaunay triangulation of p and of all the vertices v whose triangulations T v have been updated by inserting p. Then, we compute a bounding box for the current conflict zone Z p , and query the kd-tree Kd(V ) to compute all the vertices in this box. Each vertex in this box is inserted in T p if it belongs to the current conflict zone Z p .

Sliver Prevention

When the sampling is dense enough, quasi-cospherical configurations are tetrahedra with wellshaped facets (because of Rule 2) that have small volumes and are almost tangent to the surface. We call such tetrahedra slivers. We can reduce the number of slivers by enlarging the eigenvalue of the metric in the direction normal to the surface (eigenvalue e n defined in the preceding subsection). Indeed, the volume of those tetrahedra in the new metric will be increased, hence reducing the number of slivers and inconsistencies. Observe that this modification of the metric does not significantly affect the shape of the facets of the mesh since those facets are almost tangent to the surface. Experimental evidence indicates that e n should be of the same order of magnitude as the global maximum of c max . Our experiments reveal that our algorithm is roughly linear behavior with respect to the number of vertices of the final mesh. As is well-known, for most point sets and insertion orders, inserting a new point in a Delaunay triangulation takes a constant time. Moreover, our use of the conflict zones (Section 5.3) ensures that a constant number of stars are visited when inserting a new vertex and that a constant number of vertices are considered when creating the star of the new vertex.

Results

Performance and Complexity

Figure 7 shows computation times for meshing a torus defined by an implicit function. All meshes are generated with r 0 = 0.1 and for different values of γ 0 . Each curve in Figure 7 consists roughly of two linear portions with different slopes. This corresponds to the fact that, in a first stage, only Rule 1 is applied. During the second stage, Rules 2 and 3 and the Pick valid function are activated, which slows down the selection of new points.

Mesh Quality

Our meshing algorithm provides surface meshes of guaranteed quality since the algorithm controls the radius-edge ratio of each facet, as measured by the metric of any facet vertex. Figure 8 shows the histogram of the squared radius-edge ratios of the facets. For each facet, three radius-edge ratios were computed, one for each of the metrics of the vertices of the facet.

The quality of the meshes produced by our algorithm can still be improved by locally optimizing the positions of the vertices. This is done by adapting the optimization method of Alliez et al. [START_REF] Alliez | Variational tetrahedral meshing[END_REF]. Each vertex v is moved to the barycenter of the circumcenters of the facets in star S v , weighted by their respective areas. In our anisotropic setting, the circumcenters are M v -circumcenters, the weights are the M v -areas of the facets of S v , and the barycenter b v is an M v -barycenter. In order to preserve the locally uniform property of the anisotropic mesh, the motion should not introduce inconsistent, badly-shaped nor over-distortioned facet. Hence, instead of directly moving v to b v , we test new positions on the line segment (v, b v ) and check that the planned motion does not trigger the application of one of the Rules (1)-(3). We first test

v = b v . If v is a valid refinement point, we move v to v , otherwise, we set v = (v + v )/2 while d v (v, v ) > h.
Each star is optimized in turn. In our implementation, the stars are sorted by decreasing values of d v (v, b v ) and maintained in a priority queue. The first element is optimized first. The algorithm terminates when no more star can be optimized. Figure 8 shows that the distribution of radiusedge ratios is improved after the optimization step. The additional computing time required for this optimization is negligible.

Accuracy-Size Trade-off

Using an anisotropic metric field derived from the curvature tensor leads to surface meshes that enhance the trade-off between the accuracy of the approximation and the size of the mesh.

To measure the accuracy of a mesh, we use the mean error which we compute by uniformly sampling the mesh and averaging the Euclidean distance between each sample point on the mesh and its projection on the surface. Table 6.3 provides a comparison between isotropic meshes and curvature adapted anisotropic meshes of tori. Each torus has a small radius r = 1, while the big radius R takes different values. In both the isotropic and the anisotropic cases, the meshes have been obtained by a refinement process where the size criterion (r 0 ) has been replaced by a bound on the quality of the approximation over a facet. Hence, both meshes are guaranteed to approximate the surface up to the same error bound. The bound on the approximation error is 0.001. Each column of the table gives the number of Figure 8: The distributions of the squared radius-edge ratios of the mesh facets, respectively before (blue curve) and after (red curve) the optimization step. The radius-edge ratio of each mesh facet is computed three times, using each time the metric attached to one of the three vertices. Model : "fertility", r 0 = 1.0, γ 0 = 1.5.

Comparing the number of vertices between anisotropic and isotropic surface meshes with the same accuracy. vertices in the isotropic mesh of a piece of the torus, the number of vertices in the anisotropic mesh of the same piece of the torus, and the ratio κ between these two numbers. Anisotropic meshes have been computed using the following parameter values ρ 0 = 3.0, γ 0 = 1.8, β = 2.5, δ = 0.3 and r 0 = 1.0. Figure 9 shows a piece of the mesh of the torus with R = 50 that was generated by our algorithm and highlights the high anisotropy ratio that can be reached.

As can be seen, the number of vertices N a of the anisotropic mesh can be much smaller than the number of vertices N i of the corresponding isotropic mesh with the same accuracy. The table also highlights the fact that the ratio between N a and N i decreases with the ratio r/R.

Conclusion

We have introduced a novel approach to generate anisotropic meshes of surfaces. Given a surface S endowed with a metric tensor that varies over S, our algorithm produces a mesh that approximates S and whose elements conform to the metric field. The resulting mesh provides a better approximation of the surface (if one takes as metric tensor the curvature tensor of S) or a better approximation of a function f defined over S (if we take as metric tensor the hessian of f ). Our algorithm is simple and comes with theoretical guarantees without relying on heuristics. As demonstrated by our experiments, it can handle complex shapes and metric fields. The algorithm is generic in the sense that it can handle various types of surfaces such as implicit surfaces or triangulated surfaces. Although we have not presented such results, it could be easily extended to mesh isosurfaces in 3D images.

As future work, we intend to consider surfaces with sharp features following the approach of [START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF]. We also intend to extend our method to evolving surfaces [START_REF] Li | Generating well-shaped delaunay meshed in 3d[END_REF], and to 3D domains bounded by surfaces and equipped with a 3D metric field, combining the results of this paper with [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF].

More on parameters β and δ.

We provide here a few additional experimental results, to support our claim that parameters β and δ have little influence on the quality of the final mesh.

Table 7 shows computing times and numbers of vertices needed to mesh a torus (r = 1, R = 10) with r 0 = 0.1, γ 0 = 1.5 and different values of β (2,3,4,5) and δ (0.1, 0.3, 0.5). As can be seen, changing the value of β has basically no effect. The situation is different with parameter δ. As expected, choosing a small δ induces small picking regions and makes the finding of a Figure 9: Anisotropic Delaunay mesh of a piece of torus with big radius 50 and small radius 1. We took ρ 0 = 3.0, γ 0 = 1.8, r 0 = 1.0, β = 2.5, δ = 0.3. The approximation error is smaller than 0.001. valid refinement point harder. Hence, the picking process fails more often, leading to more vertex insertion trials and to an increased computing time.

Figure 10 shows the quality of the meshes obtained in the above experiment. Specifically, the figure shows the distribution of the squared radius-edge ratios of the mesh facets as in Fig. 8. We have selected three histograms to highlight the fact that the overall quality of the mesh is not impacted by varying β nor δ. The first histogram corresponds to our default values (δ = 0.3, β = 2). The two others correspond to one default value and one "large" value, specifically (δ = 0.3, β = 5), and (β = 2, δ = 0.5).

Figure 1 :

 1 Figure 1: Example of inconsistent stars in 2D : stars S v and S w are inconsistent because edge [vw] belongs to S v but not to S w .

Figure 3 :

 3 Figure 3: The picking region of triangle τ . Set σ = {q, r, s} is a hitting set for the picking region P v (τ ). It defines a forbidden region, shown with red dashes, to be avoided by the refinement point p of τ .

It follows from the

  lemma that the M p -balls B p (p, l(p)), where l(p) = C 2(1+C)Γ sf(p), have disjoint interiors. A simple volume calculation then provides an upper bound on the number of vertices of the final mesh.

Lemma 4. 4 (

 4 Number of inserted sites) If the conditions of Lemma 4.2 are satisfied, the number n of sites inserted by the algorithm is at most

Figure 5 :

 5 Figure 5: The "chair" implicit surface with a curvature-related metric field. The equation of the chair is f(x, y, z) = x 2 + y 2 + z 2 -ak 2 2 -b (z -k) 2 -2x 2 (z + k) 2 -2y 2 .In this example,

  Figure 5: The "chair" implicit surface with a curvature-related metric field. The equation of the chair isf (x, y, z) = x 2 + y 2 + z 2 -ak 2 2 -b (z -k) 2 -2x 2 (z + k) 2 -2y 2 . Inthis example, a = 0.8, b = 0.4, k = 1.0, ρ 0 = 3.0, r 0 = 0.1, γ 0 = 1.3, β = 2.5 and δ = 0.3. The mesh has 24164 vertices and 48336 facets.

Figure 7 :

 7 Figure 7: Computing time to mesh an implicit torus with r 0 = 0.1 and different values of γ 0 between 1.3 and 1.8.

  and timings to mesh a torus with different values of β and δ. number of vertices β = 2

Acknowledgments

This work has been partially supported by the Agence Nationale de la Recherche (ANR) under project GIGA (Geometric Inference and Geometric Approximation) and by the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No. 339025 GUDHI (Algorithmic Foundations of Geometry Understanding in Higher Dimensions). The second author was partially supported by Chinese 973 Program (2010CB328001), Chinese 863 Program (2012AA040902) and the NSFC (61035002, 61272235).