Physics 114: Lecture 20 2D Data Analysis

1 / 9

Physics 114: Lecture 20 2D Data Analysis - PowerPoint PPT Presentation

Physics 114: Lecture 20 2D Data Analysis. Dale E. Gary NJIT Physics Department. Reminder 1D 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

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

PowerPoint Slideshow about 'Physics 114: Lecture 20 2D Data Analysis' - dolores

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript

Physics 114: Lecture 20 2D Data Analysis

Dale E. Gary

NJIT Physics Department

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)
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)).

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)
• 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
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)
• 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))
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);
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.
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.