I am wondering how to generate uniformly distributed points on the surface of the 3-d unit sphere? Also after generating those points, what is the best way to visualize and check whether they are truly uniform on the surface $x^2+y^2+z^2=1$?
A standard method is to generate three standard normals and construct a unit vector from them.  That is, when $X_i \sim N(0,1)$ and $\lambda^2 = X_1^2 + X_2^2 + X_3^2$, then $(X_1/\lambda, X_2/\lambda, X_3/\lambda)$ is uniformly distributed on the sphere.  This method works well for $d$-dimensional spheres, too.

In 3D you can use rejection sampling: draw $X_i$ from a uniform$[-1,1]$ distribution until the length of $(X_1, X_2, X_3)$ is less than or equal to 1, then--just as with the preceding method--normalize the vector to unit length.  The expected number of trials per spherical point equals $2^3/(4 \pi / 3)$ = 1.91.  In higher dimensions the expected number of trials gets so large this rapidly becomes impracticable.

There are many ways to check uniformity.  A neat way, although somewhat computationally intensive, is with Ripley's K function.  The expected number of points within (3D Euclidean) distance $\rho$ of any location on the sphere is proportional to the area of the sphere within distance $\rho$, which equals $\pi\rho^2$.  By computing all interpoint distances you can compare the data to this ideal.

General principles of constructing statistical graphics suggest a good way to make the comparison is to plot variance-stabilized residuals $e_i(d_{[i]} - e_i)$ against $i = 1, 2, \ldots, n(n-1)/2=m$ where $d_{[i]}$ is the $i^\text{th}$ smallest of the mutual distances and $e_i = 2\sqrt{i/m}$.  The plot should be close to zero.  (This approach is unconventional.)

Here is a picture of 100 independent draws from a uniform spherical distribution obtained with the first method:



Here is the diagnostic plot of the distances:



The y scale suggests these values are all close to zero.

Here is the accumulation of 100 such plots to suggest what size deviations might actually be significant indicators of non-uniformity:



(These plots look an awful lot like Brownian bridges...there may be some interesting theoretical discoveries lurking here.)

Finally, here is the diagnostic plot for a set of 100 uniform random points plus another 41 points uniformly distributed in the upper hemisphere only:



Relative to the uniform distribution, it shows a significant decrease in average interpoint distances out to a range of one hemisphere.  That in itself is meaningless, but the useful information here is that something is non-uniform on the scale of one hemisphere.  In effect, this plot readily detects that one hemisphere has a different density than the other.  (A simpler chi-square test would do this with more power if you knew in advance which hemisphere to test out of the infinitely many possible ones.)
