Discrete K-surfaces tutorial

This tutorial explains how one may implement the discrete K-surfaces application with jReality. This tutorial is one of the longer tutorials which explains a small but still reasonably complete application. The formulas behind this application are explained in Ulrich Pinkall: Designing Cylinders with Constant Negative Curvature.

Vector operations in <math>R^3</math>

For the vector calculus in <math>R^3</math> we represent vectors as arrays of doubles of length 3 and implement operations on these vectors as static methods of a class named R3:

 
/** This class provides simple static methods to do verctor calculation in R<sup>3</sup>. 
 * The vectors of R<sup>3</sup> are represented by array of 3 doubles.
 *
 */
public class R3 {
	/** target = a + b
	 * 
	 * @param a array of 3 doubles
	 * @param b array of 3 doubles
	 * @param target array of 3 doubles
	 */
	public final static void plus(double[] a, double[] b, double[] target) {
		target[0] = a[0] + b[0];
		target[1] = a[1] + b[1];
		target[2] = a[2] + b[2];
	}
	
	/** target = a-b
	 * 
	 * @param a array of 3 doubles
	 * @param b array of 3 doubles
	 * @param target array of 3 doubles
	 */
	public final static void minus(double[] a, double[] b, double[] target) {
		target[0] = a[0] - b[0];
		target[1] = a[1] - b[1];
		target[2] = a[2] - b[2];
	}
 
	/** target = a x b, i.e., the <a href="http://en.wikipedia.org/wiki/Crossproduct">cross product or vector product</a><br>
	 * 
	 * WARNING: This method fails when the target is same as one of its arguments. 
	 *
	 * @param a array of 3 doubles
	 * @param b array of 3 doubles
	 * @param target array of 3 doubles
	 */
	public final static void cross(double[] a, double[] b, double[] target) {
		double x = -a[1]*b[2] + a[2]*b[1];
		double y = -a[2]*b[0] + a[0]*b[2];
		target[2] = -a[0]*b[1] + a[1]*b[0];
		target[0] = x;
		target[1] = y;
	}
 
	/** target = a, copy the entries of a into target
	 *  
	 * @param a array of 3 doubles
	 * @param target array of 3 doubles
	 */
	public final static void copy(double[] a, double[] target) {
		target[0] = a[0];
		target[1] = a[1];
		target[2] = a[2];
	}
 
	/** target = 0, reset the three entries of target to 0
	 *  
	 * @param target array of 3 doubles
	 */
	public final static void zero(double[] target) {
		target[0] = 0;
		target[1] = 0;
		target[2] = 0;
	}
	
	/** &lt;a,b&gt;, i.e., the <a href="http://en.wikipedia.org/wiki/Dotproduct">dot product or inner product</a>
	 * 
	 * @param a array of 3 doubles
	 * @param b array of 3 doubles
	 * @return the value of the dot product
	 */
	public final static double dot(double[] a, double[] b) {
		return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
	}
	
	/** target = s * a, product of a real with a vector
	 * 
	 * @param s the scalar real 
	 * @param a array of 3 doubles
	 * @param target array of 3 doubles
	 */
	public final static void times(double s, double[]a, double[]target) {
		target[0] = s * a[0];
		target[1] = s * a[1];
		target[2] = s * a[2];
	}
 
	/** |a|, the length of a vector.
	 * 
	 * @param a array of 3 doubles
	 * @return the length of a 
	 */
	public static double norm(double[] a) {
		return Math.sqrt(normSquared(a));
	}
 
	/** |a|^2, the square of the length of a vector.
	 * 
	 * @param a array of 3 doubles
	 * @return the square of the length of a
	 */
	public static double normSquared(double[] a) {
		return a[0]*a[0]+a[1]*a[1]+a[2]*a[2];
	}
}

Calculate the Gauss map

The Gaus map of a discrete K-surface consists of spherical parallelograms, i.e., spherical quadrilateral that are invariant under a 180° rotation. But the parallelograms are skew in our coordinates (this corresponds to the coordinate change <math>u=x+t,\ v=x-t</math> which transforms the form <math>N_{xx}-N_{tt}=0</math> of the wave equation into <math>N_{uv}=0</math>. The equation <math>f_{uv}=0</math>), i.e., a parallelogram is formed by the Gauss map at (i-2,j-1), (i-1,j-1), (i,j), and (i-1,j). Thus the Gauss map at (i,j) is 180° rotation of the Gauss map at (i-2,j-1) about the mid point of the diagonal between (i-1,j-1) and (i-1,j).

Let us implement the algorithm for computing the Gauss map as a static method of the class KSurfaces:

public class KSurfaces {
	/** Computes the Gauss map of a K-surface from initial Cauchy data, i.e., two closed curves. Assums that the points of
	* the initial data lie on the 2-sphere this method generates more points so that the quadrilaterals 
	* (i-1,j-1), (i-1,j), (i,j),(i-2,j-1) form spherical parallelograms. 
	* 
	* @param initialAnnulus the initial data, i.e., a double array double[2][n][3] containing 2 polygons with n vertices,
	*  which are both interpreted as closed curves connecting that last and the first one. 
	* @param target a double array that may be filled with the result, i.e., an array double[m][n][3], where m>1 is the number of
	* 	time steps to be calculated.
	*/
	public static void gaussMapFromInitialAnnulus(double[][][] initialAnnulus, double[][][] target) {
		final int m = target.length;
		final int n = target[0].length;
		final double[] a = new double[3];
	
		// copy first two rows
		for (int i=0; i<2; i++) {
			for (int j=0; j<n; j++) {
				R3.copy(initialAnnulus[i][j], target[i][j]);
			}
		}
		
		// compute the other rows
		for (int i=2; i<m; i++) {
			for (int j=0; j<n; j++) {
				int k = j==0 ? n-1 : j-1;
				R3.plus(target[i-1][k], target[i-1][j], a);
				double[] w = target[i-2][k];
				double s = 2 * R3.dot(a,w) / R3.dot(a,a);
				R3.times(s, a, a);
				R3.minus(a, w, target[i][j]);
			}
		}
	}
}

Display the Gauss map

Now we are able to display the Gauss map. For the moment, we choose an ad hoc initial condition: the intersection of two elliptic cones around the z-achses with the 2-sphere. Lets put the code for this in yet another class named KSurfacesApp.

The fields of KSurfacesApp

import de.jreality.geometry.IndexedFaceSetFactory;
import de.jreality.plugin.JRViewer;
import de.jreality.plugin.JRViewer.ContentType;
import de.jreality.plugin.basic.ViewShrinkPanelPlugin;
import de.jreality.plugin.content.ContentAppearance;
import de.jreality.plugin.content.ContentTools;
import de.jreality.scene.SceneGraphComponent;
import de.jtem.beans.InspectorPanel;
import de.varylab.jrworkspace.plugin.PluginInfo;
 
/** This class contains the application to display the K-surfaces and their Gauss map.
 */
public class KSurfacesApp {
	private int m=3; 
	private int n=30; 
	private double a1=1; 
	private double a2=.9;
	private double b1=1.2;
	private double b2=1.1;
	private double[][][] gaussMap;
	private double[][][] initialGaussMap;
	
	final private SceneGraphComponent root=new SceneGraphComponent("k-surfaces");
	final private SceneGraphComponent gaussMapSgc=new SceneGraphComponent("Gauss map");
	final private IndexedFaceSetFactory gaussMapIfs = new IndexedFaceSetFactory();

The constructor of KSurfacesApp

public KSurfacesApp() {
	gaussMapIfs.setGenerateEdgesFromFaces(true);
	gaussMapIfs.setGenerateFaceNormals(true);
	gaussMapSgc.setGeometry(gaussMapIfs.getGeometry());
	root.addChild(gaussMapSgc);
 
	update();
}

The method that generates the initial condition

This method generates two generations of normal vectors for use in the method gaussMapFromInitialAnnulus

private void ellipticConesInitialCondition() { 
	initialGaussMap = new double[2][n][3];
	for (int j=0; j<n; j++){
		double s=a1*Math.sin(2*Math.PI/n*j);
		double c=a2*Math.cos(2*Math.PI/n*j);
		double f=1/Math.sqrt(s*s+c*c+1);
		initialGaussMap[0][j][0]=s*f;
		initialGaussMap[0][j][1]=c*f;
		initialGaussMap[0][j][2]=1*f;
 
		s=b1*Math.sin(2*Math.PI/n*(j-.5));
		c=b2*Math.cos(2*Math.PI/n*(j-.5));
		f=1/Math.sqrt(s*s+c*c+1);
		initialGaussMap[1][j][0]=s*f;
		initialGaussMap[1][j][1]=c*f;
		initialGaussMap[1][j][2]=1*f;
	}
}

Recalculate the gauss map and update the display

private void gaussMap() {
	//calculate the gauss map
	gaussMap=new double[m][n][3];
	KSurfaces.gaussMapFromInitialAnnulus(initialGaussMap, gaussMap);
	
	//rearrange the vertices for use in an IndexedFaceSetFactory
	double[][] vertices=new double[n*m][3];
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			for (int k = 0; k < 3; k++) { 
				vertices[j+i*n][k]=gaussMap[i][j][k];
			}
		}
	}
	
	//calculate the combinatorics of the faces for use in an IndexedFaceSetFactory
	int[][] faces = new int[n*(m-2)][4];
	for (int i = 0; i < m-2; i++) {
		for (int j = 0; j < n; j++) {
			int k = j + i*n;
			faces[k][0] = j       + i*n;
			faces[k][1] = j       + (i+1)*n;
			faces[k][2] = (j+1)%n + (i+2)*n;
			faces[k][3] = (j+1)%n + (i+1)*n;
		}
	}
	//put the data into the IndexedFaceSetFactory ifs and update the Geometry
	gaussMapIfs.setVertexCount(n*m);
	gaussMapIfs.setVertexCoordinates(vertices);
	gaussMapIfs.setFaceCount(n*(m-2));
	gaussMapIfs.setFaceIndices(faces);
	gaussMapIfs.update();
}

