PRO cluster

 ; This routine contains the computation of the HOPKINS STATISTIC and the CLUSTER INDEX
 ; used in Ruthazer, E.S., and Stryker, M.P. (1996) The role of activity in the development of 
 ; long-range horizontal connections in area 17 of the ferret. J Neurosci 16: 7253-69
 ; The original paper concerning this statistic is Hopkins, B. (1954) A new method for determining 
 ; the type of distribution of plant individuals. Ann Bot (Lond) 18:213–227.
 ; It is cited in standard references like Ripley, B.D. (1981) Spatial statistics. New York: Wiley. 
 ;
 ; The program is written in the language called IDL (Interactive Data Language), 
 ; which is similar to Matlab in spirit and is now used for geological and medical image analysis
 ; www.ittvis.com/ProductServices/IDL.aspx
 ; http://en.wikipedia.org/wiki/IDL_(programming_language)
 ;
 ; The code was written a long time ago by Ed Ruthazer with some help by Michael Stryker.  
 ; It is commented fairly extensively, so it should probably be clear after you read the 
 ; Ruthazer & Stryker (1996) paper.
 ; Unfortunately, no help is available.  If you do not have IDL on your computer, you should
 ; find it fairly easy to port to a similar language like Matlab or Octave, if you use the
 ; online documentation that describes some of the more cryptic IDL functions, such as that at 
 ; http://physics.nyu.edu/grierlab/idl_html_help/home.html
 ;
 ; 
 
 
 
;THIS IS AN IDL ROUTINE WHICH:
;1) READS IN A FILE OF X,Y COORDINATES OF CELL BODY POSITIONS
;           RELATIVE TO AN INJECTION SITE AT (0,0).  THE FORMAT IS
;           ONE PAIR OF SPACE-DELIMITED X & Y VALUES (IN MICRONS) PER
;           LINE. UP TO A MAXIMUM OF 10000 LINES (10000 CELLS)
;2) SETS A REGION OF INTEREST(ROI) IN WHICH TO MAKE SPATIAL
;           STATISTICAL MEASUREMENTS THAT ENCOMPASSES THE SET OF
;           CELLS AND EXCLUDES A REGION AROUND THE INJECTION SITE.
;3) IN A SLIDING ANALYSIS WINDOW (I USED 1MM X 1MM, SLIDING
;           BY 100UM EACH TIME.), IT MEASURES THE NEAREST-NEIGHBOR
;           DISTANCES FOR 10% OF THE CELLS(W) AND THE NEAREST-NEIGHBOR
;           DISTANCES FOR AN EQUAL NUMBER OF RANDOMLY DROPPED POINTS
;           AND THEIR NEAREST CELLS(X). IT THEN GIVES A CLUSTER INDEX VALUE
;           VALUE (H) BASED ON THE FORMULA:
;           H=LOG(SUM(X^2)/SUM(W^2))
;4) IT PERFORMS THIS CALCULATION REPEATEDLY (AT LEAST 10 TIMES).
;5) IT DISPLAYS THE RESULTS IN GRAPHICAL AND NUMERICAL FORM AND
;           OFFERS VARIOUS WAYS OF SAVING THE RESULTS.
 
COMMON HOP,hopkins
 
pt_set=fltarr(2000,2)         ;array of point coordinates in quadrat
in_set=fltarr(10000,2)       ;array of point coordinates in full set
n_neigh=fltarr(2000)         ;array of nearest neighbor distances for PT_SET points
n_neigh_sq=fltarr(2000)
rand_neigh=fltarr(2000)     ;array of nearest neighbor distances for random + PT_SET points
p=bytarr(1)                     ; a hack to allow TV to plot a single pixel
select = intarr(2000)
result=fltarr(6)
filenm = ' '
read,filenm,PROMPT='Input Filename:'
in_set=reader(filenm)        ;read datafile as (x,y) pairs and divide coordinates by 10
filenm = strmid(filenm,rstrpos(filenm,'/')+1,strlen(filenm)) ; REMOVE DIRECTORY NAMES
 
alleng=non0el(in_set)       ;get length of dataset excluding (0,0) filler data
 
