This task is a straightforward generalization of Deconvolution/1D to higher dimensions. For example, the one dimensional case would be applicable to audio signals, whereas two dimensions would pertain to images. Define the discrete convolution in 





d




{\displaystyle {\mathit {d}}}

 dimensions of two functions





H
,
F
:


Z


d


→

R



{\displaystyle H,F:\mathbb {Z} ^{d}\rightarrow \mathbb {R} }

 
taking 





d




{\displaystyle {\mathit {d}}}

-tuples of integers to real numbers as the function





G
:


Z


d


→

R



{\displaystyle G:\mathbb {Z} ^{d}\rightarrow \mathbb {R} }


also taking 





d




{\displaystyle {\mathit {d}}}

-tuples of integers to reals and satisfying





G
(

n

0


,
…
,

n

d
−
1


)
=

∑


m

0


=
−
∞


∞


…

∑


m

d
−
1


=
−
∞


∞


F
(

m

0


,
…
,

m

d
−
1


)
H
(

n

0


−

m

0


,
…
,

n

d
−
1


−

m

d
−
1


)


{\displaystyle G(n_{0},\dots ,n_{d-1})=\sum _{m_{0}=-\infty }^{\infty }\dots \sum _{m_{d-1}=-\infty }^{\infty }F(m_{0},\dots ,m_{d-1})H(n_{0}-m_{0},\dots ,n_{d-1}-m_{d-1})}


for all 





d




{\displaystyle {\mathit {d}}}

-tuples of integers 



(

n

0


,
…
,

n

d
−
1


)
∈


Z


d




{\displaystyle (n_{0},\dots ,n_{d-1})\in \mathbb {Z} ^{d}}

. Assume






F




{\displaystyle {\mathit {F}}}

 and 





H




{\displaystyle {\mathit {H}}}

 (and therefore 





G




{\displaystyle {\mathit {G}}}

) are non-zero over only a finite domain bounded by the origin, hence possible to represent as finite multi-dimensional arrays or nested lists 





f




{\displaystyle {\mathit {f}}}

, 





h




{\displaystyle {\mathit {h}}}

, and 





g




{\displaystyle {\mathit {g}}}

.

For this task, implement a function (or method, procedure, subroutine, etc.) deconv to perform deconvolution (i.e., the inverse of convolution) by solving for 





h




{\displaystyle {\mathit {h}}}

 given 





f




{\displaystyle {\mathit {f}}}

 and 





g




{\displaystyle {\mathit {g}}}

. (See Deconvolution/1D for details.)

 The function should work for 





g




{\displaystyle {\mathit {g}}}

 of arbitrary length in each dimension (i.e., not hard coded or constant) and 





f




{\displaystyle {\mathit {f}}}

 of any length up to that of 





g




{\displaystyle {\mathit {g}}}

 in the corresponding dimension.
 The deconv function will need to be parameterized by the dimension 





d




{\displaystyle {\mathit {d}}}

 unless the dimension can be inferred from the data structures representing 





g




{\displaystyle {\mathit {g}}}

 and 





f




{\displaystyle {\mathit {f}}}

.
 There may be more equations than unknowns. If convenient, use a function from a library that finds the best fitting solution to an overdetermined system of linear equations (as in the Multiple regression task). Otherwise, prune the set of equations as needed and solve as in the Reduced row echelon form task.
 Debug your solution using this test data, of which a portion is shown below. Be sure to verify both that the deconvolution of 





g




{\displaystyle {\mathit {g}}}

 with 





f




{\displaystyle {\mathit {f}}}

 is 





h




{\displaystyle {\mathit {h}}}

 and that the deconvolution of 





g




{\displaystyle {\mathit {g}}}

 with 





h




{\displaystyle {\mathit {h}}}

 is 





f




{\displaystyle {\mathit {f}}}

. Display the results in a human readable form for the three dimensional case only.
dimension 1:

h: [-8, 2, -9, -2, 9, -8, -2]
f: [ 6, -9, -7, -5]
g: [-48, 84, -16, 95, 125, -70, 7, 29, 54, 10]

dimension 2:

h: [
      [-8, 1, -7, -2, -9, 4], 
      [4, 5, -5, 2, 7, -1], 
      [-6, -3, -3, -6, 9, 5]]
f: [
      [-5, 2, -2, -6, -7], 
      [9, 7, -6, 5, -7], 
      [1, -1, 9, 2, -7], 
      [5, 9, -9, 2, -5], 
      [-8, 5, -2, 8, 5]]
g: [
      [40, -21, 53, 42, 105, 1, 87, 60, 39, -28], 
      [-92, -64, 19, -167, -71, -47, 128, -109, 40, -21], 
      [58, 85, -93, 37, 101, -14, 5, 37, -76, -56], 
      [-90, -135, 60, -125, 68, 53, 223, 4, -36, -48], 
      [78, 16, 7, -199, 156, -162, 29, 28, -103, -10], 
      [-62, -89, 69, -61, 66, 193, -61, 71, -8, -30], 
      [48, -6, 21, -9, -150, -22, -56, 32, 85, 25]]

dimension 3:

