March 1, 2015 Programs No comments
sfthistogram computes a histogram for distribution of values in the input dataset.
The following example from rsf/rsf/sfnoise plots the histogram of a normally-distributed random noise: 
The output of sfhistogram contains integer values arranged in a one-dimensional array. The sampling is specified by n1=, d1=, and o1= parameters.
January 30, 2015 Programs No comments
sfmf applies a median filter on the first dimension of the input.
The size of the filter window is controlled by nfw= parameter. The following example from tccs/medianfilter/dragon shows field seismic data before and after 11-point (nfw=11) median filtering combined with bandpass filtering. 
For time-variant median filtering, see sftvmf. To simply output the median of the first axis, use sfmedian.
December 1, 2014 Programs No comments
sfbin bins traces with irregular spatial sampling to a regularly sampled 3-D cube.
The following example from sep/precon/cube shows the output of binning a common-offset seismic cube 
which corresponds to the following distribution of common midpoints 
Optionally, sfbin can also output the fold map (using fold= parameter). The fold map shows the number of input traces in each output bin. 
Parameters that control output grid sampling are nx=, dx=, x0= (for the second axis), ny=, dy=, y0= (for the third axis). Alternatively, one can specify the range values xmin=, xmax=, ymin=, ymax=.
By default, the range is determined from the input trace coordinates. The input trace coordinates can be specified in an auxiliary trace header file (head= parameter), where x and y coordinates are given in keys number xkey= and ykey= (0 and 1 by default).
By default, nearest-neighbor binning is applied. Alternatively, it is possible to use median binning by specifying interp=0 or bilinear-interpolation binning by specifying interp=2.
By default, the output values are normalized by the fold. To switch fold-normalization off, use norm=n.
November 12, 2014 Programs No comments
sfthreshold filters the input by soft thresholding (shrinkage).
Soft thresholding is a point-by-point operation, which can be described mathematically as
$T_{\mu}[u] = \left\{\begin{array}{rcl} u – \mu\,\mbox{sign}(u) & \quad & \mbox{if}\,|u| > \mu \\ 0 & \quad & \mbox{if}\,|u| \le \mu\end{array}\right.$
Soft thresholding was analyzed by Donoho (1995) and became particularly popular thanks to the iterative shrinkage-thresholding algorithm by Daubechies et al. (2004).
Donoho, D. L. (1995). De-noising by soft-thresholding. Information Theory, IEEE Transactions on, 41(3), 613-627.
Daubechies, I., Defrise, M., & De Mol, C. (2004). An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on pure and applied mathematics, 57(11), 1413-1457.
The following example from tccs/seislet/lena shows an image (Seismic Lena) and its reconstruction after soft thresholding in the seislet domain using 5% thresholding (pclip=5).
sfthreshold uses percentage parameter pclip= to set thresholding at the corresponding quantile of the data values. To do soft or hard thresholding with a fixed threshold, use sfthr.
An alternative thresholding-like operation is provided by sfsharpen.
October 8, 2014 Programs No comments
sfsigmoid generates a 2-D synthetic reflectivity model, created by Jon Claerbout.
One of the first occurrences of this model is in SEP-73 sponsor report from 1992, where it appeared in several papers:
J. F. Claerbout, 1992, Crossline Regridding by Inversion: SEP-73 report, 421-428, Stanford Exploration Project.
The model was described as “a synthetic model that illustrates local variation in bedding. Notice dipping bedding, curved bedding, unconformity between them, and a fault in the curved bedding.” Later, the sigmoid model made an appearance in Claerbout’s book Basic Earth Imaging. The following example from bei/krch/sep73 illustrates the effect of aliasing on Kirchhoff modeling and migration:
The model has appeared in numerous other tests. The following example from tccs/flat/flat shows automatic flattening of the sigmoid model by predictive painting.
sfsigmoid has several parameters that control the model. The usual n1=, n2=, o1=, o2=, d1=, d2= parameters control the mesh size and sampling, taper= indicates whether to taper the sides of the model, large= controls the length of the synthetic reflectivity series. The program takes no input.
September 24, 2014 Programs No comments
sfmax1 finds local maxima along the first axis of the input. It takes floating-point input but outputs complex numbers, where the real part stands for the location of the local minima and the imaginary part stands for the value of the input at local minima.
The number of minima to output is controlled by np= parameter. To control the range for the minima locations (in the case that it is smaller than the full range of the data), use min= and max=. The output is sorted by value so that the largest maxima appear first. Here is a quick example. Let us create some data:
bash$ sfmath n1=5 output="sin(x1)" > data.rsf
bash$ bash$ < data.rsf sfdisfil
0: 0 0.8415 0.9093 0.1411 -0.7568
Observing the data values, we can suspect that the local maximum is between 1 and 2.
sfmax1 uses local parabolic interpolation to locate the minimum at 1.581 with the value of 0.9826.
In the following example, from tccs/flat/flat, sfmax1 is used to locate the strongest-amplitude horizons for predictive painting.
August 3, 2014 Programs No comments
sfstolt implements zero-offset (post-stack) seismic migration using Stolt method.
Stolt migration is described in the classic paper
Stolt, R. H., 1978, Migration by Fourier transform: Geophysics, 43, 23-48.
The following example from gallery/french/stolt shows the result of Stolt migration in the French model:
The classic Stolt migration works for constant velocity. However, it can be extended to the case of V(z) by using Stolt stretch and cascaded migrations. See Evaluating the Stolt-stretch parameter and its references, including
Beasley, C., W. Lynn, K. Larner, and H. Nguyen, 1988, Cascaded frequency-wavenumber migration – Removing the restrictions on depth-varying velocity: Geophysics, 53, 881-893.
The following example from sep/stoltst/elfst compares the results of Stolt migration with Stolt stretch, phase-shift migration, and cascaded Stolt migration with Stolt stretch.
sfstolt2 is another version of Stolt migration, with a control on interpolation accuracy. The following example from sep/forwd/stolt compares results from two different interpolations:
July 13, 2014 Programs No comments
sfltft (Local Time-Frequency Transform) decomposes input data into frequency components.
The algorithm is described in the paper Seismic data analysis using local time-frequency decomposition and is based on regularized non-stationary regression. The following example from tccs/ltft/timefreq shows 1-D synthetic data composed of two chirp signals and the magnitude of coefficients in its time-frequency decomposition:
The frequency sampling in the output of sfltft is controled by nw=, w0=, and dw=. By default, these parameters correspond to the sampling of the discrete Fourier transform. The critical parameters for regularized regression are rect= (smoothing radius in time, in samples) and niter= (number of iterations). To output the details of iterative regularization, use verb=y. The frequency sampling and the rect= parameter provide explicit controls on time-frequency resolution. In the example above, rect=7. It is possible to change smoothing radius with frequency by using alpha= parameter. The iterative inversion can be controlled additionally by specifying a data weight (with mask=) or a model weight (with weight=). Optionally, the Fourier basis used in the decomposition can be extracted from the program by specifying basis= file.
To perform the inverse transform from time-frequency back to time domain, use inv=y. sfltft takes real-valued input and produces complex-valued output.
An analogous program for transforming complex-valued data is sfcltft. A related program is sftimefreq described in Time-frequency analysis of seismic data using local attributes.
June 11, 2014 Programs No comments
sfeikonal solves the eikonal equation using the Fast Marching Method. This computation produces first-arrival traveltimes on a fixed grid.
The following example from sep/fmeiko/fmarch shows traveltime contours for a point source inside the SEG/EAGE salt model.
The point source can be specified by its coordinates xshot=, yzhot=, and zshot= (note that the depth coordinate zshot corresponds to the first axis). In a small box around the source, the solution is computed analytically to avoid errors from the point-source singularity. The size of the box can be specified in samples (b1=, b2=, and b3=) or in physical dimensions (br1=, br2=, and br3=). For a plane-wave source instead of a point source, use plane1=y, plane2=y, or plane3=y. The plane is assumed to be aligned with the grid. For computing a traveltime table with multiple sources, the source coordinates can be specified in a file given by shotfile=. The order of accuracy in the finite-difference scheme is specified by order= parameter.
The following plot from sep/fmsec/cvel shows the error difference between the first- and second-order computations in a constant-velocity medium.
sfeikonal computes isotropic traveltimes. For an extension to VTI anisotropy, see sfeikonalvti.
May 13, 2014 Programs No comments
sfhelicon performs multidimensional convolution and inverse convolution (recursive filtering) using the helix transform.
The theory behind helical convolution is explained by Jon Claerbout in the paper
Claerbout, J., 1998, Multidimensional recursive filters via a helix: Geophysics, 63, 1532-1541
and in the corresponding chapter in his book Image Estimation by Example. The following example from gee/hlx/helicon illustrates inverse convolution on a helix.
sfhelicon works in N dimensions. The filter for is supplied by filt= parameter, which points to a real-valued RSF file with filter coefficents. The filter lags on a helix can be stored in a separate integer-value RSF file specified wih lag= parameter (which can be optionally contained inside the file with filter coefficients). The dimensions used for specifying the filter lags are not necessarily the same as the dimensions of the input data and can be specified with n= parameter (which can be optionally contained inside the file with filter lags).
The choice between forward and inverse convolution is controlled by div= parameter. The adjoint flag is supplied by adj= parameter. The following figure illustrates the act of the adjoint convolution and inverse convolution on a helix.
For an example of least-squares inversion with sfhelicon using the conjugate-gradient algorithm, see the documentation for sfconjgrad.
