# Fit image of mountain to gaussian

I want to detect the edges of an image of a mountain and then fit that to a gaussian.

I am unsure about how to only get the mountain edge, rather than every single edge in the image.

How do I detect the outline of the mountain and then fit it with a gaussian?

I think I need to

1. Import Image
2. Detect Edges
3. Binarize
4. Fit Gaussian

However I am unable to detect only the edge I want.

If I use EdgeDetect

Then I get this image

• Looks a bit like Chimborazo from one vantage point. What mountain is it? Commented Mar 14, 2016 at 9:44
• @MikeHoneychurch It IS the Chimborazo. Nice mountain:) Commented Mar 14, 2016 at 9:54
• @Dr.belisarius went mountain biking down from the refuge hut to a nearby town back in 2013. Ecuador is a nice place Commented Mar 14, 2016 at 10:20
• Yes this is Chimborazo. I grew up in Ecuador and climbed it when I was 16 :) Commented Mar 14, 2016 at 17:20

A method not needing "magic numbers":

i     = Import@"https://i.sstatic.net/1Ui83.jpg";
id    = ImageDimensions@i;
mask  = {⌊#/2⌋, ⌊#/4⌋} &@ Reverse@id;

model = a E^(-b (x - x0)^2) + c;
fit   = FindFit[p, model, {{a, Max[Last /@ p]}, b, c, {x0, First@id/2}}, x];
j     = Plot[model /. fit, {x, 0, First@id}, Axes -> False, PlotStyle -> Thick];

Show[i, j, ImageSize -> 400]


A careful method of finding the mountain boundary including what is masked by the snow but not including the clouds:

i = Import@"https://i.sstatic.net/1Ui83.jpg";
edge = Closing[
4000], {{1, 1, 1, 1, 1, 1}}]


