1    Ray Marching

In this activity,  you will march rays through a VTK render to implement various transfer and composite functions.

The outline of this activity will be:  1.  Implementing a grayscale transfer function with the max

intensity operator 2. Implementing an alpha transfer function with the over operator Starter Code:

#grade (DO NOT DELETE THIS LINE)

import random, os, gzip, numpy, vtk

from vtk.util import numpy_support as VN

from PIL import Image

VTK_FILE = "./HeadCT.vtk"

def compositeKeepA(a,b):

return a

def getVolumeScalars(volume):

scalars = VN.vtk_to_numpy(volume.GetPointData().GetArray("scalars")) return scalars.reshape(volume.GetDimensions(),order="F")

def volumeRender(displaySlice, volume, transfer, composite, yval=-1):

#get the dimensions of the volume array

xsize, ysize, zsize = volume.GetDimensions()

#get the min and max values of the scalars for volume

srange             = volume.GetScalarRange()

#get the scalars themselves

scalars           = getVolumeScalars(volume)

#determine if we're just taking a single y slice or compositing it somehow

ymin = ymax = yval

if yval == -1:

ymin = 0

ymax = ysize

#loop over every element on the z and x axes

for z in range(zsize):

for x in range(xsize):

#do transfer and composite for y values across the entire range

a = transfer(scalars[x,ymin,z],srange)

for y in range(ymin,ymax):

a = composite(a,transfer(scalars[x,y,z],srange))

ra,ga,ba,aa = a

#update displaySlice appropriately

displaySlice.SetScalarComponentFromFloat(x, z, 0, 0, ra*aa)

displaySlice.SetScalarComponentFromFloat(x, z, 0, 1, ga*aa)

displaySlice.SetScalarComponentFromFloat(x, z, 0, 2, ba*aa)

def saveSliceToImage(displaySlice,fname):

xs,zs,_ = displaySlice.GetDimensions()

scalars = VN.vtk_to_numpy(displaySlice.GetPointData().GetArray(0))

scalars = (255*scalars).astype(numpy.uint8)

scalars = scalars.reshape((xs,zs,3),order="F")

scalars = numpy.rot90(scalars,1)

im     = Image.fromarray(scalars)

im.save(fname)

display(im)

def decompressVtkData(ifname,ofname):

if not os.path.exists(ofname):

if os.path.exists(ifname):

with gzip.GzipFile(ifname, 'rb ') as fin:

with open(ofname, 'wb') as fout:

fout.write(fin.read())

#transfer = transfer function to pass to volumeRender()

#composite = composite function to pass to volumeRender()

#ymin = y value to pass to volumeRender

def render(transfer,composite,fname,yval=-1):

decompressVtkData(VTK_FILE+".gz",VTK_FILE)

reader = vtk.vtkStructuredPointsReader()

reader.SetFileName(VTK_FILE)

reader.Update()

volume = reader.GetOutput()

xsize, ysize, zsize = volume.GetDimensions()

displaySlice = vtk.vtkImageData()

displaySlice.SetExtent(0,xsize-1,0,zsize-1,0,0);

displaySlice.AllocateScalars(vtk.VTK_FLOAT,3);

volumeRender(displaySlice,volume,transfer,composite,yval)

displaySlice.Modified()

 

2    Part 1:  Maximum intensity projection

Implement transferGrayscale():

•  function should return an  (r,g,b,a) tuple

•  r, g, and b components should be equal to the linear interpolation of  scalar between the min and max values in srange, normalized to (0,1)

•  a component should always be 1

[ ]: #grade (DO NOT DELETE THIS LINE)

def transferGrayscale(scalar, srange):

###BEGIN STUDENT SOLUTION

pass

###END REFERENCE SOLUTION

Test transferGrayscale()

[ ]: #unit tests for transferGrayscale()

r        = [4,7]

scalars = [5,6,5.5,6.5]

grays   = [2/6,4/6,3/6,5/6]

for i in range(len(scalars)):

numpy .testing .assert_allclose(transferGrayscale(scalars[i],r)[0],grays[i],␣ ,→rtol=1e-2)

numpy .testing .assert_allclose(transferGrayscale(scalars[i],r)[1],grays[i],␣ ,→rtol=1e-2)

numpy .testing .assert_allclose(transferGrayscale(scalars[i],r)[2],grays[i],␣ ,→rtol=1e-2)

numpy .testing .assert_allclose(transferGrayscale(scalars[i],r)[3],1,␣ ,→rtol=1e-2)

Test render using transferGrayscale()                                                    [ ]: render(transferGrayscale,compositeKeepA,"test .png",200)                            

Your image should look like this:

Implement compositeMax():

•  a and b are both rgba tuples of the form  (red,green,blue,alpha)

•  return value should be the max of a and b using standard ordering operators (e.g., compare red, then green, then blue, then alpha)

[ ]: #grade (DO NOT DELETE THIS LINE)

def compositeMax(a,b):

###BEGIN STUDENT SOLUTION

pass

###END REFERENCE SOLUTION

Test compositeMax()

[ ]: #unit tests for compositeMax()

for i in range(1000):

a = tuple([random .random() for _ in range(4)])

b = tuple([random .random() for _ in range(4)])

if a > b:

assert(compositeMax(a,b) == a)

else:

assert(compositeMax(a,b) == b)

Test render using compositeMax()                                                         [ ]: render(transferGrayscale,compositeMax,"part1 .png")                                 

Your image should look like this:

3    Part 2:  “X-ray” projection

Implement transferAlpha():

•  function should return an  (r,g,b,a) tuple

•  r, g, and b components should always be  1

•  a component should be equal to the linear interpolation of scalar between the min and max values in srange, normalized to (0,1), and multiplied by the constant 0.0125

[ ]: #grade (DO NOT DELETE THIS LINE)

def transferAlpha(scalar, srange):

###BEGIN STUDENT SOLUTION

pass

###END REFERENCE SOLUTION

Test transferAlpha()

[ ]: #unit tests for transferAlpha()

r        = [4,7]

scalars = [5,6,5.5,6.5]

alphas  = [0.004167,0.008333,0.00625,0.010417]

for i in range(len(scalars)):

numpy .testing .assert_allclose(transferAlpha(scalars[i],r)[0],1, rtol=1e-2)

numpy .testing .assert_allclose(transferAlpha(scalars[i],r)[1],1, rtol=1e-2)

numpy .testing .assert_allclose(transferAlpha(scalars[i],r)[2],1, rtol=1e-2)

numpy .testing .assert_allclose(transferAlpha(scalars[i],r)[3],alphas[i],␣ ,→rtol=1e-2)

Implement compositeOver():