Two Simulated Annealing Optimization Schemas for Rational Bézier Curve Fitting in the Presence ...

18 Pages·2016·8 MB·English

Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2016, Article ID 8241275, 17 pages
http://dx.doi.org/10.1155/2016/8241275
Research Article
Two Simulated Annealing Optimization Schemas for
Rational Bézier Curve Fitting in the Presence of Noise
AndrésIglesias,1,2AkemiGálvez,1andCarlosLoucera1,3
1DepartmentofAppliedMathematicsandComputationalSciences,E.T.S.I.Caminos,CanalesyPuertos,
UniversityofCantabria,AvenidadelosCastros,s/n,39005Santander,Spain
2DepartmentofInformationScience,FacultyofSciences,TohoUniversity,2-2-1Miyama,Funabashi274-8510,Japan
3StateMeteorologicalAgency(AEMET),28071Madrid,Spain
CorrespondenceshouldbeaddressedtoAndre´sIglesias;[email protected]
Received31August2015;Revised18November2015;Accepted24November2015
AcademicEditor:AnnaPandolfi
Copyright©2016Andre´sIglesiasetal.ThisisanopenaccessarticledistributedundertheCreativeCommonsAttributionLicense,
whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited.
Fittingcurvestonoisydatapointsisadifficultproblemarisinginmanyscientificandindustrialdomains.Althoughpolynomial
functionsareusuallyappliedtothistask,therearemanyshapesthatcannotbeproperlyfittedbyusingthisapproach.Inthis
paper,wetacklethisissuebyusingrationalBe´ziercurves.Thisisaverydifficultproblemthatrequirescomputingfourdifferent
setsofunknowns(dataparameters,poles,weights,andthecurvedegree)stronglyrelatedtoeachotherinahighlynonlinearway.
Thisleadstoadifficultcontinuousnonlinearoptimizationproblem.Inthispaper,weproposetwosimulatedannealingschemas
(theall-in-oneschemaandthesequentialschema)todeterminethedataparameterizationandtheweightsofthepolesofthe
fittingcurve.Theseschemasarecombinedwithleast-squaresminimizationandtheBayesianInformationCriteriontocalculate
thepolesandtheoptimaldegreeofthebestfittingBe´zierrationalcurve,respectively.Weapplyourmethodstoabenchmarkof
threecarefullychosenexamplesof2Dand3Dnoisydatapoints.Ourexperimentalresultsshowthatthismethodology(particularly,
thesequentialschema)outperformspreviouspolynomial-basedapproachesforourdatafittingproblem,eveninthepresenceof
noiseoflow-mediumintensity.
1.Introduction electronic components for different purposes. For instance,
to analyze how a new device or machine in the market is
(1)Motivation.Thispaperdealswiththeproblemoffittinga builtorhowaparticularcomponentworks.Anotherrelevant
collectionofnoisydatapointsbymeansofarationalcurve. application arises in automotive and aerospace industries,
Thisproblemarisesinseveralscientificandappliedfields.It whereprototypesarebuiltonclay,foamrubber,wood,and
isanimportantissue,forinstance,inregressionanalysisfor othermaterialstohelpdesignersandengineersexplorenew
statisticsandmachinelearningandinapproximationtheory ideas about shape or size and get a visual insight of a new
for numerical analysis. It also plays a key role in several parttobedesigned.Thisisachallengingtask,sincenowadays
industrial fields, most prominently in reverse engineering. designs are becoming more and more organic in shape,
In its most comprehensive meaning, reverse engineering making them more difficult to be replicated by computer
consists of obtaining a digital replica of an already existing from scratch. A typical approach in this regard is to obtain
physicalobjectorcomponent.Thisisatypicalprocedurein a set of measurements of the object or workpiece and then
medicalandhealthareas,wherenoninvasivetechniquessuch reconstructitasa3Dmodel.Typicalwaystomeasureinclude
as magneticresonanceimaging(MRI) orcomputertomog- scanning technologies such as 3D laser scanners, touch
raphy (CT) are commonly used to visualize inner organs scanners, coordinate measuring machines, light digitizers,
or different parts of the human body for medical check, and industrial computer tomography. Typically, the process
diagnosis, and therapy. Reverse engineering is also a com- yieldsacloudofdatapoints,whichhavetobefittedinorder
mon practice in consumer products, microchips, and other torecoverthetopologicalinformationoftheoriginalmodel.
2 MathematicalProblemsinEngineering
Whendatacomefromaccuratesources(i.e.,withmini- importanttoremarkthatconicscanbecanonicallydescribed
mal measurement errors) and smooth shapes, curve recon- asrationalfunctions.Manyothercomplicatedshapescanbe
struction can be done through well-known curve interpo- dealtwithmoreefficientlyifthe(morelimited)polynomial
lation techniques such as spline and Be´zier curve inter- scheme is extended to a rational one. It is not a trivial
polation. This approach is characterized by imposing the task,however,sincetherationalschemeincludesadditional
parametric curve to pass through all data points. However, degreesoffreedomthathavealsotobecomputed.Asaconse-
most reverse engineering applications retrieve data from quence,ourpreviousmethodisnolongervalidandmustbe
scanning,digitizing[1],orsomemechanicaldevice,suchas substantiallyimproved.Thisisactuallythegoalofthepresent
pressuremachinesintheshoemakingindustry[2],thusbeing contribution.
subjected to measurement errors. In case of high-intensity
noise, interpolation methods tend to fail since they force (2)AimsandStructureofThisPaper.Inthispaperwefocuson
the curve to pass through all points, including the noisy obtaininganaccuratefittingofagivensetofnoisydatapoints
outliers.Thisundesirablesituationcanbeovercomebydata by using a rational Be´zier curve. To this aim, we consider
approximation, where the fitting curve is only required to a modification of the standard simulated annealing (SA), a
be sufficiently close to the original data according to some popular probabilistic method for global optimization based
prescribed metrics. Owing to its ability to cope with noise on thermodynamical principles. In our approach, the basic
inherenttodatainreal-worldsettings,inthispaperwewill SAalgorithmisfurtherimprovedbyconsideringtwolarger
useacurveapproximationmethod. globaloptimizationschemas.Theyareappliedtocomputeall
Several families of functions can be used for this curve relevant parameters of our curve fitting problem. A typical
approximation task. Among them, the free-form para- shortcominginthisregardisthattheaccuracyofthemodel
metric functions (such as Be´zier and B-splines) are the increasesasthenumberofpolesincreases,eventuallyleading
most common in CAD/CAM (Computer-Aided Design/ to overfitting. To prevent this undesirable effect, we also
Manufacturing), computer graphics and animation, virtual aim to seek a suitable trade-off between data fidelity and
reality, and other fields. In particular, Be´zier curves have modelcomplexity.ToattainthisgoalweapplytheBayesian
been intensively used in automotive industry for decades. Information Criterion (BIC), a model-selection technique
Mathematically, they are given by a linear combination of widelyusedinthecontextofinformationsciences.Thebasic
basisfunctionscalledtheBernsteinpolynomialswithvector idea of BIC is to combine a penalty term with the fitness
coefficientscalledpolesorcontrolpoints.Thecurvefollows function.Thispenaltytermincreaseswiththenumberoffree
approximatelytheshapeofitscontrolpolygon(thecollection parameters ofthe model, thuspenalizinganyunnecessarily
of segments joining the poles), and, hence, it reacts to largecomplexityofthemodel.
the movement of its poles by following a push-pull effect. The structure of this paper is as follows: Section2
Thisnicefeaturewasfundamentalforthepopularizationof summarizes briefly the previous work in the field. Some
free-formcurvesforinteractivedesign.Althoughnowadays basic mathematical concepts and definitions and the curve
Be´zier curves have been mostly deprecated in that field, fitting problem from noisy data points with rational Be´zier
beingovertakenbythemorepowerfulB-splines,theyarestill curvesaredescribedinSection3.Then,thesimulatedanneal-
in use in many other areas. For instance, True Type fonts ing algorithm is described in detail in Section4. Section5
usecompositecurvescomprisedofquadraticBe´ziercurves. describesallstepsofourproposedmethodologyalongwith
Similarly, all Postscript font outlines are defined in terms thetwoschemasusedinthispaper.Theexperimentalresults
ofcubicandlinearBe´ziercurves.Otherrecentapplications oftheapplicationofourmethodtothreeillustrativeexamples
include,forinstance,robotpathplanning[3]andthedeter- are reported in Section6. The section also discusses the
minationoftheairfoilgeometryfromagiven𝐶𝑝-distribution robustnessofthemethodinthepresenceofnoise.Thepaper
[4]. concludeswiththemainconclusionsofthisworkandsome
Arecentpaperbytheauthorsaddressedtheproblemof futurelinesofresearch.
curveapproximationwithpolynomialBe´ziercurves[5].The
methodwasbasedonahybridschemecombiningapopular 2.PreviousWork
single-particle metaheuristic called simulated annealing for
globalsearchandalocalsearchoptimizerforfurtherrefine- Theproblemofdataapproximationwithfree-formcurveshas
mentoftheglobalsolution.Thishybridschemewasapplied beenaclassicalsubjectofresearchformanyyears.Thefirst
to compute all relevant parameters of the approximating research works addressing this issue were published in the
functionwithsatisfactoryresults. 60sand70sandweremostlybasedonnumericalprocedures
Although the polynomial representation is very easy to [6–8]. The classical methods of this period select the free
dealwith(and,therefore,veryconvenientformanypurposes) parametersbyinferringsomegeometricpropertiesfromdata
itisstilllimitedinsomeways.Forexample,itiswellknown in order to meet specific constraints [1, 9, 10]. A major
that the polynomial-based schemes cannot represent accu- probleminthisregardistoperformdataparameterization,
ratelysomeimportantshapessuchastheconics(e.g.,circles, thatis,todeterminesuitableparametervaluesassociatedwith
ellipses, and hyperbolas). A way to overcome this draw- thedatapoints.Typicalchoicesaretheuniformparameter-
back is to consider the polynomial curve in homogeneous ization and the arc-length parameterization. Unfortunately,
coordinates, leading to rational functions (i.e., functions classical numerical methods are not able to obtain optimal
that are quotients of two polynomials). In this sense, it is solutionsinthegeneralcase.Subsequentattemptstoaddress
MathematicalProblemsinEngineering 3
thisproblemduringthenextdecadesprovidedsolutionsfor and hyperbolas). One way to overcome this limitation is to
several particular cases but were not successful in solving use homogeneous coordinates (see [10, 30] for details). The
the general problem [9, 11]. More recent methods apply basicideaistoconsidertheprojectionofthestandardpoly-
alternative approaches such as line use error bounds [12], nomialBe´ziercurveinR𝑑+1,withnewpolesbℎ.Theresulting
𝑗
curvature-basedsquareddistanceminimization[13],ordom- curveinR𝑑iscalledarationalBe´ziercurve.Mathematically,
inantpoints[14].Althoughtheyperformwell,theyarestill
thiscurvecanbedescribedasaquotientoftwopolynomials,
limitedtoparticularcasesand/orrequirestrongconstraints
orasalinearcombinationofrationalbasisfunctions:
intermsofdifferentiabilityornoiselessdata,assumptionsthat
arenotreasonableforreal-worldinstances. 𝑛
Mathematically speaking, curve fitting is usually C(𝑡)= ∑b𝑗𝑅𝑗𝑛(𝑡) with 𝑡∈[0,1], (3)
𝑗=0
addressed as a continuous optimization problem [6, 8, 9].
However, traditionalmathematical optimization techniques
wheretherationalbasisfunctionsaredefinedby
have had little success in solving this problem. In the last
few years, the field has taken a new impetus with the 𝑤 𝐵𝑛(𝑡)
applicationofseveralpowerfulartificialintelligenceandsoft 𝑅𝑗𝑛(𝑡)= ∑𝑛 𝑗𝑤𝑗𝐵𝑛(𝑡), (4)
computingtechniques.Someofthesemethodsarebasedon 𝑘=0 𝑘 𝑘
neuralnetworks,suchasthestandardneuralnetworks[15], where𝑛isthecurvedegreeand𝑤𝑗isthelastcoordinateofthe
Bernsteinbasisfunctionnetworks[16],andKohonen’sSOM
homogeneouscontrolpointbℎ.Thissetofnewscalarparam-
(Self-Organizing Maps) nets [17]. In some other cases, the 𝑗
neural approach is combined with partial differential eters {𝑤𝑗}𝑛𝑗=0, called weights, provides us with additional
equations [18]. This scheme has been extended to the more degrees of freedom for better shape approximation. They
general functional networks in [19, 20]. More recently, alsoincreasethemodelcomplexity,however,asweintroduce
powerfulmetaheuristictechniqueshavebeenappliedtothis anewsetofparametersthathavetobecomputedaswell.
problem,includingswarmintelligence[21],artificialimmune
systems [22], genetic algorithms [23, 24], physics-inspired 3.2. The Fitting Problem. Let now {Q𝑖}𝑚𝑖=1 be a set of data
thermodynamics [5], cuckoo search [25], support vector pointsinR𝑑.Theproblemconsistsofobtainingtherational
machines[26],fireflyalgorithm[27],andhybridtechniques Be´zier curve, C(𝑡), of a certain degree 𝑛 providing the best
[28, 29]. These methods obtain remarkable results for
least-squaresfittingofthedatapoints.Thisleadstoamini-
polynomialcurves,butmorecomplicateshapes(e.g.,conics) mizationproblemoftheleast-squareserrorΘdefinedasthe
are still not properly fitted. This is a good evidence of the
weightedsumofsquaresoftheresiduals:
interest and difficulty of this data fitting problem, which is
describedindetailinthenextsection. 𝑚 ∑𝑛 𝑤 b 𝐵𝑛(𝑡) 2
𝑗=0 𝑗 𝑗 𝑗 𝑖
Θ=∑𝜇 (Q − ) , (5)
𝑖 𝑖 ∑𝑛 𝑤 𝐵𝑛(𝑡)
3.DescriptionoftheProblem 𝑖=1 𝑗=0 𝑗 𝑗 𝑖
3.1.BasicConceptsandDefinitions. Inthissectionweassume where {𝜇𝑖}𝑚𝑖=1 are scalar numbers used in situations when it
maybereasonabletoassumethatsampleddatashouldnotbe
that the reader is familiar with the concept of parametric
treatedequally.Inordertoreflectfaithfullythemostcommon
curves.Be´ziercurvesareaparticularcaseoffree-formpara-
situationinreal-worldproblems,inthispaperwewillassume
metric curves, introduced in the 60s for interactive design
thatnoinformationabouttheproblemisavailablebeyondthe
intheautomotiveindustry.Anonrational (i.e.,polynomial)
parametricBe´ziercurveofdegree𝑛inR𝑑isgivenby deqautaalplyo;intthsa.tThisi,s𝜇m𝑖 e=ans1,thfaotraalllld𝑖a.tNaoptoei,nhtsomweuvsetrb,ethtaretaoteudr
𝑛 methodisindependentonthevaluesof𝜇𝑖.Torepresentthe
B(𝑡)= ∑b𝑗𝐵𝑗𝑛(𝑡) with 𝑡∈[0,1], (1) geometrical distribution of the data we need to associate a
𝑗=0 parameter𝑡𝑖foreachinputpointQ𝑖.Therefore,ourgoalisto
where {b𝑗}𝑛𝑗=0 ⊂ R𝑑 are vector coefficients (usually called oItbitsaoinbvthioeutshtrheaetsseitnscoefepaacrhambleetnedrsin{g𝑡𝑖f}u𝑚𝑖=n1,c{t𝑤io𝑗n}𝑛𝑗i=n0(,2a)nadn{db𝑗(}4𝑛𝑗)=i0s.
polesorcontrolpoints)and𝐵𝑗𝑛(𝑡)aretheBernsteinpolynomials nonlinearin𝑡,system(5)isalsononlinear.Asaconsequence,
ofindex𝑗anddegree𝑛,givenby we have to deal with a multivariate continuous nonlinear
minimization problem. In this paper we solve this problem
𝑛
𝐵𝑛(𝑡)=( )𝑡𝑗(1−𝑡)𝑛−𝑗 byapplyingtwodifferentschemasofthesimulatedannealing
𝑗 𝑗 optimizationmethod,whichisdescribedinthenextsection.
(2)
𝑛 𝑛!
with ( )= , 0!=1. 4.SimulatedAnnealing
𝑗 𝑗!(𝑛−𝑗)!
One of the major trends in global optimization during the
Notethatinthispapervectorsaredenotedinbold.Thepoly- last few decades has been to build algorithms trying to
nomialrepresentationisnotpowerfulenoughtorepresenta mimic certain efficient optimization patterns observed in
varietyofshapes,particularlytheconics(e.g.,circles,ellipses, natural processes. As a result, a series of very powerful
4 MathematicalProblemsinEngineering
nature-inspiredoptimizationalgorithms(e.g.,particleswarm physical system. Given an initial (usually random) state x0,
optimization,geneticalgorithms,orantcolonyoptimization) the goal of the SA is to obtain the state with the minimum
havebeendevised.Veryoften,theyprovidebettersolutions energy(associatedwiththebestsolutionoftheoptimization
thanprevioustraditionalmathematicalalgorithmstoseveral problem).Inourcase,wehaveareal-valuedfunction𝑓:D⊆
hardoptimizationproblems.Althoughtheyareverydiverse, R𝑑 → Rthatwewanttominimize.Thealgorithmperforms
allofthemsharetwocommonfeatures:tobeinspiredbyreal- aniterativeprocess;ateachiterationstep,anewstatex is
new
worldobservationandtosearchforsolutionsinastochastic generatedfromthecurrentone,x ,throughaneighborhood
old
way. Most of them are also derivative-free, meaning that function,denotedbyN : D → D,thatis,x = N(x ).
new old
they can be applied to problems where it is not possible to Let now 𝑓 ≡ 𝑓(x ), 𝑓 ≡ 𝑓(x ) be their associated
old old new new
compute the derivatives of the objective function (or they energies,respectively.Thealgorithmprobabilisticallydecides
areveryexpensivecomputationally).Inthispaperweapply betweenmovingthesystemtothenewstatex orstayingin
new
simulatedannealing,oneofthemostpopularnature-inspired currentstatex .Thisnewstateischosenwithaprobability
old
optimizationalgorithms. functionP:D×D → [0,1],calledtheacceptancefunction,
whichdependsontwofactorsasfollows:
4.1. Background. The simulated annealing (SA) is a nature-
inspired metaheuristic optimization algorithm introduced (i)ononehand,onthedifferencebetweentheirenergy
by Kirkpatrick et al. in the context of combinatorial opti- values,Δ=𝑓old−𝑓new,
mization [31]. The algorithm is inspired in the annealing (ii)on the other hand, on a global parameter called
processinmetallurgy,wheremetalsareheatedatveryhigh temperature,denotedby𝑇,whichvariesaccordingto
temperatures and then slowly cooled down to reach a state astrictlydecreasingfunctionT : R+ → R+ called
of lower energy. During the process, atoms tend to move thecoolingfunction.
to configurations that minimize the system energy even if
duringsuchmigrationcertainconfigurationsrisethesystem The probability function is not arbitrary, but must comply
overallenergy(whenitstabilizesforafixedtemperature,we with certain requirements. One of them is that P > 0 if
callitthermalequilibrium).Suchmovesaremoreprominent Δ < 0,meaningthatthesystemmaymovetothenewstate
at the beginning of the process than at the end, when evenifitisworse(hasahigherenergy)thanthecurrentone.
the particles loose thermal mobility in order to polish the The primary goal of this condition is to prevent stagnation
systeminnerstructuretofinallyproduceabettermetal.Asa (whenthesystemgetstrappedintheneighborhoodoflocal
result,themetalsbecomestrongerandwithbetterproperties, optima,leadingtoprematureconvergence).Underthesame
speciallyiftheprocessisconductedseveralconsecutivetimes conditionΔ < 0,wealsoimposethatP → 0as𝑇 → 0,
(aprocesscalledreannealing). while P → 𝜂 > 0 otherwise. Basically, these conditions
TheoriginalSAalgorithmisanadvancedinterpretation state that, for sufficiently small values of 𝑇, the system
of the Metropolis-Hastings algorithm [32] to generate sam- will increasingly promote “downhill” changes (i.e., changes
ple states of a thermodynamic system, showing the deep leading to lower energy values) and avoid “uphill” changes.
connectionsbetweenstatisticalmechanicsandcombinatorial Inotherwords,thelowerthetemperatureis,theeasieritis
optimization.Given an initial(usually random)state in the torejectaworsesolution.Infact,intheparticularcase𝑇 =
solutiondomain,thealgorithmiterativelyperturbsit.When- 0, the procedure will only allow downhill moves, meaning
everabettersolutionisfound,thechangeisalwaysaccepted; that the algorithm reduces to a greedy search algorithm.
otherwise, it is accepted only with a certain probability. Another desirable feature is to promote small uphill moves
Thisprobabilityishigheratthebeginning(mimickingwhat over large ones. This effect can be obtained by modulating
happensinthethermodynamicprocessathightemperatures) the probability as a function of the parameter Δ so that
thanattheend.Inotherwords,thisideaofslowcoolingis P decreases as Δ increases. Under these conditions, the
translatedasaslowdecreaseoftheprobabilityofaccepting temperature𝑇becomesacriticalparameterindescribingthe
suchworsesolutions.Soessentiallythesystemevolvesfrom evolutionofthesystem,asitsvaluedeterminesthesensitivity
a free exploration of the search space at initial stages to a ofthesystemtoenergyvariations.
stochastichill-climbingatlatterstages. Inthispaper,weconsidertheclassicalacceptancefunc-
Since its publication the algorithm has received a lot of tion derived from the Metropolis-Hastings sampling algo-
attention from the scientific community, with many real- rithm,firstintroducedbyMetropolisetal.in[32]asaMonte-
worldapplicationsinthemostdiversefields,rangingfromthe Carlo method to simulate the creation of new states in a
classicalNP-hardcombinatorialtravellingsalesmanproblem thermodynamicsystemandgivenby
[33] to the minimization of highly multimodal real-valued
functions [34]; see [35] for an in-depth review of several {{1 if Δ≤0
simulatedannealingapplications. P(Δ)={{ Δ (6)
exp(− ) otherwise,
{ 𝑇
4.2.TheAlgorithm. Thesimulatedannealingalgorithmwas
originally designed to compute a good approximation of where 𝑇 represents the system temperature at the iteration
theglobaloptimumofafitnessfunction(usuallycalledthe where x has been generated. With this function, a bet-
new
system energy) within a problem domain D, assumed to be ter candidate solution is always accepted, but even worse
continuousinthispaper.Eachpointx∈Disastateofsome solutions have a chance to be accepted. It is clear that
MathematicalProblemsinEngineering 5
Require:(InitialParameters)
(1) InitialTemperature𝑇0
(2) ScheduleCriteriaS𝑐
(3) StoppingCriteriaSC
(4) NeighbourhoodfunctionN:D → D
(5) SystemEnergy𝑓:D⊆R𝑑 → R
(6) CoolingScheduleT:R+ → R+(strictlydecreasingfunction)
(7) 𝑇←𝑇0 ⊳Temperatureinitialization
(8) Randomlyselectstartstatex0∈D
(9) xold←x0
(10)repeat ⊳Externalloop:stoppingcriteria
(11) repeat ⊳Internalloop:thermalequilibrium
(12) x ←N(x ) ⊳Newstategenerated
new old
(13) 𝑓 ←𝑓(x ) ⊳Energyevaluationofcurrentstate
old old
(14) 𝑓 ←𝑓(x ) ⊳Energyevaluationofnewstate
new new
(15) Δ←𝑓 −𝑓
old new
(16) if Δ<0then ⊳Acceptancecriterion
(17) x ←x
old new
(18) else
(19) Randomlycompute𝑢∈Rand(0,1) ⊳Metropolisprocedure
(20) if 𝑢≤exp(−Δ/𝑇)then ⊳Equation(6)
(21) x ←x
old new
(22) endif
(23) endif
(24) untilS𝑐==true
(25) 𝑇←T(𝑇) ⊳Temperatureupdate
(26)untilSC==true
(27)return x ⊳Bestfinalsolution
new
Algorithm1:Simulatedannealing.
(6) meets all required conditions for a proper acceptance dependingontheproblemathand.Inthissense,thesecriteria
functionindicatedabove.Italsoprovidesanadequatetrade- can be either variables or functions or even rules. In our
offbetweenexplorationandexploitation:atinitialstages(i.e., implementationitisassumedthattheyareBooleanfunctions
highertemperatures)thealgorithmexploresthesearchspace that return true or false depending on whether or not the
while at latter stages it resembles a hill-climbing algorithm, given conditions for each particular problem are met. The
withthedifferencethatnowthereisalwaysthepossibilityto specific conditions used for our problem are explained in
acceptaworsestate. Section5.
Thecorrespondingpseudocodeofthesimulatedanneal-
ingalgorithmusedinthispaperisshowninAlgorithm1.A 5.TheProposedMethod
typicalexecutionbeginswitharandomlychosenstate(x0,𝑓0)
and an initial (very high) temperature 𝑇0. The algorithm 5.1.OverviewoftheMethod. Asdiscussedabove,ourproblem
generates new states, according to N, at each iteration consists of reconstructing the underlying shape of a given
and probabilistically decides whether or not the new state set of noisy data points by using a rational Be´zier curve.
is accepted according to the probability P. According to This implies solving a nonlinear least-squares minimiza-
the pseudocode, the algorithm can be summarized as the tionproblemwhilesimultaneouslyminimizingtherequired
interaction between two cycles: the external one controls number of free parameters. Solving this problem requires
the temperature updating and the inner one performs the computingfourdifferentsetsofunknowns:dataparameters,
Metropolis procedure for a given outer iteration. The tem- poles, weights (represented in this section by vectors p, w,
perature 𝑇 is updated through the cooling function T only and b, resp.), and the curve degree, 𝑛. Our approach to
when the thermal equilibrium is reached. This workflow is tackle this issue is a hybrid strategy combining classical
repeateduntila stopping criterion ismet.Classicalstopping methods (least-squares minimization), modern stochastic
criteriaareasfollows:thesystemreachesastategoodenough methods (simulated annealing), and information science
forthespecificapplication,themethodreachesaprescribed metrics(BayesianInformationCriterion(BIC)).
numberofiterations,orthesolutiondoesnotimproveaftera Before explaining how our method works, let us intro-
prescribednumberofconsecutiveiterations.Notethatboth duce the following notation: from now on, we will use the
the cooling schedule criterion (line 24) and the stopping subindex (⋅)𝑤 when searching for the curve weights, (⋅)𝑝
criterion (line 26) can be defined in many different ways, when searching for the curve parameters, and (⋅)𝑤,𝑝 when
6 MathematicalProblemsinEngineering
Take next value ofnwithin its range
wR
n SA-LSQ pR BIC Θ BBeessttΘn
bR Degree
determination
wR
wR pR
n SA LSQ
pR bR
Data fitting
w0,p0 SAR(p,w) wR,pR Asllc-hinem-oene
Data parameterization and weight computation
w0,p0 SABp pB SARw wR,pB SARp wR,pR
(1) (2) (3)
(4)
Sequential schema
Figure1:Workflowoftheproposedmethoddescribedinthreelayers(fromuppertolowerlayer):thegeneralworkflow;decompositionof
theSA-LSQprocedure;thetwodifferentSAschemasintroducedinthispaper:all-in-oneschema(top)andsequentialschema(bottom).
searchingforbothparametersandweights.Thus,SA𝑝stands schema and the all-in-one schema. They are explained in
for simulated annealing applied to parameter search, 𝑓𝑤,𝑝 detailinthenextsection.
meanstheobjectivefunctionwithdomainofdefinitiony =
(w,p)∈W×D⊂R𝑛+1×[0,1]𝑚,andsoon.Thesuperindex 5.2. Data Parameterization and Weight Computation. In
(⋅)𝐵denotesthesearchforanonrationalBe´ziercurveand(⋅)𝑅 this step we perform two different (but intertwined) tasks:
standsforarationalone.Hereafter,considerW=(0,100]𝑛+1 dataparameterizationandweightcomputation.Theformer
and D = [0,1]𝑚. Without loss of generality, we can also consists of obtaining a discrete association between the set
assumethat𝑤0 =𝑤𝑛 =1,𝑡1 =0,and𝑡𝑚 =1. ofparameters{𝑡𝑖}𝑚𝑖=1 andthenoisydatapoints{Q𝑖}𝑚𝑖=1 tobe
Figure1 shows the main steps of our method. Basically, fitted, while the latter computes the weights. Both tasks are
itconsistsoffourmajortasks:dataparametrization,weight performedsimultaneouslybyusingthesimulatedannealing
computation,datafitting,anddegreedetermination.Upper approachdescribedintheprevioussection.Theinputforthe
part of this figure summarizes the method: we initially set SAmethodisgivenbythefollowing:
a range for the curve degree 𝑛; then, for each value of this
parameter 𝑛 within that range, we apply a combination of (i)thecurvedegree,𝑛,
simulatedannealingandleast-squaresoptimization(boxSA- (ii)initialrandomparametervectorsp0,w0,andb0,
LSQ) to perform the first three tasks and compute the data
parameters, weights, and poles of the best fitting rational (iii)theenergyfunction,givenby(5),
Be´zier curve for this value of 𝑛. Then, the BIC value of
(iv)aneighborhoodfunction(describedinSection5.2.1),
the resulting curve (corresponding to the last task, degree
determination)isobtained(boxBIC).Atitsturn,theSA-LSQ (v)acoolingschedule(describedinSection5.2.2),
canbedecomposedintotwosteps(middlelayerofFigure1): (vi)astoppingcriterion(describedinSection5.2.3).
SAperformsdataparameterizationandweightcomputation
(seeSection5.2fordetails),whileLSQisusedtocomputethe In this work, two different simulate annealing schemas are
poles (see Section5.3). This combination of SA and LSQ is considered:thesequentialschemaandtheall-in-oneschema.
repeatediterativelyuntilnofurtherimprovementisreached. Basically,theformercalculatesthedifferentsetsofunknowns
Inthispaper,weintroducetwodifferentSAschemas,shown of our optimization problem in a sequential way (i.e., only
graphically in the lower layer of Figure1: the sequential some parameters are initially computed and subsequently
MathematicalProblemsinEngineering 7
usedtocomputetheremainingparameters),whilethelatter where k ∈ D is a random vector holding ‖k‖2 = 1. This
computesallunknownsatonce.Letusanalyzethemindetail. is one of the tested functions in [5] but it required another
Using the notation introduced above, the sequential supportfunction(alocalsearchmethod)inordertoexploit
schema(SEQ)beginsbyfindingthebestnonrationalBe´zier the neighborhood of a solution. In this paper we remove
fitting curve through SA𝐵, then it performs the following the supporting local method by maintaining two sets of
𝑝
sequenceiteratively: controlling parameters: firstly, the global real temperature
followingtheclassicalSAtemperature.Then,thealgorithm
SA𝑅 󳨀→SA𝑅. (7) builds a second set of virtual temperatures (one for each
𝑤 𝑝
spatialdimension)thatarerestartedatthebeginningofeach
until there is no further improvement in the final solution innercycle.Theresultingequationbecomes
of either SA procedure in comparison with the previous x =x +T⊙k, (11)
one. Each routine takes the preceding solution as an input new old
parameter. The schema can be summarized as follows (see where⊙representstheHadamardproductoftwovectorsand
also the lower layer of Figure1, where the numbers in red Thasallitscomponentssettothecurrentsystemtemperature
indicatethedifferentstepsofthealgorithm): atthebeginningofeachoutercycle.Thus,wheneverx is
new
accepted through the acceptance criterion, we compute the
(1)Apply the SA𝐵𝑝 with a random initial guess p0 ∈ D absolute difference between the old and new solutions. For
to find a non-rationalBe´zier curve that fits the data eachcomponentthatbehavesbetterthanthepreviousonein
better.Letp𝐵betheresultingsolution. thecomparison,weupdatethecorrespondingcomponenton
𝑅 Taccordingtothecoolingschedule.Theproposedmethodis
(2)SearchforarationalBe´ziercurvethroughSA with
𝑤
arandominitialguessw0 ∈ Wandfixedparameters basedontheAdaptiveSimulatedAnnealing(ASA)algorithm
(p𝐵).Let(w𝑅,p𝐵)betheresultingsolution. butthereisonekeydifference:ourmethoddoesnotrequire
computing the gradient of the function. A further analysis
(3)Apply the SA𝑅𝑝 optimization with p𝐵 as the initial abouttheperformanceofthisneighborhoodfunctionforthe
guess and fixed weights w𝑅. Let (w𝑅,p𝑅) be the thirdexampleofourbenchmarkisdescribedinSection6.3.
resultingsolution. Thereaderiskindlyreferredtothatsectionforfurtherdetails
onthisissue.
(4)Repeat (2)-(3) iteratively until there is no improve-
mentintheresultingsolution. 5.2.2. Cooling Schedule. By the cooling schedule we refer to
In general, the energy functionfor the simulated annealing atriplet(𝑇0,𝑇,T)accountingfortheselectionoftheinitial
𝑅 temperature, the temperature parameter, and the cooling
procedures SA in steps (2) and (3) above is that in (5). function(T),respectively,alongwiththethermalequilibrium
However, a different energy function is required for the
criterion. The cooling schedule governs the pace at which
𝐵
nonrationalcaseinstep(1),SA ,givenby
the temperature is updated during the execution of the
2 SA. Therefore, its choice is of primary importance for the
𝑚 𝑛
Θ=∑𝜇 (Q −∑b 𝐵𝑛(𝑡)) (8) goodperformanceofthealgorithm.Becauseofthisreason,
𝑖 𝑖 𝑗 𝑗 𝑖 numerous cooling schedules have been proposed over the
𝑖=1 𝑗=0
years,suchasthelinearschedule,thelogarithmicschedule,
thatThcoerraelslp-ionn-donsetoscthheemcaase(A𝑤I𝑗O=)1coforrreaslplo𝑗n.ds to the mini- ganivdenmbanyy𝑇ontehwers=.A𝛼v𝑘e𝑇r0y,cwohmemreotnypoincealilsyt0h.e9po<we𝛼rsc<hed0u.9le9,
mizationoffunctional and 𝑘 > 0 is a scalar parameter. Another popular schedule
is the Boltzmann schedule, given by 𝑇new = 𝑇0/log(𝑘).
𝑓 :W×D󳨀→R,
𝑤,𝑝 (9) Finally, we also consider the fast schedule, governed by the
whereforeachstatevector(w,p)wecompute𝑓(w,p)asthe slaimwp𝑇linceiwty,=spe𝑇e0d/,𝑘a,nwdhgoicohdppreorvfoidrmedanacgeofoordobthaelarndcaetaafimttoinngg
𝑅
least-squaressolutionof(5)throughSA(𝑝,𝑤). problems[5].Inthispaperweconsideredinitiallythepower,
Boltzmann,andfastschedules.However,ourcomputational
5.2.1. Neighborhood Function. The neighborhood function is experiments showed that the fast schedule provides more
one of the key components of the SA algorithm. Further- visuallyappealingresultsandrunsfasterthantheothertwo
more,itsrolebecomesevenmoreimportantformultimodal alternatives.Sincethesearetwoparticularlyvaluablefeatures
optimization problems, where the objective function is of inthecontextofdatafitting,weeventuallyselectedthefast
the many local peaks surrounded by deep valleys type. This scheduleasthebestchoiceforthiswork.
is exactly what happens in this paper. There are several Thecoolingscheduleneedstoaccomplishtwoadditional
alternatives for the neighborhood function. A very popular goals. On one hand, the starting temperature needs to be
choiceisthefastneighborhoodfunction,whichbuildsthenew sufficientlyhighinordertoletthealgorithmfreelyexplorethe
solutionbymodifyingthepreviousoneinstepsproportional searchspaceDbutnotsomuchthatthesystembehaveslikea
tothesystemtemperature: randomsearchforalargenumberofiterations.Ontheother
hand, the cooling function must perform a slow reduction
x =x +𝑇⋅k, (10)
new old ofthesystemtemperaturetopreventprematureconvergence.
8 MathematicalProblemsinEngineering
Thesegoalsareachievedbylettingtheacceptancefunctionto with the role to (possibly) improve the convergence
initiallyacceptworsestatesfrequently,whichinturnrequires ratewithrespecttothestandard(nonelitist)version
hightemperatures.Thereisnotknowngeneralwaytoselect ofSA.
a good starting temperature for all problems. A common
(ii)We add a new operator, related to the domain of
practice consists of selecting a very high temperature to
theproblem,toworkincombinationwiththeneigh-
rapidly cool the system after a certain proportion of worse
borhood function. This extra functionality checks
solutions are accepted (this is called the acceptance ratio);
whether a new generated solution goes outside the
then,begintoslowlycoolthesystem,thusbeginningwiththe
search domain of the problem and sends it back
proper SA algorithm. Another common technique consists
intothesearchspacewheneveritgoesaway.Tothis
of spending most of the computation time in running fast
purpose, we apply the classical cast back operator,
sample runs of the entire algorithm for different starting
a widely accepted routine in numerical methods.
temperatures in order to determine the ideal temperature: Suppose that the SA returns a new solution x
one that promotes better solutions but still manages to outsidetheproblemdomainD,obtainedfromapnreew-
explore the space avoiding local minima. This last practice vioussolutionx withinD.Thecastbackprocedure
wasproposedin[36]asawaytoavoidwastingcomputation replacesx byoladnewvaluex givenbytheconvex
timewhenrunningtheSA.In[35]thereisabrief,although combinatinoenw:x = 𝛼Proj(x c)b+(1−𝛼)x ,where
insightful, overview of the different approaches on how to Proj(⋅)istheopcebratorthatpnreowjectsanypoinoltdoutside
selecttheinitialtemperature.
thedomainontoitsclosestpointontheboundaryof
In order to select the initial temperature we have dis- Dand𝛼isauniformrandomnumberintheinterval
carded the time-consuming restarting strategy used in our (0,1).Thisprocedurereturnsanewpointx thatis
previous work in favor of the method presented in [37]. It cb
wellwithinthesearchdomainwhilesimultaneously
consistsofsettinganacceptanceratio,𝜒0,anddetermining, ensuring that the probability of the boundary is not
throughaniterativealgorithm,acompatiblestartingtemper-
increasedbythisoperator.
ature.Amongtheseveralpossiblewaystodefine𝜒0,inthis
(iii)Finally we improve the cooling schedule with extra
workwechoosetheclassicaloneofthequotientbetweenthe
conditionsforthethermalequilibrium.Inourimple-
numberofbadtransitionsacceptedandtheattemptedones,
alreadyproposedin[38].Accordingtothisdefinition,weset mentation,theinnercycle stopsifthevalueof𝜒0 is
𝜒0 =0.8,atypicalvalueinpreviousliteratureinthefield. reachedafter𝜉iterations.
Thesenewfeaturesimprovetheperformanceofourapproach
5.2.3. Stopping Criteria. Similar to the choice of the initial
significantly in terms of computational time and quality of
temperature,thereisnotageneralsetofstoppingconditions
results.
suitable for all problems. There are, however, two common
practices: the first one is to set a maximum number of 5.3.DataFitting. Withtheparameterizationandweightscal-
iterations; the second one is to stop the execution of the culatedinprevioussteps,wecomputethecurvecoefficients
algorithm when the system is frozen, that is, when no new {b𝑗}𝑛𝑗=0.Using(3),(5)canberewrittenas
solutions (either better or worse states) are accepted for a
predefinednumberofiterations.Sincetheformercanwaste R𝑇R ⋅⋅⋅ R𝑇R b QR
0 0 𝑛 0 0 0
alotofcomputationtimewithnofurtherimprovementson [ ][ ] [ ]
[ . . . ][ . ] [ . ]
thesolution,inpracticethestoppingcriterionisoftenacom- [ .. .. .. ][ .. ]=[ .. ], (12)
[ ][ ] [ ]
bination of the two. In our implementation, the algorithm
stopswheneveroneofthefollowingconditionsismet:either [R0𝑇R𝑛 ⋅⋅⋅ R𝑛𝑇R𝑛][b𝑛] [QR𝑛]
𝑛th𝑓eevatlo>tal10n0u0m⋅b𝜉eorro♯f[𝜕e(vΔagluloabatli)o]n=s 1o0f ⋅th𝜉,ewehneerreg𝑛y𝑓feuvanlcdteinonot,e𝜉s where R𝑗 = (𝑅𝑗𝑛(𝑡1),...,𝑅𝑗𝑛(𝑡𝑚))𝑇 represents the column
denotes the number of free parameters of the system, and vector of rational basis functions given by (4) at the best
𝜕(⋅)denotesthelackofchangesofthemeanvariationofthe parameter values, Q = (Q1,...,Q𝑚) is a row vector, and
energyfunction. (⋅)𝑇 represents the transpose of a vector or matrix. System
(12) is overdetermined, meaning that no analytical solution
5.2.4. Further Improvements. In addition to the previous can be obtained. Instead, we solved it numerically through
SAcomponents,threeotherimportantmodificationsofthe least-squares minimization. If R+ denote the generalized
originalSAalgorithmhavebeenincludedinourapproach: inverse(alsoknownasMoore-Penrosepseudoinverse)ofR=
(R𝑖)𝑇, P = R+Q is the least-squares solution of this data
(i)We improve the memory capacities of the method
fitting problem. From a computational point of view, it can
through elitism: the best state from the current iter-
be obtained by either LU decomposition or singular value
ation is encoded as a vector (x𝑀 ,𝑓𝑀 ) and stored
best best decomposition(SVD). In thiswork,we have used theSVD
in a temporal buffer. Obviously, this “best so far”
solverfromthepopularnumericalprogramMatlab.
solution is updated whenever a better solution is
achieved during SA execution. We remark however 5.4.DegreeDetermination. Thethreepreviousstepsassumed
thatthis(x𝑀 ,𝑓𝑀 )vectorisnotusedtodrivetheSA a given degree for the rational fitting curve. However, the
best best
execution;instead,itisonlyusedasamemoryeffect, optimaldegreeisproblem-dependent,andthereforeweneed
MathematicalProblemsinEngineering 9
a method to compute it. This is not an easy task. It is well each component (𝑥 , 𝑦 , and 𝑧 ). When used, the
mean mean mean
known that increasing the number of poles increases the acronymN.A.standsfornotapplicable.The logo,shoe,and
accuracyofthemodel,becausewehavemoredegreesoffree- torusnamesrefertothescannedlogo,shoeprofile,andcurve
domtoadjustthedata.However,thisprocesscaneventually onatorusdatasetsdescribedbelow,respectively.Thenumber
leadtooverfitting.Inordertopreventthisundesirableeffect, after torus refers to the signal-to-noise ratio (SNR) applied.
somekindoftrade-offbetweenthesetwocompetingfactors Forthenonrationalexamples,alocalsearchwasperformed
(theaccuracyofthefittingandthecomplexityofthemodel) inordertorefinetheSAsolution.Ourresultsshowthegood
isrequired.Toaddress thisissue, in thispaper we compute performanceofthemethodeveninhighlynoisysituations,
the Bayesian Information Criterion (BIC) for the resulting whicharethecommoncaseinreal-worldapplications.
model.TheBICisaninformationcriterionprovidingagood
compromise between data fidelity and model complexity 6.1. Example 1: A Scanned Logo. The first example corre-
[39].Thisisdonebyintroducingapenaltytermforcomplex spondstotheshapeofadigitallyscannedlogo.Thedataset
modelsintothetargetfunction: consists of a set of 190 noisy 2D data points, represented
by black × symbols in Figure2. The figure shows our
BIC=𝜁log(Θ)+𝜉log(𝜁), (13)
experimentalresultsforthethreecasesanalyzed:nonrational
where Θ, 𝜁, and 𝜉 refer to the energy function, number case (top), the all-in-one rational case (middle), and the
sequentialrationalcase(bottom).Thefiguresontheleftshow
of sampled points, and the number of free parameters,
respectively.Notethat,fixing𝜁and𝜉,theBICbehaveslikethe the reconstructed points, represented as red empty circles.
errorfunctionandifwefixavalueofΘand𝜁,thecriterion Ontheright,thebestfittingcurveisdisplayedasaredsolid
line. This example has been chosen because it represents
penalizesthosemodelswithahighernumberofparameters.
a common real-world scenario: a scanned figure with the
Therefore, BIC provides us with a powerful procedure to
typical noise introduced during the scanning process. In
computetheoptimalvalueforthedegreeofthefittingcurve.
It always corresponds to the parameter value 𝑛 with the addition to the high-intensity noise (clearly visible in all
instancesofFigure2),thisshapeisalsochallengingbecause
smallestvaluefortheBIC.
it includes difficult geometric features, such as several self-
6.ExperimentalResults intersectionsandstrongchangesofslopeandcurvature.
Asthereadercanseefromthefigures,themethodisable
Inthissection,weanalyzetheperformanceofourmethodby to recover the general shape of the data points with good
applyingittoabenchmarkofthreeillustrativeexamplesof2D accuracy.Thisisaveryremarkableresultbecausetheoriginal
and3Dnoisydatapoints.Theseexampleshavebeencarefully data points are highly noisy. Best results correspond to the
chosensothattheyexhibitchallengingfeaturessuchasself- sequential rational schema, while the nonrational and the
intersections or strong changes of slopes and curvatures. AIOrationalschemasperformalmostsimilarlyinthiscase.
The first two examples correspond to real-world instances Thisfactisclearlyvisibleinthetopmostloopofthefigures
so that we can replicate the usual conditions of real-world in right column, best fitted through the sequential schema
applications,includingthepresenceofnoiseoflow-medium (bottom figure). These visual results are in good agreement
intensity. The last one is an academic example designed to withthenumericalresultsreportedinTable1.Thebestand
analyzetheeffectofdifferentlevelsofnoiseonourmethod. average BIC for the sequential schema are approximately
For each example, we report the results of three different −770and−701,respectively,whiletheyareboth−691forthe
schemasusedforcomparativepurposes:nonrational,ratio- nonrationalschemaand−661and−658fortheAIOrational
nal all-in-one, and rational sequential. Finally we analyze schema. Note also that the three methods obtain the same
the robustness of the method in the presence of noise by optimalnumberofpoles𝑛 = 13.Inotherwords,thegood
pol
comparingeachschemaagainstthesamedataperturbedwith accuracyofourmethodisnotattheexpenseofaverylarge
differentlevelsofnoiseforthelastexample. numberoffreevariables.
Regardingtheimplementationissues,alltheexperiments
were run on a AMD-FX-4100 Quad-Core Processor at 6.2. Example 2: Shoe Profile. The second example corre-
3600Mhzwith8GBDDR3RAMrunningLinux3.14.xLTS spondstoashoeprofileobtainedfromapressure-mechanical
kernelandMATLAB2012a.Foreachdatasetandschema,the method without filtering, leading to a set of 400 three-
experimenthasbeenexecuted26times.Then,20executions dimensional noisy data points. Figure3 shows our experi-
are finally selected (after removing the three bast and three mentalresults.Theinterpretationofthisfigureissimilarto
worst executions) to provide statistical evidence for the that of previous example and, hence, it is omitted here to
resultspresentedandasserttheexperimentreproducibility. avoidredundantinformation.Onceagain,thebestfittingis
Table1 reports our results for each dataset and experi- obtainedwiththesequentialrationalschema,althoughinthis
ment.Thefollowingitemsarearranged(incolumns):dataset case,thevisualandnumericalresultsarecloserforthetwo
examined,thetypeofcurvereconstructed(NR:nonrational; rationalschemasandsignificantlyworseforthenonrational
R: rational), the schema executed for rational curves (AIO: one. Note also that we obtained a similar parameter value,
all-in-one; SEQ: sequential), the total number of calls to 𝑛 =24,forthenumberofpolesinallcases.
pol
theenergyfunctioninthebestcase(representedby𝑛𝑓eval),
the best and average BIC, the number of poles for the best 6.3. Example 3: Curve on a Torus. Last instance in our
BIC (represented by 𝑛 ), and the relative mean error for benchmarkcorrespondstoanacademicexample.Ithasbeen
pol
10 MathematicalProblemsinEngineering
Table1:Experimentalresultsonthethreeexamplesinourbenchmarkforthesimulatedannealingschemasdiscussedinthispaper.
Model Statisticalresults
Dataset Type Schema 𝑛𝑓eval BIC(best) BIC(avg) 𝑛pol 𝑥mean 𝑦mean 𝑧mean
logo NR 10010 −691.913 −691.830 13 0.0005 0.0001 N.A.
logo R AIO 4208 −661.908 −658.804 13 0.0002 0.0003 N.A.
logo R SEQ 13031 −770.212 −701.338 13 0.0001 0.0001 N.A.
shoe NR 16624 2915.200 3020.101 24 0.0191 0.0152 0.0076
shoe R AIO 6648 2098.859 2470.634 24 0.0012 0.0061 0.0053
shoe R SEQ 20202 1996.226 2289.712 24 0.0024 0.0061 0.0052
torus NR 18546 −2204.821 −2301.567 20 1.06𝑒−05 1.23𝑒−05 2.12𝑒−05
torus R AIO 1529 −6067.465 −6001.324 21 6.42𝑒−06 5.39𝑒−06 4.32𝑒−06
torus R SEQ 5286 −7798.643 −7796.846 20 1.48𝑒−06 1.48𝑒−06 1.35𝑒−06
torus10 NR 13276 4307.409 4390.492 17 0.0311 0.0342 0.0458
torus10 R AIO 3768 4354.618 4414.523 19 0.0300 0.0304 0.0401
torus10 R SEQ 23119 4341.604 4381.367 17 0.0312 0.0323 0.0402
torus20 NR 12027 3887.839 4991.823 20 0.0222 0.0201 0.0364
torus20 R AIO 3778 3474.463 3479.345 20 0.0130 0.0137 0.0181
torus20 R SEQ 45775 3446.804 3465.673 20 0.0128 0.0139 0.0182
torus30 NR 16174 3033.538 3561.901 20 0.0109 0.0094 0.0172
torus30 R AIO 5546 3060.624 3066.298 20 0.0097 0.0089 0.0151
torus30 R SEQ 48511 3021.939 3035.634 20 0.0098 0.0089 0.0148
torus40 NR 11015 6662.487 6700.738 40 0.0183 0.0144 0.0346
torus40 R AIO 5111 2740.208 2788.422 20 0.0074 0.0073 0.0094
torus40 R SEQ 53640 2708.680 2729.980 20 0.0074 0.0072 0.0096
torus50 NR 14101 2386.160 2501.900 25 0.0057 0.0059 0.0076
torus50 R AIO 4205 2357.285 2408.811 21 0.0054 0.0053 0.0076
torus50 R SEQ 45206 2344.050 2378.280 20 0.0055 0.0052 0.0076
torus60 NR 12052 3574.142 4713.190 20 0.0179 0.0121 0.0340
torus60 R AIO 4012 2116.800 2145.321 20 0.0041 0.0045 0.0620
torus60 R SEQ 45191 2051.304 2080.235 20 0.0031 0.0033 0.0038
torus70 NR 8860 3574.142 3866.889 21 0.0058 0.0057 0.0116
torus70 R AIO 4205 2210.012 2284.590 21 0.0039 0.0040 0.0054
torus70 R SEQ 45797 1963.317 1995.759 20 0.0038 0.0039 0.0056
torus80 NR 14100 3574.142 3934.296 25 0.0038 0.0044 0.0052
torus80 R AIO 4002 1814.950 2008.400 21 0.0032 0.0036 0.0044
torus80 R SEQ 13366 1803.356 1931.721 20 0.0032 0.0036 0.0043
torus90 NR 14100 1828.028 1981.534 25 0.0031 0.0041 0.0048
torus90 R AIO 4204 1671.145 1799.627 21 0.0030 0.0030 0.0042
torus90 R SEQ 13350 1654.944 1786.873 20 0.0030 0.0031 0.0044
carefullydesignedtoanalyzetheperformanceofourmethod in Table1, is then perturbed with additive white noise of
againstnoiseofdifferentintensities.Tothisaim,weconsider different intensities, modulated by a signal-to-noise ratio
theparametric: (SNR) ranging from 10 (very high intensity) to 90 (low
𝑥(𝑡)=[7+2cos(5𝑡)]cos(2𝑡) intensity),withstep-size10.Thecorrespondingdatasetsare
labelled as torusN, where 𝑁 indicates the SNR intensity.
𝑦(𝑡)=[7+2cos(5𝑡)]sin(2𝑡) The simulation results with our method for the resulting
(14) 10 datasets are reported in the last 10 horizontal blocks of
𝑧(𝑡)=3sin(5𝑡) Table1.Someimportantobservationscanbeobtainedfrom
𝑡∈[0,2𝜋], the numerical data in the table. The most important one is
thatthesequentialrationalschemaoutperformstheothersin
whichcorrespondstoacurveonatorus.Weconsideraset termsofBICvalue,meaningthatitprovidesthebesttrade-
of202three-dimensionaldatapointswithuniformsampling offbetweenaccuracyandcomplexityforallinstancesinthis
intheintervaldomain[0,2𝜋].Thisdataset,labelledastorus example. According to our results, the optimal number of