;THESE VARIABLES ARE MIN AND MAX FOR ORIGINAL COORDINATE DATA(/10)(I.E., BOUNDS)
quadxmax=max(in_set(*,0))
quadymax=max(in_set(*,1))
quadxmin=min(in_set(*,0))
quadymin=min(in_set(*,1))
 
;MAKE A WINDOW THE DIMENSIONS OF THE DATASET AND PLOT DATA.
window,0,xsize=quadxmax-quadxmin+50,ysize=quadymax-quadymin+50,retain=2
plot,in_set(*,0),in_set(*,1),psym=3, title = filenm
print,quadxmax,quadxmin,quadymax,quadymin
 
yn='no'
read,yn,prompt='Save plot as postscript file?'
if yn eq'yes' then begin
  ps_open,'P'+filenm
            plot,in_set(*,0),in_set(*,1),psym=3, title = filenm
  ps_close
endif
 
;GET ANALYSIS INFO FROM USER.
read,quadradius,PROMPT='Quadrat Radius (x10):'
read,quadinc,PROMPT='Boxcar increment (x10):'
 
;CREATE TWO PIXMAP ARRAYS.  ONE OF DATA AND ONE TO BE ROI TEST ARRAY (acc_arr)
picture=bytarr(quadxmax-quadxmin+2*quadradius+quadinc,quadymax-quadymin+2*quadradius+quadinc)
acc_arr=bytarr(quadxmax-quadxmin+2*quadradius+quadinc,quadymax-quadymin+2*quadradius+quadinc)
picture(*,*)=0
 
in_set(*,0)=in_set(*,0)-quadxmin+quadradius;make all coordinates positive and
in_set(*,1)=in_set(*,1)-quadymin+quadradius;leave room for quadrat box around roi.
 
for i=0,alleng-1 do picture(in_set(i,0),in_set(i,1))=1 ; copy data into pixmap
 
;DEFINE THE ROI TEST ARRAY USING DOROI(MOUSE INTERFACE AND SHRINKS DATA
;            TO WORKABLE SIZE, THEN RE-EXPANDS IT)
window,1,xsize=quadxmax-quadxmin+50,ysize=quadymax-quadymin+50,retain=2
acc_arr=doroi(picture,quadxmax-quadxmin+2*quadradius+quadinc,quadymax-quadymin+2*quadradius+quadinc)
 
;MAKE A HOLE IN THE CENTER OF THE ROI TEST ARRAY TO EXCLUDE INJECTION SITE.
read,roe,PROMPT='Exclusion Radius (x10):'
xorigin = -quadxmin+quadradius
yorigin = -quadymin+quadradius
for x=xorigin-roe,xorigin+roe do for y=yorigin-roe,yorigin+roe do $
     if (sqrt( ((x-xorigin)^2) + ((y-yorigin)^2) ) le roe) then acc_arr(x,y)=0
picture=bytarr(quadxmax-quadxmin+2*quadradius+quadinc,quadymax-quadymin+2*quadradius+quadinc)
for i=0,alleng-1 do picture(in_set(i,0),in_set(i,1))=1+acc_arr(in_set(i,0),in_set(i,1))
picture=(acc_arr+picture)*70
 
;REDRAW DATA SET
window,1,xsize=quadxmax-quadxmin+50,ysize=quadymax-quadymin+50,retain=2
tvscl,picture
yn='no'
read,yn,prompt ='Save roi-image as gif file?'
if yn eq'yes' then write_gif,filenm+'_roi.gif',byte(picture)
 
;HOPKINS IS ARRAY OF HOPKINS NUMBER PER QUADRAT.
read,depth,prompt='iterations per quadrat (at least 10): '
hopkins = fltarr(1+((quadxmax-quadxmin)/quadinc),1+((quadymax-quadymin)/quadinc),depth+1)
sigrid = fltarr(1+((quadxmax-quadxmin)/quadinc),1+((quadymax-quadymin)/quadinc))
hopkins(*,*,*)=0
sigrid(*,*)= 0
 
;DEFINE SMALLEST RECTANGLE THAT ENCLOSES ROI ([xmin,ymin],[xmax,ymax])
 
