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$.