; MASTER COPY: /net/bali/usr1/spider/procg/cead_test.spi ; ; AUTHOR: ArDean Leith October 2006 ; ; PURPOSE: Studies influence of parameter selection on ; results of SPIDER 'CE AD' operation "HEG" ; Anisotropic Diffusion ; Operates on a single relevant slice from the volume. ; ; OUTPUT: Stack of images created with various parameter settings. ; ; USAGE: Use Web to display central slice of the original input ; volume ; Use Web -- Pixel operation to find location of a good window ; for testing. Window will be 200x200 unless you alter it below. ; Use this routine and input: ; Input volume name, ; Input window location, ; CE AD Iterations- Initial, Final, & Step ; CE AD Times- Initial, Final, & Step ; FC - Contour Levels ; While running - use Web to look at results in: cead_stk ; When finished - use Web to look at final montage in: cead_montage ; When finished - use Web to look at final imges in: cead_volume ; ; tingvol 68,229,121 10,25,5 .05,.15,.05 ; ryrvol 260,50,55 5,100,25 .00,.21,.10 20,110,90 ; benzervol 260,50,66 3,23,10 .01,.13,.06 ; Toggles to select type of input and types of outputs [FR] = 1 ; Use 'FR' Interactive input instead of local file input [PW] = 1 ; Output Power Spectrums [FC] = 1 ; Output Contoured Files [TH] = 0 ; Output Thresholded images [VOL] = 1 ; Output Volume with windows of denoised slices ; ----------------------------------- INPUT VOLUME IF ([FR] .EQ. 0) THEN FR L ; Input volume [volin]input/benzervol ENDIF ; ----------------------------------- OUTPUT FILES IF ([FR] .EQ. 0) THEN FR G ; Labeled output stack file [stkout]cead_stk FR G ; Final labeled output montage [montage]cead_montage FR G ; Final labeled denoised slice volume [volout]cead_denoisedvol ENDIF ; ----------------------------------- INPUT PARAMETERS [wnsam] = 200 ; Window X dimension [wnrow] = 200 ; Window Y dimension [wnslice] = 11 ; Window Z dimension (SHOULD BE ODD!) IF ([FR] .EQ. 0) THEN [wx] = 260 ; Window upper left X location [wy] = 50 ; Window upper left Y location [slice] = 66 ; Relevant slice from input volume [iter_go] = 3 ; Starting iteration [iter_end] = 23 ; Ending iteration [iter_step] = 10 ; Iteration step [time_go] = .01 ; Starting time-step (0... 0.25) [time_end] = .13 ; Ending time-step [time_step] = .06 ; Time-step step [fclevels] = 4 ; FC Contour Levels ENDIF [sigma_go] = 3.0 ; Starting sigma (Standard deviation) [sigma_end] = 3.0 ; Ending sigma [sigma_step] = 3.0 ; Sigma step [lambda_go] = 10 ; Starting lambda (Contrast-- integers) [lambda_end] = 10 ; Ending lambda [lambda_step] = 10 ; Lambda step [cutoff_go] = 1 ; Starting cutoff (Usually unused) [cutoff_end] = 1 ; Ending cutoff [cutoff_step] = 1 ; Cutoff step ; -------------- END BATCH HEADER ------------------------------------------- IF ([FR] .NE. 0) THEN ; Interative input parameters FR ; Input volume ? Input volume?[volin] FR ? Window location- X &Y ?[location] RR S [wx][wy] ; Split out variables [location] ; Variable 1 ; Positions FR ? Iterations- Initial, Final, & Step?[iterations] RR S [iter_go][iter_end],[iter_step]; Split out variables [iterations] ; Variable 1 ; Positions FR ? Times- Initial, Final, & Step?[times] RR S [time_go][time_end],[time_step]; Split out variables [times] ; Variable 1 ; Positions FR ? Contour levels ?[fclevels] FR ; Output volume ? Output stack?[stkout] FR ; Output volume ? Output montage?[montage] FR ; Output volume ? Output volume?[volout] VM echo " " ELSE VM echo " File: [volin]" ENDIF ; Remove any existing output files DE [stkout] ; Report input volume size, location, etc FI X [nsam],[nrow],[nslice] [volin] ; Input volume 12,2,1 ; nsam, nrow, & nslice locations IF ([slice] .EQ. 0) THEN [slice] = [nslice] / 2 ; Default is central slice ENDIF [wz] = [slice]-(INT([wnslice]/2)) ; Window upper left Z location VM echo " Size: {****[nsam]} x {****[nrow]} x {****[nslice]}" VM echo " Windowed: {****[wnsam]} x {****[wnrow]} x {***[wnslice]}" VM echo " At: {****[wx]}, {****[wy]}, {****[wz]}" ; Window the input volume WI ; Window operation [volin] ; Input _1 ; Output [wnsam],[wnrow],[wnslice] ; Window size [wx],[wy],[wz] ; Window location FS [winmax],[winmin] ; Window raw statistics _1 ; Input ; Scale image values over 0...255 AR SCA ; Scale image 0...255 _1 ; Input _2 ; Output 0 255 ; Scaling range ;NEG ; Negate the image (MAY NOT BE NECESSARY) ;_2 ;_2 ; Report windowed image range (for setting threshold) FS X [fmax],[fmin],[favg] _2 VM echo " Avg: {***[favg]}" VM echo " " [zsl] = [wnslice] / 2 [imgnum] = 0 ; Display Original Unfiltered Slice IF ([wnslice] .GT. 1) THEN ; Extract relevant slice from processed volume PS Z ; Pick slice _2 ; Input volume _3 ; Output slice [zsl] ; Slice number ELSE CP ; Copy _2 ; Input slice _3 ; Output slice ENDIF ; Increment image stack pointer [imgnum]=[imgnum]+1 ; Increment stack pointer LA ; Label the image _3 ; Input stack image [stkout]@{****[imgnum]} ; Output stack ORIGINAL IF ([PW] .GT. 0) THEN [padpx] = (4 * [wnsam]) + 1 ; Padded x dimension [padpy] = (4 * [wnrow]) + 1 ; Padded y dimension ; Power Spectrum of slice PD ; Pad into 4X larger image _3 ; Original slice input _4 ; Padded image [padpx],[padpy] ; Padded image size Y ; Use average for background PW ; Power Spectrum of slice _4 ; Padded image _5 ; Power spectrum output [px] = ([padpx] - [wnsam]) / 2 [py] = ([padpy] - [wnrow]) / 2 ; Window the PS WI ; Window operation _5 ; Input _4 ; Output [wnsam],[wnrow] ; Window size [px],[py] ; Window location [wcx] = ([wnsam]/2) + 1 [wcy] = ([wnrow]/2) + 1 ; Mask the PS MA ; Window operation _4 ; Input _5 ; Output (0,6) ; Outer & inner masking radii D ; Sharp disk E ; External background input 0 ; Background [wcx],[wcy] ; Mask center coordinates ; Increment image stack pointer [imgnum]=[imgnum]+1 ; Increment stack pointer LA ; Label the image _5 ; Input stack image [stkout]@{****[imgnum]} ; Output stack POWER SPEC ENDIF IF ([FC] .GT. 0) THEN ; Contour the output image FC ; Contour the image _3 ; Input _5 ; Output [fclevels] ; Levels Y ; Overwrite the image W ; White contours ; Increment image stack pointer [imgnum] = [imgnum] + 1 ; Put label in output image LA ; Label the image _5 ; Input image [stkout]@{****[imgnum]} ; Output stack CONTOURED ENDIF IF ([TH] .GT. 0) THEN ; Threshold the output image TH F ; Threshold the image _3 ; Input _5 ; Output B ; Below the cutoff set to zero [cutoff_go] ; Cutoff threshold (0.0) ; Value for Below the cutoff ; Increment image stack pointer [imgnum]=[imgnum]+1 ; Increment stack pointer LA ; Label the image _5 ; Input stack image [stkout]@{****[imgnum]} ; Output stack THRESHOLDED ENDIF ; ---------- process _2 Using: CE -- AD, HEG [iter_loops] = 1 ; Protect against division by zero IF ([iter_step] .GT. 0) THEN [iter_loops] = INT((([iter_end] - [iter_go]) / [iter_step]) + 0.5 ) + 1 ENDIF [time_loops] = 1 IF ([time_step] .GT. 0) THEN [time_loops] = INT((([time_end] - [time_go]) / [time_step]) + 0.5) + 1 ENDIF [sigma_loops] = 1 IF ([sigma_step] .GT. 0) THEN [sigma_loops] = INT((([sigma_end] - [sigma_go]) / [sigma_step])) + 1 ENDIF [lambda_loops] = 1 IF ([lambda_step] .GT. 0) THEN [lambda_loops] = INT((([lambda_end] - [lambda_go]) / [lambda_step])) + 1 ENDIF [cutoff_loops] = 1 IF ([cutoff_step] .GT. 0) THEN [cutoff_loops] = INT((([cutoff_end] - [cutoff_go]) / [cutoff_step])) + 1 ENDIF ; Find number of output images [frames] = [iter_loops]*[time_loops]*[sigma_loops]*[lambda_loops]*[cutoff_loops] ; Report steps VM echo " Iterations: {***[iter_loops]} Times: {***[time_loops]}" VM echo " Sigmas: {***[sigma_loops]} Lambdas: {***[lambda_loops]}" IF ([TH] .GT. 0) THEN VM echo " Cutoffs: {***[cutoff_loops]} Tests: {***[frames]}" ELSE VM echo " Tests: {***[frames]}" ENDIF VM echo " " VM echo " Output stack: [stkout]" VM echo " Output volume: [volout]" VM echo " " IF ([VOL] .GT. 0) THEN ; Want an output volume of denoised slice windows PS Z ; Pick slice [volin] ; Input volume _8 ; Output slice [slice] ; Slice number ; Scale the slice same as for denoised slice AR _8 ; Input slice _7 ; Output slice 255 * (P1-[winmin]) / ([winmax]-[winmin]) [padlocx] = [wx] + 0 ; Location for inset window [padlocy] = [wy] + 0 ; Location for inset window BL ; Blank volume [volout] [nsam],[nrow]+36,[frames]; Size adjusts for label N ; Use background 0 ; Background ENDIF ; Carry out filtering steps [vol_frame]=0 DO LB1 [iter_loop]=1,[iter_loops] [iter]= [iter_go] + ([iter_loop] -1) * [iter_step] DO LB2 [time_loop]=1,[time_loops] [time]= [time_go] + ([time_loop] -1) * [time_step] DO LB3 [sigma_loop]=1,[sigma_loops] [sigma]= [sigma_go] + ([sigma_loop] -1) * [sigma_step] DO LB4 [lambda_loop]=1,[lambda_loops] [lambda]= [lambda_go] + ([lambda_loop] -1) * [lambda_step] CE AD ; Anisotropic Diffusion _2 ; Input volume _3 ; Temp. output volume HEG ; AD filter type [iter] ; Iterations [time] ; Time step [sigma],[lambda] ; Sigma & Lambda IF ([wnslice] .GT. 1) THEN ; Extract relevant slice from processed volume PS Z ; Pick slice _3 ; Input volume _4 ; Output slice [zsl] ; Slice number ELSE CP ; Copy _3 ; Input slice _4 ; Output slice ENDIF IF ([TH] .LE. 0) THEN ; Display denoised image ; Increment image stack pointer [imgnum] = [imgnum] + 1 ; Echo image description to terminal VM ; Systems call echo " N: {***[imgnum]} Iter: {***[iter]} Time: {%f3.2%[time]} S: {%f3.1%[sigma]} L: {***[lambda]}" ; Put label in output image LA ; Label the image _4 ; Input image [stkout]@{****[imgnum]} ; Output stack I:{***[iter]} T:{%f3.2%[time]} ; Some alternatative Labels ; {***[imgnum]} ; I:{***[iter]} ; T:{%f3.2%[time]} ; S:{%f3.1%[sigma]} ; L:{**[lambda]}" ; C:{***[cutoff]} ELSE ; Display thresholded denoised image DO LB5 [cutoff_loop]=1,[cutoff_loops] [cutoff]= [cutoff_go] + ([cutoff_loop] -1) * [cutoff_step] ; Threshold the output image TH F ; Threshold the image _4 ; Input _5 ; Output B ; Below the cutoff set to zero [cutoff] ; Cutoff threshold (0.0) ; Value for Below the cutoff ; Increment image stack pointer [imgnum] = [imgnum] + 1 ; Echo image description to terminal VM ; Systems call echo " N: {***[imgnum]} TH Iter: {***[iter]} Time: {%f3.2%[time]} S: {%f3.1%[sigma]} L: {***[lambda]} C: {***[cutoff]}" ; Put label in output image LA ; Label the image _5 ; Input image [stkout]@{****[imgnum]} ; Output stack I:{***[iter]} T:{%f3.2%[time]} C: {***[cutoff]} LB5 ENDIF IF ([FC] .GT. 0) THEN ; Display contour on the denoised slice FC ; Contour the image _4 ; Input _5 ; Output 5 ; Levels Y ; Overwrite the image W ; White contours ; Increment image stack pointer [imgnum] = [imgnum] + 1 ; Echo image description to terminal VM ; Systems call echo " N: {***[imgnum]} FC Iter: {***[iter]} Time: {%f3.2%[time]} S: {%f3.1%[sigma]} L: {***[lambda]}" ; Put label in output image LA ; Label the image _5 ; Input image [stkout]@{****[imgnum]} ; Output stack I:{***[iter]} T:{%f3.2%[time]} ENDIF IF ([PW] .GT. 0) THEN ; Display Power Spectrum of denoised slice [px] = ([padpx] - [wnsam]) / 2 [py] = ([padpy] - [wnrow]) / 2 PD ; Pad into 4X larger image _4 ; Original slice input _5 ; Padded image [padpx],[padpy] ; Padded image size Y ; Use average for background [px],[py] ; Window location PW ; Power Spectrum of slice _5 ; Padded image _6 ; Power spectrum output ; Window the PS WI ; Window operation _6 ; Input _5 ; Output [wnsam],[wnrow] ; Window size [px],[py] ; Window location ; Mask the PS MA ; Window operation _5 ; Input _6 ; Output (0,6) ; Masking radii D ; Sharp disk E ; External background input 0 ; Background [wcx],[wcy] ; Mask center coordinates ; Increment image stack pointer [imgnum] = [imgnum] + 1 ; Echo image description to terminal VM ; Systems call echo " N: {***[imgnum]} PS Iter: {***[iter]} Time: {%f3.2%[time]} S: {%f3.1%[sigma]} L: {***[lambda]}" ; Put label in output image LA ; Label the image _6 ; Input image [stkout]@{****[imgnum]} ; Output stack I:{***[iter]} T:{%f3.2%[time]} ENDIF IF ([VOL] .GT. 0) THEN ; Create volume of slice with denoised areas CP _7 ; Original slice input _8 ; Padded image IN ; Pad into original slice image _4 ; Denoised slice input _8 ; Original intensity scaled slice [padlocx],[padlocy] ; Window location ; Put label in output image LA ; Label the image _8 ; Input image _9 ; Output image N: {***[imgnum]} I:{***[iter]} T:{%f3.2%[time]} S: {%f3.1%[sigma]} L: {***[lambda]} [vol_frame]=[vol_frame]+1 SK R ; Stack into labeled size volume [volout] ; Original intensity scaled slice _9 ; Denoised slice input [vol_frame] ; Slice location * ; Stop stacking ENDIF LB4 LB3 LB2 LB1 VM echo " Total Frames: {***[imgnum]}" DE [montage] MN S ; Montage the output images [stkout]@**** ; Input file template 1-[imgnum] ; File numbers 8,1 ; Images / row, margin width 0 ; Margin value [montage] ; Output file FS [montage] ; Output file VM echo " Output montage: [montage]" VM echo " " EN