Re: Clsuter analysis wiht IDL

[Craig Markwardt <craigmnet(at)cow.physics.wisc.edu>]

dw@isva.dtu.dk (Dorthe Wildenschild) writes:

> I am trying to use the cluster analysis utility in IDL, and I can't work it 
> out. Perhaps due to my lack of knowledge in the statistical field.....
> 
> I have a 2D-image (658x658) with a single band of info (CT x-ray 
> intentities) that I want to classify (into three classes) using cluster 
> analysis. I assume I have to use the CLUST_WTS function first and then the 
> CLUSTER function, but I can't work it out.
> The on-line help isn't very helpful on this topic....
> 
> If I use
> weights = clust_wts(image, n_clusters=3), with image= intarr(658x658) I get 
> alot of floating errors.
> Do I need to reform the image to (658x658,1) before using the clust_wts 
> function? Doesn't seem to work either, though. Also how do I get the 
> cluster numbers back as an overlay of my image? so that I can actually see 
> the result of the classification.

I am by no means a statistics or clusters analysis expert, but from
looking at the documentation it seems that both CLUST_WTS and CLUSTER
require the n-dimensional *positions* of the data points, and not an
intensity map.  It is pretty clear that the cluster functions are
based on an unweighted set of scatter-points.  If you want to use your
intensity information as a weighting, you may be out of luck.

Perhaps you could achieve what you desire with this code, which simply
finds the non-zero pixels:

wh = where(image GT 0, ct)
if ct EQ 0 then message, 'ERROR: the image is blank!'
x = wh MOD 658      ; form x pixel positions
y = floor(wh / 658) ; form y pixel positions

xy = transpose([[x],[y]]) ; compute the 2-d scatter positions
weights = clust_wts(xy, n_clusters=3)
etc.

I haven't tried this, so it may take some tweaking.  Good luck,
Craig

-- 
--------------------------------------------------------------------------
Craig B. Markwardt, Ph.D.         EMAIL:    craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
--------------------------------------------------------------------------