Fun with, inter alia, ENVI-IDL

Pansharpening on the Google Appengine

2012-01-22T14:03:00.001+01:00

One of the more elegant panchromatic sharpening techniques in multispectral image processing makes use of the discrete wavelet transform (DWT). I describe it in some detail in my book and provide an ENVI/IDL script for it on my software site. I also mentioned running it on the cloud with small images in an earlier blog. Using PiCloud as a backend for the App Engine, it wasn't difficult to get the algorithm to run on much larger images on my GAE App. The picture below shows a small subset of an Ikonos image with the upsampled four-meter resolution multispectral image at the center. Bands 1,2 and 3 are shown as an RGB-composite. The DWT pansharpened scene, calculated with the App using the one-meter resolution panchromatic image included with the Ikonos data, is shown on the left. For comparison, the result of applying ENVI's built-in Gram-Schmidt algorithm is shown to the right. Notice that the colors match better for DWT.

These are just 400x400 pixel cutouts but the full pansharpened scene is 2000x2000 pixels, about the largest size the App can handle just now.

Multispectral Image Analysis on the Cloud (again)

2011-12-22T13:38:00.000+01:00

Well, I've made a few fairly big jumps forward, as I continue to " muck about" on the Google Appengine. The biggest jump was to get Numpy to work, at long last. Another was to make extensive use of PiCloud as a background service, another to use SvgFig to render 2-D plots on the client side, and yet another to "mashup" with OpenLayers to show image footprints on the world map. Picloud, as I mentioned before, provides Python wrappers for GDAL, but also has the Python wrapper for the computer vision package OpenCV pre-installed. This opens up all sorts of nice tools. So far, I've only implemented the Canny edge detector, but there will be more to come.

I have not enabled billing in the GAE, so I just get 28 "instance hours" a day, each with 600 Mhz and 128 MB. Also my free storage quota is 5 GB (for me and all users together). I'm getting a much better environment on PiCloud (1.2 GHz, 20 hours per month free CPU time), but I can only send 5 MB payloads to PiCloud from the Appengine. Based on past experience, these limits will continue to get better as time goes by and so will, I hope, the App.

Getting all this stuff to work (more or less) would have been impossible for me without the very friendly help from the Appengine team, the people at PiCloud and the altruistic contributors to StackOverflow. If any of you guys are reading this blog, then

Many thanks, Season's Greetings, and I'll be bugging you again in the New Year!

Mort

Scatterplot Matching

2011-12-09T14:04:00.001+01:00

S. J. Maas and N. Rajan (2010) Remote Sensing 2, 1644-1661, describe a method to radiometrically normalize the RED and NIR bands of multitemporal Landsat TM images which uses the characteristic shape of NIR vs RED scatterplots. The so-called bare soil line (BSL) and full canopy point (FCP) derived from the scatterplots of target and reference image are used to re-scale the RED and NIR bands of the target to match those of the reference. The procedure works with other platforms too, and is illustrated below for two ASTER images (the RED band is VNIR band 2 and the NIR band is VNIR band 3):

The red lines are the derived BSLs and the red crosses mark the FCPs in the reference (left) and target (center) scatterplots. The scatterplot of the normalized target image (right) closely matches that of the reference image (left). The images below are RGB images of (R=band 2, G=band 3, B = band2) for reference (left), target (center) and normalized target (right) on the scale (0-255):

Notice that the normalized target colors are very similar to those of the reference. We can use the IR-MAD normalization method to check the result. The plots below show the reference vs. target regression lines using pseudo-invariant pixels identified by IR-MAD before (top row) and after normalization (bottom row).

The slope and intercepts are as follows:

top left (band 2): slope 1.49 intercept 11.7 top right (band 3) slope 1.40 intercept -3.44
bottom left (band 2) slope 0.92 intercept 17.6 top right (band 3) slope 1.02 intercept 1.07

So the correction seems to work somewhat better for the NIR band 3. Here is an ENVI/IDL extension for scatterplot normalization.

Still waiting for Numpy, and a bug in GDAL

2011-11-15T12:41:00.001+01:00

The latest Python software development kit for Google's App Engine (GAE) offers experimental support for Python 2.7, concurrent programming and a number of C-libraries, one of which is Numpy (Numerical Python). I've been waiting for Numpy for a long time in order to make my image analysis Webb App efficient, or more accurately, usable. Unfortunately there are problems and I still can't import Numpy into my GAE scripts. I'm hoping things will be sorted out in the next release.

Another source of trial and tribulation is GDAL (Geospatial Data Abstraction Library). I'm using it in my App to read and write raster data from different standard formats. To do this, I need to shove the imagery back and forth between GAE and PiCloud since GAE has no GDAL support. That's complicated enough, but it turns out that GDAL doesn't handle georeferenced images which are rotated away form north-up correctly. The rotation is ignored. A lot of satellite imagery is supplied path-oriented (oriented along the satellite trajectory) so that isn't very nice. So far no advice from the GDAL developers' forum or from StackOverflow either.

Free quotas on PiCloud

2011-10-26T12:27:00.000+02:00

PiCloud announced yesterday that it will provide users with 20 free c1 core hours per month. At their very reasonable prices that's only $1.00 worth of CPU time, but it means that you don't have to enable billing at all to be able to experiment with and/or use the service indefinitely!

A Python script for Kernel PCA on the GPU

2011-10-15T12:02:00.000+02:00

Just to wrap up the discussion in my last post but one: Although the pyCublas wrapper for matrix multiplication on the GPU doesn't bring much of a speedup relative to the CPU when I use Numpy together with the Math Kernel Library on 64bit Windows, it may well do so in other configurations. I only have the one machine with a CUDA-enabled graphics card, so I can't experiment. You can now download the code for both versions (with and without CUDA) from my software page.

By the way, pyCublas won't work for CUDA 4.x because nvidia has changed the CUBLAS API completely.

PiCloud adds GDAL

2011-10-11T16:39:00.000+02:00

For anyone using (or thinking of using) the Cloud Sevice PiCloud for image analysis, here is great news:

Hi Mort,

We decided to add GDAL to our systems by default. If you run a job now that uses the gdal Python library (import gdal), it should work. Let us know how it goes!

Ken Elkabany
PiCloud, Inc.

Cool!

Numpy is FAST!

2011-10-08T21:22:00.001+02:00

After upgrading my 32bit Python(XY) package to version 2.7, I discovered to my horror that Eclipse/Pydev was no longer there. It is actually the only part of the package that I regularly use. So I decided to "make nails with heads on them" as we say in German, and do my own 64bit Python 2.7 installation. In doing so, I discovered a great deal. To begin with, for anyone interested, here's how I went about it on my Windows machine:

- installed the 64bit Java Runtime (this is what Eclipse runs on)
- intalled Eclipse 3.7.1 (win32-x86_64)
- installed pydev onto Eclipse
- installed 64bit Python 2.7
- installed Numpy 1.6 (win-amd64-Py2.7)
- registered the Python interpreter with Pydev.

That's it.

Then it occured to me that I could now make use of the 64bit CUDA drivers/toolkit with my new Python installation. Not wanting to write my own CUDA kernels, there seemed to be three possibilities available to me:
PyCuda
SciKits.Cuda
PyCublas
My experience with using CUDA in IDL has been that the only really substantial speedups involve matrix multiplication on the GPU. This means that the CUBLAS library should be interfaced to whatever script language one is using in order to take full advantage of NVIDIA optimized BLAS (Basic Linear Algebra Subprograms) code, as is the case, for example, with GPULib in IDL. This eliminated PyCuda, which explicitly excludes CUBLAS. SciKits.Cuda has Python wrappers for CUBLAS, but doesn't support the Windows platform (sheesh). That left PyCublas, written by Derek Anderson, a very simple and pretty wrapper for matrix multiplication using ctypes to wrap the CUBLAS DLL.

For later comparison with Python, here is what I get for IDL with and without GPULib when multiplying two 2048x2048 matrices together:

Welcome to GPULib 1.4.4 (Revision: 2107)Graphics card: GeForce GTS 450, compute capability: 2.1, memory: 752 MB available, 993 MB totalChecking GPU memory allocation...no errors 0.10899997 2.5430000The first number is the GPU time, the second the time required by IDL on the CPU, so the speedup is about 25x.
Now for Python:

[pycublas] time (CUBLAS): 0.171000s[pycublas] time (numpy): 68.578000s[pycublas] speedup: 401.0XWow! Python/Numpy is much slower than IDL and using CUDA really makes sense. However, I then deinstalled Numpy and replaced it with Numpy compiled against the Intel Math Kernel Library (MKL). On my machine (which has an Intel Core i5 CPU) I get

[pycublas] time (CUBLAS): 0.172000s[pycublas] time (numpy): 0.218000s[pycublas] speedup: 1.3XSo, in this configuration at least, CUDA is of little or no advantage.

RADCAL Update

2011-09-06T10:41:00.000+02:00

I just realized that I forgot to include spectral wavelengths, when available, in images normalized with RADCAL. That's important if you're normalizing in order to compare vegetation indices, for example. Fixed that now.

Numpy at last

2011-08-20T15:10:00.000+02:00

"When sorrows come, they come not single spies, but in battalions." Fortunately this rule sometimes also applies to "joys". Almost simultaneously (at least within a few days of eachother) the following happened:

1. Picloud announced its REST registration interface, which allows Picloud functions to be called from the Google App Engine (GAE). So CPU-intensive calculations can be "farmed out" to Picloud servers. The servers support numpy, scipy and virually anything else you need.

2. Google announced that GAE will soon be running on Python 2.7 and, furthermore, will include numpy on its runtime.

3. Dan Sanderson announced the initial drafts of the second edition of his book Programming Google App Engine. No doubt he will discuss the new features in detail.

