Matlab : Working with Images

1.  Basic Image Display and Processing


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., 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_fig1
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
The 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:
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.

3.1.  Stretching

    The equipment related to the Mars Pathfinder in the lower area of the image is neat, but if we are interested in the Martian surface, this image appears too dark.  Hence, we can try "stretching" the range of data values so that we can see the Martian surface more clearly.  How much should we expand the range?  Let's make a guess and try 0.0  0.15 by typing:
if ishandle(h_fig1)
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

Fig. 2 - An example of a (img) stretched from 0 to 0.15
It is brighter, but if we want to be more precise, we don't have to guess!  Matlab provides several nice ways to interactively print the values of the image.  One is the Data Cursor (), which interactively displays the x-position, y-position, and pixel value at the cursor as you move it around on the image. You can also just type in the line/sample coordinates directly at the command line (e.g., img(262,211) ).    For a more hand-on approach, you can interactively save these values to arrays by using the ginput function and then addressing the img array use the 1D band sequential locations of the chosen line and sample (x and y) axes locations:
[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.

Corresponding rdpix output
x = 211, y = 262, ind = 161542, value = 0.0235857
x = 215, y = 192, ind = 164545, value = 0.0547978

    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)
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
The result is shown in Fig. 3.  Notice how the track marks are accentuated.

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:

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

Once you get profiles working correctly, you get a plot of the original data array (av_img) in its original units, either in the horizontal or vertical direction, along a line you choose with the cursor.  Also, notice that the cursor position is displayed in the plot as a large "+"symbol and its pixel coordinate are shown at the bottom of the profile window.  Please play with it until you feel comfortable.

    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:


Fig. 4 - An example of histogram plots

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:

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

Fig. 8 - Examples of gamma transformation functions created using imadjust .

Below is an example showing the results of a non-linear gamma stretch for three different values of gamma:

img_gamma1 = imadjust(img,[],[],0.4); %gamma transformation with gamma=0.4
img_gamma2 = imadjust(img,[],[],1.0); %gamma transformation with gamma=1
img_gamma3 = imadjust(img,[],[],3.0); %gamma transformation with gamma=3
h_fig6 = figure(); %make a new figure to show the gamma transformations
fig_pos = get(h_fig6,'Position'); %
fig_pos(3) = fig_pos(3).*3; %make the figure 3 times wider so we can put all three images on it
h_ax6a = axes('Position',[0.075+0*(0.05+0.25) 0.05 0.25 0.8]);
h_img6a = imagesc(img_gamma1);
h_title6a = title(h_ax6a,'Gamma = 0.4');
h_ax6b = axes('Position',[0.075+1*(0.05+0.25) 0.05 0.25 0.8]);
h_img6b = imagesc(img_gamma2);
h_title6b = title(h_ax6b,'Gamma = 1.0');
h_ax6c = axes('Position',[0.075+2*(0.05+0.25) 0.05 0.25 0.8]);
h_img6c = imagesc(img_gamma3);
h_title6c = title(h_ax6c,'Gamma = 3.0');
sz = size(img);
colormap gray

Fig. 9 - Examples of gamma stretches created using imadjust .

Gamma transformations can be exspecially useful. To learn more, feel free to read a FAQ put together by Charles Poynton here.

Logarithmic Transformations
    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

Fig. 10 - Example of a log-stretched image expressed in dB .

Contrast-Stretching Transformations
    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.

Fig. 11 - Example of a contrast-stretching transform functions .
Here is an exmaple of applying a contrast-stretching transformation with E=0.4 and m=0.4:
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

Fig. 12 - Example of a contrast-stretched image with E=0.4, m=0.4 .

3.2.  Cropping and Masking

    Suppose we are interested in only a portion of the image - say, the rock named Barnacle Bill (see the picture on the right).  You can do this easily by interactively cropping the image using imcrop:
%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;
set(h_fig7,'Name','Barnacle Bill');
colormap gray;

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

Fig. 13 - Example of a copped image of Barnacle Bill .

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)

Fig. 14 - Example of how to use imfreehand to define an ROI and use it to generaate a mask array .

3.3.  Filtering

    The capabilities described below in this section are just a few of many image processing tools available in Matlab.  Please refer to the various online text books available on the class website as well as the Matlab online help for more information (especially the documentation for the Image Processing Toolbox). 

Smoothing Images
    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));
axis equal
colormap gray
h_ax8 = gca;
h_cbar8 = colorbar('FontSize',14,'FontWeight','Bold');

Fig. 15 - Example of 5x5 Boxcar Average Applied to img. Compare with Fig. 9A. .

Compare the result to Fig. 9A.  Please refer to the Matlab online help for more information on this and also the other capabilities of the fspecial function.

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();
axis equal; colormap gray;
axis equal; colormap gray;
imagesc(edge(img,'sobel')); axis equal; colormap gray;
Frequency Domain Filtering Images
    Filtering in frequency domain is a common image and signal processing technique.  It can be used to smooth, sharpen,de-blur, and restore images.
There are three steps to frequency domain filtering:
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.
    The general form of a frequency filtering command is:
filtered_image = ifft(fft(img)*filter)
where the image is either a vector or a two-dimensional image, the filter is a vector or two-dimensional array designed to filter out certain frequencies in the image, and the image is transformed from the spatial domain to the frequency domain using fft and back using ifft.  Please refer to the Matlab online help for more detailed information.

4.  Making a TIFF File

    You can use the imwrite procedure to write images to TIFF files (as well as other formats). However, first use have to convert the doulbe img array to an 8-bit or 16-bit by array:
minval=min(equalized_img(img ~= 0));
im_out = uint8(round(imrotate( (equalized_img'-minval)./(maxval-minval).*255,90) ));

5.  More, more, more...

    There's lots more in image processing and display.  When you manipulate images mathematically, just as you can any other Matlab variables or array, please remember to work with the original array instead of its byte counterpart from above.  You can play with color tables using colormap. You can rotate images with imrotate or rot, imtranspose, imzoom, draw contours, label your images, overlay coordinates using plot after typing hold on , write text using text, and etc. 

    Please work on your Lab #1 when you are done with this tutorial!

Written by Alexander Hayes and Jim Bell