Smoothing spline fit through a set of (x,y) points

[hanshinlee77]

Hello, I have a set of (x,y) points (shown in the attached dpng) which I would like to create a smoothing spline through. Normally, I don’t have a problem doing a fit like this using Spline1D() from Dierckx.jl or interpolate() from BSplineKit.jl. The difficulty I have with this particular data set is that the points appear to follow an underlying curve that steers to the left near the top, making the left side of the profile concave. This seems to keep my go-to spline methods confused, returning a completely wrong fit, particularly that left side of the profile. I was wondering if there are ways to use the existing methods to do this fit or is there a Julia package that can handle this kind of fit?

I appreciate any help.

Data points are:
-51.88380131158897,0.5912681519319563
-48.50922599109117,0.0
-48.43237308360659,1.7738044557970056
-43.014779032608885,0.0
-42.32448040957922,0.0
-42.295198627223954,0.5912681519320131
-41.72451368217037,0.0
-40.17264832268097,0.0
-40.00866373138865,0.0
-39.9563486093366,1.7738044557970056
-36.683862952446134,1.7738044557970056
-35.82053073832424,216.995411759141
-35.82053073832424,228.22950664585403
-35.229262586392224,241.237405988364
-35.199980804036954,1.1825363038649925
-34.833205589164606,1.7738044557970056
-34.63799443446021,251.28896457121198
-34.046726282527175,202.804976112767
-33.45545813059516,1.7738044557970056
-33.426176348240006,225.27316588619306
-32.86418997866326,261.340523154061
-32.486162347560935,1.1825363038639694
-32.27292182673125,188.614540466392
-31.95371773698821,4.138877063526024
-31.78068414589893,1.7738044557969488
-31.362449585055174,224.090629582328
-31.13964145240857,1.7738044557970056
-31.087326330356746,2.956340759660975
-31.06110374051093,203.39624426469902
-31.06110374051093,214.63033915141204
-31.06110374051093,236.50726077290597
-30.77118143312316,198.66609904924098
-30.77118143312316,212.85653469561498
-30.499117370934187,176.19790927581505
-30.490418795457458,0.5912681519320131
-30.46983558857903,189.79707677025704
-30.179913281191148,187.43200416252796
-29.878567436645994,176.78917742774706
-28.997376977326212,162.00747362944003
-28.997376977326212,173.83283666808597
-28.997376977326212,232.36838370938005
-28.938553435967947,202.21370796083397
-28.938553435967947,214.03907099948003
-28.920524069841576,203.98751241663098
-28.725312915137238,164.37254623717
-28.696031132781968,163.78127808523698
-28.313876963340988,202.804976112767
-28.313876963340988,213.44780284754802
-28.313876963340988,224.68189773426104
-28.233075234298894,212.85653469561493
-28.130985570696794,245.37628305189003
-27.814840673462186,9.460290430916018
-27.737987765976584,186.24946785866297
-27.722608811408975,236.50726077290608
-27.64180708236904,225.864434038125
-27.54277661127321,3.5476089115940113
-27.54277661127321,147.81703798306603
-27.539717418766713,234.73345631710902
-27.513494828917032,243.60247859609296
-27.22357252152915,146.04323352726897
-27.164748980170884,4.730145215457981
-27.164748980170884,225.86443403812495
-27.14671961404457,7.095217823186999
-27.14671961404457,172.05903221228897
-27.14671961404457,214.630339151412
-27.13134065947696,190.38834492219007
-27.05053893043896,199.25736720117294
-26.942809883863447,7.686485975119012
-26.92222667698502,1.7738044557970056
-26.92222667698502,148.40830613499804
-26.632304369597136,132.44406603282698
-26.57348082823887,186.84073601059595
-26.45927077850888,153.13845135045602
-26.45927077850888,186.24946785866297
-26.45927077850888,201.03117165697
-26.4094962369486,4.138877063526024
-26.4094962369486,235.32472446904097
-26.360240307408162,131.85279788089497
-26.360240307408162,266.661936521451
-26.357181114896775,252.47150087507703
-26.351541731931434,215.221607303344
-26.041036217665123,16.555508254103017
-26.041036217665123,118.25363038645298
-25.98221267630595,154.32098765432102
-25.98221267630595,170.28522775649196
-25.98221267630595,233.55092001324397
-25.9641833101806,140.72182015987897
-25.94880435561197,177.38044557967999
-25.868002626569023,4.1388770635259675
-25.868002626569023,63.85696040868396
-25.868002626569023,135.99167494442094
-25.868002626569023,167.32888699683093
-25.818228085018518,190.979613074121
-25.818228085018518,222.90809327846404
-25.765912962966695,219.95175251880198
-25.76027357999942,204.578780568563
-25.76027357999942,228.22950664585397
-25.372915158247565,156.09479211011802
-25.372915158247565,225.273165886193
-25.276734474638943,79.82120051085593
-25.226959933088665,247.74135565961905
-25.174644811036615,198.07483089730795
-25.169005428066384,238.28106522870206
-25.148422221187957,120.61870299418206
-24.858499913800188,239.46360153256705
-24.68546632270909,49.07525661037795
-24.68546632270909,94.60290430916194
-24.63569178114858,88.69022278983999
-24.63569178114858,142.49562461567604
-24.63569178114858,157.86859656591497
-24.63569178114858,173.24156851615402
-24.586435851611213,9.460290430916018
-24.58337665910676,7.686485975119012
-24.58337665910676,260.74925500212794
-24.557154069255944,133.62660233669203
-24.267231761868175,104.06319474007898
-24.20840822050991,14.190435646373999
-24.20840822050991,110.56714441133295
-24.20840822050991,124.75758005770797
-24.20840822050991,136.582943096353
-24.190378854383596,26.607066836951958
-24.190378854383596,109.97587625940099
-24.190378854383596,126.53138451350497
-24.0444236292185,8.277754127051992
-24.0444236292185,125.94011636157205
-24.0444236292185,207.535121328225
-23.9951676996792,18.329312709900023
-23.9951676996792,116.47982593065603
-23.992108507166677,147.81703798306603
-23.992108507166677,183.29312709900194
-23.986469124202472,187.43200416252802
-23.96588591732393,5.321413367389994
-23.67596360993616,91.64656354950102
-23.617140068576873,44.34511139491997
-23.599110702451583,16.555508254103984
-23.599110702451583,54.396669977768
-23.599110702451583,90.46402724563694
-23.58373174788295,157.86859656591503
-23.58373174788295,245.96755120382204
-23.502930018838924,122.39250744997895
-23.502930018838924,239.463601532567
-23.418176665341946,0.14190435646401056
-23.403899547747187,88.69022278984005
-23.39520097227046,165.55508254103398
-23.374617765391008,106.42826734780806
-23.08469545800415,27.198334988884028
-23.025871916644974,58.53554704129397
-23.025871916644974,94.01163615722999
-23.007842550518603,39.61496617946199
-23.007842550518603,69.17837377607503
-22.91166186690907,16.55550825410296
-22.91166186690907,31.337212052409996
-22.861887325358566,69.17837377607503
-22.861887325358566,105.83699919587497
-22.81263139581415,33.111016508207
-22.81263139581415,72.72598268766899
-22.493427306071226,37.84116172366504
-22.493427306071226,53.21413367390403
-22.493427306071226,66.222033016414
-22.43460376471296,70.360910079939
-22.43460376471296,81.59500496665203
-22.40119544401898,139.539283856015
-22.32039371497899,107.01953549974002
-22.270619173428486,18.92058086183198
-22.270619173428486,29.56340759661299
-22.270619173428486,47.301452154581
-22.22136324388225,104.06319474007898
-22.218304051376663,133.035334184759
-22.218304051376663,162.59874178137295
-22.21266466840541,17.14677640603503
-22.21266466840541,47.301452154581
-22.21266466840541,65.63076486448102
-22.21266466840541,115.29728962679104
-22.21266466840541,129.48772527316606
-22.21266466840541,146.63450167920104
-22.192081461526982,91.05529539756901
-21.902159154139213,77.45612790312703
-21.843335612780947,27.19833498888397
-21.63009509195024,47.89272030651307
-21.63009509195024,60.30935149709103
-21.627035899436805,63.85696040868402
-21.627035899436805,98.15051322075601
-21.60081330959497,12.416631190578016
-21.25206746084791,238.281065228702
-21.21865914015399,118.84489853838602
-21.088082869558548,255.42784163473794
-21.030128364540474,32.51974835627402
-21.030128364540474,84.55134572631397
-21.030128364540474,102.28939028428198
-21.030128364540474,242.41994229222797
-21.009545157661933,23.059457925358004
-20.642769942789585,230.59457925358294
-20.444499595576644,17.73804455796801
-20.444499595576644,78.63866420699094
-20.444499595576644,112.93221701906202
-20.418277005730033,47.30145215458106
-20.418277005730033,59.12681519322598
-20.418277005730033,75.68232344732996
-20.012472110211945,7.804739605506029
-20.012472110211945,69.95884773662601
-20.012472110211945,101.46161487157701
-19.856290636153176,266.070668369519
-19.85323144364679,31.92848020434201
-19.85323144364679,49.075256610378005
-19.235740701866007,34.29355281207103
-19.235740701866007,243.60247859609296
-19.161045971429985,50.37604654462905
-19.161045971429985,84.43309209592803
-18.3545502425452,243.60247859609296
-18.309619832646945,35.05037604654501
-18.181516651448874,244.19374674802492
-17.458193693864985,23.98183624237305
-17.079782076627964,245.96755120382204
-16.308681724560188,261.340523154061
-15.714354380116674,263.70559576178994
-14.138820271534598,231.18584740551597
-13.934910541353474,243.01121044415999
-13.565581485728899,238.281065228702
-12.761072812966177,252.47150087507703
-12.73179103061102,242.41994229222803
-12.441868723223138,243.01121044416107
-12.349636861169984,235.91599262097407
-12.219060590578465,258.97545054633196
-11.200508877999937,233.55092001324397
-10.957986574813958,234.73345631710902
-9.903762524398871,244.19374674802492
-9.804732053305202,235.91599262097304
-9.80167286079677,263.114327609857
-9.40867505607656,229.41204294971902
-8.894259811629126,237.68979707677
-8.21075979764396,222.31682512653208
-8.030927597508253,219.36048436686997
-8.001645815152983,222.90809327846398
-7.0616318144739125,219.95175251880198
-6.839692718166475,241.82867414029596
-6.529187203900165,228.82077479778604
-6.227841359355921,208.12638948015706
-5.665854989779177,206.94385317629303
-5.636573207424021,196.89229459344403
-5.346650900036138,216.404143607209
-5.269797992550593,222.90809327846404
-5.269797992550593,215.221607303344
-5.071527645336801,259.566718698264
-5.045305055490985,182.70185894707004
-4.532574615458543,254.24530533087403
-4.4833186859142415,189.79707677025704
-3.4806145821859786,203.98751241663206
-3.3998128531488874,241.23740598836395
-3.3007823820502153,172.65030036422104
-3.2715005996949458,167.328886996831
-2.706455037606702,254.24530533087398
-2.3903101403741402,195.70975828958
-2.331486599015875,203.39624426469902
-2.3134572328895615,203.39624426469896
-2.1675020077286717,238.87233338063498
-2.118246078185166,157.277328413982
-1.706810126388973,192.162149377986
-1.5182793507764245,235.32472446904097
-1.5182793507764245,221.13428882266697
-1.1309209290255922,180.33678633934102
-0.9849657038585065,222.90809327846404
-0.9357097743212535,143.086892767608
-0.9064279919649607,153.138451350457
-0.6165056845782146,181.519322643205
-0.5576821432189263,190.38834492218894
-0.5396527770925559,192.75341752991795
0.03358600871308681,171.467764060357
0.051615374839400374,167.32888699683093
0.06699432940700945,175.60664112388304
0.14779605845092192,227.04697034198892
0.24682652954379591,122.98377560191102
0.24988572205325,232.95965186131195
0.2761083118990655,140.72182015987903
0.8411538739833304,247.741355659619
0.8467932569525374,203.39624426469896
0.8673764638310786,128.30518896930107
1.2495306332720588,160.233669173644
1.3303323623110828,206.35258502435994
1.3303323623110828,188.61454046639193
1.4324220259131835,220.543020670734
1.7485669231517704,164.96381438910203
1.840798785204015,144.26942907147304
2.3398350750837835,147.22576983113402
2.3398350750837835,131.26152972896307
2.416687982568419,150.18211059079493
2.611899137272758,109.38460810746903
2.6205977127496,188.02327231446003
2.6411809196280274,116.47982593065603
2.989926768374062,156.68606026204998
3.104136818111101,176.197909275815
3.104136818111101,173.83283666808592
3.104136818111101,146.04323352726902
3.153911359661379,200.439903505038
3.2324490715600405,107.01953549974007
3.5992242864334116,137.17421124828496
3.695404970040954,159.642401021711
3.797494633643282,198.07483089730795
4.113639530880846,114.70602147485897
4.205871392933034,126.53138451350503
4.205871392933034,109.38460810746903
4.336447663531317,180.33678633934102
4.388762785583367,183.88439525093395
4.394402168546549,169.69395960456
4.781760590297438,124.75758005770803
4.985670320478562,156.68606026205003
5.354999376104047,141.31308831181093
5.354999376104047,128.30518896930096
5.4692094258409725,128.305188969301
5.5682398969338465,93.42036800529797
5.597521679289002,91.64656354950102
6.060477577771053,116.47982593065598
6.110252119321331,157.27732841398205
6.168206624343611,133.03533418475905
6.478712138609808,95.19417246109504
6.5709440006620525,92.23783170143406
6.7015202712614155,135.99167494442105
6.7015202712614155,122.39250744997895
6.753835393303234,148.99957428693097
6.780057983154052,78.04739605505904
7.128803831900086,116.47982593065603
7.1468331980264,112.34094886712995
7.345103545243319,163.19000993330496
8.311340135765136,99.92431767655296
8.344748456459058,78.63866420699208
8.344748456459058,66.81330116834602
8.524580656595845,75.68232344732996
8.533279232072573,118.25363038645304
9.06659287898151,108.79333995553696
9.511905805755418,98.74178137268797
9.71017615297319,130.67026157703003
9.71581553593353,86.91641833404304
9.736398742815027,63.265692256752004
10.026321050202796,76.86485975119399
10.085144591562084,79.82120051085599
10.118552912256007,49.075256610378005
10.199354641300943,99.92431767655302
10.898351839793577,102.28939028428198
11.26768089542611,67.40456932027803
11.431665486711381,91.05529539756895
11.483980608763204,115.88855777872402
11.510203198611975,45.527647698783994
12.075248760703289,101.698122132349
12.075248760703289,76.86485975119399
12.663457720121755,56.76174258549702
12.66651691263337,88.09895463790701
12.692739502477025,29.56340759661299
12.982661809864794,59.718083345159016
13.074893671917039,34.88482096400401
13.155695400961122,87.50768648597494
13.257785064563222,62.67442410481999
13.263424447523562,69.76964192800705
13.845994023985781,40.206234331394
13.845994023985781,25.424530533087022
13.849053216493303,46.11891585071601
14.242051021213399,85.733882030178
14.257429975782031,23.65072607729104
14.445960751393613,51.440329218107024
14.979274398311418,76.86485975119405
15.03158952036324,32.519748356273965
15.406557958952135,49.075256610378005
15.424587325078448,75.09105529539801
15.520768008690993,72.13471453573601
15.520768008690993,72.13471453573601
15.649080262138,16.555508254103017
15.939002569525769,48.48398845844599
16.015855477010405,62.674424104819934
16.607123628942418,50.84906106617501
16.607123628942418,35.47608911593596
16.802334783646756,11.825363038645037
17.121538873389795,30.74594390047804
17.39666212809334,18.329312709899966
17.402301511053565,34.29355281207103
17.78965993280741,22.46818977342599
17.80503888737502,10.051558582849054
17.93561515797137,63.26569225675195
18.477108768350945,58.53554704129402
18.52688330990145,50.84906106617501
18.954166870546032,35.47608911593596
19.48661148111978,14.190435646373999
19.709419613761384,35.47608911593602
20.25091322414096,43.162575091054975
20.35864227071363,24.241994229222996
20.72797132634207,23.059457925358004
20.746000692468442,9.460290430915961
20.89195591763132,21.876921621493977
21.532479999104794,5.321413367389994
22.024717679940977,23.059457925357947
22.535184102833057,1.7738044557970056
23.309343647413243,8.869022278984005
24.275580237936083,9.460290430915961
24.497519334243634,7.686485975119012
24.51810254112206,1.7738044557970056
25.399293000441844,5.321413367389994
26.163594743471094,11.82536303864498
26.213369285021372,11.82536303864498
26.658682211791415,4.730145215457981
28.630757014803294,3.547608911592988
29.760978196611404,5.321413367389994
30.77952990919107,3.5476089115939544
30.99277043002178,1.7738044557970056
31.99547453375004,0.0
35.03261702245095,1.7738044557969488
35.14034606902351,1.7738044557970056
36.11897264270743,0.5912681519319563
36.343465579766985,0.0
37.8159241910198,1.1825363038649925
39.273583749583395,0.5912681519320131
39.64855218817513,1.7738044557970056
42.177609387191524,1.7738044557970056
43.22956942046301,0.5912681519330363
43.41810019607351,1.1825363038640262
44.98799492169144,0.0
45.77447422832779,0.0
46.26671190916113,0.0
49.700110771023105,1.1825363038639694
50.53390122614201,0.0
50.823823533529776,0.0
50.91605539558202,0.5912681519330363
51.0989467882232,0.0
53.41170427390148,1.7738044557970056
55.243463234713545,1.1825363038640262
58.587162416133424,1.1825363038649925
58.68334309974102,0.0
73.46504689805101,1.7738044557969488

