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?
Monday, 9 November 2009
Monday, 2 November 2009
Radiometric normalization revisited
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.
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.
Thursday, 8 October 2009
Shifting the mean
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:
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:
Monday, 28 September 2009
Clustering changes
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.
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.
Wednesday, 23 September 2009
Confessions of a "no-C" programmer
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.
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.
Tuesday, 8 September 2009
Clustering Images on the GPU
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.
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.
Sunday, 23 August 2009
Are you MAD because CHOLDC failed?
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.
Subscribe to:
Posts (Atom)
