After a short research, I found this site, which describes the algorithm. I didn't know it before. An important thing you didn't mention in your question is the calculation of the decision parameter $d$, which is:
$$d=f\left(N\right)+f\left(S\right)$$
It abuses the fact, that the bottom pixel will always lie on the circle or inside and the right pixel will always be on the circle or outside (for this particular segment). In conclusion, this means that $f\left(S\right) \le 0$ and $f\left(N\right) \ge 0$. If $d<=0$ you pick the right one and the left one otherwise.
Okay, we're using the circle equation to determine which is the closest point to the circle but why not simply using the euclidean distance?
Most probably for performance reasons. Using the Euclidean distance requires you to calculate 2 square roots, which is quite expensive when compared to simple multiplication ($r\cdot r$). On a Haswell CPU, you can perform 6 independent multiplications at the same time as calculating a single square root.
If you optimize the Bresenham algorithm you will probably get something like that:
$$\begin{matrix}
d &=& x_{k+1}^2 + y_k^2 - r^2 + x_{k+1}^2 + y_{k-1}^2 - r^2\\
&=& 2x_{k+1}^2- 2r^2 + y_k^2 + y_{k-1}^2 \\
&=& 2\left(x_{k+1}^2- r^2\right) + y_k^2 + y_{k-1}^2
\end{matrix}
$$
then compare
$$d \le 0$$
So in total, you end up with 9 multiplications, 2 additions, 1 subtraction and 1 comparison.
EDIT
As mentioned in the comments to this answer, the Euclidean length and the Bresenham algorithm won't give you the same results. While the Euclidean length will always give you the Pixel with the closest distance of its center to the arc, the Bresenham algorithm will favor points inside the circle, even though they are a little bit farther away. However, this error gets smaller the closer you get to the circle and shouldn't be noticeable.
Update
After rereading some tutorials like the one above or this one, I am not totally sure if they are really stating that they take the closest point. The tutorial I linked says:
So if $Dk < 0$ that means the negative $F(S)$ is bigger then the positive $F(N)$, that implies Point $N$ is closer to the circle than point $S$.
Language details are not my strong suit, but I think "implies" is meant more in the way of "indicating" something. In this case, it does not mean it is the exact result.
Anyways, after some elongated calculations, you are right, that $d=f\left(N\right)+f\left(S\right)$ will give wrong results in favor of the inner point if both distances to the circle are quite close. I did some error estimation.
With:
$$
\begin{matrix}
a &=& \sqrt{x_{k+1}^2 + y_k^2} \\
b &=& \sqrt{x_{k+1}^2 + y_{k-1}^2} \\
\end{matrix}
$$
being the distances to the points from the center, you can calculate the error $\epsilon$ as follows:
$$\begin{matrix}
a^2 - r^2 + b^2 - r^2 + \epsilon &=& ((a-r) - (r-b))^2\\
&=& (2r -a - b)^2
\end{matrix}$$
If you solve that for $\epsilon$ you will get:
$$\epsilon = 6r + 2ab - 4ra - 4 rb$$
You can get 2 interesting pieces of information from that. First, as you already claimed, the Bresenham algorithm does not take the exact distance to the circle when deciding which point is the next one, otherwise $\epsilon$ must have been 0. However, the second information is: If $a$ and $b$ are both on the circle, the error vanishes.
I was wondering how much closer $a$ must be than $b$ until it is selected. So I did the following:
$$\begin{matrix}
a &=& r+\delta + \tau\\
b &=& r-\delta
\end{matrix}$$
The reasoning here is as follows: If $\tau = 0$, then both points have the same distance $\delta$ to the circle. So the right-hand side of
$$\begin{matrix} a^2 - r^2 + b^2 - r^2 + \epsilon &=& (2r -a - b)^2 \end{matrix}$$
, the Euclidean distance, will be 0. Now I want to find an offset $\tau$, which causes the left-hand side to be 0, which is the condition where the Bresenham algorithm "assumes" both points to be equally distant to the circle. So if you say $$a^2 - r^2 + b^2 - r^2 = 0$$, replace the $a$ and $b$ with the shown formulas and solve for $\tau$ you get:
$$\tau = \sqrt{r^2 + 2\delta - \delta^2}- (r + \delta)$$
You can rewrite this a little bit so that it looks like this
$$\tau = \sqrt{(r + \delta)^2 - 2\delta^2}- \sqrt{(r + \delta)^2}$$
From this, you can see again that both approaches give the same result in the special case that both points lie on the circle since $\tau$ gets 0. Putting in some real numbers will also show, that $\tau$ does not get too big. I am no mathematician, but I think you can show, that the "error" and $\tau$ get relatively small with a growing ratio of $r:\delta$. So if you are not drawing circles with a radius of just 2-3 pixels, the differences in the visual output should be negligible.
Apart from that, taking the Euclidean distance might also not be the best way in the first place.
In the end, a rasterized circle is just an approximation and not a true circle. So there is no "correct" rasterization method for that. With enough distance, you won't notice that it contains of pixels. But in this case, you won't also notice anymore which pixels center is closer to the "true circle".
Conclusion
The Euclidean distance and Bresenham algorithm result in slightly different rasterizations (probably not noticeable most of the times).
The Euclidean distance is much more expensive to calculate due to the square roots. Additionally, as one can read in the linked sources, the Bresenham algorithm can be optimized even further than written in this answer, which makes it really fast. If I got it right, one just needs 1-2 multiplications and 2-3 additions/subtractions per pixel.
