# Image processing tips for geoscientists – 1

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

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
%
%
%% 1 - dealing with input
% once imported jpeg file as galenaAlaska
%
in=zeros(1094,1664);
n=1:1:1094; k=1:1:1664;
%
%
%% 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); 