So it looks like, at long last, I will be able to port my ENVI/IDL algorithms to GAE in an efficient, IDL-Matlab-type syntax rather than continuing to torture myself with slow pure Python versions.

Mathematica and Game Theory

2011-08-14T17:14:00.000+02:00

Inter alia = "among other things", and it's about time that I made this change to my Blog title.

Quite a few years ago I wrote a book on algorithmic game theory called Resolving Conflicts with Mathematica. It was built around the (at the time) most recent version of Wolfram Inc.'s Mathematica, namely Version 5.2. The software package on my Home Page was, till now, only compatible with Version 5.2. I think I've fixed that. I say "think" because I stopped upgading at Version 6. I believe, though, that Versions 7 and 8 are backward compatible to Version 6 (but not to 5.2).

Image Analysis on the Cloud: New App Deployed

2011-08-04T15:51:00.000+02:00

I just deployed the new, improved multispectral image analysis app to Google's cloud service (Google App Engine), see my previous blog. It still has a lot of quirks, one reason being that memory size restrictions apply on the cloud which are neither simulated on the development server nor, as far as I can see, anywhere documented. Anyway, it's fun.

Image Analysis on the Cloud (Progress report)

2011-07-24T22:06:00.004+02:00

I just spent a whole weekend modifying my Google App Engine App. I haven't released the new version yet, prior to more testing. I now use the so-called Blobstore to store the users' image data, so that larger images can in principle be processed. The image formats correspond more closely to those generated by my ENVI/IDL routines. For example, the MAD variates are now stored in floating point format with a chi-square band appended. I'd previously been byte-stretching them to save space.

I've limited the upload size to 8MB uncompressed multispectral images in byte format (formerly 1MB), allowing for up to 32MB processed images in floating point format (for instance the result of a principal components, Fourier or IR-MAD transformation). There are two main reasons for the limitation:

