Let's say that your surface is defined by $\vec{r}(\Theta) = [r_0 + s f(\Theta) ] \hat{r}$, where $\Theta$ represents the spherical angles. I'm going to conjecture that the normal vector in radial coordinates is proportional to $\vec{n} = (-s\nabla f,1)$, where the gradient is a 2-vector in the angle space, and the last component is the radial direction.
We can verify this by making sure the normal vector is orthogonal to the surface. Note that a tangent vector to the surface is $(1,s\nabla f)$. This tangent points along the gradient direction. A normal vector orthogonal to this one points along an isocontour, which by construction is $(\hat{\Theta}^\perp,0)$, where $\hat{\Theta}^\perp$ is 90 degree rotation of the unit angular vectors in the tangent plane to the sphere passing through the point. The dot product with both of these is zero; the first is obvious, while for the second, $\nabla f \cdot \hat{\Theta}^\perp = 0$.
TL;DR: In radial coordinates: $(-s\nabla f,1)$, then just normalize and convert to cartesian if needed.
Elaboration: Based on your comment, it seems that your $f(\vec{x})$ is a scalar function of a 3-vector. Then $\nabla f$ is a 3-vector. The notation I use lumps the spherical coordinates $(\phi,\theta)$ into a single abstract 2-vector-like quantity $\Theta$, so that when you say $f(\Theta)$, that means $f$ is a function of only the spherical coordinate angles (and not a function of radius). Then $\nabla f(\Theta)=\nabla_\Theta f(\Theta)$ is a gradient in this 2-dimensional angle space. Your $\nabla f$ needs to be projected onto the surface of a sphere first, since you only ever sample $f(\vec{x})$ for $\lVert \vec{x} \rVert = 1$.
So, to calculate your normals, given a point $\vec{P}(\vec{x}) = [R + s f(\vec{x})]\vec{x}$ with $\vec{x}$ such that $\lVert \vec{x} \rVert = 1$,
- Let $\vec{g}(\vec{x}) = \nabla f(\vec{x})$. Here, $\vec{g}$ is the true gradient of $f$ in 3-space, and you can calculate that in cartesian coordinates, so you get $\vec{g} = (g_x,g_y,g_z)$.
- Project out the radial component of $\vec{g}$ to get $\vec{h}(\vec{x})$. To do this, $\vec{h} = \vec{g} - \frac{\vec{g}\cdot\vec{x}}{\vec{x}\cdot\vec{x}}\vec{x}$. Note that the denominator should be 1. Here, $\vec{h}$ represents the component of $\vec{g}$ that should be tangential to a sphere centered at the origin and passing through $\vec{P}$.
- An outward normal vector to the surface is $\vec{n} = \vec{x} - s\cdot \vec{h}$. Normalize it to get a unit normal vector.
All these computations can be done in Cartesian coordinates, but notice that I never had to resort to referring to the Cartesian components of any vector; the basic operations are vector arithmetic and dot products.
Elaborating on my comment: This webpage describes a good method to interpolate and find the barycentric coordinates of a point inside a triangle, which can easily be used to calculate the 'weight' or color of the point as asked by the OP. The method involves first representing the triangle as a triplet of vectors $(\mathbf{a} = (a_x, a_y), \mathbf{b} = (b_x, b_y), \mathbf{c} = (c_x, c_y))$ and the point as another vector $\mathbf{p} = (p_x, p_y)$. (The linked webpage uses 3-dimensional vectors but the method is generalizable.)
Once the vectors are obtained, we can calculate the areas of the triangles formed by joining the vertices of the triangle to the point inside it. The area of the whole triangle is
$$\Delta = \frac{1}{2} \cdot ||\mathbf{a} - \mathbf{b}|| \cdot ||\mathbf{a} - \mathbf{c}||$$
(which is half the magnitude of the cross product of the two vectors). The areas of the triangles opposite vertices $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$, respectively, are
$$
\Delta_a = \frac{1}{2} \cdot ||\mathbf{p} - \mathbf{b}|| \cdot ||\mathbf{p} - \mathbf{c}|| \\
\Delta_b = \frac{1}{2} \cdot ||\mathbf{p} - \mathbf{a}|| \cdot ||\mathbf{p} - \mathbf{c}|| \\
\Delta_c = \frac{1}{2} \cdot ||\mathbf{p} - \mathbf{a}|| \cdot ||\mathbf{p} - \mathbf{b}||
$$
Now, we calculate the 'interpolation factors' of the point $\mathbf{p}$ as $k_a = \frac{\Delta_a}{\Delta}$, $k_b = \frac{\Delta_b}{\Delta}$, and $k_c = \frac{\Delta_c}{\Delta}$. The barycentric coordinates of $\mathbf{p}$ are finally given by
$$\mathbf{p} = \mathbf{a}k_a + \mathbf{b}k_b + \mathbf{c}k_c$$
(In the linked page, they use 'UV' coordinates, but since we are using 2 dimensions throughout, our 'UV' coordiantes are the same as our 'XY[Z]' coordinates.)
If $w_a$, $w_b$, $w_c$ denote the weights of, or colors of, or amount of cheese contained in, the points $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$, respectively, then $w_p = w_ak_a + w_bk_b + w_ck_c$ is the corresponding quantity for $\mathbf{p}$, assuming the 'linear' interpolation.
Extra: It would be interesting if this technique could be generalized to any polygon, using the areas opposite the vertices or something similar. I don't know much about this, but perhaps someone more knowledgeable can comment.
Edit: The answer to the above is very simple, just use a triangle within the polygon that contains the point. I overlooked that :)
Best Answer
You seem to be a little confused. If you know the direction of the normal, then you can determine what the order of the vertices should be. If you know the ordering of the vertices, then you can determine the direction of the normal. But if you know neither, you cannot determine either. And if you know both, what's the question?
The easy way to keep track of everything is this. When at first you are defining the triangles for your (triangulated) cube, you know which direction each face's normal should point in, so you should order the vertices of your triangles accordingly. Later, when you just have a bunch of triangles and have forgotten that they came from a cube, you can recover the direction of the normal because you know the vertices are in a consistent orientation. That's all there is to it. Leave a comment if anything's unclear.
If you have a set of vertices $\{p, q, r\}$ and you want the normal to be in the direction of some known vector $n$, then you can start by assuming the order is $(p,q,r)$, computing the corresponding normal $\tilde n = (q-p)\times(r-p)$, and checking whether the dot product $\tilde n\cdot n$ is positive; if not, take the opposite order $(r,q,p)$ instead.