Before we get started, download the LAB 1 data package from the class website here and unpack the archive. If you are a Linux or Mac user, you can do this using the commands "gunzip LAB1_DATA_Package.tar.gz" and "tar -xvf LAB1_DATA.Package.tar" from the command line. If you are a windows user, you can use your own command line environment (e.g., http://www.cygwin.com) or download and install the freeware program 7-ZIP . Once you have downloaded and extracted the data you will end up with a directory named LAB1 with sub-directories DATA, SUBROUTINES, BACKGROUND, and REPORT. Start Matlab and move to the LAB1/REPORT directory. Once there, you can add the subroutines stored in the SUBROUTINES direcory and data sored in the DATA directory to you Matlab path by typing:
Now, Matlab will know how to find the sub-routines and data required for this tutorial and LAB 1. In the LAB1/DATA directory, there is a nice image of the Martian surface in an old flood channel, Ares Vallis, obtained using the IMP (Imager for Mars Pathfinder), LAB1/DATA/imp.fits. This file is in a format called "FITS" format, which is the most popular format for astronomical images. To read FITS data files into Matlab arrays, we will use the procedure "fitsread.m", which is built into matlab. However, we have included an updated "fitsread.m" in the SUBROUTINES directory of LAB1. This is because the Matlab version of fitsread will not work for Linux or Mac environments without a small modification. We will go over this modification, which defaults all environemnts to big-endian format, during class. All you have to do is type:
img = fitsread('imp.fits');
headerinfo = fitsheader('imp.fits');%displays the header information to the screen
These commands will return the image array (img) and diplay information about the image (header). Note that the semi-colons (";") at the end of each command are important. If you don't include them, all of the values in the image and header will display to the screen.
head_date = fitsheader('imp.fits','DATE'); % returns the header keyword value "DATE"
whos img % shows basic information regarding the image array
The output from the last two commands returns a string array with the contents of the header keyword DATE and tells you that the img array is a 768 x 632 DOUBLE array. You can use the max and min functions to determine the range of the data values as shown.
minval = min(img(:)); maxval = max(img(:)); disp([minval,maxval]);
The first two commands store the minimum and maximum values while the third displays it to screen. You could have also simply left the semicolons off the end of the first two commands. The colon (":") is used to convert the image into a 470016 X 1 DOUBLE array before calculating the minimum and maximum. If we had used min(img) and max(img) , Matlab would have returned a 1 X 612 DOUBLE arrays populated by the minimum and maximum of the image on a row-by-row basis. The range of the image is from 0.0 to 0.274653. We should note that, if the img array was directly converted to a byte array, all the data values would become byte 0. This will become relevant when we attempt to export the data into a standard image format, such at JPEG or TIFF. But first, lets create an appropriately sized graphics window and try displaying the image:
h_fig1 = figure();%make the figure and place its handle in the variable h_fig1The first thing you will notice is that the image is upside down. This is because of the way that matlab reads arrays. In order to display it upright, we can either flip the display or flip the image. To flip the image we could type img = img(end:-1:1,:); . However, lets flip the display axis instead. To do this, type:
fig_pos = get(h_fig1,'Position');%return the figures position on the screen
fig_pos(3:4) = [612,768];%set the figure size in display pixels to be the size (or at lest aspect ratio) of our image
set(h_fig1,'Position',fig_pos,'Name','Martian Surface');%sets the figure size/position and puts "Martian Surface" in the title bar
h_img1 = imagesc(img);%displays the image and returns the handle for the image object
h_ax1 = gca;%gets the handle for the current axis
colormap gray;%set the colormap to grayscale (default is a blue to red jet colormap)
h_cbar1 = colorbar;%put up a colorbar to show the mapping between grayscale and pixel value
axis equal;%automatically scales the axis to be equant
set([h_ax1,h_cbar1],'FontSize',14,'FontWeight','Bold');%makes the tick marks and tick labels easier to read
h_title1 = title('Martian Surface (MPF)','FontSize',18,'FontWeight','Bold');%gives the axes a Title
set(h_ax1,'Position',[0.05,0.05,0.9,0.9]);%expand the size of the axis in the figure to remove blank space
set(h_ax,'YDir','Normal'); %change Y-Direction from "Reverse" to "Normal".
Matlab automatically inverts the y-direction of the axis when displaying images using imagesc .
Fig. 1 - An example of a (img) stretched from its minimum to maximum values
There are many ways to display images in Matlab. I prefer imagesc because it will automatically scale the colormap to the minimum and maximum pixel values in the array. You can stretch linearly stretch the displayed image using imagesc by simply passing a 1x2 vector with the desired minimum and maximum values into the command. In fact, the image looks a bit dark so lets do that. But first, we'll need to decide what the minimum and maximum values of the stretch should be.
h_fig2 = figure('Position',get(h_fig1,'Position')); %make a new figure the same size and location as the first
h_fig2 = figure(); % make a new figure
h_img2 = imagesc(img,[0,0.15]); %display image stretched from 0 to 0.15
colormap gray; %set the colormap
h_cbar2 = colorbar; %open the colorbar
h_ax2 = gca; %save the handle for the new axis
set(h_ax2,'YDir','Normal'); %flip the image as before
axis equal; %automatically scales the axis to be equant
set([h_ax2,h_cbar2],'FontSize',14,'FontWeight','Bold'); %makes the tick marks and tick labels easier to read
h_title2 = title('Martian Surface (MPF Stretched Image)','FontSize',18,'FontWeight','Bold'); %add a title
if ishandle(h_ax1) ;%
set(h_ax2,'Position',get(h_ax1,'Position'));%if you still have the other figure open, set the axis to the same size
[x,y] = ginput(4);%grabs four positions on the image by left-clicking on them (if you did not include a #, the function would continue grabbing values until you right-clicked your mouse).
ind = sub2ind(size(img),round(y),round(x));%turns the line/sample coordinates into 1D band sequential indices for addressing.
z = img(ind);%address the image using 1D band sequential indices
Alternatively, you could use the built in impixel function to do much the same.
Suppose we want to focus on the rover track marks in the image, then use the Data Cursor () or ginput to find the data values at the spots representative of the lowest and highest values associated with this region - for example, the locations marked respectively 1 and 2 as shown on the right. These selected locations produced the following output.
Now, use the imagesc command to display the img array again as above, but limit the range using the minimum and maximum data values in the vicinity of the rover tracks. This will scale the colormap data range between 0.0236 and 0.0548. This will display any original data values that exceed 0.0548 as white (or whatever the largest value in the colormap is set to), and any that are smaller than 0.0236 equal to black (or the lowest value in the colormap), obliterating information on the stronger and weaker features in the displayed image:
if ishandle(h_fig2)The result is shown in Fig. 3. Notice how the track marks are accentuated.
h_fig3 = figure('Position',get(h_fig2,'Position')); %make a new figure the same size and location as the first
h_fig3 = figure(); % make a new figure
h_img3 = imagesc(img,[0.0236,0.0548]); %display image stretched from 0.0236 to 0.548
colormap gray; %set the colormap
h_cbar3 = colorbar; %open the colorbar
h_ax3 = gca; %save the handle for the new axis
set(h_ax3,'YDir','Normal'); %flip the image as before
axis equal; %automatically scales the axis to be equant
set([h_ax3,h_cbar3],'FontSize',14,'FontWeight','Bold');%makes the tick marks and tick labels easier to read
h_title3 = title('Martian Surface (MPF Stretched Image)','FontSize',18,'FontWeight','Bold'); %add a title
if ishandle(h_ax2) ;%
set(h_ax3,'Position',get(h_ax2,'Position'));%if you still have the other figure open, set the axis to the same size
Fig. 3 - An example of stretched image set to maximum contrast in the rover tracks
Nicer than the brute force ginput algorithm or even the more elegant impixel function is improfile, which interactively draws profiles of an image in a separate newly-created window. With some inginuitiy, it can be used to show the values of either the data array you are actually seeing (img) or any other data array of the same dimensions (img). You will be using th is capability in a later lab to extract elevation values for bedding planes that you'll identify from an image. To see how it works, type:
improfileafter setting the current figure to the one you want to use to identify a profile. Then, select two or more points to define the profile on the image, right-clicking when you are done. You can also define the number of click points when you call the function (see the matlab help). Note that the new window will show a 2D plot if you select a straight line and a 3D plot if your profile defines a more complex path.
Often even better is to make a histogram of the image pixel values to see their distribution. To make and display a hitogram of the pixel values in img with 500 bins use the following commands:
h_fig4 = figure();%open a new figure for the histogram
h_hist4a = histogram(img,500);%creates a histogram of the image
h_ax4a = gca;%saves the axis handle
set(h_ax4a,'FontSize',14,'FontWeight','Bold');%cleans up the fonts for better viewing
h_ax4a_xlab = xlabel(h_ax4a,'Pixel Value','FontSize',16,'FontWeight','Bold');%inserts an x-label
h_ax4a_ylab = ylabel(h_ax4a,'# of Pixels','FontSize',16,'FontWeight','Bold');%inserts a y-label
You'll notice the the resulting histogram shows a significant number of pixels with a value of exactly 0. These represent regions of the array that do not have valid data. To better see the pixel distribution, you can generate another histogram that ignores these values:
h_ax4b = axes();%generates a new axes for the next histogram
set(h_ax4a,'Position',[0.14,0.6,0.7750,0.35]); %moves the h_ax4a axes out of the way
set(h_ax4b,'Position',[0.14,0.1,0.7750,0.35]); %moves the h_ax4b axes into the now empty space
h_hist4b = histogram(img(img ~= 0),500);%creates a histogram of the image using pixel values not equal to zero
set(h_ax4b,'FontSize',14,'FontWeight','Bold');%cleans up the fonts for better viewing
h_ax4b_xlab = xlabel(h_ax4b,'Pixel Value','FontSize',16,'FontWeight','Bold');%inserts an x-label
h_ax4b_ylab = ylabel(h_ax4b,'# of Pixels','FontSize',16,'FontWeight','Bold');%inserts a y-label
It appears that most of the information is contained within pixel values between 0 and 0.1. We can use the 'XLim' parameter of the axes handles to zoom our plots into this region:
Looking at the histogram, it looks like a decent stretch would be from 0 to 0.08. Using the prctile or quantile functions in matlab, which tell you percentile (0-100) or quantile (0-1) of particular values in a given dataset, we can find that 85% of the images values are contained betwen 0 and 0.08. Another way to determine an optimum linear stretch would be to calculate the mean (or median) and standard deviation of valid pixels in img and then set the stretch to be between the mean (or median) plus/minus n standard deviations . In our case, a 1 std contrast stretch would be from 0.01 to 0.1. As long as you still have h_fig3 open, we can change the range to stretch the displayed image by changing the "CLim" paramters in the axes handle:
set(h_ax3,'CLim',[0,0.08])The resulting display is shown in Fig.6 below.
Fig. 5 - An example of an image scaled from 0 to 0.08 using histogram
So far, we have only concerned ourselves with linear stretching. Thre are also a host of non-linear stretching algorithms. Non-linear stretching is useful to extract detail from one end of the pixel distribution (either bright or dark values) without eliminating all contrast amongst pixel values on the opposite end of the distribution. Examples of non-linear stretching include Histogram Equilization , Logarithimic Stretching , Gamma Transformations , and Contrast-Stretching Transformations.
Histogram Equalization is commonly used nonlinear mapping technique, in which the pixel value color mapping is modified so that all of the colormap values are used in an equal number of pixels. As in our example above (Fig. 4), the distribution of pixel values in an image are often found to be clustered, giving the image has a narrow dynamic color/intensity range. If you spread the pixel distribution out, so that each range of values has approximately the same number of pixels, the information content of the image can sometimes be enhanced. This process of spreading the pixel distribution over the entire dynamic color/intensity range is known as histogram equalization, and the histeq function can be used to do this easily. Try the following example.
equalized_img = histeq(img./max(img(:)),1000); %created a histogram equalized image by spreading the pixel values out into 100 evenly spaced bins (note histeq expects an image with values between 0 and 1 for a double array
set(h_img3,'CData',equalized_img); %replace the data array used to display h_img3 in h_fig3 with the equlized image
set(h_ax3,'CLim',[min(equalized_img(:)),max(equalized_img(:))]);%reset the stretch in h_ax3 to match the minimum and maximum of the equalized image
set(h_title3,'String','Martian Surface (MPF Histogram Equalized Image)');%change the image title
Fig. 5 - An example of a hitogram equalized image
In the above example, instead of creating a new figure for the histogram equialized image, we have simply replaced the data array (parameter CData in the axes object h_ax3 ) used to display the image in h_fig3. This is one of the benefits of Matlab's object based graphics environment.
We can also compare the histograms between the original and histogram equalized images:
h_fig5 = figure();%open a new figure for the histogram
h_hist5a = histogram(img(img ~= 0)./max(img(:)),100);%creates a histogram of the image scaled from 0-1
h_ax5a = gca;%saves the axis handle
set(h_ax5a,'FontSize',14,'FontWeight','Bold');%cleans up the fonts for better viewing
h_ax5a_xlab = xlabel(h_ax4,'Pixel Value','FontSize',16,'FontWeight','Bold');%inserts an x-label
h_ax5a_ylab = ylabel(h_ax4,'# of Pixels','FontSize',16,'FontWeight','Bold');%inserts a y-label
h_ax5b = axes();%generates a new axes for the next histogram
set(h_ax5a,'Position',[0.14,0.6,0.7750,0.35]); %moves the h_ax5a axes out of the way
set(h_ax5b,'Position',[0.14,0.1,0.7750,0.35]); %moves the h_ax5b axes into the now empty space
h_hist5b = histogram(equalized_img(img ~= 0),100);%creates a histogram of the equalized image
set(h_ax5b,'FontSize',14,'FontWeight','Bold');%cleans up the fonts for better viewing
h_ax5b_xlab = xlabel(h_ax5b,'Pixel Value','FontSize',16,'FontWeight','Bold');%inserts an x-label
h_ax5b_ylab = ylabel(h_ax5b,'# of Pixels','FontSize',16,'FontWeight','Bold');%inserts a y-label
h_title5a = title(h_ax5a,'Original Image','FontSize',18,'FontWeight','Bold');
h_title5b = title(h_ax5b,'Histogram Equalized Image','Fontsize',18,'FontWeight','Bold'); set([h_ax5a],'XLim',[0,1]);
Fig. 7 - Histogram of a histogram equalized image.
Gamma Transformations can curve the grayscale components to either brighten the intensity (gamma parameter < 1) or darkent the intensity (gamma parameter > 1) using a non-linear power law relationship between input and output pixel values. In this way, the gamma transformation can enhanced contrast for either the darker or brighter parts of the intensity distribution. The Matlab funciton that creates these non-linear gamma transformations is imadjust . The imadjust function can also linearly scale or clip the image prior to performing the gamma transformation by providing pixel value ranges for the input and output images (see online Matlab help for imadjust). For the example below, we will ignore this capability and focus on the gamma transformation itself, mapping the full range of the input into the full range of the output. Fig. 8 below shows the shape of the gamma transformation for varying values of gamma. Notice tha the red line has gamma=0.4, which creates an upward curve and will brighten the image by expanding contrast at the lower pixel instensities.
Below is an example showing the results of a non-linear gamma stretch for three different values of gamma:
Gamma transformations can be exspecially useful. To learn more, feel free to read a FAQ put together by Charles Poynton here.
Logarithmic Transformations can be used to brighten the intensities of an image (like the Gamma Transformation for gamma less than 1) to increase the contrast of lower pixel values at the expense of the high value pixels. The general equation of a logarithmic transform is img_log = c*log(img)/log(d), where the paramter c sets the magnitude of image brightness and d determines base of the logarithim. A standard logarithmic transform is img_log = 10*log10(img), which is expressed in units of decibels. Note that image values below 1 will become negative and image values less than or equal to zero are non-physical in the logarithmic transform. As a result, it is important to either remove these values from your image ahead of time or use the find function to only perform the logarithmic transform on the valid portions of the image. Below is an example that converts our MPF image into decibels of I/F and displays the result in h_fig3.
bad_ind = find( img == 0);%find the band sequential indices of the invalid pixels
good_ind = find( img > 0);%find the band sequential indices of the valid pixels above 0
if numel(bad_ind)+numel(good_ind) == numel(img); display('No Pixels < 0'); end;
img_log = img;%define an array for the logarithmic image
img_log(good_ind) = 10.*log10(img(good_ind));%take the logarithm
img_log(bad_ind)=nan;%set the invalid pixels to nan or "Not A Number"
set(h_img3,'CData',img_log)%replace the data array in h_img3 with the logarithmic image
set(h_ax3,'CLim',[-25,-5])%reset the stretch range of the imagesc command
set(h_title3,'String','Martian Surface (Logarithimic Transform) [dB]');%change the image title
Contrast-stretching transformations increase the contrast between the dark and light pixel values. The process is similar to the histogram equialization algorithm in that you can stretch the histogram to fill up the image's intensity domain, although contrast-stretching allows you to stretch the intensity values around a certain level. You end up with darker darks and lighter lights, with only a few gray levels around the central region of the stretch. The general equation for the contrast-stretching transformation is img_constretch = 1./(1 + (m./(img + eps)).^E). E controls the slope of the function and m is the mid-line where you want to switch the dark values to light values. eps is a small constant (~1e-10) that is used to avoid dividing by zero. Fig. 11 shows the transform function for several values of E and m.
img_constretch = 1./(1+(0.4./(img+1e-10)).^(0.4));%applies the constrast-stretch transformation
set(h_img3,'CData',img_constretch)%replace the data array in h_img3 with the contrast-stretched image
set(h_ax3,'CLim',[0.11,0.45])%reset the stretch range of the imagesc command
set(h_title3,'String','Martian Surface (Contrast-Stretch Transform) [E=0.4 m=0.4]');%change the image title
%img2 = imcrop(h_fig3);%will allow you to interactively define a region of interest (ROI) on h_fig3
img2 = imcrop(img,[65 327 186 125]);%crops img at specific locations around Barnicle Bill
h_fig7 = figure();
h_img7 = imagesc(img2);
h_ax7 = gca;
Note that the imcrop function will also work by passing it an array. In that case, the function with either create a figure or simply overplot the array into the currently selected figure. The downside of this option, however, is that no stretching will be performed so you may not be able to easily see you region of interest (ROI).
Matlab also has a seriews of powerful tools for defining regions of interest that you can use to define areas in which you want to apply analysis algorithms. The functions imfreehand, imellipse, and imrectangle allow you to interactively define regions of interest on a figure. Once the ROI is defined, you can use it to perform a host of useful tasks. One of the most useful is generating a mask array populated by ones in the ROI and zeros outside of it. The following example shows you how to use imfreehand to define an ROI around Barnacle Bill and then generate an image mask. When image masks are coupled to Matlab's find command, it is a powerful combination that will be used repeatedly over the course of the semester (including in LAB 1).
roi = imfreehand;%interactively draw the roi on h_ax3
mask=createMask(roi);%creates the image mask, NOTE THAT THIS FUNCTION IS CASE SENSITIVE
roi_ind = find(mask==1);%find the pixels in the mask that are in the ROI (i.e., have pixel value 1)
disp(mean(img(roi_ind)));%displays the mean value of the img in the vicinity of Barnacle Bill (0.0328 for my ROI)
Images can be smoothed by averaging each pixel value with the values of its surrounding neighbors. This is known as mean or boxcar smoothing, and it is accomplished in Matlab using the the fspecial and imfilter functions. fspecial defines a 2-D filter and imfilter applies it via convolution. fspecial can generate any number of filters, but the simplist is an equally weighted window using a square neighborhood of a given odd width. For example, if the neighborhood is 3 by 3, thent the filter is generated by the command filt = fspecial('average',3); and the subsequent command img_smooth = imfilter(img,filt); would replace each pixel by the mean value of it and its eight surrounding nearest neighbors (except at the edges). Let's try the fspecial and imfilter functions on img using a 5 by 5 boxcar average.
filt = fspecial('average',5);%defines a 5x5 uniform weight boxcar filter
smoothed_img = imfilter(img,filt);%apply the filter to img
h_fig8 = figure('Position',get(h_fig3,'Position'));
h_img8 = imagesc(imadjust(smoothed_img,,,0.4));
h_ax8 = gca;
h_cbar8 = colorbar('FontSize',14,'FontWeight','Bold');
Removing Noise From Images
A common form of noise is salt and pepper noise, in which random pixels in the image have extreme values. The medfilt2 function is an excellent choice for removing this kind of noise from images. It is similar to the fspecial and imfilter functions discussed above, except that it calculates the median value of the pixel neighborhood instead of the weighted mean value. This has two important effects. First, it can eliminate extreme values in images. Second, it does not blur the edges or features of images whose sizes are larger than the neighborhood. Please refer to the Matlab online help for more details.
Enhancing the Edges of Images
An image may be sharpened or have its edges enhanced by differentiation. Matlab provides a built-in edge enhancement function edge that can use several different edge detecting algorithms including Canny, Log (Laplacian of Gaussian), Prewitt, Roberts, Sobel, and Zero Cross. Try the following commands, and see the Matlab online help for more information.
h_fig9 = figure();Frequency Domain Filtering Images
axis equal; colormap gray;
axis equal; colormap gray;
imagesc(edge(img,'sobel')); axis equal; colormap gray;
There are three steps to frequency domain filtering:The general form of a frequency filtering command is:
1. The image is transformed from the spatial domain into the frequency domain using the Fast Fourier Transform (FFT).
2. The transformed image is multiplied by a frequency filter.
3. The filtered image is transformed back to the spatial domain.
minval=min(equalized_img(img ~= 0));
im_out = uint8(round(imrotate( (equalized_img'-minval)./(maxval-minval).*255,90) ));
Please work on your Lab #1 when you are done with this tutorial!