Python3, 336 bytes

E=enumerate
def f(m):
 S,M=(0,0),[-1,0,1]
 J,K=len(m),len(m[0])
 q,d,F=[(*S,V:=m[0][0])],{S:V},[]
 for x,y,c in q:
  if(x,y)==(J-1,K-1):F+=[c];continue
  for X in M:
   for Y in M:
    Q,W=x+X,y+Y
    if[X,Y]!=[0,0]and 0<=Q<J and 0<=W<K and((V:=d.get((Q,W),-1))==-1 or V>c+m[Q][W]):q+=[(Q,W,c+m[Q][W])];d[(Q,W)]=c+m[Q][W]
 return min(F)

Try it online!

APL (Dyalog Classic), 33 bytes

{⊃⌽,(⊢⌊⍵+(⍉3⌊/⊣/,⊢,⊢/)⍣2)⍣≡+\+⍀⍵}

Try it online!

{ } function with argument ⍵

+\+⍀⍵ take partial sums by row and by column to establish a pessimistic upper bound on path distances

( )⍣≡ repeat until convergence:

• (⍉3⌊/⊣/,⊢,⊢/)⍣2 min of distances to neighbours, i.e. do twice (( )⍣2): prepend leftmost column (⊣/,) to self (⊢) and append rightmost columns (,⊢/), find minima in horizontal triples (3⌊/) and transpose (⍉)

• ⍵+ add each node's value to its min of distances to neighbours

• ⊢⌊ try to beat the current best distances

⊃⌽, finally, return the bottom right cell

Octave + Image Processing package, 175 162 157 151 142 139 bytes

Saved 14 bytes thanks to @Luis Mendo and 1 byte thanks to @notjagan

function P(G)A=inf(z=size(G));A(1)=G(1);for k=G(:)'B=im2col(padarray(A,[1,1],inf),[3,3])+G(:)';B(5,:)-=G(:)';A=reshape(min(B),z);end,A(end)

Uses the Image Processing package, because why not? Isn't that how everybody solves graph problems?

Try it online!

Exploded

function P(G)
   A=inf(z=size(G));         % Initialize distance array to all Inf
   A(1)=G(1);                % Make A(1) = cost of start cell
   for k=G(:)'               % For a really long time...
      B=im2col(padarray(A,[1,1],inf),[3,3])+G(:)';
       %  B=padarray(A,[1,1],inf);     % Add border of Inf around distance array
       %  B=im2col(B,[3,3]);           % Turn each 3x3 neighborhood into a column
       %  B=B+G(:)';                   % Add the weights to each row
      B(5,:)-=G(:)';         % Subtract the weights from center of neighborhood
      A=reshape(min(B),z);   % Take minimum columnwise and reshape to original
   end
   A(end)                    % Display cost of getting to last cell

Explanation

Given an array of weights:

7   12    6    2    4
5   13    3   11    1
4    7    2    9    3
4    2   12   13    4
9    2    7    9    4

Initialize a cost array so that the cost to reach every element is Infinity, except the starting point (the upper left element) whose cost is equal to its weight.

  7   Inf   Inf   Inf   Inf
Inf   Inf   Inf   Inf   Inf
Inf   Inf   Inf   Inf   Inf
Inf   Inf   Inf   Inf   Inf
Inf   Inf   Inf   Inf   Inf

This is iteration 0. For each subsequent iteration, the cost to reach a cell is set to the minimum of:

• the current cost to reach that element, and
• the current cost to reach the element's neighbors + the weight of the element

After the first iteration, the cost of the path to element (2,2) (using 1-based indexing) will be

minimum([  7   Inf   Inf]   [13  13  13]) = 20
        [Inf   Inf   Inf] + [13   0  13]
        [Inf   Inf   Inf]   [13  13  13]

The full cost array after the first iteration would be:

  7    19   Inf   Inf   Inf
 12    20   Inf   Inf   Inf
Inf   Inf   Inf   Inf   Inf
Inf   Inf   Inf   Inf   Inf
Inf   Inf   Inf   Inf   Inf

After iteration k, each element will be the lowest cost of reaching that element from the start taking at most k steps. For example, the element at (3,3) can be reached in 2 steps (iterations) for a cost of 22:

  7    19    25   Inf   Inf
 12    20    22   Inf   Inf
 16    19    22   Inf   Inf
Inf   Inf   Inf   Inf   Inf
Inf   Inf   Inf   Inf   Inf

But on the 4th iteration, a path of 4 steps is found with a cost of 20:

 7   19   25   24   28
12   20   22   32   25
16   19   20   30   34
20   18   30   34   35
27   20   25   40   39

Since no path through the mxn matrix can be longer than the number of elements in the matrix (as a very loose upper bound), after m*n iterations every element will contain the cost of the shortest path to reach that element from the start.

Mathematica 279 Bytes

Basic idea is to create a graph with vertices corresponding to matrix entries and directed edges between any two vertices separated by a ChessboardDistance greater than zero but less than or equal to 1. Incidentally, this happens to be known as a King graph, since it corresponds to the valid moves of a king on a chessboard.

FindShortestPath is then used to get the minimal path. It works on EdgeWeight, not VertexWeight, so there is some extra code to define the EdgeWeight as the matrix entry corresponding to the destination of each directed edge.

Code:

(m=Flatten[#];d=Dimensions@#;s=Range[Times@@d];e=Select[Tuples[s,2],0<ChessboardDistance@@(#/.Thread[s->({Ceiling[#/d[[1]]],Mod[#,d[[1]],1]}&/@s)])≤1&];Tr[FindShortestPath[Graph[s,#[[1]]->#[[2]]&/@e,EdgeWeight->(Last@#&/@Map[Extract[m,#]&,e,{2}])],1,Last@s]/.Thread[s->m]])&

Note that the  character is the transpose symbol. It will paste into Mathematica as-is.

Usage:

%@{{2, 55, 5, 3, 1, 1, 4, 1},
  {2, 56, 1, 99, 99, 99, 99, 5},
  {3, 57, 5, 2, 2, 2, 99, 1},
  {3, 58, 4, 2, 8, 1, 99, 2},
  {4, 65, 66, 67, 68, 3, 99, 3},
  {2, 5, 4, 3, 3, 4, 99, 5},
  {75, 76, 77, 78, 79, 80, 81, 2},
  {5, 4, 5, 1, 1, 3, 3, 2}}

Output:

67

If you set g=Graph[...,GraphLayout->{"GridEmbedding","Dimension"->d},VertexLabels->Thread[s->m] and p=FindShortestPath[... then the following graphic will visually display the solution (top of the matrix corresponds to the bottom of the graph):

HighlightGraph[g,PathGraph[p,Thread[Most@p->Rest@p]]]

Python 2, 356 320 bytes

s=input()
r=lambda x:[x-1,x,x+1][-x-2:]
w=lambda z:[z+[(x,y)]for x in r(z[-1][0])for y in r(z[-1][1])if x<len(s)>0==((x,y)in z)<len(s[0])>y]
l=len(s)-1,len(s[0])-1
f=lambda x:all(l in y for y in x)and x or f([a for b in[l in z and[z]or w(z)for z in x]for a in b])
print min(sum(s[a][b]for(a,b)in x)for x in f([[(0,0)]]))

Try it here!

-36 bytes thanks to notjagan!

Receives a list of lists as an input, and outputs the lowest cost when navigating the matrix from the upper left to the bottom right.

Explanation

Find every possible route from the upper left to the bottom right of the matrix, creating a list of x,y coordinates for each route. The routes cannot backtrack, and they must end at (len(s)-1,len(s[0])-1).

Sum the integers on each path of coordinates, and return the minimum cost.

The print can be easily changed to output the list of coordinates for the shortest route.

Haskell, 228 bytes

Positions are lists of two elements, because those are easy to generate with sequence and just as easy to pattern match as 2-tuples.

h=g[[-1,-1]]
g t@(p:r)c|p==m=0|1<2=minimum$(sum$concat c):(\q@[a,b]->c!!a!!b+g(q:t)c)#(f(e$s$(\x->[0..x])#m)$f(not.e t)$zipWith(+)p#s[[-1..1],[-1..1]])where m=[l(c)-1,l(head c)-1]
(#)=map
f=filter
e=flip elem
s=sequence
l=length

Start at -1,-1 and count the cost of each steps destination field.

Alternative first two lines: start at 0,0, count the departure fields, terminate at the coordinates equal to the matrix dimensions (so down-right from the goal, which needs to be added to the list of legal destinations) - exact same length but slower:

j=i[[0,0]]
i t@(p@[a,b]:r)c|p==m=0|1<2=c!!a!!b+(minimum$(sum$concat c):(\q->i(q:t)c)#(f(e$m:(s$(\x->[0..x-1])#m))$f(not.e t)$zipWith(+)p#s[[-1..1],[-1..1]]))where m=[l c,l$head c]

Using an infix for map does not save bytes here but I substitute it as soon as it doesn't cost one, because it can only get better with more uses, and sometimes with other restructurings as well which shave off another pair of parentheses.

To be improved: Redundant filters. Merging/in-lining them to filter(flip elem$(s$(\x->[0..x])#m)\\p) with import Data.List for \\ costs 3 bytes.

Also, too bad (fromEnumTo 0) is 2 bytes longer than (\x->[0..x]).

sum$concat c is all fields' cost summed up and thus a concisely expressible upper bound on the path cost which is given to the minimum to avoid an empty list (my type checker has already determined the whole thing to work on Integers, so no hard-coding the maximum, hehe). No matter how I restrict steps based on the previous one made (which would speed up the algorithm a lot, but also cost bytes), I can not avoid the dead ends that make this fall-back necessary.

• One filter idea was ((not.e n).zipWith(-)(head r)) with extracting n=s[[-1..1],[-1..1]], which necessitates adding ,[-1,-1] to the initial path. The algorithm then avoids going where it could already have gone in the previous step, which makes stepping on an edge field orthogonally to that edge a dead end.

• Another was ((>=0).sum.z(*)d) with extracting z=zipWith, which introduces a new argument d to the recursive function that is given as (z(-)p q) in the recursion and [1,1] in the initial case. The algorithm avoids successive steps with a negative scalar product (d being the previous step), which means no sharp 45°-turns. This still narrows down the choices considerably, and avoids the previous trivial dead end, but there are still paths that end up enclosed in already-visited fields (and possibly an 'escape' which however would be a sharp turn).

JavaScript, 442 412 408 358 bytes

This is my first PPCG submission. Feedback would be appreciated.

(m,h=m.length,w=m[0].length)=>{for(i=0;i<h*w;i++)for(x=0;x<w;x++){for(y=0;y<h;y++){if(m[y][x]%1==0)m[y][x]={c:m[y][x],t:m[y][x]};for(X=-1;X<=1;X++)for(Y=-1;Y<=1;Y++){t=x+X;v=y+Y;if((X==0&&Y==0)||t<0||t>=w||v<0||v>=h)continue;if(m[v][t]%1==0)m[v][t]={c:m[v][t],t:null};c=m[y][x].t+m[v][t].c;if (c<m[v][t].t||m[v][t].t==null)m[v][t].t=c}}}return m[h-1][w-1].t}

This takes a multi-dimensional array as input.

Explanation

Basically, loop through all of the cells over and over adjusting the lowest known cost to get to each of the neighbors. Eventually, the grid will reach a state where the total cost to reach the bottom right is the lowest cost to get there.

Demo

f=(m,h=m.length,w=m[0].length)=>{for(i=0;i<h*w;i++)for(x=0;x<w;x++){for(y=0;y<h;y++){if(m[y][x]%1==0)m[y][x]={c:m[y][x],t:m[y][x]};for(X=-1;X<=1;X++)for(Y=-1;Y<=1;Y++){t=x+X;v=y+Y;if((X==0&&Y==0)||t<0||t>=w||v<0||v>=h)continue;if(m[v][t]%1==0)m[v][t]={c:m[v][t],t:null};c=m[y][x].t+m[v][t].c;if (c<m[v][t].t||m[v][t].t==null)m[v][t].t=c}}}return m[h-1][w-1].t}

//Tests
console.log(f([[1,1,1],[1,1,1]])===3);
console.log(f([[7,9,6,6,4],[6,5,9,1,6],[10,7,10,4,3],[4,2,2,3,7],[9,2,7,9,4]])===28);
console.log(f([[2,42,6,4,1],[3,33,1,1,1],[4,21,7,59,1],[1,7,6,49,1],[1,9,2,39,1]])===27);
console.log(f([[5,6,7,4,4],[12,12,25,25,25],[9,4,25,9,5],[7,4,25,1,12],[4,4,4,4,4]])===34); 
console.log(f([[1,1,1,1],[9,9,9,1],[1,9,9,9],[1,9,9,9],[1,1,1,1]])===15)
console.log(f([[2,55,5,3,1,1,4,1],[2,56,1,99,99,99,99,5],[3,57,5,2,2,2,99,1],[3,58,4,2,8,1,99,2],[4,65,66,67,68,3,99,3],[2,5,4,3,3,4,99,5],[75,76,77,78,79,80,81,2],[5,4,5,1,1,3,3,2]])===67);

Edit: Special thanks to @ETHproductions for helping me shave dozens of tasty bytes.

More thanks to @Stewie Griffin for your tips that knocked off 50 bytes.

Python 3 + numpy + scipy, 239 222 186 bytes

from numpy import*
from scipy.sparse.csgraph import*
def f(M):m,n=s=M.shape;x,y=indices(s);return dijkstra([(M*(abs(i//n-x)<2)*(abs(i%n-y)<2)).flatten()for i in range(m*n)])[0,-1]+M[0,0]

Try it online!

JavaScript, 197 bytes

a=>(v=a.map(x=>x.map(_=>1/0)),v[0][0]=a[0][0],q=[...(a+'')].map(_=>v=v.map((l,y)=>l.map((c,x)=>Math.min(c,...[...'012345678'].map(c=>a[y][x]+((v[y+(c/3|0)-1]||[])[x+c%3-1]||1/0)))))),v.pop().pop())

Prettify:

a=>(
  // v is a matrix holds minimal distance to the left top
  v=a.map(x=>x.map(_=>1/0)),
  v[0][0]=a[0][0],
  q=[
     // iterate more than width * height times to ensure the answer is correct
    ...(a+'')
  ].map(_=>
    v=v.map((l,y)=>
      l.map((c,x)=>
        // update each cell
        Math.min(c,...[...'012345678'].map(
          c=>a[y][x]+((v[y+(c/3|0)-1]||[])[x+c%3-1]||1/0)
        ))
      )
    )
  ),
  // get result at right bottom
  v.pop().pop()
)