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.

 0  0  0
 0  1  0
 0 -1  0

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