Repairing a panorama

Question

I was on a cruise and the panorama I took was warped due to the ship rocking:

I'm hoping to correct this so that the horizon is level. Here's the full resolution image:

pano = CloudImport @ CloudObject[
 "https://www.wolframcloud.com/objects/505ab4f5-2f3b-4284-b159-23b650ec45e6"];
{w, h} = ImageDimensions@pano

Here's what I've tried so far: basically making an Interpolation function for ImageTransformation to use:

(* manually chosen points *)
horizon = {{10824.347571942448`, 1828.3412145283764`}, {8600.42723321343`, 1926.247726818544`}, {6562.498576139089`, 2068.3502822741793`}, {4752.2000024980025`, 2189.78709532374`}, {2206.937599920064`, 2306.927070843324`}};
stable = {{16036.104578836932`, 117.17682104316555`}, {16057.281849520385`, 3465.1293964828137`}, {8767.18545163869`, 3474.4392111310954`}, {8663.24290567546`,  173.2403202438045`}, {349.88534172661866`, 160.86133593125487`}, {343.2354741207032`, 3377.453449740208`}};
fixed = Thread @ {horizon[[All, 1]], horizon[[All, 2]] // Mean};

HighlightImage[pano, {Red, Point[horizon], Blue, Point[stable], 
                      Green, Line[Thread[{horizon, fixed}]]}]

t = Thread[{horizon, fixed}]; s = Thread[{stable, stable}];
inter = Interpolation[Join[t, s], InterpolationOrder -> All]
f[{x_: 0, y_: 0}] :=  With[{v = inter[x, y]}, {Clip[v[[1]], {0, w}], Clip[v[[2]], {0, h}]}]

So it looks like it should be correcting it:

Finally, then this returns an image with a flat horizon, but it looks a bit too pinched:

ImageForwardTransformation[pano, f, DataRange -> Full]

Update:

Reference links and more examples:

• https://m.youtube.com/watch?v=QZFfJgRVqrE
• https://fstoppers.com/education/quick-fix-correcting-distortion-stitched-panoramas-190686
• https://graphicdesign.stackexchange.com/questions/98005/how-to-undistort-this-panorama-image-in-adobe-photoshop-elements
• https://www.youtube.com/watch?v=QZFfJgRVqrE

Answer 1

Here's an approach I took, but it's slightly manual and far from perfect.

First I worked with a smaller image and then applied the function to the original dataset.

Original data:

pano = CloudImport @ CloudObject[
 "https://www.wolframcloud.com/objects/505ab4f5-2f3b-4284-b159-23b650ec45e6"];

Smaller dataset to work on:

panosmall = ImageResize[pano, Scaled[1/8]];
{dimx, dimy} = ImageDimensions[panosmall];

I manually set the range of the horizon, edge detected, and extracted:

{xmin, xmax} = {290, -600};
edges = EdgeDetect[ImageTake[panosmall, All, {xmin, xmax}]];

cands = ImageSubtract[edges, DeleteBorderComponents[edges]];
horizon = SelectComponents[cands, "Elongation", -1];

Now we find a smooth curve to represent the horizon:

hline = Sort[PixelValuePositions[horizon, 1]][[All, 2]];

cpts = Join[
 ConstantArray[First[hline], xmin - 1], 
 hline, 
 hline[[-1]] - Range[1, -xmax]/18.5
][[1 ;; -1 ;; 1]];

Here I manually extended the x-range into the areas I cropped out. On the left I extrapolated with a constant value and on the right a linear one:

ListLinePlot[cpts, GridLines -> {{xmin, ImageDimensions[panosmall][[1]] + xmax}, None}]

A smoother version:

bf = BSplineFunction[cpts, SplineDegree -> Length[cpts] - 1];

Plot[bf[t], {t, 0, 1}, GridLines -> {{xmin/dimx, 1 + xmax/dimx}, None}]

The approach now is to shift each column with this curve. We can speed things up by 'caching' floats with a NearestFunction.

With[{fac = 1.0/dimx, center = {.9, .1}.{bf[0], bf[1]}},
  shifts = fac*(bf /@ Range[.5/dimx, 1, 1./dimx] - center);
]

nf = Nearest[Range[.5/dimx, 1, 1./dimx] -> shifts];

imfunc[{x_, y_}] := {x, y + nf[x][[1]]}

ImageTransformation[panosmall, imfunc]

And 10 min later we can get the full res version:

dimxl = ImageDimensions[pano][[1]];

With[{fac = 1.0/dimx, center = {.9, .1}.{bf[0], bf[1]}},
  shiftslarge = fac*(bf /@ Range[.5/dimxl, 1, 1./dimxl] - center);
]

nflarge = Nearest[Range[.5/dimxl, 1, 1./dimxl] -> shiftslarge];

imfunclarge[{x_, y_}] := {x, y + nflarge[x][[1]]}

ImageTransformation[pano, imfunclarge]