#!/scisoft/bin/python # #====================================================================================================================== # # 2009-12-03: Mark Neeser # # A demo that shows the following python attributes: # 1. how to create an image built from a (500x500) grid and then assign it to a Gaussian-weighted sinusoid # 2. how to perform a 2-D FFT # 3. how to easily mask the image of part 1. # # #====================================================================================================================== from numpy.fft import fft # importing fast fourier transform modules from scipy import * # importing scipy modules from pylab import * # importing plotting modules #============ # Here we go: #============ lim_radius = 100 # define the radius of the mask x0 = 150 # define the center of the radius y0 = 150 x = y = arange(1,501,1) # create a grid image 500 x 500 pixels in size X, Y = meshgrid(x, y) Z = sin(X/pi)*exp(-((X-250)^2 + (Y-250)^2)/2500.) # the image now = a Gaussian-weighted sinusoid #=================================================== # This is simply for defining plot intensity limits: #=================================================== zmin = Z.min() zmax = Z.max() zmean = Z.mean() zstd = Z.std() z1 = zmean - 3.0*zstd z2 = zmean + 10.0*zstd levels = arange(z1, z2, 0.5) Z_fft = abs(fft2(Z)) # this is the real part of the 2-D FFT #========================================================== # Again, this is simply for defining plot intensity limits: #========================================================== zFmin = Z_fft.min() zFmax = Z_fft.max() zFmean = Z_fft.mean() zFstd = Z_fft.std() zF1 = zFmean - 0.02*zFstd zF2 = zFmean + 0.1*zFstd levels = arange(zF1, zF2, 0.5) figure(1,figsize=(20,8)) # define the figure and its size (in inches) subplot(131) imshow(Z, vmin=z1, vmax=z2, interpolation='nearest', cmap=cm.jet, origin='lower') title(r'$sin(x/\pi ) e^{-((x-250)^2 + (y-250)^2)/2500.}$', fontsize=18) xlabel('X (pixels)') ylabel('Y (pixels)') #contour(Z, levels, hold='on', colors = 'k', origin='lower') subplot(132) imshow(Z_fft, vmin=zF1, vmax=zF2, interpolation='nearest', cmap=cm.jet, origin='lower') #colorbar() title(r'$\rm{Real\ part\ of\ FFT\ [\ } sin(x/\pi ) e^{-((x-250)^2 + (y-250)^2)/2500.}\rm{]}$', fontsize=14) xlabel('X (pixels)') #============================================================ # To show off python's masking abilities, let us create a # mask of the first image with a radius of lim_radius around # a point centered on (x0, y0). # All this will be done WITHOUT do loops! #============================================================ y, x = indices(Z.shape, dtype=float32) # make an image space holder the same size as image 1 (note order of x and y!) x = x - x0 # let's put our center around the point (x0, y0) y = y - y0 radius = sqrt(x**2 + y**2) mask = radius < lim_radius # this "mask" is an image with values="False" outside of lim_radius, and "True inside subplot(133) imshow(mask*Z, vmin=z1, vmax=z2, interpolation='nearest', cmap=cm.jet, origin='lower') title1 = r'$\rm{Mask\ where\ [\ } sin(x/\pi ) e^{-((x-250)^2 + (y-250)^2)/2500.}\rm{]}\ <\$' + str(lim_radius) title2 = 'radius centered on the point (' + str(x0) + ',' + str(y0) + ')' title(title1, fontsize=14) xx1, xx2 = xlim() # extract the limits of the plot yy1, yy2 = ylim() ylim() text(xx2+10,yy1+5, title2, color='r', size=12, rotation='vertical') xlabel('X (pixels)') xlim(1,500) ylim(1,500) #============================================================================= # Let's now sum up some stuff within the masked area (look Mom, no do loops!): #============================================================================= sum_pixels = mask.sum() # this gives you the number of pixels contained within the radius < lim_radius flux = (mask*Z).sum() # Sum up all the pixel values in the image that are within the radius < lim_radius # sum_mask = Z[where(mask)].sum() would give you the SAME answer netflux = flux/sum_pixels # this gives you the number of pixels contained within the radius < lim_radius title3 = 'The total number of pixels within the mask = %5i'%(sum_pixels) title4 = 'The summed up flux contained within the mask = %6.3f counts'%(flux) title5 = 'The average flux contained within the mask = %8.5f counts'%(netflux) text(xx1, yy1-70, title3, color='b', size=12) text(xx1, yy1-90, title4, color='b', size=12) text(xx1, yy1-110, title5, color='b', size=12) show()