{ Sobel gradient kernels (3x3,5x5,7x7,9x9) A 3x3 Sobel looks like this: 1 0 -1 1 2 1 D1 = 2 0 -2 D2 = 0 0 0 1 0 -1 -1 -2 -1 Magnitude is sqrt(D1^2 + D2^2) and the phase is atan(D2 / D1). Note: Due to a numerical overflow problem 9x9 is the largest Sobel filter computable via this method. A 11x11 kernel or larger will result in erroneous results. v1.0 - 17 April 1994 v1.1 - 21 April 1994 Changed to operate only on an ROI if selected. A timing procedure has been added to 1) estimate the amount of time necessary and 2) show the elapsed time in a status window. The status window shows the size of the ROI if one is selected. Calibration of estimated timing is performed internal to the macro. If your not using a Centris 650 use the Calibrate Machine Timing macro to produce a calibration factor. Then change the parameter mach in the procedure slowsobel. Perform a search for *Machine to find the correct line. Doug Morris University of Illinois Biomedical Magnetic Resonance Lab Urbana IL 61801 dmorris@bmrl.med.uiuc.edu Robert Rohland University of Minnesota Department of Pharmacology Mpls, MN 55455 bonzo@lenti.med.umn.edu } var ndim,tstatus,dtime,cal:integer; { global variables } procedure Timeit; var year, month,day,hour, minute,second,dayofweek:integer begin if tstatus=1 then begin GetTime(year,month,day,hour,minute,second,dayofweek); dtime := (hour*3600 + 60*minute + second) - dtime; SelectWindow('Sobel Status'); Writeln('Elapsed Time = ',dtime:10:0,' s'); if cal=1 then begin PutMessage('Calibration factor :',dtime/98:10:5); cal := 0; end; tstatus := 0; end; if tstatus=0 then begin GetTime(year,month,day,hour,minute,second,dayofweek); dtime := (hour*3600 + 60*minute + second); tstatus := 1; end; end; procedure hsobel; { Horizontal Filter } begin RequiresVersion(1.53); NewTextWindow('horizontal slow sobel',150,140); if ndim=3 then begin writeln('1 0 -1'); writeln('2 0 -2'); writeln('1 0 -1'); end; if ndim=5 then begin writeln('2 1 0 -1 -2'); writeln('3 2 0 -2 -3'); writeln('4 3 0 -3 -4'); writeln('3 2 0 -2 -3'); writeln('2 1 0 -1 -2'); end; if ndim=7 then begin writeln('3 2 1 0 -1 -2 -3'); writeln('4 3 2 0 -2 -3 -4'); writeln('5 4 3 0 -3 -4 -5'); writeln('6 5 4 0 -4 -5 -6'); writeln('5 4 3 0 -3 -4 -5'); writeln('4 3 2 0 -2 -3 -4'); writeln('3 2 1 0 -1 -2 -3'); end; if ndim=9 then begin writeln('4 3 2 1 0 -1 -2 -3 -4'); writeln('5 4 3 2 0 -2 -3 -4 -5'); writeln('6 5 4 3 0 -3 -4 -5 -6'); writeln('7 6 5 4 0 -4 -5 -6 -7'); writeln('8 7 6 5 0 -5 -6 -7 -8'); writeln('7 6 5 4 0 -4 -5 -6 -7'); writeln('6 5 4 3 0 -3 -4 -5 -6'); writeln('5 4 3 2 0 -2 -3 -4 -5'); writeln('4 3 2 1 0 -1 -2 -3 -4'); end; Convolve(''); Dispose; end; procedure vsobel; { Vertical sobel filter } begin RequiresVersion(1.53); NewTextWindow('vertical slow sobel',150,140); if ndim=3 then begin writeln(' 1 2 1'); writeln(' 0 0 0'); writeln('-1 -2 -1'); end; if ndim=5 then begin writeln('2 3 4 3 2'); writeln('1 2 3 2 1'); writeln('0 0 0 0 0'); writeln('-1 -2 -3 -2 -1'); writeln('-2 -3 -4 -3 -2'); end; if ndim=7 then begin writeln('3 4 5 6 5 4 3'); writeln('2 3 4 5 4 3 2'); writeln('1 2 3 4 3 2 1'); writeln('0 0 0 0 0 0 0'); writeln('-1 -2 -3 -4 -3 -2 -1'); writeln('-2 -3 -4 -5 -4 -3 -2'); writeln('-3 -4 -5 -6 -5 -4 -3'); end; if ndim=9 then begin writeln('4 5 6 7 8 7 6 5 4'); writeln('3 4 5 6 7 6 5 4 3'); writeln('2 3 4 5 6 5 4 3 2'); writeln('1 2 3 4 5 4 3 2 1'); writeln('0 0 0 0 0 0 0 0 0'); writeln('-1 -2 -3 -4 -5 -4 -3 -2 -1'); writeln('-2 -3 -4 -5 -6 -5 -4 -3 -2'); writeln('-3 -4 -5 -6 -7 -6 -5 -4 -3'); writeln('-4 -5 -6 -7 -8 -7 -6 -5 -4'); end; Convolve(''); Dispose; end; procedure slowsobel; var xsize,ysize,pid0,pid1,pid2,pid3,pid4,i,j,left,top,width,height:integer; n,mode,min,max,mean1,mean2,dev1,dev2,mag,phase,esttime:integer; mach:real; begin tstatus := 0; { *Machine factors for timing loop. Mac Centris 650, mach := 1.0 ; Mac IIci , mach := 2.044 ; Mac IIx , mach := 3.195 ; The procedure scales linearly with number of pixels. Local Calibration can be accomplished by executing the included macro " Calibrate Machine Timing". } mach := 1.0; ndim:=GetNumber('Filter Dimension? (3,5,7 or 9 )',3); if (((ndim<>3)AND(ndim<>5))AND((ndim<>7)AND(ndim<>9))) then begin PutMessage('Supported Kernel Dimensions are 3,5,7 and 9.'); Exit; end; SetCursor('watch'); GetRoi(left,top,width,height); if width <> 0 then begin Copy; SetNewSize(width,height); MakeNewWindow('ROI'); Paste; end; ScaleConvolutions(true); GetPicSize(xsize,ysize); esttime := trunc(mach*(3.9247 + 0.0057984*(xsize*ysize))); PutMessage(' Estimated time: ',esttime:6:0,' s'); pid0:=PicNumber; NewTextWindow('Sobel Status',200,100); Writeln('Estimated time: ',esttime:6:0,' s'); Timeit; SelectPic(pid0); SelectAll; Duplicate('D1'); hsobel; pid1:=PicNumber; SelectAll; Measure; GetResults(n,mean1,mode,min,max); SelectPic(pid0); SelectAll; Duplicate('D2'); vsobel; pid2:=PicNumber; SelectAll; Measure; GetResults(n,mean2,mode,min,max); SelectPic(pid0); Duplicate('Magnitude'); pid3:=PicNumber; SelectPic(pid0); Duplicate('Phase'); pid4:=PicNumber; SelectPic(pid3); for j:=0 to ysize do begin for i:=0 to xsize do begin ChoosePic(pid1); dev1:=GetPixel(i,j)-mean1; ChoosePic(pid2); dev2:=GetPixel(i,j)-mean2; mag:=sqrt(dev1*dev1+dev2*dev2); phase:=0; if dev1<>0 then phase:=128+80*arctan(dev2/dev1); ChoosePic(pid3); PutPixel(i,j,trunc(mag)); ChoosePic(pid4); PutPixel(i,j,trunc(phase)); end; end; SelectPic(pid2); Dispose; SelectPic(pid1); Dispose; if width <> 0 then begin SelectWindow('Sobel Status'); Writeln('ROI size: ',width:6:0,height:6:0); SelectWindow('ROI'); Dispose; end; Timeit; SetCursor('arrow'); end; macro 'Slow Sobel' begin slowsobel; end; procedure makeramp; var width,height,i,j,max,min,diff,mw,mh:integer; begin; height:= 128; width:=128; max := 127; min := 0; diff := max-min; mw := diff/width; mh := diff/height; SetNewSize(width,height); MakeNewWindow('Ramp'); SetCursor('watch'); for i:=0 to width do begin for j:=0 to height do begin Putpixel(i,j,i*mw+j*mh+min); end; end; SetCursor('arrow'); end; macro 'Calculate Machine Timing' var cal:integer; begin cal :=1; makeramp; SelectWindow('Ramp'); slowsobel; end;