Rendering planets accurately
Table of Contents

Rendering planets accurately

Today I have finished the revision of planetary rendering in JPARSEC. Initial implementation was quite good, but not fully accurate. After improving performance it was clear that a redo of some parts was necessary. Rotations around position angle of axis, bright limb angle, and phase angle were mixed everywhere without a good organization. Now I have implemented two core functions to rotate positions as seen from the Sun to positions as seen from the Earth and the opposite, that are used instead of the previous algorithms to bring more clean and accurate code. The functions are:

double[] fromEarthTofromSun(double[] xyz, double brightLimbAngle, double phaseAngle, double positionAngleOfAxis) {
	double px = xyz[0], py = xyz[1], pz = xyz[2];
	double pr = Math.sqrt(px*px+py*py);
	double pang = FastMath.atan2(-px, py);
 
	// From the Earth: derotate along the position angle of axis
	px = pr * FastMath.sin(pang + positionAngleOfAxis);
	py = pr * FastMath.cos(pang + positionAngleOfAxis);
	pang = FastMath.atan2(-px, py);
 
	// From the Sun: rotate to the bright limb angle, rotate along the
	// phase angle, and back
	px = pr * FastMath.sin(pang + brightLimbAngle);
	py = pr * FastMath.cos(pang + brightLimbAngle);
 
	double r = Math.sqrt(py * py + pz * pz);
	double ang = FastMath.atan2(pz, py) - Math.abs(phaseAngle);
	double x = -px;
	double y = r * FastMath.cos(ang);
	double z = r * FastMath.sin(ang);
 
	ang = positionAngleOfAxis - brightLimbAngle;
	pang = FastMath.atan2(-x, y);
	pr = Math.sqrt(x * x + y * y);
	px = pr * FastMath.sin(pang + ang);
	py = pr * FastMath.cos(pang + ang);
	pz = z;
 
	return new double[] {px, py, pz};
}
double[] fromSunTofromEarth(double[] xyz, double brightLimbAngle, double phaseAngle, double positionAngleOfAxis) {
	double px = xyz[0], py = xyz[1], pz = xyz[2];
 
	double ang = positionAngleOfAxis - brightLimbAngle;
	double pang = -FastMath.atan2(-px, py);
	double pr = Math.sqrt(px * px + py * py);
	px = pr * FastMath.sin(pang - ang);
	py = pr * FastMath.cos(pang - ang);
 
	pr = Math.sqrt(py * py + pz * pz);
	ang = (FastMath.atan2(pz, py) + Math.abs(phaseAngle));
	px = -px;
	py = pr * FastMath.cos(ang);
	pz = pr * FastMath.sin(ang);
 
	pr = Math.sqrt(px*px+py*py);
	pang = FastMath.atan2(-px, py);
 
	// From the Sun: rotate to the bright limb angle, rotate along the
	// phase angle, and back
	px = pr * FastMath.sin(pang - brightLimbAngle);
	py = pr * FastMath.cos(pang - brightLimbAngle);
 
	pang = -FastMath.atan2(-px, py);
 
	// From the Earth: rerotate along the position angle of axis
	px = pr * FastMath.sin(pang + positionAngleOfAxis);
	py = pr * FastMath.cos(pang + positionAngleOfAxis);
 
	return new double[] {px, py, pz};
}

An application of the previous code to all operations in planetary rendering, not only to draw planets and rings themselves, but also to project shadows on the surface of the planets gives very nice results, both visually and from the implementation point of view. The new algorithm to draw shadows on the planet is here.

