You want to favor edges that are "close" to the line segment joining the endpoints of the path (its "axis," let's say). One direct way to do this is to weight the edges accordingly. How you weight them will determine what kind of path is "nicest." Just make sure that edges further from the axis get proportionately greater weight.
As an illustration, I have removed three points from the lattice shown in the question:
To obtain a "nicest" path from, say, (1,2) to (6,6), let's transform the coordinates to stretch the points quadratically relative to their distance to the axis:
(The quadratic stretching ensures that edge weights increase linearly with distance from the axis. Stronger transformations, such as exponential stretching, could be used if this does not appear sufficient to achieve the intended "niceness" of the ensuing solutions.)
Using the new distances for the edge weights, find the shortest path (using any efficient algorithm; Dijkstra's will do fine):
Although in the original lattice metric (where all edges have unit length) there are many shortest paths between (1,2) and (6,6)--such as the one going from (1,2) east to (6,2) and thence north to (6,6)--the quadratic distortion causes edges closer to the axis to be favored.
Note, please, that the lattice points themselves do not need to be changed: once the "nicest" path is found based on the reweighted edges, draw it in the original lattice (using the original vertex coordinates).
Evidently, each new origin-destination pair will require a reweighting of the edges. Although this may seem inefficient for large lattices, its signal advantage is that it automatically handles the second (more difficult) part of the question, where constraints (in the form of vertices to be avoided) are included: no special new algorithm is needed.
For completeness, here is the Mathematica code used to generate examples like this: it shows how to carry out the metric distortion.
m = 30; n = 20; (* Lattice dimensions *)
s = {m - 1, 7}; t = {2, n - 2}; (* Path endpoints *)
v = s - t;
v = {-v[[2]], v[[1]]} / Norm[v];
a = {{v[[2]], -v[[1]]}, v};
o = v.t;
vertices = Flatten[Table[{i, j}, {i, 1, m}, {j, 1, n}], 1];
vertices = Sort[RandomSample[vertices,
Length[vertices] - 5 Floor[Sqrt[Length[vertices]]]]];
With[{\[Epsilon] = 10/(m + n)}, coords = Transpose[a] . {#[[1]], (#[[2]] +
Sign[#[[2]] - o] \[Epsilon] (#[[2]] - o)^2)} & /@ (a.# & /@ vertices);
f[i_Integer, j_Integer] :=
If[Total[Abs[vertices[[i]] - vertices[[j]]]] != 1, Infinity, Norm[coords[[i]] - coords[[j]]]];
Show[ContourPlot[(#[[2]] - o + Sign[#[[2]] - o] \[Epsilon] (#[[2]] - o)^2) & @ (a.{x, y}),
{x, 1, m}, {y, 1, n}, AspectRatio -> n/m],
ListPlot[vertices, PlotStyle -> White], Epilog -> {Red, PointSize[0.015], Point[{s, t}]}]]
adj = Outer[f, Range[Length[vertices]], Range[Length[vertices]], 1];
h = WeightedAdjacencyGraph[adj , VertexCoordinates -> vertices]
g = WeightedAdjacencyGraph[adj , VertexCoordinates -> coords]
{v1, v2} = Flatten[Position[vertices, #] & /@ {s, t}];
Show[HighlightGraph[g, PathGraph[FindShortestPath[g, v1, v2]], GraphHighlightStyle -> "Thick"],
Graphics[{Gray, Dashed, Line[{coords[[v1]], coords[[v2]]}]}]]
Another example of its output, showing the grid and the solution in weighted and original coordinates:
