Detecting Changes

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

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

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

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

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

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

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.