February 2006 Archives
I needed to plot the intensity profile of some test images for my thesis. It's really useful to be able to visualise 2D grayscale image data by treating the intensity as a height-field and displaying it in 3D. There's a variety of ways to do this, but I wanted something that would produce good printed results, so EPS was the best output option.
So here's a recipe for using GNU R to produce such a plot.
First, generate your image data and crop it appropriately. For example, I used ITK to generate a Gaussian field.
Then convert it to the pgm file format. The most excellent ImageMagick tools make this a snap (note the '%' is the shell prompt):
% convert -compress none inputfile.tif inputfile.pgm
The compress option ensures the PGM file is written in ASCII format, for easy reading in.
Then, fire up R and do the following (note the '>' character is the default R prompt):
> sz <- scan("gaussian.pgm", what=integer(0), skip=1, comment.char="#", nmax=2)
Read 2 items
> g <- scan("gaussian.pgm", what=integer(0), skip=3, comment.char="#")
Read 65536 items
> gi <- matrix(g, sz[1], sz[2])
> rm(g)
The first line reads the image size from the PGM header into the variable sz, the second reads in the actual data as an array. The third line creates a matrix from the data, with the appropriate size. The last line removes the temporary array data.
Now you have a matrix object with your image data in it, you can do all sorts of groovy things. Such as:
> image(gi)
> contour(gi, add=TRUE)
[img_assist|fid=13|thumb=1|alt=Gaussian 2D Image Contour plot]
will display it as an image, with a scalar mapping function (which you can customize), along with contour lines showing the isolevels. But more interesting is this:
> persp(gi, theta=25, phi=30)
This will produce a 3D contour plot in perspective (the theta and phi parameters rotate the viewing angle).
[img_assist|fid=16|thumb=1|alt=Gaussian 2D Image Perspective plot]
You can also do other handy stuff, like pick out a scanline and do a 2D intensity plot (or a series).
