Detection of geological faults
3 views (last 30 days)
Show older comments
Hello. I am trying to use this code in a satellite image (SRTM DEM) for detection of geolofical faults:
close all;
clear all;
%%
%Automatic enhancement and identification of geological fault
%
%This code is a simple(not speed optimized) implementation of fault detection method based on the paper [1], [2].
%If you want to use the software for commercial purpose you have to contact with the authors of [1], [2].
%You can find more details in https://sites.google.com/site/costaspanagiotakis/research/faults
%We will appreciate if you cite our papers [1] and [2] in your work:
%[1] C. Panagiotakis and E. Kokinou, Automatic enhancement and detection of active sea faults from bathymetry,International
%Conference on Pattern Recognition, 2014.
%[2] C. Panagiotakis and E. Kokinou, Linear Pattern Detection of Geological Faults via a Topology and Shape Optimization Method,
%IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of Selected Topics in Applied Earth Observations and
%Remote Sensing, vol. 8, no. 1, pp. 3-11, 2015.
%%
%Reading data variables
%DataSR: Bathymatric image
%Y_cord,X_cord: Coordinates of the map
load('data.mat');
figure;
colormap('gray');
imagesc(Y_cord,X_cord,DataSR);
title('Selected Area');
axis xy;
fsize = 8;
H = fspecial('gaussian', fsize, 2);
blurredDataSR = imfilter(DataSR,H,'replicate');
[Asp,Slope] = getSlopeAndAspect(blurredDataSR,5,5);
[~,DSlope] = getSlopeAndAspect(Slope,5,5);
[DAspect] = getSlopeOfAspect(Asp,5,5);
FDSlope = imfilter(DSlope,H,'replicate');
FDAspect = imfilter(DAspect,H,'replicate');
%Fault Enhancement image
I = (Slope.^2.*FDSlope.^1.*FDAspect.^1).^(1/4);
figure;
imagesc(Y_cord,X_cord,Slope);
title('Slope');
axis xy;
figure;
imagesc(Y_cord,X_cord,Asp);
title('Aspect');
axis xy;
figure;
imagesc(Y_cord,X_cord,I);
colormap('gray');
title('Fault Enhancement image');
axis xy;
%parameters for step Filtering
grammes = 20;
Step_Filter_step = pi/18; %angle step
Step_Filter_DTheta = [0:Step_Filter_step:pi-(Step_Filter_step)];
Step_Filter_Dbw = [2 4 6]; %vector of width of main wave
I = Slope;
%step Filtering
[Y]= step_Filtering(I,grammes, Step_Filter_DTheta,Step_Filter_Dbw,1);
%detection of all faults from Y - BW may contains fault faults
[BW,plotSkel] = getFaultDetection(Y,DataSR);
%detection of strong faults
BWorig = BW;
C = bwconncomp(BW);
RS = regionprops(C,'Area','Orientation','Centroid');
L = length(RS);
area = [RS(1:L).Area];
Map = bwlabel(BW, 8);
RS_Y = regionprops(C,Y,'PixelValues');
feat = zeros(L,2);%features
for i=1:L,
vec = RS_Y(i).PixelValues;
feat(i,1) = sum(vec);
end
feat(:,2) = feat(:,1);
sv = sort(feat(:,1),'descend');
%
% parameter for selection the NUMBER OF Faults
%
NUMBER_OF_Faults = round(0.3*length(sv));
Thresh = sv(NUMBER_OF_Faults);
Labels = find(feat(:,1) >= Thresh);
%Plots Detected regions on Depth
[plotH,~] = plotNewColorMap2(Y,DataSR,Map,Labels,'Strong faults',Y_cord,X_cord);
After I tried to run the code, that was the result:
Error using load
Unable to read file 'DEM_clip_01.tif'. Input must be a MAT-file or an ASCII file containing numeric data with same number of columns in each row.
Error in untitled (line 19)
load('DEM_clip_01.tif');
0 Comments
Answers (0)
See Also
Categories
Find more on Geology in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!