Approximate a shape outline using constrained B-splines

Question

I'm looking for a possibility to generate a constrained spline in order to approximate a shape (in my case, a footprint outline). As raw data, I have a table with several hundred xy-coordinate pairs, which have been collected from the boundary of the footprint. The spline should only approximate the data points (the spline does not need to pass the data points). I want to be able to smooth the spline to certain degrees. Also, I need to be able to constrain the spline: Defining several critical data points which the spline has to pass.

The R package "cobs" (COnstrained B-Splines, https://cran.r-project.org/web/packages/cobs/index.html) comes very close to providing a solution, offering parameters to constrain the spline as wanted. However, this package does not spline through an ordered sequence of data points, which of course is crucial when you want the spline to follow the boundary of a shape. I tried to spline x and y coordinates separately, but after recombining them two distinct shapes appear in the plot, so this does not seem to work (Or I got something wrong?). Is anybody aware of a solution?

Update: working example (dinosaur footprint outline)

data.txt:

structure(list(V1 = c(124.9, 86.44, 97.22, 81.34, 49.09, 57.18, 
-77.6, -191.95, -284.67, -383.18, -379.27, -492.85, -547.72, 
-600.67, -713.29, -814.36, -868.27, -926.99, -958.76, -1025.18, 
-1077.16, -1105.07, -1126.25, -1112.77, -1087.74, -989.56, -911.59, 
-859.61, -745.06, -656.5, -682.01, -637.25, -601.71, -539.09, 
-394.79, -219.17, -170.17, -201.48, -122.52, -43.56, 127.97, 
344.42, 539.09, 686.11, 987.63, 1253.31, 1283.15, 1536.32, 1741.14, 
1832.35, 1700.3, 1787.43, 1911.31, 2017.49, 2097.81, 2135.93, 
2093.73, 2066.96, 2063.78, 2022.94, 1978.69, 1919.44, 1904.03, 
1895.37, 1854.22, 1810.23, 1771.09, 1741.48, 1642.45, 1553.96, 
1472.96, 1396.04, 1141.65, 1085.82, 1055.02, 1358.24, 1325.94, 
1031.91, 1287.14, 1265.36, 931.15, 872.12, 811.48, 755.65, 738.32, 
697.41, 682.49, 647.35, 628.25, 620.09, 629.62, 675.22, 709.25, 
718.78, 717.42, 551.09, 535.21, 540.98, 534.73, 546.76, 811.96, 
823.03, 822.07, 607.4, 626.18, 637.73, 659.87, 756.13, 753.72, 
735.91, 720.99, 676.71, 576.6, 508.26, 339.8, 179.53, 121.16, 
45.6, 12.93, -9.87, -12.59, 16, 27.91, 37.78, 49.35, 8.51, 2.72, 
-1.02, 59.22, 58.2, 51.73, 54.45, 0.96, 10.59, 138.62, 149.69, 
144.87, 142.26, 146.34, 125.24, 124.9, 86.44, 97.22, 81.34, 49.09, 
57.18, -77.6, -191.95, -284.67, -383.18, -379.27, -492.85, -547.72, 
-600.67, -713.29, -814.36, -868.27, -926.99, -958.76, -1025.18, 
-1077.16, -1105.07, -1126.25, -1112.77, -1087.74, -989.56, -911.59, 
-859.61, -745.06, -656.5, -682.01, -637.25, -601.71, -539.09, 
-394.79, -219.17, -170.17, -201.48, -122.52, -43.56, 127.97, 
344.42, 539.09, 686.11, 987.63, 1253.31, 1283.15, 1536.32, 1741.14, 
1832.35, 1700.3, 1787.43, 1911.31, 2017.49, 2097.81, 2135.93, 
2093.73, 2066.96, 2063.78, 2022.94, 1978.69, 1919.44, 1904.03, 
1895.37, 1854.22, 1810.23, 1771.09, 1741.48, 1642.45, 1553.96, 
1472.96, 1396.04, 1141.65, 1085.82, 1055.02, 1358.24, 1325.94, 
1031.91, 1287.14, 1265.36, 931.15, 872.12, 811.48, 755.65, 738.32, 
697.41, 682.49, 647.35, 628.25, 620.09, 629.62, 675.22, 709.25, 
718.78, 717.42, 551.09, 535.21, 540.98, 534.73, 546.76, 811.96, 
823.03, 822.07, 607.4, 626.18, 637.73, 659.87, 756.13, 753.72, 
735.91, 720.99, 676.71, 576.6, 508.26, 339.8, 179.53, 121.16, 
45.6, 12.93, -9.87, -12.59, 16, 27.91, 37.78, 49.35, 8.51, 2.72, 
-1.02, 59.22, 58.2, 51.73, 54.45, 0.96, 10.59, 138.62, 149.69, 
144.87, 142.26, 146.34, 125.24), V2 = c(-446.8, -415.83, -394.43, 
-259.19, -104.69, -4.03, 58.59, -80.26, 52.11, -48.33, -142.23, 
-176.89, -233.68, -321.28, -416.57, -457.97, -458.93, -429.09, 
-422.35, -450.27, -431.98, -379.03, -260.63, -123.94, -2.65, 
269.76, 455.55, 548.92, 616.3, 691.38, 756.84, 888.72, 1016.97, 
1157.18, 1198.02, 1101.37, 1025.14, 929.84, 852.25, 766.48, 717.47, 
733.81, 784.18, 835.91, 1225.63, 1198.68, 925.3, 742.4, 814.13, 
732.45, 586.79, 394.84, 212.42, 28.64, -111.58, -337.56, -490.03, 
-526.07, -528.82, -547.2, -551.97, -552.3, -585.51, -551.34, 
-543.16, -526.1, -494.11, -466.88, -355.93, -274.94, -215.04, 
-114.3, -194.21, -103.73, -3.62, 104.2, 230.8, 154.25, 380.55, 
416.62, 260.07, 295.75, 295.75, 251.47, 220.67, 225.96, 180.72, 
121.52, 4.14, -127.23, -176.24, -332.11, -408.35, -494.11, -573.75, 
-582.62, -678.88, -730.38, -788.62, -831.94, -846.38, -895.95, 
-934.46, -968.15, -1033.12, -1097.62, -1150.08, -1157.3, -1254.04, 
-1340.2, -1441.75, -1500.47, -1550.52, -1605.39, -1681.44, -1709.84, 
-1715.22, -1672.34, -1607, -1522.59, -1440.57, -1421.18, -1345.62, 
-1247.95, -1190.77, -1181.58, -1071.65, -1037.62, -1010.39, -998.82, 
-986.57, -937.9, -887.29, -842.05, -831.46, -774.66, -703.91, 
-573.75, -533.59, -448.16, -446.8, -415.83, -394.43, -259.19, 
-104.69, -4.03, 58.59, -80.26, 52.11, -48.33, -142.23, -176.89, 
-233.68, -321.28, -416.57, -457.97, -458.93, -429.09, -422.35, 
-450.27, -431.98, -379.03, -260.63, -123.94, -2.65, 269.76, 455.55, 
548.92, 616.3, 691.38, 756.84, 888.72, 1016.97, 1157.18, 1198.02, 
1101.37, 1025.14, 929.84, 852.25, 766.48, 717.47, 733.81, 784.18, 
835.91, 1225.63, 1198.68, 925.3, 742.4, 814.13, 732.45, 586.79, 
394.84, 212.42, 28.64, -111.58, -337.56, -490.03, -526.07, -528.82, 
-547.2, -551.97, -552.3, -585.51, -551.34, -543.16, -526.1, -494.11, 
-466.88, -355.93, -274.94, -215.04, -114.3, -194.21, -103.73, 
-3.62, 104.2, 230.8, 154.25, 380.55, 416.62, 260.07, 295.75, 
295.75, 251.47, 220.67, 225.96, 180.72, 121.52, 4.14, -127.23, 
-176.24, -332.11, -408.35, -494.11, -573.75, -582.62, -678.88, 
-730.38, -788.62, -831.94, -846.38, -895.95, -934.46, -968.15, 
-1033.12, -1097.62, -1150.08, -1157.3, -1254.04, -1340.2, -1441.75, 
-1500.47, -1550.52, -1605.39, -1681.44, -1709.84, -1715.22, -1672.34, 
-1607, -1522.59, -1440.57, -1421.18, -1345.62, -1247.95, -1190.77, 
-1181.58, -1071.65, -1037.62, -1010.39, -998.82, -986.57, -937.9, 
-887.29, -842.05, -831.46, -774.66, -703.91, -573.75, -533.59, 
-448.16)), .Names = c("V1", "V2"), class = "data.frame", row.names = c(NA, 
-280L))

require(cobs)

xy <- dget(data.txt)

#Cumchord function (from Claude, 2008): Cumulative chordal distance vector
cumchord<-function(M)
{cumsum(sqrt(apply((M-rbind(M[1,],
M[-(dim(M)[1]),]))^2,1,sum)))}

z <- cumchord(xy)

#Calculating B-spline for x and y values separately
x <- cobs(z,xy[,1],nknots=50)
y <- cobs(z,xy[,2],nknots=50)

#Plot spline
plot(xy)
lines(x$fitted,y$fitted)

Image of resulting plot

I think this is an interesting question, but it would benefit from a [reproducible-example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). Can you share a little data, showing a shape that people can work with please. — user20650, Apr 07 '16 at 11:13
Just wondering, are you trying to reproduce/improve the example in Claude _Morphometrics with R_ (2008) pp. 208 ? — Vincent Bonhomme, Apr 07 '16 at 11:32
@Vincent: Cool thanks for the hint, that's exactly what I tried to, but way more elegant. In the working example above, I now followed Claude 2008, but changed it to use "cobs" instead of the "spline" function. Works perfectly with the "spline" function, but with cobs I get this wired double outline. But why? — tretelrusch, Apr 07 '16 at 19:15
@user20650: I shared my original dinosaur footprint data, the coordinates list is a bit long I hope this is ok! — tretelrusch, Apr 07 '16 at 19:16
your outline is the problem (seems that there are two of them). I have code that shows and makes all of this shorter, but it uses Momocs. If you want it, I turn it into an answer, cant fit here. By the way why not elliptical Fourier trasnforms? — Vincent Bonhomme, Apr 07 '16 at 19:26
@VincentBonhomme; as an interested observer, I would like to see your answer. — user20650, Apr 07 '16 at 19:38

score 5 · Accepted Answer · answered Apr 07 '16 at 20:14

Following the comments thread, here are some graphs. I use Momocs since I'm familiar with it and it will shorten examples.

I brief, your problem is that there are two outlines in your outline.

I also include original use of spline by Julien Claude, and two additional examples with bezier curves and elliptic Fourier transforms. All 4 can be used to describe an outline (and reconstruct it) and it's probably worth gathering them here.

The picture gathers the original shape and these 4 methods

Now the code. It's not particularly long but quite repetitive.

# devtools::install_github("vbonhomme/Momocs")
library(Momocs) # version 1.0.3

xy <- structure(list(
  V1 = c(124.9, 86.44, 97.22, 81.34, 49.09, 57.18, -77.6, -191.95, -284.67, -383.18, -379.27, -492.85, -547.72, -600.67, -713.29, -814.36, -868.27, -926.99, -958.76, -1025.18, -1077.16, -1105.07, -1126.25, -1112.77, -1087.74, -989.56, -911.59, -859.61, -745.06, -656.5, -682.01, -637.25, -601.71, -539.09, -394.79, -219.17, -170.17, -201.48, -122.52, -43.56, 127.97, 344.42, 539.09, 686.11, 987.63, 1253.31, 1283.15, 1536.32, 1741.14, 1832.35, 1700.3, 1787.43, 1911.31, 2017.49, 2097.81, 2135.93, 2093.73, 2066.96, 2063.78, 2022.94, 1978.69, 1919.44, 1904.03, 1895.37, 1854.22, 1810.23, 1771.09, 1741.48, 1642.45, 1553.96, 1472.96, 1396.04, 1141.65, 1085.82, 1055.02, 1358.24, 1325.94, 1031.91, 1287.14, 1265.36, 931.15, 872.12, 811.48, 755.65, 738.32, 697.41, 682.49, 647.35, 628.25, 620.09, 629.62, 675.22, 709.25, 718.78, 717.42, 551.09, 535.21, 540.98, 534.73, 546.76, 811.96, 823.03, 822.07, 607.4, 626.18, 637.73, 659.87, 756.13, 753.72, 735.91, 720.99, 676.71, 576.6, 508.26, 339.8, 179.53, 121.16, 45.6, 12.93, -9.87, -12.59, 16, 27.91, 37.78, 49.35, 8.51, 2.72, -1.02, 59.22, 58.2, 51.73, 54.45, 0.96, 10.59, 138.62, 149.69, 144.87, 142.26, 146.34, 125.24, 124.9, 86.44, 97.22, 81.34, 49.09, 57.18, -77.6, -191.95, -284.67, -383.18, -379.27, -492.85, -547.72, -600.67, -713.29, -814.36, -868.27, -926.99, -958.76, -1025.18, -1077.16, -1105.07, -1126.25, -1112.77, -1087.74, -989.56, -911.59, -859.61, -745.06, -656.5, -682.01, -637.25, -601.71, -539.09, -394.79, -219.17, -170.17, -201.48, -122.52, -43.56, 127.97, 344.42, 539.09, 686.11, 987.63, 1253.31, 1283.15, 1536.32, 1741.14, 1832.35, 1700.3, 1787.43, 1911.31, 2017.49, 2097.81, 2135.93, 2093.73, 2066.96, 2063.78, 2022.94, 1978.69, 1919.44, 1904.03, 1895.37, 1854.22, 1810.23, 1771.09, 1741.48, 1642.45, 1553.96, 1472.96, 1396.04, 1141.65, 1085.82, 1055.02, 1358.24, 1325.94, 1031.91, 1287.14, 1265.36, 931.15, 872.12, 811.48, 755.65, 738.32, 697.41, 682.49, 647.35, 628.25, 620.09, 629.62, 675.22, 709.25, 718.78, 717.42, 551.09, 535.21, 540.98, 534.73, 546.76, 811.96, 823.03, 822.07, 607.4, 626.18, 637.73, 659.87, 756.13, 753.72, 735.91, 720.99, 676.71, 576.6, 508.26, 339.8, 179.53, 121.16, 45.6, 12.93, -9.87, -12.59, 16, 27.91, 37.78, 49.35, 8.51, 2.72, -1.02, 59.22, 58.2, 51.73, 54.45, 0.96, 10.59, 138.62, 149.69, 144.87, 142.26, 146.34, 125.24),
  V2 = c(-446.8, -415.83, -394.43, -259.19, -104.69, -4.03, 58.59, -80.26, 52.11, -48.33, -142.23, -176.89, -233.68, -321.28, -416.57, -457.97, -458.93, -429.09, -422.35, -450.27, -431.98, -379.03, -260.63, -123.94, -2.65, 269.76, 455.55, 548.92, 616.3, 691.38, 756.84, 888.72, 1016.97, 1157.18, 1198.02, 1101.37, 1025.14, 929.84, 852.25, 766.48, 717.47, 733.81, 784.18, 835.91, 1225.63, 1198.68, 925.3, 742.4, 814.13, 732.45, 586.79, 394.84, 212.42, 28.64, -111.58, -337.56, -490.03, -526.07, -528.82, -547.2, -551.97, -552.3, -585.51, -551.34, -543.16, -526.1, -494.11, -466.88, -355.93, -274.94, -215.04, -114.3, -194.21, -103.73, -3.62, 104.2, 230.8, 154.25, 380.55, 416.62, 260.07, 295.75, 295.75, 251.47, 220.67, 225.96, 180.72, 121.52, 4.14, -127.23, -176.24, -332.11, -408.35, -494.11, -573.75, -582.62, -678.88, -730.38, -788.62, -831.94, -846.38, -895.95, -934.46, -968.15, -1033.12, -1097.62, -1150.08, -1157.3, -1254.04, -1340.2, -1441.75, -1500.47, -1550.52, -1605.39, -1681.44, -1709.84, -1715.22, -1672.34, -1607, -1522.59, -1440.57, -1421.18, -1345.62, -1247.95, -1190.77, -1181.58, -1071.65, -1037.62, -1010.39, -998.82, -986.57, -937.9, -887.29, -842.05, -831.46, -774.66, -703.91, -573.75, -533.59, -448.16, -446.8, -415.83, -394.43, -259.19, -104.69, -4.03, 58.59, -80.26, 52.11, -48.33, -142.23, -176.89, -233.68, -321.28, -416.57, -457.97, -458.93, -429.09, -422.35, -450.27, -431.98, -379.03, -260.63, -123.94, -2.65, 269.76, 455.55, 548.92, 616.3, 691.38, 756.84, 888.72, 1016.97, 1157.18, 1198.02, 1101.37, 1025.14, 929.84, 852.25, 766.48, 717.47, 733.81, 784.18, 835.91, 1225.63, 1198.68, 925.3, 742.4, 814.13, 732.45, 586.79, 394.84, 212.42, 28.64, -111.58, -337.56, -490.03, -526.07, -528.82, -547.2, -551.97, -552.3, -585.51, -551.34, -543.16, -526.1, -494.11, -466.88, -355.93, -274.94, -215.04, -114.3, -194.21, -103.73, -3.62, 104.2, 230.8, 154.25, 380.55, 416.62, 260.07, 295.75, 295.75, 251.47, 220.67, 225.96, 180.72, 121.52, 4.14, -127.23, -176.24, -332.11, -408.35, -494.11, -573.75, -582.62, -678.88, -730.38, -788.62, -831.94, -846.38, -895.95, -934.46, -968.15, -1033.12, -1097.62, -1150.08, -1157.3, -1254.04, -1340.2, -1441.75, -1500.47, -1550.52, -1605.39, -1681.44, -1709.84, -1715.22, -1672.34, -1607, -1522.59, -1440.57, -1421.18, -1345.62, -1247.95, -1190.77, -1181.58, -1071.65, -1037.62, -1010.39, -998.82, -986.57, -937.9, -887.29, -842.05, -831.46, -774.66, -703.91, -573.75, -533.59, -448.16)), .Names = c("V1", "V2"), class = "data.frame", row.names = c(NA, -280L))


### First thing first: double outline ---------------------
coo_plot(xy)
ldk_labels(xy) # blurry since superimposed
coo_plot(xy[1:140, ], lwd=3) # first shape
coo_draw(xy[-(1:140), ], border="white") # second shape

# so from here, we will use the first 140th points from xy and name it 'shp' to avoid confusion
# if you're bored with dinos footprints, you can use beer bottles with shp <- bot[9] for a guinness

shp <- xy[1:140, ]

### 1.Natural splines ---------------------

shp_cumchord <- coo_perimcum(shp) # cumchord equivalent
shp_spline <- cbind(spline(shp_cumchord, shp[, 1], method="natural", n=120)$y,
                    spline(shp_cumchord, shp[, 2], method="natural", n=120)$y)

coo_plot(shp, main="natural spline", zoom=1.2)
coo_draw(shp_spline, border="blue", lwd=2)

### 2. B-splines with cobs ---------------------
library(cobs)

shp_bspline <- cbind(cobs(shp_cumchord, shp[, 1], nknots=50)$fitted,
                     cobs(shp_cumchord, shp[, 2], nknots=50)$fitted)

coo_plot(shp, main = "bspline", zoom=1.2)
coo_draw(shp_bspline, border="blue")

### 3. Bezier curves  ---------------------
# built in function so it's shorter
shp_bezier <- shp %>% bezier() %>% bezier_i()

coo_plot(shp, main = "bezier", zoom=1.2)
coo_draw(shp_bezier, border="blue")


### 4. elliptic Fourier transforms ---------------------
# another built in function

shp_eft <- shp %>% efourier() %>% efourier_i()

coo_plot(shp, main = "bspline", zoom=1.2)
coo_draw(shp_eft, border="blue")

### 5. A panel of original shape and 4 methods ---------
Out(list(original=shp,
         nat_spline=shp_spline, bspline=shp_bspline,
         bezier=shp_bezier, eft=shp_eft)) %>% 
panel(names=TRUE, dim=c(1, 5))

also, If you think splines should be included in Momocs I would be very happy to discuss how and why. — Vincent Bonhomme, Apr 07 '16 at 20:15
Awesome, thank you for the comprehensive overview, and for finding the error with the double outline! Elliptical Fourier with Momocs was in fact the first I tried, and worked quite well. The reason I wanted to try cobs is its ability to constrain the spline (i.e., define points which the spline has to pass). Is this possible with elliptical Fourier transforms as well? — tretelrusch, Apr 07 '16 at 21:36
It's possible to align them using landmarks, not particularly to constraint them. But the outlines will be very very close (closer than the pixel sometime, so technically passing by these points) anyway after, say 7-8 harmonics. — Vincent Bonhomme, Apr 08 '16 at 06:27

Approximate a shape outline using constrained B-splines

1 Answers1

Linked