-
Notifications
You must be signed in to change notification settings - Fork 0
/
fortfilt.f90
194 lines (153 loc) · 6.75 KB
/
fortfilt.f90
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
185
186
187
188
189
190
191
192
193
194
module fortfilt
implicit none
private ! everything private by default
public :: naiveGauss
public :: fastGauss
public :: BoxBlur
public :: naiveBox
contains
subroutine naiveGauss (source, filtered, r)
implicit none
integer, intent(in) :: r
double precision, intent(in) :: source(:,:)
double precision, intent(out) :: filtered(:,:)
double precision, parameter :: PI = 4.*atan(1.)
integer :: i, j, k, l
integer :: ii, jj ! inside the array
integer :: nx, ny
integer :: dim(2)
double precision :: val, wsum
double precision :: dsq ! distance squared
double precision :: wght ! weight
dim = shape(source)
nx = dim(1)
ny = dim(2)
do i = 1, nx
do j = 1, ny
val = 0
wsum = 0
do k = i-r, i+r ! inner loop over kernel
do l = j-r, j+r ! inner loop over kernel
ii = min(nx, max(1,k)) ! make sure i is always inside the grid (this implies values are extendet (stretched at the boundaries))
jj = min(ny, max(1,l)) ! make sure j is always inside the grid (this implies values are extendet (stretched at the boundaries))
dsq = (j-l)**2 + (i-k)**2
wght = exp(-dsq / (2.d0*r**2)) / (2.d0*PI*r**2)
val = val + source(ii,jj) * wght
wsum = wsum + wght
! print *, i,j, k, l, ii, jj, dsq
end do
end do
filtered(i,j) = val / wsum
end do
end do
end subroutine naiveGauss
subroutine BoxBlurH (source, filtered, r)
! computes horizontal blur
implicit none
integer, intent(in) :: r
double precision, intent(in) :: source(:,:)
double precision, intent(out) :: filtered(:,:)
integer :: nx, ny
integer :: dim(2)
double precision :: wght ! weight
double precision :: sum
integer :: i,j
integer :: il ! leftmost pixel which should be removed from the accumulator
integer :: ir ! rightmost pixel which should be added to the accumulator
dim = shape(source)
nx = dim(1)
ny = dim(2)
wght = 1.d0 / (2.d0*r+1.d0)
do j = 1, ny ! loop over all rows
! compute sum at first pixel
sum = source(1,j)
do i = 1, r ! loop over box kernel
sum = sum + source(1,j) + source(i+1,j) ! always take 1 as left pixel, to not get out of grid
! print *, sum
end do
! generate output pixel, then update running sum
do i = 1, nx
! print *, j, i, sum
filtered(i,j) = sum * wght
il = max(i-r, 1) ! make sure we dont get off the grid
ir = min(i+r+1, nx) ! make sure we dont get off the grid
sum = sum + source(ir,j) - source(il,j)
end do
end do
end subroutine BoxBlurH
subroutine BoxBlur (source, filtered, r)
! computes box blur, by calling horizontal boxblur twice, the second time with tansposed matrix
implicit none
integer, intent(in) :: r
double precision, intent(in) :: source(:,:)
double precision, intent(out) :: filtered(:,:)
integer :: nx, ny
integer :: dim(2)
double precision, allocatable :: source_t(:,:) ! source transposed
double precision, allocatable :: filtered_t(:,:) ! filtered transposed
dim = shape(source)
nx = dim(1)
ny = dim(2)
allocate(source_t(ny,nx))
allocate(filtered_t(ny,nx))
! first horizontal blur
call BoxBlurH (source, filtered, r)
! then transpose result and call horizontal blur again, now really doing vertical blur
source_t = transpose (filtered)
call BoxBlurH (source_t, filtered_t, r)
! transpose again, to get back initial shape
filtered = transpose (filtered_t)
end subroutine BoxBlur
subroutine fastGauss (source, filtered, r)
! computes a fast approximation to a Gaussian filter, by applying a series of box filters
implicit none
integer, intent(in) :: r
double precision, intent(in) :: source(:,:)
double precision, intent(out) :: filtered(:,:)
integer :: nx, ny
integer :: dim(2)
double precision, allocatable :: tmpfilter1(:,:) ! tmp array to store intermediate results
double precision, allocatable :: tmpfilter2(:,:) ! tmp array to store intermediate results
dim = shape(source)
nx = dim(1)
ny = dim(2)
allocate(tmpfilter1(nx,ny))
allocate(tmpfilter2(nx,ny))
call BoxBlur (source, tmpfilter1, r)
call BoxBlur (tmpfilter1, tmpfilter2, r)
call BoxBlur (tmpfilter2, filtered, r)
end subroutine fastGauss
subroutine naiveBox(source, filtered, r)
! computes a box filter as in the original code from Ann Lebrocq to compare against
implicit none
integer, intent(in) :: r
double precision, intent(in) :: source(:,:)
double precision, intent(out) :: filtered(:,:)
integer :: nx, ny
integer :: dim(2)
integer :: i,j,ii,jj !,nlen,xlen,ylen,jjl,jjh,iil,iih
double precision :: dsum, dcount
dim = shape(source)
nx = dim(1)
ny = dim(2)
do j=r+1,ny-r-1
do i=r+1,nx-r-1
dsum = 0.0
dcount = 0.0
if ( source(i,j) >= 0 ) then ! was imask before, but we want to avoid additional arrays for comparison
do jj=j-r,j+r
do ii=i-r,i+r
if ( source(ii,jj) >= 0) then
dsum = dsum + source(ii,jj)
dcount = dcount + 1.0
filtered(i,j) = dsum/dcount
end if
end do
end do
else
filtered(i,j) = source(i,j)
end if
end do
end do
end subroutine naiveBox
end module fortfilt