Construction of 3D Surfaces from Image Regions of Interest

QUESTION: I have 50 MRI images in TIFF files in a sub-directory. They represent serial sections of a head. I'd like to open these files, draw a region of interest (ROI) around the interesting feature on each of them, and then re-construct a 3D surface from the ROIs that I can view as a polygon mesh. Can you explain how to do this in IDL?

ANSWER: I hope so. :-)

You probably don't have 50 medical image TIFF files laying around in a sub-directory, but if you would like to follow along with this example, use this IDL program, named MakeMRIImageFiles, to construct the files. You may give it the name of a directory to use, or the program will simply create a sub-directory in your current directory named "mri_images" if you don't pass it a directory name.

   IDL> MakeMRIImageFiles, 'mri_images'

The files that were just created are named image01.tif, image02.tif, through image50.tif. They are serial MRI images of a head. Each image is a 2D array of 80 columns and 100 rows. The idea here is to read these files into IDL, draw or create a region of interest on each image, and then use the ROIs to construct a 3D representation of the surface as a polygon mesh or isosurface.

Suppose we just wanted to read the data files, one at a time, and create a volumetric data set from them. With such a data set, we could create an isosurface (a surface of constant value). First, we need a list of the data files in the correct order. We could obtain such a list by moving into the image subdirectory and finding all the TIFF files with Findfile. The code might look like this:

   CD, 'mri_images', Current=thisDirectory
   files = FindFile('*.tif')

You might expect the files would be ordered correctly in the files variable. But your expectations could easily be wrong. A good rule of thumb is to sort your files yourself if you expect them to be sorted. The IDL Sort command can be used, but many IDL users are surprised to learn that doesn't do a good job of sorting on all platforms. Almost everyone I know uses the BSORT routine from the NASA Astronomy IDL Library, since this gives consistent results on all platforms.

   index = BSort(files, sortedFiles)

Next, we simply read each image data file in a loop and store it in a 3D volume array. Each image is 80 columns by 100 rows, and there are 50 of them, so the code looks like this:

   volume = BytArr(80,100,50)
   FOR j=0,49 DO BEGIN
      image = Read_Tiff(sortedFiles[j])
      volume[0, 0, j] = image
   ENDFOR
   CD, thisDirectory

Note the use of the zeros in the left-hand expression. This is several orders of magnitude faster than using asterisks to indicate the columns and rows of the volume array. This is a trick worth keeping in mind while you are writing IDL programs!

Finally, let's calculate an isosurface for this volume. We will set a 3D coordinate system up with the Scale3 command. The PolyShade command will use this 3D coordinate system to orient the polygons in the final rendering step. The Shade_Volume command will compute the polygons, along with their vertices, that make up the isosurface at the pixel value of 50. The Low keyword indicates that we will display the low side of the contour surface. In other words, the surface will enclose the volume values greater than 50. Finally, the PolyShade command constructs the shaded polygon surface.

   Window, 6, XSize=300, YSize=300, Title='MRI Head Data IsoSurface'
   Scale3, XRange=[0,80], YRange=[0,100], ZRange=[0,50], AZ=-150
   Shade_Volume, volume, 50, vertices, polygons, /Low
   theHead = PolyShade(vertices, polygons, /T3D)
   TV, theHead

The results are seen in the figure below.

Creating an IsoSurface from a Stack of ROIs

The result is not so bad. Can we get a similar result if we use regions of interest (ROIs) and stack those together? Let's find out.

First, let's have a look at what the image slices themselves. Let's put them all in one window. Type this:

   Window, Title='Original Image Slices', XSize=800, YSize=500, 1, XPos=0, YPos=0
   LoadCT, 0
   Device, Decomposed=0
   FOR j=0,49 DO TV, volume[*,*,j], j

The results are seen in the figure below.

What I would like to do is create an outline, or region of interest, about the light area in each of these 50 images. I can do this in various ways. For example, I could use the IDL tool XROI to carefully draw a free-hand line around the light area in each image. When I had completely outlined the region of interest, I could save the outline as an ROI. For example, XROI will allow me to save the outline as an IDLgrROI object. For example, I might draw the outline like this on image 15 in my sample:

   XROI, volume[*,*,14], Regions_Out=thisROI, /Block

Notice the Block keyword. This is essential for getting the ROI out of the program! The XROI program looks like this when I am about half-way finished drawing the ROI.

When I am finished drawing the ROI, I exit the XROI tool and retrieve the IDLgrROI object. Because XROI is an object graphics program, the ROI is reported in XYZ coordinates. The Z coordinate, of course, is 0, since we are drawing on a flat XY plane. But we wish to create an ROI stack, eventually, so we have to artificially assign a Z value to our ROI. We can obtain the ROI data, and modify the Z values to reflect the image's place in the stack, like this:

   thisROI -> GetProperty, Data=theROIData
   numElements = N_Elements(theROIData[0,*])
   theROIData[2,*] = Replicate(14, numElements)

We would do this for each image in the stack, capture an ROI for each image, and then store each ROI in an IDLgrROIGroup object. Doing so would be tedious, but necessary for some types of scientific work (as you well know).

