How to straighten a curve?

Question

What's the easiest way to get the following animation for a given parabola where the arc length is preserved (that is, for two points in the parabola, I straighten it as follows)

My main goal is the following, given a set of random points and a set of points following a parabola, I want to "straighten" the parabola and all the surrounding points according to the straightening of that parabola. Consider the following schematic

For instance, we could consider the points as a starter

f[x_, a_, b_] := a x^2 + b;
pts1 = RandomReal[{-1, 1}, {150, 2}];
ptsf = RandomReal[{-1, 1}, 50];
pts2 = Map[{#, f[#, .5, -.5]} &, ptsf];
Graphics[{Gray, Point /@ pts1, Red, Point /@ pts2}]

Where the red dots follow a parabola. How do change pts1 based on a and b? From such a pattern, I'd expect, after the transformation, to get something like

where the dots are sparser above the line and denser below it. I tried to make it work by using preservation of the arc length somehow, but it's too slow. For more details on the math, here's a more general question.

Edit: I'm pretty satisfied with the answers provided and they do answer part of the question. As for the plane transformation, this seems to be trickier. Consider the following (very messy) image.

What I want, as suggested by @I.M., is to transform the plane such that we get two focusing and defocusing regions. The tricky part is how to "straighten" the surrounding points. If we track three line segments with end points $\{1,2,3,4,5,6\}$, the line segment will change size as the curve changes, and I suspect the new points will follow new parabolas (see thiner coloured lines). Not sure if I can prove this. Ultimately, what I'm asking is how to implement the following plane transformation, where the x-axis becomes a parabola

What do you think?

Answer 1

Here is a "homotopy style" transformation between functions $h(x, q) = (1-q) f(x) + q g(x)$.

(* set initial data points *)
ClearAll[x] ;
x = Subdivide[-1, 1, 14] ;

(* replace x^2 with q/2 + (1-q) x^2, i.e q = 0 -- initial function, q = 1 -- final function (line) *)
ClearAll[fun];
fun[q_][x_] := 0.5*q + (1.0 - q)*x^2 ;

(* arc length *)
ClearAll[arc];
arc[q_][x_] := 1/4 (2 x Sqrt[1 + 4 (-1 + q)^2 x^2] + ArcSinh[2 (-1 + q) x]*1/(-1 + q)) ;

(* given x and q,find new x with same arc length *)
ClearAll[fuc];
fuc[q_][x_] := Catch[
    Block[
        {y},
        If[q == 1.0, Throw[{arc[0.0][x], fun[q][x]}]] ;
        y = y /. FindRoot[arc[0.0][x] == arc[q][y], {y, x}] ;
        {y, fun[q][y]}
    ]
] ;

(* test arc conservation *)
{arc[0.0][1.0], arc[0.25][First[fuc[0.25][1.0]]]}

(* plot *)
Manipulate[
    ListLinePlot[
        Map[fuc[q], x],
        Mesh -> All, AspectRatio -> 1, PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}, Background -> Black, PlotStyle -> Cyan, Axes -> False,
        PlotLabel -> Style[StringTemplate["q = ``"][q], White]
    ],
    {q, 0.0, 1.0, 0.01}
]
(* {1.4789428575445975`,1.478942857544597`} *)

Answer 2

The arclength is given by this and is the big expression containing ArcSinh later on:

Integrate[Sqrt[1 + (2 a x)^2], {x, -k, k}, 
 Assumptions -> a > 0 && k > 0]

For a given a and len (arc length) find the best k (integration range) that gives that length by doing an NMinimize (for whatever reason, NSolve didn't work that well so I minimized a square error). The error in the NMinimize is extremely small.

Manipulate[
 {err, sol} = 
  NMinimize[(k Sqrt[1 + 4 a^2 k^2] + ArcSinh[2 a k]/(2 a) - len)^2, k];
 k1 = k /. sol;
 Plot[a t^2, {t, -k1, k1}, Background -> Black, 
  PlotRange -> {{-1, 1}, {-0.2, 1}}, AspectRatio -> 1, Axes -> False]
 , {a, 0.001, 4}, {len, .5, 3}]

Answer 3

f[t_] = t^2;
{l, u} = {-1, 1};

To transition to a line, curves that preserve the arc length and scales the signed curvature by q are used. First such curves are made for a variety of q values:

sols = <|# -> First[{x, y} /. NDSolve[##2]] & @@@ Table[Evaluate[{1 - q,
  {D[x'[t]^2 + y'[t]^2 == Total[D[{t, f[t]}, t]^2], t],
  (x'[t] y''[t] - y'[t] x''[t])/(x'[t]^2 + y'[t]^2)^(3/2) == (q f''[t])/(1 + f'[t]^2)^(3/2),
  x[l] == l, y[l] == f[l], x'[l] == 1, y'[l] == f'[l]}, {x, y}, {t, l, u}}], {q, Subdivide[25]}]|>;

The curves are heuristically adjusted to be similar to the initial curve:

rigids = Last[FindGeometricTransform[#, #2, TransformationClass -> "Rigid"] & @@
      Transpose[Table[{{t, f[t]}, Through[#[t]]}, {t, Subdivide[l, u, 15]}]]] & /@ sols;

Linear interpolation between the solved q values is used to allow an arbitrary q:

getKeys = Nearest[Keys[sols]];
straightening[q_, t_] := Interpolation[{#, rigids[#][Through[sols[#][t]]]} & /@ getKeys[q, 2], q, InterpolationOrder -> 1]

Manipulate[ParametricPlot[straightening[q, t], {t, -1, 1},
  PlotRange -> {{-2, 2}, {-2, 2}}, AspectRatio -> Automatic], {q, 0, 1}]

Answer 4

Clear[a, x, y, x1, x2]
SeedRandom[2]
npts = 200;
r = 3;
pts1 = RandomReal[{-1.5, 1.5}, {npts, 2}];
parms = {a -> 1/2, b -> -1/2};
f[x_, y_] := y - a x^2 - b /. parms
px[x_] := 1/2 x Sqrt[1 + 4 a^2 x^2] + ArcSinh[2 a x]/(4 a) /. parms
vecn = Grad[f[x, y], {x, y}];
GRS = {};
INTS = {};

For[k = 1, k <= npts, k++,
  {x0, y0} = pts1[[k]];
  equs = Join[Thread[{x, y} == {x0, y0} + lambda vecn], {f[x, y] == 0}] /. parms;
  sols = Quiet@Solve[equs, {x, y, lambda}, Reals];
  {x1, y1, l1} = {x, y, lambda} /. sols[[1]];
  AppendTo[INTS, {px[x1], Sign[f[x0, y0]] Sqrt[(x1 - x0)^2 + (y1 - y0)^2]}];
  vn1 = vecn /. sols[[1]] /. parms;
  AppendTo[GRS, ParametricPlot[Thread[{x0, y0} + l vn1], {l, 0, l1}, PlotStyle -> {Thin, Black}]]
  ];

gr1 = ContourPlot[(f[x, y] /. parms) == 0, {x, -r, r}, {y, -r, r}, ContourStyle -> Red];
gr2 = Graphics[{Black, PointSize[0.01], Point /@ pts1}];
xk = Transpose[INTS][[1]];
yk = Transpose[INTS][[2]] + b /. parms;
INTS0 = Join[{xk}, {yk}] // Transpose;
gr3 = ListPlot[INTS0, PlotStyle -> {Red, PointSize[0.01]}, Filling -> (b /. parms)];
gr4 = Plot[(b /. parms), {x, -r, r}, PlotStyle -> Red];
Show[gr1, gr2, GRS, PlotRange -> {{-r, r}, {-r, r}}, AspectRatio -> 1]
Show[gr4, gr3, PlotRange -> {{-r, r}, {-r, r}}, AspectRatio -> 1]
Show[gr1, gr2, gr3, gr4, GRS, PlotRange -> {{-r, r}, {-r, r}}, AspectRatio -> 1]

EDIT

Included a Manipulateversion. Run at low speed.

Clear["Global`*"]
SeedRandom[5]
npts = 10; r = 3; a = 2; b = -1/2;
pts1 = RandomReal[{-2, 2}, {npts, 2}];
L = (x0 - x)^2 + (y0 - y)^2 + lambda (y - a x^2 - b);
grad = Grad[L, {x, y, lambda}];
sollambda1 = Solve[grad[[1]] == 0, lambda][[1]];
sollambda2 = Solve[grad[[2]] == 0, lambda][[1]];
equ = (lambda /. sollambda1) == (lambda /. sollambda2);
solsdist = Solve[{equ, y - a x^2 - b == 0}, {x, y}];
f[x_, y_, a_] := y - a x^2 - b
px[x_, a_] := 1/2 x Sqrt[1 + 4 a^2 x^2] + ArcSinh[2 a x]/(4 a)
vecn[x_, a_] := {-((2 a x)/Sqrt[1 + 4 (a x)^2]), 1/Sqrt[1 + 4 (a x)^2]}

nearpt[x0_, y0_] := Module[{sols, ptsp, dists, p0, pnear, near}, sols = solsdist // N // Chop;
  ptsp = {x, y} /. sols;
  p0 = {x0, y0};
  dists = Table[{Norm[ptsp[[k]] - p0], k}, {k, 1, 3}];
  pnear = Sort[dists][[1]];
  near = ptsp[[pnear[[2]]]];
  {x1, y1} = near;
  Return[{x1, y1, Norm[near - p0] Sign[f[p0[[1]], p0[[2]], a]]}]
  ]

Manipulate[
 GRS1 = {};
 INTS1 = {};
 For[k = 1, k <= npts, k++,
  {x0, y0} = pts1[[k]];
  {x1, y1, l1} = nearpt[x0, y0];
  dist = px[x1, a];
  x2 = x /. Quiet@FindRoot[px[x, a1] - dist == 0, {x, x1}];
  y2 = y /. Quiet@Solve[f[x2, y, a1] == 0, y][[1]];
  vn2 = vecn[x2, a1];
  AppendTo[GRS1, ParametricPlot[{x2, y2} - l vn2, {l, 0, l1}, PlotStyle -> {Thin, Black}]];
  AppendTo[INTS1, {x2, y2} - l1 vn2]
  ];
 gr11 = ContourPlot[f[x, y, a1] == 0, {x, -r, r}, {y, -r, r}, PlotPoints -> 30, ContourStyle -> Red];
 Show[gr11, ListPlot[INTS1, PlotStyle -> {Black, PointSize[0.02]}], GRS1, PlotRange -> {{-r, r}, {-r, r}}, AspectRatio -> 1],
 {{a1, 0.01}, 0.01, a, 0.01}
 ]