pro prep_vectmag, imgq, imgu, imgv, imgi, bx, by, bz, bpx, bpy, qq = qq, uu = uu, vv = vv, $ deriv_pol = deriv_pol, ct_cutoff = ct_cutoff, plot_ct = plot_ct, $ lscale = lscale, lzero = lzero, tscale = tscale, tzero = tzero, $ aoffset = aoffset, qreverse = qreverse, ureverse = ureverse ;+ ;PURPOSE: ; derive Bx, By, Bz from Q/U/V and display the vector map ; the procedure calls sub-routines FIT_LINEAR.PRO to do ; least-square linear fit; XYOFF1.PRO to obtain flat shift ; between images; LFF_FIELD.PRO to obtain linear-force-free magnetic ; field extrapolation. ;INPUT: ; imgq: Q component, flatfielding done; ; imgu: U component, flatfielding done; ; imgv: V component, flatfielding done; ; imgi: I (intensity) ;OUTPUT: ; bz: longitudinal field component ; bx: X-component of the trainsverse field. ; by: Y-component of the trainsverse field. ;OPTIONAL OUTPUTS: ; qq: calibrated Q component ; uu: calibrated U component ; vv: calibrated V component ;KEYWORDS INPUT: ; deriv_pol: if set, Q/U/V are not converted to degree of polarization yet, so do it. ; ct_cutoff: the lower/upper cutoff to derive cross-talk, in the same units of V. ; default = [0.001, 0.005] ; plot_ct: if set, plot the crosstalk fit ; lscale: scaling factor for longitudinal field calibration; default = 1. ; (the default can be modified later once we get the calibration done). ; (tentatively, I have got this value as 17769.5+/-1347.8 as of 2001 March) ; lzero: offset of longitudinal field calibration; default = 0. ; tscale: scaling factor for transverse field calibration; default = 1. ; (the default can be modified later once we get the calibration done). ; (tentatively, I have got this value as sqrt(152056.9+/-15%) as of 2001 March) ; tzero: offset of transverse field calibration; default = 0. ; aoffset: offset angle of the Zeiss pre-filter; default = 22 degree. ; qreverse: if set, q = - q; ; ureverse: if set, u = - u; according to my experience, we should set this at this moment. ;HISTORY: ; February 24 2001: simplified version of a messy package, QiuJ @BBSO ; 2001/06/07: slight change and more documentation. QiuJ@NJIT ;- ; nx = (size(imgq))(1) ny = (size(imgq))(2) qq = imgq uu = imgu vv = imgv ii = imgi if (keyword_set(deriv_pol) or (size(imgq))(0) gt 2) then begin ; if Q/U/V are not already obtained as degree of polarization: ; align Q/U/V components to I image, though usually we do not have to. ; first get the rigid shift by cross-correlation method. shift_q = xyoff1(ii, reform(qq(*, *, 0) + qq(*, *, 1)), nx - 20, ny - 20) shift_u = xyoff1(ii, reform(uu(*, *, 0) + uu(*, *, 1)), nx - 20, ny - 20) shift_v = xyoff1(ii, reform(vv(*, *, 0) + vv(*, *, 1)), nx - 20, ny - 20) ; now correct the shift for each component: for i = 0, 1 do begin qq(*, *, i) = poly_2d(qq(*, *, i),[-shift_q(0),0,1,0],[-shift_q(1),1,0,0],2, cubic = -0.5) uu(*, *, i) = poly_2d(uu(*, *, i),[-shift_u(0),0,1,0],[-shift_u(1),1,0,0],2, cubic = -0.5) vv(*, *, i) = poly_2d(vv(*, *, i),[-shift_v(0),0,1,0],[-shift_v(1),1,0,0],2, cubic = -0.5) endfor ; now derive Q/U/V as degree of polarization: qq = (qq(*, *, 0) - qq(*, *, 1))/(qq(*, *, 0)+qq(*, *, 1)) uu = (uu(*, *, 0) - uu(*, *, 1))/(uu(*, *, 0)+uu(*, *, 1)) vv = (vv(*, *, 0) - vv(*, *, 1))/(vv(*, *, 0)+vv(*, *, 1)) endif ; if keyword_set(deriv_pol) ; now reduce bias ; vv = vv - median (vv) qq = qq - median (qq) uu = uu - median (uu) ; now correct cross-talk: if n_elements(ct_cutoff) eq 0 then ct_cutoff = [0.001, 0.005] ; this is the lower/upper cutoff ; to obtain CT fit ss = where(abs(vv) gt ct_cutoff(0) and abs(vv) lt ct_cutoff(1)) fit_linear, vv(ss), qq(ss), alp_q, bet_q, sig_q, plot = plot_ct fit_linear, vv(ss), uu(ss), alp_u, bet_u, sig_u, plot = plot_ct print, '' print, 'Percentage of V into Q: '+string(alp_q*100, format = '(f5.1)') print, 'Percentage of V into U: '+string(alp_u*100, format = '(f5.1)') print, '' qq = qq - alp_q * vv uu = uu - alp_u * vv ; now obtain the logitudinal field component: if n_elements(lscale) eq 0 then lscale = 1. ; calibration scaling factor if n_elements(lzero) eq 0 then lzero = 0. ; calibration offset bz = vv * lscale + lzero ; now obtain the transverse field components: if n_elements(tscale) eq 0 then tscale = 1. ; calibration scaling factor if n_elements(tzero) eq 0 then tzero = 0. ; calibration offset. bt = sqrt(float(qq)^2+float(uu)^2) bt = sqrt(bt) * tscale + tzero ; determine the intensity ; get the orientation: if keyword_set(qreverse) then qq = - qq ; so far I recommend to take -qq if keyword_set(ureverse) then uu = - uu if n_elements(aoffset) eq 0 then aoffset = -22. ; the angle offset in the pre-filter. ;if n_elements(aoffset) eq 0 then aoffset = 0. ; somehow I find 22 ;does not work! Why?? tangle = atan(uu, qq+1.0e-5) tangle = tangle tangle = tangle*0.5 + aoffset*!dtor ; solve the 180degree ambiguity using simple potential field approximation: lff_field, bz, 0, 0, bpx, bpy, /noplot ss = where((cos(tangle)*bpx+sin(tangle)*bpy) lt 0.) tangle(ss) = tangle(ss) + !pi ; eventual output: bx = bt * cos(tangle) by = bt * sin(tangle) ; here to make the figures look not too small; be careful that this ; gives the polarization degree in percentage. If proper calibration ; is done, the output is in Gauss and the following should NOT be done. if (fix(lscale+0.1) eq 1) then begin ; ok here we know no calib. is given, so give an arbitrary. bz = bz * 100. bx = bx * 100. by = by * 100. endif end ;************************************************************************ pro lff_field, bz, ha, nz, nbx, nby, nbz, noplot=noplot ;+ ;NAME: ; lff_field.pro ;PURPOSE: ; to get the LFF field from the boundary longitudinal field. ;INPUTS: ; bz: the photospheric longitudinal field. ; nz: the height at which the field is calculated ; ha: the linear force-free factor ;OUTPUTS: ; nbx, nby: the calculated transverse field ; nbz: the calculated longitudinal field ;KEYWORD PARAMETERS: ; noplot: if set, do not display ;OPTIONAL INPUTS: ; ;OPTIONAL OUTPUTS: ; ;HISTORY: ; 1998 June, first edition at BBSO ;- sz=size(bz) & nx=sz(1) & ny=sz(2) nbx=fltarr(nx, ny) & nby=nbx i=complex(0,1) sub=indgen(nx, ny) n=sub/nx m=(sub-n*nx) kx=2*!pi*(m)/float(nx) ky=2*!pi*(n)/float(ny) kk=sqrt((sin(ky))^2+(sin(kx))^2)+0.0001 ;IFFT to get Bk first bk=fft(bz, -1) bk(0,0)=0. ;hlp,float(bk) ;window,1,r=2 ;tvscl, float(bk) ikk=i/(kk^2) kh=sqrt(abs(kk^2-ha^2)) ;now to get the Fourier components of the field components. fx=ikk*(ha*sin(ky)-kh*sin(kx))*exp(-kh*nz)*(abs(kk) ge (ha)) $ +ikk*(-kh*kx*sin(kh*nz)+ha*ky*cos(kh*nz))*(abs(ha) gt abs(kk)) fy=ikk*(-ha*sin(kx)-kh*sin(ky))*exp(-kh*nz)*(abs(kk) ge (ha)) $ +ikk*(-kh*ky*sin(kh*nz)-ha*kx*cos(kh*nz))*(abs(ha) gt abs(kk)) fz=exp(-kh*nz)*(abs(kk) ge abs(ha)) $ +cos(kh*nz)*(abs(ha) gt abs(kk)) ;FFT to get the field components at the height of nz. nbx=float(fft(fx*bk,1)) nby=float(fft(fy*bk,1)) nbz=float(fft(fz*bk,1)) ;hlp,nbx, nby, nbz ;let's check a bit? if not keyword_set(noplot) then begin xwin=800. & ywin=800. window, 0, xsize=xwin, ysize=ywin, r=2 key=(xwin/nx ge ywin/ny) zm=key*ywin/ny+(1-key)*xwin/nx print, 'Max and Min of Bz: ' hlp,nbz p_lev=fltarr(5) read,'Give positive contour levels(5) of Bz: ', p_lev p_lev=p_lev(sort(p_lev)) n_lev=-p_lev(sort(-p_lev)) contour, nbz, pos=[0,0,nx*zm, ny*zm],level=p_lev, /dev, xst=1, yst=1 contour, nbz, pos=[0,0,nx*zm, ny*zm],level=n_lev, /dev, xst=1, yst=1,/noera, c_lin=1 contour, nbz, pos=[0,0,nx*zm, ny*zm],level=[0], /dev, xst=1, yst=1,/noera, thick=2 vect, nbx, nby, 1, 40, 40, pos=[0,0,nx*zm, ny*zm],/dev,xst=1, yst=1, /noera endif end ;********************************************************************** ;------------------------------------------------------------- ;+ ; 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 THE IT IS A SMOOTH CURVE AND THE MAXIMUM IS THE GLOBALMAXIMUM, THEN ; XMAX AND YMAX OBTAINED ABOVE ARE EXACT; OTHERWISE, ???? @QJ 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 endif rmax=[xmax,ymax] finishup: return,rmax end ;********************************************************************** ;------------------------------------------------------------- ;+ ; NAME: ; XYOFF1 ; PURPOSE: ; This function, given two images, calculates the ; correlation function between those 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, [wdow,mask,hx,hy,fna]) ; INPUTS: ; a = The reference image. ; type: array,any type,arr(nx,ny) ; b = The image that you want to track with ; respect to a. ; type: array,any type,arr(nx,ny) ; bx = The size of the box in the x-dimension that you wish ; to use to generate the cross correlation function. ; type: scalar,integer ; by = The size of the box in the y-dimension that you wish ; to use to generate the cross correlation function. ; type: scalar,integer ; wdow = A 2-d apodization window. (Optional) ; type: array,any type,arr(nx,ny) ; mask = A 2-d mask to be applied in fourier space. (Optional) ; type: array,any type,arr(nx,ny) ; hx = The x position in the image that corresponds to the center ; of the box. (Optional) ; type: scalar,integer ; hy = The y position in the image that corresponds to the center ; of the box. (Optional) ; type: scalar,integer ; fna = The fourier transform of the 256x256 (Optional) ; subtended image from the center of a. ; type: array,floating point,fltarr(nx,ny) ; KEYWORD PARAMETERS: ; /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 opetrator) ; OUTPUTS: ; s = A vector containing the x and y offsets that ; should be applied to an image to maximize ; the cross correlation. ; type: vector,floating point,fltarr(2) ; s(0) = The x offset. ; type: scalar,floating point ; s(1) = The y offset. ; type: scalar,floating point ; COMMON BLOCKS: ; NOTES: ; The optional input parameters allow the user to ; speed up long jobs that, for instance, use ; the same reference. Simply add the 3 variable names ; hx, hy and fna to the calling of xyoff.pro in your driver ; routine. The first time it is called, those 3 variables ; won't be defined...so, xyoff will define them, ; the next time, xyoff, will see that they are defined ; and will not define them again, therefore removing ; that part of the algorithm and speeding up the routine. ; (If you don't want either a mask or a window, simply ; define them as 1.) ; ; If you wan't to change the reference some time ; later you must remove those 3 variables from the ; calling of xyoff.pro or make sure that the 3 ; variables are undefined. ; 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. ; Qiu J. 17 Jun, 1998 --- Just want to kick off the wdow mask stuff. ;- ;------------------------------------------------------------- function xyoff1,a,b,bx,by,hx, hy, wdow,mask,fna,$ sobel=sobel,help=help ;Display idl header if help is required. if keyword_set(help) or n_params() lt 4 then begin get_idlhdr,'xyoff1.pro' s=-1 goto,finishup endif ;Set correlation box size. mx=bx/2 & my=by/2 ;Test to see if window is specified. if n_elements(wdow) eq 0 then wdow=1. if n_elements(mask) eq 0 then mask=1. ;Calculate size of image. ssz=size(a) ;Determine whether hx exists if not, then define it. if n_elements(hx) eq 0 then begin nx=ssz(1) hx=nx/2 endif ;Determine whether hy exists if not, then define it. if n_elements(hy) eq 0 then begin ny=ssz(2) hy=ny/2 endif ;Speedup section. if n_elements(fna) eq 0 then begin na=a(hx-mx:hx+mx-1,hy-my:hy+my-1) if n_elements(sobel) eq 0 then na=(na-mean(na))*wdow $ else na=sobel(na-mean(na)) fna=fft(na*mask,-1) endif ;Subset correltion box around the center of image. nb=b(hx-mx:hx+mx-1,hy-my:hy+my-1) if n_elements(sobel) eq 0 then nb=(nb-mean(nb))*wdow $ else nb=sobel(nb-mean(nb)) ;Compute the fourier transform of nb. fnb=fft(nb*mask,-1) ;Compute the correlation function between na and nb. ccf=fft(fna*conj(fnb),1) ;Rename the cross correlation function as it's real part. ccf=float(ccf) ;Shift correlation function. ccf=shift(ccf,mx,my) ;Find maximum of correlation function. ps=max_pos(ccf) ;Calculate shifts. s=fltarr(2) s(0)=ps(0)-mx s(1)=ps(1)-my print,' Shift in X - ',s(0) print,' Shift in Y - ',s(1) finishup: return,s end ;*************************************************************** pro fit_linear, x, y, alpha, beta, sigma, s_alpha, s_beta, plot = plot, error_plot = error_plot ;+ ;PURPOSE: ; least-square linear fit of y = alpha*x+beta ;INPUT: ; x: array of independent variable ; y: array of dependent variable ;OUTPUT: ; alpha: parameter from the fit ; beta: parameter from the fit ; sigma: standard deviation of the fit ; s_alpha: fitting error in alpha ; s _beta: fitting error in beta ;KEYWORDS: ; plot: if set, plot the fit result; ; error_plot: if set, add the error bars in the fit plot. ;HISTORY: ; 1998: first version; QiuJ. ; 2000: add keyword error_plot; QiuJ ;- sumx=0. sumy=0. sumxx=0. sumxy=0. sumn=0. sumyy=0. sigma = 0. alpha = 0. beta = 0. nx=n_elements(x) ny=n_elements(y) if nx ne ny then stop else sumn=nx sumx=total(x) sumy=total(y) sumxx=total(x*x) sumyy=total(y*y) sumxy=total(x*y) determinant=sumxx*sumn-sumx*sumx if (abs(determinant) ge 0.0001) then begin alpha=(sumxy*sumn-sumx*sumy) / determinant beta=(sumxx*sumy-sumx*sumxy) / determinant sigma=total((y-alpha*x-beta)*(y-alpha*x-beta)) sigma=sqrt(sigma/(sumn-2)) s_alpha = sqrt(sumn*sigma^2/determinant) s_beta = sqrt(sumxx*sigma^2/determinant) ;print, 'y = alpha * x + beta' ;print, 'Alpha, beta, sigma, s_alpha, s_beta' ;print, Alpha, beta, sigma, s_alpha, s_beta if keyword_set(plot) then begin plot, x, y, psym = 1 oplot, x, x*alpha+beta if keyword_set(error_plot) then for i = 0, nx - 1 do plots, $ [x(i), x(i)], [y(i)-sigma, y(i)+sigma], /data endif endif else message, 'Cannot fit!' end