Physics 114: Lecture 20 2D Data Analysis

1. Physics 114: Lecture 20 2D Data Analysis Dale E. Gary NJIT Physics Department

2. Reminder1D Convolution and Smoothing • Let’s create a noisy sine wave: • u = -5:.1:5; • w = sin(u*pi)+0.5*randn(size(u)); • plot(u,w) • We can now smooth the data by convolving it with a vector [1,1,1], which does a 3-point running sum. • wsm = conv(w,[1,1,1]); • whoswsm Name Size Bytes Class Attributes wsm 1x103 824 double • Notice wsm is now of length 103. That means we cannot plot(u,wsm), but we can plot(u,wsm(2:102)). Now we see another problem. • Try this: • plot(u,w) • hold on; plot(u,wsm(2:102)/3,’r’,’linewidth’,2)

3. 2D Convolution • To do 2D smoothing, we will use a 2D kernel k = ones(3,3), and use the conv2() function. So to smooth the residuals of our fit, we can use • zsm = conv2(ztest-z,k)/9.; • imagesc(x,y,zsm(2:102,2:102)) • Now we can see the effect of missing the value of cx by 0.05 due to our limited search range. • There are other uses for convolution, such as edge detection. For example, we can convolve with a kernel k = [1,1,1,1,1,1]. • Or a kernel k = [1,1,1,1,1,1]’. Or even a kernel k = eye(6). Or k = rot90(eye(6)).

4. Convolution and Resolution • Convolution can be used for smoothing data, but it is also something that happens naturally whenever measurements are made, due to instrumental resolution limitations. • For example, an optical system (telescope or microscope, etc.) has an aperture size that limits the resolution due to diffraction (called the diffraction limit). Looking at a star with a telescope, assuming no other effects like atmospheric turbulence, results in a star image of a certain size, surrounded by an “airy disk” with diffraction rings. • This shape is mathematically just the sinc() function we introduced last time: • x = -5:0.1:5; y = -5:0.1:5; • [X, Y] = meshgrid(x,y); • Z = sinc(sqrt(X.^2 + Y.^2)); • imagesc(x,y,Z); • In fact, this is the electric field pattern, and to get the intensity we need to square the electric field: imagesc(x,y,Z.^2)

5. Point Spread Function • To show this point better, consider a “perfect” instrument that perhaps has noise, but shows stars as perfect point sources. Let’s generate an image of some stars: • stars = randn(size(X))*0.1; • stars(50,50) = 1; stars(20,37) = 4; stars(87,74) = 2; stars(45,24) = 0.5; • imagesc(stars) • To see the effect of observing such a star pattern with an instrument, convolve the star image with the sinc function representing the diffraction pattern of the instrument (the point spread function or PSF): • Z = sinc(sqrt(X.^2 + Y.^2)*5).^2; % the *5 makes it smaller/sharper • imagesc(conv2(stars,Z)) • You see that the result is twice as large due to the way convolution works. Try • fuzzy = conv2(stars,Z); colormap(gray(256)); • imagesc(stars); axis square • imagesc(fuzzy(51:150,51:150)); axis square

6. Deconvolution • It is actually possible to do the inverse of convolution, called deconvolution. Let’s read in an image and fuzz it up (download fruit.gif from course web pg) • [img map] = imread(‘fruit.gif’); • fuzzy = conv2(single(img),Z)/sum(sum(Z); • image(img) % Original image—observe the sharpness • image(fuzzy(51:515,51:750)) % fuzzy image • Now let’s sharpen it again. MatLAB has a family of deconvolution routines. The standard one is deconvreg(): • image(deconvreg(fuzzy,Z)) • The image is dark, because we have to normalize the convolving function: • image(deconvreg(fuzzy,Z)*sum(sum(Z))) • This looks pretty good, but note the edge effects. Try another routine • image(deconvlucy(fuzzy,Z)*sum(sum(Z))) • This one looks almost perfect. However, if you compare images you do see differences • sharp = deconvlucy(fuzzy,Z)*sum(sum(Z)) • imagesc(sharp(51:515,51:750) – single(img))

7. Deconvolution Problems • Any time you do an inversion of data, the result can be unstable. Success depends critically on having the correct point spread function. • The deconvolution we just did was after convolving the image with a “perfect” instrument and neglecting atmospheric turbulence. Further blurring by the atmosphere acts to increase the size of the “airy disk” and smear out the diffraction rings. • With some time average, the above pattern smears out into an equivalent gaussian. The equivalent gaussian to • Zsinc = sinc(sqrt(X.^2 + Y.^2)*5).^2; is • Zgaus = exp(-(X.^2 + Y.^2)*(5*1.913)^2);

8. Incorrect PSF • Let’s convolve the image with the Gaussian (i.e. instrument plus atmospheric turbulence), creating a larger PSF • Zgaus = exp(-(X.^2 + Y.^2)*(3*1.913)^2); % Note use of 3 to enlarge Gaussian • Convolve with this blurred PSF • fuzzy = conv2(double(img),Zgaus)/sum(sum(Zgaus)); • image(fuzzy) • Now deconvolve with the instrumental PSF • dconl = deconvlucy(fuzzy,Zsinc)*sum(sum(Zsinc)); • image(dconl) • We see that we cannot recover the original instrumental resolution. The clarity is lost due to atmospheric turbulence. • However, if we measure the PSF of instrument plus atmosphere, we CAN recover the blurring due to the atmosphere.

9. Laser Guide Stars • Astronomers now use a laser to create a bright, nearby “guide star” near the region of the sky of interest. • By imaging the laser scintillation pattern instantaneously, they can freeze the atmosphere and correct the images in real time.