In the interest of time, I prefer to do this another way. I'm simply going to brute force my way through this. I will set any value in our volume data set less than 50 to 0. Then I will set any pixel with a value greater than 50 to 255, perform a very heavy dilation operation (followed by an erosion) to close any reasonable small connections in the bright area of the image, and find the outline of the result with the Find_Boundary program from the Coyote Library. The code will look like this:

      ; Create some windows for displaying the output and creating the ROIs.

   Window, Title='ROI of Each Image', XSize=800, YSize=500, 2, XPos=0, YPos=530
   Window, XSize=800, YSize=500, /Pixmap, 3
   Window, XSize=80, YSize=100, 4, /Pixmap

      ; Create a binary volume image with all pixel values greater than 50 turned on.

   index = Where(volume LT 50, count)
   IF count GT 0 THEN volume[index] = 0
   volume = (volume GT 50) * 255B

      ; Create an ROI group object for holding all the ROI objects we are going to create.
      ; (The ROI objects will be stored in an object array to facilitate object clean-up later.)

    roiGroup = Obj_New('IDLgrROIGroup')
    theROIs = ObjArr(50)
   
      ; For each image slice, extract it, find its boundary (ROI), save the ROI in an
      ; IDLgrROI object, and add each ROI object to an IDLgrROIGroup object.

   FOR j=0,49 DO BEGIN

	 ; Extract the image.

      image = Reform(volume[*,*,j])

	 ; Close open areas with a healthy dilation kernel.

      kernel = Replicate(1, 5, 5)
      image = Morph_Close(image, kernel)
      image = BytScl(Temporary(image))

         ; Draw in a large window so the coordinates are right for viewing.
         ; Later, draw into a small window so the coordinates are right for
         ; the isosurface I am trying to build.

      WSet, 3
      Erase
      TV, image, j
      i = Where(TVRD() GT 0)
      b = Find_Boundary(i, XSize=800, YSize=500)

         ; Draw the image (temporary) in the ROI window. Wait one second,
         ; then replace the image with the ROI boundary. This gives the user
         ; some sense of what is going on.

      WSet,2
      TV, image, j
      Wait, 1
      TV, BytArr(80,100), j
      PlotS, b, /Device

         ; That was for show. This is for dough. :-)

      WSet, 4
      Erase
      TV, image
      i = Where(TVRD() GT 0)
      b = Find_Boundary(i, XSize=80, YSize=100)

	 ; Create the ROI object. Take care to create proper Z values in the ROI. Each
         ; Z value is the position of the slice in the stack.

      theROIs[j] = Obj_New('IDLgrROI', b[0,*], b[1,*], Replicate(j, N_Elements(b[0,*]) ) )
      roiGroup -> Add, theROIs[j]

   ENDFOR

You can see the collection of ROIs in the figure below.

Creating the Surface Mesh

The next step is to create a surface mesh from this stack of ROIs. Fortunately, we don't have to do that ourselves, because the hard work has been done for us. The IDLgrROIGroup object has a ComputeMesh method that is written specifically for this purpose. The code looks like this. The variables vertices and polygons contain a list of the polygon vertices, and the polygons that make up the mesh surface. This is exactly the information we need to create a shaded polygon surface. (See PolyShade above.)

numTriangles = roiGroup -> ComputeMesh(vertices, polygons)

We have a couple of choices. We can make the shaded polygon surface in direct graphics (similar to the isosurface above), or in object graphics, which would make the surface easier to rotate in 3D space. Let's do both and compare the results.

To create the direct graphics representation, we need to see up a 3D coordinate system, as before. The code will look like this:

   Window, 5, XSize=300, YSize=300, Title='The Direct Graphics Polygon Mesh Surface'
   T3D, /RESET
   Scale3, XRange=[0,80], YRange=[0,100], ZRange=[0,50], AZ=-30, AX=-30
   theHead = PolyShade(vertices, polygons, /T3D)
   TV, theHead

The result is shown in the figure below.

Notice the quality is not quite as high as the original isosurface, which you might expect. But a reasonable result, I think. There do seem to be some artifacts from the brute force way I extracted the ROIs. I think we could do better if we took more time with the ROI step.

But how about viewing the same data in object graphics? We will need to create a polygon object, add that to a model, scale and translate the model into a reasonable 3D view, and then display the model with XObjView. The code looks like this:

   theModel = Obj_New('IDLgrModel')
   theModel -> Scale, 1./80, 1./100, 1./50.
   theModel -> Translate, -0.5, -0.5, -0.5
   theModel -> Rotate, [1,0,0], -90
   theModel -> Rotate, [0,1,0], 210
   theModel -> Rotate, [1,0,0], 30
   numTriangles = roiGroup -> ComputeMesh(vertices, polygons)
   thePolygon = Obj_New('IDLgrPolygon', vertices, Polygons=polygons, Color=[128,128,128], Shading=1)
   theModel -> Add, thePolygon 
   XObjView, theModel, /Block, Scale=1.0, XSize=300, YSize=300

The results are shown in the figure below. The head in the XObjView program can be rotated in 3D space by dragging the object with your mouse.

To complete your program, be sure to clean up all your objects and pixmaps. Here is the code:

   Obj_Destroy, [theROIs, roiGroup, theModel , thePolygon]
   WDelete, 3
   WDelete, 4

If you like, you can download the MRI_MESH program that has all the above commands in it.

Like this program? Want to encourage this sort of thing? Now you can make it happen! Just click the PayPal button below to contribute whatever amount you like. Coyote thanks you!

[Return to IDL Programming Tips]