forked from avikhlinin/wvdecomp
-
Notifications
You must be signed in to change notification settings - Fork 0
/
gausswv.f
190 lines (154 loc) · 4.14 KB
/
gausswv.f
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
174
175
176
177
178
179
180
181
182
183
184
subroutine gausswv (img,z,bg,nx,ny,scale)
c
implicit none
integer nx,ny
real img(nx,ny)
real z(nx,ny)
real bg(nx,ny)
integer scale
c real w
integer w11(7),w12(7),w13(7),w21(7),w22(7),w23(7)
integer ww11,ww12,ww13,ww21,ww22,ww23
integer iscale
c scale 1 2 3 4 5 6 7
c sig1 0.0 1.0 2.0 4.0 8.0 16. 32.
c sig2 1.0 2.0 4.0 8.0 16. 32. 64.
data w11 / 1, 3, 3, 9, 15, 31, 65 /
data w12 / 1, 3, 5, 7, 17, 35, 53 /
data w13 / 1, 1, 3, 7, 15, 31, 71 /
data w21 / 3, 3, 9, 15, 31, 65, 127 /
data w22 / 3, 5, 7, 17, 35, 53, 123 /
data w23 / 1, 3, 7, 15, 31, 71, 139 /
integer wvecx(3),wvecy(3)
call mcopy(z,img,nx,ny)
call mcopy(bg,z,nx,ny)
if (scale.le.7) then
ww11 = w11(scale)
ww12 = w12(scale)
ww13 = w13(scale)
ww21 = w21(scale)
ww22 = w22(scale)
ww23 = w23(scale)
else
ww11 = w11(7)
ww12 = w12(7)
ww13 = w13(7)
ww21 = w21(7)
ww22 = w22(7)
ww23 = w23(7)
do iscale=8,scale
ww11 = ww11 * 2
ww12 = ww12 * 2
ww13 = ww13 * 2
ww21 = ww21 * 2
ww22 = ww22 * 2
ww23 = ww23 * 2
enddo
ww11 = ww11 - 1
ww12 = ww12 - 1
ww13 = ww13 - 1
ww21 = ww21 - 1
ww22 = ww22 - 1
ww23 = ww23 - 1
endif
wvecx(1)=ww11
wvecx(2)=ww12
wvecx(3)=ww13
wvecy(1)=ww11
wvecy(2)=ww12
wvecy(3)=ww13
call conv_rect1_n(z,nx,ny,wvecx,wvecy,3)
* w=float(ww11**2)*float(ww12**2)*float(ww13**2)
* w=1.0/w
* call msarith (z,"*",w,nx*ny)
wvecx(1)=ww21
wvecx(2)=ww22
wvecx(3)=ww23
wvecy(1)=ww21
wvecy(2)=ww22
wvecy(3)=ww23
call conv_rect1_n(bg,nx,ny,wvecx,wvecy,3)
* w=float(ww21**2)*float(ww22**2)*float(ww23**2)
* w=1.0/w
* call msarith (bg,"*",w,nx*ny)
return
end
subroutine conv_rect1_n(d,nx,ny,wvecx,wvecy,nconv)
implicit none
integer nx,ny
real d(nx,ny)
integer nconv
integer wvecx(nconv),wvecy(nconv)
integer nmax
parameter (nmax=16384)
double precision row(nmax), drow(nmax)
integer i,j
double precision w
integer iconv,nrecx,nrecy
c==Checks
do iconv=1,nconv
nrecx=wvecx(iconv)
nrecy=wvecy(iconv)
if ((nrecx/2)*2.eq.nrecx.or.(nrecy/2)*2.eq.nrecy) then
print*,iconv,nrecx,nrecy
call exiterror('must be odd wvec in conv_rect1_n')
endif
enddo
c== Columns
do i=1,nx
c..copy column
do j=1,ny
row(j)=d(i,j)
drow(j)=row(j)
enddo
do iconv = 1,nconv
nrecy = wvecy(iconv)
call conv_rect1_1D_d (row,drow,ny,nrecy)
c..Copy the result of convolution to the 1-D data
if (iconv.lt.nconv) then
do j=1,ny
row(j)=drow(j)
enddo
endif
enddo
c..Copy results back to the array
do j=1,ny
d(i,j)=sngl(drow(j))
enddo
enddo
c== Rows
do j=1,ny
c..copy row
do i=1,nx
row(i)=dble(d(i,j))
drow(i)=row(i)
enddo
do iconv = 1,nconv
nrecx = wvecx(iconv)
call conv_rect1_1D_d (row,drow,nx,nrecx)
c..Copy the result of convolution to the 1-D data
if (iconv.lt.nconv) then
do i=1,nx
row(i)=drow(i)
enddo
endif
enddo
c..Copy results back to the array
do i=1,nx
d(i,j)=sngl(drow(i))
enddo
enddo
c==Normalize
w=wvecx(1)*wvecy(1)
do iconv=2,nconv
w = w*wvecx(iconv)*wvecy(iconv)
enddo
w=1.0/w
do i=1,nx
do j=1,ny
d(i,j)=sngl(d(i,j)*w)
enddo
enddo
return
end
*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx