' 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