[JM_Beckers]

My guess is that you have sorted the point on the x-coordinate ? If you have the points in their “natural” order, I think a parametric analysis could do the trick (fit x(s) and y(s) interdependently and then plot the result). An example where I tried to sort the points to get a “natural order” is included, but if you have their original order, the idea might work better

using DelimitedFiles
using Plots
using DIVAnd
using Statistics
xy=readdlm("xyt.txt",',')
xd=xy[:,1]
yd=xy[:,2]
# Sort the points 
xscale=maximum(xd)
yscale=maximum(yd)

xd=xd./xscale
yd=yd./yscale
dist=(xd.+3).^2 +(yd.+0.15).^2 
indexes=sortperm(dist)

trange=size(xy,1)
NP=10*trange
tgrid=ndgrid(range(0,stop=trange+1,length=NP+2))

mask = trues(NP+2);
mask[1]=false
mask[end]=false
pmn=ones(Float64,NP+2)
t=float.(collect(1:size(xd,1)))
xmean=mean(xd)
ymean=mean(yd)
xgrid,s=DIVAndrun(mask,(pmn,),tgrid,(t,),xd[indexes].-xmean,250.0,1.1;alpha=[0,2,1]);
ygrid,s=DIVAndrun(mask,(pmn,),tgrid,(t,),yd[indexes].-ymean,250.0,1.1;alpha=[0,2,1]);

