; $Id: filter.pro,v 1.10 1995/07/17 22:49:01 ali Exp $ ; Copyright (c) 1992-1995, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ;------------------------------------------------------ ; procedure SetSliders ;------------------------------------------------------ ; This routine sets the sliders that are data dependent. ;------------------------------------------------------ PRO SetSliders, sliderBase COMMON filterCommon, filterbase, image, imagetrans, view, $ viewsize, order, cutoff, bandwidth, filter, filteredimage, $ filteredspectrum, nviews, bandtype, wcutoffslider, $ wbandwidthslider, filtertype, orderslider ; Determine the dimensions of the selected data. dimensions = SIZE(image) IF (N_PARAMS() NE 0) THEN BEGIN ; Define the width of the sliders as a function of hardware display ; size (VIEWSIZE) and number of plot windows (NVIEWS). VIEWSIZE is ; set at 192 for horizontal hardware display sizes less than 1155 ; and 256 otherwise. sliderWidth = (NVIEWS * VIEWSIZE) + ((NVIEWS + 1) * 20) - 670 ; Define a floating-point slider to adjust the frequency cutoff value. wCutOffSlider = CW_FSLIDER(sliderBase, XSIZE=sliderWidth, VALUE=20.0, $ MINIMUM=0.0, MAXIMUM=dimensions(1)/2., $ UVALUE="FILTERCUTOFF", TITLE="Frequency Cutoff",$ FORMAT='(F6.2)') ; Define a floating-point slider to adjust the bandwidth value. wBandWidthSlider = CW_FSLIDER(sliderBase, XSIZE=sliderWidth, VALUE=10., $ MINIMUM=0, MAXIMUM=dimensions(1)/4., $ UVALUE="FILTERBANDWIDTH", TITLE="Bandwidth", $ FORMAT='(F6.2)') ; Desensitize the frequency bandwidth slider. WIDGET_CONTROL, wBandWidthSlider, SENSITIVE=0 ENDIF $ ELSE BEGIN ; Initialize the frequency cutoff value. cutoff = 20. WIDGET_CONTROL, wCutOffSlider, SET_VALUE=cutoff ; Initialize the frequency bandwidth value. WIDGET_CONTROL, wBandWidthSlider, SET_VALUE=10.0 ENDELSE END ;------------ end of procedure SetSliders ------------- ;------------------------------------------------------ ; procedure LoadData ;------------------------------------------------------ ; This routine gets the name of the new data file to be ; loaded using the procedure getfile. It then loads the ; data, resets the controls, and computes the new FFT of ; the new data. It also erases the old filtered image ; and the old filtered power spectrum. ; ; KEYWORDS: ; NAME - if name is specified, the user is not ; prompted with a choice and name is used in its ; place. ;------------------------------------------------------ PRO LoadData, NEW=newFile COMMON filterCommon ; Load a new signal or image from the "data" subdirectory. GetImage, newdata, DESCRIPTION=datadesc, DIMENSIONS=datadim, /ONE_DIM, $ /TWO_DIM, TITLE="Please Select Sample For Filtering", $ SUBDIRECTORY = ["examples", "data"] ; Was the GetImage successful fileOK = (SIZE(newdata))(0) GT 0 ; If GetImage was successful IF fileOK THEN BEGIN IF (KEYWORD_SET(newdata)) THEN BEGIN IF (datadim(1) NE 1) THEN BEGIN IF (datadim(0) LT (VIEWSIZE / 2)) THEN BEGIN datadim(0) = datadim(0) * (VIEWSIZE / datadim(0)) datadim(1) = datadim(1) * (VIEWSIZE / datadim(1)) image = REBIN(newdata, datadim(0), datadim(1)) ENDIF ELSE $ IF (datadim(0) GT VIEWSIZE) THEN BEGIN datadim = datadim < [ viewsize, viewsize] image = CONGRID(newdata, datadim(0), datadim(1)) ENDIF ELSE $ image = newdata ENDIF ELSE $ image = newdata ENDIF IF !d.n_colors le 256 then image = BYTSCL(image, TOP=!d.n_colors-1) IF (NOT(KEYWORD_SET(newFile))) THEN $ SetSliders ; Compute the forward fast Fourier transform. IF (KEYWORD_SET(image)) THEN $ imagetrans = FFT(image, -1) filteredimage = 0 filteredspectrum = 0 ENDIF END ;------------ end of procedure LoadData ------------- ;------------------------------------------------------ ; procedure ViewEventHndlr ;------------------------------------------------------ PRO ViewEventHndlr, event COMMON filterCommon WIDGET_CONTROL, event.id, GET_UVALUE=viewindex IF (view(viewindex).viewtype NE event.index) THEN BEGIN view(viewindex).viewtype = event.index ResetView, viewindex ENDIF END ;------------- end of procedure ViewEventHndlr -------------- ;------------------------------------------------------ ; procedure ResetFilter ;------------------------------------------------------ ; This procedure computes the correct filter depending ; on the bandpass type. Here a high pass filter is ; one minus the low pass filter. ;------------------------------------------------------ PRO ResetFilter COMMON filterCommon COMMON filterCommon_dist, distfun ; Reset the filtered image and the filtered spectrum since ; changing the filter affects them. filteredimage = 0 filteredspectrum = 0 imagesize = SIZE(image) IF (N_ELEMENTS(distfun) NE N_ELEMENTS(image)) THEN BEGIN ; If single dimensional... IF (imagesize(0) EQ 1) THEN BEGIN ; Compute the Euclidean dist function. distfun = FINDGEN(imagesize(1)) distfun = distfun < (imagesize(1) - distfun) ENDIF $ ELSE $ ; Else 2D case... distfun = DIST(imagesize(1)) distfun(0) = 1d-4 ; Avoid division by 0 ENDIF ; New dist fcn ; Display busy while computing filter. WIDGET_CONTROL, filterbase, /HOURGLASS CASE filtertype OF ; Define Butterworth filter types. 0: BEGIN IF (bandtype EQ 0) THEN $ ; Compute lowpass filter. filter = 1.0 / (1.0d + (distfun / cutoff)^(2 * order)) $ ELSE $ IF (bandtype EQ 1) THEN $ ; Compute highpass filter. filter = 1.0 / (1.0 + (cutoff / distfun)^(2 * order)) $ ELSE BEGIN filter = distfun * distfun - cutoff^2 ;Dist squared zeroes = WHERE(filter EQ 0.0, count) IF count NE 0 THEN $ filter(zeroes) = 1e-4 ;Avoid divide by 0 ; Compute bandreject filter. filter = 1.0 / (1.0 + $ ((distfun * bandwidth) / filter) ^ (2 * order)) IF (bandtype EQ 2) THEN $ ; Compute bandpass filter. filter = 1.0 - filter ENDELSE ENDCASE ; Define Exponential filter types. 1: BEGIN IF (bandtype EQ 0) THEN $ ; Compute lowpass filter. filter = EXP((-(distfun / cutoff)^order)) $ ELSE $ IF (bandtype EQ 1) THEN $ ; Compute highpass filter. filter = 1.0 - EXP((-(distfun / cutoff)^order)) $ ELSE BEGIN ;Bandpass / reject, avoid underflow in EXP fcn ; Compute bandpass filter. filter = (distfun * distfun - cutoff^2) / (distfun * bandwidth) filter = EXP(-(filter ^ (2 * order) < 25)) IF (bandtype EQ 3) THEN $ ; Compute bandreject filter. filter = 1.0 - filter ENDELSE ENDCASE ; Define Ideal Mathematical filter types. 2: BEGIN IF (bandtype EQ 0) THEN $ ; Compute lowpass filter. filter = (distfun LE cutoff) $ ELSE $ IF (bandtype EQ 1) THEN $ ; Compute highpass filter. filter = 1.0 - (distfun LE cutoff) $ ELSE $ IF (bandtype EQ 3) THEN $ ; Compute bandreject filter. filter = ((distfun LT (cutoff - bandwidth / 2.0)) or $ (distfun GT (cutoff + bandwidth / 2.0))) $ ELSE $ ; Compute bandpass filter. filter = ((distfun GE (cutoff - bandwidth / 2.0)) and $ (distfun LE (cutoff + bandwidth / 2.0))) ENDCASE ELSE: MESSAGE, "ResetFilter: Bad Filter Type" ENDCASE FOR i = 0, NVIEWS - 1 DO $ IF (view(i).viewtype GT 1) THEN $ ; Redraw all the windows that depend upon the filter. ResetView, i END ;---------------- end of ResetFilter ------------------ ;------------------------------------------------------ ; procedure ResetView ;------------------------------------------------------ ; This procedure resets the view number passed in. It ; is assumed that the structure associated with the ; passed in view has been altered in view type so it ; needs to be redrawn. In the case of one dimensional ; data, there is no filter intensity image and cross - ; sections of the filter are the same as filter plots ; so they just call those routines. Once the view has ; been updated, the button that controls that view is ; renamed to reflect the new view type. ;------------------------------------------------------ PRO ResetView, windowIndex COMMON filterCommon ; Set the active draw window. WSET, view(windowIndex).windownum ; If the view is not "saved" then erase the contents of the draw window. IF (view(windowIndex).viewtype NE 8) THEN $ ERASE, COLOR=100 imagesize = SIZE(image) CASE view(windowIndex).viewtype OF ; Redraw original signal or image data. 0: BEGIN IF (imagesize(0) GT 1) THEN $ TVSCL, image, (VIEWSIZE - imagesize(1)) / 2, $ (VIEWSIZE - imagesize(2)) / 2 $ ELSE $ PLOT, image, /XSTYLE, /YSTYLE ENDCASE ; Redraw logarithmic power spectrum. 1: BEGIN IF (imagesize(0) GT 1) THEN $ TVSCL, SHIFT(ALOG(ABS(imagetrans)), $ imagesize(1) / 2, imagesize(2) / 2),$ (VIEWSIZE - imagesize(1)) / 2, $ (VIEWSIZE - imagesize(2)) / 2 $ ELSE $ PLOT_IO, SHIFT(ABS(imagetrans), imagesize(1)/2), /XSTYLE ENDCASE ; Redraw filter. 2: BEGIN grid_size = 32 IF (imagesize(0) GT 1) THEN BEGIN thinfilter = CONGRID(filter, grid_size, grid_size) SURFACE, SHIFT(thinfilter, grid_size/2, grid_size/2),$ XSTYLE=4, YSTYLE=4, ZSTYLE=4 ENDIF ELSE $ PLOT, SHIFT(filter, imagesize(1) / 2), /XSTYLE EMPTY ENDCASE ; Redraw filter intensity. 3: BEGIN IF (imagesize(0) GT 1) THEN $ TVSCL, SHIFT(filter, imagesize(1) / 2, imagesize(2) / 2), $ (VIEWSIZE - imagesize(1)) / 2, $ (VIEWSIZE - imagesize(2)) / 2 $ ELSE BEGIN view(windowIndex).viewtype = 2 ResetView, windowIndex ENDELSE ENDCASE ; Redraw filter cross section. 4: BEGIN IF (imagesize(0) GT 1) THEN BEGIN PLOT, SHIFT(filter(*, 0), imagesize(1) / 2), /XSTYLE EMPTY ENDIF ELSE BEGIN view(windowIndex).viewtype = 2 ResetView, windowIndex ENDELSE ENDCASE ; Redraw filter shaded surface. 5: BEGIN grid_size = 48 IF (imagesize(0) GT 1) THEN BEGIN thinfilter = CONGRID(filter, grid_size, grid_size) SHADE_SURF, SHIFT(thinfilter, grid_size/2, grid_size/2), $ XSTYLE=4, YSTYLE=4, ZSTYLE=4 ENDIF ELSE $ PLOT, SHIFT(filter, imagesize(1) / 2), /XSTYLE EMPTY ENDCASE ; Redraw filtered logarithmic power spectrum. 6: BEGIN IF (N_ELEMENTS(filteredspectrum) LE 1) THEN $ filteredspectrum = imagetrans * filter IF (imagesize(0) GT 1) THEN $ TVSCL, SHIFT(ALOG(ABS(filteredspectrum)> 1e-10), $ imagesize(1) / 2, imagesize(2) / 2), $ (VIEWSIZE - imagesize(1)) / 2, $ (VIEWSIZE - imagesize(2)) / 2 $ ELSE $ PLOT_IO, SHIFT(ABS(FLOAT(filteredspectrum)) > 1.0E-12, $ imagesize(1)/2), XSTYLE=1, YSTYLE=1 ENDCASE ; Redraw filtered signal or image. 7: BEGIN IF (N_ELEMENTS(filteredimage) LE 1) THEN $ filteredimage = FLOAT(FFT(imagetrans * filter, 1)) IF (imagesize(0) GT 1) THEN $ TVSCL, filteredimage, (VIEWSIZE - imagesize(1)) / 2, $ (VIEWSIZE - imagesize(2)) / 2 $ ELSE $ PLOT, filteredimage ENDCASE ; Saved view. 8: ELSE: ENDCASE END ;------------------- end of ResetView ----------------- ;------------------------------------------------------ ; procedure FilterEventHndlr ;------------------------------------------------------ ; This is the main event handler for the Filter demo. ;------------------------------------------------------ PRO FilterEventHndlr, event COMMON filterCommon WIDGET_CONTROL, event.id, GET_UVALUE=retval if n_elements(retval) eq 0 then retval = event.value ; Take the following action based on the corresponding event. CASE retval OF "Help" : BEGIN infoText = [ $ "This example demonstrates some of the filtering "+ $ "capabilities of IDL using fourier transformations. "+ $ "The original signals are filtered to remove certain "+ $ "frequencys while letting the others pass through. "+ $ "This filtering takes place in the frequency domain. "+ $ "The fast fourier transform converts signals between "+$ "the time and frequency domains. ", "", $ "In the example, there are multiple views that let "+ $ "you see the filtering process at its various stages. "+ $ "The button on top of each view lets you pick which "+ $ "type of view you want. Pressing the button displays "+ $ "a menu that lets you choose between views. ", "", $ "Below the views are the different settings "+ $ "that specify what type of filter is to be used and the "+ $ "cutoff frequency for filtering and the width of the "+ $ "filter when using bandpass or bandreject filtering. ", "", $ "Different signals are available and can be "+ $ "loaded. Some of the signals are two "+ $ "dimensional and others are single dimensional. ", "", $ "You will notice that when viewing single "+ $ "dimensional data, the filter cross section, filter "+ $ "image and shaded surface of filter views are all "+ $ "plots of the filter. This is because these views "+ $ "do not apply to single dimensional filters."] ShowInfo, TITLE="Filtering Example Information", GROUP=event.top, $ WIDTH=80, HEIGHT=24, INFOTEXT=infoText ENDCASE ; Load new signal or image data. "Open..." : BEGIN LoadData ResetFilter FOR i = 0, NVIEWS - 1 DO $ IF (view(i).viewtype LT 2) THEN $ ResetView, i ENDCASE ; Load pre-defined color table tool. "Colors" : XLOADCT, /SILENT, GROUP=event.top ; Reset the filter if the order value changed. "ORDERSELECTION" : BEGIN WIDGET_CONTROL, ORDERSLIDER, GET_VALUE=temp IF (temp NE order) THEN BEGIN order = temp ResetFilter ENDIF ENDCASE ; Reset the filter if frequency cutoff has changed. "FILTERCUTOFF" : BEGIN WIDGET_CONTROL, wCutOffSlider, GET_VALUE=cutoff ResetFilter ENDCASE ; Reset the filter if frequency bandwidth has changed and also if ; bandwidth matters with current filter type. "FILTERBANDWIDTH" : BEGIN WIDGET_CONTROL, wBandWidthSlider, GET_VALUE=bandwidth ResetFilter END "FILTERSELECTION": BEGIN IF (filtertype NE event.value) THEN BEGIN filtertype = event.value ResetFilter ENDIF ENDCASE "BANDSELECTION": BEGIN IF (bandtype NE event.value) THEN BEGIN bandtype = event.value ResetFilter ; If not bandpass and not bandreject... IF (bandtype NE 2) AND (bandtype NE 3) THEN $ WIDGET_CONTROL, wBandWidthSlider, SENSITIVE=0 $ ELSE $ WIDGET_CONTROL, wBandWidthSlider, /SENSITIVE ENDIF ENDCASE "Exit": BEGIN ; deallocate variables in common. ; Destroy widget hierarchy. WIDGET_CONTROL, event.top, /DESTROY ENDCASE ELSE: ENDCASE END PRO CleanUpFilter, wFilterWindow COMMON filterCommon ; Get the color table saved in the window's user value WIDGET_CONTROL, wFilterWindow, GET_UVALUE=previousState ; Restore the previous color table. TVLCT, previousState.colorTable ; Restore previous margins !X.MARGIN = previousState.xMargins !Y.MARGIN = previousState.yMargins ; Clean up images in common. image = 0 & imagetrans=0 & filter = 0 & filteredimage = 0 filteredspectrum = 0 END ;------------------------------------------------------ ; procedure Filter ;------------------------------------------------------ ; This is the main procedure for the filter demo. ;------------------------------------------------------ PRO Filter, GROUP=group COMMON filterCommon IF XRegistered("Filter") THEN $ ;only one instance at a time GOTO, FINISH ; Get the current color vectors to restore when this application is exited. TVLCT, savedR, savedG, savedB, /GET ; Build color table from color vectors colorTable = [[savedR],[savedG],[savedB]] ; Save the current margins in order to restore it when ; the spiro application is exited xMargins = !X.MARGIN yMargins = !Y.MARGIN ; Save items to be restored on exit in a structure previousState = {colorTable:colorTable, xMargins:xMargins, yMargins:yMargins} ; Load a particular default color table loadct, 0, /SILENT !X.MARGIN = [6, 2] !Y.MARGIN = [2, 2] ; Initialize filter states. filteredimage = 0 filteredspectrum = 0 ; Initialize the filter's characteristics. ;The filter's order. order = 2 ; Butterworth filter. filtertype = 0 ; The cutoff frequency. cutoff = 20. ; The bandwidth. bandwidth = 10. ; Lowpass filter. bandtype = 0 ; Determine the hardware device size in pixels. DEVICE, GET_SCREEN=screendims if screendims(0) ge 1100 then i = 2 $ ;General screen sizing else if screendims(0) ge 800 then i = 1 $ else i = 0 ; Number of views and view size in pixels (device). viewsize = ([ 192, 192, 256])(i) nviews = ([ 2, 3, 3 ])(i) ; Load the initial data. LoadData, /NEW ; If no data was selected, drop out of routine. IF (NOT(KEYWORD_SET(image))) THEN $ GOTO, FINISH ; Define an array of view structures. view = REPLICATE({vinst, viewtype:0b, windowid:0L, windownum:0}, NVIEWS) ; The default view settings that are used for differing numbers of views view.viewtype = ([0,7,2,1,6])(0:nviews-1) ; The types of views we can show viewoptions = [ "Original Signal", "Log Power Spectrum", "Filter Plot", $ "Filter Intensity Image", "Filter Cross Section", $ "Shaded Surface of Filter", "Log Filtered Power Spectrum", $ "Filtered Signal", "Save View" ] ; Define a main filter base. filterBase = WIDGET_BASE(TITLE="Filtering Example", /COLUMN, $ MBAR=mbar, UVALUE=previousState) menubar = CW_PDMENU(mbar, /RETURN_NAME, /MBAR, /HELP, $ ['1\File', '0\Open...', '2\Exit', '1\Edit', '2\Colors', $ '1\Help', '0\Help']) ; Define a sub-base belonging to filterBase. viewbase = WIDGET_BASE(filterBase, SPACE=5, /ROW) ; Define the view windows. FOR viewindex = 0, NVIEWS - 1 DO BEGIN tempbase = WIDGET_BASE(viewbase, /FRAME, /COLUMN, /BASE_ALIGN_CENTER) ; Create the filter method drop list methodDropList = WIDGET_DROPLIST(tempbase, VALUE=viewoptions, $ TITLE='', EVENT_PRO='ViewEventHndlr', $ UVALUE=viewindex) tmp = widget_info(methodDropList, /GEOMETRY) if (tmp.scr_xsize lt VIEWSIZE) then $ WIDGET_CONTROL, SCR_XSIZE=VIEWSIZE ; Set the drop list to the current method WIDGET_CONTROL, methodDropList, $ SET_DROPLIST_SELECT = view(viewindex).viewtype view(viewindex).windowid = WIDGET_DRAW(tempbase, $ XSIZE=VIEWSIZE, YSIZE=VIEWSIZE, RETAIN=2) ENDFOR controlbase = WIDGET_BASE(filterBase, /ROW) ; Define a button group to select filter function. filters = cw_bgroup(controlbase, LABEL_TOP='Filter Function', $ /EXCLUSIVE, /FRAME, /NO_RELEASE, UVALUE='FILTERSELECTION', $ ["Butterworth Filter", "Exponential Filter", "Ideal Filter"]) orderbase = WIDGET_BASE(controlbase, /COLUMN) orderlabel = WIDGET_LABEL(orderbase, VALUE="Filter Order") ; Define a widget slider to select the order of the filter. ORDERSLIDER = WIDGET_SLIDER(orderbase, /VERTICAL, /FRAME, $ VALUE=order, MINIMUM=1, MAXIMUM=6, UVALUE='ORDERSELECTION') ; Define a button group to select the filter type. filtertypes = cw_bgroup(controlbase, LABEL_TOP='Filter Type', $ /EXCLUSIVE, /FRAME, /NO_RELEASE, UVALUE='BANDSELECTION', $ ["Low Pass", "High Pass", "Band Pass", "Band Reject"]) sliderBase = WIDGET_BASE(controlbase, /COLUMN) ; Define the floating-point sliders for cutoff freq. and bandwidth. SetSliders, sliderBase ; Display HOURGLASS during setup. WIDGET_CONTROL, filterBase, /REALIZE, /HOURGLASS ; Select Butterworth as default filter function. WIDGET_CONTROL, filters, SET_VALUE=0 ; Select lowpass as default filter type. WIDGET_CONTROL, filtertypes, SET_VALUE=0 ; Get the drawing area numbers and set the view structure array so ; these values can be referenced when drawing. FOR i = 0, NVIEWS - 1 DO BEGIN WIDGET_CONTROL, view(i).windowid, GET_VALUE=temp view(i).windownum = temp ENDFOR ; Draw the default views that don't rely on the filter and then reset ; the filter which automatically draws the filter's dependent views. ResetFilter ResetView, 0 XManager, "Filter", filterBase, EVENT_HANDLER="FilterEventHndlr", $ GROUP_LEADER=group, CLEANUP="CleanUpFilter" FINISH: END ;================ procedure Filter ====================