Re: Clsuter analysis wiht IDL

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

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

> --=====================_524776338==_.ALT
> Content-Type: text/plain; charset="us-ascii"; format=flowed
> 
> <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
> 
> I'm such a beginner at this I don't know what the MOD function does? (no 
> on-line help listing for it)
> When trying to transpose, IDL corrects me with
> Arrays are allowed 1 - 8 dimensions

MOD is documented as the "modulo" function.  You may have to scroll to
see the documentation entry in the index.  Since WHERE "pretends" that
your array is 1-d, you have to reconstruct the x and y positions by
dividing by the number of columns, and taking the remainder (=x) and
quotient (=y).

Here is an example with some simulated data.  The TRANSPOSE command
works as advertised.  The GAUSS2 simply generates a 2-d gaussian
function, and is available from my web page
http://cow.physics.wisc.edu/~craigm/idl/idl.html

This script finds all three clusters successfully, but it shows that
the routine is not optimal since it doesn't find the cluster centers.

Would you consider conferring with Mr. Rojas on his question about
clusters and kmeans?

Craig


;; Simulate some data.  Three clusters at (2,3), (-3,1), (1,-2)
x = findgen(100)*0.1 - 5. & y = x
xx = x # (y*0 + 1)
yy = (x*0 + 1) # y
z = 30 * gauss2(xx, yy, [2D, 3D, .2, 1]) + $
    10 * gauss2(xx, yy, [-3D, 1D, .2, 1]) + $
    20 * gauss2(xx, yy, [1D, -2D, .2D, 1]) 
zi = floor(z)   ;; Convert to integer

;; Find the positions of significant data points
wh = where(z GT 5, ct)
if ct EQ 0 then message, 'ERROR: no signif points!'
xi = x(wh MOD 100)
yi = y(floor(wh/100))
xy = transpose([[xi],[yi]])
weights = clust_wts(xy, n_clusters=3)

plot, xi, yi, psym=3
oplot, weights(0,*), weights(1,*), psym=1, symsize=3


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