fixed curvature

src/main/java/bdv/jogl/VolumeRenderer/ShaderPrograms/ShaderSources/functions/accumulator/MaxCurvatureDifference.java

@@ -3,11 +3,26 @@ package bdv.jogl.VolumeRenderer.ShaderPrograms.ShaderSources.functions.accumulat
 import static bdv.jogl.VolumeRenderer.ShaderPrograms.ShaderSources.MultiVolumeRendererShaderSource.scvMaxNumberOfVolumes;
 import static bdv.jogl.VolumeRenderer.ShaderPrograms.ShaderSources.MultiVolumeRendererShaderSource.sgvRayPositions;
 import static bdv.jogl.VolumeRenderer.ShaderPrograms.ShaderSources.MultiVolumeRendererShaderSource.suvVolumeTexture;
+import static bdv.jogl.VolumeRenderer.ShaderPrograms.ShaderSources.MultiVolumeRendererShaderSource.sgvVolumeNormalizeFactor;
 import static bdv.jogl.VolumeRenderer.utils.ShaderSourceUtil.addCodeArrayToList;
 import static bdv.jogl.VolumeRenderer.utils.ShaderSourceUtil.appendNewLines;
+import static bdv.jogl.VolumeRenderer.utils.VolumeDataUtils.calculateCurvatureOfVolume;
+import static bdv.jogl.VolumeRenderer.utils.VolumeDataUtils.createVolumeTexture;
 
+import java.nio.FloatBuffer;
 import java.util.ArrayList;
+import java.util.HashMap;
 import java.util.List;
