FUNCTION MEAN, A B = TOTAL(A)/N_ELEMENTS(A) RETURN, B & END ;The executing program is xyoff_tex.pro. It performs a 2D ;cross-correlation, yielding the lag of maximum correlation, the ;coefficient of maximum correlation and the coefficient at lag 0. In ;practice, I usually used lag 0. Within the limits of floating point ;computation, the cross correlation at 0 shift of an image and itself is ;1.0 and the shifts are 0.0 and 0.0 as would be expected. xyoff_tex.pro ;calls max_pos.pro to find the maximum position within an array...the ;cross-correlation function in this case. rca_2D.pro is a setup and ;calling program for xyoff_tex.pro. You may not need it. ; ;If in running these, you find that they call other programs that are ;not included, let Tim know and he will supply you with them. Likewise, ;if you have any other problems. ;xyoff_tex.pro: REMOVE THIS LINE ;+ ; NAME: ; XYOFF ; PURPOSE: ; Calculates the correlation function between two images using a box of ; size (bx,by) surrounding a point in the image. By default this point ; is the center of the image, (nx/2,ny/2). The routine calculates the ; offsets by finding the maximum of the correlation fucntion. ; CATEGORY: ; IMAGE PROCESSING ; CALLING SEQUENCE: ; s = xyoff(a,b,bx,by,hx,hy,xx,yy,amax,fna,wdow=wdow,mask=mask,/sobel,ccf=ccf ; one=one,xy,amax) ; INPUTS: ; a 2D reference image. ; b 2D image that you want to track with respect to a ; bx size of the box in the x-dimension that you wish ; to use to generate the cross correlation function ; by size of the box in the y-dimension that you wish ; to use to generate the cross correlation function. ; hx x position in the image that corresponds to the center ; of the box. (Optional) ; hy y position in the image that corresponds to the center ; of the box. (Optional) ; fna fourier transform of the 256x256 (Optional) ; subtended image from the center of a. ; KEYWORD PARAMETERS: ; wdow=wdow 2-d apodization window. (Optional; default=1) ; mask=mask 2-d mask to be applied in fourier space. (Optional;default=1) ; /SOBEL = If specified, then SOBEL keyword causes this routine ; to track on the gradient of the image instead of the ; image itself. (It uses a Sobel edge enhancement operator) ; OUTPUTS: ; s = A vector containing the x and y offsets that should be applied ; to an image to maximize the cross correlation. ; s(0) = The x offset. ; s(1) = The y offset. ; s(2) = max value of the correlation ; s(3) = cross correlation at 0 shift ; s(4) = maximum of the real part of absolute value of cross correlation ; NOTES: ; The executing program is xyoff.pro. It performs a 2D ; cross-correlation, yielding the lag of maximum correlation, the ; coefficient of maximum correlation and the coefficient at lag 0. In ; practice, I usually used lag 0. Within the limits of floating point ; computation, the cross correlation at 0 shift of an image and itself is ; 1.0 and the shifts are 0.0 and 0.0 as would be expected. xyoff_tex.pro ; calls max_pos.pro to find the maximum position within an array...the ; cross-correlation function in this case. ; MODIFICATION HISTORY: ; H. Cohl, 25 Mar, 1992 --- Made documentation more understandable. ; H. Cohl, 13 Feb, 1992 --- Made it so that you can enter in hx or hy. ; H. Cohl, 13 Dec, 1991 --- Added SOBEL keyword. ; H. Cohl, 7 Jun, 1991 --- Initial programming. ; T. Henry December, 1993 --- some changes made to size ;- function xyoff, a,b,bx,by,hx,hy,xx,yy,amax,fna,wdow=wdow,mask=mask,sobel=sobel, $ ccf=ccf,one=one sobel = keyword_set(sobel) ssz = size(a) ; Size of first image if ssz(0) ne 2 then message, 'Main image must be 2D array nx = ssz(1) ny = ssz(2) ;message, /info, 'Main image is A('+ $ ; strcompress(nx,/rem)+','+strcompress(ny,/rem)+')' ssz = size(b) if ssz(0) ne 2 then message, 'Secondary image must be 2D image if nx ne ssz(1) or ny ne ssz(2) then begin message, /info, 'Secondary image is A('+ $ strcompress(ssz(1),/rem)+','+strcompress(ssz(2),/rem)+')' message, 'Both images must have same array dimensions endif if n_elements(hx) eq 0 then hx = nx/2 ; Check center of image if n_elements(hy) eq 0 then hy = ny/2 if n_elements(bx) eq 0 then bx = nx-1 if n_elements(by) eq 0 then by = ny-1 mx = bx/2 & my = by/2 ; Set correlation box size. if not keyword_set(wdow) then wdow=1. ; Test to see if window is specified. if not keyword_set(mask) then mask=1. xlo = hx-mx & xhi = hx+mx ylo = hy-my & yhi = hy+my ; In the original XYOFF program the values of xhi and yhi were one less ; than in the current version. Activate the next line to go back to the ; original version ; xhi = xhi-1 & yhi = yhi-1 ;message, /info, 'Center of correlation window: ('+ $ ; strcompress(hx,/rem)+','+strcompress(hy,/rem)+')' ;message, /info, 'Window used for correlation: ' + $ ; 'A('+strcompress(xlo,/rem)+':'+strcompress(xhi,/rem)+','+ $ ; strcompress(ylo,/rem)+':'+strcompress(yhi,/rem)+')' if xlo lt 0 or xhi gt nx-1 or $ ylo lt 0 or yhi gt ny-1 then message, 'Invalid subarray index range na = a(xlo:xhi,ylo:yhi) if min(na) eq max(na) then message, 'First image is flat na = na-total(na)/n_elements(na) ; Subtract average if sobel then na = 1.*sobel(na) else na = na*wdow fna = fft(na*mask,-1) ; Fourier transform rms_na = sqrt(total(na*na)/n_elements(na)) nb = b(xlo:xhi,ylo:yhi) if min(na) eq max(na) then message, 'Second image is flat nb = nb-total(nb)/n_elements(nb) if sobel then nb = 1.*sobel(nb) else nb = nb*wdow fnb = fft(nb*mask,-1) ; Fourier transform rms_nb = sqrt(total(nb*nb)/n_elements(nb)) ccf = fft(fna*conj(fnb),1) ; Correlation fnc between na and nb ccf = float(ccf) ; Get real part of cross correlation fnc ccf = ccf/(rms_na*rms_nb) ; Normalize ccf = shift(ccf,mx,my) ; Shift correlation function. ps = max_pos(ccf) ; Find maximum of ccf s = [ps(0)-mx, ps(1)-my, max(ccf), ccf(bx/2,by/2), max(abs(ccf)) ] ;help, ccf ; xx = s(0) yy = s(1) amax = s(2); ; printf,one, xy,amax ; ;print, 'RMS IMAGE A & B : ', rms_na, rms_nb ;print, 'Position of Maximum(CCF)= ', [mx,my]+s(0:1) ;print, 'Shift in X - ', s(0) ;print, 'Shift in Y - ', s(1) ;print, 'Maximum(CCF) - ', s(2) ;print, 'Center CCF - ', s(3) ; ss=fltarr(9) ; ss(0) = ccf(bx/2-4,by/2) ; ss(1) = ccf(bx/2-3,by/2) ; ss(2) = ccf(bx/2-2,by/2) ; ss(3) = ccf(bx/2-1,by/2) ; ss(4) = ccf(bx/2,by/2) ; ss(5) = ccf(bx/2+1,by/2) ; ss(6) = ccf(bx/2+2,by/2) ; ss(7) = ccf(bx/2+3,by/2) ; ss(8) = ccf(bx/2+4,by/2) ;print,'x cor values ',ss(0),ss(1),ss(2),ss(3),ss(4),ss(5),ss(6),ss(7),ss(8) ;print, 'Maximum( ABS(CCF) )= ', s(4) ;print, 'Position of Maximum( ABS(CCF) )= ',max_pos(abs(ccf)) if n_elements(one) ne 0 then begin ;printf,one, 'RMS IMAGE A & B : ', rms_na, rms_nb ;printf,one, 'Position of Maximum(CCF)= ', [mx,my]+s(0:1) ;printf,one, 'Shift in X - ', s(0) ;printf,one, 'Shift in Y - ', s(1) ;printf,one, 'Maximum(CCF) - ', s(2) ;printf,one, 'Center CCF - ', s(3) ;printf,one,'x cor values ',ss(0),ss(1),ss(2),ss(3),ss(4),ss(5),ss(6),ss(7),ss(8) ;printf,one, 'Maximum( ABS(CCF) )= ', s(4) ;printf,one, 'Position of Maximum( ABS(CCF) )= ',max_pos(abs(ccf)) endif return, s & end ;max_pos.pro: REMOVE THIS LINE ;------------------------------------------------------------- ;+ ; NAME: ; MAX_POS ; PURPOSE: ; This function, given an array with a single well ; defined maximum, finds the maximum position ; of that array to decimal accuracy using an interpolation ; technique to find the maximum of an array from, Niblack, W., ; "An Introduction to Digital Image Processing", p 139. ; CATEGORY: ; IMAGE PROCESSING ; CALLING SEQUENCE: ; max = max_pos(f) ; INPUTS: ; f = An array. ; type: array,any type ; KEYWORD PARAMETERS: ; OUTPUTS: ; max = An array containing the position of the maximum of an array. ; type: vector,floating point,fltarr(2) ; max(0) = The decimal position of the x-maximum. ; type: scalar,floating point ; max(1) = The decimal position of the y-maximum. ; type: scalar,floating point ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; H. Cohl, 20 Nov, 1990 --- Initial implementation. ;- ;------------------------------------------------------------- function max_pos,f,help=help ;Display idl header if help is required. if keyword_set(help) or n_params() lt 1 then begin get_idlhdr,'max_pos.pro' max=-1 goto,finishup endif ssz=size(f) fmax=max(f,loc) xmax=loc mod ssz(1) ymax=loc/ssz(1) ;A more complicated interpolation technique to find the maximum of an array. ;from, Niblack, W., "An Introduction to Digital Image Processing", p 139. if (xmax*ymax gt 0) and (xmax lt (ssz(1)-1)) and (ymax lt (ssz(2)-1)) then begin denom=fmax*2-f(xmax-1,ymax)-f(xmax+1,ymax) xfra=(xmax-.5)+(fmax-f(xmax-1,ymax))/denom denom=fmax*2-f(xmax,ymax-1)-f(xmax,ymax+1) yfra=(ymax-.5)+(fmax-f(xmax,ymax-1))/denom xmax=xfra ymax=yfra endif rmax=[xmax,ymax] finishup: return,rmax end ; function bound,ix,iy,xc,yc,scx,scy,az,ht ; ; Determine the azimuth and height from the input ix and iy ; radtd = 57.29578 ; scx2 = scx*scx scy2 = scy*scy az = atan((iy - yc),(ix - xc))*radtd ht = sqrt((iy - yc)*(iy - yc)*scy2 + (ix - xc)*(ix - xc)*scx2) if az lt 0.0 then az = 360.0 + az ; END ; function bestv,sp,amp,ang,spn,ampn,angn,nx,ny,xc,yc,scx,scy, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,bestsp0,bestsp1,bestsp2, $ bestamp0,bestamp1,bestamp2,bestang0,bestang1,bestang2 minxx = 2 minyy = 2 ivv = miny ivp = maxy ihh = minx ihp = maxx if ivv lt 0 then ivv = 0 if ivp ge ny then ivp = ny - 1 if ihh lt 0 then ihh = 0 if ihp ge nx then ihp = nx - 1 bx0 = -1 by0 = -1 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin if amp(ixx,iyy) gt ampn(ixx,iyy) then begin ampn(ixx,iyy) = amp(ixx,iyy) spn(ixx,iyy) = sp(ixx,iyy) angn(ixx,iyy) = ang(ixx,iyy) endif if amp(ixx,iyy) gt bestamp0 then begin bestamp0 = amp(ixx,iyy) bestsp0 = sp(ixx,iyy) bestang0 = ang(ixx,iyy) bx0 = ixx by0 = iyy endif endif endif ENDFOR ENDFOR if bx0 eq -1 then goto, endit bx1 = -1 by1 = -1 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin xdif = abs(ixx - bx0) ydif = abs(iyy - by0) if amp(ixx,iyy) gt bestamp1 $ and xdif gt minxx and ydif gt minyy then begin bestamp1 = amp(ixx,iyy) bestsp1 = sp(ixx,iyy) bestang1 = ang(ixx,iyy) bx1 = ixx by1 = iyy endif endif endif ENDFOR ENDFOR if bx1 eq -1 then goto, endit bx2 = -1 by2 = -1 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin if bx1 ne -1 and by1 ne -1 then begin xdif = abs(ixx - bx0) ydif = abs(iyy - by0) xdif1 = abs(ixx - bx1) ydif1 = abs(iyy - by1) xdif2 = abs(bx0 - bx1) ydif2 = abs(by0 - by1) if amp(ixx,iyy) gt bestamp2 $ and xdif gt minxx and ydif gt minyy $ and xdif1 gt minxx and ydif1 gt minyy $ and xdif2 gt minxx and ydif2 gt minyy then begin bestamp2 = amp(ixx,iyy) bestsp2 = sp(ixx,iyy) bestang2 = ang(ixx,iyy) bx2 = ixx by2 = iyy endif endif endif endif ENDFOR ENDFOR endit: ; end ; function cimage,i1,i2,nh1,nv1,nh2,nv2,xc1,yc1,xc2,yc2, $ scx1,scy1,scx2,scy2,pspeed,dtime,az1,az2,ht1,ht2, $ minx,maxx,miny,maxy, $ pspeedxy,xcent,ycent,kk,itrigger,imgg, $ az11,az22,ht11,ht22,xsize,ysize,one, $ angl,sspeed,fspeed,smax,badspeed,badangle,badcorr, $ cspeed,cangl,cmax ; ; This function obtains sets of arrays of the proper size and displacement ; to be used in a correlation ; radtd = 57.29578 ; istep = 1 ; step increment between correlated data sets in x jstep = 1 ; step increment between correlated data sets in y xsize2p = 2*xsize+1 ; 2 times the x size of the array to correlate ysize2p = 2*ysize+1 ; 2 times the y size of the array to correlate jsizeh = ysize/2 ; half of maximum y value of correlation isizeh = xsize/2 ; half of maximum x value of correlation scxsq = scx1*scx1 scysq = scy1*scy1 ; azb1 = az1 htb1 = ht1 azb2 = az2 htb2 = ht2 azb11 = az11 htb11 = ht11 azb22 = az22 htb22 = ht22 ; ; Correlation arrays ; i11 = fltarr(xsize2p,ysize2p) i22 = fltarr(xsize2p,ysize2p) ; print, ' Enter correlation' ; ivvv = ysize if miny gt ivvv then ivvv = miny ; ivvp = nv1m - ysize if maxy lt ivvp then ivvp = maxy ; ihhh = xsize if minx gt ihhh then ihhh = minx ; ihhp = nh1m - xsize if maxx lt ihhp then ihhp = maxx nxxx = 0L nxxxx = 0L kk = 0L FOR iy=ivvv,ivvp,jstep do BEGIN ixxx = 0L FOR ix=ihhh,ihhp,istep do BEGIN ; ; Determine the boundary within which to do the correlation on image 1 ; bd1 = bound(ix,iy,xc1,yc1,scx1,scy1,azb,htb) if azb ge azb2 and azb le azb1 then begin if htb ge htb1 and htb le htb2 then begin ; ivv = iy-ysize ivp = iy+ysize ihh = ix-xsize ihp = ix+xsize ; print, '1st ivv ivp ihh ihp', ivv,ivp,ihh,ihp ; ; Load up the first array to use in the correlation. ; iyv = 0 FOR iyy=ivv,ivp do BEGIN ixh = 0 FOR ixx=ihh,ihp do BEGIN i11(ixh,iyv) = i1(ixx,iyy) ixh = ixh + 1 ENDFOR iyv = iyv + 1 ENDFOR ; iflag1 = 0 iflag2 = 0 iflag3 = 0 iflag4 = 1 ; ; Determine a new x, y position of the center of the second window to ; correlate by assuming radial outward motion at the preferred speed. ; phid = atan((iy - yc1)*scx1,(ix - xc1)*scy1) ; ; Preferred x displacement to next image ; ixd = (pspeed*dtime)/scx1*cos(phid) ; ; Preferred y displacement to next image ; iyd = (pspeed*dtime)/scy1*sin(phid) ; inx = ix + ixd inx = fix(inx) iny = iy + iyd iny = fix(iny) ; ; Determine if the boundary is within the boundary correlation of image 2 ; bd2 = bound(inx,iny,xc2,yc2,scx2,scy2,azb,htb) if azb ge azb22 and azb le azb11 then begin if htb ge htb11 and htb le htb22 then begin ; ; Load up the second array to use in the correlation. ; ivv = iny-ysize ivp = iny+ysize ihh = inx-xsize ihp = inx+xsize ; print, '2nd ivv ivp ihh ihp', ivv,ivp,ihh,ihp iyv = 0 FOR iyy=ivv,ivp do BEGIN ixh = 0 FOR ixx=ihh,ihp do BEGIN if iyy lt 0 or iyy ge nv2 or ixx lt 0 or ixx ge nh2 $ then iflag4 = 0 else BEGIN i22(ixh,iyv) = i2(ixx,iyy) ; if i11(ixh,iyv) ne i22(ixh,iyv) then iflag1 = 1 if i11(ixh,iyv) ne i11(0,0) then iflag2 = 1 if i22(ixh,iyv) ne i22(0,0) then iflag3 = 1 ; ixh = ixh + 1 endelse ENDFOR iyv = iyv + 1 ENDFOR ; endif endif ; ; The load-up of arrays to use in the correlation is now finished. ; Now, correlate the arrays as long as there is data in the arrays and ; as all the data in one of the arrays is not flat (without structure). ; if iflag1 eq 1 then begin if iflag2 eq 1 then begin if iflag3 eq 1 then begin if iflag4 eq 1 then begin ; i1Xi2 = xyoff(i11,i22,xsize,ysize,xsize,ysize,xxi,yyj, $ amax,ccf=ccf,one=one) ; nxxx = nxxx + 1L ; ; print, 'xxi yyj amax', xxi,yyj,amax if abs(xxi) gt xsize/3 or abs(yyj) gt xsize/3 then begin xys = badspeed dang = badangle amax = badcorr endif else begin xx = xxi - ixd yy = yyj - iyd phi = atan((iy - yc1)*scy1,(ix - xc1)*scx1)*radtd ipflag = 0 if phi lt 0.0 then begin phi = phi + 360.0 ipflag = 1 endif the = atan(-yy,-xx)*radtd if the lt 0.0 then begin the = the + 360.0 if ipflag eq 0 then phi = phi + 360.0 endif dang = (the - phi) xy = sqrt(xx*xx*scxsq + yy*yy*scysq) xys = xy/dtime ; mini1 = min(i11(*,*)) ; maxi1 = max(i11(*,*)) ; print, 'mini1 maxi1 ix iy', mini1, maxi1, ix, iy ; mini2 = min(i22(*,*)) ; maxi2 = max(i22(*,*)) ; print,' mini2 maxi2 inx iny', mini2, maxi2, inx, iny ; print, 'xys dang amax',xys,dang,amax endelse ; print, yy,xx,xy,xys,iy,yc1,ix,xc1,phi*radtd,the*radtd ;if amax gt .9 then print, yy,xx,xy,xys,iy,yc,ix,xc,phi*radtd,the*radtd,'angle diff =',dang if abs(dang) gt angl then begin xys = badspeed dang = badangle amax = badcorr endif if amax lt smax then begin xys = badspeed dang = badangle amax = badcorr endif if xys lt sspeed then begin xys = badspeed dang = badangle amax = badcorr endif if xys gt fspeed then begin xys = badspeed dang = badangle amax = badcorr endif endif endif endif endif if iflag1 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag2 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag3 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag4 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif ; if xys ne badspeed and dang ne badangle and amax ne badcorr then begin nxxxx = nxxxx + 1L ; *** if imgg eq 1 then begin if amax gt .20 and amax lt .25 then begin if iy gt (yc1 + 30) and ix gt xc1 then begin if itrigger eq 0 and xys gt 150 and xys lt 200 and dang gt -2.0 and $ dang lt 2.0 then begin print, 'set of images correlated and correlation', ixd, xx, xxi, iyd, $ yy, yyj, xy, xys, phi, the, $ dang, amax printf, one, 'set of images correlated and correlation', ixd, xx, xxi, iyd, $ yy, yyj, xy, xys, phi, the, $ dang, amax ; print, 'Contouring first image' TWIN, /new, xsize = 200, ysize=200 mini1 = min(i11(*,*)) maxi1 = max(i11(*,*)) print, '1st image mini1 maxi1 ix iy', mini1, maxi1, ix, iy printf, one, '1st image mini1 maxi1 ix iy', mini1, maxi1, ix, iy i11(*,*) = i11(*,*) - mini1 i11(*,*) = I11(*,*)*(0.5)/(maxi1-mini1) contour, i11, levels = [0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35, $ 0.4,0.45,0.5], $ c_label = [1,1,1,1,1,1,1,1,1,1,1,1], $ c_colors = [81,84,88,92,96,100,104,108,112,116,120], /fill, $ xstyle=1, ystyle=1, pos=[0.1,0.1,0.9,0.9] ; print, 'Contouring second image' TWIN, /new, xsize = 200, ysize=200 mini2 = min(i22(*,*)) maxi2 = max(i22(*,*)) print,'2nd image mini2 maxi2 inx iny', mini2, maxi2, inx, iny printf, one, '2nd image mini2 maxi2 inx iny', mini2, maxi2, inx, iny i22(*,*) = i22(*,*) - mini2 i22(*,*) = I22(*,*)*(0.5)/(maxi2-mini2) contour, i22, levels = [0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35, $ 0.4,0.45,0.5], $ c_label = [1,1,1,1,1,1,1,1,1,1,1,1], $ c_colors = [81,84,88,92,96,100,104,108,112,116,120], /fill, $ xstyle=1, ystyle=1, pos=[0.1,0.1,0.9,0.9] ; print, 'Contouring correlation height for a sample' TWIN, /new, xsize = 200, ysize=200 minic = min(ccf(*,*)) maxic = max(ccf(*,*)) print,'Correlation minic maxic inx iny', minic, maxic, xxi, yyj printf, one, 'Correlation minic maxic inx iny', minic, maxic, xxi, yyj cc = fltarr(xsize,ysize) cc(*,*) = ccf(*,*) - minic cc(*,*) = cc(*,*)*(0.5)/(maxic-minic) contour, cc, levels = [0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35, $ 0.4,0.45,0.5], $ c_label = [1,1,1,1,1,1,1,1,1,1,1,1], $ c_colors = [121,124,128,132,136,140,144,148,152,156,160], /fill, $ xstyle=1, ystyle=1, pos=[0.1,0.1,0.9,0.9] ; itrigger = 1 endif endif endif endif ; ; *** ; endif cspeed(ix,iy) = xys cangl(ix,iy) = dang cmax(ix,iy) = amax ixxx = ixxx + 1L ; ; if iflag1 eq 1 and iflag2 eq 1 and iflag3 eq 1 then begin if xys ne badspeed and dang ne badangle and amax ne badcorr then begin pspeedxy(kk) = cspeed(ix,iy) xcent(kk) = ix ycent(kk) = iy kk = kk + 1L endif ; endif endif ENDFOR ENDFOR ; if nxxx gt 0L then begin print, ' success ', nxxx, ' pixels ',xsize,' by ',ysize,' were correlated' printf,one,' success ', nxxx, ' pixels ',xsize,' by ',ysize,' were correlated' if nxxxx gt 0L then begin print, ' more success ', nxxxx, ' pixels past all criterion' printf, one, ' more success ', nxxxx, ' pixels past all criterion' endif endif if nxxx eq 0L then begin print, ' failure - no pixels correlated printf, one,' failure - no pixels correlated endif END ; ; function findmx, az1,az2,ht1,ht2,xc,yc,scx,scy,xsize,ysize,iave, $ minx,maxx,miny,maxy ; ; This function finds the minimum and maximum x and y values from the values ; of azimuth and height. ; radtd = 57.29578 ; naz = abs(az1-az2) nazp = naz+1 nazp2 = nazp*2 htx = fltarr(nazp2) hty = fltarr(nazp2) htx(*) = -999 hty(*) = -999 ; for i=0L,nazp do begin si = sin((az1-i)/radtd) co = cos((az1-i)/radtd) htx(i) = xc + (ht1/scx)*co hty(i) = yc + (ht1/scy)*si htx(nazp2-1-i) = xc + (ht2/scx)*co hty(nazp2-1-i) = yc + (ht2/scy)*si endfor iav = xsize if ysize gt iav then iav = ysize if iave gt iav then iav = iave minx = min(htx) - iav minx = fix(minx) maxx = max(htx) + iav + 1.0 maxx = fix(maxx) miny = min(hty) - iav miny = fix(miny) maxy = max(hty) + iav + 1.0 maxy = fix(maxy) ; END ; function filtr,ar,nh,nv,iave,xsize,ysize,xc,yc,scx,scy, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,i1 ; ; This function filters an image file a within the boundaries of the windowed ; area ; radtd = 57.29578 ; i1 = fltarr(nh,nv) ; FOR iyy=0L,nv-1 do BEGIN FOR ixx=0L,nh-1 do BEGIN i1(ixx,iyy) = -5.0 ENDFOR ENDFOR ; nxxx = 0L ; iavx = iave if xsize gt iave then iavx = xsize iavy = iave if ysize gt iave then iavy = ysize iavs = iavx*scx if iavy*scy gt iavs then iavs = iavy*scy ; dazbd = abs(atan(iavs,ht1))*radtd ; azb1 = az1 + dazbd htb1 = ht1 - iavs azb2 = az2 - dazbd htb2 = ht2 + iavs ; print, 'iavs azb1 azb2 htb1 htb2', iavs,azb1,azb2,htb1,htb2 ; ivv = iave-1 if miny gt ivv then ivv = miny ; ivp = nv - iave if maxy lt ivp then ivp = maxy ; ihh = iave-1 if minx gt ihh then ihh = minx ; ihp = nh - iave if maxx lt ihp then ihp = maxx ; print, ' Enter filtering' FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge azb2 and azb le azb1 then begin if htb ge htb1 and htb le htb2 then begin nxxx = nxxx + 1L xyave = ar(ixx,iyy) xxyy = sqrt((iyy -yc)*(iyy -yc) + (ixx -xc)*(ixx -xc)) if xxyy ne 0.0 then BEGIN ax = (ixx - xc)/xxyy endif else BEGIN ax = 0.0 endelse if xxyy ne 0.0 then BEGIN ay = (iyy - yc)/xxyy endif else BEGIN ay = 0.0 endelse FOR i=0L,iave-1 do BEGIN iix = i*ax iiy = i*ay ixm = ixx - iix iym = iyy - iiy ixp = ixx + iix iyp = iyy + iiy xyave = xyave + ar(ixm,iym) + ar(ixp,iyp) ENDFOR xyave = xyave/(2.0*iave-1) i1(ixx,iyy) = ar(ixx,iyy) - xyave endif endif ENDFOR ENDFOR print, ' success - filtering ',nxxx,' points' end ; function mxmnbv, a,az1,az2,ht1,ht2,xc,yc,scx,scy, $ minx,maxx,miny,maxy,mxbv,mnbv ; ; This function finds the minimum and maximum image values within the area ; bound by azimuth and height. ; radtd = 57.29578 ; mxbv = -100000.0 mnbv = 100000.0 for i=minx,maxx do begin for j=miny,maxy do begin bd = bound(i,j,xc,yc,scx,scy,azb,htb) if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin if a(i,j) gt mxbv then mxbv = a(i,j) if a(i,j) lt mxbv then mnbv = a(i,j) endif endif endfor endfor END ; function smoo,ar,nh,nv,xsize,ysize,xc,yc,scx,scy, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,i1 ; ; This function smooths an image file within the boundaries of the windowed ; area so that point glitches, cosmic rays, especially are removed from ; the images. ; radtd = 57.29578 ; i2 = fltarr(nh,nv) i1 = fltarr(nh,nv) ; i1(*,*) = -5.0 i2(*,*) = ar(*,*) ; difmax = 800. ; iavx = xsize iavy = ysize ; iavs = iavx*scx if iavy*scy gt iavs then iavs = iavy*scy ; dazbd = abs(atan(iavs,ht1))*radtd ; azb1 = az1 + dazbd htb1 = ht1 - iavs azb2 = az2 - dazbd htb2 = ht2 + iavs ; ivv = iavy-1 if miny gt ivv then ivv = miny ; ivp = nv - iavy if maxy lt ivp then ivp = maxy ; ihh = iavx-1 if minx gt ihh then ihh = minx ; ihp = nh - iavx if maxx lt ihp then ihp = maxx ; print, ' Enter Smoothing' ; ; First pass ; nx1 = 0L FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) if azb ge azb2 and azb le azb1 then begin if htb ge htb1 and htb le htb2 then begin arsur = (ar(ixx+1)*2.0 + ar(ixx-1)*2.0 + ar(iyy+1)*2.0 $ + ar(iyy-1)*2.0 + ar(ixx+1,iyy+1) + ar(ixx-1,iyy-1) $ + ar(ixx+1,iyy-1) + ar(ixx-1,iyy+1))/12.0 ardiff = abs(ar(ixx,iyy) - arsur) if ardiff gt difmax then begin i2(ixx,iyy) = arsur nx1 = nx1 + 1L endif endif endif ENDFOR ENDFOR if nx1 gt 0 then print, 'First pass found', nx1, ' points to smooth' ar(*,*) = i2(*,*) ; ; Second pass ; nx2 = 0L FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) if azb gt azb2 and azb lt azb1 then begin if htb gt htb1 and htb lt htb2 then begin arsur = (ar(ixx+1)*2.0 + ar(ixx-1)*2.0 + ar(iyy+1)*2.0 $ + ar(iyy-1)*2.0 + ar(ixx+1,iyy+1) + ar(ixx-1,iyy-1) $ + ar(ixx+1,iyy-1) + ar(ixx-1,iyy+1))/12.0 ardiff = abs(ar(ixx,iyy) - arsur) if ardiff gt difmax then begin i2(ixx,iyy) = arsur nx2 = nx2 + 1L endif endif endif ENDFOR ENDFOR if nx2 gt 0 then print, 'Second pass found', nx2, ' points to smooth' ; i1(*,*) = i2(*,*) ; end ; function getcp,xc,yc,scx,scy,az1,az2,ht1,ht2,iok ; ; Determine an image area to analyze in the study. ; radtd = 57.29578 ; if iok eq 0 then begin ; pcursor, xp11, yp11, /data pcursor, xp22, yp22, /data ; xp1 = xp11 yp1 = yp11 xp2 = xp22 yp2 = yp22 scx2 = scx*scx scy2 = scy*scy az11 = atan((yp1 - yc),(xp1 - xc))*radtd if az11 lt 0.0 then az11 = az11 + 360.0 az22 = atan((yp2 - yc),(xp2 - xc))*radtd if az22 lt 0.0 then az22 = az22 + 360.0 print, az11, az22 yp1m = yp1 - yc yp2m = yp2 - yc if az11 lt az22 then az11 = 360.0 + az11 ht11 = sqrt((yp1 - yc)*(yp1 - yc)*scy2 + (xp1 - xc)*(xp1 - xc)*scx2) ht22 = sqrt((yp2 - yc)*(yp2 - yc)*scy2 + (xp2 - xc)*(xp2 - xc)*scx2) az1 = az11 az2 = az22 ht1 = ht11 ht2 = ht22 if ht1 gt ht2 then begin ht2 = ht11 ht1 = ht22 endif endif ; END ; function getok, ih, iv, okok ; ; If the cursor is clicked within the dot on the overplot box of imagea, then ; a "true" is returned for the function. ; getok = 1 ; the only time getok is false is if ; you say so with the cursor if okok eq 0 then begin pcursor, xp, yp, /data getok = xp gt ih and xp lt ih+20 and yp gt 0 and yp lt 20 if getok eq 1 then print, ' Keep this one!' if getok eq 0 then print, ' Try again!' endif return, getok END ; function getspeed,xc,yc,scx,scy,ht11,ht22,dt,psp,iok ; ; Determine an image area to analyze in the study. ; radtd = 57.29578 ; if iok eq 0 then begin ; ; Ht from cursor mark ; print, 'define speed - click on center of new box - speed is now',psp pcursor, xp11, yp11, /data ; xp1 = xp11 yp1 = yp11 scx2 = scx*scx scy2 = scy*scy ht2 = sqrt((yp1 - yc)*(yp1 - yc)*scy2 + (xp1 - xc)*(xp1 - xc)*scx2) ; ; Ht of earlier box center ; ht1 = (ht11 + ht22)/2.0 ; ; Preferred speed ; psp = (ht2 - ht1)/dt ; endif ; END ; function image,name,ih,iv,mxbv,mnbv,imgg,a ; ; Image reader ; this function reads in a LASCO difference image and plots it ; a = readfits(name,hdr) ; mxdif = mxbv - mnbv mxd3 = mxdif/3.0 mx = mxbv + mxd3 mn = mnbv - mxd3 print, 'mxbv mnbv mx mn', mxbv, mnbv, mx, mn if imgg eq 1 then begin TWIN, /new, xsize=ih+20, ysize=iv MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Reddish(0) + bytscl( a, max = mx, $ min = mn, top = Reddish(1)-Reddish(0)) ; print, reddish(0),reddish(1),greenish(0),greenish(1),blueish(0),blueish(1) TV, a0 endif ; END ; function image3,spv,ih,iv,sspeed,fspeed,minc,maxc ; ; This function images the speed two-dimensional image inside the boundary. ; All other values are set to the slowest speed. ; TWIN, /new, xsize=ih, ysize=iv a0 = minc + bytscl( spv, max = fspeed, $ min = sspeed, top = maxc-minc) TV, a0 ; END ; function image4,ang,ih,iv,angl,sspeed,fspeed,smax,az1,az2,ht1,ht2 ; ; This function images the angle two-dimensional image ; TWIN, /new, xsize=ih, ysize=iv MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Greenish(0) + bytscl( a, max = angl, $ min = -angl, top = Greenish(1)-Greenish(0)) TV, a0 ; END ; function image5,sma,ih,iv,angl,sspeed,fspeed,smax,az1,az2,ht1,ht2 ; ; This function images the correlation maximum two-dimensional image ; TWIN, /new, xsize=ih, ysize=iv MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Reddish(0) + bytscl( a, max = 1.0, $ min = 0.0, top = Reddish(1)-Reddish(0)) TV, a0 ; END ; function image9,name,spv,ang,sma,ih,iv,badspv,badang,badsma,one,imgg, $ xsiz,ysiz,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,mxbv,mnbv ; ; This function reads in a LASCO difference image, plots it and overplots ; arrows which map the outward speed of features within the window. ; if imgg eq 1 then begin radtd = 57.29578 ; xsz = xsiz/3 ; nearest x corr. step to arrow base ysz = ysiz/3 ; nearest y corr. step to arrow base ihm = ih - 1 ivm = iv - 1 ; a = readfits(name,hdr) ; twin, /new, xsize = ih, ysize = iv mxdif = mxbv - mnbv mxd3 = mxdif/3.0 mx = mxbv + mxd3 mn = mnbv - mxd3 print, 'mxbv mnbv mx mn', mxbv, mnbv, mx, mn TWIN, /new, xsize=ih+20, ysize=iv MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Reddish(0) + bytscl( a, max = mx, $ min = mn, top = Reddish(1)-Reddish(0)) TV, a0 ; pix = 10.0/100.0 ; arrow length scale (pixels per km/s) ahead = 5.0 ; length of arrowhead dphi = 25.0/radtd ; angular shift for arrowhead spvs = 0.0 angs = 0.0 smas = 0.0 nsums = 0.0 arx = fltarr(2) ary = fltarr(2) arhx = fltarr(3) arhy = fltarr(3) FOR i=0L, ihm, xsz DO BEGIN FOR j=0L, ivm, ysz DO BEGIN if spv(i,j) ne badspv and ang(i,j) ne badang and sma(i,j) ne badsma $ then begin print, ' success', i,j,spv(i,j),ang(i,j),sma(i,j) printf, one,' success', i,j,spv(i,j),ang(i,j),sma(i,j) spvs = spvs + spv(i,j) angs = angs + ang(i,j) smas = smas + sma(i,j) nsums = nsums + 1.0 phid = atan((j - vcent1)*scy1,(i - hcent1)*scx1)*radtd if phid lt 0.0 then phid = phid + 360.0 phi = (phid + ang(i,j))/radtd ; phi in radians xd = spv(i,j)*pix*cos(phi) yd = spv(i,j)*pix*sin(phi) xp = i + xd yp = j + yd xpm = xp - ahead*cos(phi-dphi) ypm = yp - ahead*sin(phi-dphi) xpp = xp - ahead*cos(phi+dphi) ypp = yp - ahead*sin(phi+dphi) ; ; set up plotting arrays ; arx(0) = i ary(0) = j arx(1) = xp ary(1) = yp arhx(0) = xpm arhy(0) = ypm arhx(1) = xp arhy(1) = yp arhx(2) = xpp arhy(2) = ypp ; oplot, arx, ary, thick=2, /noclip oplot, arhx, arhy, thick=2, /noclip endif ENDFOR ENDFOR if nsums ne 0.0 then begin spvs = spvs/nsums angs = angs/nsums smas = smas/nsums print, ' summary of successes ', nsums,spvs,angs,smas printf, one, ' summary of successes ', nsums,spvs,angs,smas print, ' ' print, ' ' printf, one, ' ' printf, one, ' ' endif endif ; END ; function image10,name,ih,iv,az1,az2,ht1,ht2,xc,yc,scx,scy,imgg, $ minx,maxx,miny,maxy,value,minva,maxva,mxbv,mnbv,minc,maxc ; ; Image reader ; ; This function reads in a LASCO difference image and plots it outside the ; boundary defined before. Inside the boundary is plotted the appropriate ; parameter - pixel by pixel. ; if imgg eq 1 then begin a = readfits(name,hdr) ; azb1 = az1 htb1 = ht1 azb2 = az2 htb2 = ht2 ; ivv = miny ivp = maxy ihh = minx ihp = maxx ; TWIN, /new, xsize=ih, ysize=iv mxdif = mxbv - mnbv mxd3 = mxdif/3.0 mx = mxbv + mxd3 mn = mnbv - mxd3 MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Reddish(0) + bytscl( a, max = mx, $ min = mn, top = Reddish(1)-Reddish(0)) TV, a0 dif = maxc - minc FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) if azb gt azb2 and azb lt azb1 then begin if htb gt htb1 and htb lt htb2 then begin va = value(ixx,iyy) if dif gt 39 then begin if va lt 0 then begin vlu = minc + bytscl( [va], max = 0.0, min = minva, $ top = 39) endif if va ge 0 then begin vlu = maxc - bytscl( [va], max = maxva, min = 0.0, $ top = 39) endif endif else begin vlu = minc + bytscl( [va], max = maxva, min = minva, $ top = maxc-minc) endelse TV, vlu, ixx, iyy endif endif ENDFOR ENDFOR ; endif ; END ; function imagea,ih,iv,xc,yc,scx,scy,az1,az2,ht1,ht2,imgg, $ xsiz,ysiz,az1i,ht2i,ti1,ti2,dtime,ps ; ; this overplots the "area to be analyzed lines" on the image as well as an ; "OK" box on the bottom right of the image. ; if imgg eq 1 then begin ; radtd = 57.29578 ; naz = abs(az1-az2) htx1 = fltarr(naz+1) hty1 = fltarr(naz+1) htx2 = fltarr(naz+1) hty2 = fltarr(naz+1) ; for i=0L,naz do begin si = sin((az1-i)/radtd) co = cos((az1-i)/radtd) htx1(i) = xc + (ht1/scx)*co hty1(i) = yc + (ht1/scy)*si htx2(i) = xc + (ht2/scx)*co hty2(i) = yc + (ht2/scy)*si endfor ; azxp1 = fltarr(2) azyp1 = fltarr(2) azxp2 = fltarr(2) azyp2 = fltarr(2) ; abox = fltarr(5) aboy = fltarr(5) ; aboxs = fltarr(5) aboys = fltarr(5) ardisx = fltarr(2) ardisy = fltarr(2) ; azxp1(0) = htx1(0) azyp1(0) = hty1(0) azxp1(1) = htx2(0) azyp1(1) = hty2(0) azxp2(0) = htx1(naz) azyp2(0) = hty1(naz) azxp2(1) = htx2(naz) azyp2(1) = hty2(naz) ; abox(0) = ih+2 abox(1) = ih+2 abox(2) = ih+18 abox(3) = ih+18 abox(4) = ih+2 aboy(0) = 2 aboy(1) = 18 aboy(2) = 18 aboy(3) = 2 aboy(4) = 2 ; dsp = ps*(ti2-ti1) ; preferred distance travelled ; dsp100 = 100.*(dtime) ; 100 km/s distance travelled pdistx = xc + (ht2i+dsp)/scx*cos(az1i/radtd) ; preferred x distance pdisty = yc + (ht2i+dsp)/scy*sin(az1i/radtd) ; preferred y distance aboxs(0) = pdistx aboxs(1) = pdistx xsi = xsiz if az1i gt 90. then xsi = -xsiz aboxs(2) = pdistx + xsi aboxs(3) = pdistx + xsi aboxs(4) = pdistx aboys(0) = pdisty ysi = ysiz if az1i gt 180. then ysi = -ysi aboys(1) = pdisty + ysi aboys(2) = pdisty + ysi aboys(3) = pdisty aboys(4) = pdisty ; ardisx(0) = pdistx + xsi ; ardisx(1) = pdistx + xsi + dsp100/scx*cos(az1i/radtd) ; ardisy(0) = pdisty + ysi ; ardisy(1) = pdisty + ysi + dsp100/scy*sin(az1i/radtd) ; oplot, azxp1, azyp1, thick = 2, /noclip oplot, azxp2, azyp2, thick = 2, /noclip oplot, htx1, hty1, thick = 2, /noclip oplot, htx2, hty2, thick = 2, /noclip oplot, abox,aboy, thick = 2, /noclip oplot, aboxs,aboys, thick = 2, /noclip ; oplot, ardisx,ardisy, thick = 2, /noclip ; endif ; END ; function imageb,ih,iv,xc,yc,scx,scy,az1,az2,ht1,ht2,imgg, $ xsiz,ysiz,az1i,ht2i,ti1,ti2,dtime,ps, $ minv,maxv,minc,maxc ; ; this overplots the "area to be analyzed lines" on the image as well as a ; color label box at the bottom right of the image. ; if imgg eq 1 then begin radtd = 57.29578 ; naz = abs(az1-az2) htx1 = fltarr(naz+1) hty1 = fltarr(naz+1) htx2 = fltarr(naz+1) hty2 = fltarr(naz+1) ; for i=0L,naz do begin si = sin((az1-i)/radtd) co = cos((az1-i)/radtd) htx1(i) = xc + (ht1/scx)*co hty1(i) = yc + (ht1/scy)*si htx2(i) = xc + (ht2/scx)*co hty2(i) = yc + (ht2/scy)*si endfor ; azxp1 = fltarr(2) azyp1 = fltarr(2) azxp2 = fltarr(2) azyp2 = fltarr(2) ; abox = fltarr(5) aboy = fltarr(5) ; aboxs = fltarr(5) aboys = fltarr(5) ardisx = fltarr(2) ardisy = fltarr(2) ; azxp1(0) = htx1(0) azyp1(0) = hty1(0) azxp1(1) = htx2(0) azyp1(1) = hty2(0) azxp2(0) = htx1(naz) azyp2(0) = hty1(naz) azxp2(1) = htx2(naz) azyp2(1) = hty2(naz) ; ht = maxc - minc abox(0) = ih-18 abox(1) = ih-18 abox(2) = ih-2 abox(3) = ih-2 abox(4) = ih-18 aboy(0) = 20 aboy(1) = ht + 20 aboy(2) = ht + 20 aboy(3) = 20 aboy(4) = 20 ; scale1 = fltarr(16,ht) FOR j=0L,ht-1 do BEGIN FOR i=0L,15 do BEGIN scale1(i,j) = minc + j if j gt 39 then scale1(i,j) = maxc - j + 39 endfor ; print, ht, j, minc + j, maxc - j + 39 endfor ; dsp = ps*(ti2-ti1) ; preferred distance travelled pdistx = xc + (ht2i+dsp)/scx*cos(az1i/radtd) ; preferred x distance pdisty = yc + (ht2i+dsp)/scy*sin(az1i/radtd) ; preferred y distance ; aboxs(0) = pdistx aboxs(1) = pdistx xsi = xsiz if az1i gt 90. then xsi = -xsiz aboxs(2) = pdistx + xsi aboxs(3) = pdistx + xsi aboxs(4) = pdistx aboys(0) = pdisty ysi = ysiz if az1i gt 180. then ysi = -ysi aboys(1) = pdisty + ysi aboys(2) = pdisty + ysi aboys(3) = pdisty aboys(4) = pdisty ; scl = minc + bytscl( scale1, max = maxc, $ min = minc, top = maxc-minc) TV, scl, ih-18,20 oplot, azxp1, azyp1, thick = 2, /noclip oplot, azxp2, azyp2, thick = 2, /noclip oplot, htx1, hty1, thick = 2, /noclip oplot, htx2, hty2, thick = 2, /noclip oplot, abox,aboy, thick = 2, /noclip oplot, aboxs,aboys, thick = 2, /noclip ; endif ; END function imagec,hsize,vsize,minv,maxv,minc,maxc,imgg ; ; This function plots the scale on the box at the left of the image. ; if imgg eq 1 then begin ; if minv ge 0.0 then min = fix(minv) if minv lt 0.0 then min = - fix(abs(minv)) max = fix(maxv) dif = maxc - minc x1 = hsize - 85 y1 = 17 y2 = dif + 17 xyouts, x1, y1, min xyouts, x1, y2, max ; endif ; END ; function nextim,az1,az2,ht1,ht2,filta,filtb, $ azd,htd,az11,az22,ht11,ht22 ; ; This function determines the "area to be analyzed locations" on the second ; image ; az11 = az1 + azd + filta az22 = az2 + azd - filta ht11 = ht1 + htd - filtb ; ht in km, time in s. ht22 = ht2 + htd + filtb ; ht in km, time in s. ; end ; function param,name,img,istep,fini1,name1,h1,v1,hc,vc,scx,scy,doy,tim ; ; This function reads the IDL data file and for this image, determines ; its important parameters and returns the image. ; ; name of differenced image ; ; ir = lasco_readfits(name1,hdr) ; namep = namep+strcompress(image,/rem)+finif ; namep = name image = img+istep name1 = namep+strcompress(image,/rem)+fini1 a = lasco_readfits(name1,hdr) ; h1 = hdr.NAXIS1 ; NAXIS1 horizontal size of differenced image v1 = hdr.NAXIS2 ; NAXIS2 vertical size of differenced image sun = getsuncen(hdr) hc = sun.xcen vc = sun.ycen scxas = hdr.PLATESCL ; CDELT1 horizontal scale of image (arcsec/pixel) scyas = hdr.PLATESCL ; CDELT2 vertical scale of image (arcsec/pixel) doy = hdr.MID_DATE ; MID_DATE time of image in Modified Julian Day tim = hdr.MID_TIME ; MID_TIME time of day of image in miliseconds ; ; Determine rsun = sun size in arc seconds ; date = {mjd:hdr.MID_DATE,time:hdr.MID_TIME*1000.} rsu = sun_ephem(date) rsun = rsu(4)*60.0 ; convert to size in arcseconds print, rsun ; rs = 695980.0 ; Sun radius in km kmpas = rs/rsun ; Km on sun per arcsec scx = kmpas*scxas scy = kmpas*scyas print, name1,h1,v1,hc,vc,scx,scy,doy,tim return, a ; end ; pro cori33,name,img,cor=cor,istep=istep,psp=psp, $ sspeed=sspeed,fspeed=fspeed,smax=smax,angl=angl, $ az1=az1,az2=az2,ht1=ht1,ht2=ht2,ok=ok,imgg=imgg, $ xsize=xsize,ysize=ysize,iave=iave,nimage=nimage,nsmo=nsmo ; ; This program takes image sets from LASCO, uses the best current technique ; to difference the data from a base, and then displays this data. Using a ; mouse cursor the researcher defines an area on the image to track. This ; area is drawn on the next image and is also defined with outward motion on ; this image. The program figures out how many different preferred speeds to ; try within the limits of sspeed and fspeed and then trys them on a single ; image set before continuing. The best three answers within the given area ; are averaged together to determine the correct answer to track to the next ; image. The program redefines this box on the next image and it is ; redrawn. The LASCO data is filtered using the best available filter, and ; the following images within that same sector in the data set are, too. A ; two-D correlation program is used to follow the best correlation in that ; box outward from that area in following images. ; if n_elements(name) eq 0 then name = 'jan15_' ; Jan. 15 date default if n_elements(cor) eq 0 then cor = 3 ; C3 coronagraph default if n_elements(img) eq 0 then img = 2444 ; Image 2447 default if n_elements(istep) eq 0 then istep = 4 ; Default image increment if n_elements(psp) eq 0 then psp = 300.0 ; Pref. default speed ; ; Size of the correlation window used (must be an odd number) ; if n_elements(xsize) eq 0 then xsize = 9 ; horizontal size if n_elements(ysize) eq 0 then ysize = xsize ; vertical size ; if n_elements(angl) eq 0 then angl = 20.0 ; max angle from radial if n_elements(sspeed) eq 0 then sspeed = 50.0 ; slowest speed allowed if n_elements(fspeed) eq 0 then fspeed = 1000.0 ; fastest speed allowed if n_elements(smax) eq 0 then smax = 0.1 ; lowest corr. allowed ; if n_elements(iave) eq 0 then iave = 5 ; filter tangential ; pixel average away ; from center ; if n_elements(az1) eq 0 then az1 = 50 ; bigger azimuth (1) if n_elements(az2) eq 0 then az2 = 40 ; lesser azimuth (2) if n_elements(ht1) eq 0 then ht1 = 4.0 ; smaller height (1) if n_elements(ht2) eq 0 then ht2 = 5.0 ; greater height (2) if n_elements(nimage) eq 0 then nimage = 33 ; number of image sets if n_elements(nsmo) eq 0 then nsmo = 0 ; image smoothing factor rs = 695980.0 ; Sun radius in km ht1 = ht1*rs ; ht1 in km ht2 = ht2*rs ; ht2 in km ; if n_elements(ok) eq 0 then ok = 0 ; if 1, automatic run if n_elements(imgg) eq 0 then imgg = 1 ; if 1, images plotted iok = 0 if imgg ne 1 then iok = 1 badspv = -5.0 ; value in array if speed correlation is out of limit badang = 999.0 ; value in array if speed correlation is out of limit badamp = -5.0 ; value in array if speed correlation is out of limit ; dayins = 24.*60.*60. ; number of seconds in a day radtd = 57.29578 ; degrees in a radian itrigger = 0 ; dir1 = '[bjackson.lasco.dat]' dir2 = '[bjackson.lasco.results]' fini1 = '.dif' fini2 = '.dat' openw, one, dir2+name+strcompress(img,/rem)+fini2,/get_lun print, ' RUN ON ', name+strcompress(img,/rem)+fini2 printf, one, ' RUN ON ', name+strcompress(img,/rem)+fini2 ; name = dir1+name ; npspn = 0 pcos = 1.0 psin = 0.0 pang = 0.0 htm = (ht1 + ht2)/2.0 ; start: ; ; sspeed = psp - psp/2.0 ; if sspeed lt 50.0 then sspeed = 50 ; fspeed = psp+psp ; ; Determine different parameters for this specific image. ; is = 0 ; par = param(name,img,is,fini1,namea,hsize1,vsize1,hcent1,vcent1,scx1,scy1, $ day1,time1) ; ; Plot the differenced image on a window ; ; May 24 mxbv = 611 mnbv = 499 ; Jan 15 mxbv = 400 mnbv = 0 ; up: img1 = image(namea,hsize1,vsize1,mxbv,mnbv,imgg,a) ; ; Determine the area of the image to analyze (azimuth and height values) ; if imgg eq 1 then begin getc = getcp(hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,ok) endif ; ; Determine the area of the second image to analyze. ; ; Plot the first differenced image on a window with the area ; superimposed. ; dtime = 0.0 ; fd = findmx(az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1,xsize,ysize,iave, $ minx,maxx,miny,maxy) ; mxmn = mxmnbv(a,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1, $ minx,maxx,miny,maxy,mxbv,mnbv) print, 'az1 az2',az1,az2 print, mxbv,mnbv ; imga=imagea(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time1,dtime,psp) print, 'az1 az2',az1,az2 if imgg eq 1 then iok = getok(hsize1,vsize1,ok) if not iok then goto, up ; ; Determine the parameters for the second image. ; par = param(name,img,istep,fini1,nameb,hsize2,vsize2,hcent2,vcent2, $ scx2,scy2,day2,time2) add = 0.0 if day2 lt day1 then add = (day1-day2)*dayins if day2 gt day1 then add = (day2-day1)*dayins time2 = time2 + add ; time of image 2 in seconds dtime = time2-time1 ; time difference between images ; ; ; determine the updated height and azimuth differences to be used on ; the second image for the next correlations ; htd = psp*dtime*pcos ; html = htm htm = (ht1 + ht2)/2.0 htmmm = (html + htm)/2.0 ; azd = atan((psin*htd),(htm + htd))*radtd ; ; Determine the next image area to filter and analyze ; filta = abs(atan((1.5*xsize),ht1/scx1))*radtd ; extra azimuth value in deg filtb = 1.5*xsize*scx1 ; extra height value in km ; up2: ; next = nextim(az1,az2,ht1,ht2,filta,filtb, $ azd,htd,az11,az22,ht11,ht22) print, 'az11 az22',az11,az22 ; ; Find the min and max x and y values for the area analyzed on the second image ; fd=findmx(az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2,xsize,ysize,iave, $ minx11,maxx11,miny11,maxy11) ; spcw = xsize*scx1/dtime ; ; Plot the second differenced image on a window with the new area ; superimposed. ; ; May 24 mxbv = 611 mnbv = 499 ; Jan 15 mxbv = 40 mnbv = 20 ; img2 = image(nameb,hsize2,vsize2,mxbv,mnbv,imgg,b) ; ; locate the updated height and the azimuth location on the second image ; ht1d = ht1 + htd ht2d = ht2 + htd az1d = az1 + azd az2d = az2 + azd ; imga=imagea(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az1d,az2d,ht1d,ht2d,imgg, $ xsize,ysize,az1,ht2,time1,time2,dtime,psp) imga=imagea(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time2,dtime,psp) imga=imagea(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az11,az22,ht11,ht22,imgg, $ xsize,ysize,az1,ht2,time1,time2,dtime,psp) print, 'az1 az2',az1,az2 print, 'az11 az22',az11,az22 print, 'Speed and angle on the image',psp,pang ; if imgg eq 1 then iok = getok(hsize2,vsize2,ok) if not iok and imgg eq 1 then begin gsp = getspeed(hcent2,vcent2,scx2,scy2,ht1,ht2,dtime,psp,ok) goto, up2 endif ok = 1 ; ; Filter the first and second image areas. The area beyond this region is ; not filtered. ; ; First image ; smo = smoo(a,hsize1,vsize1,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,ar) if nsmo ge 3 then begin ar = smooth(ar,nsmo) print, 'Image 1 smoothed',nsmo,' by',nsmo printf, one, 'Image 1 smoothed',nsmo,' by',nsmo endif ; ft1 = filtr(ar,hsize1,vsize1,iave,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,i1) ; ; Second image ; smo = smoo(b,hsize2,vsize2,xsize,ysize,hcent2,vcent2,scx2,scy2, $ az11,az22,ht11,ht22,minx11,maxx11,miny11,maxy11,br) if nsmo ge 3 then begin br = smooth(br,nsmo) print, 'Image 2 smoothed',nsmo,' by',nsmo printf, one, 'Image 2 smoothed',nsmo,' by',nsmo endif ; ft2 = filtr(br,hsize2,vsize2,iave,xsize,ysize,hcent2,vcent2,scx2,scy2, $ az11,az22,ht11,ht22,minx11,maxx11,miny11,maxy11,i2) ; ; Correlate from one location to the other in the two images using the ; nominal speed value. ; nkx = long(maxx-minx+1) nky = long(maxy-miny+1) nk = nkx*nky print, 'nk', nk spv = fltarr(hsize1,vsize1) ampv = fltarr(hsize1,vsize1) angv = fltarr(hsize1,vsize1) spv(*,*) = sspeed ampv(*,*) = smax angv(*,*) = 0.0 spvn = fltarr(hsize1,vsize1) ampvn = fltarr(hsize1,vsize1) angvn = fltarr(hsize1,vsize1) spvn(*,*) = sspeed ampvn(*,*) = smax angvn(*,*) = 0.0 pspeedxy = fltarr(nk) xcent = fltarr(nk) ycent = fltarr(nk) ; pspeedxy(*) = badspv npspn = npspn + 1 ; hh1 = time1/(60.*60.) hh1 = fix(hh1) mm1f = time1/60. - hh1*60. mm1 = fix(mm1f) ss1 = mm1f*60. - mm1*60. hh2 = time2/(60.*60.) hh2 = fix(hh2) mm2f = time2/60. - hh2*60. mm2 = fix(mm2f) ss2 = mm2f*60. - mm2*60. print, ' ' print, 'first image is at ', hh1, mm1, ss1 printf, one, 'first image is at', hh1, mm1, ss1 print, 'second image is at', hh2, mm2, ss2 printf, one, 'second image is at', hh2, mm2, ss2 print, 'correlation window width is equivalent to ',spcw,' km/s' print, 'azimuths ', az1, az2,' hts ', ht1/rs, ht2/rs printf, one, 'correlation window width is equivalent to ',spcw,' km/s' printf, one, 'azimuths', az1, az2,' hts', ht1/rs, ht2/rs ; bestsp0 = 0 bestsp1 = 0 bestsp2 = 0 bestamp0 = 0 bestamp1 = 0 bestamp2 = 0 bestang0 = 0 bestang1 = 0 bestang2 = 0 ; ; Determine the best preferred speed values to use to propagate to the next ; image (psp is cycled from sspeed to fspeed at distances of 1/3 xsize). ; sxsize3 = spcw/3.0 ; 1/3 the speed equivalent to the length of xsize psp = sspeed again: print, ' psp', psp ; cx = cimage(i1,i2,hsize1,vsize1,hsize2,vsize2,hcent1,vcent1,hcent2,vcent2, $ scx1,scy1,scx2,scy2,psp,dtime,az1,az2,ht1,ht2, $ minx,maxx,miny,maxy, $ pspeedxy,xcent,ycent,kk,itrigger,imgg, $ az11,az22,ht11,ht22,xsize,ysize,one, $ angl,sspeed,fspeed,smax,badspv,badang,badamp, $ spv,angv,ampv) ; Determine the highest amplitude correlation of speed and angle value ; within the area defined ; bv = bestv(spv,ampv,angv,spvn,ampvn,angvn,hsize1,vsize1,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,bestsp0,bestsp1,bestsp2, $ bestamp0, bestamp1, bestamp2, bestang0, bestang1, bestang2) ; if psp gt fspeed then goto, next psp = psp + sxsize3 goto, again next: print, 'speed0 =', bestsp0 printf, one, 'speed0 =', bestsp0 print, 'speed1 =', bestsp1 printf, one, 'speed1 =', bestsp1 print, 'speed2 =', bestsp2 printf, one, 'speed2 =', bestsp2 ; print, 'amplitude0 =', bestamp0 printf, one, 'amplitude0 =', bestamp0 print, 'amplitude1 =', bestamp1 printf, one, 'amplitude1 =', bestamp1 print, 'amplitude2 =', bestamp2 printf, one, 'amplitude2 =', bestamp2 ; print, 'angle0 =', bestang0 printf, one, 'angle0 =', bestang0 print, 'angle1 =', bestang1 printf, one, 'angle1 =', bestang1 print, 'angle2 =', bestang2 printf, one, 'angle2 =', bestang2 ; ; Plot the best speed values within the correlated area superposed on the ; image ; if imgg eq 1 then begin MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; mincs = Blueish(0) maxcs = Blueish(1) mincg = Reddish(0) maxcg = Greenish(1) mincp = Reddish(0) maxcp = Reddish(1) endif ; ; May 24 mxbv = 611 mnbv = 499 ; Jan 15 mxbv = 40 mnbv = 20 ; img10 = image10(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1, $ imgg,minx,maxx,miny,maxy,spvn,sspeed,fspeed,mxbv,mnbv,mincs,maxcs) ; ; Overplot the boundary and the speed color scale box ; imgb = imageb(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time1,dtime,psp, $ sspeed,fspeed,mincs,maxcs) ; ; Overplot the speed scale on the box ; imgc = imagec(hsize1,vsize1,sspeed,fspeed,mincs,maxcs,imgg) ; ; Plot the best amplitude values within the correlated area superposed on the ; image ; img10 = image10(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1,imgg, $ minx,maxx,miny,maxy,ampvn,0.0,1.0,mxbv,mnbv,mincp,maxcp) ; ; Overplot the boundary and the amplitude color scale box ; imgb = imageb(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time1,dtime,psp, $ 0.0,1.0,mincp,maxcp) ; ; Overplot the amplitude scale on the box ; imgc = imagec(hsize1,vsize1,0.0,1.0,mincp,maxcp,imgg) ; ; Plot the best angle values within the correlated area superposed on the ; image ; img10 = image10(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1,imgg, $ minx,maxx,miny,maxy,angvn,-angl,angl,mxbv,mnbv,mincg,maxcg) ; ; Overplot the boundary and the angle color scale box ; imgb = imageb(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time1,dtime,psp, $ -angl,angl,mincg,maxcg) ; ; Overplot the angle scale on the box ; imgc = imagec(hsize1,vsize1,-angl,angl,mincg,maxcg,imgg) ; print, 'npspn nimage', npspn, nimage ; ; Do the correlations once more if npspn lt nimage ; psp = (bestsp0*bestamp0 + bestsp1*bestamp1 + bestsp2*bestamp2)/ $ (bestamp0 + bestamp1 + bestamp2) pang = (bestang0*bestamp0 + bestang1*bestamp1 + bestang2*bestamp2)/ $ (bestamp0 + bestamp1 + bestamp2) ; img = img + istep ; New starting image ; ; determine the updated height and azimuths to be used on the first image ; of the next set for the next correlations ; pcos = cos(pang/radtd) psin = sin(pang/radtd) ; if pcos ne 0 then begin htd = psp*dtime*pcos endif ; htm = (ht1 + ht2)/2.0 ; ht1 = ht1 + htd ; New starting height 1 ht2 = ht2 + htd ; New starting height 2 ; azd = atan((psin*htd),(htm + htd))*radtd ; az1 = az1 + azd ; New starting azimuth 1 az2 = az2 + azd ; New starting azimuth 2 print, 'Best speed', psp,' at angle', pang,' at', htmmm/rs,' Rs' printf, one, 'Best speed', psp,' at angle', pang,' at', htmmm/rs,' Rs' print, 'New img az1 az2 ht1 ht2 psp pang',img,az1,az2,ht1/rs,ht2/rs,htm/rs,psp,pang ; if npspn ge nimage then goto, end1 ; last time through goto, start ; another time ; end1: ; free_lun, one ; end