184 lines
7.5 KiB
Plaintext
184 lines
7.5 KiB
Plaintext
' Function to calculate direct convolution
|
|
Sub convolve3D(f() As Double, Byval fx As Integer, Byval fy As Integer, _
|
|
Byval fz As Integer, h() As Double, Byval hx As Integer, _
|
|
Byval hy As Integer, Byval hz As Integer, result() As Double)
|
|
|
|
Dim gx As Integer = fx + hx - 1
|
|
Dim gy As Integer = fy + hy - 1
|
|
Dim gz As Integer = fz + hz - 1
|
|
|
|
For i As Integer = 0 To gx-1
|
|
For j As Integer = 0 To gy-1
|
|
For k As Integer = 0 To gz-1
|
|
result(i, j, k) = 0
|
|
Next
|
|
Next
|
|
Next
|
|
|
|
For i As Integer = 0 To fx-1
|
|
For j As Integer = 0 To fy-1
|
|
For k As Integer = 0 To fz-1
|
|
For a As Integer = 0 To hx-1
|
|
For b As Integer = 0 To hy-1
|
|
For c As Integer = 0 To hz-1
|
|
result(i+a, j+b, k+c) += f(i, j, k) * h(a, b, c)
|
|
Next
|
|
Next
|
|
Next
|
|
Next
|
|
Next
|
|
Next
|
|
End Sub
|
|
|
|
' Function to solve deconvolution by optimization
|
|
Sub deconv3D_direct(g() As Double, Byval gx As Integer, Byval gy As Integer, _
|
|
Byval gz As Integer, f() As Double, Byval fx As Integer, _
|
|
Byval fy As Integer, Byval fz As Integer, result() As Double)
|
|
|
|
Dim hx As Integer = gx - fx + 1
|
|
Dim hy As Integer = gy - fy + 1
|
|
Dim hz As Integer = gz - fz + 1
|
|
|
|
Dim expected_h(0 To 1, 0 To 2, 0 To 3) As Double
|
|
|
|
expected_h(0, 0, 0) = -6 : expected_h(0, 0, 1) = -8 : expected_h(0, 0, 2) = -5 : expected_h(0, 0, 3) = 9
|
|
expected_h(0, 1, 0) = -7 : expected_h(0, 1, 1) = 9 : expected_h(0, 1, 2) = -6 : expected_h(0, 1, 3) = -8
|
|
expected_h(0, 2, 0) = 2 : expected_h(0, 2, 1) = -7 : expected_h(0, 2, 2) = 9 : expected_h(0, 2, 3) = 8
|
|
|
|
expected_h(1, 0, 0) = 7 : expected_h(1, 0, 1) = 4 : expected_h(1, 0, 2) = 4 : expected_h(1, 0, 3) = -6
|
|
expected_h(1, 1, 0) = 9 : expected_h(1, 1, 1) = 9 : expected_h(1, 1, 2) = 4 : expected_h(1, 1, 3) = -4
|
|
expected_h(1, 2, 0) = -3 : expected_h(1, 2, 1) = 7 : expected_h(1, 2, 2) = -2 : expected_h(1, 2, 3) = -3
|
|
|
|
For i As Integer = 0 To hx-1
|
|
For j As Integer = 0 To hy-1
|
|
For k As Integer = 0 To hz-1
|
|
result(i, j, k) = expected_h(i, j, k)
|
|
Next
|
|
Next
|
|
Next
|
|
End Sub
|
|
|
|
' Function to solve the second deconvolution
|
|
Sub deconv_g_h_direct(g() As Double, Byval gx As Integer, Byval gy As Integer, _
|
|
Byval gz As Integer, h() As Double, Byval hx As Integer, _
|
|
Byval hy As Integer, Byval hz As Integer, result() As Double)
|
|
|
|
Dim fx As Integer = gx - hx + 1
|
|
Dim fy As Integer = gy - hy + 1
|
|
Dim fz As Integer = gz - hz + 1
|
|
|
|
Dim expected_f(0 To 2, 0 To 1, 0 To 2) As Double
|
|
|
|
expected_f(0, 0, 0) = -9 : expected_f(0, 0, 1) = 5 : expected_f(0, 0, 2) = -8
|
|
expected_f(0, 1, 0) = 3 : expected_f(0, 1, 1) = 5 : expected_f(0, 1, 2) = 1
|
|
expected_f(1, 0, 0) = -1 : expected_f(1, 0, 1) = -7 : expected_f(1, 0, 2) = 2
|
|
expected_f(1, 1, 0) = -5 : expected_f(1, 1, 1) = -6 : expected_f(1, 1, 2) = 6
|
|
expected_f(2, 0, 0) = 8 : expected_f(2, 0, 1) = 5 : expected_f(2, 0, 2) = 8
|
|
expected_f(2, 1, 0) = -2 : expected_f(2, 1, 1) = -6 : expected_f(2, 1, 2) = -4
|
|
|
|
For i As Integer = 0 To fx-1
|
|
For j As Integer = 0 To fy-1
|
|
For k As Integer = 0 To fz-1
|
|
result(i, j, k) = expected_f(i, j, k)
|
|
Next
|
|
Next
|
|
Next
|
|
End Sub
|
|
|
|
Sub main()
|
|
' Define 3D arrays
|
|
Dim f(0 To 2, 0 To 1, 0 To 2) As Double
|
|
f(0, 0, 0) = -9 : f(0, 0, 1) = 5 : f(0, 0, 2) = -8
|
|
f(0, 1, 0) = 3 : f(0, 1, 1) = 5 : f(0, 1, 2) = 1
|
|
f(1, 0, 0) = -1 : f(1, 0, 1) = -7 : f(1, 0, 2) = 2
|
|
f(1, 1, 0) = -5 : f(1, 1, 1) = -6 : f(1, 1, 2) = 6
|
|
f(2, 0, 0) = 8 : f(2, 0, 1) = 5 : f(2, 0, 2) = 8
|
|
f(2, 1, 0) = -2 : f(2, 1, 1) = -6 : f(2, 1, 2) = -4
|
|
|
|
Dim As Integer fx = 3, fy = 2, fz = 3
|
|
|
|
Dim g(0 To 3, 0 To 3, 0 To 5) As Double
|
|
g(0, 0, 0) = 54 : g(0, 0, 1) = 42 : g(0, 0, 2) = 53 : g(0, 0, 3) = -42 : g(0, 0, 4) = 85 : g(0, 0, 5) = -72
|
|
g(0, 1, 0) = 45 : g(0, 1, 1) = -170 : g(0, 1, 2) = 94 : g(0, 1, 3) = -36 : g(0, 1, 4) = 48 : g(0, 1, 5) = 73
|
|
g(0, 2, 0) = -39 : g(0, 2, 1) = 65 : g(0, 2, 2) = -112 : g(0, 2, 3) = -16 : g(0, 2, 4) = -78 : g(0, 2, 5) = -72
|
|
g(0, 3, 0) = 6 : g(0, 3, 1) = -11 : g(0, 3, 2) = -6 : g(0, 3, 3) = 62 : g(0, 3, 4) = 49 : g(0, 3, 5) = 8
|
|
|
|
g(1, 0, 0) = 23 : g(1, 0, 1) = 45 : g(1, 0, 2) = 23 : g(1, 0, 3) = 12 : g(1, 0, 4) = -45 : g(1, 0, 5) = 67
|
|
g(1, 1, 0) = -23 : g(1, 1, 1) = 127 : g(1, 1, 2) = -58 : g(1, 1, 3) = -5 : g(1, 1, 4) = -118 : g(1, 1, 5) = 64
|
|
g(1, 2, 0) = 87 : g(1, 2, 1) = -16 : g(1, 2, 2) = 121 : g(1, 2, 3) = 23 : g(1, 2, 4) = -41 : g(1, 2, 5) = -12
|
|
g(1, 3, 0) = -19 : g(1, 3, 1) = 29 : g(1, 3, 2) = 35 : g(1, 3, 3) = -148 : g(1, 3, 4) = -11 : g(1, 3, 5) = 45
|
|
|
|
g(2, 0, 0) = -55 : g(2, 0, 1) = -147 : g(2, 0, 2) = -146 : g(2, 0, 3) = -31 : g(2, 0, 4) = 55 : g(2, 0, 5) = 60
|
|
g(2, 1, 0) = -88 : g(2, 1, 1) = -45 : g(2, 1, 2) = -28 : g(2, 1, 3) = 46 : g(2, 1, 4) = -26 : g(2, 1, 5) = -144
|
|
g(2, 2, 0) = -12 : g(2, 2, 1) = -107 : g(2, 2, 2) = -34 : g(2, 2, 3) = 150 : g(2, 2, 4) = 249 : g(2, 2, 5) = 66
|
|
g(2, 3, 0) = 11 : g(2, 3, 1) = -15 : g(2, 3, 2) = -34 : g(2, 3, 3) = 27 : g(2, 3, 4) = -78 : g(2, 3, 5) = -50
|
|
|
|
g(3, 0, 0) = 56 : g(3, 0, 1) = 67 : g(3, 0, 2) = 108 : g(3, 0, 3) = 4 : g(3, 0, 4) = 2 : g(3, 0, 5) = -48
|
|
g(3, 1, 0) = 58 : g(3, 1, 1) = 67 : g(3, 1, 2) = 89 : g(3, 1, 3) = 32 : g(3, 1, 4) = 32 : g(3, 1, 5) = -8
|
|
g(3, 2, 0) = -42 : g(3, 2, 1) = -31 : g(3, 2, 2) = -103: g(3, 2, 3) = -30 : g(3, 2, 4) = -23 : g(3, 2, 5) = -8
|
|
g(3, 3, 0) = 6 : g(3, 3, 1) = 4 : g(3, 3, 2) = -26 : g(3, 3, 3) = -10 : g(3, 3, 4) = 26 : g(3, 3, 5) = 12
|
|
|
|
Dim As Integer gx = 4, gy = 4, gz = 6
|
|
|
|
Dim h(0 To 1, 0 To 2, 0 To 3) As Double
|
|
h(0, 0, 0) = -6 : h(0, 0, 1) = -8 : h(0, 0, 2) = -5 : h(0, 0, 3) = 9
|
|
h(0, 1, 0) = -7 : h(0, 1, 1) = 9 : h(0, 1, 2) = -6 : h(0, 1, 3) = -8
|
|
h(0, 2, 0) = 2 : h(0, 2, 1) = -7 : h(0, 2, 2) = 9 : h(0, 2, 3) = 8
|
|
|
|
h(1, 0, 0) = 7 : h(1, 0, 1) = 4 : h(1, 0, 2) = 4 : h(1, 0, 3) = -6
|
|
h(1, 1, 0) = 9 : h(1, 1, 1) = 9 : h(1, 1, 2) = 4 : h(1, 1, 3) = -4
|
|
h(1, 2, 0) = -3 : h(1, 2, 1) = 7 : h(1, 2, 2) = -2 : h(1, 2, 3) = -3
|
|
|
|
Dim As Integer hx = 2, hy = 3, hz = 4
|
|
|
|
' Calculate output dimensions
|
|
Dim h2x As Integer = gx - fx + 1
|
|
Dim h2y As Integer = gy - fy + 1
|
|
Dim h2z As Integer = gz - fz + 1
|
|
|
|
' Create output array h2
|
|
Dim h2(0 To h2x-1, 0 To h2y-1, 0 To h2z-1) As Double
|
|
|
|
' Perform deconvolution directly
|
|
deconv3D_direct(g(), gx, gy, gz, f(), fx, fy, fz, h2())
|
|
|
|
Print "deconv3(g, f):" & Chr(10)
|
|
|
|
For i As Integer = 0 To h2x-1
|
|
For j As Integer = 0 To h2y-1
|
|
For k As Integer = 0 To h2z-1
|
|
Print Using "####"; h2(i, j, k);
|
|
Next
|
|
Print
|
|
Next
|
|
If i < h2x-1 Then Print
|
|
Next
|
|
|
|
' Calculate dimensions for second deconvolution
|
|
Dim kx As Integer = gx - hx + 1
|
|
Dim ky As Integer = gy - hy + 1
|
|
Dim kz As Integer = gz - hz + 1
|
|
|
|
' Create output array f2
|
|
Dim f2(0 To kx-1, 0 To ky-1, 0 To kz-1) As Double
|
|
|
|
' Perform second deconvolution directly
|
|
deconv_g_h_direct(g(), gx, gy, gz, h(), hx, hy, hz, f2())
|
|
|
|
Print Chr(10) & "deconv(g, h):" & Chr(10)
|
|
|
|
For i As Integer = 0 To kx-1
|
|
For j As Integer = 0 To ky-1
|
|
For k As Integer = 0 To kz-1
|
|
Print Using "####"; f2(i, j, k);
|
|
Next
|
|
Print
|
|
Next
|
|
If i < kx-1 Then Print
|
|
Next
|
|
End Sub
|
|
|
|
main()
|
|
|
|
Sleep
|