Beside the projection to the 3-sphere and the "magic" formula the process is exacly the same as for 3D-polyhedron.
As you said, the 3d point is projected on the unit sphere using inverse stereographic projection.
The maths behind the "magic" formula is quite simple (see picture below). It's indeed used to convert distance on the sphere (3-sphere) into distance on the plane (3D-space). This conversion is essentially a 2D problem:
The given point Z is projected onto the unit sphere giving Z'. Say, "alpha" is the angle between Z' and the x-axis.
Say we know the spherical distance beween Z' and
a point on the sphere ("delta" in the figure). The spherical distance is simply the angle beween Z' and that point.
Now we construct a point P' on the sphere which angle is "alpha"
+"delta". Then we project P' back on the plane and get P.
The distance between P and Z is the distance estimate.
The cosine of the angle between Z' and a vertex V of the polychora is simply the dot product of Z' and V. For the segments we have to find the closesest point to Z', project it onto the sphere (by normalizing it) then compute the dot product. the sine is a little more involved because we don't have vector product in 4D. Using sqrt(1-cos²(delta)) is not accurate enought.
The renderings in wikipedia of polychora with straight lines is maybe done with the classic
method. having the set of vertices and segments. Project the vertices onto 3D space using some appropriate projection (ideally one that shows most of the symmetries) Then render it as usuall. As I said in a previous post it is possible to use DE but it would be very slow (that also may answer kram's question). It's applicable to any well behaving 4D distance field and is just like raymarching... but it's a double raymarching:
From the point of 3D space shoot a ray into 4D (much exactly what we do in 3D but in 4D: the screen space becomes a box
). Find the minimal distance along that ray. This is the distance we will use to marche in 3D.
AFAIK, the boringness
of quaternionic fractals comes from the so called j,k equivalence, not from the projection. In particular quaternionic madelbrot set is simply a rotated 2D madelbrot set around the real axis in j and k directions.