pro undistort, fitsfile ;Will read in the fits file and create an undistorted version. The ;algorithm uses nearest neighbor sampling. ;Note that these coeffs are based on a pixel origin that is not 1,1 ;(like IRAF). The center of the origin pixel is actually 0.5,0.5 (so ;that the lower left corder of the origin pixel is 0,0). image = mrdfits(fitsfile, 0, header, /dscale, /silent) xsize = n_elements(image[*,0]) ysize = n_elements(image[0,*]) x = dindgen(xsize) + 0.5d x = x#transpose(replicate(1d,ysize)) y = dindgen(ysize) + 0.5d y = replicate(1d,xsize)#y ;undistort Altair/NIRI images p = dblarr(7) p[0] = 512d p[1] = 512d p[2] = 0d p[3] = 1d p[4] = 1.32d-5 ;readcol, 'coeffs2', co1, dummy, p, perr, format='a,a,d,d' ;coeff file should contain coeffs 4..14: ;4: x center of distortion ;5: y center of distortion ;6: zeroth order offset (defined as zero) ;7-14: first..eighth order xp = x - p[0] yp = y - p[1] r = sqrt(xp*xp + yp*yp) theta = atan(yp, xp) rp = 0d ;correct for distortion ndegree = n_elements(p)-3 for i=0l,ndegree do begin rp = rp + p[i+2]*(r^i) endfor ;convert back into cartesian coordinates xp = rp*cos(theta) yp = rp*sin(theta) xp = xp + p[0] yp = yp + p[1] ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Do a linear solution print, "Doing a bilinear solution" print, "running triangulate..." triangulate, xp, yp, tri gs = [1d,1d] limits = [0d,0d,double(xsize-1),double(ysize-1)] + 0.5d print, "computing solution..." result = trigrid(xp, yp, image, tri, gs, limits) help, result ; write the new file newfile = 'l' + fitsfile mwrfits, result, newfile, header, /create ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Do a quintic solution print, "Doing a quintic solution" print, "running triangulate..." triangulate, xp, yp, tri gs = [1d,1d] limits = [0d,0d,double(xsize-1),double(ysize-1)] + 0.5d print, "computing solution..." result = trigrid(xp, yp, image, tri, gs, limits, /quintic) help, result ; write the new file newfile = 'q' + fitsfile mwrfits, result, newfile, header, /create ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Do a nearest neighbor solution print, "Doing a nearest neighbor solution" newimage = image - image ;for each of the pixels in the new image, find the pixel in the old ;image that is closest oldj = -1 maxdev = 10l ;the maximum distortion neccessary (could compute this in real time) for j=0l,ysize-1 do begin newj = j/10 if(newj ne oldj) then begin ;print, "j: ", j oldj = newj endif for i=0l,xsize-1 do begin ;don't need to compute the whole image, just a subset xmin = (i - maxdev) > 0 xmax = (i + maxdev) < xsize - 1 ymin = (j - maxdev) > 0 ymax = (j + maxdev) < ysize - 1 sxp = xp[xmin:xmax,ymin:ymax] syp = yp[xmin:xmax,ymin:ymax] simage = image[xmin:xmax,ymin:ymax] xdev = x[i,j] - sxp ydev = y[i,j] - syp dev = sqrt(xdev*xdev + ydev*ydev) m = min(dev, index) newimage[i,j] = simage[index] endfor endfor ;write out the new file newfile = 'n' + fitsfile mwrfits, newimage, newfile, header, /create end