Fun with IDLgrShader

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
  ; load Maroon Bells binary data, using dimensions and data type from index.txt
  m = 350
  n = 450
  heights = Float(Read_Binary(FilePath('surface.dat', ROOT=!Dir, $
                                       SUBDIR=['examples', 'data']), $
                              DATA_TYPE=2, DATA_DIMS=[m, n]))
  ; scale height data to make commensurate with horizontal scale
  heights *= 0.1
  minHeight = Min(heights, MAX=maxHeight)
  meanHeight = (minHeight + maxHeight) * 0.5
  ; Basic scene graph
  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
  ; build model and add to view
  oModel = IDLgrModel()
  oView->Add, oModel
  ; create surface and add to model
  oSurf = IDLgrSurface(DATAX=Findgen(m), DATAY=Findgen(n), $
                       DATAZ=heights, STYLE=2, USE_TRIANGLES=1)
  oModel->Add, oSurf
  ; create image to use as texture in fragment shader
  oHeightImage = IDLgrImage(DATA=heights, INTERNAL_DATA_TYPE=3, INTERPOLATE=0)
  textureSize = [1.0*pow2Ceil(m), 1.0*pow2Ceil(n)]
  ; create shader
  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
  ; create mouse event observer
  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.