Internal: A Modal Filter Function for Classic ENVI
Anonym
Modal filtering is commonly used to clean up classification results in an attempt to make the results more spatially coherent and realistic. In this type of filter, a rectangular window is passed over each pixel and the DN of the window's central pixel is replaced with the window's statistical mode (i.e., the most common value).
This User Function is a good example of how to use ENVI's status reporting mechanism without tiling the processing. It also provides one method for how to handle image borders for spatial filtering. The status reporting routines used are:
ENVI_REPORT_INIT
ENVI_REPORT_INC
ENVI_REPORT_STAT
Click here to download the user function code modal_filter.pro as an ASCII file.
;--------------------------------------------------------+
; NAME
;
; MODAL_FILTER
;
; TYPE
;
; ENVI User Function
;
; CALLING SEQUENCE
;
; MODAL_FILTER, event
;
; PURPOSE
;
; Modal Filtering replaces the central pixel
; within the filtering window with the
; neighborhood mode (i.e., most common value).
; This is typically done to clean up classification
; results, although this function allows any type
; of file as input.
;
; MODIFICATION HISTORY
;
; July, 1998
; Writen under IDL 5.0.3a, ENVI 3.0
;
; December 1999
; Updated to handle image borders and to
; be more memory efficient
;
; AUTHOR INFO
;
; David Gorodetzky
; RSI Technical Support and Profession Services Group
; goro@rsinc.com
;
; PROGRAM FILE ORGANIZATION
;
; This file contains the following procedures and
; functions in this order:
;
; NEIGHBORHOOD
;
; A function that returns the (sample, line)
; pixel addreses of the filer window pixels
; in an array of (2,n_pixels). This function
; also removes invalid addresses -- those that
; are outside the image boundaries.
;
; MODAL_FILTER
;
; The widget creation routine (that uses ENVI's
; auto-managed widget events) and the image
; processing routine.
;
; PROGRAM NOTES
;
; Currently input is limited to only a single image band,
; and tiling is not used so the entire input image must
; fit into memory.
;
;---------------------------------------------------------------+
function NEIGHBORHOOD, WinX, WinY, samples, lines, S, L
; returns the filter window's pixel addresses
;
; WinX and WinY: the X and Y filter window sizes
; samples and lines: 2D arrays that are used to
; compute the pixel addresses
; S and L: the address of the filter window's
; central pixel
;
n_addr = LONARR(2,winX*winY)
n_samples = samples + S
n_lines = lines + L
counter=0L
FOR i=0,(winX*winY)-1 DO BEGIN
IF (n_samples[i] ge 0 AND n_lines[i] ge 0) THEN BEGIN
n_addr[*,counter] = [ n_samples[i], n_lines[i] ]
counter=counter+1
ENDIF
ENDFOR
; trim any extra space off the array
;
n_addr = n_addr[*,0:counter-1]
; return the valid neighborhood addresses
;
return, n_addr
END
;---------------------------------------------------------------------+
pro MODAL_FILTER, event
; select a file for processing
;
ENVI_SELECT, title='modal filter input file', fid=fid, $
pos=pos, dims=dims, /band_only
IF fid EQ -1 THEN return
; get required file info
;
ENVI_FILE_QUERY, fid, fname=fname, data_type=dt, $
file_type=FT, bnames=orig_bname, ns=orig_ns, $
nl=orig_nl, xstart=xs, ystart=ys, sname=sname
; if the image was spatially subsetted,
; update the xstart and ystart values
;
IF dims(1) NE 0 THEN xs = xs+dims(1)
IF dims(3) NE 0 THEN ys = ys+dims(3)
; compute actual dimensions which may differ from
; the header info b/c of a spatial subset
;
ns = dims[2]-dims[1] +1
nl = dims[4]-dims[3] +1
; default text string for WIDGET_OUTFM below
; (must be defined before the catch statement)
;
txt = 'modal_f_'+sname
; catch any errors to prevent a crash
;
Catch, error
IF (error NE 0) THEN BEGIN
ok = DIALOG_MESSAGE(!error_state.msg, /cancel)
IF n_elements(lun) ne 0 THEN Free_LUN, lun
IF n_elements(base) ne 0 THEN $
ENVI_REPORT_INIT, base=base, /finish
IF (STRUPCASE(ok) EQ 'CANCEL') THEN return
ENDIF
; set up the User Interface for collecting
; the rest of the required parameters
;
TLB = WIDGET_AUTO_BASE(title='Modal Filter Parameters')
sbase1 = WIDGET_BASE(TLB, /row, /frame)
; determine the max window size to
; be used as the 'ceiling' for WIDGET_PARAM,
; and make sure they're odd numbers
;
max_sam = ns/10
max_line =nl/10
IF (max_sam MOD 2 EQ 0) THEN max_sam = max_sam-1
IF (max_line MOD 2 EQ 0) THEN max_line = max_line-1
sam = WIDGET_PARAM(sbase1, prompt='Size of Filter Window: samples', $
ceil=max_sam, floor=3, dt=2, default=5, xsize=3, $
uvalue='samples', /auto)
lin = WIDGET_PARAM(sbase1, prompt='lines', ceil=max_line, $
floor=3, dt=2, default=5, xsize=3, uvalue='lines', /auto)
; only add this section to the interface
; if the data type of the input band is
; some sort of floating point
;
IF (dt GT 3) THEN BEGIN
sbase2 = WIDGET_BASE(TLB, /col, /frame)
lable2 = WIDGET_LABEL(sbase2, value='Set Increment for Modal Precision',$
/align_left)
lable3 = WIDGET_LABEL(sbase2, /align_left, $
value='(i.e., the binsize for determining the mode)')
sbase3 = WIDGET_BASE(sbase2, /row)
bin = WIDGET_PARAM(sbase3, prompt='binsize', field=3, default=.1, $
xsize=4, uvalue='binsize', /auto)
ENDIF
file = WIDGET_OUTFM(TLB, uvalue='outfile', /frame, default=txt, /auto)
; let ENVI auto-manage the widgets
;
result = AUTO_WID_MNG(TLB)
IF (result.accept EQ 0) THEN return
; make sure the filter sizes are odd --
; if not, invoke in error, but keep
; the previously defined output filename
; to use as the default for WIDGET_OUTFM
;
IF ( result.samples MOD 2 EQ 0 OR $
result.lines MOD 2 EQ 0 ) THEN BEGIN
IF result.outfile.in_memory NE 1 THEN txt=result.outfile.name
MESSAGE, 'The modal filter window sizes must be odd.'
ENDIF
; if results are being saved to memory, make sure
; you're not going to exceed the cache limit
;
mem_res = MAGIC_MEM_CHECK(fid=fid, out_dt=dt, nb=1, dims=dims, $
out_name=result.outfile.name, in_memory=result.outfile.in_memory)
IF (mem_res.cancel eq 1) THEN return
; just to make the code below more readable
;
winX = FIX(result.samples) ; x window size
hwX = winX/2 ; the 'odd' half in the x direction
winY = FIX(result.lines) ; y window size
hwY = WinY/2 ; the 'odd' half in the y direction
; set up the INHERIT structure so that
; if the input band is a class image,
; the filtered image will retain the
; class mapping info; otherwise assume
; its normal image data, so inherit
; the band and wavelength information
;
inherit = {fid:fid, flag:0B, POS:pos}
IF (FT EQ 3) THEN inherit.flag=4 ELSE inherit.flag=3
; get the input data
;
input = ENVI_GET_DATA(dims=dims, fid=fid, pos=pos[0])
; define the modal filter binsize
;
IF (dt LT 4) THEN bin=1 ELSE bin=result.binsize
; initilize the status report
;
IF (mem_res.in_memory EQ 1) THEN oname = 'In Memory' $
ELSE oname = mem_res.out_name
outstr = ['Input Band :' +fname+'('+orig_bname(pos)+')',$
'Output File: '+oname]
ENVI_REPORT_INIT, outstr, title='Modal Filter Progress', $
base=base, /interupt
; initilize the increment counter
;
ENVI_REPORT_INC, base, nl
; open the output file for writting,
; or make an array for the filtered results
; (make the temp array for writting to file
; b/c its an easy way to match the data
; type of the input image)
;
IF (mem_res.in_memory eq 0) THEN BEGIN
temp = MAKE_ARRAY(ns, type=dt)
OpenW, lun, mem_res.out_name, /Get_LUN
ENDIF ELSE outarray = MAKE_ARRAY(ns, nl, type=dt, /no)
; In order to filter the entire image, you
; have to pay special attention to the edges
; since the filter window is going to include
; pixels outsize of the image boundary. For
; the image bounradies, use the NEIGHBORHOOD
; function included here to get filter window's
; pixel adresses, since it excludes the invlaid
; addresses. So, prepare the inputs for
; using this function.
;
samples = INTARR(WinX,WinY)
lines = samples
FOR i=0,winY-1 DO samples[*,i] = indgen(winX) - hwX
FOR i=0,WinY-1 DO lines[*,i] = (indgen(winY) - hwY)[i]
; filter the image -- this is the image processing
;
FOR l=0,nl-1 DO BEGIN ; for each line
;
; update the status report dialog
;
ENVI_REPORT_STAT, base, l, nl, cancel=cancel
IF (cancel) THEN BEGIN
ENVI_REPORT_INIT, base=base, /finish
IF n_elements(lun) ne 0 THEN Free_LUN, lun
return
ENDIF
;
; Replace each sample in the line with its
; neighborhood mode. For the edge pixels
; use the neighborhood fxn to get the valid
; adresses. For all other pixels directly
; subscript the input data array since its
; much faster. Data is written to file only
; once per line since excessive file I/O is
; pretty slow.
;
FOR s=0,ns-1 DO BEGIN
IF (l lt hwY OR abs(l-(nl-1)) lt hwY OR $
s lt hwX OR abs(s-(ns-1)) lt hwX) THEN BEGIN
hood = neighborhood(WinX, WinY, samples, lines, s, l)
filter_data = input[ hood[0,*], hood[1,*] ]
ENDIF ELSE $
filter_data = input[s-hwX:s+hwX, l-hwY:l+hwY]
hist = HISTOGRAM(filter_data, binsize=bin, omin=omin)
mode = WHERE(hist EQ MAX(hist))
IF (mem_res.in_memory eq 0) THEN $
temp[s] = omin+(mode[0]*bin) $
ELSE outarray[s,l] = omin+(mode[0]*bin)
ENDFOR
IF (mem_res.in_memory eq 0) THEN WriteU, lun, temp
ENDFOR
; set up some output variables
;
descrip = 'modal filter result for '+sname+', band = '+orig_bname[pos]
outbnames = 'modal filter ('+sname+')'
; take care of the output
;
IF (mem_res.in_memory EQ 1) THEN $
ENVI_ENTER_DATA, outarray, bnames=outbnames, $
descrip=descrip, inherit=inherit, xstart=xs, ystart=ys $
ELSE BEGIN
Free_LUN, lun
ENVI_SETUP_HEAD, fname=mem_res.out_name, ns=ns, nl=nl, $
nb=1, data_type=dt, interleave=0, descrip=descrip, $
bnames=outbnames, xstart=xs, ystart=ys, inherit=inherit, $
/open, /write
ENDELSE
; delete the status reporting dialog
;
ENVI_REPORT_INIT, base=base, /finish
END