h: [
      [[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]], 
      [[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]]
f: [
      [[-9, 5, -8], [3, 5, 1]], 
      [[-1, -7, 2], [-5, -6, 6]], 
      [[8, 5, 8],[-2, -6, -4]]]
g: [
      [
         [54, 42, 53, -42, 85, -72], 
         [45, -170, 94, -36, 48, 73], 
         [-39, 65, -112, -16, -78, -72], 
         [6, -11, -6, 62, 49, 8]], 
      [
         [-57, 49, -23, 52, -135, 66], 
         [-23, 127, -58, -5, -118, 64], 
         [87, -16, 121, 23, -41, -12], 
         [-19, 29, 35, -148, -11, 45]], 
      [
         [-55, -147, -146, -31, 55, 60], 
         [-88, -45, -28, 46, -26, -144], 
         [-12, -107, -34, 150, 249, 66], 
         [11, -15, -34, 27, -78, -50]], 
      [
         [56, 67, 108, 4, 2, -48], 
         [58, 67, 89, 32, 32, -8], 
         [-42, -31, -103, -30, -23, -8],
         [6, 4, -26, -10, 26, 12]]]

package require Tcl 8.5
namespace path {::tcl::mathfunc ::tcl::mathop}
 
# Utility to extract the number of dimensions of a matrix
proc rank m {
    for {set rank 0} {[llength $m] > 1} {incr rank} {
	set m [lindex $m 0]
    }
    return $rank
}
 
# Utility to get the size of a matrix, as a list of lengths
proc size m {
    set r [rank $m]
    set index {}
    set size {}
    for {set i 0} {$i<$r} {incr i} {
	lappend size [llength [lindex $m $index]]
	lappend index 0
    }
    return $size
}
 
# Utility that iterates over the space of coordinates within a matrix.
#
# Arguments:
#   var   The name of the variable (in the caller's context) to set to each
#         coordinate.
#   size  The size of matrix whose coordinates are to be iterated over.
#   body  The script to evaluate (in the caller's context) for each coordinate,
#         with the variable named by 'var' set to the coordinate for the particular
#         iteration.
proc loopcoords {var size body} {
    upvar 1 $var v
    set count [* {*}$size]
    for {set i 0} {$i < $count} {incr i} {
	set coords {}
	set j $i
	for {set s $size} {[llength $s]} {set s [lrange $s 0 end-1]} {
	    set dimension [lindex $s end]
	    lappend coords [expr {$j % $dimension}]
	    set j [expr {$j / $dimension}]
	}
	set v [lreverse $coords]
	uplevel 1 $body
    }
}
 
# Assembles a row, which is one of the simultaneous equations that needs
# to be solved by reducing the whole set to reduced row echelon form. Note
# that each row describes the equation for a single cell of the 'g' function.
#
# Arguments:
#   g	The "result" matrix of the convolution being undone.
#   h	The known "input" matrix of the convolution being undone.
#   gs	The size descriptor of 'g', passed as argument for efficiency.
#   gc	The coordinate in 'g' that we are generating the equation for.
#   fs	The size descriptor of 'f', passed as argument for efficiency.
#   hs	The size descriptor of 'h' (the unknown "input" matrix), passed
#	as argument for efficiency.
proc row {g f gs gc fs hs} {
    loopcoords hc $hs {
	set fc {}
	set ok 1
	foreach a $gc b $fs c $hc {
	    set d [expr {$a - $c}]
	    if {$d < 0 || $d >= $b} {
		set ok 0
		break
	    }
	    lappend fc $d
	}
	if {$ok} {
	    lappend row [lindex $f $fc]
	} else {
	    lappend row 0
	}
    }
    return [lappend row [lindex $g $gc]]
}
 
# Utility for converting a matrix to Reduced Row Echelon Form
# From http://rosettacode.org/wiki/Reduced_row_echelon_form#Tcl
proc toRREF {m} {
    set lead 0
    set rows [llength $m]
    set cols [llength [lindex $m 0]]
    for {set r 0} {$r < $rows} {incr r} {
	if {$cols <= $lead} {
	    break
	}
	set i $r
	while {[lindex $m $i $lead] == 0} {
	    incr i
	    if {$rows == $i} {
		set i $r
		incr lead
		if {$cols == $lead} {
		    # Tcl can't break out of nested loops
		    return $m
		}
	    }
	}
	# swap rows i and r
	foreach j [list $i $r] row [list [lindex $m $r] [lindex $m $i]] {
	    lset m $j $row
	}
	# divide row r by m(r,lead)
	set val [lindex $m $r $lead]
	for {set j 0} {$j < $cols} {incr j} {
	    lset m $r $j [/ [double [lindex $m $r $j]] $val]
	}
 
	for {set i 0} {$i < $rows} {incr i} {
	    if {$i != $r} {
		# subtract m(i,lead) multiplied by row r from row i
		set val [lindex $m $i $lead]
		for {set j 0} {$j < $cols} {incr j} {
		    lset m $i $j \
			[- [lindex $m $i $j] [* $val [lindex $m $r $j]]]
		}
	    }
	}
	incr lead
    }
    return $m
}
 
# Deconvolve a pair of matrixes. Solves for 'h' such that 'g = f convolve h'.
#
# Arguments:
#   g     The matrix of data to be deconvolved.
#   f     The matrix describing the convolution to be removed.
#   type  Optional description of the type of data expected. Defaults to 32-bit
#         integer data; use 'double' for floating-point data.
proc deconvolve {g f {type int}} {
    # Compute the sizes of the various matrixes involved.
    set gsize [size $g]
    set fsize [size $f]
    foreach gs $gsize fs $fsize {
	lappend hsize [expr {$gs - $fs + 1}]
    }
 
    # Prepare the set of simultaneous equations to solve
    set toSolve {}
    loopcoords coords $gsize {
	lappend toSolve [row $g $f $gsize $coords $fsize $hsize]
    }
 
    # Solve the equations
    set solved [toRREF $toSolve]
 
    # Make a dummy result matrix of the right size
    set h 0
    foreach hs [lreverse $hsize] {set h [lrepeat $hs $h]}
 
    # Fill the results from the equations into the result matrix
    set idx 0
    loopcoords coords $hsize {
	lset h $coords [$type [lindex $solved $idx end]]
	incr idx
    }
 
    return $h
}