[Question]
  [
*This is the first in a series of Island Golf challenges. [Next challenge](https://codegolf.stackexchange.com/q/114053/16766)*
Given an island in ASCII-art, output an optimal path to circumnavigate it.
## Input
Your input will be a rectangular grid consisting of two characters, representing land and water. In the examples below, land is `#` and water is `.`, but you may substitute any two distinct characters you wish.
```
...........
...##......
..#####....
..#######..
.#########.
...#######.
...#####.#.
....####...
...........
```
There will always be at least one land tile. The land tiles will all be contiguous (i.e. there's only one island). The water tiles will also be contiguous (i.e. there are no lakes). The outer border of the grid will all be water tiles. Land tiles will **not** be connected diagonally: i.e., you will never see something like
```
....
.#..
..#.
....
```
## Output
Your code must output the same grid, with a *shortest circumnavigation* drawn on it. In the examples below, the circumnavigation path is drawn with `o`, but you may substitute any character as long as it is distinct from your land and water characters.
A *circumnavigation* is a simple closed curve, drawn entirely on water tiles, that fully encircles all land tiles on the grid. Diagonal connections **are** allowed. For instance, this is a circumnavigation of the above island (but not a shortest one):
```
.ooooo.....
o..##.oo...
o.#####.o..
o.#######o.
o#########o
ooo#######o
..o#####.#o
..oo####..o
....oooooo.
```
The length of a circumnavigation is computed as follows: For every pair of adjacent tiles on the path, if they are connected horizontally or vertically, add 1; if they are connected diagonally, add √2. The length of the above path is 22 + 7√2 (≈ 31.9).
A *shortest* circumnavigation is a circumnavigation with the shortest possible length. Your program should output any one path that satisfies this condition. For most islands, there will be multiple possible solutions. Here is one solution for the above island, with length 10 + 13√2 (≈ 28.4):
```
...oo......
..o##oo....
.o#####oo..
.o#######o.
o#########o
.o.#######o
..o#####.#o
...o####.o.
....ooooo..
```
## Details
Your solution may be a [full program or a function](https://codegolf.meta.stackexchange.com/a/2422/16766). Any of the [default input and output methods](https://codegolf.meta.stackexchange.com/q/2447/16766) are acceptable.
Your input and output may be a multiline string or a list of strings. **If** your language has a character type distinct from single-character strings, you may substitute "list of characters" for "string" in the previous sentence. If your language needs to input the height and/or width of the grid, you may do so. Your output may (optionally) have a single trailing newline. As mentioned above, you may use any three distinct characters in place of `#.o` (please specify in your submission which characters you're using).
## Test cases
**A.** Islands with unique shortest circumnavigations:
```
...
.#.
...
.o.
o#o
.o.
......
.####.
......
.oooo.
o####o
.oooo.
......
......
..##..
...#..
......
......
......
..oo..
.o##o.
..o#o.
...o..
......
.......
.#####.
...#...
...#...
.#####.
.......
.ooooo.
o#####o
o..#..o
o..#..o
o#####o
.ooooo.
.......
...#...
...#...
.#####.
...#...
...#...
.......
...o...
..o#o..
.o.#.o.
o#####o
.o.#.o.
..o#o..
...o...
.......
.#####.
.##..#.
..#..#.
.......
.ooooo.
o#####o
o##..#o
.o#..#o
..oooo.
```
**B.** Example of an island with multiple possible solutions:
```
........
....##..
...####.
..###...
.#####..
.#####..
..##....
........
```
Possible outputs:
```
....oo..
...o##o.
..o####o
.o###.o.
o#####o.
o#####o.
.o##oo..
..oo....
....oo..
...o##o.
..o####o
.o###.o.
o#####o.
o#####o.
.o##.o..
..ooo...
....oo..
...o##o.
..o####o
.o###..o
o#####.o
o#####o.
.o##oo..
..oo....
....oo..
...o##o.
..o####o
.o###..o
o#####.o
o#####o.
.o##.o..
..ooo...
```
**C.** [Large test case as a Gist](https://gist.github.com/dloscutoff/6bc4a16c1ee9b7932aef9c00ea99cd3e)
---
This is [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'"): the shortest code in each language wins.
  
[Answer]
## Mathematica (version 9), 165 bytes
The nice, short `ConvexHullMesh` function that Greg Martin used was only introduced in Mathematica version 10, so I thought I'd make an attempt without it, using my ancient Mathematica version 9. I managed to get it slightly shorter, though! It's a function that takes and returns a list of strings (with `.`, `#` and `o` as the symbols).
```
""<>#&/@("o"MorphologicalTransform[MorphologicalComponents[#,Method->"ConvexHull"],Max[#(1-#[[2,2]])CrossMatrix@1]&]+"#"#/.{0->"."})&[Characters@#/.{"."->0,"#"->1}]&
```
Explanation:
* First, `Characters@# /. {"."->0, "#"->1}` turns the input into a matrix of `0`s and `1`s.
* `"o"MorphologicalTransform[MorphologicalComponents[#,Method->"ConvexHull"],Max[#(1-#[[2,2]])CrossMatrix@1]&]+"#"#` then uses Mathematica's powerful (but extremely byte-heavy…) image processing capabilities to first fill in the island's [convex hull](https://en.wikipedia.org/wiki/Convex_hull) (which is the shape you'd get if you stretched a piece of string around it), and then take its boundary. We then multiply this matrix by the string `"o"` to get a matrix of `0`s and `"o"`s (thanks to Mathematica's impressive adaptability about types), and add this to `"#"` times the matrix of the island.
* Finally, `""<>#& /@ (... /. {0->"."})` turns this matrix of `"o"`s, `"#"`s and `0`s into a matrix of `"o"`s, `"#"`s and `"."`s, and joins each row into a string.
When we test this on example **B**, we get the output
```
{"....oo..",
 "...o##o.",
 "..o####o",
 ".o###..o",
 "o#####o.",
 "o#####o.",
 ".o##oo..",
 "..oo...."}
```
[Edit, thanks to Greg Martin:] If we're allowed to use arrays of characters instead of lists of strings, we can trim this down to 144 bytes:
```
"o"MorphologicalTransform[MorphologicalComponents[#,Method->"ConvexHull"],Max[#(1-#[[2,2]])CrossMatrix@1]&]+"#"#/.{0->"."}&[#/.{"."->0,"#"->1}]&
```
[Answer]
(But go upvote [Notatree's solution](https://codegolf.stackexchange.com/a/113637/56178), it's better!)
# Mathematica, 168 bytes
```
(x_~n~y_:=Min[Norm[x-#]&/@y];i=#;p=i~Position~#&;t=p["#"|"."]~Select~#&;(i~Part~##="o")&@@@t[#~n~t[ConvexHullMesh[Join[s=p@"#",#+{.1,.1}&/@s]]~RegionMember~#&]==1&];i)&
```
Pure function taking a 2D array of characters as input and returning a 2D array of characters. An easier-to-read version:
```
1  (x_~n~y_ := Min[Norm[x - #] & /@ y];
2  i = #; p = i~Position~# &; 
3  t = p["#" | "."]~Select~# &;
4  (i~Part~## = "o") & @@@ 
5    t[#~n~
6      t[ConvexHullMesh[
7        Join[s = p@"#", # + {.1, .1} & /@ s]]
8      ~RegionMember~# &] == 1 &];
9  i) &
```
Line 1 defines a function `n` that produces the (smallest) distance between a point `x` in the plane and a set `y` of other points. Line 2 initializes the variable `i` to the input, both to resolve a currying ambiguity later and to be able to modify it to produce the eventual output; line 2 also defines a function `p` that returns the coordinates of all occurrences of its input in `i`.
On line 3, `p["#" | "."]` represents every single coordinate from the input map (since all its characters are either `"#"` or `"."`), so `t` is a function that selects only those coordinates that satisfy a yet-unspecified property. On line 4, `i~Part~## = "o"` is going to change a bunch of entries of `i` to the character `"o"`; those characters will be selected from the set of possible coordinates according to the stuff on lines 5-8. And line 9 just returns the answer once it's computed.
Okay, infrastructure done, now to the actual computation. `ConvexHullMesh` is Mathematica's built-in for computing the convex hull of a set of points (the smallest convex polygon containing those points). Morally speaking, this should "fill in" the coves and fjords of the island (which is `s = p@"#"`), to rule them out of our cicrumnavigation. There's a little problem with `ConvexHullMesh` when that set of points is all in a line (thank you, test case #2), which we resolve by appending a slightly offset version of `s` to itself in line 7. This output is a polygon, so lines 7-9 (`t[...~RegionMember~# &]`) produces a list of the points with integer coordinates in that polygon. Finally, line 5 and the end of line 9 calculate all points that are at a distance exactly 1 (hence not 0) from this set of integer points; that set becomes the circumnavigating path.
Below is the output for the large test case at the OP's link. Notice at the upper left, the unusual choices of when to go west versus southwest hint at the fact that it's shadowing an invisible line of slope -2/3 between two peninsulas (said line segment being part of the convex hull's boundary).
```
........................
.............o..........
...........oo#ooooooo...
..........o#.#.##...#o..
........oo.#.#.###.##o..
.......o..########.##o..
.....oo...############o.
...oo#....############o.
..o#.###.##############o
.o##.##################o
.o####################o.
.o.##################.o.
.o##################..o.
.o..################..o.
o###################..o.
o#####################o.
o.##################.o..
o####################o..
o#...##############.o...
o##...#############o....
o#.....###....#oooo.....
.oooooo#ooooooo.........
.......o................
```
[Answer]
# Python 3, 779 bytes (indenting with tabs)
This is the whole program. It reads input from stdin and prints it to stdout. Stdin must end with EOF.
Example run with the big input: <https://ideone.com/XIfYY0>
```
import itertools,sys
L=list
C=itertools.count
d=L(map(L,filter(None,sys.stdin.read().split('\n'))))
X=len(d[0])
Y=len(d)
R=range
def P(r):return all(d[y][x]=='.'for x,y in r)
def S(f):
    for n in C(0):
        if P(f(n)):l=n
        else:break
    for n in C(l+1):
        if P(f(n)):return l,n
def f(V,a,*b):return L(eval('lambda '+a+':('+i+')',V)for i in b)
V=locals()
def D(n):
    y=min(n,Y-1);x=n-y
    while y>=0and x<X:yield(x,y);x+=1;y-=1
def E(n):
    x=max(0,n-Y);y=x+Y-n
    while y<Y and x<X:yield(x,y);x+=1;y+=1
F=f(V,'p','(p,y)for y in R(0,Y)','(x,p)for x in R(0,X)')+[D,E]
r=f(V,'x,y','x','y','x+y','x-y+Y')
B=L(map(S,F))
for x in R(0,X):
    for y in R(0,Y):
        z=L(zip(r,B))
        if all(g(x,y)in R(a,b+1)for g,(a,b)in z)and any(g(x,y)in e for g,e in z):d[y][x]='o'
print('\n'.join(''.join(x)for x in d))
```
The idea is simple: it computes smallest octagonal bounds, and draws cells which are inside all of the computed bounds and intersect at least one of the edges.
[Answer]
## JavaScript (ES6), ~~369~~ 343 bytes
```
f=s=>(a=s.split`
`.map(s=>[...s]),m=Array(8),a.map((b,i)=>b.map((c,j)=>c>'#'||[i-j,i,j+i,j,j-i,-i,-i-j,-j].map((d,k)=>d>m[k]||(m[k]=d-1)))),[o,p,q,r,t,u,v,w]=m,g=(i,j,k,l,...p)=>i-k|j-l?a[i][j]=g(i+(k>i)-(k<i),j+(l>j)-(l<j),k,l,...p):1/p[0]?g(k,l,...p):'o',g(p,p-o,p,q-p,q-r,r,r-t,r,-u,t-u,-u,u-v,w-v,-w,o-w,-w,p,p-o),a.map(b=>b.join``).join`
`)
```
Explanation: The string is split into a character array (I'm unclear as to whether character array input is acceptable). The array is then iterated over and the positions of all the land squares are located. The bounding lines given by the equations `x - y = o`, `x = p`, `x + y = q`, `y = r`, `y - x = t`, `-x = u`, `-x - y = v`, `-y = w` are determined such that the maximum possible parameter is chosen where all the land lies beyond the line. This has the effect of enclosing the island in an octagon. The coordinates of the corners of the octagon are readily calculated from the parameters and the cells on its edge are filled in. The array is then joined back into a string. The reason an octagon suffices is as follows:
```
   /o#     /o#     /o#
 |/o #   |/o #   |/ o#
 *o###   * o #   *  o#
/|o #   /|o #   /| o#
 |o#     |o#     |o#
```
Consider a corner of the octagon. At some point along the two edges the path will be limited by the land because we constructed the octagon to touch the land as closely as possible. If there is no land at the corner itself, the path could take the alternate routes as shown on the right, but that is still the same number of orthogonal and diagonal steps, so the distance is unchanged.
[Answer]
# Python 3.5, ~~224, 263, 234~~ 218 bytes
Golfed off another 16 bytes by getting rid of the nested function and making it a one-liner.
```
def h(s,k=0,i=0):w=s.find('\n')+1;x=s.find('X')-w;k=k or x;d=[1,w+1,w,w-1,-1,-w-1,-w,-w+1]*2;u=s[:k]+'o'+s[k+1:];return['X'>s[k]and i<8and(h(u,k+d[i+2],i+2)or h(u,k+d[i+1],i+1)or h(u,k+d[i],i))or'',s][s[k]>'X'and k==x]
```
Golfed off 29 bytes:
```
def f(s):
 w=s.find('\n')+1;x=s.find('X')-w;d=[1,w+1,w,w-1,-1,-w-1,-w,-w+1]*2
 def h(s,k,i):u=s[:k]+'o'+s[k+1:];return['X'>s[k]and i<8and(h(u,k+d[i+2],i+2)or h(u,k+d[i+1],i+1)or h(u,k+d[i],i))or'',s][s[k]>'X'and k==x]
 return h(s,x,0)
```
Input is a single string using '~' for ocean, 'X' for land, and 'o' for the boundary. (Using 'X' saves a byte for '>' instead of '==')
Less golfed version with comments:
```
def f(s):
    w=s.find('\n')+1                         # width of one row
    x=s.find('X')-w                          # starting point
    d=[1,w+1,w,w-1,-1,-w-1,-w,-w+1]*2        # delta to add to current index to move in 
                                             # the 8 directions: E, SE, S, SW, W, NW, 
                                             # N, NE. Make it long to avoid
                                             # lots of modulo operations in 
                                             #    the recursive calls
    def h(s,k,i):                            # s is the island string, k the current
                                             # position, i the direction index
        if s[k]>'X'and k==x:                 # if back at the begining,
            return s                         #   return the map
        elif 'X'>s[k] and i<8:               # if there is water here, and haven't
                                             #  looped around,
            u=s[:k]+'o'+s[k+1:]              #  make a new map with an 'o' in the 
                                             #  current spot
            r = h(u,k+d[i+2],i+2)            # try a 90 degree right turn
            if r: return r
            r = h(u,k+d[i+1],i+1)            # try a 45 degree turn
            if r: return r
            r= h(u,k+d[i],i)                 # try straight ahead
            if r: return r
        return ''                            # this path failed
    return h(s,x,0)
```
[Answer]
# C# 7 - ~~414~~ ~~369~~ 327 bytes
***Edit***: switched to 1D looping, computing `i` and `j` on the fly
***Edit***: changed input method, collapsed lookup table, and switched to well defined initial bounds ...and removed the pointless space in the last outer for-loop
```
using C=System.Console;class P{static void Main(){var D=C.In.ReadToEnd().Replace("\r","");int W=D.IndexOf('\n')+1,H=D.Length,z=H,k,q,c;int P()=>z%W*(k%3-1)+z/W*(k/3-1)+H;var B=new int[9];for(;z-->0;)for(k=9;k-->0&D[z]%7<1;)if(B[k]<=P())B[k]=P()+1;for(;++z<H;C.Write(q>9?'o':D[z]))for(q=k=9;k-->0;)q*=(c=P()-B[k])>0?0:c<0?1:2;}}
```
[Try It Online](https://tio.run/nexus/cs-core#lZJdT4MwFIbv9ysIZNKOD0EvzFbKkjETTDQuarKLbRcEukmYRT6cyrLfPimD6TaIWm7enue87wltt2@JTxechR8/k5S8qFZIk3BJkLt0koQbrZPUSX2XW4W@x905PgVwvXJibogt9YaqD8TxnsJr6gGY69el4xLAT2Ne5nmIfJpyYzzM@zzycT8H4pSKUNJlO6/dErpIn@UM23IgR7JbNI8AxGbWHndA0L5UdChl50yfF9pGbO4AU/LO5c2T7gzNwxigTFFMDUGmA9xFAdueDSfZrH1l6Aj6czCYBDMD5@GQKSYkfeeVpMywkaWOYz8lIDK7fTEUe8wMi8AI7yMRjDoYuMytsBhoan2t5xpaX@9doM1mu1UbVutvQKgHAvuEkp@CEv4AQrkaQQFLINSDKvmAtFjnUfEb1NR3UTX1esf@z@uB0OQQGof/x7EDJ9PVynFMKlAd7eEVtk4v9ZfH8AU)
Complete program, takes input to standard in, prints it to standard out, uses `#`, `.`, and `o`. For each cell, it computes a 'profile' (which is the distance over 8 directions (it appears to compute a ninth for convenience, but this is always `0`), and records a maximum of each of these. It then writes out the whole map again, and replaces any cell which is both on a boundary and not outside any with a 'o'. The commented code below explains how it all works.
As per my answer to [Save the Geese from Extinction](https://codegolf.stackexchange.com/questions/50829/save-the-geese-from-extinction), this produces the smallest octagon (valid circumnavigation with largest area) which bounds the island.
***Note***: that for once in my life I'm using something from the current decade, and this code requires C# 7 to compile. If you do not have C# 7, there is one line that will need to be replaced, which is clearly marked in the code.
Example usage and output:
```
type t7.txt | IslandGolf1.exe
.........ooooooooooo....
........o....#......o...
.......o...#.#.##...#o..
......o....#.#.###.##.o.
.....o....########.##..o
....o.....############.o
...o.#....############.o
..o#.###.##############o
.o##.##################o
o.####################.o
o..##################..o
o.##################...o
o...################...o
o###################...o
o#####################.o
o.##################..o.
o####################o..
o#...##############.o...
o##...#############o....
o#.....###....#...o.....
.o.....#.........o......
..ooooooooooooooo.......
```
Formatted and commented code:
```
using C=System.Console;
class P
{
    static void Main()
    {
        // \n 10
        // # 35
        // . 46
        // o 111
        var D=C.In.ReadToEnd().Replace("\r",""); // map
        int W=D.IndexOf('\n')+1, // width
            H=D.Length, // length
            z=H, // position in map (decomposed into i and j by and for P)
            k, // bound index
            q, // bound distance, and later cell condition (0 -> outside, 8 -> inside, >8 -> on boudary)
            c; // (free), comparison store
        // 'indexes' into a profile for the point z at index k
        // effectively {i=z%W,j=z/W,-i,-j,i+j,j-i,-i-j,i-j,0}[k] (inside order is a bit different) (0 const is always treated as 'inside bounds')
        // each non-zero-const entry describes the distance in one of the 8 directions: we want to maximise these to find the 'outer bounds'
        // the non-zero-const bounds describe 8 lines, together an octogen
        int P()=>z%W*(k%3-1)+z/W*(k/3-1)+H; // new C#7 local method syntax (if you don't have C#7, you can test this code with the line below instead)
        //k=0;System.Func<int>P=()=>z%W*(k%3-1)+z/W*(k/3-1)+H; // old lambda syntax (must pre-assign k to make static checker happy)
        var B=new int[9]; // our current bounds, each is initially null (must only call P() when on a #)
        // B[k] starts off a 0, P() has a +H term, and W+(H/W)<H for W >= 3, so B[k] is assigned the first time we compare it (H-i-j always > 0)
        for(;z-->0;) // for each cell
            for(k=9;k-->0& // for each bound
                D[z]%7<1;) // if this cell is #
                if(B[k]<=P())B[k]=P()+1; // update bound if necessary (add one so that we define the bound _outside_ the hashes)
        // z=-1
        for(;++z<H; // for each cell
                C.Write(q>9?'o':D[z])) // print the cell (if q > 9, then we are on the bounds, otherwise, spit out whatever we were before)
            // check we are not 'outside' any of the bounds, and that we are 'on' atleast one of them
            for(q=k=9;k-->0;) // for each bound
                q*=(c=P()-B[k])>0?0: // outside bound (q=0)    (??0 is cheaper than (int) or .Value)
                    c<0?1: // inside (preserve q)
                    2; // on bound (if q != 0, then q becomes > 9)
    }
}
```
]