xmin=0  ;COUNT UP FROM LEFT
repeat begin
xmin=xmin+1
endrep until (max(acc_arr(xmin,*)) ne 0)
 
ymin=0  ;COUNT UP FROM BOTTOM
repeat begin
ymin=ymin+1
endrep until (max(acc_arr(*,ymin)) ne 0)
 
ymax=quadymax-quadymin            ;COUNT DOWN FROM TOP
repeat begin
ymax=ymax-1
endrep until (max(acc_arr(*,ymax)) ne 0)
 
xmax=quadxmax-quadxmin            ;COUNT DOWN FROM RIGHT
repeat begin 
xmax=xmax-1
endrep until (max(acc_arr(xmax,*)) ne 0)
 
;QUADRAT CENTERS AT [quadx,quady]
for quady=ymin,ymax,quadinc do begin
print,string(fix((quady-ymin)/(ymax-ymin)*100))+'%'
 for quadx=xmin,xmax,quadinc do begin
    ;pt_set IS SUBSET OF in_set, CONTAINED IN RECTANGLE DEFINED BY FUNCTION CALL.
    pt_set=all2quad(alleng,in_set,[quadx-quadradius,quady+quadradius],$
                                                [quadx+quadradius,quady-quadradius],acc_arr)
 
    leng=non0el(pt_set)      ;LENG IS NUMBER OF NON-ZERO COORDINATES, REST IS IGNORED.
           
