Volume Rendering in Cesium

1. A concise explanation of the problem you’re experiencing.

I want to integrate the volume raycaster into Cesium,I’ve been doing it these days.I use the Volume Ray Casting (http://www.lebarba.com/) for creating the rays and raymarching ,but it still failed,I doubt whether the rays cooridnate is incorrect in my way.

2. A minimal code example. If you’ve found a bug, this helps us reproduce and repair it.

class CubePrimitive {

constructor(options) {

this.geometry = options.geometry;

this.uniformMap = options.uniformMap;

this.backVertexShader = `precision highp float;

attribute vec3 position;

varying vec3 v_pos;

void main(){

v_pos=position;

gl_Position = czm_modelViewProjection * vec4(position,1.0);

}`;

this.backFragmentShader = `precision highp float;

const float xyzLimit=30000000.0;

varying vec3 v_pos;

void main(){

gl_FragColor = vec4((v_pos+xyzLimit/2.0)/xyzLimit,1.0); //move it to [0,1]

}`;

this.frontVertexShader = `precision highp float;

const float xyzLimit=30000000.0;

attribute vec3 position;

varying vec3 v_pos;

varying vec4 v_prjPos;

void main(){

v_pos=position;

v_prjPos=czm_modelViewProjection * vec4(position,1.0);

gl_Position = v_prjPos;

}`;

this.frontFragmentShader = `precision highp float;

const float xyzLimit=30000000.0;

uniform sampler2D u_volumeTexture;

uniform sampler2D u_colorTexture;

uniform sampler2D u_backTexture;

uniform sampler2D u_depthTexture;

uniform vec3 u_minPos;

uniform vec3 u_maxPos;

varying vec3 v_pos;

varying vec4 v_prjPos;

const float zdef=63.0;

const float zslice=8.0;

float sampleAs3DTexture(vec3 texCoord){

float colorSlice1, colorSlice2;

vec2 texCoordSlice1, texCoordSlice2;

float zSliceNumber1 = floor(texCoord.z * zdef);

float zSliceNumber2 = min( zSliceNumber1 + 1.0, zdef);

float zDifference = mod(texCoord.z * zdef, 1.0);

texCoord.xy /= zslice;

texCoordSlice1 = texCoordSlice2 = texCoord.xy;

texCoordSlice1.x += mod(zSliceNumber1, zslice ) / zslice;

texCoordSlice1.y += floor((zdef - zSliceNumber1) / zslice) / zslice;

texCoordSlice1.y=1.0-texCoordSlice1.y;

texCoordSlice2.x += mod(zSliceNumber2, zslice ) / zslice;

texCoordSlice2.y += floor((zdef - zSliceNumber2) / zslice) / zslice;

texCoordSlice2.y=1.0-texCoordSlice2.y;

colorSlice1 = texture2D(u_volumeTexture, texCoordSlice1).r;

colorSlice2 = texture2D(u_volumeTexture, texCoordSlice2).r;

return mix(colorSlice1, colorSlice2, zDifference);

}

const vec3 oneOverRadii=vec3(1.567855942887398e-7,1.567855942887398e-7,1.573130351105623e-7);

const vec3 oneOverRadiiSquared=vec3(2.458172257647332e-14,2.458172257647332e-14,2.4747391015697002e-14);

const float centerToleranceSquared=0.1;

// from cesium source code

vec3 cartesianToCartographic(vec3 cartesian){

//scaleToGeodeticSurface

float x2 = cartesian.x * cartesian.x * oneOverRadii.x * oneOverRadii.x;

float y2 = cartesian.y * cartesian.y * oneOverRadii.y * oneOverRadii.y;

float z2 = cartesian.z * cartesian.z * oneOverRadii.z * oneOverRadii.z;

float squaredNorm = x2 + y2 + z2;

float ratio = sqrt(1.0 / squaredNorm);

vec3 intersection = cartesian*ratio;

// If the position is near the center, the iteration will not converge.

if (squaredNorm < centerToleranceSquared) {

}

vec3 gradient=intersectiononeOverRadiiSquared2.0;

float lambda = (1.0 - ratio) * length(cartesian) / (0.5 * length(gradient));

float correction = 0.0;

float func;

float denominator;

float xMultiplier;

float yMultiplier;

float zMultiplier;

float xMultiplier2;

float yMultiplier2;

float zMultiplier2;

float xMultiplier3;

float yMultiplier3;

float zMultiplier3;

const int loopCount = 100;

for (int i=0; i < loopCount; ++i) {

lambda -= correction;

xMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquared.x);

yMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquared.y);

zMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquared.z);

xMultiplier2 = xMultiplier * xMultiplier;

yMultiplier2 = yMultiplier * yMultiplier;

zMultiplier2 = zMultiplier * zMultiplier;

xMultiplier3 = xMultiplier2 * xMultiplier;

yMultiplier3 = yMultiplier2 * yMultiplier;

zMultiplier3 = zMultiplier2 * zMultiplier;

func = x2 * xMultiplier2 + y2 * yMultiplier2 + z2 * zMultiplier2 - 1.0;

denominator = x2 * xMultiplier3 * oneOverRadiiSquared.x +

y2 * yMultiplier3 * oneOverRadiiSquared.y +

z2 * zMultiplier3 * oneOverRadiiSquared.z;

float derivative = -2.0 * denominator;

correction = func / derivative;

if (abs(func) > 1e-12) break;

}

vec3 p=cartesian*vec3(xMultiplier,yMultiplier,zMultiplier);

//geodeticSurfaceNormal

vec3 n=p*oneOverRadiiSquared;

n=normalize(n);

vec3 h = cartesian-p;

float longitude = degrees(atan(n.y, n.x));

float latitude = degrees(asin(n.z));

float height = sign(dot(h, cartesian)) * length(h);

return vec3(longitude-u_minPos.x,latitude-u_minPos.y,height-u_minPos.z)/(u_maxPos-u_minPos);

//return (cartesian+xyzLimit/2.0)/xyzLimit;

}

void main(){

vec2 texc = vec2(((v_prjPos.x / v_prjPos.w) + 1.0 ) / 2.0,

((v_prjPos.y / v_prjPos.w) + 1.0 ) / 2.0 );

vec3 backPos = texture2D(u_backTexture, texc).xyz;

backPos=backPos*xyzLimit-xyzLimit/2.0;

vec3 frontPos = v_pos;

vec3 dir = backPos - frontPos;

// raycast sampling properties

const int numSamples = 50;

float stepSize = 1.0/float(numSamples);

// Perform the ray marching:

vec3 pos = frontPos;

vec3 step = normalize(dir) * stepSize;

float travel = length(dir);

vec4 color;

for (int i=0; i < numSamples; ++i) {

if (travel <= 0.0) break;

vec3 volumePos=cartesianToCartographic(pos);

float val = sampleAs3DTexture(volumePos);

vec4 val_color = vec4(texture2D(u_colorTexture, vec2(val, 0.0)).rgb, val);

color.rgb += (1.0 - color.a) * val_color.a * val_color.rgb;

color.a += (1.0 - color.a) * val_color.a;

if (color.a >= 0.95) break;

pos += step;

travel -= stepSize;

}

gl_FragColor = texture2D(u_backTexture, texc);

}`;

this.framebuffer = options.framebuffer;

this.autoClear = Cesium.defaultValue(options.autoClear, false);

this.show = true;

this.backCommand = undefined;

this.frontCommand = undefined;

this.clearCommand = undefined;

if (this.autoClear) {

this.clearCommand = new Cesium.ClearCommand({

color: new Cesium.Color(0.0, 0.0, 0.0, 0.0),

depth: 1.0,

framebuffer: this.framebuffer,

pass: Cesium.Pass.OPAQUE

});

}

}

createCommand(context, isback) {

var vertexArray = Cesium.VertexArray.fromGeometry({

context: context,

geometry: this.geometry,

bufferUsage: Cesium.BufferUsage.STATIC_DRAW,

});

var shaderProgram;

var renderState;

if (isback) {

shaderProgram = Cesium.ShaderProgram.fromCache({

context: context,

vertexShaderSource: this.backVertexShader,

fragmentShaderSource: this.backFragmentShader

});

renderState = Cesium.RenderState.fromCache({

depthTest: {

enabled: true

},

cull: {

enabled: true,

face: Cesium.CullFace.FRONT

},

depthMask: false,

blending: Cesium.BlendingState.ALPHA_BLEND

});

return new Cesium.DrawCommand({

owner: this,

vertexArray: vertexArray,

modelMatrix: Cesium.Matrix4.IDENTITY,

shaderProgram: shaderProgram,

framebuffer: this.framebuffer,

renderState: renderState,

pass: Cesium.Pass.OPAQUE

});

} else {

shaderProgram = Cesium.ShaderProgram.fromCache({

context: context,

vertexShaderSource: this.frontVertexShader,

fragmentShaderSource: this.frontFragmentShader

});

renderState = Cesium.RenderState.fromCache({

depthTest: {

enabled: true

},

cull: {

enabled: true,

face: Cesium.CullFace.BACK

},

depthMask: false,

blending: Cesium.BlendingState.ALPHA_BLEND

});

return new Cesium.DrawCommand({

owner: this,

vertexArray: vertexArray,

uniformMap: this.uniformMap,

modelMatrix: Cesium.Matrix4.IDENTITY,

shaderProgram: shaderProgram,

renderState: renderState,

pass: Cesium.Pass.OPAQUE

});

}

}

update(frameState) {

if (!this.show) {

return;

}

if (!Cesium.defined(this.backCommand)) {

this.backCommand = this.createCommand(frameState.context, true);

}

if (!Cesium.defined(this.frontCommand)) {

this.frontCommand = this.createCommand(frameState.context, false);

}

if (Cesium.defined(this.clearCommand)) {

frameState.commandList.push(this.clearCommand);

}

frameState.commandList.push(this.backCommand); //First render pass

frameState.commandList.push(this.frontCommand); //Second render pass

}

isDestroyed() {

return false;

}

destroy() {

if (Cesium.defined(this.backCommand)) {

this.backCommand.shaderProgram = this.backCommand.shaderProgram && this.backCommand.shaderProgram.destroy();

}

if (Cesium.defined(this.frontCommand)) {

this.frontCommand.shaderProgram = this.frontCommand.shaderProgram && this.frontCommand.shaderProgram.destroy();

}

return Cesium.destroyObject(this);

}

}

function render(data) {

// data neghip_64x64x64_uint8.raw

var minLon = 70;

var maxLon = 138.39092;

var minLat = 3;

var maxLat = 54.39436;

var minHeight = 0;

var maxHeight = 600000;

var position_min = new Cesium.Cartesian3.fromDegrees(minLon, minLat, minHeight);

var position_max = new Cesium.Cartesian3.fromDegrees(maxLon, maxLat, maxHeight);

var ellipsoid = viewer.scene.globe.ellipsoid;

var cartesian3 = new Cesium.Cartesian3(position_max.x, position_max.y, position_max.z);

var cartographic = ellipsoid.cartesianToCartographic(cartesian3);

var lat = Cesium.Math.toDegrees(cartographic.latitude);

var lng = Cesium.Math.toDegrees(cartographic.longitude);

var alt = cartographic.height;

var topBound = ;

var botBound = ;

var stepCount = 100;

var lonStep = (maxLon - minLon) / stepCount;

for (var i = 0; i <= stepCount; i++) {

var lon = minLon + i * lonStep;

topBound.push(lon, maxLat);

lon = maxLon - i * lonStep;

botBound.push(lon, minLat)

}

var degreesArray = [minLon, minLat,

minLon, maxLat

];

degreesArray = degreesArray.concat(topBound);

degreesArray.push(maxLon, minLat);

degreesArray = degreesArray.concat(botBound);

var polygon = new Cesium.PolygonGeometry({

polygonHierarchy: new Cesium.PolygonHierarchy(

Cesium.Cartesian3.fromDegreesArray(degreesArray)

),

height: maxHeight,

extrudedHeight: 0

});

var geometry = Cesium.PolygonGeometry.createGeometry(polygon);

var geoInstance = new Cesium.GeometryInstance({

geometry: geometry

});

var context = viewer.scene.context;

var colorRender = [

107, 226, 16, 255,

228, 214, 8, 255,

231, 210, 8, 255,

234, 202, 8, 255,

242, 128, 7, 255,

243, 113, 7, 255,

243, 102, 7, 255,

244, 83, 7, 255,

244, 70, 7, 255,

244, 60, 7, 255,

244, 43, 7, 255,

244, 29, 7, 255,

146, 5, 196, 255,

97, 5, 138, 255,

97, 5, 138, 255,

91, 5, 120, 255,

85, 5, 102, 255,

79, 5, 83, 255,

73, 5, 65, 255

];

var colorTexture = new Cesium.Texture({

context:context,

width: colorRender.length / 4,

height: 1,

pixelFormat: Cesium.PixelFormat.RGBA,

pixelDatatype: Cesium.PixelDatatype.UNSIGNED_BYTE,

source: {

arrayBufferView: new Uint8Array(colorRender)

},

sampler: new Cesium.Sampler({

wrapS: Cesium.TextureWrap.CLAMP_TO_EDGE,

wrapT: Cesium.TextureWrap.CLAMP_TO_EDGE,

minificationFilter: Cesium.TextureMinificationFilter.LINEAR,

magnificationFilter: Cesium.TextureMagnificationFilter.LINEAR

})

});

colorTexture.type = ‘sampler2D’;

var volumnTextureData = makeTexture(data);

var volumeTexture = new Cesium.Texture({

context: context,

width: volumnTextureData.width,

height: volumnTextureData.height,

pixelFormat: Cesium.PixelFormat.RGBA,

pixelDatatype: Cesium.PixelDatatype.UNSIGNED_BYTE,

source: {

arrayBufferView: volumnTextureData.data

},

sampler: new Cesium.Sampler({

wrapS: Cesium.TextureWrap.CLAMP_TO_EDGE,

wrapT: Cesium.TextureWrap.CLAMP_TO_EDGE,

minificationFilter: Cesium.TextureMinificationFilter.LINEAR,

magnificationFilter: Cesium.TextureMagnificationFilter.LINEAR

})

});

volumeTexture.type = ‘sampler2D’;

var f_backTexture = new Cesium.Texture({

context: context,

width: context.drawingBufferWidth,

height: context.drawingBufferHeight,

pixelFormat: Cesium.PixelFormat.RGBA,

pixelDatatype: Cesium.PixelDatatype.UNSIGNED_BYTE

});

f_backTexture.type = ‘sampler2D’;

var f_depthTexture = new Cesium.Texture({

context: context,

width: context.drawingBufferWidth,

height: context.drawingBufferHeight,

pixelFormat: Cesium.PixelFormat.DEPTH_COMPONENT,

pixelDatatype: Cesium.PixelDatatype.UNSIGNED_INT,

sampler: new Cesium.Sampler({

wrapS: Cesium.TextureWrap.CLAMP_TO_EDGE,

wrapT: Cesium.TextureWrap.CLAMP_TO_EDGE,

minificationFilter: Cesium.TextureMinificationFilter.LINEAR,

magnificationFilter: Cesium.TextureMagnificationFilter.LINEAR

})

});

f_depthTexture.type = ‘sampler2D’;

var framebuffer = new Cesium.Framebuffer({

context: context,

colorTextures: [f_backTexture],

depthTexture: f_depthTexture

});

var primitive = new CubePrimitive({

geometry: geometry,

primitiveType: Cesium.PrimitiveType.TRIANGLES,

uniformMap: {

u_backTexture:function(){

return f_backTexture;

},

u_depthTexture: function(){

return f_depthTexture;

},

u_volumeTexture:function(){

return volumeTexture;

},

u_colorTexture:function(){

return colorTexture;

},

u_minPos:function(){

return new Cesium.Cartesian3(minLon,minLat,minHeight);

},

u_maxPos:function(){

return new Cesium.Cartesian3(maxLon,maxLat,maxHeight);

}

},

framebuffer: framebuffer,

autoClear: true

});

viewer.scene.primitives.add(primitive);

}

3. Context. Why do you need to do this? We might know a better way to accomplish your goal.

4. The Cesium version you’re using, your operating system and browser.

Cesium1.56

Chrome74.0.3729.157

The correct result image like this(without cesium)

在 2019年5月20日星期一 UTC+8下午5:18:27,pollux写道:

What does the image look like in Cesium? Can you post a screenshot of that? Or is it not rendering at all? Are you getting any errors?

This looks like a pretty cool application - what kind of project are you working on?

Follow these steps:

1、in first render pass,I store back faces position in a framebuffer texture which called u_backTexture,
2、in second render pass, only draw the front faces,where each point of the front face will be the ray starting point.and then sample the u_backTexture to get the world space position as ending position.

3、calc the ray from the starting position to ending position vec3 dir = backPos - frontPos;

4、rendering the volume

I doubt whether the result ray is incorrect in step1 and step2.and cann’t get a right result in .cesium.

在 2019年5月21日星期二 UTC+8上午6:37:45,Omar Shehata写道:

a clearer picture, but still incorrect.It’s been bothering me for days.

在 2019年5月20日星期一 UTC+8下午5:18:27,pollux写道:

Hello,friend. a year has passed. Have you solved your problem now? Since I also need to do volume rendering in Cesium recently, I would like to get your help very much.
I’ve read your sample code and I don’t know why it’s height: maxheight, extrudedheight: 0, not height: 0, extrudedheight: maxheight when creating polygon geometry. ( I’m not sure if you want to create a simple polygon or cuboid geometry)
By the way, are you Chinese?(because your use name"liu liang" is like Chinese pinyin),If you are really Chinese, can you give me your QQ number or wechat number and other contact information……
I look forward to your reply.