diff options
Diffstat (limited to 'lib/stwcs/distortion/models.py')
| -rw-r--r-- | lib/stwcs/distortion/models.py | 361 |
1 files changed, 361 insertions, 0 deletions
diff --git a/lib/stwcs/distortion/models.py b/lib/stwcs/distortion/models.py new file mode 100644 index 0000000..96d5d72 --- /dev/null +++ b/lib/stwcs/distortion/models.py @@ -0,0 +1,361 @@ +from __future__ import division # confidence high + +import types +# Import PyDrizzle utility modules +import mutil +import numpy as np +import mutil +from mutil import combin + +yes = True +no = False + +################# +# +# +# Geometry/Distortion Classes +# +# +################# + +class GeometryModel: + """ + Base class for Distortion model. + There will be a separate class for each type of + model/filetype used with drizzle, i.e., IDCModel and + DrizzleModel. + + Each class will know how to apply the distortion to a + single point and how to convert coefficients to an input table + suitable for the drizzle task. + + Coefficients will be stored in CX,CY arrays. + """ + # + # + # + # + # + # + # + NORDER = 3 + + def __init__(self): + " This will open the given file and determine its type and norder." + + # Method to read in coefficients from given table and + # populate the n arrays 'cx' and 'cy'. + # This will be different for each type of input file, + # IDCTAB vs. drizzle table. + + # Set these up here for all sub-classes to use... + # But, calculate norder and cx,cy arrays in detector specific classes. + self.cx = None + self.cy = None + self.refpix = None + self.norder = self.NORDER + # Keep track of computed zero-point for distortion coeffs + self.x0 = None + self.y0 = None + + # default values for these attributes + self.direction = 'forward' + + self.pscale = 1.0 + + def shift(self,cx,cy,xs,ys): + """ + Shift reference position of coefficients to new center + where (xs,ys) = old-reference-position - subarray/image center. + This will support creating coeffs files for drizzle which will + be applied relative to the center of the image, rather than relative + to the reference position of the chip. + """ + + _cxs = np.zeros(shape=cx.shape,dtype=cx.dtype) + _cys = np.zeros(shape=cy.shape,dtype=cy.dtype) + _k = self.norder + 1 + # loop over each input coefficient + for m in xrange(_k): + for n in xrange(_k): + if m >= n: + # For this coefficient, shift by xs/ys. + _ilist = range(m, _k) + # sum from m to k + for i in _ilist: + _jlist = range(n, i - (m-n)+1) + # sum from n to i-(m-n) + for j in _jlist: + _cxs[m,n] += cx[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n))) + _cys[m,n] += cy[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n))) + return _cxs,_cys + + def convert(self, tmpname, xref=None,yref=None,delta=yes): + """ + Open up an ASCII file, output coefficients in drizzle + format after converting them as necessary. + First, normalize these coefficients to what drizzle expects + Normalize the coefficients by the MODEL/output plate scale. + + 16-May-2002: + Revised to work with higher order polynomials by John Blakeslee. + 27-June-2002: + Added ability to shift coefficients to new center for support + of subarrays. + """ + cx = self.cx/self.pscale + cy = self.cy/self.pscale + x0 = self.refpix['XDELTA'] + cx[0,0] + y0 = self.refpix['YDELTA'] + cy[0,0] + #xr = self.refpix['XREF'] + #yr = self.refpix['YREF'] + xr = self.refpix['CHIP_XREF'] + yr = self.refpix['CHIP_YREF'] + + + + ''' + if xref != None: + # Shift coefficients for use with drizzle + _xs = xref - self.refpix['XREF'] + 1.0 + _ys = yref - self.refpix['YREF'] + 1.0 + + + if _xs != 0 or _ys != 0: + cxs,cys= self.shift(cx, cy, _xs, _ys) + cx = cxs + cy = cys + + # We only want to apply this shift to coeffs + # for subarray images. + if delta == no: + cxs[0,0] = cxs[0,0] - _xs + cys[0,0] = cys[0,0] - _ys + + # Now, apply only the difference introduced by the distortion.. + # i.e., (undistorted - original) shift. + x0 += cxs[0,0] + y0 += cys[0,0] + ''' + self.x0 = x0 #+ 1.0 + self.y0 = y0 #+ 1.0 + + # Now, write out the coefficients into an ASCII + # file in 'drizzle' format. + lines = [] + + + lines.append('# Polynomial distortion coefficients\n') + lines.append('# Extracted from "%s" \n'%self.name) + lines.append('refpix %f %f \n'%(xr,yr)) + if self.norder==3: + lines.append('cubic\n') + elif self.norder==4: + lines.append('quartic\n') + elif self.norder==5: + lines.append('quintic\n') + else: + raise ValueError, "Drizzle cannot handle poly distortions of order %d"%self.norder + + str = '%16.8f %16.8g %16.8g %16.8g %16.8g \n'% (x0,cx[1,1],cx[1,0],cx[2,2],cx[2,1]) + lines.append(str) + str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[2,0],cx[3,3],cx[3,2],cx[3,1],cx[3,0]) + lines.append(str) + if self.norder>3: + str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[4,4],cx[4,3],cx[4,2],cx[4,1],cx[4,0]) + lines.append(str) + if self.norder>4: + str = '%16.8g %16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[5,5],cx[5,4],cx[5,3],cx[5,2],cx[5,1],cx[5,0]) + lines.append(str) + lines.append("\n") + + str = '%16.8f %16.8g %16.8g %16.8g %16.8g \n'% (y0,cy[1,1],cy[1,0],cy[2,2],cy[2,1]) + lines.append(str) + str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[2,0],cy[3,3],cy[3,2],cy[3,1],cy[3,0]) + lines.append(str) + if self.norder>3: + str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[4,4],cy[4,3],cy[4,2],cy[4,1],cy[4,0]) + lines.append(str) + if self.norder>4: + str = '%16.8g %16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[5,5],cy[5,4],cy[5,3],cy[5,2],cy[5,1],cy[5,0]) + lines.append(str) + + output = open(tmpname,'w') + output.writelines(lines) + output.close() + + + def apply(self, pixpos,scale=1.0,order=None): + """ + Apply coefficients to a pixel position or a list of positions. + This should be the same for all coefficients tables. + Return the geometrically-adjusted position + in arcseconds from the reference position as a tuple (x,y). + + Compute delta from reference position + """ + + """ + scale actually is a ratio of pscale/self.model.pscale + what is pscale? + """ + if self.cx == None: + return pixpos[:,0],pixpos[:,1] + + if order is None: + order = self.norder + + # Apply in the same way that 'drizzle' would... + _cx = self.cx / (self.pscale * scale) + _cy = self.cy / (self.pscale * scale) + _convert = no + _p = pixpos + + # Do NOT include any zero-point terms in CX,CY here + # as they should not be scaled by plate-scale like rest + # of coeffs... This makes the computations consistent + # with 'drizzle'. WJH 17-Feb-2004 + _cx[0,0] = 0. + _cy[0,0] = 0. + + if isinstance(_p,types.ListType) or isinstance(_p,types.TupleType): + _p = np.array(_p,dtype=np.float64) + _convert = yes + + dxy = _p - (self.refpix['XREF'],self.refpix['YREF']) + # Apply coefficients from distortion model here... + c = _p * 0. + for i in range(order+1): + for j in range(i+1): + c[:,0] = c[:,0] + _cx[i][j] * pow(dxy[:,0],j) * pow(dxy[:,1],(i-j)) + c[:,1] = c[:,1] + _cy[i][j] * pow(dxy[:,0],j) * pow(dxy[:,1],(i-j)) + xc = c[:,0] + yc = c[:,1] + + # Convert results back to same form as original input + if _convert: + xc = xc.tolist() + yc = yc.tolist() + # If a single tuple was input, return just a single tuple + if len(xc) == 1: + xc = xc[0] + yc = yc[0] + + return xc,yc + + def setPScaleCoeffs(self,pscale): + self.cx[1,1] = pscale + self.cy[1,0] = pscale + + self.refpix['PSCALE'] = pscale + self.pscale = pscale + + +class IDCModel(GeometryModel): + """ + This class will open the IDCTAB, select proper row based on + chip/direction and populate cx,cy arrays. + We also need to read in SCALE, XCOM,YCOM, XREF,YREF as well. + """ + def __init__(self, idcfile, date=None, chip=1, direction='forward', + filter1='CLEAR1',filter2='CLEAR2',offtab=None, binned=1): + GeometryModel.__init__(self) + # + # Norder must be derived from the coeffs file itself, + # then the arrays can be setup. Thus, it needs to be + # done in the sub-class, not in the base class. + # Read in table. + # Populate cx,cy,scale, and other variables here. + # + self.name = idcfile + self.cx,self.cy,self.refpix,self.norder = mutil.readIDCtab(idcfile, + chip=chip,direction=direction,filter1=filter1,filter2=filter2, + date=date, offtab=offtab) + + if self.refpix.has_key('empty_model') and self.refpix['empty_model']: + pass + else: + self.refpix['PSCALE'] = self.refpix['PSCALE'] * binned + self.cx = self.cx * binned + self.cy = self.cy * binned + self.refpix['XREF'] = self.refpix['XREF'] / binned + self.refpix['YREF'] = self.refpix['YREF'] / binned + self.refpix['XSIZE'] = self.refpix['XSIZE'] / binned + self.refpix['YSIZE'] = self.refpix['YSIZE'] / binned + + self.pscale = self.refpix['PSCALE'] + + +class WCSModel(GeometryModel): + """ + This class sets up a distortion model based on coefficients + found in the image header. + """ + def __init__(self,header,rootname): + GeometryModel.__init__(self) + + + if header.has_key('rootname'): + self.name = header['rootname'] + else: + self.name = rootname + # Initialize all necessary distortion arrays with + # default model... + #self.cx,self.cy,self.refpix,self.order = mutil.defaultModel() + + # Read in values from header, and update distortion arrays. + self.cx,self.cy,self.refpix,self.norder = mutil.readWCSCoeffs(header) + + self.pscale = self.refpix['PSCALE'] + + + +class DrizzleModel(GeometryModel): + """ + This class will read in an ASCII Cubic + drizzle coeffs file and populate the cx,cy arrays. + """ + + def __init__(self, idcfile, scale = None): + GeometryModel.__init__(self) + # + # We now need to read in the file, populate cx,cy, and + # other variables as necessary. + # + self.name = idcfile + self.cx,self.cy,self.refpix,self.norder = mutil.readCubicTable(idcfile) + + # scale is the ratio wcs.pscale/model.pscale. + # model.pscale for WFPC2 is passed from REFDATA. + # This is needed for WFPC2 binned data. + + if scale != None: + self.pscale = scale + else: + self.pscale = self.refpix['PSCALE'] + + """ + The above definition looks wrong. + In one case it's a ratio in the other it's pscale. + + """ + +class TraugerModel(GeometryModel): + """ + This class will read in the ASCII Trauger coeffs + file, convert them to SIAF coefficients, then populate + the cx,cy arrays. + """ + NORDER = 3 + + def __init__(self, idcfile,lam): + GeometryModel.__init__(self) + self.name = idcfile + self.cx,self.cy,self.refpix,self.norder = mutil.readTraugerTable(idcfile,lam) + self.pscale = self.refpix['PSCALE'] + # + # Read in file here. + # Populate cx,cy, and other variables. + # + + |