I think a large part of the difficulty here is the variation in
the magnitudes of your "normals" in general.
When you give a "normal" in the form $\langle -f_x, - f_y, 1\rangle$,
basically what you have is a radial projection of
the hemisphere parameterized by
$\langle \sin\phi \cos\theta, \sin\phi \sin\theta, \cos\phi\rangle$
(for $0 \leq \phi \leq \frac\pi2$) onto the plane $z = 1$.
That is, you have a projection that takes a point $P$ on the surface
of the hemisphere to a projected point $P'$ on the plane such that $P$
and $P'$ are collinear with the center of the sphere.
(This is also called a gnomonic projection of the sphere.)
The coordinates of your "normals" are the coordinates of the
projected points.
As long as $0 \leq \phi < \frac\pi2$,
each point of your sphere does in fact project onto that plane,
although as $\phi$ approaches $\frac\pi2$ the magnitudes of your
"normals" grow without bound.
And obviously the boundary of your hemisphere, where $\phi = \frac\pi2$,
does not project on to the plane $z = 1$ at all.
If you change your choice of $f$ as recommended in the answer by H.R.,
this problem goes away. All your normals then will have the same magnitude,
and they will be defined at all points of the hemisphere.
When you give your "normals" in the form
$\langle\sin^2\phi\cos\theta, \sin^2\phi\sin\theta, \sin\phi\cos\phi\rangle$,
you again have non-uniform magnitudes, but this time the magnitudes
go to zero as $\phi$ approaches zero.
In effect, you are radially projecting the hemisphere
onto a kind of degenerate torus given by $r = \sin\phi$.
But notice that the three components of your "normals" all have the
common factor $\sin\phi$. You can normalize the magnitudes of all your
"normals" (except for the case $\phi = 0$)
by multiplying by the scalar $1/\sin\phi$.
If you do this, you get the vectors
$\langle\sin\phi\cos\theta, \sin\phi\sin\theta, \cos\phi\rangle$,
that is, the coordinates of each vector (for $0 < \phi \leq \frac\pi2$)
are simply the coordinates of the point on the sphere.
If you define these vectors as the normals for all points such
that $0 < \phi \leq \frac\pi2$, you may see that you can use the same
formula to define the normal for $\phi = 0$ as well, and it works very nicely.
The formula ${\bf r}_\phi \times {\bf r}_\theta$
works nicely for the measurement of area (which is what it is used for
in the example you quoted) precisely because its magnitude does go to zero
as $\phi$ goes to zero and does so in the same
way that the area element $r \,d\phi\,d\theta$ goes to zero.
But as you observed, this property is not so desirable when you are
trying to construct a set of normal vectors rather than trying to
integrate some scalar function over area.
You just need to solve your inequality for $r$:
\begin{align*}
a^2 &\geq 2 r^2 \cos^2(\theta) + 5 r^2 \sin^2(\theta) + 4 r^2 \sin(\theta) \cos(\theta)\\
&= r^2(2 \cos^2(\theta) + 5 \sin^2(\theta) + 4 \sin(\theta) \cos(\theta))\\
&= r^2(2 + 3 \sin^2(\theta) + 4 \sin(\theta) \cos(\theta)) \, .
\end{align*}
Then
$$
r^2 \leq \frac{a^2}{2 + 3 \sin^2(\theta) + 4 \sin(\theta) \cos(\theta)}
$$
(Note that the denominator is always $\geq 0$, so the inequality is preserved.) Thus your bounds of integration should be
\begin{align*}
&0 \leq \theta \leq 2 \pi\\
&0 \leq r \leq \frac{a}{\sqrt{2 + 3 \sin^2(\theta) + 4 \sin(\theta) \cos(\theta)}} \, .
\end{align*}
Best Answer
Your parametrization is valid.
The arc-length integral is not done by finding an antiderivative, but by numerical methods.$$ ||r'(\theta)|| =\sqrt{(-a\rho \sin{\theta})^2 + (b\rho \cos{\theta})^2} $$
Does not have a closed form antiderivative.
That is why we do not have a simple formula for the circumference of an ellipse.