1. My free Blobstore quota is 1GB, and a single user can fill up to 5 image slots. So in the worst case, if every user fills her slots to the maximum 5 x 32MB = 160MB, I can only serve 1000 mod 160 = 6 users. (Unless I start paying Google for more resources, which I don't intend to do.)

2. There is still no Numpy, so the algorithms are very slow and would be hopelessly slow for larger images. IR-MAD runs on the App Engine task queue and calls itself recursively for each iteration. Since the task response time deadline is 10 minutes, this allows for up to 10 minutes per iteration. With Numpy this would be plenty, but in pure Python the itertations are very slow. For two 1000 x 1000 x 6 byte LANDSAT images, each iteration takes around 4 minutes. I have the App send the user an email notification when the calculation finishes.

I have enquired with Picloud about the possibility of "farming out" computationally intensive functions from the Google App Engine. Since Picloud has Numpy/Scipy, this would solve my problem. The Picloud team is looking into it and I may be asked to beta-test. Of course, they eventually want to see cash, too.

IDLe moments

2011-07-02T21:27:00.009+02:00

With nothing better to do this evening I decided to add some gadgets to the blog design. The "SEARCH THIS BLOG" box on the right is probably useful now that there are almost 60 pages of my blatherings. I also took the opportunity to update my photo with one that more closely corresponds to my not inconsiderable age.

Still waiting for the Google people to enable scientific computing in Python on their fantastic App Engine platform, which someone went so far as to call an "operating system for the internet". My begging on the App Engine/Python discussion group the other day drew a friendly response from Greg D'Alesandre, App Engine Senior Product Manager:

Hi Mort,

Your time (and, um, begging) isn't wasted, we hear you and we know this is a very important ask and are prioritizing it with our other work. Unfortunately at this time I don't have anything to announce about it...

So maybe one day ...

Modern IDL

2011-06-24T22:04:00.005+02:00

I ordered Michael Galloy's new book Modern IDL last week. The hardcopy hasn't arrived yet (I'm in Germany) but Mike sent me a link to the PDF version and I've browsed through it. In my own work I use IDL mainly for algorithm implementation and leave most of the graphics user interface, data I/O, map transforms and other pizazz to ENVI. But I've always wanted a thorough, concise, up-to-date overview of the the IDL language and its vast capabilities. This is exactly what Mike's book provides in 464 very informative pages. Can't wait to have the "real thing" on my bookshelf. Highly recommended!

Scatterplots for the Google Cloud App

2011-06-23T10:20:00.006+02:00

If you run automatic radiometric normalization on my GAE App you can now plot the orthogonal regression lines next to the table of regression coefficients (click to enlarge):

The App uses the HTML5 CANVAS object to draw the plots, so it should run on any up-to-date browser, with the exception (of course) of IE9.

IR-MAD on the Google App Engine

2011-05-17T21:32:00.004+02:00

I just taught myself how to use background tasks in the App Engine (a lot easier than I thought). The IR-MAD algorithm now will run completely in background in my App, so you no longer have to iterate manually. Things are happening at last with the GAE, see John K. Waters' blog. I'm still hoping that Numpy will be made available, and if it is, then the sky's the limit.

iMAD on Mac OS

2011-05-09T10:11:00.005+02:00

Gillian Galford from Woods Hole Research Center figured out how to build the PROVMEANS DLM on a Mac. This code will place the library PROVMEANS.SO in the default DLM path (and it couldn't be simpler):

pro make_provmeans

PRINT, !MAKE_DLL.COMPILER_NAME

INDIR = '/Applications/itt/canty/auxil'
OUTDIR = $
'/Applications/itt/idl/idl80/bin/bin.darwin.x86_64'

MAKE_DLL, 'provmeans', 'IDL_Load', $
INPUT_DIRECTORY=INDIR, $
OUTPUT_DIRECTORY=OUTDIR, /SHOW_ALL_OUTPUT, $
compile_directory=outdir

end

After running that, just copy the description file PROVMEANS.DLM from my auxil package to the same OUTDIR and you're all set. (The DLM is used by iMAD for calculating weighted means and variance-covariance matrices.)

Python script for kernel PCA

2011-04-12T20:44:00.002+02:00

The script has pretty much the same functionality as the ENVI/IDL program described in my book, except that it doesn't take advantage of CUDA. This is not because Python has no wrappers for CUDA, but because I'm programming on x64 Windows with a 64-bit driver/CUDA toolkit. Being a green novice at Python, I have no idea how to set up my own scientific Python computing environment, so I'm working with the (very nice) Pythonxy distribution. It has everything from Numpy to Scipy to GDAL to Eclipse/Pydev but is, alas, 32 bit.

The MADMAN Rides Again

2011-03-29T11:34:00.007+02:00

I've just finished modifiying the Python scripts for IR-Mad and radiometric normalization to allow for spatial and spectral subsetting as well as for radiometric normalization of large scenes. I have now included their documentation in the MAD MAN which originally described the ENVI/IDL extensions only.

Radiometric Normalization with Python

2011-03-19T17:42:00.016+01:00

(Maybe I should start calling this blog "Fun Without ENVI-IDL".) I've just improved the Python RADCAL script to allow a user to normalize large images. You can get both iMAD and RADCAL together with a readme here. The syntax is now as follows:

To run iMad on images fname1 and fname2, which have the same spatial and spectral dimensions and are co-registered:

python iMad.py fname1 fname2

This generates the output file

MAD[fname1-fname2]

To normalize fname2 to fname1, (optionally) setting a minimum no-change probability threshold of 0.9 (default 0.95) and then (optionally) using the regression coefficients to normalize a larger file called bigfname2 having the same spectral dimension (e.g. a full scene) to fname1
python radcal.py MAD[fname1-fname2]
-ncpt 0.9 -fsfn bigfname2
This creates the normalized output files

fname2_norm
bigfname2_norm

Several image formats are allowed and can even be mixed, see my earlier blogs. Here is an example with ENVI standard files n1small, n2small and n2, all located in the directory \imagery\gdal\ (click to enlarge):

Before you can use the scripts, you must have the right Python environment installed. I'm running Windows 7 x64 with Python(xy) Version 2.6.5.6. You can get it here, and it's childsplay to install. Just be sure to include the GDAL plugin.

Porting IR-MAD to Python (continued)

2011-02-23T12:29:00.010+01:00

I described a Python port of IR-MAD in an earlier blog. The IR-MAD algorithm is used by quite a few people to do automatic radiometric normalization, so I have added a Python script called radcal.py which can be used for normalization. In order to keep things simple (if a bit inflexible) I changed the file naming conventions slightly from the previous post:

If you have two co-registered images having the same spatial/spectral dimensions, say image1.tif in GeoTiff format and image2.img in ERDAS IMAGINE format, then you can run IR-MAD on them as follows:

python iMad.py path1/image1.tif path2/image2.img

The MAD result will be written to the file

path1/MAD(image1.tif-image2.img).tif

in other words, in the same directory and with the same format as the first image.

Now, to radiometrically normalize image2 to image1, you can type

python radcal.py path1/MAD(image1.tif-image2.img).tif [noChangeThresh]

where the optional parameter noChangeThresh is the desired minimum probability of no change for selecting the invariant pixels to be used in the normalization (default 0.95). The script will write the normalized image to

path1/image2_norm.img

thus placing it in the same directory as image1 and maintaining its original format.

If you want to normalize image1 to image2, just run IR-MAD the other way round, i.e.,

python iMad.py path2/image2.img path1/image1.tif

You can download both iMad.py and radcal.py from my software page.

Face Recognition with IDL

2011-02-07T11:44:00.026+01:00

Google's Picasa has a very impressive face recognition function. Just naming a few of the hundreds or thousands of faces that the program extracts out of a photo collection is sufficient to categorize a large fraction of the remaining faces correctly. The more unknown faces the user labels, the more exhaustive Picasa's classification becomes.

I don't have any idea how Picasa face recognition works and I expect they are not going to tell me. But I do know that one standard technique makes use of linear and nonlinear principal components analysis (PCA), subjects which I deal with in some detail in my book. The original suggestion goes back to Kirby and Sirovich (1990) and is, basically, as follows:

Given, say, M labeled faces (the training set) each consisting of N pixel intensity values, find the directions in N-dimensional "face space" which maximize the variance (spread) in faces. This is an eigenvalue problem with M solutions (for M less than N). The directions (eigenvectors) are called in this case "eigenfaces". Next, project each of the M labeled faces along the first k eigenfaces, where k is usually smaller than M. The k-component vector of projections comprises the principal components of each face. Now, given an unlabeled face which is, however, known to be one of the persons in the training set, determine its principal component vector as well. Finally, give it the label of the training set face whose principal component vector is nearest. By using the so-called kernel trick, the whole procedure can be carried out in a nonlinear face space which (perhaps) can characterize the faces more efficiently. For the mathematically minded, here is a more detailed explanation.

I wrote a short face recognition routine in IDL in the form of an object class which uses both linear and nonlinear (kernel) PCA as described above. There are several public domain face recognition databases on the net, but I decided it would be more fun to use my own, as generated by Picasa. I constructed training sets of (on average) about 130 examples drawn randomly from Picasa face images of 30 relatives and friends (so each person is represented by about 4 different faces). I resampled them all to 100x120-pixel b/w images. Here are the first 25 (ghostly and ghastly-looking) eigenfaces generated from one of the training sets:

Then I tested with (on average) about 170 new faces for the same 30 persons (so 5-6 faces per subject), counting the number of correct assignments. For linear PCA the hit rate was 13%, for kernel PCA with a Gaussian kernel it jumped to 36%. That may not sound very impressive, but the probability of getting it right by pure chance is one in thirty or about 3%. It's also interesting that nonlinear PCA does so much better a job.

SVD on the GPU Revisited

2011-01-28T10:03:00.009+01:00

A while ago I wrote an IDL DLM wrapper for CULA Tools Singular Value Decomposition called CUDA_SVD, see my previous blog. I just tested it with an inexpensive, and therefore somewhat castrated, Fermi GPU (Asus GeForce GTS 450, 1Gb, €125.00) running on a Core i5 PC with 64b Windows 7. Instead of the earlier 7-fold acceleration for a 2000x2000 matrix relative to the IDL procedure LA_SVD, the speedup is now almost 11-fold (Abscissa: Number of matrix elements, Ordinate: Speedup):

Compiling an IDL DLM on Windows x64

2011-01-24T09:30:00.009+01:00

After installing 64 bit ENVI/IDL on my new Windows 7 PC, I had to face the problem of recompiling my DLMs to the x64 platform. It turns out to be quite easy to do with Microsoft Visual C++ 2010 Express if you know which settings to change, which I didn't. After two frustrating days reading MS documentation and Visual Studio forums, I got the decisive tip from, of all places, Wikipedia. Basically you have to:

- install the full Visual C++ Express (free from MS)

- install the latest Microsoft Windows SDK (also free, includes the 32bit-64bit cross compiler called from Visual C++, which in turn runs as a 32bit application on 64bit Windows)

- set the IDL include and library paths of the Visual C++ project (the DLM) to the corresponding locations for the 64bit IDL version

- from the Visual C++ Configuration Manager, choose the x64 target platform (this becomes available only after installing the SDK)

- from the project's property sheet, change the Platform Toolset to the Windows SDK instead of the built-in v100 (this essential little point was in Wikipedia)

So the first thing I did was to make a 64-bit version of the provmeans DLM, which is required if you want to run my IR-MAD ENVI extension on a 64bit Windows OS. There are now two versions in the AUXIL.ZIP package:

provmeans.x86.dll - for 32bit Windows
provmeans.x86_64.dll - for 64bit Windows

With this naming convention IDL will automatically load the appropriate version.

My shiny new PC also has a GeForce GTS 450 graphics card so I can finally start playing around with the FERMI GPU. More on that later.

Porting IR-MAD to Python

2011-01-01T19:51:00.018+01:00

If I wait until Google App Engine allows Python scripts to use numpy/scipy I'll be old and gray. Come to think of it, I'm already old and gray. So I've done what I guess is the next best thing, namely a direct port of the ENVI/IDL IR-MAD (iMAD) script to Python/numpy/scipy. My main obstacle was how to read in different image formats, since if the user has ENVI/IDL or Matlab anyway, he or she won't need a Python version. I knew about GDAL and the GDAL Python bindings, but the documentation is, shall we say, sparse? What got me off the ground was Chris Gerrard's very nice Geoprocessing with Python programming class material.

The script requires the Python packages numpy, scipy and osgeo/gdal and assumes that the GDAL binaries are correctly installed, see the links above. I'm using Pythonxy, which has everything you need apart from GDAL.

Usage: iMad infilename1 infilename2

where infilename1 and infilename2 are fully qualified file names referring to geometrically registered, multispectral images having the same spatial and spectral dimensions. They can be in any of the following formats:

ENVI Standard e.g. image1
image1.hdr

PCI Geomatics e.g. image1.pix

ERDAS IMAGINE e.g. image1.img

ARCVIEW Raster e.g. image1.bil
image1.hdr

TIFF/GeoTIFF e.g. image1.tif

These formats can be mixed, for example the first image might be GeoTIFF and the second image ERDAS IMAGINE, but the output file will be written in the same format as that of the first image, with filename in this case

MAD(image1-image2).tif

The iMAD algorithm iterates to convergence or terminates after 50 iterations. The output file is stored in BSQ floating point format and consists of the stacked MAD variates in order of decreasing variance followed by the chi-square image, i.e., same as with my ENVI/IDL script.

Next project will be a port of automatic radiometric normalization.

Happy New Year!

Pan-sharpening on the GAE

2010-12-05T15:46:00.019+01:00

Improving the spatial resolution of multispectral (MS) images via fusion with a higher resolution panchromatic (PAN) image is a common processing technique. The goal is to sharpen the MS image while at the same time preserving its spectral properties as much as possible. One of the best methods I know to achieve this involves the discrete wavelet transform (DWT). You can perform DWT pansharpening of small images on my Google App Engine app. The routine is a port to Python of the IDL program described in my book.

The screenshot below shows, on the left, a 1 m resolution panchromatic Ikonos image and, on the right, the same image after two applications of a DWT filter.

The small square in the upper left hand corner on the right has been compressed by a factor of 4, so its spatial resolution is 4 m, the rest of the right hand image consists of the so-called wavelet coefficients, i.e., the high spatial frequency details. These allow a lossless reconstruction of the panchromatic image. We can now fuse the pan image with one of the 4 m resolution MS bands by substituting the latter for the small square and inverting the transform (there is also a normalization step required), thereby mixing the high frequency PAN details with the MS band. This can be done for all of the MS bands successively. The image below shows the original MS image (bands 1,2,3 in RGB) magnified by a factor of four (left) and the pan-sharpened image (right).

For more details see Section 5.3 in my book.

GPULib 1.4.2

2010-11-24T20:30:00.005+01:00

I recently installed ENVI 4.8/IDL 8.0 and, along with them, the latest release of GPULib. The Windows installation is now very slick and taylored to Windows nerds like me: You just run GPULib-setup.exe and say whether your CUDA graphics card supports single precision, double precision or (lucky you!) has a Fermi GPU. That's it. Then, provided you have your Nvidia driver and CUDA Toolkit correctly configured, you can open IDL, point your IDL and DLM paths to the GPULib /lib directory and start parallel programming.

With the variable objects and operator overloading introduced in this version of GPULib, GPU programming is really getting nice. The code below transfers two large arrays a and b to the GPU, matrix multiplies the transpose of a with b, returns the result to the CPU and frees the graphics memory:

a_gpu = gpuPutArr(randomu(seed,1000,1000))
b_gpu = gpuPutArr(randomu(seed,1000,1000))
c_gpu = gpuTranspose(a_gpu,LHS=a_gpu)##b_gpu
c = gpuGetArr(c_gpu)
gpuFree,[a_gpu,b_gpu,c_gpu]

I'm now prettying up my CUDA-enabled ENVI extensions and thinking of new applications.

Singular Value Decomposition on the GPU

2010-10-23T11:33:00.020+02:00

EM Photonics Inc. provide a cost-free basic version of CULA Tools which ports a number of single-precision LAPACK routines to CUDA and the GPU. One of these is SVD (singular value decomposition). I wrapped it up into a DLM called CUDA_SVD which has similar functionality to IDL's LA_SVD procedure (without the DOUBLE keyword). On my old GeForce 8600GT graphics card I see a 7-fold speedup relative to LA_SVD for a 1000x2000 matrix, see the Figure below.

You can get the source together with the Windows compiled DLL here. To run, you need a CUDA-enabled Nvidia GPU, its latest graphics driver, CUDA Toolkit 3.0 or later, the CULA 2.0 runtime and, of course, IDL.

CUDA By Example

2010-08-28T18:14:00.003+02:00

I just got hold of Sanders' and Kandrot's "CUDA By Example". It is a great introduction for beginners like me. A lot of fuzziness and confusion in my mind about what's going on has been blown away just reading the first few chapters. One small criticism: The authors state right at the beginning that CUDA C is a (relatively minor) extension of standard C. But in their examples they switch back and forth between C and C++ in a very cavalier manner. For instance, a code snippet given as part of an explanation may use C++ structures, whereas the final listing later in the chapter switches back to standard C. They even use C++ constructs in __device__ code. Is CUDA C is really CUDA C++? So I'm still a bit confused.

Cloud computing (another update)

2010-08-06T16:16:00.003+02:00

Still no sign that Google will allow Numpy to be imported from Python scripts on their App Engine servers, and hence permit efficient IDL/Matlab-like applications. Still, I find it a nice (and cheap) playground to learn web programming. I've added fast Fourier image transforms and frequency domain filtering to my App

Image Analysis on the Google Cloud

and jazzed it up with a bit of Web 2.0 stuff, i.e., Ajax.

Cloud computing (an update)

2010-06-29T16:00:00.006+02:00

Well, now it's even possible to run IR-MAD (iterated MAD) on my App. All you need is a bit of patience, see the Help Page. It has to be done this way because the App Server will time out after 30 seconds, and that's not long enough to iterate the algorithm automatically. It would be plenty of time if Numpy were available! (Are you listening, bosses of the Google App Engine?!)

A GPU "gem" that didn't make it

2010-06-19T13:08:00.006+02:00

I wrote a detailed description of my GPULib implementation of the Kernel K-Means algorithm which didn't quite make the wire for the next edition of NVIDIA's "GPU Computing Gems". The reviewers' main objection was lack of detail regarding GPU computing considerations, e.g.,

"data about how the algorithm behaves on the GPU such as # of coalesced memory reads & writes achieved, # branches, instructions/clock achieved, observed pcie bandwidth achieved with the mapping, whether the library achieved good occupancy or not ..."

This struck me as a bit silly, considering that GPULib is, among other things, just wrapping CUBLAS, which in turn is written by NVIDIA's own software engineers. Anyway, here is the paper:

Kernel K-Means Clustering of Remote Sensing Imagery with CUDA

MAD on the Cloud (cont'd)

2010-05-13T15:16:00.013+02:00

I must admit, this stuff gets addicting after a while. I now have a rather slick-looking interface to the Google App Engine application

Image Analysis on the Google Cloud

which allows you to upload (small) multispectral images in about any format, view them, display their histograms, do various spectral transformations (PCA, NDVI, Tasseled Cap), perform change detection and/or radiometric normalization with MAD, and then download the results in ENVI standard format.

Severe size restrictions on data structures and the lack of Python C-library support make it heavy going, but some things are maybe useful. For example, if you have ENVI EX or ENVI Runtime but no IDL (and hence no access to my programs), to radiometrically normalize spatial subsets of two co-registered scenes you can

- upload them to the cloud in image slots 1 (reference image) and 2 (target image)
- run MAD, putting the result in image slot 3
- calculate a change probability image from the MAD image in image slot 4
- run the normalization, overwriting image slot 2 with the normalized target
- download the normalized image from image slot 2

The screenshot below shows the radiometric normalization result page for two ASTER VNIR images. The regression coefficients could be used, eg, with ENVI's Band Math to normalize the full scenes.

This screenshot shows the reference image (center), the target image (left) and the normalized target image (right).

MAD on the Cloud

2010-04-05T15:12:00.006+02:00

After "mucking about" over the Easter weekend, my App Engine website Image Analysis on the Google Cloud now includes change detection with the MAD transformation.

Caveats:

1. Because of Datastore restrictions, uploaded images must be 1MB or less. For instance, a 400x400 pixel 6-band LANDSAT image (= 960,000 bytes) will just fit. Using the Blobstore might fix this, but I'm not sure yet.

2. This is not IR-MAD, just plain, ordinary, uniterated MAD. Iteration is my next project.

3. You can view the result in RGB, but you can't download it. This will be fixed if and when I use the Blobstore (might cost money).

4. CPU usage is high since everything is programmed in pure Python. This may cause the App Engine to throttle the service if the CPU quota (whatever it is) is exceeded:

Anyway, feedback always welcome.

Python the hard way? (Cont'd)

2010-04-02T10:13:00.004+02:00

My what a beautiful language. Here is a Python implementation of IDL's CHOLDC procedure for Cholesky decomposition:

from paida.tools import Matrix
from math import sqrt

def choldc(A):
# Cholesky-Banachiewicz algorithm
L = A - A
for i in range(len(L)):
for j in range(i):
sum = 0.0
for k in range(j):
sum += L[i,k]*L[j,k]
L[i,j] = (A[i,j]-sum)/L[j,j]
sum = 0.0
for k in range(i):
sum += L[i,k]*L[i,k]
L[i,i] = sqrt(A[i,i]-sum)
return L

The parameter A is a Paida Matrix instance for which +/- and * are matrix addition/subtraction and multiplication. If A is symmetric positive definite, the routine returns the unique lower triangular matrix L such that

A = L*L.transpose() # in Python/Paida
A = L##transpose(L) ; in IDL

With this routine it is easy to solve the generalized eigenvalue problem which is needed to calculate MAD transformation coefficients. Whether or not all this will run on the App Engine within my free CPU quotas remains to be seen.

Python the hard way? (Cont'd)

2010-03-18T21:07:00.004+01:00

OK, mucking about with Python on the Google App Engine is a bit of fun after all. I discovered, next to PyPng, another pure Python toolbox called PAIDA, which has a matrix object, can solve eigenvalue problems and implements BLAS (Basic Linear Algebra Subprograms). Also I fugured out how to upload and store not just ENVI images, but JPEG, PNG, GIF, TIFF, BMP and ICO in ENVI standard format, see the App.

I discovered in the process that there are some very nasty differences between the App Engine SDK and the App Engine proper (i.e. Google's Cloud).

For example, there is a memory limit on the App Engine, but none on the SDK (others have remarked with dismay on this). So you get a nice application running on your own machine with the SDK, deploy it and ... April Fool!

Here's a subtle one: when you convert, say, a JPEG to a PNG with App Engine's Image class, it hangs an alpha channel onto the result. When you do it with the SDK's Image class, it doesn't.

Next thing on the agenda: principal components analysis. Then, if I can use PAIDA to write a pure Python routine for Cholesky decomposition, I'll try to implement the MAD transformation.

Python the hard way? (Cont'd)

2010-03-15T20:50:00.004+01:00

Maybe I should start a new Blog: Fun With Python/Google App Engine. Except ... it's not fun. I thought I could at least play around with the PIL on my App. But the App Engine clips the PIL's wings so severely that about all you can do is make thumbnails. Cool.

Then there's the datastore. It supplies the only persistence, since you're in a sandbox and can't write to the App server's file system. It allows the user to store and retrieve a maximum of 1 MB at a time. Cooler.

Only pure Python modules can be used. So no Numpy, no Scipy, no Matplotlib. Coolest.

Fortunately, PyPng is written in pure Python, so I included it in my deployed app. I can take an uploaded ENVI image from the datastore and, after doing something trivial like calculating an NDVI, serve the result as an image. Can't for the life of me figure out how to insert the PNG image into an HTML document, though. But it must be possible.

I guess scientific computing on Google's cloud has a long way to go.

Python the hard way?

2010-03-10T18:03:00.007+01:00

Most introductory books I've seen on Python emphasize how easy it is to learn, and then attempt to illustrate that fact in the next 400 to 1000 pages. OK, Python is elegant and great to work with. But my experience is that, when learning any new computer language, it is always best to have a concrete, nontrivial task in mind and then just dive into it.

So I'm diving into the Python API for the Google App Engine (see an earlier Blog), never having before programmed a single line of Python in my life. Here you can watch me mucking about.

Kernel PCA

2010-02-06T13:14:00.012+01:00

In a paper on kernel principal component analysis (KPCA) (Kim et al., Transactions on Pattern Analysis and Machine Intelligence 27(9) 2005, 1351-1366) it is stated:

"In contrast to linear PCA, KPCA is capable of capturing part of the higher-order statistics which are particularly important for encoding image structure."

I offer a (CUDA-enabled) script for doing kernel PCA on my software page. Is it useful for remote sensing image analysis? Here is band 3 (near infrared) of an ASTER image, as usual over my home town Juelich:

I ran a kernel principal component analysis on all three visual/near infrared bands. Here is the fourth kernel principal component:

If you click on these two images to enlarge them, you will see that KPCA can act as a pretty good edge detector, something that ordinary PCA certainly can't.

Learning CUDA

2010-01-14T21:35:00.017+01:00

I've been experimenting a lot with GPULib (see many earlier blogs), a library which shields you nicely from the details of CUDA programming. But to learn more about what goes on behind the scenes, I decided to try my hand at writing my own CUDA DLM for IDL. The function is pretty pointless: it does image contrast stretching on the GPU with a pre-calculated lookup table, something that plain IDL does very efficiently anyway. But it was a start.

For example it was interesting to time what's going on. The function, called CUDA_STRETCH(LUT, IMG), takes a 256-byte lookup table LUT and an image band IMG as input and returns the stretched result. So both arrays have to be passed to and from the GPU each time. Here is the calling program:

pro STRETCH
envi_select, title='Choose MS image band', $
fid=fid, dims=dims, pos=pos, /band_only
if (fid eq -1) then return
img = bytscl(envi_get_data(fid=fid,dims=dims,pos=pos))
h = histogram(img,min=0,max=255,nbins=256)
; make a histogram equalization lookup table
LUT = byte(255*(total(h, /cumulative)/total(h)))
; call the CUDA DLM once to get it loaded
imgout = cuda_stretch(LUT, img)
; now time it
start = systime(2)
imgout = cuda_stretch(LUT, img)
print, 'CUDA: ',systime(2)-start
; return result to ENVI
envi_enter_data, imgout
; now time the IDL version
start = systime(2)
imgout = LUT[img]
print, 'IDL: ',systime(2)-start
end

The console output for a 4000x4000 image is

CUDA: 0.15700006
IDL: 0.015000105

(I said it was pointless.) My first thought was that most of the overhead was due to host - device transfers, but here is what the internal timer in the DLM says:

cuda_stretch ---------------------
clocks per second 1000
transferring arrays to device ...
cols = 4000 rows = 4000
clocks required 15
launching kernel ...
kernel done ...
clocks required 110
transferring result to host ...
clocks required 15

So the kernel is the problem, and I guess it is just programmed too naively.

Anyone else who would like to experiment along these lines can get the C code, DLM description file, Visual Studio 2008 project file, etc., etc here.

The Madman

2009-12-06T15:27:00.012+01:00

My book documents the ENVI/IDL extensions for MAD change detection and radiometric normalization pretty thoroughly, but after answering several emails with the same or similar questions from users, I thought either a short manual or a Frequently Asked Questions compilation might be in order. I opted for the manual because, to my childlike (childish?) brain, MAD MAN sounded cooler than MAD FAQ.

I hope this helps, as they say on the newsgroups.

GPULib 1.2.2

2009-11-18T18:37:00.007+01:00

Just downloaded and installed the latest version of GPULib from the Tech-X website.

I noticed in their announcement that they have "added several missing Windows build files." This caught my interest, as I haven't yet been able to get a build to work on my Windows machine.

So I loaded the Visual C++ solution file provided in the distribution into Visual C++ 2008 Express Edition and let the Wizard convert it from 2005 to 2008 format.

The first build attempt complained about not finding idl_export.h and idl.lib. After copying these from my IDL installation to the gpulib-1.2.2\idl directory, Lo and Behold -- it built!

However, I got 466 warning messages, so I didn't really believe it until I replaced the distributed DLL with the newly built one. All of my IDL/GPULib programs run as before. Very nice.

Cloudy thoughts

2009-11-09T19:00:00.012+01:00

Browsing (physically) through a book store a few days ago I came across Charles Severance's Using Google App Engine. I took it home and devoured it in one evening. It's absolutely fascinating what Google allows you to do on the web for free. I downloaded the App Engine simulator to play around with and started dreaming of putting MAD/RADCAL on the "cloud", programmed in Python. Of course number crunching is not exactly what Google App Engine is for, but it would be a great platform for demonstration. Small images wouldn't use much CPU time.

Anyway, I sobered up pretty quickly when I learned that NumPy (the Python extension that allows IDL/Matlab-style programming) isn't supported. A lot of users are asking for NumPy support in the Google forums but, if I understand correctly, its reliance on C-libraries more or less disqualifies it from deployment on the Google infrastucture (the "cloud"). But maybe I'm wrong?

Radiometric normalization revisited

2009-11-02T18:51:00.028+01:00

Judging from the emails I've received lately, quite a few people are using the ENVI/IDL extension RADCAL.PRO for radiometric normalization. In an earlier blog I gave an example, taken from Canty and Nielsen (2008), where the method gave a nice result. Now I'll show a situation (from the same paper) where it fails rather miserably, and then indicate how it might be made to work after all.

The example again involves two LANDSAT 7 ETM+ images acquired over Juelich, Germany, this time in May and late August of the same year. The amount of vegetation change in the scene is extreme, so that when IR-MAD is run on the six non-thermal LANDSAT bands it does not converge to a well-defined background of no-change pixels:

The above is an RGB composite of MAD bands 4,5 and 6. It is very noisy, and a clean no-change background is difficult to make out. Using the associated Chi-Square image with RADCAL to perform automatic normalization doesn't lead to a reasonable result. Here for example is the regression line for band 5:

The test statistics from RADCAL confirm the failure. The p-values for the F-test for equal variances of the reference and normalized target images are zero or very small for all six bands (Blogger is truncating the table at three bands, but you get the idea):

RADCAL statistics: 4245 training and 2123 test pixels
Variances
target 165.2327 297.3767 895.3079 401.0657 991.1505 1274.5616
reference 76.3946 138.7130 381.9422 347.6269 410.2062 399.8563
normalized 28.5304 43.2249 16.4233 301.8247 56.1722 1.2324
F-stat 2.6777 3.2091 23.2562 1.1518 7.3027 324.4514
p-value 0.0000 0.0000 0.0000 0.0011 0.0000 0.0000

But now we'll choose a spatial subset of the images, centered on one of the mining areas, where we might expect to find a higher fraction of invariant pixels and therefore hope for better convergence. The IR-MAD image certainly looks much cleaner:

We're also encouraged by the fact that the no-change pixels (middle gray) correspond both to forest canopy and to built-up areas. Running RADCAL again, this time on the chosen subset, we get good regression lines. For instance for Band 5:

and the statistical test also tells us that we can't reject the hypothesis of equal variances for the reference and normalized images:

RADCAL statistics: 453 training and 227 test pixels
Variances
target 45.7903 69.3390 182.2367 435.2744 487.0745 352.9176
reference 32.0876 49.1145 122.2010 139.4574 270.9004 197.6096
normalized 34.9058 53.8019 131.8851 137.8543 275.6107 203.5141
F-stat 1.0878 1.0954 1.0792 1.0116 1.0174 1.0299
p-value 0.5274 0.4938 0.5670 0.9308 0.8970 0.8251

So now we can tell RADCAL to use the regression coefficients obtained from the spatial subset to calibrate the entire scene, and we may have a useful normalization after all.

Shifting the mean

2009-10-08T21:08:00.015+02:00

One of the subjects that I take up in the second edition of my book is image segmentation. A very popular segmetation algorithm is the so-called mean shift and it works like this: Suppose we want to segment an image like the one below, which consists of the first 3 principal components of a 9-band ASTER image:

Take all of the 3-dimensional pixel vectors and extend them to 5-dimensional vectors by adding their x and y spatial coordinates, normalized to the 3 spectral intensity coordinates in some convenient way. Choose any pixel and, in the 5-dimensional feature space, center a hypersphere of radius r on that pixel. Calculate the (5-dimensional) mean of all of the pixels in the hypersphere and then shift its center to that mean. Then recalculate the mean and shift again. Keep doing this until convergence. The point of convergence is called a "mode", a local maximum in the density distribution of the pixels in the feature space. Replace the spectral components of the original pixel with those of the mode but retain the pixel's original spatial coordinates. Now do the same for all of the other pixels in the image. The result is a segmented image:

There are about 740 modes or segments in the above picture, the spatial coordinates of the modes being roughly at the segment centroids. The positions of the modes in the spectral subspace look like this:

Here I've clustered them for fun, using 5 classes. If we re-lable the segments with the corresponding class color, we get - in a roundabout way - a rather nice unsupervised classification:

Here's a blow-up of a spatial subset of the above with the original segment boundaries overlayed:

Clustering changes

2009-09-28T10:55:00.012+02:00

I've discussed the generation of multi-band change images with the iterated MAD method and, more recently, the clustering of multispectral images with a Gaussian mixture model. So why not combine the two? Here's an example:

The image below is an RGB composite of three of six MAD components showing changes in the region around Juelich, Germany, where I live. The MAD components were derived from two Landsat 7 ETM scenes acquired in the months of June and August, 2001.

The blue-gray (featureless) regions correspond to built-up areas and to forest canopy, neither of which have changed significantly between acquisitions. The significant changes, indicated by the various colored areas, are quite heterogeneous and not easy to interpret, but from their form they are clearly associated with cultivated fields. The main crops around here are sugar beets, corn and cereal grain. The former two are still maturing in August, whereas the grain fields have been harvested. So we might postulate two main classes of change: maturing crops and harvested crops. Allowing for a no-change class and a "catch-all" class for other kinds of change, we can try to classify the MAD image with the Guassian mixture clustering algorithm by assuming four classes in all. The resulting clusters, in the feature space of the three MAD components, are shown below with to help of ENVI's n-d visualiser.

There are three well-defined clusters (the central (white) one is no change) and a diffuse cluster (blue). Identifying the red and green classes as the two postulated ones and overlaying them onto band 4 of the August Landsat image we get:

The red regions correspond to low reflectance in band 4 (near infrared) and therefore (presumably) to harvested grain, the green regions (presumably) to maturing crops.

Confessions of a "no-C" programmer

2009-09-23T14:58:00.005+02:00

I wrote my first computer program way back in 1961 in a weird machine language on an even weirder machine called, as I recall, the Bendix G15D. From thence I graduated to Fortran on an IBM 1620 and, over many subsequent years, made sacrificial offerings of boxes upon boxes of punched cards to ever larger and hungrier IBM mainframes.

But then a wonderful/horrible invention called the IBM PC suddenly appeared on the scene. So I had to learn to use BASIC for a while, followed by a love affair with Turbo Pascal which continued up until about Delphi 3.0. After two illicit affairs with Prolog and Mathematica I finally landed in the here and now with good old IDL.

Something essential is missing, right? Yes, over all these years I have actually succeeded not to learn how to use C (or its derivatives C++ and Java). So why am I confessing this shameful fact? Because, at 66, I'm living proof that "you're never too old to learn". I proudly announce that I've just succeeded in writing my very first DLM for IDL, using IDL_VPTR's, stuff like "X->value.arr->data", the whole, cryptic lot!

The DLM replaces the old CALL_EXTERNAL() DLL for the provisional means algorithm which is used in the MAD program to collect image statistics. I haven't put it on my website along with a version of MAD which uses it yet, because experience shows that that can cause no end of confusion (see earlier blog). I'll think of something.

Clustering Images on the GPU

2009-09-08T16:37:00.024+02:00

Well, I seem at last to have a version of Gaussian mixture clustering working with GPULib. Unfortunately, it crashes my system when I run it on large images (> around 1000x1000 pixels) because of limitations on the GPU memory. It gives a speedup factor of about four relative to the original CPU version described in Chapter 7 of my book (same as for the neural net I talked dabout a few blogs ago). Among some other cool features, the algorithm uses multiresolution image pyramids created with the discrete wavelet transformation to speed up clustering.

Click on the image above to see successive stages for a 3-level pyramid created from a Landsat 7 ETM image. The red cluster is mixed forest. It becomes more sharply defined from left (120m resolution) to right (30m resolution). Clustering only at the native 30m resolution (no pyramid) takes five times longer.

You can get the routine here.

Are you MAD because CHOLDC failed?

2009-08-23T15:45:00.004+02:00

Been getting several emails from users of my ENVI/IDL extension for MAD, saying that they are getting the message "CHOLDC failed". The reason is that I updated the external routine PROV_MEANS.DLL some while back. The current version of the MAD program also requires the current version of the DLL. If the old DLL is anywhere on the OS path, it still may be picked up, so be sure to eradicate it completely when installing the new version. Humble apologies to all affected.

GPULib 1.2

2009-07-16T14:43:00.006+02:00

I downloaded and installed the latest version of GPULib on my windows machine. It gave me wierd results at first, till I realized that I had to point the DLM path to

\gpulib-1.2\bindist\idl\x86\single

as my antiquated nvidia card doesn't support double precision.

Version 1.2 caters mostly to Matlab, see the Tech-X blog. But some bugs in IDL have also been fixed, in particular the dimension parameter in GPUTOTAL() which was giving me headaches.

For anyone else who might be developing ENVI extensions with GPULib, I find the following little program to be useful when placed in SAVE_ADD:

PRO aaa_run_define_buttons, buttonInfo
ENVI_DEFINE_MENU_BUTTON, buttonInfo, $
VALUE = 'GPULib API', $
REF_VALUE = 'About ENVI', $
EVENT_PRO = 'aaa_run', $
UVALUE = 'AAA',$
POSITION = 'after'
gpuinit
END

PRO aaa_run,event
mg_open_url, $
'http://www.txcorp.com/products/GPULib/'+ $
'idl_docs/index.html'
END

ENVI picks it up first (I think it scans the SAVE_ADD directory alphabetically) and runs the define_buttons procedure to place a link to the GPULib API in the help menu. (You can get MG_OPEN_URL from Mike Galloy.) As a side effect, GPUINIT is called so that any other GPULib code in SAVE_ADD will compile properly. Maybe it only works on windows.

Normalizing Images

2009-07-07T20:43:00.023+02:00

I said in my previous blog that the MAD components are insensitive to uninteresting effects like solar illumination, differing sensor gains or atmospheric absorption. In fact, if the "uninteresting effects" are linear in nature, then the MAD components are strictly invariant.

Mathematically, let X be any pixel vector in an ASTER VNIR image at time T1 and Y any pixel vector in a second image taken at a later time T2. These vectors each have three components corresponding to the three spectral bands. Suppose that A and B are two arbitrary 3x3 matrices and that C and D are two equally arbitrary 3-component vectors. Now look at the linear transformations X' = AX+C, and Y' = BY+D. Performing the MAD procedure on X' and Y' gives exactly the same result as performing it on X and Y!

What this means is that the MAD algorithm will find the same change pixels and, by implication, the same no-change pixels irrespective of any linear perturbations, attenuations or shifts of the image data.

Suppose one image is acquired under "favorable" conditions, e.g., high solar elevation, clear visibility, optimal sensor gain, and the second image under poorer conditions. Then it would be nice to radiometrically normalize the second image to the first one. Then one could, for example, compare vegetation indices in the two images in a meaningful way. To do the normalization properly, we need to base it on the no-change pixels. Including pixels which correspond to real, phyisical changes on the ground would invalidate the normalization. The MAD algorithm will provide us with the no-change pixels.

An example: The two ASTER VNIR images below were acquired over an area east of Esfahan, Iran in April 2001 (left) and September 2005 (right). The pixel intensities are displayed as digital numbers in the range 0-255 and are proportional to the radiance measured at the satellite sensor.

There's alot of cloud cover in the left hand image and big differences in vegetation cover (indicated by the red pixels). We'll normalize the second image to the first image. After running MAD to obtain the pixels with highest probability of no change (greater than 90% in this example), a band-by-band orthogonal linear regression of the first (reference) image against the second (target) image is performed:

The regression coefficients are then used to scale the target image intensities to those of the reference image:

Click on the above picture to enlarge. On the right hand image the pixels that took part in the normalization are visible as white dots. There's only a few hundred of them, but they suffice. The best part is, the whole procedure is automatic.

Detecting Changes

2009-06-16T12:29:00.000+02:00

I'll leave the topic of GPU processing in IDL for a while and blog a little on another of my favorite subjects, change detection with multispectral imagery. Allan Nielsen and I have collaborated on several papers dealing with his MAD (multivariate alteration detection) algorithm. Here is an unpublished example of MAD in action:

The disasterous earthquake of October 8, 2005 had a magnitude of 7.6 and was centered near the city of Muzaffarabad in Kashmir. The figure below is an ASTER image over the area acquired on Sept. 7, 2005. The image size is 1000 by 1000 pixels. Ground resolution is 15 meters and the image consists of three bands in the VNIR (visual and near infra-red) spectrum. Since the terrain is rough, the image was orthorectified using ASTER's stereo capability.

Here is the same scene acquired with the ASTER satellite on October 27, 2005, similarly orthorectified and carefully registered to the first image:

We could now look for changes caused by the earthquake just by subtracting the two images from eachother, band-for-band, and this is in fact what many people do. The MAD algorithm is more sophisticated. It generates three new bands (for each of the two images) which are linear combinations of the original bands. Whereas the original bands were sorted by spectral wavelength, the new bands of the first image are sorted according to their maximum similarity or correlation with the corresponding new bands of the second image. This is called canonical correlation and is a very old and well-known statistical device. The new bands of the second image are then subtracted from those of the first image to give the MAD components (in this case three of them) which contain the change information.

Forcing the two images in this way to be as similar as possible before subtracting them makes sense. Many differences between the original images will be due to uninteresting effects like solar illumination, differing sensor gains or atmospheric absorption and the canonical correlation procedure will render the MAD components insensitive to such effects.

However what we really wish to do is to force the invariant pixels to be as similar as possible, not the pixels which signal true changes on the ground. This is where the MAD algorithm really excels, because it can be iterated. Having obtained the MAD components, we use them to make an initial estimate of the no-change probability of each pixel in the image. This is possible because of the nice statistical properties of the MAD components. Then we repeat the canonical correlation, this time weighting each pixel with its probability of no change, and recalculate the no-change probabilities. The procedure continues until the similarity measures (the correlations) cease to change. The figure below shows the values of the correlations for each successive iteration. They increase steadily as the invariant pixels become better and better discriminated.

Have a look at the first MAD component projected onto Google Earth using Mike Galloy's IDL program MG_WRITE_KML.PRO. Bright and dark pixels signify change, middle gray no change. Use the slider to vary the overlay's transparency. Tilt the view and get in close to see that, at this resolution, it is the landslides that are most evident.

Another interesting use of MAD is for automatic radiometric normalization of multitemporal images. Stay tuned.

A neural network on the GPU - Concluded

2009-05-31T21:06:00.000+02:00

OK, it's "finished". Here's a summary of the features:

- Two-layer, feed-forward neural network for supervised
land cover/land use classification of multispectral imagery,
fully integrated into ENVI.

- Makes use of the cross entropy cost function and softmax
outputs to model posterior class membership probabilities.

- Trained with the scaled conjugate gradient algorithm,
guaranteed to converge monotonically to a (local) minimum
in the cost function, see the Figure. The algorithm is
described in detail in Appendix B of my book.

- Apart from the number of neurons in the first (hidden) layer,
there are no adjustable parameters whatsoever.

- Uses the R-operator method (Bishop 1995, Neural
Networks for Pattern Recognition) to evaluate the Hessian
(matrix of second derivatives of the cost function) efficiently.

- All matrix operations are carried out on the GPU using
GPULib IDL bindings to CUDA. CUDA is giving me a speedup
factor of around 3-4 on my system (see earlier blogs)
relative to the CPU version.

- One heck of a lot faster than ENVI's built-in network
classifier. Also much faster that ENVI's support vector machine,
with comparable accuracy.

- You can get the code here.

A neural network on the GPU

2009-04-29T20:36:00.000+02:00

In my spare time I'm trying to put a neural
network classifier onto the GPU with the help
of GPULib. It's a lot of work because the training
algorithm, scaled conjugate gradient, is fairly
complicated. I'm converting the original CPU version
that I wrote for my book piece-by-piece, so as
not to get completely frustrated.

Here is an chunk of code (which is now working)
which propagates an input vector through the network
to give an output signal that is supposed to reflect
its class membership. (The network is programmed as
an IDL object class, hence all the SELF-references.)

Pro GPUFFNCG::forwardPass
gpuView,*self.N_gpu,self.np,self.LL*self.np,Nv_gpu
gpuReform,Nv_gpu,self.np,self.LL
;logistic output of hidden layer
expnt_gpu = gpuMatrix_Multiply(*self.Gs_gpu,$
*self.Wh_gpu,/btranspose)
tmp_gpu = gpuExp(1.0,-1.0,expnt_gpu,0.0,1.0)
gpuPutArr,fltarr(self.np,self.LL)+1.0,onesL_gpu
gpuDiv,onesL_gpu,tmp_gpu,Nv_gpu
;softmax network output
Io_gpu = gpuMatrix_Multiply(*self.N_gpu,$
*self.Wo_gpu,/btranspose)
maxIo_gpu = gpuMax(Io_gpu,dimension=2)
for k=0,self.KK-1 do begin
gpuView,Io_gpu,k*self.np,self.np,Iov_gpu
gpuSub,Iov_gpu,maxIo_gpu,Iov_gpu
endfor
A_gpu = gpuExp(Io_gpu)
sum_gpu = gpuTotal(A_gpu,2)
for k=0,self.KK-1 do begin
gpuView,*self.M_gpu,k*self.np,self.np,Mv_gpu
gpuView,A_gpu,k*self.np,self.np,Av_gpu
gpuDiv,Av_gpu,sum_gpu,Mv_gpu
endfor
;cleanup
gpuFree, [tmp_gpu,expnt_gpu,onesL_gpu,Io_gpu,$
maxIo_gpu,A_gpu,sum_gpu]
End

I try to make as much use of GPUVIEW as possible
in order to avoid shifting arrays around, and of
course I try to avoid transfers back to the CPU.

The hardest part is calculating the matrix of second
derivatives of the cost function with respect to the
synaptic weights (the Hessian). This is where I hope
for a big speedup using GPULib. Still working on that.

More to come, I hope.

Home page update

2009-02-24T21:03:00.000+01:00

Just re-vamped my home page a bit, making terribly efficient use of all 10 (count them, ten) MB of webspace that Deutsche Telekom generously provides me free of charge with my call and surf "comfort package", as they humorously call it.

I've put up the most recent versions of my ENVI extensions, some of which actually work. The kernel PCA and kernel k-means routines, both of which use GPULib, seem stable enough. I'm still stuggling with a CUDA version of the Gaussian mixture classifier, unsure as usual if there's a bug in GPULib or if I'm to blame.

All of the programs and their subroutines have been churned through IDLDOC, but the list of warnings about partial or missing documentation is still very long. I want to fill in the gaps gradually as I have time.

Learning IDLDOC

2009-01-19T18:06:00.000+01:00

I just installed Mike Galloy's IDLDOC plugin for the IDL Workbench and I'm learning how to use it. My first feeble attempts can be scoffed at on my software page, where I've started to document the latest versions of the ENVI extensions in my book. By the time the Revised Second Edition comes out later this year, I hope to have mastered the techniques. Already, though, I think it's a fantastic tool.

One problem I have is that the ".\" prefixed to the filenames in ALL-FILES.HTML generated by IDLDOC is confusing my (much-loved) Firefox browser. IE6 has no difficulty. Just now I'm removing them by hand.

Using GPULib 1.0

2008-12-29T11:24:00.000+01:00

Now that I've got GPULib 1.0 up and running, I'm making use
of the new functional forms of the library routines to write
kernel-based algorithms. Given a set of m n-dimensional
observation vectors G_i

G_i = ( g_1, g_2 ... g_n )_i, i=1,2,...m,

in the form of a m x n data matrix:

G = ( G_1, G_2, ... G_m)^T

where ^T denotes transposition, I want to calculate the symmetric
m x m Gaussian kernel matrix

K_ij = exp(-gma ||G_i-G_j||^2 ), i,j = 1,2,...m.

In ordinary IDL:

function kernel_matrix, G, gma=gma
m = n_elements(G[0,*])
K = transpose(total(G^2,1))##(intarr(m)+1)
K = K + transpose(K) - 2*G##transpose(G)
return, exp(-gma*K)
end

WIth GPULib:

function gpukernel_matrix, G_gpu, gma=gma
m = G_gpu.dimensions[1]
n = G_gpu.dimensions[0]
onesm_gpu = gpuputarr(fltarr(m)+1)
onesn_gpu = gpuputarr(fltarr(1,n)+1)
G2_gpu = gpumult(G_gpu,G_gpu)
tmp1_gpu = gpumatrix_multiply(onesn_gpu,G2_gpu)
gpufree,G2_gpu
K_gpu = gpumatrix_multiply(onesm_gpu,tmp1_gpu)
tmp2_gpu = gpumatrix_multiply(tmp1_gpu, $
onesm_gpu,/atranspose,/btranspose)
gpufree,[onesm_gpu,onesn_gpu,tmp1_gpu]
K_gpu = gpuadd(K_gpu,tmp2_gpu,lhs=K_gpu)
gpufree,tmp2_gpu
tmp3_gpu = gpumatrix_multiply(G_gpu,G_gpu, $
/atranspose)
K_gpu = gpuadd(1.0,K_gpu,-2.0,tmp3_gpu,0.0,$
lhs=K_gpu)
gpufree,tmp3_gpu
K_gpu = gpuexp(1.0,-gma,K_gpu,0.0,0.0,$
lhs=K_gpu)
return, K_gpu
end

Since this routine is called many times in the
processing of a multispectral image, a memory
leak would be disasterous. Note the careful use
of the LHS keyword and of the GPUFREE procedure
to avoid this.

The code is now more readable than for the
procedural forms, but would be more readable
still if there were a GPUTRANSPOSE() function and
if GPUTOTAL() had a dimension parameter.

More on GPULib and Kernel Methods

2008-12-12T17:16:00.000+01:00

My last post was way back in August and I'm writing this one really just to keep the blog alive, so to speak. So far I haven't been able to get GPULib 1.0 going on my system, so I have no idea if the programs on my website will still run under 1.0. I'm getting help from Tech-X and haven't given up hope.

Anyway, I'm very keen to apply CUDA (via GPULib) to the kind of nonlinear algorithms (kernel methods) I've described previously. These algorithms involve two steps: training and generalization. In remote sensing applications, the training data set can be kept fairly small (order 10^3), but generalization invariably means processing all the pixel vectors in an image, typically of the order 10^6-10^8. Since generalization can be safely done in single precision and involves manipulation of large matrices, it is a natural task for the GPU.

GPULib: Kernel K-means

2008-08-28T18:58:00.000+02:00

Just like principal components analysis, the subject of previous posts, the well-known K-means clustering algorithm can be kernelized. Shawe-Taylor and Cristianini, in their excellent book Kernel Methods for Pattern Analysis, give Matlab code. Based on that I've written an ENVI extension for kernel K-means using the GPULib library. Here are some times for the multispectral image clustering shown above:

Without GPULib
-------------------
Kernel K-means
Thu Aug 28 19:44:12 2008
-------------------
input file juelich_august_29_pca
Number of Classes 8
elapsed time: 32.500000
Result written to memory

With GPULib

-------------------
Kernel K-means
Thu Aug 28 19:48:16 2008
-------------------
input file juelich_august_29_pca
Number of Classes 8
elapsed time: 4.0160000
Result written to memory

Speedup factor: 8 .

Kernelized PCA with GPULib

2008-08-10T12:35:00.000+02:00

I haven't had much time lately to play with GPULib, although there are plenty of ENVI/IDL apps that I want to experiment with. An ENVI extension for kernelized PCA (nonlinear principal components analysis using a radial basis kernel, mentioned in an earlier blog) can be downloaded from my software website. It is pretty rudimentary, but illustrates the use of quite a few GPULib functions. See also this short ppt on nonlinear PCA with GPULib.

GPULIB (Clustering speedup)

2008-07-19T13:35:00.000+02:00

The clustering program that I've been re-working with GPULib fits a mixture of K multivariate Gaussians to a cloud of n N-dimensional image pixel vectors using the Expectation Maximization algorithm:

1. Initialize n x K cluster membership
probability matrix U

2. Maximization step: calculate
K N-component mean vectors,
K NxN covariance matrices,
K mixture coefficients

3. Expectation step: recalcuate
U and normalize.

4. If U hasn't changed significantly,
stop, else go to 2.

Step 3 requires inversion of the cluster covariance matrices, which has to be done on the CPU. Apart from that, I still get into trouble every time I use gpuWhere (don't know why yet). So, since I have to check U for zeroes, I have to do its normalization on the CPU as well. Here is my attempt at normalization on the GPU, which usually crashes IDL:

; normalize U
gpuMatrix_Multiply,U_gpu,onesKxK_gpu,den_gpu
gpuEq,den_gpu,fltarr(n,K),zeroes_gpu
gpuWhere,zeroes_gpu,ind_gpu,count
if count gt 0 then $
gpuSubscript,ind_gpu,den_gpu,$
fltarr(count)+1,/LHS
gpuDiv,U_gpu,den_gpu,U_gpu
gpuFree,[den_gpu,ind_gpu,zeroes_gpu]

Speedup (so far) on a 1000 x 1000 x 4 image clustered into 8 classes (i.e. n=1,000,000, N=4, K=8) is a factor of 2.8.

GPULIB (trip to never-neverland)

2008-07-10T16:51:00.000+02:00

I seem to be making a little progress with the clustering algorithm. One problem is that I have to invert a covariance matrix on each iteration, which means jumping back to the CPU to do it. This probably slows things down. Doesn't CUBLAS provide matrix inversion?

Anyway, while coding I stumbled into a horrible trap:

gpuPutArr,A,A_gpu
gpuView,A_gpu,n,m,Av_gpu
;...
;lots of code, so I kind of forgot
;what Av_gpu is
;...
gpuFree,A_gpu
gpuFree,Av_gpu ; very, very bad thing to do

This put me Somewhere Over the Rainbow. Nothing, but nothing works anymore on the GPU until I exit IDL and restart. Took a long time to get back to earth. Incidently, another way to get up there with Judy Garland is to forget that

gpuWhere,A_gpu,ind_gpu,count

overwrites A_gpu.

GPULIB succes story (cont'd)

2008-07-08T16:42:00.000+02:00

Regarding my last post, Peter Messmer at Tech-X pointed out to me that I should use gpuView rather than shifting subarrays around with gpuSubArr. Doesn't make things much faster in this case, but the code is nicer:

while i lt num_rows do begin
; ggi = GG[*,i*num_cols:(i+1)*num_cols-1]
gpuView,GG_gpu,i*num_cols*num_bands,$
num_bands*num_cols,ggi_gpu
gpuReform,ggi_gpu,num_bands,num_cols
; KK = transpose(G2)##(intarr(num_cols)+1)
gpuMatrix_Multiply,onesnc_gpu,G2_gpu,$
KK_gpu,/btranspose
; KK = KK + transpose(intarr(m)+1)##$
; GG2[i*num_cols:(i+1)*num_cols-1]
gpuView,GG2_gpu,i*num_cols,num_cols,gg2i_gpu
gpuMatrix_Multiply,gg2i_gpu,$
onesm_gpu,res_gpu,/btranspose
gpuAdd,KK_gpu,res_gpu,KK_gpu
; KK = KK - 2*G##transpose(ggi)
gpuMatrix_Multiply,ggi_gpu,G_gpu,$
res_gpu,/atranspose
gpuAdd,1.0,KK_gpu,-2.0,res_gpu,0.0,KK_gpu
; KK = exp(-gma*KK)
gpuExp,1.0,-gma,KK_gpu,0.0,0.0,KK_gpu
; image[*,i,*] = alpha[*,0:npc-1]##KK
gpuMatrix_Multiply,KK_gpu,alpha_gpu,res1_gpu
gpuGetArr,res1_gpu,res
image[*,i,*]=res
i++ ; next row
endwhile

With my new-found skills ;-) I want next to have a go at image clustering. If you cluster all of the pixel vectors in a multispectral image (typically millions) then things can get pretty slow.

GPULIB success story (at last)

2008-07-01T15:56:00.000+02:00

In kernelized PCA (nonlinear principal components analysis) you have to diagonalize a large (say 1000x1000) Gram matrix (real, symmetric kernel matrix) and then project the original image along the most significant eigenvectors. There is no way to do diagonalizations on the GPU (yet?), but the projection part is also very time consuming. Here is my solution using GPULIB. The commented statements are the original code on the CPU:

start_time=systime(2)
while i lt num_rows do begin
; ggi = GG[*,i*num_cols:(i+1)*num_cols-1]
gpuSubArr,GG_gpu,[0L,-1L],$
[i*num_cols,(i+1)*num_cols-1],$
ggi_gpu,[0L,-1L],[0L,-1L]
; KK = transpose(G2)##(intarr(num_cols)+1)
gpuMatrix_Multiply,onesnc_gpu,$
G2_gpu,KK_gpu,/btranspose
; KK = KK + transpose(intarr(m)+1)##$
; GG2[i*num_cols:(i+1)*num_cols-1]
gpuSubArr,GG2_gpu,[i*num_cols,(i+1)*num_cols-1],$
gg2i_gpu,[0L,-1L]
gpuMatrix_Multiply,gg2i_gpu,onesm_gpu,res_gpu,$
/btranspose
gpuAdd,KK_gpu,res_gpu,KK_gpu
; KK = KK - 2*G##transpose(ggi)
gpuMatrix_Multiply,ggi_gpu,G_gpu,res_gpu,$
/atranspose
gpuAdd,1.0,KK_gpu,-2.0,res_gpu,0.0,KK_gpu
; KK = exp(-gma*KK)
gpuExp,1.0,-gma,KK_gpu,0.0,0.0,KK_gpu
; image[*,i,*] = alpha[*,0:npc-1]##KK
gpuMatrix_Multiply,KK_gpu,alpha_gpu,res1_gpu
gpuGetArr,res1_gpu,res
image[*,i,*]=res
i++
endwhile
print,'elapsed time: ',systime(2)-start_time

Here is the time taken to project a 1000x1000x6 image on the GPU

elapsed time: 12.093000

and on the CPU

elapsed time: 197.39100

for a speedup of 16. Not bad.

GPULIB Where-Woes

2008-06-30T15:57:00.000+02:00

No one can live without the where() function. Is it faster on the GPU? Here is my attempt to time it:

pro test_gpuWhere
gpuinit
A = float(round(randomu(s,1000000)))
st = systime(2)
for i=0,99 do begin
gpuPutArr, A, A_gpu
ind = where(A,count)
endfor
print,count
CPUtime = systime(2)-st
st = systime(2)
_ = temporary(count)
for i=0,99 do begin
gpuPutArr, A, A_gpu
gpuWhere,A_gpu,ind_gpu,count
endfor
print,count
GPUtime = systime(2)-st
print, 'Speedup: ',CPUtime/GPUtime
gpuFree,A_gpu
gpuFree,ind_gpu
end

Several things to notice here:

A is an array of "float flags", i.e. 1.0s and 0.0s. Everything, repeat, everything in GPULIB is float. Also the count parameter in gpuWhere returns the sum of A, not the number of nonzero elements.

Before reusing count in gpuWhere I have to undefine it, because it is returned as longint from where(). GPULIB will report an error otherwise.

In the loop, I have to use putArr,A,A_gpu each time, because gpuWhere overwrites its first argument.

So is all this worth the trouble?

% Compiled module: TEST_GPUWHERE.
500933
500933.
Speedup: 1.3521798

Evidently not.

GPULIB rude awakening

2008-06-30T15:32:00.000+02:00

Turns out I had set the GPU into some kind of undefined state by incorrectly using the gpuWhere procedure in the session in the last post. Here's the correct output from test_gpuMatrix_Multiply:

Length : 1000
CPU time: 0.047000170
GPU time: 0.030999899
Speedup : 1.5161395
Length : 10000
CPU time: 0.34399986
GPU time: 0.26500010
Speedup : 1.2981122
Length : 100000
CPU time: 3.0780001
GPU time: 2.6250000
Speedup : 1.1725715
Length : 1000000
CPU time: 30.734000
GPU time: 26.282000
Speedup : 1.1693935

It was way too good to be true. Ah well.

OK, what about simple array products?

pro test_gpuMult
; initialize
gpuinit
; array of 1000000 3-element observation vectors (rows)
A = randomu(s,3,1000000)
gpuPutArr,A,A_gpu
; calculate square of A on the CPU
start = systime(2)
for j=0L,99 do C = A*A
CPUtime = systime(2)-start
print, 'CPU time: ',CPUtime
; now on the GPU
start = systime(2)
for j=0L,99 do gpuMult,A_gpu,A_gpu,C_gpu
GPUtime = systime(2)-start
print,'GPU time: ', GPUtime
print,'Speedup : ', CPUtime/GPUtime
gpuFree,A_gpu
; check that the results are the same
gpuGetArr,C_gpu,C1
print,'Check:'
print, total(C-C1)
gpuFree,C_gpu
end

This gives a speedup of 10:

% Compiled module: TEST_GPUMULT.
CPU time: 2.1100001
GPU time: 0.20300007
Speedup : 10.394086
Check:
0.000000

GPULIB continued

2008-06-28T17:30:00.000+02:00

It is often necessary to calculate the covariance matrix of a multispectral image from the array of observation vectors (pixels). The program

pro test_gpuMatrix_Multiply
; initialize
gpuinit
for i=3,6 do begin
; array of 6-element observation vectors (rows)
A = randomu(s,6,10L^i)
gpuPutArr,A,A_gpu
; calculate covariance matrix C on the CPU
start = systime(2)
for j=0L,99 do C = $
Matrix_Multiply(A,A,/btranspose)
CPUtime = systime(2)-start
print, 'Length : ',10L^i
print, 'CPU time: ',CPUtime
; now on the GPU
start = systime(2)
for j=0L,999 do $
gpuMatrix_Multiply,A_gpu,A_gpu,C_gpu,/btranspose
GPUtime = systime(2)-start
print,'GPU time: ', GPUtime/10
print,'Speedup : ', 10*CPUtime/GPUtime
gpuFree,A_gpu
endfor
; check that the results are the same
gpuGetArr,C_gpu,C1
print,'Check:'
print, determ(C),determ(C1)
gpuFree,C_gpu
end

produces

% Compiled module: TEST_GPUMATRIX_MULTIPLY.
Length : 1000
CPU time: 0.030999899
GPU time: 0.0062000036
Speedup : 4.9999808
Length : 10000
CPU time: 0.25000000
GPU time: 0.0046999931
Speedup : 53.191567
Length : 100000
CPU time: 0.96900010
GPU time: 0.0046999931
Speedup : 206.17054
Length : 1000000
CPU time: 9.5620000
GPU time: 0.0032000065
Speedup : 2988.1190
Check:
6.33251e+030 6.33251e+030

So for instance the covariance matrix of a 1000x1000x6 multispectral image could be calculated 3000 times faster on the GPU!? My ENVI/IDL change detection extension iterates on the weighted covariance matrix of a bitemporal image. Sooo ...

GPULIB Trials (and Tribulations)

2008-06-26T19:16:00.000+02:00

With active and friendly help from Tech-X Corporation I've been experimenting with GPULIB, an IDL/Matlab/Python ... interface to NVIDIA's CUDA. I was (and still am) hoping to speed up some of the ENVI/IDL extensions for remote sensing image analysis described in my book.

So I've started blogging (at my age!) to document my experiences for anyone who wants to take a similar plunge into the world of "massively parallel computing" on his or her home PC.

Mine is:

Intel Pentium 4 650, 3400 MHz (17 x 200)
Intel Lakeport-G i945G
2048 MB (DDR2 SDRAM)
NVIDIA GeForce 8600 GT (512 MB)

Having gotten a Windows XP build from Tech-X (I evidently don't have the right C++ compiler to do my own builds), I did the following

-Copied GPULIB.DLM and GPULIB.DLL to C:\Programme\ITT\IDL70\bin\bin.x86 (my DLL/DLM Path)

-Placed GPUINIT.PRO in my IDL path

-Downloaded and installed the latest NVIDIA drivers for my graphics card

-Downloaded the CUDA Toolkit version 1.1 for Windows XP 32bit

-Copied the DLLs (including CUDART.DLL) to C:\Programme\ITT\IDL70\bin\bin.x86

-Started IDL 7.0 (with ENVI 4.5)

-Ran the benchmark program BENCH.PRO

Here's what it said:

% Compiled module: BENCH.
0.756607 2.33993 0.196372 0.516154 0.0442747 0.839950
0.756607 2.33993 0.196372 0.516154 0.0442747 0.839950
N iter = 50
CPU Time = 9.2029998
GPU Time = 0.10900021
Speedup = 84.431032 (emphasis added)

Wow! If that isn't exciting I don't know what is. This with a cheap, passively-cooled graphics card!

Of course, what BENCH.PRO actually does is calculate the log gamma function of a 1000000-element array 50 times, something I seldom have occasion to do in my ENVI extensions.

And that's where the story starts. Stay tuned.