; docformat = 'rst' ; ; NAME: ; shapeRmsImg ; PURPOSE: ; Compute the RMS deviation between a shape model and imaging 2-D profile ;+ ; :Description: ; Compute the RMS deviation between a shape model and imaging 2-D profile ; ; :Categories: ; Shape Models ; ; :Returns: The (model-profile) RMS ; ; INPUTS: ; XY: in, required, type=fltarr ; Cartesian coordinates of the profile ; img: in, required, type=structure ; The structure of the imaging data (single epoch). See ; koalaReadAO for more details on the fields. ; ; :Keywords: ; ; :Examples: ; Read a model profile, and an imaging data, and compute the rms:: ; IDL> img = koalaReadAO( 'my_AO_file' ) ; IDL> model = readVer( 'my_model_file' ) ; IDL> spin = readSpin( 'my_spin_file' ) ; IDL> ssp = subObserver_coord( img[0].sun, *spin[0],'iso_date_of_occultation', PA=PA, /LTC ) ; IDL> sep = subObserver_coord( img[0].obs, *spin[0],'iso_date_of_occultation', PA=PA, /LTC ) ; IDL> profile = shapeProfile( *model[0], sep=sep, ssp=ssp, pa=pa ) ; IDL> rms = shapeRmsImg( profile, img[0] ) ; ; :Uses: ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in August 2016, B. Carry (OCA) ; 2017 Dec - B.Carry (OCA) - Corrected bug in resampling ;- function shapeRmsImg, xy, img COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Check Input Parameters -------------------------------------------------------------- if N_PARAMS() LT 2 then begin message, /ioError, 'Syntax : rms = shapeRmsImg( xy, img )' return, -1 endif ;--I.2-- Profile Structure ------------------------------------------------------------------- if not keyword_set(img.polar) then begin endif ;--I.x-- Useful Constants -------------------------------------------------------------------- AUtoKM = 149597870.691D ;-km ARCtoRAD = !DTOR/3600. ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Check and Homogeneise Profile Sampling -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; resample=0 ;--II.1-- Profiles Number of Points ---------------------------------------------------------- nbImg = n_elements(*img.polar) nbMod = n_elements(xy)/2. if nbMod ne nbImg then resample=1 ;--II.2-- Convert Cartesian Profile into Polar ----------------------------------------------- pol = cv_coord( from_rect=xy, /to_polar ) theta = (*img.polar).ang mod !PI ;--II.3-- Further Azimuth Test --------------------------------------------------------------- if resample eq 0 then begin ;--II.3.1-- Sort Model Azimuth sortMod = sort( pol[0,*] ) ordMod = reform( pol[0,sortMod] ) ;--II.3.2-- Sort Image Azimuth sortImg = sort( (*img.polar).ang ) ordImg = theta[sortImg] ;--II.3.3-- Resample if Difference is notable totAngDiff = total( (ordMod-ordImg)^2 )/nbMod if totAngDiff ge 1e-5 then resample=1 endif ;--II.4-- Resample Model if Needed ----------------------------------------------------------- if resample ne 0 then begin ;--II.4.1-- Replicate Contour for Interpolation lonRad = [ reform(pol[1,sortMod]), reform(pol[1,sortMod]), reform(pol[1,sortMod]) ] lonAng = [ reform(pol[0,sortMod])-2*!PI, reform(pol[0,sortMod]), reform(pol[0,sortMod])+2*!PI] ;--II.4.2-- Contour Smoothing modKM = interpol( lonRad, lonAng, theta ) endif else modKM = reform(pol[1,*]) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Image and Model Unit Conversion -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Range to Observer ----------------------------------------------------------------- range = sqrt( total( img.obs^2 ) ) * AUtoKM ;--III.2-- Model Conversion ------------------------------------------------------------------ modAS = atan( modKM / range ) / ARCtoRAD ;-Model radius in arcsec modPIX = modAS / img.scale ;-Model radius in pixel ;--III.3-- Image Conversion ------------------------------------------------------------------ imgPIX = (*img.polar).rad ;-Image radius in pixel imgAS = (*img.polar).rad * img.scale ;-Image radius in arcsec imgKM = range* tan( imgAS* arctoRAD) ;-Image radius in KM ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Compute RMS and Return -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Model-Image Differences ----------------------------------------------------------- rmsKM = sqrt( mean( (modKM-imgKM)^2 ) ) rmsAS = sqrt( mean( (modAS-imgAS)^2 ) ) rmsPIX = sqrt( mean( (modPIX-imgPIX)^2 ) ) rmsSIG = sqrt( mean( ((modPIX-imgPIX)/ (*img.polar).unc )^2 ) ) ;--III.2-- Compute Root Mean Square ---------------------------------------------------------- rms = [rmsPIX, rmsAS, rmsKM, rmsSIG] return, rms end