;+ ; $Id: scc_normalize.pro,v 1.9 2013/06/05 17:47:32 mcnutt Exp $ ; ; Project : STEREO - SECCHI ; ; Name : SCC_normalize ; ; Purpose : Create standard normalized image for A-B comparison ; ; Category : STEREO, SECCHI, Registration, display, movie ; ; Explanation : This procedure creates an image with standard alignment from one ; STEREO spacecraft. An image is read in, processed, and rescaled to ; a common center and image size. The orientation of the image is ; aligned for proper stereo viewing. ; ; Syntax : Output = SCC_normalize( filename ) ; ; or ; ; Output = SCC_normalize( image, header ) ; ; Examples : file = '20070321_122811_s4euA.fts' ; Output = SCC_normalize(file) ; EXPTV, Output, TRUE=1 ; ; Inputs : filename ; ; Opt. Inputs : Instead of passing the filename, one can pass in image and ; index which have already been read in. In that case, the ; calling parameters are: ; ; image = Image(s) for either observatory ; header = FITS index structure(s) for image ; ; To return corrected wcshdr then input image = -1 wcs = SCC_normalize( -1, bheader, ref_hdr=aheader ,/wcscorr) ; ; Outputs : The result of the function is an image with ; dimensions [Nx,Ny] where Nx, Ny are the dimensions of the image ; ; Opt. Outputs: INDEX = FITS index structures for the image ; ; Keywords : INTERP = Interpolation type, 0=nearest neighbor, 1=bilinear ; interpolation (default), 2=cubic interpolation. ; ; CUBIC = Interpolation parameter for cubic interpolation. ; See the IDL documentation for POLY_2D for more ; information. ; ; NOMISSING = If set, then don't filter out missing pixels around ; the edge of the CCD. ; ; SECCHIPREP= If set, call SECCHI_PREP ; ; REF_HDR= Header to use as reference ; ; ROTATE_BY= Specify how much to rot image (degrees) instead of to STEREO mission plane ; ; WCSCORR = return corrected wcs header instead of image ; ; /NOROTATE Do not correct roll angle ; ; Calls : SCCREADFITS, FITSHEAD2WCS, WCS_GET_PIXEL, GET_STEREO_ROLL, ; WCS2FITSHEAD, SECCHI_PREP ; ; Common : None. ; ; Restrictions: Currently only works for STEREO COR1, COR2 and EUVI images. Does ; not (yet) work for combinations of STEREO and other ; observatories (e.g. SOHO). ; ; Side effects: The returned index structures are updated to reflect the ; image manipulations. CRVALi is set to zero and CRPIX becomes the sun center. ; ; Prev. Hist. : This is the guts of scc_stereopair.pro ; ; Written : 29-Aug-2006, William Thompson, GSFC ; ; $Log: scc_normalize.pro,v $ ; Revision 1.9 2013/06/05 17:47:32 mcnutt ; added adjustment to stereo_roll system STPLN need after 2011-02-06 when spacecraft seperation reasched 180 ; ; Revision 1.8 2013/05/22 17:07:08 nathan ; do not stop of wcs.time.observ_avg not present ; ; Revision 1.7 2011/12/22 20:40:07 nathan ; add /norotate and allow non-STEREO input ; ; Revision 1.6 2011/02/17 21:21:33 nathan ; do not change pixel0 or dsun if USE_ROLL set ; ; Revision 1.5 2010/12/09 22:30:29 nathan ; pixel0+1 for crpix; move poly_2d ; ; Revision 1.4 2010/01/15 20:09:52 mcnutt ; added keyword wcscorr to return corrected wcs hdr only ; ; Revision 1.3 2008/06/18 18:23:09 nathan ; added ROTATE_BY= ; ; Revision 1.2 2008/04/04 18:53:46 nathan ; omit MISSING keyword for poly_2d because getting NaN for mask ; ; Revision 1.1 2008/04/02 16:50:52 nathan ; standardized version of scc_stereopair.pro for scc_mkframe.pro ; ; Revision 1.6 2007/09/04 19:19:37 thompson ; Fixed bug ; ; Revision 1.5 2007/08/03 22:25:54 thompson ; Added REVERSE keyword. Treat secondary coordinate system. ; ; Revision 1.4 2007/08/03 19:55:36 thompson ; Corrected bug with wrong roll angle being returned in header structures. ; Take telescope roll relative to spacecraft frame into account ; ; Revision 1.3 2007/07/02 22:32:51 thompson ; Add SECCHIPREP keyword, option of passing in images and headers ; ; Revision 1.2 2007/06/28 22:09:03 thompson ; *** empty log message *** ; ; Revision 1.1 2006/09/27 17:04:06 nathan ; moved from dev/display ; ; Revision 1.5 2006/09/26 21:45:18 nathan ; add CVS keywords, Written label ; ; Version 3, 25-Sep-2006, William Thompson, GSFC ; Split INDEX into INDEX_BEHIND and INDEX_AHEAD ; Version 2, 19-Sep-2006, William Thompson, GSFC ; Added output parameter INDEX ; ; Contact : WTHOMPSON ;- ; function scc_normalize, files_ahead, index_ahead, interp=interp, REF_HDR=ref_hdr, $ secchiprep=secchiprep, MISSING=missing, ROTATE_BY=rotate_by, $ WCSCORR=WCSCORR, debug=debug, NOROTATE=norotate, _extra=_extra ; ; Check whether the filename or image array is passed. ; if keyword_set(WCSCORR) then begin n_files=0 header=index_ahead goto,wcsonly endif names_passed = datatype(files_ahead,1) eq 'String' if names_passed then begin if n_params() lt 1 then message, $ 'Syntax: Result = SCC_normalize(FILEname) end else begin if n_params() lt 2 then message, $ 'Syntax: Result = SCC_normalize(IMAGE, INDEX) sza = size(files_ahead) case sza[0] of 2: n_files = 1 3: n_files = sza[3] else: message, 'IMAGE_AHEAD must have 2 or 3 dimensions' endcase endelse ; ; Get the INTERP value to pass to POLY_2D. ; if n_elements(interp) eq 0 then interp = 1 ; ; Read the image. ; if names_passed then begin if keyword_set(secchiprep) then begin secchi_prep, files_ahead, header_ahead, image_ahead, $ _extra=_extra end else begin image_ahead = sccreadfits(files_ahead[*], header_ahead) endelse header=header_ahead endif ELSE header=index_ahead ; if names_passed then image = image_ahead else image = files_ahead ; ; Create the output array. ; if names_passed then sz = size(image_ahead) else sz = size(files_ahead) type = sz[sz[0]+1] if keyword_set(anaglyph) then dim = [3,sz[1:2]] else dim = [2,sz[1:2]] if n_files gt 0 then dim = [dim, n_files] output = make_array(type=type, dimension=dim) ; ; Select a missing value, based on the datatype. ; IF ~keyword_set(MISSING) THEN $ case type of 4: missing = !values.f_nan 5: missing = !values.d_nan else: missing = 0 endcase wcsonly: ; Get the solar distance, the pixel spacing, and the pixel location of Sun ; center. All the other images will be matched to this. ; IF keyword_set(REF_HDR) THEN BEGIN wcs0 = fitshead2wcs(ref_hdr) dsun_obs0 = wcs0.position.dsun_obs cdelt0 = wcs0.cdelt[0] pixel0 = wcs_get_pixel(wcs0, [0,0]) crota0 = wcs0.roll_angle ENDIF ELSE BEGIN dsun_obs0 = header.dsun_obs ;oneau('KM')*1000. cdelt0 = header.cdelt1 ;getsccsecpix(header,'a') pixel0 = [header.naxis1/2.-0.5, header.naxis2/2.-0.5] ; center of image crota0 = 0 ENDELSE ; ; ; Resample each image so that each has the same Sun center, scale. Rotate ; the image by the roll angle to the STEREO Mission Plane coordinate system. ; p = dblarr(2,2) q = dblarr(2,2) wcs = fitshead2wcs(header) ;if ~tag_exist(wcs.time,'observ_avg') then stop IF keyword_set(NOROTATE) THEN roll = 0 ELSE $ IF keyword_set(ROTATE_BY) THEN BEGIN roll = rotate_by ;pixel0=[header.crpix1,header.crpix2]-1 ;dsun_obs0=header.dsun_obs ENDIF ELSE IF strmid(header.obsrvtry,0,3) EQ 'STE' THEN BEGIN stpln=get_stereo_roll(wcs.time.observ_date, header.obsrvtry, system='STPLN') stplnadj=0.0 if stpln lt -90 then stplnadj=180. if stpln gt 90 then stplnadj=(-180.) roll = (stpln+stplnadj) - get_stereo_roll(wcs.time.observ_date, header.obsrvtry, system='RTN') - wcs.roll_angle ENDIF ELSE $ roll = crota0-wcs.roll_angle message,header.filename+' roll change of '+trim(roll)+' deg',/info cosr = cos(roll * !dpi / 180.d0) sinr = sin(roll * !dpi / 180.d0) pixel = wcs_get_pixel(wcs, [0,0]) scale = (cdelt0 / wcs.cdelt) * (dsun_obs0 / wcs.position.dsun_obs) ; ; Modify the WCS structure to reflect the new image scale. ; wcs.crpix = pixel0+1 wcs.crval[*] = 0 wcs.cdelt = wcs.cdelt * scale wcs.roll_angle = wcs.roll_angle + roll cos_a = cos(wcs.roll_angle * !dpi / 180.d0) sin_a = sin(wcs.roll_angle * !dpi / 180.d0) lambda = wcs.cdelt[1] / wcs.cdelt[0] wcs.pc[0,0] = cos_a wcs.pc[0,1] = -lambda * sin_a wcs.pc[1,0] = sin_a / lambda wcs.pc[1,1] = cos_a if ~keyword_set(WCSCORR) then begin ; skip secondary coordinate system and return corrected wcshdr used in scc_playmovie call when defing wcshdrs for ab movies ; ; Call POLY_2D to resample the image. ; p[0,0] = pixel[0] - scale[0]*pixel0[0]*cosr + scale[1]*pixel0[1]*sinr q[0,0] = pixel[1] - scale[0]*pixel0[0]*sinr - scale[1]*pixel0[1]*cosr p[0,1] = scale[0]*cosr p[1,0] = -scale[1]*sinr q[0,1] = scale[1]*sinr q[1,0] = scale[0]*cosr if n_files gt 0 then output = poly_2d(image, p, q, interp, _extra=_extra) iindex = wcs2fitshead(wcs, old=header, /add_xcen, /add_roll, $ /structure) temp = index_ahead struct_assign, iindex, temp ; ; Update the secondary coordinate system. ; wcsa = fitshead2wcs(header, system='a') origina = wcs_get_coord(wcsa, pixel) temp.crpix1a = temp.crpix1 temp.crpix2a = temp.crpix2 temp.crval1a = origina[0] temp.crval2a = origina[1] temp.cdelt1a = temp.cdelt1a * scale[0] temp.cdelt2a = temp.cdelt2a * scale[1] rolla = wcsa.roll_angle - roll cos_a = cos(rolla * !dpi / 180.d0) sin_a = sin(rolla * !dpi / 180.d0) temp.pc1_1a = cos_a temp.pc1_2a = -lambda * sin_a temp.pc2_1a = sin_a / lambda temp.pc2_2a = cos_a ; ; Store the modified header. ; index_ahead = temporary(temp) endif else output=wcs ; ; return, output end