Image processing tips for geoscientists – 1

Today I would like to show a way to quickly create a pseudo-3D display from this map:

Original image

The map is a screen capture of a meandering river near Galena, Alaska, taken in Google Earth. I love this image; it is one of my favorite maps for several reasons. First of all it is just plainly and simply a stunningly beautiful image. Secondly, and more practically, the meanders look not too dissimilar to what they would appear on a 3D seismic time slice displayed in grayscale density which is great because it is difficult to get good 3D seismic examples to work with. Finally, this is a good test image from the filtering standpoint as it has a number of linear and curved features of different sizes, scales, and orientation.The method I will use to enhance the display is the shift and subtract operation illustrated in The Scientist and Engineer’s Guide to Digital Signal Processing along with other 3×3 edge modification methods. The idea is quite simple, and yet extremely effective – we convolve the input image with a filter like this one:

0  0  0
0  1  0
0  0 -1

which will create a shifted version of the image (at the -1 location) and subtract it from the original (at the 1). Quoting the author:” this processing produces the optical illusion that some objects are closer or farther away than the background, making a 3D or embossed effect. The brain interprets images as if the lighting is from above, the normal way the world presents itself. If the edges of an object are bright on the top and dark on the bottom, the object is perceived to be poking out from the background”.

This is the same result you would achieve using Emboss filter in Photoshop and other similar programs. The only time I saw embossing used on seismic data was a few years ago when I took a  Seismic Geomorphology course taught by Henry Posamentier for the Canadian Society of Exploration Geophysicists (CSEG). So thank you for the idea, Henry.

Like with derivative maps and with illumination, the outcome of this operation is not unique, in the sense that with changing the direction (or azimuth) of the process, i.e. the position of the -1, different features will be enhanced, and some in fact will be obscured as a result of the subtraction. Let’s look at a few examples below.

Here is the shift and subtract to the South, using this kernel:
 0  0  0
 0  1  0
 0 -1  0

It looks like a good result, with most linear and curved features are enhanced. It has however suppressed features oriented N-S. For example many of the channels in the NE corner, and the “l” in Google, which is all but gone. Notice also that the pseudo-3D effect is noticeable but not as pronounced as in the shift and subtract to the Southwest below. Look at the main river on the left for instance. The reason is that when applied along the diagonals, the copied image is shifted both horizontally and vertically, as opposed to just vertically or just horizontally:
0  0  0      0  0  0
0  1  0      0  1  0
0 -1  0     -1  0  0

The shift and subtract to the West delivers the least pseudo-3D effect. This is because not only are we applying the shift in one direction, but also that happens to be the direction a lot of the channels in the center are oriented. So directionality is a drawback of the method, but on the other hand I find very often that going through the exercise of comparing the alternatives I will learn a lot more about my data and the geological features in it than if I hadn’t, so I do it if I have the opportunity.

Try it yourself, and let me know if you agree with my conclusions. Did you find any of the remaining 5 direction better overall? The original image is here, and the code to create all filters and output displays is below:

%% Shift and subtract script
% A simple script to apply shift and subtract edge enhancement
% Additional functionality to display and QC results
%
% Reference: The Scientist and Engineer's Guide to Digital Signal Processing
% Chapter 24 - Linear Image Processing
% 3x3 Edge Modification
% http://www.dspguide.com/ch24/2.htm
%
% Acknowledgements: linkzoom function from the Matlab File Exchange
% www.mathworks.com/matlabcentral/fileexchange/21414-linkzoom-m-v1-3-aug-2009
%
%
%% 1 - dealing with input
% once imported jpeg file as galenaAlaska
%
in=zeros(1094,1664);
n=1:1:1094; k=1:1:1664;
in(n,k)=0.2989 * galenaAlaska(n,k,1) + 0.5870* galenaAlaska(n,k,2) ...
 + 0.1140 * galenaAlaska(n,k,3);
%
%
%% 2 - initializing
% this section not only create and applies the filters through convolution
% but also stores input, filters, and output in portable multi-dimentional
% arrays for ease of access
%
% input image array
IN=repmat(in,[1 1 8]);
% filtered output array
[mi,ni]=size(in);
[mf,nf]=size(flt);
OUT=zeros(mi+mf-1,ni+nf-1,8);
% filters array
flt=zeros(3,3); flt(2,2)=1;
FLT=zeros(3,3,8);
% this is the main part, populates filters, performs the filtering,
% and outputs filtered images
idx=[1,2,3,4,6,7,8,9];
for i=1:8
 s=idx(i);
 flt1=flt;flt1(s)=-1;
 FLT(:,:,i)=flt1;
 out=conv2(flt1,in);
 OUT(:,:,i)=out;
end
%
% to give you an idea of how it is done, just type at the prompt:
% flt
% which will print:
% flt =
%
%  0  0  0
%  0  1  0
%  0  0  0
% which is the initial kernel without shifted copies
%
% and then type:
% FLT
% which will print:
% FLT(:,:,1) =
%
% -1  0  0
%  0  1  0
%  0  0  0
%
%
%FLT(:,:,2) =
%
%  0  0  0
% -1  1  0
%  0  0  0
%
%
%FLT(:,:,3) =
%
%  0  0  0
%  0  1  0
% -1  0  0
%
% ...
% ...
% FLT(:,:,8) =
%
%  0  0  0
%  0  1  0
%  0  0 -1
% which are all the filters to apply shift and subtract
% in the 8 neighbors of the central pixel

For plotting you try use something like these few lines, which is what I used for the three examples above:

imagesc(OUT(:,:,1),[-25 25]);colormap(gray);
axis tight
axis equal
axis off

and for the original:

imagesc(in,[0 255]);colormap(bone);
axis tight
axis equal
axis off

