-
-
Notifications
You must be signed in to change notification settings - Fork 2
/
Mparsmoother.c
173 lines (149 loc) · 5.72 KB
/
Mparsmoother.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
/* CRS Parameter smoother - Use moving average in a user defined window in space and time
This program is a parameter smoother for the parameter data cube generated by sfvfsacrsnh. It uses the parameter data cube organized in (t0,m0,parameter) dimensions. The parameters are stored in RN, RNIP, BETA, Semblance order for each (t0, m0) pair, and they are obtained using VFSA global optimization.
So, this program moves a space x time window in each parameter section and it calculates the average for data samples into this window. The value at the window center will be the average.
This program uses the following smoothing algorithm (Mann and Duveneck, 2003):
The basic processing steps to calculate the smoothed wavefield attributes for a given ZO location read as
follows:
- reject all points within the window with S < Smin to avoid the use of unreliable attributes
- reject all points within the window with dbeta > dbetaMax to avoid a mixing of attributes associated with independent reflection events
- separately compute the average of each attribute for a fraction f of the remaining data points centered around the respective median of its distribution.
If no data points remain for averaging after the application of the above-mentioned criteria, the original wavefield attributes are used. Note that for the normal wavefront.
*/
#include <rsf.h>
void sortingXinAscendingOrder(
float *x, /* x vector to sort */
int n /* Vectors dimension */)
/*< x vector sorting in ascending order >*/
{
int i; // Loop counter
float tmpx; // Temporary variables
int k; // Sorting key (number of changes)
do{
k=0;
for(i=1;i<n;i++){
if(x[i-1]>x[i]){
tmpx=x[i-1];
x[i-1]=x[i];
x[i]=tmpx;
k++;
}
} // Loop vector samples
}while(k!=0);
}
int main(int argc, char *argv[]){
float ***p; // Parameters data cube RN, RNIP, BETA, SEMBLANCE
int n1, n2, n3; // Parameters data cube dimension
float d1, d2, d3; // Parameters data cube sampling
float o1, o2, o3; // Parameters data cube axis origin
int rect[2]; // Window dimensions
char key[6]; // String to get parameters through cmd
int i; // Loop counter
float smin; // Minimum semblance
float dbeta; // Minimun delta beta
float f; // Control fraction
int ip, im0, it0, ix, it; // Loop counter
float sum; // Samples sum
float ***pp; // Parameters after smoothing
bool **mask; // Samples mask
float ***filemask;
int imm, itt; // Samples index in CMP and Time
int ns, nss; // Number of samples
float *v; // Samples vector to sorting and median calculation
float median; // Median of the distribution
sf_file in;
sf_file out;
sf_file mask_file;
sf_init(argc,argv);
in = sf_input("in");
out = sf_output("out");
mask_file = sf_output("mask");
for(i=0;i<2;i++){
snprintf(key,6,"rect%d",i+1);
if (!sf_getint(key,rect+i)) rect[i]=1;
/*( rect#=(1,1,...) smoothing radius on #-th axis )*/
}
if(!sf_getfloat("smin",&smin)) smin=0.;
/* Minimun semblance */
if(!sf_getfloat("dbeta",&dbeta)) dbeta=1.;
/* beta interval */
if(!sf_getfloat("f",&f)) f=1.;
/* Control fraction */
/* Read parameters file */
if(!sf_histint(in,"n1",&n1)) sf_error("No n1= in input file");
if(!sf_histfloat(in,"o1",&o1)) sf_error("No o1= in input file");
if(!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input file");
if(!sf_histint(in,"n2",&n2)) sf_error("No n2= in input file");
if(!sf_histfloat(in,"o2",&o2)) sf_error("No o2= in input file");
if(!sf_histfloat(in,"d2",&d2)) sf_error("No d2= in input file");
if(!sf_histint(in,"n3",&n3)) sf_error("No n3= in input file");
if(!sf_histfloat(in,"o3",&o3)) sf_error("No o3= in input file");
if(!sf_histfloat(in,"d3",&d3)) sf_error("No d3= in input file");
p = sf_floatalloc3(n1,n2,n3);
pp = sf_floatalloc3(n1,n2,n3);
filemask = sf_floatalloc3(n1,n2,n3);
sf_floatread(p[0][0],n1*n2*n3,in);
mask = sf_boolalloc2(rect[0],rect[1]);
// Smoothing
for(ip=0;ip<n3-1;ip++){
for(im0=0;im0<n2;im0++){
for(it0=0;it0<n1;it0++){
pp[ip][im0][it0] = p[ip][im0][it0];
filemask[ip][im0][it0]=1.;
sum = 0.; ns=0;
for(ix=0;ix<rect[1];ix++){
for(it=0;it<rect[0];it++){
imm = im0+ix-rect[1]/2; itt=it0+it-rect[0]/2;
if(imm >= 0 && imm<n2 && itt >= 0 && itt < n1 && p[3][imm][itt]>=smin && fabs(p[2][imm][itt]) <= dbeta){
ns++;
mask[ix][it]=true;
}else{
mask[ix][it]=false;
}
} // Loop window time
} // Loop window space
if(p[3][im0][it0]<smin) filemask[ip][im0][it0]=0.;
if(ns<=1){
pp[ip][im0][it0]=p[ip][im0][it0];
}else{
i=0;
v = sf_floatalloc(ns);
for(ix=0;ix<rect[1];ix++){
for(it=0;it<rect[0];it++){
imm = im0+ix-rect[1]/2; itt=it0+it-rect[0]/2;
if(mask[ix][it]){
v[i] = p[ip][imm][itt]; i++;
}
}
}
// sort v
sortingXinAscendingOrder(v,ns);
if(ns%2==0){
median = (v[ns/2]+v[ns/2-1])/2.;
}else{
median = v[ns/2];
}
nss=0;
for(i=(1-f)*ns/2;i<(1+f)*ns/2;i++){
sum += v[i]; nss++;
}
pp[ip][im0][it0] = sum/(nss);
free(v);
}
} // Loop for t0s
} // Loop for m0s
} // Loop for parameters
// Semblance should not be smoothed
// Just keep semblance values
for(ip=n3-1;ip<n3;ip++){
for(im0=0;im0<n2;im0++){
for(it0=0;it0<n1;it0++){
pp[ip][im0][it0]=p[ip][im0][it0];
filemask[ip][im0][it0]=1.;
}
}
}
//sf_setformat(mask_file,"native_int");
//sf_putint(mask_file,"n3",1);
sf_floatwrite(filemask[0][0],n1*n2*n3,mask_file);
sf_floatwrite(pp[0][0],n1*n2*n3,out);
}