First you need to order your edge location (xy) coordinates along the arc of the curve; this can be achieved with the image processing toolbox in-built "bwboundaries", so you can do that before you extract the edge (because the edge image is not going to be all that helpful to get your ordered coordinates along perimeter).
In principle, you can then calculate the arc coordinate for each xy using pythagoras on each distance between each consecutive pair: 
% generate a curve
x = transpose(linspace(0,2,7));
y = x.^2;
% pack into a coordinates array
xy = [x,y];
% compute [pseudo] arc-length
d = diff(xy); % nearest neighbor component diff
r = sqrt(sum(d.^2,2)); % nearest neighbor distance
s = cumsum([0;r]); % arc length
You can then think about several ways to achieve the derivatives: finite differences, splines, etc., to implement your formula above.
Here's a simple way with 1st order finite differences using diff(), but you'd have to modify to deal with edges (you keep shortening the coordinate lists by 1 each application of diff)
dxds = diff(xy(:,1))./diff(s); % this also equals cos(phi) where phi is inclination angle
dyds = diff(xy(:,2))./diff(s);
phi = unwrap(atan2(dyds,dxds);
curv = diff(phi)./diff(s(2:end)); % curvature is just the derivative of inclination angle
ddxdss = diff(dxds)./diff(s(2:end));
% etc
You've already linked to another answer involving spline interpolation and differentiating on the spline.
You could imagine other ways to take derivatives, like higher order finite differences (can do on unevenly sampled arc length) or convolutions (need evenly sampeld arc length).
In practice, the problem is that coordinates you get from images will be aliased so derivatives of your curve, parameterized or not, are going to be really bad - useless if you are only using a few neighboring pixels.
So your real problem is actually about how to pre-process your pixel coordinates and/or how to sample your coordinates to avoid aliasing problems.