(UPDATE: If one considers that at the left bottom corner we have not a snow but clouds, one can completely remove them beforehand with the following code:

dims = ImageDimensions[i];
mask = Graphics[{White, Rectangle[{0., 364}, {70, 431}]}, Background -> Black,
ImageSize -> dims, PlotRange -> Transpose[{{0, 0}, dims}]];
ImageApply[If[#[[1]] > .7, {.6, .8, .9}, #] &, i, Masking -> mask]


).

Extract the upper pixels and compare with the original image:

data = SortBy[Max /@ Transpose[#] & /@ GatherBy[ImageValuePositions[edge, White], First],
First];
Show[i, Graphics[{Red, Thick, Line[data]}]]


Sequentially fit the data by sums of Gaussians:

gauss[i_] := s[i]*PDF[NormalDistribution[a[i], b[i]], x]
model[n_] := base + Sum[gauss[i], {i, n}];
style[n_] :=
Sequence[Evaluated -> True,
PlotStyle -> {{Red, Thick, Dashing[{}]},
Sequence @@
Take[{{Green, Dashed}, {Blue, Dashed}, {Yellow, Dashed}, {Magenta, Dashed}},
n], {Red, Thick, Dashed}}];
toPlot[n_] := Join[{model[n]}, Table[base + gauss[i], {i, n}], {base}] /. fit[n]
fit[1] = FindFit[data,
base + gauss[1], {{base, 300}, {a[1], 1000}, {b[1], 100}, {s[1], 2*10^5}}, x]
Show[i, Plot[{model[n] /. fit[1], base /. fit[1]}, {x, 0, 1900},
PlotStyle -> {{Red, Thick, Dashing[{}]}, {Red, Thick, Dashed}}], ImageSize -> 600]

{base -> 426.193, a[1] -> 983.737, b[1] -> 365.883, s[1] -> 359741.}


fit[2] = FindFit[
data, {base + gauss[1] + gauss[2], base >= 0 && And @@ Table[b[i] > 0, {i, 2}]},
List @@@ Join[fit[1], Rest[fit[1]] /. (p_[1] -> v_) :> {p[2], .95 v}], x]
Show[i, Plot[toPlot[2], {x, 0, 1900}, Evaluate[style[2]]], ImageSize -> 600]

{base -> 3.98312*10^-7, a[1] -> 1061.57, b[1] -> 1171.57, s[1] -> 1.7309*10^6,
a[2] -> 946.356, b[2] -> 265.231, s[2] -> 161146.}


fit[3] = FindFit[
data, {base + gauss[1] + gauss[2] + gauss[3],
base >= 0 && And @@ Table[b[i] > 0, {i, 3}]},
List @@@ Join[fit[2], Rest[fit[1]] /. (p_[1] -> v_) :> {p[3], v}], x]
Show[i, Plot[toPlot[3], {x, 0, 1900}, Evaluate[style[3]]], ImageSize -> 600]

{base -> 320.683, a[1] -> 1792.98, b[1] -> 267.814, s[1] -> 87111.7, a[2] -> 961.519,
b[2] -> 371.918, s[2] -> 471238., a[3] -> 233.468, b[3] -> 190.528, s[3] -> 37203.8}


fit[4] = FindFit[
data, {base + gauss[1] + gauss[2] + gauss[3] + gauss[4],
base >= 0 && And @@ Table[b[i] > 0, {i, 4}]},
List @@@ Join[fit[3], Rest[fit[1]] /. (p_[1] -> v_) :> {p[4], v}], x]
Show[i, Plot[toPlot[4], {x, 0, 1900}, Evaluate[style[4]]], ImageSize -> 600]

{base -> 57.7774, a[1] -> 1865.96, b[1] -> 342.875, s[1] -> 289761., a[2] -> 4.27395,
b[2] -> 883.184, s[2] -> 626611., a[3] -> 254.158, b[3] -> 159.113, s[3] -> 22215.1,
a[4] -> 992.112, b[4] -> 385.078, s[4] -> 582369.}


• data = HitMissTransform[edge, Join[ConstantArray[{-1}, 20], {{1}}]] Commented Mar 14, 2016 at 9:51
• @belisarius Nice, but the edge is shifted: HighlightImage[i, HitMissTransform[edge, Join[ConstantArray[{-1}, 20], {{1}}]]]. Commented Mar 14, 2016 at 10:00
• Oh, yes,sorry. It is 20 pixels up due to the ConstantArray. Let me see if there is an easy fix :) Commented Mar 14, 2016 at 10:08
• ImagePad[HitMissTransform[edge, Join[ConstantArray[{-1}, 20], {{1}}]], 10 {{0, 0}, {-1, 1}}] but not so terse Commented Mar 14, 2016 at 10:13

A rough draft of a solution:

img = Import["https://i.sstatic.net/1Ui83.jpg"];
data = ImageValuePositions[
EdgeDetect@
DeleteSmallComponents[
MorphologicalBinarize[ImageTake[img, -840], {0.595, 0.9958}]],
White];
nlm = NonlinearModelFit[data,
c + s*PDF[NormalDistribution[a, b], x], {{a, 1000}, {b, 100}, {c, 300}, {s, 2*10^5}}, x];
Show[ListPlot[data, PlotStyle -> Black],
Plot[nlm[x], {x, 0, 1900}, PlotStyle -> Red], Frame -> True]


A image-processing difference from other's

pic = ImageAdjust@Import["https://i.sstatic.net/1Ui83.jpg"];
mountain =
SelectComponents[
Binarize[pic, 0.71] // DeleteSmallComponents // ColorNegate,
"Count", -1]


Get the point of the top in every column and hightlight it.

point = First@*MaximalBy[Last] /@
GatherBy[ImageValuePositions[mountain, 1], First];
HighlightImage[pic, point]


The last work have some boring and you can refer to other's.

Using a modification of Karsten 7.'s code to fit a sum of Gaussians I got this approximation:

Here is the code:

img = Import["https://i.sstatic.net/1Ui83.jpg"];
data = ImageValuePositions[
EdgeDetect@
DeleteSmallComponents[
MorphologicalBinarize[ImageTake[img, -840], {0.595, 0.9958}]],
White];
gaussianSum =
c + Total[Flatten@Table[ s[m, sd]*PDF[NormalDistribution[a + u*m,
b + v*sd], x], {m, -2, 2, 1}, {sd, -2, 2, 1}]];
svars = Flatten@Table[s[m, sd], {m, -2, 2, 1}, {sd, -2, 2, 1}];
nlm = NonlinearModelFit[data,
gaussianSum,
{{a, 1000}, {b, 100}, {c, 300},
Sequence @@ Thread[{svars, 5*10^4}], {u, 20}, {v, 5}}, x];
Show[ListPlot[data, PlotStyle -> Black],
Plot[nlm[x], {x, 0, 1900}, PlotStyle -> Red], Frame -> True]


Here is the result function: