Deriving the non-paraxial form of Rayleigh criterion Background
The Rayleigh criterion of imaging resolution says that two incoherent point sources are barely resolved by a diffraction-limited system with a circular aperture where the center of the Airy intensity pattern generated by one point source falls exactly on the first zero of the Airy pattern generated by the second. Usually the Airy pattern of a circular aperture is derived as the Fraunhofer diffraction of the aperture, with the following result:

$I(r) = \left(\frac{A}{\lambda z}\right)^2 \left[2 \frac{J_1(kwr/z)}{kwr/z} \right]^2$ 

where $A$ is the area of the aperture, $\lambda$ and $k$ are wavelength and wavenumber respectively, $w$ is the radius of the aperture, $r$ is the radius coordinate of the observation plane, and $z$ is the distance between the aperture and the observation plane. From this expression, one can derive that the Rayleigh resolution is 

$\delta = 0.61 \frac{\lambda z}{w}$

Question
Fraunhofer diffraction or Fourier optics in general is based on the paraxial approximation of the Kirchhoff diffraction formula. So you cannot trust the expression above to be valid in the non-paraxial situation often encountered in microscopy. In fact, "Introduction to Fourier Optics" by Goodman says (pg 158) the corresponding result in the nonparaxial case can be shown to be 

$\delta = 0.61 \frac{\lambda}{\sin\theta} = 0.61 \frac{\lambda}{NA}$

where $NA$ is the numerical aperture. But the textbook does not provide any reference for this derivation. If you look at Wikipedia article on Rayleigh criterion, you will see that $z/w$ is related to  F-number, and $\text{F-number}\approx 1/(2NA)$, but again that relation only holds up on the paraxial case. I would like to see a full derivation of the non-paraxial expression for the Rayleigh resolution. 
Personal attempt at solution
Often, a high NA objective is paired with a long focal length tube lens. So when we measure the point spread function (PSF) of a point source with this (objective + tube lens) system, we are measuring the magnified version of the true PSF. You can say that in most cases, the paraxial approximation is valid on the tube lens side. That is, the PSF we measure is described by Rayleigh resolution of $0.61 \frac{\lambda}{NA_{\text{imaging space}}}$. If we can say that magnification of the imaging system is $M = \frac{NA_{\text{objective space}}}{NA_{\text{imaging space}}}$, then we can say that the Rayleigh resolution of the high NA objective is indeed $0.61 \frac{\lambda}{NA}$. The fact that the magnification is described as a ratio of two numerical apertures is described by the Abbe sine condition, which must be obeyed by diffraction-limited imaging system. 
I think the above argument is sound, but it relies on the fact that you use a long focal lens tube lens paired with an infinity-corrected objective. I would like to see a more general derivation.
 A: The reason for long focal length tube lens is to achieve magnification, usually considered a plus for a microscope. You could use a shorter focal length tube lens and have less magnification. The sine condition is valid independent of the relative focal lengths and does not rely on paraxial rays. It was introduced by Abbe and is about as first principle as possible. 
A: In all the diagrams (usually the same diagram repeated) used to explain the derivation of the Rayleigh criterion, I find that the geometry is explained incorrectly, but no one ever seem to comment on this or even notice it. From the diagram shown, for example, here:

(Source)
it is clear to me that $(0.5D)/d = \tan(a)$, not $\sin(a)$.  The approximation that $\tan(a)=\sin(a)$ is only valid for small angles. The derivation of the equation $r=0.61(\lambda)/NA$ is only possible from Rayleigh's angular resolution equation, $\sin(\theta R)=1.22(\lambda)/D$, if $\lambda\ll D\ll d$, meaning that both $a$ and $\theta R$ are small.
