M. Aurand4.21.2010EE552FFT-Based Astronomical Image Registration and Stacking using GPUThe productive imaging of faint astronomical targets mandates vanishingly low noise due to the small amount of signal amplitude available. The readout noise of the camera is gaussian in nature and dominant in presence; hence a popular method of reducing its influence is image stacking. Several images can be summed to minimize this white noise as a function of 1/√N. This project explores the use of MATLAB in combination with GPUMAT to create a functional program to register and stack a series of images to reduce total noise contributions.IntroductionThe images generated by the telescopes I use are typically a few megapixels in size, at 16-bit intensity depth. Each image in a series is dithered by a small amount to reduce the influence of pixel defects on the final, stacked image. Any program to process the image series must precisely identify the offsets and apply spatial transformations necessary to produce a high-resolution, low-noise final photograph.My program performs several steps to ensure adequate and reliable operation under real-world circumstances. The FFT cross-correlation is used to measure the offset of any image with respect to a user-selected “reference” image. FFT correlation is not very sensitive to sub-pixel movements, but works well as a coarse centering tool. Once the image has been roughly centered (within an integer pixel value), several stars are measured and compared to the “reference”, quantifying sub-pixel residual offsets. The rotation angle of the image is also measured to determine whether or not the image should be rotated. Scale from image to image is invariant, and needs no correction.Below is a sequential list of operations performed by my software.A. Read in reference image and image series to be registered/stackedB. Pre-process reference image1. Automatically generate valid list of reference stars2. Quantify position and angular moment of each reference starC. FFT-register (to integer pixel values)D. Centroid-register1. Quantify residual offset/angular moments relative to B.2E. Image transformation1. Generate transformation matrix2. Apply/interpolateF. Image summationA. Image openingThe images are opened using the MATLAB fitsread( ) function, which supports the 16-bit data format of scientific instruments. Alternatively, more popular image formats such as .jpg and .bmp can be imported by the imread( ) function. These images are allocated within GPU memory as single-precision floating points.B. Reference Image Pre-ProcessingAfter coarse FFT registration, star images will be used as tie points for fine adjustment of position and rotation. In order to maximize solution robustness, only the center 50% of the reference image area is utilized to generate tie stars. If one was to use stars near the edges of the reference frame, it is potentially unlikely that they will appear in the off-center images. This center 50% of the image implies that even with an input image shifted by M/2 or N/2, at least 25% of the centroid points should still be usable during fine adjustment.In order to automate the selection of reference stars, the usable center 50% is thresholded to generate a binary image. This binary image is sorted to identify and label the points. Next, some criteria are applied to select the best stars. In order to preserve the accuracy of a center-of-mass centroiding algorithm, the stars should not be saturated nor should they be too spatially close to an adjacent star. The stars identified in the binary image are sequentially examined and the ones that are saturated or within the exclusion radius (user-adjustable) of another star are discarded. This leaves a list of only the favorable tie-point stars.Once the list of valid tie stars has been generated, the center-of-mass centroiding algorithm is applied, and a reference feature list created. The reference feature list parametrizes each star's position and angle relative to the center of the reference image.C. FFT RegistrationFFT registration provides an exceedingly convenient way of obtaining a test image's Cartesian X and Y displacement relative to the reference image by means of correlation. Suppose we have two identical images (f1 & f2), one with a spatial offset of x0 and y0. The spatially offset image (f2) can be represented in terms of the other image:f2(x, y) = f1( x – x0 , y – y0 )Upon conversion to the frequency domain by Fourier transformation, a separable term A is created.F2(u, v) = F1(u , v) * e-2jπ ( ux0 + vy0 )A = e-2jπ ( ux0 + vy0 ) : F2(u, v) = F1(u , v) * ABy taking the cross-spectrum of the two transforms, F1 and F2, A is readily solved.A = F1(u , v) F2* (u, v) / | F1(u, v) F2(u, v) |Where A is visualized as an impulse in the cross-spectrum. It is located at the offset: u, v = (x0, y0 )In order to obtain this correlation numerically, the GPU functions fft2( ) and ifft2( ) are employed to produce a 2D result array. This result array features an impulse centered about the measured X and Y offset of the test image, which is recorded as a rough offset. This is a coarse result (integer pixel value) because the FFT is critically sampled. Better resolution could be obtained by padding the input images to larger size, at the expense of memory and computation time. In this case, however, any extra resolution is of no utility because the centroiding method is next employed. The centroiding algorithm is significantly less computationally intensive than the FFT method, but requires some degree of preregistration to ensure trustworthy output.D. Centroid RegistrationTaking into account the FFT-derived offsets, a test image is shifted and the areas where the reference stars should lie are input to a center-of-mass (COM) centroiding routine. Each reference star location is expanded by the measurement radius to form a region of interest for each star point. Providing the rotation is within reasonable limits (±5º), the stars will now be within the measurement radius specified during the reference star selection routine (measurement radius < exclusion radius). Each measured test star is compared in X and Y to the parameters found during reference pre-processing, and angular rotation/offset residuals are averaged to generate final rotation and fine offset figures. E. Image TransformationUsing the total linear offset and rotation figures,