From: thompson@stars.gsfc.nasa.gov (William Thompson, code 682.1, x2040)
Newsgroups: comp.lang.idl-pvwave
Subject: Re: HELP: Density of Points on Scatterplot
Date: 27 Aug 92 13:56:00 GMT
Organization: NASA/GSFC-Laboratory for Astronomy and Solar Physics
News-Software: VAX/VMS VNEWS 1.4-b1

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