However, it is far better for comparison and QC to display all the directions at once together with the original, like this:

I did that with this additional code below. Notice that to run it you will need to download from the Matlab File Exchange the function linkzoom by Carlos Vargas Aguilera. Linkzoom allows synchronizing of all subplots so that when zooming or panning one of the images, all other images will zoom or pan accordingly.

%% 3 - displaying
% displays the input image and all results as a matrix
% uses function linkzoom.m from Matlab File Exchange to allow
% simultaneous zooming and panning of all subplots at once
%
figure;
hold;
ax= zeros(3,3);
ax(1)=subplot(3,3,1);imagesc(OUT(:,:,1),[-25 25]);colormap(gray);
set(ax(1),'Position',[0.13,0.7093,0.2774,0.2298]);
axis tight
axis equal
axis off
title('NW','fontsize',14);
ax(2)=subplot(3,3,2);imagesc(OUT(:,:,4),[-25 25]);colormap(gray);
set(ax(2),'Position',[0.4108,0.7093,0.2774,0.2298]);
axis tight
axis equal
axis off
title('N','fontsize',14);
ax(3)=subplot(3,3,3);imagesc(OUT(:,:,6),[-25 25]);colormap(gray);
set(ax(3),'Position',[0.6916,0.7093,0.2774,0.2298]);
axis tight
axis equal
axis off
title('NE','fontsize',14);
ax(4)=subplot(3,3,4);imagesc(OUT(:,:,2),[-25 25]);colormap(gray);
set(ax(4),'Position',[0.13,0.4096,0.2774,0.2298]);
axis tight
axis equal
axis off
title('W','fontsize',14);
ax(5)=subplot(3,3,5);imagesc((in),[0 255]);colormap(bone);
set(ax(5),'Position',[0.4108,0.4096,0.2774,0.2298]);
axis tight
axis equal
axis off
title('ORIGINAL','fontsize',14);
ax(6)=subplot(3,3,6);imagesc(OUT(:,:,7),[-25 25]);colormap(gray);
set(ax(6),'Position',[0.6916,0.4096,0.2774,0.2298]);
axis tight
axis equal
axis off
title('E','fontsize',14);
ax(7)=subplot(3,3,7);imagesc(OUT(:,:,3),[-25 25]);colormap(gray);
set(ax(7),'Position',[0.13,0.11,0.2774,0.2298]);
axis tight
axis equal
axis off
title('SW','fontsize',14);
ax(8)=subplot(3,3,8);imagesc(OUT(:,:,5),[-25 25]);colormap(gray);
set(ax(8),'Position',[0.4108,0.11,0.2774,0.2298]);
axis tight
axis equal
axis off
title('S','fontsize',14);
ax(9)=subplot(3,3,9);imagesc(OUT(:,:,8),[-25 25]);colormap(gray);
set(ax(9),'Position',[0.6916,0.11,0.2774,0.2298]);
axis tight
axis equal
axis off
title('SE','fontsize',14);
linkzoom('xy');
zoom on;

In future posts I will show how to get non-directional dependent results with several methods:

– using gradient to render shading (similar to bump mapping) – UPDATE: see my post visualization tips for geoscientists – Matlab

local binary patterns

– using high pass filtered version of the data to render shading

10 responses to “Image processing tips for geoscientists – 1

  1. Here’s the ImageJ code for that. I like the stack presentation better than the montage for seeing the differences between the various kernel orientations.

    run(“Duplicate…”, “title=NW”);
    run(“Northwest”);
    selectWindow(“Galena.jpg”);
    run(“Duplicate…”, “title=N”);
    run(“North”);
    selectWindow(“Galena.jpg”);
    run(“Duplicate…”, “title=NE”);
    run(“Northeast”);
    selectWindow(“Galena.jpg”);
    run(“Duplicate…”, “title=W”);
    run(“West”);
    selectWindow(“Galena.jpg”);
    run(“Duplicate…”, “title=E”);
    run(“East”);
    selectWindow(“Galena.jpg”);
    run(“Duplicate…”, “title=SW”);
    run(“Southwest”);
    selectWindow(“Galena.jpg”);
    run(“Duplicate…”, “title=S”);
    run(“South”);
    selectWindow(“Galena.jpg”);
    run(“Duplicate…”, “title=SE”);
    run(“Southeast”);
    run(“Tile”);
    run(“Images to Stack”, “name=[Stop Animation with Back Slash] title=[] use”);
    run(“Make Montage…”, “columns=3 rows=3 scale=0.5 first=1 last=9 increment=1 border=0 font=12 label”);
    selectWindow(“Stop Animation with Back Slash”);
    doCommand(“Start Animation [\\]”);

    • Thank you for your comment Ron.
      And for the script. This is a great way for me to get started with code in ImageJ, so far I’ve only used plugins. Look forward to trying it.
      Cheers
      Matteo

      • Matteo

        save the code to a text file
        remove my comments at the top, or comment them out with //
        drag the text file and drop it on the ImageJ interface…an editor will open it.
        open your image that you would like to process in imageJ from the File menu or drag and drop the image
        select Macro>Run Macros, or press Ctrl+R in the editor

        Ron

      • Great
        Now I only need to find some time to try it. I may also look for a great new image to work with, maybe a medical image or some forensic stuff.
        Thanks Ron

  2. Pingback: Visualization tips for geoscientists: Surfer «·

  3. Pingback: Visualization tips for geoscientists – Matlab «·

  4. Pingback: Visualization tips for geoscientists: Matlab, part II «·

  5. Pingback: Visualization tips for geoscientists: Matlab, part III « MyCarta·

  6. Pingback: Visualization tips for geoscientists: Matlab, part III | MyCarta·

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s