vol*: rendering oblique plane through volumes

[alisterburt]

This is a high value feature for volume visualisation identified in #1086

@arose said in #1086

For the volume image slice we currently do the interpolation in the shader. For arbitrary planes we probably want to do it on the CPU side (or just disable?).

The answer this depends on whether we think significant computation will be necessary for the slice generation that benefits from it being in the shader.

• for simple slices you can get away with doing this on the CPU
• for slices generated by projecting data normal to the slice you may benefit from the GPU

The need for projection to be part of the slicing mechanism enables

• overcoming SNR constraints in individual slices
• visualising larger volume regions on the slice whilst maintaining the defined 3D->2D mouse->plane intersection for interaction

I don't know that doing this in the shader is strictly necessary for real-time performance, I just know that my implementation in the shader in napari is significantly faster than the CPU implementation of the same behaviour in dynamo_tomoslice, a viewer based on similar principles.

My shader implementation can be found here
https://github.com/vispy/vispy/blob/9b0f58907b151a1114fce708212ced59ddb7c1cc/vispy/visuals/volume.py#L326-L363

Looking at the direct volume fragment shader it should be relatively quick to get the same logic going using this shader if we agree on that strategy

molstar/src/mol-gl/shader/direct-volume.frag.ts

Lines 8 to 364 in 10a3f92

	export const directVolume_frag = `
	precision highp float;
	precision highp int;
	#include common
	#include light_frag_params
	#if dClipObjectCount != 0
	uniform int uClipObjectType[dClipObjectCount];
	uniform bool uClipObjectInvert[dClipObjectCount];
	uniform vec3 uClipObjectPosition[dClipObjectCount];
	uniform vec4 uClipObjectRotation[dClipObjectCount];
	uniform vec3 uClipObjectScale[dClipObjectCount];
	#endif
	#include common_clip
	#include read_from_texture
	#include texture3d_from_1d_trilinear
	#include texture3d_from_2d_nearest
	#include texture3d_from_2d_linear
	uniform mat4 uProjection, uTransform, uModelView, uModel, uView;
	uniform vec3 uCameraDir;
	uniform sampler2D tDepth;
	uniform vec2 uDrawingBufferSize;
	varying vec3 vOrigPos;
	varying float vInstance;
	varying vec4 vBoundingSphere;
	varying mat4 vTransform;
	uniform mat4 uInvView;
	uniform vec3 uGridDim;
	uniform vec3 uBboxSize;
	uniform sampler2D tTransferTex;
	uniform float uTransferScale;
	uniform float uStepScale;
	uniform float uJumpLength;
	uniform int uObjectId;
	uniform int uVertexCount;
	uniform int uInstanceCount;
	uniform int uGroupCount;
	#if defined(dColorMarker)
	uniform vec3 uHighlightColor;
	uniform vec3 uSelectColor;
	uniform vec3 uDimColor;
	uniform float uHighlightStrength;
	uniform float uSelectStrength;
	uniform float uDimStrength;
	uniform int uMarkerPriority;
	uniform float uMarkerAverage;
	uniform float uMarker;
	uniform vec2 uMarkerTexDim;
	uniform sampler2D tMarker;
	#endif
	uniform float uMetalness;
	uniform float uRoughness;
	uniform bool uFog;
	uniform float uFogNear;
	uniform float uFogFar;
	uniform vec3 uFogColor;
	uniform float uAlpha;
	uniform bool uTransparentBackground;
	uniform float uXrayEdgeFalloff;
	uniform float uExposure;
	uniform int uRenderMask;
	uniform float uNear;
	uniform float uFar;
	uniform float uIsOrtho;
	uniform vec3 uCellDim;
	uniform vec3 uCameraPosition;
	uniform mat4 uCartnToUnit;
	#if __VERSION__ != 100
	// for webgl1 this is given as a 'define'
	uniform int uMaxSteps;
	#endif
	#if defined(dGridTexType_2d)
	precision highp sampler2D;
	uniform sampler2D tGridTex;
	uniform vec3 uGridTexDim;
	#elif defined(dGridTexType_3d)
	precision highp sampler3D;
	uniform sampler3D tGridTex;
	#endif
	#if defined(dColorType_uniform)
	uniform vec3 uColor;
	#elif defined(dColorType_texture)
	uniform vec2 uColorTexDim;
	uniform sampler2D tColor;
	#endif
	#ifdef dOverpaint
	#if defined(dOverpaintType_groupInstance) || defined(dOverpaintType_vertexInstance)
	uniform vec2 uOverpaintTexDim;
	uniform sampler2D tOverpaint;
	#endif
	#endif
	#ifdef dUsePalette
	uniform sampler2D tPalette;
	#endif
	#if defined(dGridTexType_2d)
	vec4 textureVal(vec3 pos) {
	return texture3dFrom2dLinear(tGridTex, pos + (vec3(0.5, 0.5, 0.0) / uGridDim), uGridDim, uGridTexDim.xy);
	}
	vec4 textureGroup(vec3 pos) {
	return texture3dFrom2dNearest(tGridTex, pos + (vec3(0.5, 0.5, 0.0) / uGridDim), uGridDim, uGridTexDim.xy);
	}
	#elif defined(dGridTexType_3d)
	vec4 textureVal(vec3 pos) {
	return texture(tGridTex, pos + (vec3(0.5) / uGridDim));
	}
	vec4 textureGroup(vec3 pos) {
	return texelFetch(tGridTex, ivec3(pos * uGridDim), 0);
	}
	#endif
	float calcDepth(const in vec3 pos) {
	vec2 clipZW = pos.z * uProjection[2].zw + uProjection[3].zw;
	return 0.5 + 0.5 * clipZW.x / clipZW.y;
	}
	float transferFunction(float value) {
	return texture2D(tTransferTex, vec2(value, 0.0)).a;
	}
	float getDepth(const in vec2 coords) {
	#ifdef depthTextureSupport
	return texture2D(tDepth, coords).r;
	#else
	return unpackRGBAToDepth(texture2D(tDepth, coords));
	#endif
	}
	const float gradOffset = 0.5;
	vec3 v3m4(vec3 p, mat4 m) {
	return (m * vec4(p, 1.0)).xyz;
	}
	float preFogAlphaBlended = 0.0;
	vec4 raymarch(vec3 startLoc, vec3 step, vec3 rayDir) {
	mat3 normalMatrix = transpose3(inverse3(mat3(uModelView * vTransform)));
	mat4 cartnToUnit = uCartnToUnit * inverse4(vTransform);
	#if defined(dClipVariant_pixel) && dClipObjectCount != 0
	mat4 modelTransform = uModel * vTransform * uTransform;
	#endif
	mat4 modelViewTransform = uModelView * vTransform * uTransform;
	vec3 scaleVol = vec3(1.0) / uGridDim;
	vec3 pos = startLoc;
	vec4 cell;
	float prevValue = -1.0;
	float value = 0.0;
	vec4 src = vec4(0.0);
	vec4 dst = vec4(0.0);
	float fragmentDepth;
	vec3 posMin = vec3(0.0);
	vec3 posMax = vec3(1.0) - vec3(1.0) / uGridDim;
	vec3 unitPos;
	vec3 nextPos;
	float nextValue;
	vec4 material;
	vec4 overpaint;
	float metalness = uMetalness;
	float roughness = uRoughness;
	vec3 gradient = vec3(1.0);
	vec3 dx = vec3(gradOffset * scaleVol.x, 0.0, 0.0);
	vec3 dy = vec3(0.0, gradOffset * scaleVol.y, 0.0);
	vec3 dz = vec3(0.0, 0.0, gradOffset * scaleVol.z);
	float maxDist = min(vBoundingSphere.w * 2.0, uFar - uNear);
	float maxDistSq = maxDist * maxDist;
	for (int i = 0; i < uMaxSteps; ++i) {
	// break when beyond bounding-sphere or far-plane
	vec3 distVec = startLoc - pos;
	if (dot(distVec, distVec) > maxDistSq) break;
	unitPos = v3m4(pos, cartnToUnit);
	// continue when outside of grid
	if (unitPos.x > posMax.x || unitPos.y > posMax.y || unitPos.z > posMax.z ||
	unitPos.x < posMin.x || unitPos.y < posMin.y || unitPos.z < posMin.z
	) {
	prevValue = value;
	pos += step;
	continue;
	}
	cell = textureVal(unitPos);
	value = cell.a; // current voxel value
	if (uJumpLength > 0.0 && value < 0.01) {
	nextPos = pos + rayDir * uJumpLength;
	nextValue = textureVal(v3m4(nextPos, cartnToUnit)).a;
	if (nextValue < 0.01) {
	prevValue = nextValue;
	pos = nextPos;
	continue;
	}
	}
	vec4 mvPosition = modelViewTransform * vec4(unitPos * uGridDim, 1.0);
	if (calcDepth(mvPosition.xyz) > getDepth(gl_FragCoord.xy / uDrawingBufferSize))
	break;
	#if defined(dClipVariant_pixel) && dClipObjectCount != 0
	vec3 vModelPosition = v3m4(unitPos * uGridDim, modelTransform);
	if (clipTest(vec4(vModelPosition, 0.0))) {
	prevValue = value;
	pos += step;
	continue;
	}
	#endif
	vec3 vViewPosition = mvPosition.xyz;
	material.a = transferFunction(value);
	#ifdef dPackedGroup
	float group = unpackRGBToInt(textureGroup(floor(unitPos * uGridDim + 0.5) / uGridDim).rgb);
	#else
	vec3 g = floor(unitPos * uGridDim + 0.5);
	// note that we swap x and z because the texture is flipped around y
	#if defined(dAxisOrder_012)
	float group = g.z + g.y * uGridDim.z + g.x * uGridDim.z * uGridDim.y; // 210
	#elif defined(dAxisOrder_021)
	float group = g.y + g.z * uGridDim.y + g.x * uGridDim.y * uGridDim.z; // 120
	#elif defined(dAxisOrder_102)
	float group = g.z + g.x * uGridDim.z + g.y * uGridDim.z * uGridDim.x; // 201
	#elif defined(dAxisOrder_120)
	float group = g.x + g.z * uGridDim.x + g.y * uGridDim.x * uGridDim.z; // 021
	#elif defined(dAxisOrder_201)
	float group = g.y + g.x * uGridDim.y + g.z * uGridDim.y * uGridDim.x; // 102
	#elif defined(dAxisOrder_210)
	float group = g.x + g.y * uGridDim.x + g.z * uGridDim.x * uGridDim.y; // 012
	#endif
	#endif
	#if defined(dColorType_direct) && defined(dUsePalette)
	material.rgb = texture2D(tPalette, vec2(value, 0.0)).rgb;
	#elif defined(dColorType_uniform)
	material.rgb = uColor;
	#elif defined(dColorType_instance)
	material.rgb = readFromTexture(tColor, vInstance, uColorTexDim).rgb;
	#elif defined(dColorType_group)
	material.rgb = readFromTexture(tColor, group, uColorTexDim).rgb;
	#elif defined(dColorType_groupInstance)
	material.rgb = readFromTexture(tColor, vInstance * float(uGroupCount) + group, uColorTexDim).rgb;
	#elif defined(dColorType_vertex)
	material.rgb = texture3dFrom1dTrilinear(tColor, unitPos, uGridDim, uColorTexDim, 0.0).rgb;
	#elif defined(dColorType_vertexInstance)
	material.rgb = texture3dFrom1dTrilinear(tColor, unitPos, uGridDim, uColorTexDim, vInstance * float(uVertexCount)).rgb;
	#endif
	#ifdef dOverpaint
	#if defined(dOverpaintType_groupInstance)
	overpaint = readFromTexture(tOverpaint, vInstance * float(uGroupCount) + group, uOverpaintTexDim);
	#elif defined(dOverpaintType_vertexInstance)
	overpaint = texture3dFrom1dTrilinear(tOverpaint, unitPos, uGridDim, uOverpaintTexDim, vInstance * float(uVertexCount));
	#endif
	material.rgb = mix(material.rgb, overpaint.rgb, overpaint.a);
	#endif
	#ifdef dIgnoreLight
	gl_FragColor.rgb = material.rgb;
	#else
	if (material.a >= 0.01) {
	#ifdef dPackedGroup
	// compute gradient by central differences
	gradient.x = textureVal(unitPos - dx).a - textureVal(unitPos + dx).a;
	gradient.y = textureVal(unitPos - dy).a - textureVal(unitPos + dy).a;
	gradient.z = textureVal(unitPos - dz).a - textureVal(unitPos + dz).a;
	#else
	gradient = cell.xyz * 2.0 - 1.0;
	#endif
	vec3 normal = -normalize(normalMatrix * normalize(gradient));
	#include apply_light_color
	} else {
	gl_FragColor.rgb = material.rgb;
	}
	#endif
	gl_FragColor.a = material.a * uAlpha * uTransferScale;
	#if defined(dColorMarker)
	float marker = uMarker;
	if (uMarker == -1.0) {
	marker = readFromTexture(tMarker, vInstance * float(uGroupCount) + group, uMarkerTexDim).a;
	marker = floor(marker * 255.0 + 0.5); // rounding required to work on some cards on win
	}
	#endif
	#include apply_marker_color
	preFogAlphaBlended = (1.0 - preFogAlphaBlended) * gl_FragColor.a + preFogAlphaBlended;
	fragmentDepth = calcDepth(mvPosition.xyz);
	#include apply_fog
	src = gl_FragColor;
	if (!uTransparentBackground) {
	// done in 'apply_fog' otherwise
	src.rgb *= src.a;
	}
	dst = (1.0 - dst.a) * src + dst; // standard blending
	// break if the color is opaque enough
	if (dst.a > 0.95)
	break;
	pos += step;
	}
	return dst;
	}
	// TODO: support float texture for higher precision values???
	// TODO: support clipping exclusion texture support
	void main() {
	if (gl_FrontFacing)
	discard;
	vec3 rayDir = mix(normalize(vOrigPos - uCameraPosition), uCameraDir, uIsOrtho);
	vec3 step = rayDir * uStepScale;
	float boundingSphereNear = distance(vBoundingSphere.xyz, uCameraPosition) - vBoundingSphere.w;
	float d = max(uNear, boundingSphereNear) - mix(0.0, distance(vOrigPos, uCameraPosition), uIsOrtho);
	vec3 start = mix(uCameraPosition, vOrigPos, uIsOrtho) + (d * rayDir);
	gl_FragColor = raymarch(start, step, rayDir);
	float fragmentDepth = calcDepth((uModelView * vec4(start, 1.0)).xyz);
	float preFogAlpha = clamp(preFogAlphaBlended, 0.0, 1.0);
	#include wboit_write
	}
	`;

@arose also mentioned that for controlling the position/orientation of this plane we need a dedicated ParameterControl - I don't yet fully understand what this is but it looks like it's implemented here

molstar/src/mol-plugin-ui/controls/parameters.tsx

Line 43 in 10a3f92

export class ParameterControls<P extends PD.Params> extends React.PureComponent<ParameterControlsProps<P>> {

It is important that we can obtain the intersection of the mouse with these planes to enable annotation tools which work on the plane

@arose said in #1086

We can obtain the grid cell of the original volume when mouse over a slice pixel which should give us the data to build a UI and storage for it.

I think it's important that this intersection be available in continuous coordinates rather than giving the center of the nearest grid cell - sometimes pixel sizes get larger than ideal with the data I work with and we are doing lots of things on this data which require subpixel precision