void drawOvalShadow(double posx, double sx, double posy, double sy, MoonEphemElement m, EphemElement e, double incl_up, double scale, double dr, int pixel, int[][] screen_out,
		double sun_dr)
{		
	if (isInTheScreen((int) (posx + sx + 0.5), (int) (posy + sy + 0.5), (int) (4 * dr + 1)))
	{
		// Get position of the center of the shadow respect to planet axis
		double salfa = FastMath.atan2(sy, sx);
		double srr = Math.sqrt(sx*sx+sy*sy);
		double sdx0 = srr * FastMath.cos(salfa + incl_up);
		double sdy0 = srr * FastMath.sin(salfa + incl_up);
 
		double sampling = (int) (dr + sun_dr * 2); // size of shadow + penumbra
		double step = 20.0/scale; // step when sampling points according to planet size in pixels
		if (step > 1) step = 1;
		// test points
		int[][] screenCopy = clone(screen_out);
		for (double i = -sampling; i <= sampling; i = i + step)
		{
			double ssx = sx + i;
			for (double j = -sampling; j <= sampling; j = j + step)
			{
				double ssy = (sy + j);
 
				// Get position around the center of the shadow respect to planet axis
				salfa = FastMath.atan2(ssy, ssx);
				srr = Math.sqrt(ssx*ssx+ssy*ssy);
				double sdx = (srr * FastMath.cos(salfa + incl_up) - sdx0) / scale;
				double sdy = (srr * FastMath.sin(salfa + incl_up) - sdy0) / scale;
				double sdr = Math.sqrt(Math.pow(sdx, 2.0)+Math.pow(sdy, 2.0))*scale;
 
				// if we are inside shadow + penumbra ...
				if (sdr <= sampling) {
					// OK, now get position of satellite respect to sun and project it onto
					// planet surface, just calculating the adequate value of z = cos beta (r = 1 on planet).
					double x = m.xPositionFromSun+sdx;
					double y = -m.yPositionFromSun+sdy;
					// not using x^2+y^2+z^2 = 1 neither obtaining cosb from y since y depends on oblateness...
					double a = FastMath.atan2(y, x);
					double sinb = x / FastMath.cos(a);
					double cosb = Math.sqrt(1.0 - sinb * sinb) * FastMath.sign(m.zPositionFromSun);
 
					// Now pass the vision from Sun to from Earth to see where is this point around the shadow center.
					// From Sun the apparent shape is round, but from Earth not necessarily.
					double xyz[] = this.fromSunTofromEarth(new double[] {x, y, cosb}, e.brightLimbAngle, e.phaseAngle, e.positionAngleOfAxis);
 
					// Rotate back again from planet axis to current orientation, and set position in pixels
					srr = Math.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1])*scale;
					salfa = FastMath.atan2(xyz[1], xyz[0]);
					sdx = srr * FastMath.cos(salfa - incl_up);
					sdy = srr * FastMath.sin(salfa - incl_up);
 
					// Get shadow position
					int shadowX = (int) (posx + sdx + 0.5);
					int shadowY = (int) (posy + sdy + 0.5);
 
					// After all this fun, now simply draw shadow and penumbra
					if (screenCopy[shadowX][shadowY] == screen_out[shadowX][shadowY]) { // dont process a pixel twice
						if (sdr<dr) {
							screen_out[shadowX][shadowY] = pixel;						
							if (texture_step == 1) { // More points to avoid gaps when rendering to a big size
								screen_out[shadowX-1][shadowY] = pixel;						
								screen_out[shadowX][shadowY-1] = pixel;						
								screen_out[shadowX-1][shadowY-1] = pixel;						
							}
						} else {
							// Reduce intensity from penumbra to edge softly
							Color col = new Color(screen_out[shadowX][shadowY]);
							float penumbraIntensity = (0.3f + 0.7f * (float) ((sdr-dr)/(sampling-dr))) / 255f;
							int ppixel = new Color(penumbraIntensity*col.getRed(), penumbraIntensity*col.getGreen(), penumbraIntensity*col.getBlue()).getRGB();
							screen_out[shadowX][shadowY] = ppixel;
						}
					}
				}					
			}
		}
	}
}

Some comments are given in the code, although I have no time to explain everything here. Anyway, the result after applying all that to, for example, the triple eclipse on Jupiter on March 28, 2004, is:

 

At the right is the same image taken with HST (APOD November 11, 2004, http://apod.nasa.gov/apod/ap041111.html). I don't know if it was taken exactly at the same time, anyway it is hard to see any difference. All maths are done using FastMath class, so some lose of precision exists. Anyway, is a nice result taking into account that everything is done pixel by pixel, without any 3d library. I'm that brute.

The topic covered in this post was revised and extended here.

 
blog/mail_name_date_rendering_planets.txt · created: 2010/02/11 20:18 (Last modified 2016/05/17 12:51) by Tom├ís Alonso Albi
 
Recent changes RSS feed Creative Commons License Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki