Fun with IDLgrShader
Anonym
In my previous job I did a lot of OpenGL programming, particularly with GLSL shaders. So I was keen to see what the capabilities were for IDL's IDLgrShader class, and how to use it with the other graphics objects. This post discusses a shader I prototyped that visualizes an isocontour line on a surface, which can be very quickly moved or animated without having to recalculate its geometry.
Let's start with the PRO code that loads one of the example data files, builds a basic scene graph for visualization with a surface, and a custom shader:
functionpow2Ceil, x
Pow2 = IShft(1, Lindgen(16))
w = Where(x and Pow2, count)
return, (count eq 1) ? Pow2[w[0]] : Pow2[w[-1]+1]
end
pro SurfaceIsoContour
m = 350
n = 450
heights = Float(Read_Binary(FilePath('surface.dat', ROOT=!Dir, $
SUBDIR=['examples', 'data']), $
DATA_TYPE=2, DATA_DIMS=[m, n]))
heights *= 0.1
minHeight = Min(heights, MAX=maxHeight)
meanHeight = (minHeight + maxHeight) * 0.5
oWin = IDLitWindow(DIMENSIONS=[m, n])
oView = IDLgrView(COLOR=[150,150,255], VIEWPLANE_RECT=[-m*0.5, -n*0.5, m, n], $
ZCLIP=[n*0.5, -n*0.5], EYE=n)
oWin.SetProperty, GRAPHICS_TREE=oView
oModel = IDLgrModel()
oView->Add, oModel
oSurf = IDLgrSurface(DATAX=Findgen(m), DATAY=Findgen(n), $
DATAZ=heights, STYLE=2, USE_TRIANGLES=1)
oModel->Add, oSurf
oHeightImage = IDLgrImage(DATA=heights, INTERNAL_DATA_TYPE=3, INTERPOLATE=0)
textureSize = [1.0*pow2Ceil(m), 1.0*pow2Ceil(n)]
oShader = Obj_New('IDLgrShader')
baseDir = File_DirName((Scope_TraceBack(/STRUCTURE))[-1].filename) + Path_Sep()
vertFile = baseDir + 'HeightMapShader.vert'
fragFile = baseDir + 'SurfaceIsocontourShader.frag'
oShader->SetProperty, VERTEX_PROGRAM_FILENAME=vertFile, $
FRAGMENT_PROGRAM_FILENAME=fragFile
oShader->SetUniformVariable, 'heightMap', oHeightImage
oShader->SetUniformVariable, 'uImageDim', [textureSize, 1.0/textureSize]
oShader->SetUniformVariable, 'uContourValue', 0.0
oShader->SetUniformVariable, 'uHeightOffsetScale', [-Min(heights), $
1.0/(Max(heights) - Min(heights))]
oSurf->SetProperty, SHADER=oShader, TEXTURE_MAP=oHeightImage
oObs = Obj_New('surfaceObserver')
oObs->SetProperty, MODEL=oModel, SHADER=oShader, $
ORIGIN=[m/2.0, n/2.0, meanHeight], $
DATA_RANGE=maxHeight-minHeight, DATA_MEAN=meanHeight
oWin->SetEventMask, /TIMER_EVENTS, /BUTTON_EVENTS, /MOTION_EVENTS
oWin->AddWindowEventObserver, oObs
oWin->SetTimerInterval, 1.0/30.0
end
The SurfaceIsoContour function loads the example data file surface.dat, which contains a chunk of the Maroon Bells terrain in a regular grid. This height data is scaled down by a factor of 0.1 so that the vertical scale is comparable to the horizontal scale. The we go about building a bare bones IDLitWindow object to hold a basic scene graph, with one model that has the one surface. The IDLgrView is configured to create an orthographic view of the surface, though some minor clipping does occur at certain azimuth rotations. The surface is constructed with some Findgen arrays for the x- and y-coordinates, to make it a regular grid. The STYLE property was set to 2 for a filled surface, though this prototype will work for wire mesh (1), RuledXZ (3), and RuledYZ (4) modes as well. It works equally well whether USE_TRIANGLES is set to 0 or 1 for quads or triangles. The height array is also loaded into an IDLgrImage object, which will be passed to the shader as a texture. The shader initialization and configuration will be discussed below after I show the source code.
The SurfaceObserver class is a very basic mouse event handler, which only allows rotation around the origin with the left mouse button. It also has a timer callback, which updates the contour value over the height extents in a sinusoidal pattern. The timer is set at 30 Hz, so it should take about 12 seconds for the contour loop to repeat.
functionSurfaceObserver::init
self.elevation = 90.0
return, 1
end
pro SurfaceObserver::cleanup
end
pro SurfaceObserver::SetProperty, MODEL=model, SHADER=shader, ORIGIN=origin, $
DATA_RANGE=dataRange, DATA_MEAN=dataMean
if (N_ELEMENTS(model) gt 0) then self.oModel = model
if (N_ELEMENTS(shader) gt 0) then self.oShader = shader
if (N_ELEMENTS(dataRange) gt 0) then self.dataRange = dataRange
if (N_ELEMENTS(dataMean) gt 0) then self.dataMean = dataMean
if (N_ELEMENTS(origin) gt 0) then begin
self.origin = origin
self->UpdateEye
endif
end
pro SurfaceObserver::OnMouseDown, oWin, x, y, ButtonMask, Modifiers, NumClick
self.mouseDown = ButtonMask
self.lastMousePos = [x, y]
end
pro SurfaceObserver::OnMouseMotion, oWin, x, y, Modifiers
if ((self.mouseDown AND 1) gt 0) then begin
delta = [x, y] - self.lastMousePos
self.lastMousePos = [x, y]
self.azimuth += delta[0]
self.elevation += delta[1]
self->UpdateEye
oWin->Draw
endif
end
pro SurfaceObserver::OnMouseUp, oWin, x, y, ButtonMask
self.mouseDown = 0
end
pro SurfaceObserver::OnTimer, oWin
self.time = (self.time + 1) mod 360
contourVal = self.dataMean + 0.5 * self.dataRange * sin(self.time * !DTOR)
self.oShader->SetUniformVariable, 'uContourVal', contourVal
oWin->Draw
end
pro SurfaceObserver::UpdateEye
self.oModel->Reset
self.oModel->Translate, -self.origin[0], -self.origin[1], -self.origin[2]
self.oModel->Rotate, [0, 0, 1], self.azimuth
self.oModel->Rotate, [1, 0, 0], 90.0 - self.elevation
end
pro SurfaceObserver__define
void = { SurfaceObserver, $
oModel: Obj_New(), $
oShader: Obj_New(), $
mouseDown: 0L, $
lastMousePos: [0L, 0L], $
origin: [0.0, 0.0, 0.0], $
azimuth: 0.0, $
elevation: 0.0, $
time: 0L, $
dataRange: 0.0, $
dataMean: 0.0 $
}
end
Now let's take a look at the shader source code, since that's where the magic happens. Keep in mind that the IDLgrShader class uses the OpenGL 2.1 API, so it has all the built-in uniforms like gl_ModelViewProjectionMatrix. The vertex shader is a very basic passthrough shader, transforming the input vertex coordinate to normalized device coordinates and copying the texture coordinates, which are vital to the fragment shader. The vertex shader is stored in the file named HeightMapShader.vert, and is in the same folder as the main PRO file.
void main(void)
{
// get the height value from the texture
gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;
gl_TexCoord[0] = gl_MultiTexCoord0;
}
The main work is done in the fragment shader, which is stored in the file named SurfaceIsocontour.frag, and is in the same folder as the main PRO file. This shader defines four uniform variables:
- heightMap- This is the height map as a texture, marshalled for us by the IDLgrImage class. Since we're using OpenGL 2.0, and the sampler2D interface, the texture must be a power of 2 in both dimensions, and the sampling coordinates are in the [0, 1] range. The texture padding is automatically handled for us by IDLgrImage.
- uImageDim- This is a 4-element array [texture width, texture height, 1/texture width, 1/texture height]. Since the texture has to be a power of 2 in both dimensions, we can't use the height array's real dimensions, we need to round the dimensions up to the next power of 2, which is what the pow2Ceil function does for us.
- uHeightOffsetScale- This is a 2-element array [-min height value, 1/height range], which is used to calculate a smooth black->white color ramp for the fragment based on their interpolated height value.
- uContourVal- This is the height value at which to draw the isocontour line.
uniform sampler2D heightMap;
uniform vec4 uImageDim;
uniform vec2 uHeightOffsetScale;
uniform float uContourVal;
void main(void)
{
// Get the fractional texel coordinate for this fragment, offset by 0.5
vec2 fractST = fract(gl_TexCoord[0].st * uImageDim.st) - vec2(0.5);
vec2 dsdt = sign(fractST) * uImageDim.zw;
// use the dsdt values, which are +1 or -1, to sample four closest neighbor texels
float texVal1 = texture2D(heightMap, gl_TexCoord[0].st).r;
float texVal2 = texture2D(heightMap, gl_TexCoord[0].st + vec2(dsdt.s, 0.0)).r;
float texVal3 = texture2D(heightMap, gl_TexCoord[0].st + vec2(0.0, dsdt.t)).r;
float texVal4 = texture2D(heightMap, gl_TexCoord[0].st + dsdt).r;
// use the fractional texel coordinate values to linearly interpolate all four edges
float ave12 = mix(texVal1, texVal2, abs(fractST.s));
float ave34 = mix(texVal3, texVal4, abs(fractST.s));
float ave13 = mix(texVal1, texVal3, abs(fractST.t));
float ave24 = mix(texVal2, texVal4, abs(fractST.t));
// bilinearly interpolate fragment's value
float myValue = mix(ave12, ave34, abs(fractST.t));
// use minVal/maxVal uniforms to perform linear bytescale of fragment value
vec3 baseColor = mix(vec3(0.0), vec3(1.0), (myValue + uHeightOffsetScale.x) * uHeightOffsetScale.y);
// construct texture gradient vector
vec2 valGrad = vec2(ave24 - ave13, ave34 - ave12);
// use derivative operators to construct Jacobian matrix, scaled by image size
vec2 xgrad = uImageDim.st * dFdx(gl_TexCoord[0].st);
vec2 ygrad = uImageDim.st * dFdy(gl_TexCoord[0].st);
// mat2 constructor is column major, so xgrad and ygrad are its columns
mat2 jacob = mat2(xgrad, ygrad);
// transform texture gradient by Jacobian into screen space
valGrad = jacob * valGrad;
// use gradient length as line width to perform antialiased blend of line color with height color ramp
float dist = abs(myValue - uContourVal);
gl_FragColor = mix(vec4(1.0, 0.0, 0.0, 1.0), vec4(baseColor, 1.0), smoothstep(0.0, length(valGrad), dist));
}
The first thing the shader has to do is identify which quadrant of the source texel the outgoing fragment is located in. It takes the fragment's texture coordinate, and multiplies that by the texture dimensions, so that each texel occupies [m, m+1) x [n, n+1), centered at (m+0.5, n+0.5). Calling fract() drops the integer offset, clamping it to [0, 1) x [0, 1). Then offsetting this by -0.5 in both dimensions allows us to use the sign() function to identify the texel quadrant. This is then scaled by the size of one pixel in each dimension, which gives us the offsets needed to get the neighboring pixels to perform bilinear interpolation and gradient estimates with.
It then uses those offsets, to get the 2x2 neighborhood of pixels surrounding the fragment. These pixels are then pairwise interpolated with the mix() function, using the texel quadrant coordinates as the interpolation weight. The top and bottom edges are then interpolated in the vertical direction to give us the bilinearly interpolated value for the fragment. While this could be done for us automatically by setting the IDLgrImage to have its INTERPOLATION property set to 1, we still need those linearly interpolated edge values to approximate the surface gradient for us.
We use the fragment's interpolated height value to calculate the color ramp, which is stored in the baseColor variable. Next comes the linear algebra lesson, where we construct the surface gradient, which is in the coordinate system of the model. We need to transform this to screen coordinates, so we construct the Jacobian matrix to do just this. There are the very powerful dFdx and dFdy functions in GLSL, which can calculate the rate of change of any varying variable in screen coordinates. Here we calculate the rates of change of the texture coordinates, and again scale them by the texture dimensions so that they are change per pixel. Once we've transformed the gradient to screen coordinates, we use that as the upper threshold to the smoothstep function, which give us a cubic approximation to the step function. The value tested by smoothstep is the absolute value of the different between the interpolated value and the isocontour value, and the output is used as the interpolation factor between the baseColor value and bright red, which is the color of the isocontour line.
When the surface is viewed from above, the contour lines will look perfect, but at oblique angles you can see some smearing artifacts of the contour line on the vertical surfaces. This is most likely due to the way the dFdx and dFdy functions do their calculations. They operate on 2x2 neighborhoods of screen pixels in order to calculate finite differences, so there is no guarantee that all four pixels in a given neighborhood are from the same triangle, or even are from the same part of the surface when self occlusion occurs. So the Jacobian matrix is an approximation, not an exact representation of the partial derivatives, and we get these artifacts. But the complexity of this approach is based on the number of pixels on the screen, not the number of datapoints or triangles being rendered, so it can be used for fast evaluation of the isocontours when exploring new data.