scatter(xd,yd)
plot!(xgrid.+xmean,ygrid.+ymean)

[rafael.guerra]

Please give it a try to Dierckx’s smooth parametric splines. See this example.

[stevengj]

@hanshinlee77

Hello, I have a set of (x,y) points (shown in the attached dpng) which I would like to create a smoothing spline through.

Many of these methods (e.g. Dierckx) require you to supply a parameterization (t_i,x_i,y_i) of your points. If your points are roughly “in order”, then the simplest parameterization is (i,x_i,y_i). This recent paper Zhong (2022) references a number of algorithms for this problem (and proposes a new algorithm), though it doesn’t provide code.

If you don’t have a parameterization, or even an ordering, it seems like a much harder problem. See this stackoverflow question. There are some suggestions of using a nearest-neighbor ordering. Fang and Gossard (1995) is an older paper addressing precisely the problem of fitting unordered & noisy points, and there are more recent papers on similar topics such as Dinh et al (2002) and Wang et al (2014) that may be worth looking at.

[dlakelan]

Interesting problem. I’m thinking create a 2D log density function and maximize the log-likelihood of the points. Then create a spline/approxfun or similar and maximize the integral of the above log-likelihood along the contour.

[hanshinlee77]

Thank you all for chiming in on this problem. I ended up ditching spline fit and fitting a generalized gaussian function with a skew angle to the data points to get a reasonable fit. The plot below is upside down compared to the original plot. The fit parameters were center, width, depth, \gamma (profile shape), and a skew angle.

