Necessary and sufficient condition for quadrilateral to be cyclic Can you provide a proof for the following proposition:

Proposition. Given any quadrilateral $ABCD$. Let $P,Q,R,S$ be nine-point centers of triangles $\triangle ABD$,$\triangle ABC$,$\triangle BCD$ and $\triangle ACD$ respectively. Then, the quadrilateral $ABCD$ is cyclic if and only if $P,Q,R,S$ are concyclic.


GeoGebra applet that demonstrates this proposition can be found here.
 A: Regarding the "only if" part, a stronger result actually holds: if the quadrilateral $ABCD$ is cyclic, then $PQRS$ is similar to $ABCD$, and so it is cyclic, too.
This is stated, without proof, in many elementary geometry textbooks, see for instance at p. 44 of
D. G. Wells: The Penguin dictionary of curious and interesting geometry, New York, NY: Penguin Books. xiv, 285 p. (1991). ZBL0856.00005.
A proof can be found at p. 36 of the paper
F. V. Morley: Notes on the cyclic quadrilateral, Annals of Math. (2) 22, 35-42, 43 (1920). ZBL47.0566.01.
A: This is not an answer but a comment. However, it will be too long and it would be awkward to break it up.  It is of course natural to ask what is special about the nine point centre here.  One can put the question in the following context.  Given the shape of a quadrilateral $ABCD$, what is the shape of  $PQRS$ (constructed as above but using any of the triangle centres protocolled in the online Encyclopedia of Triangle Centers)?
The shape of a quadrilateral is the unique  pair $(p_1,q_1)$ and $(p_2,q_2)$ for which it is similar to the one with $(0,0)$, $(1,0)$ and these two points as vertices (in terms of complex numbers they are just
$$   \frac{z_C-z_A}{z_B-z_A},\frac{z_D-z_A}{z_B-z_A}.      $$
Many structural properties of a quadrilateral (in particular cyclicity) can be expressed as a simple equation in the $p$'s and $q$'s.  For cyclicity, one equates the coordinates of the circumcentres of $ABC$ and $ABD$.
It is then a simple, if usually tedious, task to compute the shape of $PQRS$ (for a given centre function) in terms of that of $ABCD$  (easily automatised using Mathematica) and so provide far-reaching generalisations of Morley's result.
To be explicit, if $Z_1$ and $Z_2$ are $p_1+iq_1$, resp. $p_2+iq_2$, then we can easily compute the shape of $PQRS$ by computing the complex numbers which specify its vertices (using the centre function) and then forming the corresponding quotients as above.