;MEASURE NEAREST NEIGHBOR DISTANCES FOR EACH POINT (IF THERE'S MORE THAN 1 POINT 
    if (leng/10 ge 3) then begin      
            n_neigh(*)=0
            rand_neigh(*)=0
           
            ;THIS FILLS THE NEAREST NEIGHBOR ARRAY n_neigh
            for i=0,leng-1 do begin
                        n_neigh(i)=nearest(pt_set,i,leng)
                        ;print,'real:',pt_set(i,0),pt_set(i,1),n_neigh(i)
                        end;for i
           
            ;NOW SELECT leng RANDOM COORDINATES IN QUADRAT, CONFIRM IN ROI, FILL ARRAY WITH DISTANCES.
            for h=0,depth-1 do begin 
               for i=0,(leng/10)-1 do begin
                        repeat begin
                                    pt_set(leng,0)=quadx-quadradius + randomu(seed)*quadradius*2
                                    pt_set(leng,1)=quady-quadradius + randomu(seed)*quadradius*2
                        endrep until (acc_arr(pt_set(leng,0),pt_set(leng,1)) ne 0)
                        rand_neigh(i)=nearest(pt_set,leng,leng)
                        end;for i
 
                        ;print,quadx,quady,total(n_neigh)/leng,total(rand_neigh)/leng,leng
;RANDOMLY SELECT 1 OUT OF 10 POINTS
                        select(*)=0
                        for i = 0,(leng/10)-1 do select(i)= fix(randomu(seed)*leng)
;CALCULATE LOG(HOPKINS NUMBER) WHICH IS: LOG(SUM(RAND_NEIGH^2)/SUM(N_NEIGH^2) )
                        n_neigh_sq=n_neigh^2
                        rand_neigh=rand_neigh^2
                        hopkins(quadx/quadinc,quady/quadinc,h)=alog(total(rand_neigh)/total(n_neigh_sq(select(0:(leng/10)-1))))
                        ;print,'Hopkins:',hopkins(quadx/quadinc,quady/quadinc,h),(leng/10)
                        ;wset,1
                        ; xyouts,quadx+quadradius,quady+quadradius,fix(hopkins(quadx/quadinc,quady/quadinc,h)*100),$
                        ;                       color=hopkins(quadx/quadinc,quady/quadinc,h)*25,alignment=0.5
           
                        ;wset,0
                        ;p(0)=(hopkins(quadx/quadinc,quady/quadinc,h)*50)+128 ; TV NEEDS AN ARRAY NOT A SCALAR
                        ;if (p(0) ge 255) then p(0)=254
                        ;for i=0,leng-1 do tv,p(*),pt_set(i,0),pt_set(i,1);  PLOT PIXELS OF QUADRAT  
            end;for h
            result = moment(hopkins(quadx/quadinc,quady/quadinc,0:depth-1),sdev = temp)
            hopkins(quadx/quadinc,quady/quadinc,depth)=result(0)
            sigrid(quadx/quadinc,quady/quadinc) = temp
            print,"H=",result(0),"+/-",temp
   endif
end;for quadx
end; for quady
 
;CREATE PIXEL ARRAY OF HOPKINS NUMBERS BY QUADRAT POSITION and SCALE SIZE & INTENSITY
picture=congrid(hopkins(*,*,depth),quadxmax-quadxmin+2*quadradius,quadymax-quadymin+2*quadradius)
picture = (picture *50) +128
index=where(picture ge 254,count)
if count NE 0 then picture(index)=254         ;PREVENT TV ROUTINE FROM WRAPPING VALUES >255
read,w,PROMPT='display cluster index in window #:'            ;OPEN NEW WINDOW TO PLOT IMAGE AND THEN PLOT IT
window,w,xsize=quadxmax-quadxmin+100,ysize=quadymax-quadymin+100,retain=2
tv,picture
 
yn='no'
read,yn,prompt ='Save cluster index array as gif file?'
if yn eq'yes' then write_gif,filenm+'_hop.gif',byte(picture(*,*))
 
yn='no'
read,yn,prompt='Save cluster index array as raw data?'
if yn eq 'yes' then begin
            openw,1,filenm+'.idl'
            printf,1,1+((quadxmax-quadxmin)/quadinc),1+((quadymax-quadymin)/quadinc),depth+1
            printf,1,hopkins(*,*,*)
            close,1
            endif
 
;CREATE PIXEL ARRAY OF SIGNIFICANCE BY QUADRAT POSITION and SCALE SIZE & INTENSITY
for quady=ymin,ymax,quadinc do for quadx=xmin,xmax,quadinc do begin
   if (sigrid(quadx/quadinc,quady/quadinc) ne 0) then sigrid(quadx/quadinc,quady/quadinc) = t_pdf($
             hopkins(quadx/quadinc,quady/quadinc,depth)/(sigrid(quadx/quadinc,quady/quadinc)/sqrt(depth)),depth-1)
   if (hopkins(quadx/quadinc,quady/quadinc,depth) le 0) then sigrid(quadx/quadinc,quady/quadinc) = .5
   ;xyouts,quadx,quady,string(sigrid(quadx/quadinc,quady/quadinc))
end ;for quadx
 
picture=congrid(sigrid(*,*),quadxmax-quadxmin+2*quadradius,quadymax-quadymin+2*quadradius)
picture= (picture)
index1=where(picture ge 0.5,count1)
index2=where(picture ge 0.95,count2)
if count1 NE 0 then picture(index1)=100
if count2 NE 0 then picture(index2)=200     
read,w,PROMPT='plot clustered sites (p=0.95)in window #';OPEN NEW WINDOW TO PLOT IMAGE AND THEN PLOT IT
window,w,xsize=quadxmax-quadxmin+100,ysize=quadymax-quadymin+100
tv,picture
 
read,yn,prompt='Save Significance plot as gif file?'
if yn eq'yes' then write_gif,filenm+'_sig.gif',byte(picture(*,*))
 
;DO SOME STATISTICS ON ARRAY AND PRINT OUT
tmp = hopkins(*,*,depth)
result=moment(tmp(where(tmp ne 0,count)))
Print,'mean:',result(0),' ± ',string(sqrt(result(1))),'stdev     N=',count
print,'median:', median(tmp(where(tmp ne 0)))
print,'p > ',1-t_pdf(result(0)/sqrt(result(1)/count),count)
 
;DRAW A CUMULATIVE HISTOGRAM OF QUADRAT HOPKINS VALUES
read,w,PROMPT='histogram in window #:'
window,w,retain=2
c_hist=histogram(tmp(where(tmp ne 0)),min = -4, max =4,binsize=.01)
for i=1,n_elements(c_hist)-1 do c_hist(i)=c_hist(i)+c_hist(i-1)
plot,float(c_hist)/float((c_hist(n_elements(c_hist)-1))),xtickname=[-4,-2,0,2,4]
xyouts,5,0.8,filenm
xyouts,5,0.75,'mean:'+strcompress(result(0))+'+/-'+strcompress(sqrt(result(1)))
xyouts,5,0.7,'median:'+strcompress(median(tmp(where(tmp ne 0))))+' N='+strcompress(count)
;xyouts,5,0.65,'p ='+strcompress(1-t_pdf(result(0)/sqrt(result(1)/count),count))
 
read,yn,prompt='Save histogram as postscript file?'
if yn eq'yes' then begin
  ps_open,'H'+filenm
   plot,float(c_hist)/float((c_hist(n_elements(c_hist)-1))),xtickname=[-4,-2,0,2,4]
   xyouts,5,0.8,filenm
   xyouts,5,0.75,'mean:'+strcompress(result(0))+'+/-'+strcompress(sqrt(result(1)))
   xyouts,5,0.7,'median:'+strcompress(median(tmp(where(tmp ne 0))))+' N='+strcompress(count)
  ps_close
endif
 
END ;pro ripley
 
;  THIS FUNCTION READS IN LISTS OF X,Y COORDINATES FROM THE FILENAME
;  IN THE ARGUMENTS AND RETURNS A 10000
;  ELEMENT ARRAY OF X,Y PAIRS
 
function reader,filenm
   tmp=fltarr(2)
   a=fltarr(10000,2)
   openr,1,filenm
        print,filenm
        for i=0,9999 do $
        if not eof(1) then begin
                readf,1,tmp
                a(i,0)=tmp(0)/10
                a(i,1)=tmp(1)/10
            endif
        close,1
return,a
end
 
;THIS FUNCTION RESCALES THE X,Y PLOT OF POINTS IF IT IS LARGER THAN
;THE NORMAL SIZE OF AN IDL WINDOW (FEED IT XSIZE AND YSIZE) AND THEN
;ASKS THE USER TO DRAW A REGION OF INTEREST ON THE PLOT.  THE RIGHT
;BUTTON CLOSES THE ROI.  A BYTE ARRAY IS RETURNED WHICH HAS THE SAME
;DIMENSTIONS AS THE XY PLOT AND CONTAINS 0 OUTSIDE THE ROI AND 1
;INSIDE THE ROI.
 
function doroi,picture,xsize,ysize
            rebxs=fix(xsize+.5)
            rebys=fix(ysize+.5)
 
            if ((xsize ge 400) or (ysize ge 400)) then $
             begin
             rebxs=fix((xsize/2)+.5)
             rebys=fix((ysize/2)+.5)
             picture = congrid(picture,rebxs,rebys)
            end
 
            tvscl,picture
            x=defroi(rebxs,rebys)
            roi_arr=bytarr(rebxs,rebys)
            roi_arr(*,*)=0
            roi_arr(x)=1
            if (xsize ne rebxs) then begin
              tmp=bytarr(xsize,ysize)
              tmp=congrid(roi_arr,xsize,ysize)
              roi_arr=tmp
            end
  return,roi_arr
end   
 
;THIS FUNCTION TAKES THE FULL DATASET ARRAY (in_set) OF leng ELEMENTS,
;AND RETURNS A SMALLER ARRAY CONTAING ONLY THOSE MEMBERS WITHIN THE ROI
;INSIDE THE SQARE DEFFINED BY up_left AND low_right.
;(ROI IS DEFINED BY acc_arr)
function all2quad,leng,in_set,up_left,low_right,acc_arr
  out_set=fltarr(2000,2)
  count=0
  for i=0,leng-1 do $
    if (in_set(i,0) GE up_left(0)) and (in_set(i,0) LE low_right(0)) and $
       (in_set(i,1) LE up_left(1)) and (in_set(i,1) GE low_right(1)) then $
            if (acc_arr(fix(in_set(i,0)),fix(in_set(i,1))) ne 0) then $
            begin
           out_set(count,0) = in_set(i,0)
           out_set(count,1) = in_set(i,1)
           count=count+1
    end
 return,out_set
end