00001 ##@package Msmoothgradient 00002 # Return gradient of input image with smoothing orthogonal to gradient direction 00003 # 00004 # @code 00005 # import rsf as sf 00006 # import numpy as np 00007 # @endcode 00008 # Import RSF and NumPy packages as they will be used later in the code 00009 # 00010 # @code 00011 # from Sadiabatic import * 00012 # @endcode 00013 # Sadiabatic is a package containing a function of the same name used for 00014 # padding the boundaries of the image with repetition (Neumann boundary conditions) 00015 # 00016 # Get the RSF parameters for this function 00017 # @code 00018 # par = sf.Par() 00019 # @endcode 00020 # 00021 # Get the name of the file for which we want to calculate the gradients from the RSF input. 00022 # As well, output1 is the name of the RSF file to which we write the gradient values in the 00023 # x direction, and output2 is the RSF filename for writing the gradient values in the y 00024 # direction. 00025 # @code 00026 # input = sf.Input() 00027 # output1 = sf.Input(par.string("dx_outputfilename")) 00028 # output2 = sf.Input(par.string("dy_outputfilename")) 00029 # @endcode 00030 # 00031 # n1 is the width of the image and n2 is the height of the image 00032 # @code 00033 # n1 = input.int("n1") 00034 # n2 = input.int("n2") 00035 # @endcode 00036 # 00037 # Make sure the input is 2-D. Otherwise, the gradients cannot be calculated with this code. 00038 # @code 00039 # ni = input.size(2) 00040 # assert ni == 1,"sfsmoothgradient needs 2D input" 00041 # @endcode 00042 # 00043 # Read in the input image 00044 # @code 00045 # inputimage = np.zeros((n1,n2),'f') 00046 # input.read(inputimage) 00047 # @endcode 00048 # 00049 # Pad the input image with repetition so that gradients can be calculated on the boundaries. 00050 # @code 00051 # input_padded = Sadiabatic(inputimage) 00052 # @endcode 00053 # 00054 # Used for indexing of image for computation of gradient 00055 # @code 00056 # I = array(range(1,n1+1)) 00057 # J = array(range(1,n2+1)) 00058 # @endcode 00059 # 00060 # Calculate gradient in x direction 00061 # @code 00062 # Imx = (3*input_padded[I+1][:,J-1]+10*input_padded[I+1][:,J]+3*input_padded[I+1][:,J+1]-3*input_padded[I-1][:,J-1]-10*input_padded[I-1][:,J]-3*input_padded[I-1][:,J+1])/32 00063 # @endcode 00064 # 00065 # Calculate gradient in y direction 00066 # @code 00067 # Imy = (3*input_padded[I-1][:,J+1]+10*input_padded[I][:,J+1]+3*input_padded[I+1][:,J+1]-3*input_padded[I-1][:,J-1]-10*input_padded[I][:,J-1]-3*input_padded[I+1][:,J-1])/32 00068 # @endcode 00069 # 00070 # Write directional gradients (in x and y directions) to appropriate output files 00071 # @code 00072 # output1.put("n1",n1) 00073 # output1.put("n2",n2) 00074 # output1.write(Imx) 00075 # 00076 # output2.put("n1",n1) 00077 # output2.put("n2",n2) 00078 # output2.write(Imy) 00079 # @endcode 00080 00081 #!/usr/bin/env python 00082 # Author: R. Shahidi 00083 # Seismic Laboratory for Imaging and Modeling 00084 # Department of Earch & Ocean Sciences 00085 # The University of British Columbia 00086 # 00087 # Date : July, 2008 00088 00089 # Copyright (C) 2008 The University of British Columbia at Vancouver 00090 # 00091 # This program is free software; you can redistribute it and/or modify 00092 # it under the terms of the GNU General Public License as published by 00093 # the Free Software Foundation; either version 2 of the License, or 00094 # (at your option) any later version. 00095 # 00096 # This program is distributed in the hope that it will be useful, 00097 # but WITHOUT ANY WARRANTY; without even the implied warranty of 00098 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00099 # GNU General Public License for more details. 00100 # 00101 # You should have received a copy of the GNU General Public License 00102 # along with this program; if not, write to the Free Software 00103 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00104 00105 import rsf as sf 00106 import numpy as np 00107 import sys 00108 from Sadiabatic import * 00109 00110 par = sf.Par() 00111 00112 input = sf.Input() 00113 print>>sys.stderr,par 00114 print>>sys.stderr,par.string("dx_outputfilename") 00115 print>>sys.stderr,par.string("dy_outputfilename") 00116 output1 = sf.Output(par.string("dx_outputfilename")) 00117 output2 = sf.Output(par.string("dy_outputfilename")) 00118 00119 n1 = input.int("n1") 00120 n2 = input.int("n2") 00121 ni = input.size(2) 00122 assert ni == 1,"sfsmoothgradient needs 2D input" 00123 00124 inputimage = np.zeros((n1,n2),'f') 00125 input.read(inputimage) 00126 input_padded = Sadiabatic(inputimage) 00127 00128 # Used for indexing of image for computation of gradient 00129 I = array(range(1,n1+1)) 00130 J = array(range(1,n2+1)) 00131 00132 # Gradient in x direction 00133 Imx = (3*input_padded[I+1][:,J-1]+10*input_padded[I+1][:,J]+3*input_padded[I+1][:,J+1]-3*input_padded[I-1][:,J-1]-10*input_padded[I-1][:,J]-3*input_padded[I-1][:,J+1])/32 00134 00135 # Gradient in y direction 00136 Imy = (3*input_padded[I-1][:,J+1]+10*input_padded[I][:,J+1]+3*input_padded[I+1][:,J+1]-3*input_padded[I-1][:,J-1]-10*input_padded[I][:,J-1]-3*input_padded[I+1][:,J-1])/32 00137 00138 # Write directional gradients to output files 00139 output1.put("n1",n1) 00140 output1.put("n2",n2) 00141 output1.write(Imx) 00142 00143 output2.put("n1",n1) 00144 output2.put("n2",n2) 00145 output2.write(Imy)
1.5.6