+import java.util.Map;
+
+import bdv.jogl.VolumeRenderer.Scene.Texture;
+import bdv.jogl.VolumeRenderer.ShaderPrograms.MultiVolumeRenderer;
+import bdv.jogl.VolumeRenderer.utils.CurvatureContainer;
+import bdv.jogl.VolumeRenderer.utils.VolumeDataBlock;
+import bdv.jogl.VolumeRenderer.utils.VolumeDataManager;
+
+import com.jogamp.common.nio.Buffers;
+import com.jogamp.opengl.GL4;
 
 /**
  * Calculates the maximum difference of the curvature values
@@ -19,123 +34,111 @@ public class MaxCurvatureDifference extends AbstractVolumeAccumulator {
 		super("curvature_difference");
 	}
 	
+	private boolean needsReset = true;
+	
+	private final Map<Integer, Texture> curvatureTexture = new HashMap<Integer, Texture>();
+	
+	private final Map<Integer, CurvatureContainer> evaluatedCurvatures = new HashMap<Integer, CurvatureContainer>();
+	
+	private final static String suvCurvatureTexture = "inCurvatureTexture"; 
+	
+	private final static String suvCurvatureMax = "inCurvatureMax";
+
+	private final static String suvCurvatureMin = "inCurvatureMin";
+	@Override
+	public void init(GL4 gl) {
+		getParent().mapUniforms(gl,new String[]{suvCurvatureTexture});
+		getParent().mapUniforms(gl,new String[]{suvCurvatureMax});
+		getParent().mapUniforms(gl,new String[]{suvCurvatureMin});
+		super.init(gl);
+	}
+
+	@Override
+	public void disposeGL(GL4 gl2) {
+		for(Texture t: curvatureTexture.values()){
+			t.delete(gl2);
+		}
+		curvatureTexture.clear();
+		evaluatedCurvatures.clear();
+		super.disposeGL(gl2);
+	}
+	
+	@Override
+	public void updateData(GL4 gl) {
+		VolumeDataManager dataManager = getParent().getDataManager();
+		float globalMin = Float.MAX_VALUE;
+		float globalMax = Float.MIN_VALUE;
+		
+		//create and update laplace textures 
+		for(Integer key:dataManager.getVolumeKeys()){
+			VolumeDataBlock data = dataManager.getVolume(key);
+			Texture t;
+			
+			//create texture objects
+			if(!curvatureTexture.containsKey(key)||needsReset){
+				t = createVolumeTexture(gl, getParent().getArrayEntryLocation(gl, suvCurvatureTexture,key)) ;
+				curvatureTexture.put(key, t);
+			}
+			t = curvatureTexture.get(key);
+			
+			//update laplacians
+			if(data.needsUpdate()||!evaluatedCurvatures.containsKey(key)||needsReset){
+				CurvatureContainer container = calculateCurvatureOfVolume(data); 
+				evaluatedCurvatures.put(key,container);
+				
+			
+				//fill textures
+				FloatBuffer buffer = Buffers.newDirectFloatBuffer(container.valueMesh3d);
+				t.update(gl, 0, buffer, container.dimension);
+				
+			}
+			CurvatureContainer tmp = evaluatedCurvatures.get(key);
+			
+			//get global min max value of laplace
+			globalMax = Math.max(globalMax, tmp.maxValue);
+			globalMin = Math.min(globalMin, tmp.minValue);
+		}
+		//update globals
+		gl.glUniform1f(getParent().getLocation(suvCurvatureMin), globalMin);
+		gl.glUniform1f(getParent().getLocation(suvCurvatureMax), globalMax);
+		
+		needsReset = false;
+		super.updateData(gl);
+	}
+	
+	
+	
+	@Override
+	public void setParent(MultiVolumeRenderer parent) {
+		needsReset = true;
+		super.setParent(parent);
+	}
+	
 	@Override
 	public String[] declaration() {
 		List<String> code = new ArrayList<String>();
 		String[] dec= new String[]{
 				"#line "+Thread.currentThread().getStackTrace()[1].getLineNumber()+ " 6666",
-				"const float curvoffset = 0.2;",
-				"mat3x3 getIdentity(){",
-				"	mat3x3 identity = mat3x3(0.0);",
-				"	for(int i =0; i < 3; i++){",
-				"		identity[i][i] = 1.0;",
-				"	}",
-				"	return identity;",
-				"}",
-				"",
-				"float frobeniusNorm(mat3x3 m){",
-				"	float f = 0.0;",
-				"	for(int i = 0; i< 3; i++){",
-				"		for(int j = 0; j < 3; j++){",
-				"			f+= m[i][j]*m[i][j];",
-				"		}",
-				"	}",
-				"	f= sqrt(f);",
-				"	return f;",
-				"}",
-				"",
-				"float trace(mat3x3 m){",
-				"	float t = 0.0;",
-				"	for(int i =0; i< 3; i++){",
-				"		t+= m[i][i];",
-				"	}",
-				"	return t;",
-				"}",
-				"",
-				"vec3 gradCentral(mat3x3 v[3]){",
-				"	return vec3(v[2][1][1]-v[0][1][1],v[1][2][1]-v[1][0][1],v[1][1][2]-v[1][1][0])/(2.0*curvoffset);",
-				"}",
-				"",
-				"mat3x3[3] getValuesForGradient(vec3 pos, int volume){",
-				"	mat3x3 values[3];",
-				"	vec3 offset = vec3(curvoffset);",
-				"	for(int z = 0; z < 3; z++){",
-				"		for(int y = 0; y < 3; y++){",
-				"			for(int x = 0; x < 3; x++){",
-				"				vec3 offsetFactor= vec3(x-1,y-1,z-1);",
-				"				vec3 texC = getCorrectedTexturePositions(pos + offsetFactor*offset, volume);",
-				"				values[x][y][z] = texture("+suvVolumeTexture+"[volume],texC).r;",	
-				"			}",
-				"		}",
-				"	}",
-				"	return values;",
-				"}",
-				"",
-				"float dfdn(mat3x3 v[3],int n, int off, int offn){",
-				"	if(n==0){",
-				"		if(offn == 1){",
-				"			return (v[2][off][1]-v[0][off][1])/(2.0*curvoffset);",
-				"		}else{",
-				"			return (v[2][1][off]-v[0][1][off])/(2.0*curvoffset);",
-				"		}",
-				"	}",
-				"	if(n==1){",
-				"		if(offn == 0){",
-				"			return (v[off][2][1]-v[off][0][1])/(2.0*curvoffset);",
-				"		}else{",
-				"			return (v[1][2][off]-v[1][0][off])/(2.0*curvoffset);",
-				"		}",
-				"	}",
-				"	if(n==2){",
-				"		if(offn == 0){",
-				"			return (v[off][1][2]-v[off][1][0])/(2.0*curvoffset);",
-				"		}else{",
-				"			return (v[1][off][2]-v[1][off][0])/(2.0*curvoffset);",
-				"		}",
-				"	}",
-				"	return 0.0;",
-				"}",
-				"",
-				"mat3x3 getHessianMatrix(mat3x3 v[3]){",
-				"	mat3x3 hessian;",
-				"	vec3 offset = vec3(curvoffset);",
-				"	for(int i = 0; i< 3; i++){",
-				"		for(int j = 0; j < 3; j++){",
-				"			hessian[i][j]=(dfdn(v,i,2,j)-dfdn(v,i,0,j))/(2.0*curvoffset);",
-				"		}",
-				"	}",
-				"	return hessian;",
-				"}",
-				"",
-				"float[2] getCurvatureFactors(vec3 position,int volume){",
-				"	//neighborhood for gradient",
-				"	mat3x3 values[3]=getValuesForGradient("+sgvRayPositions+",volume);",
-				"	mat3x3 h = getHessianMatrix(values);",
-				"	vec3 gradient = gradCentral(values);",	
-				"	mat3x3 n = mat3x3(-1.0* gradient/length(gradient),vec3(0.0),vec3(0.0));",
-				"	mat3x3 p = getIdentity() - n * transpose(n);",
-				"	mat3x3 g = -p*h*p/length(gradient);",
-				"	float t = trace(g);",
-				"	float f = frobeniusNorm(g);",
-				"	float s = sqrt(2.0* f*f - t*t);",
-				"	float k[2];",
-				"	k[0] = (t+s)/2.0;",
-				"	k[1] = (t-s)/2.0;",
-				"	return k;",
-				"}",
+				"uniform sampler3D "+suvCurvatureTexture+"["+scvMaxNumberOfVolumes+"];",
+				"uniform float "+suvCurvatureMax+";",
+				"uniform float "+suvCurvatureMin+";",
+				"float curveNormalizeFactor = 1.0/("+suvCurvatureMax+"-"+suvCurvatureMin+");",
 				"",
 				"float "+getFunctionName()+"(float densities["+scvMaxNumberOfVolumes+"]) {",
 				"	float difference = 0.0;",		
 				"	for(int n = 0; n < "+scvMaxNumberOfVolumes+"; n++){",
+			    "		vec3 texCN = getCorrectedTexturePositions("+sgvRayPositions+", n);",
 				"		for(int m = 0; m < "+scvMaxNumberOfVolumes+";m++){",
 				"			if(densities[n]<0 || densities[m]<0){",
 				"				continue;",
 				"			}",	
-				"			float km[2] = getCurvatureFactors("+sgvRayPositions+",m);",
-				"			float kn[2] = getCurvatureFactors("+sgvRayPositions+",n);",
+				"			vec3 texCM = getCorrectedTexturePositions("+sgvRayPositions+", m);",
+				"			float cn = texture("+suvCurvatureTexture+"[n],texCN).r;",
+				"			float cm = texture("+suvCurvatureTexture+"[m],texCM).r;",
 				"",
+				"			"+sgvVolumeNormalizeFactor+" = curveNormalizeFactor;",
 				"			//manhatten distance",
-				"			float currentDifference =abs(3.0*length(vec2(km[0],km[1])))-abs(3.0*length(vec2(kn[0],kn[1]))); // max(km[0]-kn[0]) + abs(km[1]-kn[1]) ;",
+				"			float currentDifference = cn-cm;",
 				"			difference = max(difference,currentDifference);",	
 				"		}",
 				"	}",

src/main/java/bdv/jogl/VolumeRenderer/utils/CurvatureContainer.java

+package bdv.jogl.VolumeRenderer.utils;
+
+import java.util.TreeMap;
+
+public class CurvatureContainer {
+	
+	public float[] valueMesh3d;
+	public final TreeMap<Float,Integer> distribution = new TreeMap<Float, Integer>(); 
+	public float minValue = Float.MAX_VALUE;
+	public float maxValue = Float.MIN_VALUE;
+	public int dimension[]; 
+}

src/main/java/bdv/jogl/VolumeRenderer/utils/GradientContainer.java

+package bdv.jogl.VolumeRenderer.utils;
+
+public class GradientContainer {
+	public float[][] valueMesh3d;
+	public int dimension[]; 
+}

src/main/java/bdv/jogl/VolumeRenderer/utils/HessianContainer.java

+package bdv.jogl.VolumeRenderer.utils;
+
+public class HessianContainer {
+	public float[][] valueMesh3d;
+	public int dimension[]; 
+}

src/main/java/bdv/jogl/VolumeRenderer/utils/MatrixUtils.java

@@ -89,6 +89,98 @@ public class MatrixUtils {
 		return transformation;
 	}
 	
+	/**
+	 * 
+	 * @return a Matrix filled with zeros
+	 */
+	public static Matrix4 getNewNullMatrix(){
+		Matrix4 nullMatrix = getNewIdentityMatrix();
+		for(int i =0; i < 4; i++){
+			nullMatrix.getMatrix()[i*4+i]=0;
+		}
+		return nullMatrix;
+	}
+	
+	/**
+	 * adds two matrices and returns a new matrix object containing the result
+	 * @param a
+	 * @param b
+	 * @return
+	 */
+	public static Matrix4 addMatrix(Matrix4 a, Matrix4 b){
+		Matrix4 ret = new Matrix4();
+		for(int i =0; i < 16; i++){
+			ret.getMatrix()[i] = a.getMatrix()[i] + b.getMatrix()[i];
+		}
+		return ret;
+	}
+	
+	/**
+	 * subs two matrices and returns a new matrix object containing the result
+	 * @param a
+	 * @param b
+	 * @return
+	 */
+	public static Matrix4 subMatrix(Matrix4 a, Matrix4 b){
+		Matrix4 binv = new Matrix4();
+		for(int i = 0; i < 16; i++ ){
+			binv.getMatrix()[i] = -b.getMatrix()[i];
+		}
+		
+		return addMatrix(a, binv);
+	}
+	
+	/**
+	 * Creates a matrix instance from given data
+	 * @param data 
+	 * @return
+	 */
+	public static Matrix4 createMatrix4X4(float data[]){
+		Matrix4 m = new Matrix4();
+		for(int i =0; i < 16; i++){
+			m.getMatrix()[i] = data[i];
+		}
+		return m;
+	}
+	
+	/**
+	 * scales the given matrix by n and returns a new matrix of its content
+	 * @param n
+	 * @return
+	 */
+	public static Matrix4 scale(Matrix4 m,  float n){
+		Matrix4 k = copyMatrix(m);
+		for(int i=0; i < 16; i++){
+			k.getMatrix()[i]*=n;
+		}
+		return k;
+	}
+	
+	/**
+	 * Does matmul without altering a and b
+	 * @param a
+	 * @param b
+	 * @return
+	 */
+	public static Matrix4 matMul(Matrix4 a, Matrix4 b){
+		Matrix4 c = copyMatrix(a);
+		c.multMatrix(b);
+		return c;
+	}
+	
+	/**
+	 * Calculates the trace of a given Matrix 
+	 * @param m
+	 * @return
+	 */
+	public static float trace(Matrix4 m){
+		float t = 0;
+		for(int i =0; i < 4; i++){
+			t += m.getMatrix()[i*4+i];
+		}
+		return t;
+	}
+	
 	/**
 	 * Calculates a closely fitting bounding box of transformations 
 	 * @param transformations The list of transformations of 0-1 space coordinates

src/main/java/bdv/jogl/VolumeRenderer/utils/VolumeDataUtils.java

@@ -71,13 +71,13 @@ public class VolumeDataUtils {
 		tmp = Views.flatIterable(Views.interval(dataField, minMax[0], minMax[1]));
 
 		String name = source.getSpimSource().getName();
-		
+
 		int maxOcc =0;
 		// copy values 
 		int i = 0;
 		float minValue = Float.MAX_VALUE;
 		float maxValue = Float.MIN_VALUE; 
-		
+
 		@SuppressWarnings("unchecked")
 		Iterator<UnsignedShortType> values = (Iterator<UnsignedShortType>) tmp.iterator();
 		for(long z = minMax[0][2]; z <=  minMax[1][2]; z++){
@@ -100,7 +100,7 @@ public class VolumeDataUtils {
 				}
 			}
 		}
-		
+
 		tmp.min(data.memOffset);
 		for(int d = 0;d < 3 ;d++){
 			data.memSize[d] = minMax[1][d]+1-minMax[0][d];
@@ -223,7 +223,7 @@ public class VolumeDataUtils {
 	public static Matrix4 fromVolumeToGlobalSpace(final VolumeDataBlock block){
 		return copyMatrix( block.getLocalTransformation());
 	}
-	
+
 	public static Matrix4 calcScaledVolumeTransformation(final VolumeDataBlock block){
 		Matrix4 trans = getNewIdentityMatrix();
 
@@ -234,7 +234,7 @@ public class VolumeDataUtils {
 	}
 
 	public static Matrix4 fromCubeToVolumeSpace(final VolumeDataBlock block){
-		
+
 		Matrix4 tmp = copyMatrix(block.getLocalTransformation());
 		//tmp.scale(block.dimensions[0], block.dimensions[1], block.dimensions[2]);
 		tmp.invert();
@@ -245,7 +245,7 @@ public class VolumeDataUtils {
 		return tmp;
 		//return trans;
 	}
-	
+
 
 	/**
 	 * convolves the data in x and y dimension since these dimension have the highest resolution
@@ -262,10 +262,10 @@ public class VolumeDataUtils {
 											  {-1,-1,-1}};*/
 		int xyz[] = new int[]{(int)data.memSize[0],(int)data.memSize[1],(int)data.memSize[2]};
 		float convolved[] = new float[(int)(xyz[2]*xyz[1]*xyz[0])];
