Detection of geological faults

2 visualizzazioni (ultimi 30 giorni)
Arlete Conde
Arlete Conde il 25 Lug 2022
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');

Risposte (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by