Getters and Setters

 
	public int getM() { return m; }
	public void setM(int m) {
		if (m<2) return;
		this.m = m;
	}
 
	public int getN() { return n; }
	public void setN(int n) {
		if (n<2) return;
		this.n = n;
	}
 
	public double getA1() {	return a1; }
	public void setA1(double a1) {	this.a1 = a1; 	}
 
	public double getA2() { return a2; }
	public void setA2(double a2) { this.a2 = a2; }
 
	public double getB1() {	return b1; }
	public void setB1(double b1) {	this.b1 = b1; }
 
	public double getB2() {	return b2; }
	public void setB2(double b2) { this.b2 = b2; }

/**
* @return the root of the scene graph with the Gauss map of the k-surface
 */
public SceneGraphComponent getRoot() {
	return root;
}

Update method

/** This method needs to be called whenever the changes of the parameters should have an effect on the
* scene graph one gets from {@link #getRoot()}. 
*/
public void update() {
	ellipticConesInitialCondition();
	gaussMap();
}

The main method

public static void main(String[] args) {
	KSurfacesApp ksa = new KSurfacesApp();
	
	//initialize the ViewerVR
	JRViewer v = new JRViewer();
	v.addBasicUI();
	v.addVRSupport();
	v.addContentSupport(ContentType.TerrainAligned);
	v.registerPlugin(new ContentAppearance());
	v.registerPlugin(new ContentTools());
	v.setContent(ksa.getRoot());
 
	//add an Inspector for the KSurfacesApp to the viewer
	InspectorPanel inspector = new InspectorPanel();
	inspector.setObject(ksa, "update");
	ViewShrinkPanelPlugin plugin = new ViewShrinkPanelPlugin() {
		public PluginInfo getPluginInfo() {
			return new PluginInfo("Domain");
		}
	};
	plugin.getShrinkPanel().add(inspector);
	v.registerPlugin(plugin);
 
	v.startup();
}

Editor for the initial curve and velocity

Have a look at the Polygon editor tutorial and use this for the KSurfacesApp class.

Here is the essential code you may add to the constructor

//initial condition at time 1 draggable points
a1=0;a2=0;
b1=.2;b2=.2;
ellipticConesInitialCondition();
//to keep the points on the unit sphre override the transform method
dragPointSet=new DragPointSet(initialGaussMap[1]) {
	@Override
	public void transform(double[] vertex) {
		double l=R3.norm(vertex);
		vertex[0]=vertex[0]/l;
		vertex[1]=vertex[1]/l;
		vertex[2]=vertex[2]/l;
	}
};
	dragPointSet.getBase().getAppearance().setAttribute(CommonAttributes.POINT_SHADER+"."+CommonAttributes.POINT_RADIUS, .16);
gaussMapSgc.addChild(dragPointSet.getBase());
dragPointSet.addChangeListener(new ChangeListener() {
	public void stateChanged(ChangeEvent e) {
		KSurfacesApp.this.update();
	}
});
		
 
//initial condition at time 1 subdivided Polygon points
time1SubdivPoly = new SubdividedPolygon(dragPointSet);
time1SubdivPoly.setSubdivisionLevel(0);
PointSequenceView time1CurveView = new PointSequenceView(time1SubdivPoly);
time1CurveView.setPointRadius(0.01);
time1CurveView.setLineRadius(0.04);
time1CurveView.setPointColor(Color.yellow);
time1CurveView.setLineColor(Color.green);
gaussMapSgc.addChild(time1CurveView.getBase());

And two lines two the update() method

n=time1SubdivPoly.getPoints().length;
ellipticConesInitialCondition();
initialGaussMap[1]=time1SubdivPoly.getPoints();

Integrate the Gauss map to get the surface

