The calculation follow the steps of @mixedmath:
Let's denote the points of the rectangle with
\begin{align*} A&=(625.49939,632.40015)=(x_A,y_A)\\ B&=(636.49939,632.40015)=(x_B,y_A)\\ C&=(636.49939,672.40015)=(x_B,y_C)\\ D&=(625.49939,672.40015)=(x_A,y_C)\\ \end{align*}
1. Step: Center $M$
We calculate the center $M$ of the rectangle $\Box(A,B,C,D)$: \begin{align*} M=\frac{1}{2}(A+C)=\frac{1}{2}(B+D)=\frac{1}{2}(x_A+x_B,y_A+y_C) \end{align*}
2. Step: Shift $M$ to $(0,0)$
We shift the rectangle so that $M$ becomes the origin and obtain a rectangle $\Box(A^\prime, B^\prime, C^\prime, D^\prime)$ \begin{align*} A^\prime&=A-M=\frac{1}{2}(x_A-x_B,y_A-y_C)\\ B^\prime&=B-M=\frac{1}{2}(-x_A+x_B,y_A-y_C)\\ C^\prime&=C-M=\frac{1}{2}(-x_A+x_B,-y_A+y_C)\\ D^\prime&=D-M=\frac{1}{2}(x_A-x_B,-y_A+y_C)\\ \end{align*}
Intermezzo: Rotation matrix
We consider the rotation matrix $\Phi$ clockwise through $90^\circ$ which corresponds to $-\frac{\pi}{2}$ \begin{align*} \Phi=\begin{pmatrix} \cos \phi&-\sin \phi\\ \sin \phi& \cos \phi \end{pmatrix} = \begin{pmatrix} \cos\left(-\frac{\pi}{2}\right)&-\sin\left( -\frac{\pi}{2}\right)\\ \sin\left( -\frac{\pi}{2}\right)& \cos\left( -\frac{\pi}{2}\right) \end{pmatrix} = \begin{pmatrix} 0&1\\ -1& 0 \end{pmatrix} \end{align*}
A point $P=(x_P,y_P)$ will be transformed via $\Phi$ to
\begin{align*} \Phi P=\begin{pmatrix} 0&1\\ -1& 0 \end{pmatrix} \begin{pmatrix} x_P\\ y_P \end{pmatrix} = \begin{pmatrix} y_P\\ -x_P \end{pmatrix} \end{align*}
3. Step: Rotate rectangle clockwise through $-\frac{\pi}{2}$
We now rotate the rectangle $\Box(A^\prime, B^\prime, C^\prime, D^\prime)$ and obtain a rectangle $\Box(\Phi A^\prime, \Phi B^\prime, \Phi C^\prime ,\Phi D^\prime)$ \begin{align*} \Phi A^\prime&=\frac{1}{2}(y_A-y_C,-x_A+x_B)\\ \Phi B^\prime&=\frac{1}{2}(y_A-y_C,x_A-x_B)\\ \Phi C^\prime&=\frac{1}{2}(-y_A+y_C,x_A-x_B)\\ \Phi D^\prime&=\frac{1}{2}(-y_A+y_C,-x_A-x_B)\\ \end{align*}
4. Step: Shift rectangle back from origin to $M$
We reverse step 2 by adding $M$ to the points of $\Box(\Phi A^\prime, \Phi B^\prime, \Phi C^\prime ,\Phi D^\prime)$ \begin{align*} \Phi A^\prime+M&=\frac{1}{2}(x_A+x_B+y_A-y_C,-x_A+x_B+y_A+y_C)=(610.99939,657.90015)\\ \Phi B^\prime+M&=\frac{1}{2}(x_A+x_B+y_A-y_C,x_A-x_B+y_A+y_C)=(610.99939,646.90015)\\ \Phi C^\prime+M&=\frac{1}{2}(x_A+x_B-y_A+y_C,x_A-x_B+y_A+y_C)=(650.99939,646.90015)\\ \Phi D^\prime+M&=\frac{1}{2}(x_A+x_B-y_A+y_C,-x_A-x_B+y_A+y_C)=(650.99939,657.90015)\\ \end{align*}
Note: A code can be easily derived from the calculations
[2015-07-28] Add-on: According to @BROYs comment here's the rotational displacement of the rectangle using homogeneous matrices.
A translational displacement by a vector $U=\begin{pmatrix}x_U\\y_U\end{pmatrix}$ can be represented by the homogeneous matrix $T_U$ with \begin{align*} T_U=\begin{pmatrix} 1&0&x_U\\ 0&1&y_U\\ 0&0&1\\ \end{pmatrix} \end{align*}
The point $A=(x_A,y_A)$ in homogeneous coordinates is $\begin{pmatrix}x_A\\y_A\\1\end{pmatrix}$. The translation of $A$ in direction $U$ is \begin{align*} T_UA=\begin{pmatrix} 1&0&x_U\\ 0&1&y_U\\ 0&0&1\\ \end{pmatrix} \begin{pmatrix} x_A\\ y_A\\ 1 \end{pmatrix} = \begin{pmatrix} x_A+x_U\\ y_A+y_U\\ 1 \end{pmatrix} \end{align*} Observe $T_UA$ corresponds to $A+U$.
A rotational displacement by a vector $\Phi=-\frac{\pi}{2}$ can be represented by the homogeneous matrix $R_\Phi=R_{-\frac{\pi}{2}}$ with \begin{align*} T_\Phi&=\begin{pmatrix} \cos(\Phi)&-\sin(\Phi)&0\\ \sin(\Phi)&\cos(\Phi)&0\\ 0&0&1\\ \end{pmatrix}\\ T_{-\frac{\pi}{2}}&= \begin{pmatrix} 0&1&0\\ -1&0&0\\ 0&0&1\\ \end{pmatrix} \end{align*}
The rotation of $A$ through an angle $\Phi=-\frac{\pi}{2}$ \begin{align*} T_{-\frac{\pi}{2}}A= \begin{pmatrix} 0&1&0\\ -1&0&0\\ 0&0&1\\ \end{pmatrix} \begin{pmatrix} x_A\\ y_A\\ 1 \end{pmatrix} = \begin{pmatrix} y_A\\ -x_A\\ 1 \end{pmatrix} \end{align*}
Note: In the following it's more convenient to use vector notation for the points. I will omit the transpose symbol.
Homogeneous matrices in action
The players: The four points of the rectangle $A,B,C,D$ and the center $M$ of the rectangle
\begin{align*} A&=\begin{pmatrix} x_A\\ y_A\\ 1 \end{pmatrix} \qquad B=\begin{pmatrix} x_B\\ y_B\\ 1 \end{pmatrix} \qquad C=\begin{pmatrix} x_C\\ y_C\\ 1 \end{pmatrix} \qquad D=\begin{pmatrix} x_D\\ y_D\\ 1 \end{pmatrix}\\ M&=\frac{1}{2}(A+C)=\frac{1}{2}(B+D)= \frac{1}{2}\begin{pmatrix} x_A+x_B\\ y_A+y_C\\ 1 \end{pmatrix} \end{align*}
First we want to move $M$ to the origin. This is a translational displacement by $-M$ represented by the homogeneous matrix $T_{-M}$. Next we rotate through $-\frac{\pi}{2}$ represented by $R_{-\frac{\pi}{2}}$ and at the end we shift back by $T_M$. The complete displacement is given by
\begin{align*} T_M R_{-\frac{\pi}{2}} T_{-M}&= \begin{pmatrix} 1&0&\frac{1}{2}(x_A+x_B)\\ 0&1&\frac{1}{2}(y_A+y_C)\\ 0&0&1 \end{pmatrix} \begin{pmatrix} 0&1&0\\ -1&0&0\\ 0&0&1\\ \end{pmatrix} \begin{pmatrix} 1&0&-\frac{1}{2}(x_A+x_B)\\ 0&1&-\frac{1}{2}(y_A+y_C)\\ 0&0&1 \end{pmatrix}\\ &= \begin{pmatrix} 1&0&\frac{1}{2}(x_A+x_B)\\ 0&1&\frac{1}{2}(y_A+y_C)\\ 0&0&1 \end{pmatrix} \begin{pmatrix} 0&1&-\frac{1}{2}(y_A+y_C)\\ -1&0&-\frac{1}{2}(x_A+x_B)\\ 0&0&1 \end{pmatrix}\\ &= \begin{pmatrix} 0&1&-\frac{1}{2}(x_A+x_B-y_A-y_C)\\ -1&0&-\frac{1}{2}(x_A+x_B+y_A+y_C)\\ 0&0&1 \end{pmatrix} \end{align*}
The complete displacement of the point $A$ is obtained by
\begin{align*} T_MR_{-\frac{\pi}{2}}T_{-M}A&=T_MR_{-\frac{\pi}{2}}T_{-M}\begin{pmatrix} x_A\\ y_A\\ 1 \end{pmatrix}= \frac{1}{2}\begin{pmatrix} x_A+x_B+y_A-y_C\\ -x_A+x_B+y_A+y_C\\ 2 \end{pmatrix}= \begin{pmatrix} 610.99939\\ 657.90015\\ 1 \end{pmatrix} \end{align*}
and similarly for $B,C,D$.