julia, 104 97 bytes

This solution uses the quotient space metric. For the projective plane, if \$\sim\$ is the equivalence relation induced by the edge identification (i.e. \$(p, 0) \sim (N-p, N); (0, p) \sim (N, N-p)\$) and \$d\$ is the Euclidean metric, then the projective plane metric is $$ d'(a, b) = \min_{w \sim w'} [d(a, w) + d(b, w')]. $$

To implement this, we create two grids \$x\$ and \$y\$ spanning \$[0, N]^2\$ such that \$x_{ij} \sim y_{ij}\$ and compute \$d(a, x_{ij}) + d(b, y_{ij})\$ for each corresponding pair of grid entries.

We make the vertical spacing between grid points \$0.005\$ to make sure that the result is correct to two decimal places.

u&v=((h=u-v.-.5)'h)^.5 g(N,a,b)=minimum(w->a&w+b&((N.-2w)any(w.%N.==0)+w),vcat.((r=0:.005:N)',r)) 

Attempt This Online!

Thanks to MarcMush for a 7-byte improvement!

How?

# u&v is the distance between u-0.5 (centre of cell u) and v. u&v = ((h = u - v .- .5)'h)^.5 # Minimise the function w -> d(a, w) + d(b, w') g(N, a, b) = minimum( # (N .- 2w)any(w .% N .== 0) + w is N-w on the edges of x and w otherwise. w -> a&w + b&((N .- 2w)any(w .% N .== 0) + w), # Map over a grid of pairs [x1, x2] spanning [0, N]^2. vcat.((r = 0:.005:N)', r) ) 

Jelly, 20 bytes

_~jUŒHƲ€N⁹ṭ"_þ/ẎÆḊ€Ṃ 

A dyadic Link that accepts the size on the left and a pair of pairs (the points) on the right and yields the Euclidean distance of the projective plane.

Try it online! Or see the test-suite.

How?

_~jUŒHƲ€N⁹ṭ"_þ/ẎÆḊ€Ṃ - Link: Size = S; Points = [[x1, y1], [x2, y2]] _ - {Size} subtract {Points} -> [[S-x1, S-y1], [S-x2, S-y2]] € - for each of these pairs, [a, b]: Ʋ - last four links as a monad: ~ - bitwise NOT (vectorises) -> [-1-a, -1-b] U - upend -> [b,a] j - join -> [-1-a, b, a, -1-b] ŒH - halve -> [[-1-a, b], [a, -1-b]] N - negate -> [[a+1, -b], [-a, b+1]] ⁹ - chain's right argument -> [[x1, y1], [x2, y2]] " - zip with: ṭ - tack -> [[[S-x1+1, y1-S], [x1-S, S-y1+1], [x1, y1]], [[S-x2+1, y2-S], [x2-S, S-y2+1], [x2, y2]], ] / - reduce by: þ - table with: _ - subtraction Ẏ - tighten ÆḊ€ - vector norm of each Ṃ - minimum 

Python 3, 94 bytes

lambda n,X,Y,x,y:min(abs(X-a+(Y-b)*1j)for k in[-n,n]for a,b in[(x,y),(n+1-x,y+k),(x+k,n+1-y)]) 

Try it online!

Easier to understand as:

lambda n,X,Y,x,y:min(abs(X-a+(Y-b)*1j)for a,b in[(x,y),(n+1-x,y-n),(n+1-x,y+n),(x-n,n+1-y),(x+n,n+1-y)]) 

Try it online!

Fixes one point in place, and picks the closest version of the other point out of 5 options -- the original and the four reflections across the edges of the square. There's no need to consider diagonally-adjacent copies because they're always further from the fixed point than the original.

Python, 83 bytes (-35 thanks to @xnor)

lambda N,A,B:min(map(abs,(B-A,C:=B-1j*N-~N-2*B.real-A,D:=1+1j-C-2*A,C+2j*N,D+2*N))) 

Attempt This Online!

Obsolete Python, 118 bytes

lambda N,A,B,o=.5+.5j:min(map(lambda x:abs(A-o-x),[B:=B-o,C:=N/o+B-2*B.real,-B,2*N-B,2j*N-B,4*N*o-B,-C,C+2j*N,2*N-C])) 

Attempt This Online!

Expects input coordinates as complex numbers.

Python3, 458 bytes

R=range E=enumerate def T(d,x,y,X,Y): if y%2:d=[[(*a,W+Y*len(d),k)for k,(*a,_,_)in E(i[::-1])]for W,i in E(d)] if x%2:d=[[(*a,W,k+X*len(i))for k,(*a,_,_)in E(i)]for W,i in E(d[::-1])] return[j for k in d for j in k] def f(n,p,P): d=[[(i,j,i-1,j-1)for j in R(1,n+1)]for i in R(1,n+1)];O=[*T(d,0,0,0,0),*T(d,0,1,0,-1),*T(d,1,0,1,0),*T(d,0,1,0,1),*T(d,1,0,-1,0)] return min(((X-x)**2+(Y-y)**2)**0.5 for a,b,x,y in O for A,B,X,Y in O if[a,b]==p and[A,B]==P) 

Try it online!

Charcoal, 41 bytes

⊞υηＦ⟦θ±θ⟧Ｆ²⊞υＥη⎇⁼κμ⁻⊕θλ⁺λιＩ₂⌊ＥυΣＸＥι⁻λ§ζμ² 

Try it online! Link is to verbose version of code. Takes the size and the two points as input. Explanation:

⊞υη 

Push the first point to the predefined empty list.

Ｆ⟦θ±θ⟧Ｆ² 

Loop over each orthogonal direction.

⊞υＥη⎇⁼κμ⁻⊕θλ⁺λι 

Push the transform of the point in that direction to the list.

Ｉ₂⌊ＥυΣＸＥι⁻λ§ζμ² 

Output the minimum Euclidean distance from the second point to any of the five points inthe list.

I tried writing it as a single expression but it still came out to 41 bytes:

Ｉ₂⌊Ｅ⊞ＯＥ⁴Ｅη⎇﹪⁺ιμ²⁻⊕θλ⎇÷ι²⁺λθ⁻λθηΣＸＥι⁻λ§ζμ² 

Try it online! Link is to verbose version of code.

JavaScript (ES6), 104 bytes

Saved 9 bytes by refactoring in a way similar to what xnor did

Expects (N,x1,y1,x2,y2).

(N,x,y,X,Y)=>Math.min(...(g=q=>[x,x+q,N-x].map((v,i)=>Math.hypot(X-v,Y-[y,N-y,y+q][i])))(N++),...g(1-N)) 

Try it online!

Scala, 168 bytes

Golfed version. Try it online!

def f(N:Int,x:Int,y:Int,X:Int,Y:Int)={val(a,b)=(x-X,Y+y-N-1);(Seq(a,a-N,X+x-N-1,X+x-N-1,X+x-N-1,a+N)zip Seq(y-Y,b,y-Y-N,N-b,y-Y+N,b)).map{case(i,j)=>sqrt(i*i+j*j)}.min} 

Ungolfed version. Try it online!

import scala.math.{sqrt, min} object Main { def f(N: Int, x: Int, y: Int, X: Int, Y: Int): Double = { val a = Array(x - X, x - X - N, X + x - N - 1, X + x - N - 1, X + x - N - 1, x - X + N) val b = Array(y - Y, Y + y - N - 1, y - Y - N, N - (Y + y - N - 1), y - Y + N, Y + y - N - 1) val distances = (a zip b).map{ case (i, j) => sqrt(i*i + j*j) } distances.min } def main(args: Array[String]): Unit = { println(f(5, 1, 1, 1, 1)) // 0.0 println(f(5, 2, 4, 5, 3)) // 2.24 println(f(5, 2, 4, 4, 1)) // 2.0 println(f(10, 2, 2, 8, 8)) // 4.12 println(f(10, 1, 1, 10, 10)) // 1.0 println(f(10, 4, 4, 6, 6)) // 2.83 } }