physics 114 lecture 20 2d data analysis n.
Skip this Video
Loading SlideShow in 5 Seconds..
Physics 114: Lecture 20 2D Data Analysis PowerPoint Presentation
Download Presentation
Physics 114: Lecture 20 2D Data Analysis

Loading in 2 Seconds...

play fullscreen
1 / 9

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

  • Uploaded on

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.
Download Presentation

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

An Image/Link below is provided (as is) to download presentation

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

Physics 114: Lecture 20 2D Data Analysis

Dale E. Gary

NJIT Physics Department

reminder 1d convolution and smoothing
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
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 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)
point spread function
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
  • 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))
deconvolution problems
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;


    • Zgaus = exp(-(X.^2 + Y.^2)*(5*1.913)^2);
incorrect psf
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
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.