pro point_filter2,indata,bw,tol,niter,outdata,outpts
;+
; $Id: point_filter2.pro,v 1.1 2008/03/31 14:36:39 nathan Exp $
; NAME:
;	point_filter2
;
; PURPOSE:
;	This procedure filters out the pixels which are brighter than
;	the area in which they are found. It differs from point_filter.pro
;   	in that jfilter_image.pro is called instead of filter_image.pro.
;
; CALLING SEQUENCE:
;	point_filter,indata,bw,tol,niter,outdata,outpts
;
; INPUTS:
;	indata  = 2-dim array of data to be filtered
;	bw	= width of square filter box - must be ODD, recommended
;		  value 5.
;	tol	= scaling factor to controll how bright the point
;		  is before being replaced, recommended value 7.
;	niter	= number of time to repeat the process
;
; OUTPUTS:
;	outdata = the filtered array
;	outpts	= the value and location of the points removed from
;		  the indata array.  outpts is an (n,3) array.
;		  n is the number of points subtracted; ,0 is the data
;		  value (indata - outdata) with pixel coordinates of
;		  x=,1, and y=,2.
; PROCEDURE:
;	Uses filter_image to calculate the mean, median, and variance
;	over a square box of size bw.  For those point which differ from
;       the mean by more than tol*standard deviation, replace the value
;       with the median.  Finally, fill the outpoint array with the 
;       subtracted values and pixel coordinates.
;
;  MODIFICATION HISTORY:
;	Written by Mike Andrews LASCO/NRL/HUGHES STX  12 Feb 1996.
;       Modified by MDA 20 Mar 1996 to reduce memory requirement by
;		reusing the hold array.
;	Modified by MDA 26 Feb 97 to calculate x and y correctly.
;-
; $Log: point_filter2.pro,v $
; Revision 1.1  2008/03/31 14:36:39  nathan
; moved from adam/programs, was dated Nov 7, 2007
;
if ( 2*(bw/2) eq bw ) then bw = bw+1
;
bw2=bw*bw
fact=float(tol^2)/(bw2-2)
outdata=float(indata)
;
for i=1,niter do begin
;
	med=jfilter_image(outdata,MEDIAN=bw,/ALL)
; mean
	hold=( jfilter_image(outdata,SMOOTH=bw,/ALL) - outdata )/(bw2-1)
; image deviation
	hold=(outdata - temporary(hold))^2
; image variance
	hold=abs( ( jfilter_image(hold,SMO=bw,/ALL)*bw2 - $
		hold ) )
; standard deviation
	hold=sqrt(fact*temporary(hold))
;
;   Select the data where outdata - med >0 (bright points only) 
;   and where this difference squared is larger than hold.
;
	q=where( ( (outdata-med) gt 0)  and ( (outdata-med)^2 gt hold),nq)
;
	if nq lt 1 then goto,done
;
;print, 'At iteration ',i, ' Replacing ',nq,' Points'
	outdata(q)=med(q)
;
    endfor
;
done:
;
;  Fill the outpt array
;
hold=indata-outdata
s=size(hold)
q=where(hold gt 0,nq)
;
	if nq gt 0 then begin
		y= q/s(1) 
		x= q-  y*s(1) 
		outpts=fltarr(nq,3)
		outpts(*,0)=hold(q)
		outpts(*,1)=x
		outpts(*,2)=y
	   endif
;
return
;
end