The first row, i.e., to calculate f(0,j) we use the formula <math>f_{\tilde d}=f_d +(N_d+N_{\tilde d})\times N_{r}</math>, where in our xt-coordinates <math>d=(0,j),\ \tilde d=(0,j+1),\ r=(1,j+1)</math>. For the other rows we use <math>f_r=f_d+N_d\times N_r</math>, where <math>d=(i-1,j),\ r=(i,j)</math>. Let us put this algorithm as a static method into the class KSurfaces:

	/** Calculates the K-surface from the given Gauss map, assuming that the Gauss map consists
	 * of spherical parallelograms as described in {@link #gaussMapFromInitialAnnulus(double[][][], double[][][])}. 
	 * 
	 * @param gaussMap the given Gauss map.
	 * @param target a double array with enough space to hold the resulting surface. A quad mesh where the
	 * quadrilaterals consist of (i-1,j-1), (i-1,j), (i,j),(i-2,j-1).
	 */
	public static void kSurfaceFromGaussMap(double[][][] gaussMap, double[][][] target) {
		final int m = gaussMap.length;
		final int n = gaussMap[0].length;
		final double[] a = new double[3];
 
		// compute first row
		R3.zero(target[0][0]);
		for (int j=0; j<n-1; j++) {
			int k = (j+1)%n;
			R3.plus(gaussMap[0][j], gaussMap[0][k], a);
			R3.cross(gaussMap[1][k], a, a);
			R3.plus(target[0][j], a, target[0][k]);
		}
			
		// compute the other rows
		for (int i=1; i<m; i++) {
			for (int j=0; j<n; j++) {
				R3.cross(gaussMap[i-1][j], gaussMap[i][j], a);
				R3.plus(target[i-1][j], a, target[i][j]);
			}
		}
	}

Display the surface

We add to the class KSurfacesApp two new imports

import de.jreality.math.MatrixBuilder;
import de.jreality.scene.Transformation; 

Then a field for the map, its scene graph component and its factory:

private double[][][] map;
final private SceneGraphComponent mapSgc=new SceneGraphComponent("The K-surface");
final private IndexedFaceSetFactory mapIfs = new IndexedFaceSetFactory();

Then the following code to the constructor KSurfacesApp()

mapIfs.setGenerateEdgesFromFaces(true);
mapIfs.setGenerateFaceNormals(true);
mapSgc.setGeometry(mapIfs.getGeometry());
mapSgc.setTransformation(new Transformation(
				MatrixBuilder.euclidean().translate(3,0,0).getArray()));
	root.addChild(mapSgc);

The we refactor the methode gaussMap() to avoid code doubling. Just replace the old version by

//calculate the gauss map and update it via gaussMapIfs
private void gaussMap() {
	//calculate the gauss map
	gaussMap=new double[m][n][3];
	KSurfaces.gaussMapFromInitialAnnulus(initialGaussMap, gaussMap);
	
	double[][] vertices = rearrangeVertices(gaussMap);
	int[][] faces = calcFaces();
 
	//put the data into the IndexedFaceSetFactory and update the Geometry
	gaussMapIfs.setVertexCount(n*m);
	gaussMapIfs.setVertexCoordinates(vertices);
	gaussMapIfs.setFaceCount(n*(m-2));
	gaussMapIfs.setFaceIndices(faces);
	gaussMapIfs.update();
	}
 
	//calculate the map of K-surface and update it via gaussMapIfs
	private void map() {
	//calculate the map
	map=new double[m][n][3];
	KSurfaces.kSurfaceFromGaussMap(gaussMap, map);
		
	double[][] vertices = rearrangeVertices(map);
	int[][] faces = calcFaces();
 
	//put the data into the IndexedFaceSetFactory and update the Geometry
	mapIfs.setVertexCount(n*m);
	mapIfs.setVertexCoordinates(vertices);
	mapIfs.setFaceCount(n*(m-2));
	mapIfs.setFaceIndices(faces);
	mapIfs.update();
}
 
//calculate the combinatorics of the faces for use in an IndexedFaceSetFactory
private int[][] calcFaces() {
	int[][] faces = new int[n*(m-2)][4];
	for (int i = 0; i < m-2; i++) {
		for (int j = 0; j < n; j++) {
			int k = j + i*n;
			faces[k][0] = j       + i*n;
			faces[k][1] = j       + (i+1)*n;
			faces[k][2] = (j+1)%n + (i+2)*n;
			faces[k][3] = (j+1)%n + (i+1)*n;
		}
	}
	return faces;
}
 
//rearrange the vertices for use in an IndexedFaceSetFactory
private double[][] rearrangeVertices(double[][][] xtVertices) {
	double[][] vertices=new double[n*m][3];
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			for (int k = 0; k < 3; k++) { 
				vertices[j+i*n][k]=xtVertices[i][j][k];
			}
		}
	}
	return vertices;
}

Finally call the method map() from update()

map();