| Re: HELP: Density of Points on Scatterplot [message #404] |
Thu, 27 August 1992 06:56 |
thompson
Messages: 584
Registered: August 1991
|
Senior Member |
|
|
In article <1992Aug27.075540@highwire.gsfc.nasa.gov>,
burel@highwire.gsfc.nasa.gov (Jonathan Burelbach) writes...
> Help! I am trying to create a scatterplot of 2 Landsat Image bands (ie
> band 4 vs band 5) for clustering purposes and I need to be able to determine
> the density of points on the plot. I have tried to create a 256x256 z array
> by looping through the images, but this takes forever with images of any
> size. Does anyone have any suggestions?
This may do what you want. It uses HISTOGRAM to cut down on some of the
looping.
Bill Thompson
------------------------------------------------------------ -------------------
PRO FORM_HISTO2,X,Y,HISTO,XSTEPS,YSTEPS,XDELTA,YDELTA
;+
; NAME:
; FORM_HISTO2
; PURPOSE:
; Forms a two-dimensional histogram from a set of X,Y points.
; CALLING SEQUENCE:
; FORM_HISTO2, X, Y, HISTO, XSTEPS, YSTEPS [, XDELTA, YDELTA ]
; INPUT PARAMETERS:
; X, Y = Arrays giving the X,Y points to form the histogram from.
; OPTIONAL INPUT PARAMETERS:
; XDELTA = Spacing between histogram bins in the X direction. If not
; passed, then a suitable value is selected automatically.
; YDELTA = Spacing between histogram bins in the Y direction.
; OUTPUT PARAMETERS:
; HISTO = Two-dimensional array containing the histogram values.
; XSTEPS = Bin coordinate values along the X axis.
; YSTEPS = Bin coordinate values along the Y axis.
; OPTIONAL KEYWORD PARAMETERS:
; None.
; COMMON BLOCKS:
; None.
; SIDE EFFECTS:
; None.
; RESTRICTIONS:
; X and Y must have the same number of points.
; PROCEDURE:
; The number of points within bins bounded by XSTEP(i:I+1), YSTEP(I:I+1)
; are counted and stored into HISTO.
; MODIFICATION HISTORY:
; William Thompson, April 1992, incorporated into SERTS library.
;-
;
ON_ERROR,2
;
; Check the number of parameters.
;
IF (N_PARAMS(0) NE 5) AND (N_PARAMS(0) NE 7) THEN BEGIN
PRINT,'*** FORM_HISTO2 must be called with 5 or 7 parameters:'
PRINT,' X, Y, HISTO, XSTEPS, YSTEPS [, XDELTA, YDELTA ]'
RETURN
ENDIF
;
; Get the minimum and maximum values of X and Y.
;
BANG_C = !C
XMIN = MIN(X,MAX=XMAX)
YMIN = MIN(Y,MAX=YMAX)
!C = BANG_C
;
; If passed, then check the value of XDELTA.
;
IF N_PARAMS(0) EQ 7 THEN BEGIN
IF N_ELEMENTS(XDELTA) NE 1 THEN BEGIN
MESSAGE,'XDELTA must be scalar'
END ELSE IF XDELTA LE 0 THEN BEGIN
MESSAGE,'XDELTA must be positive'
ENDIF
;
; If XDELTA was not passed, then determine the approximate number of histogram
; levels from the number of elements of X.
;
END ELSE BEGIN
NX = FLOAT(N_ELEMENTS(X))
NX = NX < 100. < (7.*ALOG10(NX) + NX/8.)
;
; Use NX to determine the spacing of the histogram levels. Break this number
; down into mantissa and exponent.
;
XDELTA = (XMAX - XMIN) / (NX - 1)
XPOWER = FIX(ALOG10(XDELTA))
IF XPOWER GT ALOG10(XDELTA) THEN XPOWER = XPOWER - 1
XDELTA = XDELTA / 10.^XPOWER
;
; Ensure that the spacing of the histogram levels is either 1,2 or 5 times
; some power of ten.
;
XVAL = [10,5,2]
XVALUE = 1
FOR I = 0,2 DO IF XVAL(I) GT XDELTA THEN XVALUE = XVAL(I)
XDELTA = XVALUE * 10.^XPOWER
;
; If X is of some integer type (byte, integer or long), then ensure that
; XDELTA is at least one.
;
TYPE = SIZE(X)
TYPE = TYPE(TYPE(0) + 1)
IF ((TYPE EQ 2) OR (TYPE EQ 4) OR (TYPE EQ 16)) THEN $
XDELTA = XDELTA > 1
ENDELSE
;
; Find the nearest multiple of XDELTA which is LE the minimum of X.
; Do the same for the maximum of X.
;
IXMIN = LONG(XMIN / XDELTA)
IXMAX = LONG(XMAX / XDELTA)
IF IXMIN*XDELTA GT XMIN THEN IXMIN = IXMIN - 1
IF IXMAX*XDELTA LE XMAX THEN IXMAX = IXMAX + 1
XMIN = IXMIN * XDELTA
XSTEPS = XMIN + XDELTA * INDGEN(IXMAX - IXMIN + 1)
;
; If passed, then check the value of YDELTA.
;
IF N_PARAMS(0) EQ 7 THEN BEGIN
IF N_ELEMENTS(YDELTA) NE 1 THEN BEGIN
MESSAGE,'YDELTA must be scalar'
END ELSE IF YDELTA LE 0 THEN BEGIN
MESSAGE,'YDELTA must be positive'
ENDIF
;
; If YDELTA was not passed, then determine the approximate number of histogram
; levels from the number of elements of Y.
;
END ELSE BEGIN
NY = FLOAT(N_ELEMENTS(Y))
NY = NY < 100. < (7.*ALOG10(NY) + NY/8.)
;
; Use NY to determine the spacing of the histogram levels. Break this number
; down into mantissa and exponent.
;
YDELTA = (YMAX - YMIN) / (NY - 1)
YPOWER = FIX(ALOG10(YDELTA))
IF YPOWER GT ALOG10(YDELTA) THEN YPOWER = YPOWER - 1
YDELTA = YDELTA / 10.^YPOWER
;
; Ensure that the spacing of the histogram levels is either 1,2 or 5 times
; some power of ten.
;
YVAL = [10,5,2]
YVALUE = 1
FOR I = 0,2 DO IF YVAL(I) GT YDELTA THEN YVALUE = YVAL(I)
YDELTA = YVALUE * 10.^YPOWER
;
; If Y is of some integer type (byte, integer or long), then ensure that
; YDELTA is at least one.
;
TYPE = SIZE(Y)
TYPE = TYPE(TYPE(0) + 1)
IF ((TYPE EQ 2) OR (TYPE EQ 4) OR (TYPE EQ 16)) THEN $
YDELTA = YDELTA > 1
ENDELSE
;
; Find the nearest multiple of YDELTA which is LE the minimum of Y.
; Do the same for the maximum of Y.
;
IYMIN = LONG(YMIN / YDELTA)
IYMAX = LONG(YMAX / YDELTA)
IF IYMIN*YDELTA GT YMIN THEN IYMIN = IYMIN - 1
IF IYMAX*YDELTA LE YMAX THEN IYMAX = IYMAX + 1
YMIN = IYMIN * YDELTA
YSTEPS = YMIN + YDELTA * INDGEN(IYMAX - IYMIN + 1)
;
; Form the histogram.
;
HISTO = FLTARR( N_ELEMENTS(XSTEPS), N_ELEMENTS(YSTEPS) )
FOR J = 0, N_ELEMENTS(YSTEPS)-1 DO BEGIN
Y1 = YSTEPS(J)
Y2 = Y1 + YDELTA
W = WHERE( (Y GE Y1) AND (Y LT Y2) ,N_FOUND)
IF N_FOUND EQ 1 THEN BEGIN
I = LONG( (X(W) - XMIN) / XDELTA )
HISTO(I,J) = 1
END ELSE IF N_FOUND GT 1 THEN BEGIN
XX = X(W)
XMAX = MAX(XX)
IF XMAX EQ XMIN THEN BEGIN
HISTO(0,J) = N_ELEMENTS(XX)
END ELSE BEGIN
HISTO(0,J) = HISTOGRAM(XX,MIN=XMIN,BINSIZE=XDELTA )
ENDELSE
ENDIF
ENDFOR
;
RETURN
END
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