-		
+
 		LaplaceContainer ret = new LaplaceContainer();
 		ret.dimension = xyz.clone();
-		
+
 		//folding
 		for(int z = 0; z < xyz[2]; z++){
 			int zOffset = z*xyz[0]*xyz[1];
@@ -273,7 +273,7 @@ public class VolumeDataUtils {
 				float kernelValue = 0f;
 				int yOffset = y*xyz[0];
 				int vertOffset = xyz[0];
-				
+
 
 				for(int x=0; x < xyz[0]; x++){
 					//simple kernel evaluation
@@ -289,7 +289,7 @@ public class VolumeDataUtils {
 								dataValue = data.data[zOffset + yOffset  + vertOffset*y1 +x +x1];
 							}
 							float koeffizient;
-							
+
 							if(y1 == 0 && x1 == 0){
 								koeffizient = -8;
 							}else{
@@ -307,7 +307,7 @@ public class VolumeDataUtils {
 		ret.valueMesh3d = convolved;
 		return ret;
 	}
-	
+
 	/**
 	 * Creates a default volume texture (interpolation, clamping , etc,)
 	 * @param gl The gl context.
@@ -325,8 +325,8 @@ public class VolumeDataUtils {
 		volumeTexture.setTexParameteri(gl, GL4.GL_TEXTURE_WRAP_R, GL4.GL_CLAMP_TO_BORDER);
 		return volumeTexture;
 	}
-	
-	
+
+
 	/**
 	 * Calculates eye and center positions of a camera defined by a hull volume and a view direction
 	 * @param hull the hull volume
@@ -339,13 +339,199 @@ public class VolumeDataUtils {
 		float normal[] = new float[3];
 		float step [] = new float[3];
 		float[] eye = {0,0,0};
-		
-		
+
+
 		VectorUtil.normalizeVec3(normal, viewDirection.clone());
 		VectorUtil.scaleVec3(step, normal.clone(), distanceFromCenter);
 		VectorUtil.subVec3(eye, center.clone(), step);
-		
+
 		return new float[][]{eye.clone(),center.clone()}; 
 	}
+
+	private static float minStep = 0.1f;
+
+	public static CurvatureContainer calculateCurvatureOfVolume(final VolumeDataBlock data){
+		CurvatureContainer cc = new CurvatureContainer();
+		GradientContainer cg= calculateGradientOfVolume(data);
+		HessianContainer ch = calculateHessianOfGradients(cg, data);
+
+		int xyz[] = new int[]{(int)data.memSize[0],(int)data.memSize[1],(int)data.memSize[2]};
+		float convolved[] = new float[(int)(xyz[2]*xyz[1]*xyz[0])];
+
+		cc.dimension = xyz.clone();
+		//folding
+		for(int z = 0; z < xyz[2]; z++){
+			int zOffset = z*xyz[0]*xyz[1];
+			for(int y = 0; y < xyz[1]; y++){
+				int yOffset = y*xyz[0];
+
+
+
+				for(int x=0; x < xyz[0]; x++){
+					float totalCurvature = 0;
+					float g[] = cg.valueMesh3d[zOffset+yOffset+x].clone();
+					float h[] = ch.valueMesh3d[zOffset+yOffset+x].clone();
+					float l = VectorUtil.normVec3(g);
+
+					if(l>0.00001){
+						//get the normal to column matrix
+						Matrix4 n = getNewNullMatrix();
+						for(int d = 0; d < 3; d++){
+							n.getMatrix()[d*4]=-1f*g[d]/l;
+						}
+						//create nt
+						Matrix4 nt = copyMatrix(n);
+						nt.transpose();
+
+						Matrix4 H = createMatrix4X4(new float[]{
+								h[0],h[1],h[2],0,
+								h[3],h[4],h[5],0,
+								h[6],h[7],h[8],0,
+								0,0,0,0
+						});
+
+						//mat3x3 p = getIdentity() - n * transpose(n);
+						Matrix4 P = subMatrix(getNewIdentityMatrix(), matMul(n, nt));
+
+						//mat3x3 g = -p*h*p/length(gradient);
+						Matrix4 G =scale(matMul(scale(P, -1), matMul(H,P)),1f/l);
+						Matrix4 Gt = copyMatrix(G);
+						Gt.transpose();
+
+						totalCurvature = (float)Math.sqrt( trace(matMul( G,Gt)));
+					}else{
+						totalCurvature = 0;
+					}
+					//distribution
+					int count = 0;
+					if(cc.distribution.containsKey(totalCurvature)){
+						cc.distribution.get(totalCurvature);
+					}
+					count ++;
+					cc.distribution.put(totalCurvature, count);
+
+					cc.maxValue = Math.max(cc.maxValue, totalCurvature);
+					cc.minValue = Math.min(cc.minValue, totalCurvature);
+					convolved[zOffset+yOffset+x] = totalCurvature;
+				}
+			}
+		}
+		cc.valueMesh3d = convolved;
+		/*"	mat3x3 n = mat3x3(-1.0* gradient/length(gradient),vec3(0.0),vec3(0.0));",
+		"	mat3x3 p = getIdentity() - n * transpose(n);",
+		"	mat3x3 g = -p*h*p/length(gradient);",
+		"	float t = trace(g);",
+		"	float f = frobeniusNorm(g);",
+		"	float s = sqrt(2.0* f*f - t*t);",
+		"	float k[2];",
+		"	k[0] = (t+s)/2.0;",
+		"	k[1] = (t-s)/2.0;",*/
+		return cc;
+	}
+
+
+	private static HessianContainer calculateHessianOfGradients(
+			GradientContainer cg, VolumeDataBlock data) {
+		//get voxel distances per dim
+		float voxelDist[] = calculateVoxelDistance(data);		
+		HessianContainer hc = new HessianContainer();
+
+		int xyz[] = new int[]{(int)data.memSize[0],(int)data.memSize[1],(int)data.memSize[2]};
+		hc.valueMesh3d = new float[(int)(xyz[2]*xyz[1]*xyz[0])][9];
+		hc.dimension = xyz.clone();
+		for(int z = 0; z < xyz[2]; z++){
+			int zOffset = z*xyz[0]*xyz[1];
+			for(int y = 0; y < xyz[1]; y++){
+				int yOffset = y*xyz[0];
+				for(int x = 0; x < xyz[0]; x++ ){
+					int cPos[] = {x,y,z};
+					int index = zOffset+yOffset+x;
+					int offsets[] = {1,xyz[0],xyz[0]*xyz[1]};
+					//central diff optimizable through symetry
+					for(int i = 0; i < 3; i++ ){
+						for(int j = i; j< 3; j++){
+
+							float[] gpdj = (cPos[j] < xyz[j]-1)?cg.valueMesh3d[index + offsets[j]]:cg.valueMesh3d[index];
+							float[] gmdj = (cPos[j] > 0)?cg.valueMesh3d[index - offsets[j]]:cg.valueMesh3d[index];
+							float[] gpdi = (cPos[i] < xyz[i]-1)?cg.valueMesh3d[index + offsets[i]]:cg.valueMesh3d[index];
+							float[] gmdi = (cPos[i] > 0)?cg.valueMesh3d[index - offsets[i]]:cg.valueMesh3d[index];
+
+							float partialDerivative=
+									(gpdj[i] - gmdj[i])/4f*voxelDist[j]+
+									(gpdi[j] - gmdi[j])/4f*voxelDist[i];
+
+							hc.valueMesh3d[index][i*3+j] = partialDerivative;
+							//symetry
+							hc.valueMesh3d[index][j*3+i] = partialDerivative;
+
+						}
+					}
+				}
+			}
+		}
+
+		return hc;
+	}
+
+
+	/**
+	 * Calculates for each data point a gradient vector
+	 * @param data
+	 * @return
+	 */
+	public static GradientContainer calculateGradientOfVolume(
+			VolumeDataBlock data) {
+		//get voxel distances per dim
+		float voxelDist[] = calculateVoxelDistance(data);
+
+		GradientContainer gc = new GradientContainer();
+
+		int xyz[] = new int[]{(int)data.memSize[0],(int)data.memSize[1],(int)data.memSize[2]};
+		gc.valueMesh3d = new float[(int)(xyz[2]*xyz[1]*xyz[0])][3];
+
+		gc.dimension = xyz.clone();
+		//folding
+		for(int z = 0; z < xyz[2]; z++){
+			int zOffset = z*xyz[0]*xyz[1];
+			for(int y = 0; y < xyz[1]; y++){
+				int yOffset = y*xyz[0];
+				for(int x = 0; x < xyz[0]; x++ ){
+					int cPos[] = {x,y,z};
+					int index = zOffset+yOffset+x;
+					int offsets[] = {1,xyz[0],xyz[0]*xyz[1]};
+					//central diff
+					for(int d = 0; d < 3; d++ ){
+
+						float higherPosValue = (cPos[d] < xyz[d]-1)?data.data[index + offsets[d]]:data.data[index];
+						float lowerPosValue = (cPos[d] > 0)?data.data[index - offsets[d]]:data.data[index];
+						gc.valueMesh3d[index][d] = (higherPosValue - lowerPosValue)/(2.f*voxelDist[d]);
+					}
+
+				}
+			}
+		}
+		return gc;
+	}
+
+	/**
+	 * Interpolates linear between 2 points 
+	 * @param v1 startvalue
+	 * @param v2 endvalue
+	 * @param dist distance between start and end
+	 * @param minStep distance from start to interpolate
+	 * @return 
+	 */
+	private static float interpolate(float v1, float v2, float dist, float minStep) {
+		float m = (v2 -v1) /dist;
+		return v1+ m*minStep;
+	}
+
+
+	private static float[] calculateVoxelDistance(VolumeDataBlock data) {
+		// TODO Auto-generated method stub
+		// TODO correct dim
+		float distances[] = new float[]{1,1,1}; 
+		return distances;
+	}
 }
 

src/test/java/bdv/jogl/VolumeRenderer/tests/MatrixUtilsTest.java

 package bdv.jogl.VolumeRenderer.tests;
 
+import java.beans.Expression;
+
 import net.imglib2.realtransform.AffineTransform3D;
 
 import org.junit.Test;
@@ -16,26 +18,26 @@ import com.jogamp.opengl.math.Matrix4;
  */
 public class MatrixUtilsTest {
 	private float floatAccuracy = 0.0001f;
-	
+
 	@Test
 	public void identityMatrixTest(){
 		Matrix4 identity = new Matrix4();
 		identity.loadIdentity();
-		
+
 		Matrix4 testMatrix = MatrixUtils.getNewIdentityMatrix();
 		assertArrayEquals(identity.getMatrix(),testMatrix.getMatrix(),floatAccuracy);
 	}
-	
+
 	@Test
 	public void copyMatrixTest(){
 		Matrix4 matrixToCopy = new Matrix4();
 		matrixToCopy.loadIdentity();
 		matrixToCopy.makePerspective(23, 42, 0.1f, 1000);
-		
+
 		Matrix4 testMatrix = MatrixUtils.copyMatrix(matrixToCopy);
 		assertArrayEquals(matrixToCopy.getMatrix(),testMatrix.getMatrix(),floatAccuracy);
 	}
-	
+
 	@Test
 	public void convertMatrixTest(){
 		float[] matrixContent = {
@@ -44,27 +46,27 @@ public class MatrixUtilsTest {
 				12,45,43,12,
 				0,0,0,1
 		};
-		
+
 		float[] s = {
 				10,0,0,0,
 				0,20,0,0,
 				0,0,30,0,
 				0,0,0,1
 		};
-		
+
 		float[] t = {
 				1,0,0,10,
 				0,1,0,20,
 				0,0,1,30,
 				0,0,0,1
 		};
-		
+
 		//check matrix
 		Matrix4 matrixWanted = new Matrix4();
 		matrixWanted.loadIdentity();
 		matrixWanted.multMatrix(matrixContent);
 		matrixWanted.transpose();
-		
+
 		//test matrix 4*3
 		AffineTransform3D matrixToConvert = new AffineTransform3D();
 		matrixToConvert.set(
@@ -72,11 +74,11 @@ public class MatrixUtilsTest {
 				matrixContent[4],matrixContent[5],matrixContent[6],matrixContent[7],
 				matrixContent[8],matrixContent[9],matrixContent[10],matrixContent[11]
 				);
-		
-		
+
+
 		Matrix4 testMatrix = MatrixUtils.convertToJoglTransform(matrixToConvert);
 		assertArrayEquals(matrixWanted.getMatrix(), testMatrix.getMatrix(), floatAccuracy);
-		
+
 		float[]testVec = {1,1,1,1};
 		float[]testTrans = new float[4];
 		AffineTransform3D sa = new AffineTransform3D();
@@ -84,17 +86,74 @@ public class MatrixUtilsTest {
 				s[0],s[1],s[2],s[3],
 				s[4],s[5],s[6],s[7],
 				s[8],s[9],s[10],s[11]
-			);
+				);
 		MatrixUtils.convertToJoglTransform(sa).multVec(testVec, testTrans);
 		assertArrayEquals(new float[]{10,20,30,1}, testTrans,0.01f);
-		
+
 		AffineTransform3D ta = new AffineTransform3D();
 		ta.set(
 				t[0],t[1],t[2],t[3],
 				t[4],t[5],t[6],t[7],
 				t[8],t[9],t[10],t[11]
-			);
+				);
 		MatrixUtils.convertToJoglTransform(ta).multVec(testVec, testTrans);
 		assertArrayEquals(new float[]{11,21,31,1}, testTrans,0.01f);		
 	}
+
+	@Test
+	public void zeroMatrixTest(){
+		Matrix4 zeroMatrix = MatrixUtils.getNewNullMatrix();
+		for(int i = 0; i < 16;i++){
+			assertEquals(0, zeroMatrix.getMatrix()[0], 0.000f);
+		}
+	} 
+
+	@Test
+	public void addMatrixTest(){
+		Matrix4 a = MatrixUtils.createMatrix4X4(new float []{
+				1,1,1,1,
+				2,2,2,2,
+				3,3,3,3,
+				4,4,4,4});
+		Matrix4 b = MatrixUtils.createMatrix4X4(new float []{
+				1,2,3,4,
+				1,2,3,4,
+				1,2,3,4,
+				1,2,3,4});
+		Matrix4 expected = MatrixUtils.createMatrix4X4(new float []{
+				2,3,4,5,
+				3,4,5,6,
+				4,5,6,7,
+				5,6,7,8
+		});
+
+		Matrix4 r = MatrixUtils.addMatrix(a, b);
+		
+		assertArrayEquals(expected.getMatrix(), r.getMatrix(), 0.0001f);
+	}
+	
+	@Test
+	public void subMatrixTest(){
+		Matrix4 a = MatrixUtils.createMatrix4X4(new float []{
+				1,1,1,1,
+				2,2,2,2,
+				3,3,3,3,
+				4,4,4,4});
+		Matrix4 b = MatrixUtils.createMatrix4X4(new float []{
+				1,2,3,4,
+				1,2,3,4,
+				1,2,3,4,
+				1,2,3,4});
+		Matrix4 expected = MatrixUtils.createMatrix4X4(new float []{
+				0,-1,-2,-3,
+				1,0,-1,-2,
+				2,1,0,-1,
+				3,2,1,0
+		});
+
+		Matrix4 r = MatrixUtils.subMatrix(a, b);
+		
+		assertArrayEquals(expected.getMatrix(), r.getMatrix(), 0.0001f);
+	}
+	
 }

src/test/java/bdv/jogl/VolumeRenderer/tests/VolumeDataUtilsTest.java

+package bdv.jogl.VolumeRenderer.tests;
+
+import static org.junit.Assert.*;
+
+import org.junit.Test;
+
+import com.jogamp.opengl.math.VectorUtil;
+
+import bdv.jogl.VolumeRenderer.utils.CurvatureContainer;
+import bdv.jogl.VolumeRenderer.utils.GradientContainer;
+import bdv.jogl.VolumeRenderer.utils.MatrixUtils;
+import bdv.jogl.VolumeRenderer.utils.VolumeDataBlock;
+import bdv.jogl.VolumeRenderer.utils.VolumeDataUtils;
+
+public class VolumeDataUtilsTest {
+
+	@Test
+	public void testCalculateCurvatureOfVolume() {
+		VolumeDataBlock testBlock = new VolumeDataBlock();
+		testBlock.data = new float[]{0,0,0,
+									 0,0,0,
+									 0,0,0,
+			
+									 1,1,1,
+									 1,1,1,
+									 1,1,1,
+									 
+									 2,2,2,
+									 2,2,2,
+									 2,2,2};
+		testBlock.memSize = new long[]{3,3,3};
+		CurvatureContainer c = VolumeDataUtils.calculateCurvatureOfVolume(testBlock);
+		//test center gradient zero curvature on homogene data
+		assertEquals(0, c.valueMesh3d[13],0.01);
+	}
+
+	@Test
+	public void testCalculateGradientOfVolume() {
+		VolumeDataBlock testBlock = new VolumeDataBlock();
+		testBlock.data = new float[]{0,0,0,
+									 0,0,0,
+									 0,0,0,
+			
+									 1,1,1,
+									 1,1,1,
+									 1,1,1,
+									 
+									 2,2,2,
+									 2,2,2,
+									 2,2,2};
+		testBlock.memSize = new long[]{3,3,3};
+		float wantedDir[] = new float[]{0,0,1};
+		testBlock.setLocalTransformation(MatrixUtils.getNewIdentityMatrix());
+		
+		GradientContainer c = VolumeDataUtils.calculateGradientOfVolume(testBlock);
+		
+		//test center gradient
+		assertEquals(1f, Math.abs(VectorUtil.dotVec3(wantedDir, c.valueMesh3d[13])),0.01);
+	}
+
+}