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,nx,ny,xc,yc,scx,scy,xsize,ysize, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,bestsp0,bestsp1,bestsp2, $ bestamp0,bestamp1,bestamp2,bestang0,bestang1,bestang2 minxx = xsize/3.0 minyy = ysize/3.0 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 ; print, 'bestv az1 az2 ht1 ht2',az1,az2,ht1,ht2 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) if azb lt az2 then azb = azb + 360.0 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 bestamp0 then begin bestamp0 = amp(ixx,iyy) bestsp0 = sp(ixx,iyy) bestang0 = ang(ixx,iyy) bx0 = ixx by0 = iyy ; print, '1st ixx iyy azb htb', ixx,iyy,azb,htb,amp(ixx,iyy),bestamp0 endif endif endif ENDFOR ENDFOR 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, '2nd ixx iyy azb htb', ixx,iyy,azb,htb,amp(ixx,iyy),bestamp1 if azb lt az2 then azb = azb + 360.0 if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin if bx0 ne -1 and by0 ne -1 then begin xdif = abs(ixx - bx0) ydif = abs(iyy - by0) ; print, 'ixx iyy xdif ydif bestamp1',ixx,iyy,xdif,ydif,bestamp1 if amp(ixx,iyy) gt bestamp1 then begin if xdif gt minxx or ydif gt minyy then begin bestamp1 = amp(ixx,iyy) bestsp1 = sp(ixx,iyy) bestang1 = ang(ixx,iyy) bx1 = ixx by1 = iyy endif endif endif endif endif ENDFOR ENDFOR 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 lt az2 then azb = azb + 360.0 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 then begin if xdif gt minxx or ydif gt minyy then begin if xdif1 gt minxx or ydif1 gt minyy then begin if xdif2 gt minxx or 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 endif endif endif ENDFOR ENDFOR ; end ; function cimage4,i1,i2,nh1,nv1,nh2,nv2,xc1,yc1,xc2,yc2, $ scx1,scy1,scx2,scy2,pspeed,dtime,xpm,ypm, $ minx,maxx,miny,maxy, $ 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 ; azb11 = az11 htb11 = ht11 azb22 = az22 htb22 = ht22 ; ; ; Arrays into which to place correlated parameters ; cspeed = fltarr(nh2,nv2) cangl = fltarr(nh2,nv2) cmax = fltarr(nh2,nv2) ; ; Zero arrays of correlated values ; nh2m = nh2 - 1 nv2m = nv2 - 1 FOR i=0L,nh2m do BEGIN FOR j=0L,nv2m do BEGIN cspeed(i,j) = badspeed ; no correlation speed cangl(i,j) = badangle ; no angle between radial and corr. cmax(i,j) = badcorr ; no correlation value ENDFOR ENDFOR ; ; Correlation arrays ; ; print, xsize2p,ysize2p i11 = fltarr(xsize2p,ysize2p) i22 = fltarr(xsize2p,ysize2p) ; print, ' Enter correlation' ; nxxx = 0L nxxxx = 0L ix=fix(xpm) iy=fix(ypm) ; print, 'ix iy',ix,iy ivv = iy-ysize ivp = iy+ysize ihh = ix-xsize ihp = ix+xsize ; ; 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 ; ; print, miny,maxy,minx,maxx,jstep,istep FOR iyyy=miny,maxy,jstep do BEGIN FOR ixxx=minx,maxx,istep do BEGIN iflag1 = 0 iflag2 = 0 iflag3 = 0 iflag4 = 1 ; ; Determine if the boundary is within the boundary correlation of image 2 ; bd2 = bound(ixxx,iyyy,xc2,yc2,scx2,scy2,azb,htb) ; print, 'azb azb11 azb22 htb htb11 htb22', azb,azb11,azb22,htb,htb11,htb22 if azb lt azb22 then azb = azb + 360.0 if azb ge azb22 and azb le azb11 then begin if htb ge htb11 and htb le htb22 then begin ; print, 'azb azb11 azb22 htb htb11 htb22', azb,azb11,azb22,htb,htb11,htb22 ; ; Load up the second array to use in the correlation. ; ; Calculate angle of center of new array from center of old ; xphi = (ix-xc1) yphi = (iy-yc1) phin = atan(yphi*scy1,xphi*scx1)*radtd phinn = phin ithpflag = 0 if phin lt 0.0 then begin phin = phin + 360.0 ithpflag = 1 endif xxphi = ixxx-ix yyphi = iyyy-iy thep = atan((yyphi)*scy1,(xxphi)*scx1)*radtd if thep lt 0.0 then begin thep = thep + 360.0 if ithpflag eq 0 then phin = phin + 360.0 endif dangn = (thep - phin) xyspnn = sqrt(xxphi*xxphi*scxsq + yyphi*yyphi*scysq) xyspn = xyspnn/dtime ; print, '2nd ixxx iyyy ivv ivp ihh ihp nxxx dangn xyspn',ixxx,iyyy,ivv,ivp,ihh,ihp,nxxx,dangn,xyspn ; ; Do not accept array unless it is within +-angl of old or within xsize/3 ; pixels of the old and within mionimum and maximum speed limits ; sdangn = abs(sin(phin/radtd))*xyspnn ; spdi = (((xsize/3.0)*scx1)/dtime) fspeedp = fspeed + spdi sspeedm = sspeed - spdi ; if abs(dangn) lt angl or sdangn lt (xsize/3.0)*scx1 then begin if xyspn lt fspeedp and xyspn gt sspeedm then begin ; ivv = iyyy-ysize ivp = iyyy+ysize ihh = ixxx-xsize ihp = ixxx+xsize ; print, '2nd ixxx iyyy ivv ivp ihh ihp nxxx dangn xyspn',ixxx,iyyy,ivv,ivp,ihh,ihp,nxxx,dangn,xyspn 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 ixh lt (xsize+isizeh) then begin ; if iyv lt (ysize+jsizeh) then begin ; if i11(ixh,iyv) ne i11(isizeh,jsizeh) then iflag2 = 1 ; if i22(ixh,iyv) ne i22(isizeh,jsizeh) then iflag3 = 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 nxxx', xxi,yyj,amax,nxxx ; if abs(xxi) gt xsize/3 or abs(yyj) gt xsize/3 then begin xys = badspeed dang = badangle amax = badcorr endif else begin phi = phinn ithpflag = 0 if phi lt 0.0 then begin phi = phi + 360.0 ithpflag = 1 endif xx = xxphi - xxi yy = yyphi - yyj the = atan(yy,xx)*radtd if the lt 0.0 then begin the = the + 360.0 if ithpflag 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 endelse ; print, 'xys dang amax',xys,dang,amax ; 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 ; print, 'xys dang amax ixxx iyyy',xys,dang,amax,ixxx,iyyy ; ; **************************************************************************** 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 ; ; **************************************************************************** ; cspeed(ixxx,iyyy) = xys cangl(ixxx,iyyy) = dang cmax(ixxx,iyyy) = amax endif endif ; end of boundary if endif ; end of boundary if ; 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 lt azb2 then azb = azb + 360.0 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,ih,iv,az1,az2,ht1,ht2,xc,yc,scx,scy, $ minx,maxx,miny,maxy,mxbv,mnbv,imgg ; ; This function finds the minimum and maximum image values within the area ; bound by azimuth and height. ; if imgg le 4 then begin radtd = 57.29578 ; if minx lt 0 then minx=0 if maxx ge ih then maxx = ih - 1 if miny lt 0 then miny=0 if maxy ge iv then maxy = iv - 1 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 lt az2 then azb = azb + 360.0 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 endif 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 lt azb2 then azb = azb + 360.0 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 lt azb2 then azb = azb + 360.0 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,xp1,yp1,azm,htm,iok ; ; Determine an image area to analyze in the study. ; radtd = 57.29578 ; if iok eq 0 then begin ; pcursor, xp11, yp11, /data ; xp1 = xp11 yp1 = yp11 scx2 = scx*scx scy2 = scy*scy azm = atan((yp1 - yc),(xp1 - xc))*radtd htm = sqrt((yp1 - yc)*(yp1 - yc)*scy2 + (xp1 - xc)*(xp1 - xc)*scx2) print, azm, htm 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 le 5 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 le 4 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 4 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 lt azb2 then azb = azb + 360.0 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 va lt minva or va gt maxva then begin ; use default if O.O.B vdef = Greenish(0) TV, [vdef], ixx, iyy endif else begin 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 ; print, va, vlu, ixx, iyy TV, vlu, ixx, iyy endelse 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,ibound ; ; this overplots the "area to be analyzed lines" (if ibound = 1) on the image ; as well as a square the size of correlation box and an "OK" box on the ; bottom right of the image. ; if imgg eq 5 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 pdistx = xc + (ht2i+dsp)/scx*cos(az1i/radtd) ; preferred x distance pdistx = pdistx - xsiz/2.0 pdisty = yc + (ht2i+dsp)/scy*sin(az1i/radtd) ; preferred y distance pdisty = pdisty - ysiz/2.0 ; aboxs(0) = pdistx aboxs(1) = pdistx xsi = xsiz aboxs(2) = pdistx + xsi aboxs(3) = pdistx + xsi aboxs(4) = pdistx aboys(0) = pdisty ysi = ysiz aboys(1) = pdisty + ysi aboys(2) = pdisty + ysi aboys(3) = pdisty aboys(4) = pdisty ; if ibound eq 1 then begin oplot, azxp1, azyp1, thick = 2, /noclip oplot, azxp2, azyp2, thick = 2, /noclip oplot, htx1, hty1, thick = 2, /noclip oplot, htx2, hty2, thick = 2, /noclip endif oplot, abox,aboy, thick = 2, /noclip oplot, aboxs,aboys, 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 le 4 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 le 4 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,azm1,htm1,azm2,htm2,hc,vc,scx,scy, $ filta,filth1,filth2,az11,az22,ht11,ht22 ; ; This function determines the "area to be analyzed locations" on the second ; image ; radtd = 57.29578 ; az11 = azm1 + filta az22 = azm2 - filta ht11 = htm1 + filth1 ; ht in km, time in s. ht22 = htm2 + filth2 ; ht in km, time in s. ; ; ypmm = sin(azm/radtd)*htm/scy + vc ; xpmm = cos(azm/radtd)*htm/scx + hc ; print, ypmm,xpmm,azm/radtd,htm/scy,vc ; 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 cori44,name,img,cor=cor,istep=istep,psp=psp,spfrac=spfrac, $ smax=smax,angl=angl, $ azm=azm,htm=htm,htmax=htmax,ok=ok,imgg=imgg, $ xsize=xsize,ysize=ysize,iave=iave,nimage=nimage,nsmo=nsmo ; This program IS a correlation tracker. ; ; 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 a correlation area on the image to track. ; This xsize by ysize area is drawn on the next image and is also defined with ; outward motion on the next image. The program tracks this box to the next ; image as best it can by looking within the defined area on the next image. ; The program then 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 2444 default img11 = img if n_elements(istep) eq 0 then istep = 4 ; Default image increment if n_elements(psp) eq 0 then psp = 350.0 ; Pref. default speed if n_elements(spfrac) eq 0 then spfrac = 0.2 ; Pref. speed fraction ; ; 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 = 5.0 ; max angle from radial 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(azm) eq 0 then azm = 345. ; mean starting azimuth azm = float(azm) if n_elements(htm) eq 0 then htm = 4.5 ; mean starting height if n_elements(htmax) eq 0 then htmax = 30.0 ; maximum image height if n_elements(nimage) eq 0 then nimage = 32 ; number of image sets if n_elements(nsmo) eq 0 then nsmo = 0 ; image smoothing factor ; 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 ; 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 ; rs = 695980.0 ; Sun radius in km htm = float(htm) htm = htm*rs ; htm in km htmax = float(htmax) htmax = htmax*rs ; htm in km 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' nam1 = name 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 ; pang = 0.0 azmd = 0.0 sspeed = psp - psp*spfrac fspeed = psp + psp*spfrac npspn = 0 maxini = 1000 imgs1 = fltarr(maxini) imgs2 = fltarr(maxini) psps = fltarr(maxini) pangs = fltarr(maxini) pamps = fltarr(maxini) htms = fltarr(maxini) ini = 0 ; start1: ; ; Determine different parameters for this specific image. ; is = 0 ; par = param(name,img,is,fini1,namea,hsize1,vsize1,hcent1,vcent1,scx1,scy1, $ day1,time1) ; ; Determine the default x and y pos of the azimuths and heights on the images ; if azm eq 180. or azm eq 360. or azm eq 0. then begin ypm = vcent1 xxxx = htm/scx1 if azm eq 180. then xxxx = -xxxx xpm = hcent1 + xxxx endif else begin ypm = sin(azm/radtd)*htm/scy1 + vcent1 xpm = cos(azm/radtd)*htm/scx1 + hcent1 endelse if ok eq 1 then print, 'xpm ypm azm htm',xpm,ypm,azm,htm/rs ; ; Plot the differenced image on a window ; ; May 24 if nam1 eq 'may24_' then begin mxbv1 = 611 mnbv1 = 499 endif ; Jan 15 if nam1 eq 'jan15_' then begin mxbv1 = 400 mnbv1 = 0 endif ; up: ; img1 = image(namea,hsize1,vsize1,mxbv1,mnbv1,imgg,a) ; ; ; Determine the area of the first image to analyze (azimuth and height value) ; getc1 = getcp(hcent1,vcent1,scx1,scy1,xpm,ypm,azm,htm,ok) if azm lt 0.0 then begin azm = azm + 360.0 endif ; ; Determine ht1 and ht2 from htm. htm is in km! ; ht1 = htm - (xsize/1.3)*scx1 ht2 = htm + (xsize/1.3)*scx1 az1 = azm + abs(atan((xsize/1.3)*scx1,htm)*radtd) az2 = azm - abs(atan((xsize/1.3)*scx1,htm)*radtd) if az2 lt 0.0 then begin az1 = az1 + 360.0 az2 = az2 + 360.0 endif ;; ; Determine the minimum and maximum values of the boundary for the first image ; fd = findmx(az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1,xsize,ysize,iave, $ minx,maxx,miny,maxy) ; ; Plot the first differenced image on a window with the area ; superimposed. ; dtime = 0.0 ; mxmn = mxmnbv(a,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1, $ minx,maxx,miny,maxy,mxbv1,mnbv1,imgg) print, ' 1st mxbv1 mnbv1 az1 az2 ht1 ht2',mxbv1,mnbv1,az1,az2,ht1/rs,ht2/rs ; azc = azm htc = htm ibound = 1 imga=imagea(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,azc,htc,time1,time1,dtime,psp,ibound) iok = getok(hsize1,vsize1,ok) if not iok then goto, up ; ; Determine the area of the second image to analyze. ; ; 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 next image area to filter and analyze ; dhtf = cos(pang/radtd)*psp*dtime filtax = abs(atan((1.5*xsize),htm/scx1))*radtd ; extra azimuth in deg filta = abs(atan(cos(angl/radtd)*dhtf,htm))*radtd if filtax gt filta then filta = filtax filthx = 1.5*xsize*scx1 ; extra height value in km filth1 = sspeed*dtime - filthx filth2 = fspeed*dtime + filthx ; print, 'Nextim azm htm', azm, htm/rs next = nextim(azm,htm,azm,htm,hcent2,vcent2,scx2,scy2, $ filta,filth1,filth2,az11,az22,ht11,ht22) if az22 lt 0.0 then begin az11 = az11 + 360.0 az22 = az22 + 360.0 endif az11n = az11 az22n = az22 if az11 lt az22 then begin az11 = az22n az22 = az11n endif print, 'Nextim az11 az22 ht11 ht22',az11,az22,ht11/rs,ht22/rs ; fd=findmx(az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2,xsize,ysize,iave, $ minx11,maxx11,miny11,maxy11) ; spcw = xsize*scx1/dtime ; ; May 24 if nam1 eq 'may24_' then begin mxbv2 = 611 mnbv2 = 499 endif ; Jan 15 if nam1 eq 'jan15_' then begin mxbv2 = 400 mnbv2 = 0 endif ; ; Plot the second differenced image on a window with the new area ; superimposed. ; img2 = image(nameb,hsize2,vsize2,mxbv2,mnbv2,imgg,b) ; mxmn = mxmnbv(b,hsize2,vsize2,az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2, $ minx11,maxx11,miny11,maxy11,mxbv2,mnbv2,imgg) ; ; locate the height and the azimuth locations on the second image and the ; maximum boundary. ; azc = azm htc = htm ibound = 0 imga=imagea(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az11,az22,ht11,ht22,imgg, $ xsize,ysize,azc,htc,time1,time1,dtime,psp,ibound) azc = azm + azmd htc = htm + dhtf ibound = 1 imga=imagea(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az11,az22,ht11,ht22,imgg, $ xsize,ysize,azc,htc,time1,time1,dtime,psp,ibound) ; iok = getok(hsize2,vsize2,ok) if not iok then goto, up 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(maxx11-minx11+1) nky = long(maxy11-miny11+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 ; 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, 'speed shift this time', psp print, 'correlation window width is equivalent to ',spcw,' km/s' print, '1st and 2nd image azimuth-hts',azm,htm/rs,az11,ht11/rs,az22,ht22/rs ; printf, one, 'correlation window width is equivalent to ',spcw,' km/s' ; printf,one,'1st and 2nd image azimuth-hts',azm,htm/rs,az11,ht11/rs,az22,ht22/rs ; start2: ; bestsp0 = 0 bestsp1 = 0 bestsp2 = 0 bestamp0 = 0 bestamp1 = 0 bestamp2 = 0 bestang0 = 0 bestang1 = 0 bestang2 = 0 ; cx = cimage4(i1,i2,hsize1,vsize1,hsize2,vsize2,hcent1,vcent1,hcent2,vcent2, $ scx1,scy1,scx2,scy2,psp,dtime,xpm,ypm, $ minx11,maxx11,miny11,maxy11, $ 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,hsize2,vsize2,hcent2,vcent2,scx2,scy2,xsize,ysize, $ az11,az22,ht11,ht22,minx11,maxx11,miny11,maxy11,bestsp0,bestsp1,bestsp2, $ bestamp0, bestamp1, bestamp2, bestang0, bestang1, bestang2) ; 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 ; 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 ; ; Plot the best speed values within the correlated area superposed on the ; image ; img10 = image10(nameb,hsize2,vsize2,az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2, $ imgg,minx11,maxx11,miny11,maxy11,spv,sspeed,fspeed,mxbv2,mnbv2,mincs,maxcs) ; ; ; Overplot the boundary and the speed color scale box ; imgb = imageb(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az11,az22,ht11,ht22,imgg, $ xsize,ysize,az11,ht22,time1,time1,dtime,psp, $ sspeed,fspeed,mincs,maxcs) ; ; Overplot the speed scale on the box ; imgc = imagec(hsize2,vsize2,sspeed,fspeed,mincs,maxcs,imgg) ; ; Plot the best amplitude values within the correlated area superposed on the ; image ; img10 = image10(nameb,hsize2,vsize2,az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2, $ imgg,minx11,maxx11,miny11,maxy11,ampv,0.0,1.0,mxbv2,mnbv2,mincp,maxcp) ; ; Overplot the boundary and the amplitude color scale box ; imgb = imageb(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az11,az22,ht11,ht22,imgg, $ xsize,ysize,az11,ht22,time1,time1,dtime,psp, $ 0.0,1.0,mincp,maxcp) ; ; Overplot the amplitude scale on the box ; imgc = imagec(hsize2,vsize2,0.0,1.0,mincp,maxcp,imgg) ; ; Plot the best angle values within the correlated area superposed on the ; image ; img10 = image10(nameb,hsize2,vsize2,az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2, $ imgg,minx11,maxx11,miny11,maxy11,angv,-angl,angl,mxbv2,mnbv2,mincg,maxcg) ; ; Overplot the boundary and the angle color scale box ; imgb = imageb(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az11,az22,ht11,ht22,imgg, $ xsize,ysize,az11,ht22,time1,time1,dtime,psp, $ -angl,angl,mincg,maxcg) ; ; Overplot the angle scale on the box ; imgc = imagec(hsize2,vsize2,-angl,angl,mincg,maxcg,imgg) ; pspa = (bestsp0*bestamp0 + bestsp1*bestamp1 + bestsp2*bestamp2)/ $ (bestamp0 + bestamp1 + bestamp2) panga = (bestang0*bestamp0 + bestang1*bestamp1 + bestang2*bestamp2)/ $ (bestamp0 + bestamp1 + bestamp2) ; ; (Generally only one value can be best, since the defined first area is in a ; very specific location. Thus, the average speeds determined are not used ; except as a check.) ; if bestsp0 ne 0 then begin psp = bestsp0 pang = bestang0 pamp = bestamp0 endif else begin psp = psp ; if no best speed, use the last one for the projection pang = 0.0 pamp = 0.0 endelse dpspa = 50.0 if abs(psp-pspa) gt dpspa then begin print, 'Solution not stable this time' ; printf, one, 'Solution not stable this time' endif ; ; Determine the updated height and azimuths to be used on the first image ; of the set for the next correlations. ; ; New Height ; dht = cos(pang/radtd)*psp*dtime ; htmean = htm + dht/2.0 ; Mean height for speed determination. tmean = time1 + dtime/2.0 ; Mean time for speed determination. hrmean = tmean/(60.*60.) hrmean = fix(hrmean) ; Mean hr mntmean = tmean/60. - hrmean*60. mnmean = fix(mntmean) ; Mean minute scmean = mntmean*60. - mnmean*60. ; Mean second ; htm = htm + dht ; New starting height ; ; ; New azimuth ; azmd = atan(sin(pang/radtd)*psp*dtime,htm)*radtd azmean = azm + azmd/2.0 ; Mean azimuth for speed determination. azm = azm + azmd ; ; Determine the updated fast and slow speeds to use for limits on next ; projection ; sspeed = psp - psp*spfrac fspeed = psp + psp*spfrac ; print, 'speed =',psp,' at mean height =',htmean/rs,' and mean time =',$ tmean,hrmean,mnmean,scmean,' (sec of day - hr - min - sec)' print, 'mean azimuth =',azmean ; if img eq img11 then begin ; printf, one,'image1 image2 mean time speed mean height, azimuth' printf, one,' (sec of day - hr - min - sec)' endif ; printf, one, format="(I5,3X,I5,2X,F10.2,I6,I6,F7.1,F9.2,F9.3,F13.2)",$ img,img+istep,tmean,hrmean,mnmean,scmean,psp,htmean/rs,azmean ; print, 'New image azm htm psp pang',img,azm,htm/rs,psp,pang imgs1(ini) = img imgs2(ini) = img + istep psps(ini) = psp pangs(ini) = pang pamps(ini) = pamp htms(ini) = htm/rs ini = ini + 1 if ini ge maxini then goto, endit ; total possible values ; if htm gt htmax then goto, endit ; end of track ; npspn = npspn + 1 ; print, 'npspn nimage', npspn, nimage ; if npspn ge nimage then goto, endit ; last image specified ; ; Do the correlations again if npspn lt nimage and all else is within limits. ; img = img + istep ; New starting image ; goto, start1 ; another time ; endit: ; print, 'Image and speed summary' ; printf, one, 'Image and speed summary' for i=0,ini-1 do begin print, imgs1(i), imgs2(i), psps(i), pangs(i), pamps(i), htms(i) ; printf, one, imgs1(i), imgs2(i), psps(i), pangs(i), pamps(i), htms(i) endfor; ; free_lun, one ; end