I originally wanted to use somewhat non-parametric approach to this fit and think that the suggestion by @JM_Beckers would work. In the past, I had a similar problem, for which I created a 2D histogram of the points and ended up tracing the peaks of that histogram (or density map) to build a 2D curve which I then spline fit. I am guessing that this might be similar to what @dlakelan suggested. Unlike the behavior of the points in the plot of the original post though, that trace did not curve back in x so a spline fit was just fine.

Since spline is a collection of basis functions, I was thinking that this is a matter of finding coefficients of those basis functions that optimally fit the given (x,y) coordinate points. I still need to review the paper suggested by @stevengj. Hopefully something is there!

Thanks all again. I will get back to this post if I have a break through!

[dlakelan]

@hanshinlee77

tracing the peaks of that histogram (or density map) to build a 2D curve which I then spline fit. I am guessing that this might be similar to what @dlakelan suggested.

That’s very similar to what I suggested, except my suggestion basically generalized that a bit.

Since your thing isn’t a function of x you’d want to use a parametric function (x(t),y(t)) each of the x(t) and y(t) functions would be their own spline… then the question is “how do you pick one”. You want it to be “near” the points… you could do something like an “orthogonal” distance but it’s not trivial to figure out which point on the curve to compare each data point to. Hence my suggestion to make a density function p(x,y) which is itself a parameterized function and chosen to maximize the likelihood of the data points then choose (x(t),y(t)) on your spline to maximize the integrated p() over the path.

[hanshinlee77]

Yes, it was indeed a tricky problem I encountered while thinking about this: how to measure the distance (as a merit of a fit) from a spline curve to those points. Your density function approach makes sense and is the one I will be trying. Appreciate it!