Deconvolution/2D+: Difference between revisions

Content added Content deleted

Inline

Revision as of 18:21, 21 March 2010

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 ${\mathit {d}}$ dimensions of two functions

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

taking ${\mathit {d}}$ -tuples of integers to real numbers as the function

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

also taking ${\mathit {d}}$ -tuples of integers to reals and satisfying

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 ${\mathit {d}}$ -tuples of integers $(n_{0},\dots ,n_{d-1})\in \mathbb {Z} ^{d}$ . Assume ${\mathit {F}}$ and ${\mathit {H}}$ (and therefore ${\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 ${\mathit {f}}$ , ${\mathit {h}}$ , and ${\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 ${\mathit {h}}$ given ${\mathit {f}}$ and ${\mathit {g}}$ . (See Deconvolution/1D for details.)

The function should work for ${\mathit {g}}$ of arbitrary length in each dimension (i.e., not hard coded or constant) and ${\mathit {f}}$ of any length up to that of ${\mathit {g}}$ in the corresponding dimension.
The deconv function will need to be parameterized by the dimension ${\mathit {d}}$ unless the dimension can be inferred from the data structures representing ${\mathit {g}}$ and ${\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 ${\mathit {g}}$ with ${\mathit {f}}$ is ${\mathit {h}}$ and that the deconvolution of ${\mathit {g}}$ with ${\mathit {h}}$ is ${\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]]]

J

Actually it is a matter of setting up the linear equations and then solving them.

Solution: <lang j>deconv3 =: 4 : 0

sz  =. x >:@-&$ y                                      NB. shape of z
poi =.  ,<"1 ($y) ,"0/&(,@i.) sz                       NB. pair of indexes
t=. /: sc=: , <@(+"1)/&(#: ,@i.)/ ($y),:sz             NB. order of ,y
T0=. (<"0,x) ,:~ (]/:"1 {.)&.> (<, y) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi   NB. set of boxed equations
T1=. (,x),.~(<0 #~ */sz) (({:@])`({.@])`[})&> {.T0     NB. set of linear equations
sz $ 1e_8 round ({:"1 %. }:"1) T1

)

  h3 -: g3 deconv3 f3                      NB. h3 matches g3 deconv3 f3

1</lang>

Tcl

The trick to doing this (without using a library to do all the legwork for you) is to recast the higher-order solutions into solutions in the 1D case. This is done by regarding an n-dimensional address as a coding of a 1-D address.

<lang tcl>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 f {

   set r [rank $f]
   set index {}
   set size {}
   for {set i 0} {$i<$r} {incr i} {

lappend size [llength [lindex $f $index]] lappend index 0

   }
   return $size

}

Utility that iterates over the space of coordinates within a matrix

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

}

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

}

proc deconvolve {g f {type int}} {

   set gsize [size $g]
   set fsize [size $f]
   foreach gs $gsize fs $fsize {

lappend hsize [expr {$gs - $fs + 1}]

   set toSolve {}
   loopcoords coords $gsize {

lappend toSolve [row $g $f $gsize $coords $fsize $hsize]

   set solved [toRREF $toSolve]

   set h 0
   foreach hs [lreverse $hsize] {set h [lrepeat $hs $h]}
   set idx 0
   loopcoords coords $hsize {

lset h $coords [$type [lindex $solved $idx end]] incr idx

   return $h

}</lang> Demonstrating how to use for the 3-D case: <lang># A pretty-printer proc pretty matrix {

   set size [rank $matrix]
   if {$size == 1} {

return \[[join $matrix ", "]\]

   } elseif {$size == 2} {

set out "" foreach row $matrix { append out " " [pretty $row] ",\n" } return \[[string trimleft [string trimright $out ,\n]]\]

   }
   set rowout {}
   foreach row $matrix {append rowout [pretty $row] ,\n}
   set rowout2 {}
   foreach row [split [string trimright $rowout ,\n] \n] {

append rowout2 " " $row \n

   }
   return \[\n[string trimright $rowout2 \n]\n\]

}

The 3D test data

set f {

   {{-9 5 -8} {3 5 1}} 
   {{-1 -7 2} {-5 -6 6}} 
   {{8 5 8} {-2 -6 -4}}

} set 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}} }

Now do the deconvolution and print it out

puts h:\ [pretty [deconvolve $g $f]]</lang> Output:

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]]

Ursala

This is done mostly with list operations that are either primitive or standard library functions in the language (e.g., zipp, zipt, and pad). The equations are solved by the dgelsd function from the Lapack library. The break function breaks a long list into a sequence of sublists according to a given template, and the band function is taken from the Deconvolution/1D solution. <lang Ursala>#import std

import nat

break = ~&r**+ zipt*+ ~&lh*~+ ~&lzyCPrX|\+ -*^|\~&tK33 :^/~& 0!*t

band = pad0+ ~&rSS+ zipt^*D(~&r,^lrrSPT/~&ltK33tx zipt^/~&r ~&lSNyCK33+ zipp0)^/~&rx ~&B->NlNSPC ~&bt

deconv = # takes a natural number n to the n-dimensional deconvolution function

~&?\math..div! iota; ~&!*; @h|\; (~&al^?\~&ar break@alh2faltPrXPRX)^^(

  ~&B->NlC~&bt*++ gang@t+ ~~*,
  lapack..dgelsd^^(
     (~&||0.!**+ ~&B^?a\~&Y@a ^lriFhNSS2iDrlYSK7LS2SL2rQ/~&alt band@alh2faltPrDPMX)^|\~&+ gang,
     @t =>~&l ~&L+@r))</lang>

The equations tend to become increasingly sparse in higher dimensions, so the following alternative implementation uses the sparse matrix solver from the UMFPACK library instead of Lapack, which is also callable in Ursala, adjusted as shown for the different calling convention. <lang Ursala>deconv = # takes a number n to the n-dimensional deconvolution function

~&?\math..div! iota; ~&!*; @h|\; -+

  //+ ~&al^?\~&ar @alh2faltPrXPRX @liX ~&arr2arl2arrh3falrbt2XPRXlrhPCrtPCPNfallrrPXXPRCQNNCq,
  ^^/-+~&B->NlC~&bt*+,gang@t,~~*+- (umf..di_a_trp^/~&DSLlrnPXrmPXS+num@lmS ^niK10mS/num@r ~&lnS)^^(
     gang; ^|\~&; //+ -+
        ^niK10/~& @NnmlSPASX ~&r->lL @lrmK2K8SmtPK20PPPX ^/~&rrnS2lC ~&rnPrmPljASmF@rrmhPSPlD,
        num+ ~&B^?a\~&Y@a -+
           ~&l?\~&r *=r ~&K7LS+ * (*D ^\~&rr sum@lrlPX)^*D\~&r product^|/~& successor@zhl,
           ^/~&alt @alh2faltPrDPMX -+
              ~&rFS+ num*rSS+ zipt^*D/~&r ^lrrSPT/~&ltK33tx zipt^/~&r ~&lSNyCK33+ zipp0,
              ^/~&rx ~&B->NlNSPC ~&bt+-+-+-,
     @t =>~&l ~&L+@r)+-</lang>

UMFPACK doesn't solve systems with more equations than unknowns, so the system is pruned to a square matrix by first selecting an equation containing only a single variable, then selecting one from those remaining that contains only a single variable not already selected, and so on until all variables are covered, with any remaining unselected equations discarded. A random selection is made whenever there is a choice. This method will cope with larger data sets than feasible using dense and overdetermined matrices, but is less robust in the presence of noise. However, some improvement may be possible by averaging the results over several runs. Here is a test program. <lang Ursala>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.>>>

cast %eLLLm

test =

<

  'h': deconv3(g,f),
  'f': deconv3(g,h)>

</lang> output:

<
   'h': <
      <
         <
            -6.000000e+00,
            -8.000000e+00,
            -5.000000e+00,
            9.000000e+00>,
         <
            -7.000000e+00,
            9.000000e+00,
            -6.000000e+00,
            -8.000000e+00>,
         <
            2.000000e+00,
            -7.000000e+00,
            9.000000e+00,
            8.000000e+00>>,
      <
         <
            7.000000e+00,
            4.000000e+00,
            4.000000e+00,
            -6.000000e+00>,
         <
            9.000000e+00,
            9.000000e+00,
            4.000000e+00,
            -4.000000e+00>,
         <
            -3.000000e+00,
            7.000000e+00,
            -2.000000e+00,
            -3.000000e+00>>>,
   'f': <
      <
         <-9.000000e+00,5.000000e+00,-8.000000e+00>,
         <3.000000e+00,5.000000e+00,1.000000e+00>>,
      <
         <-1.000000e+00,-7.000000e+00,2.000000e+00>,
         <-5.000000e+00,-6.000000e+00,6.000000e+00>>,
      <
         <8.000000e+00,5.000000e+00,8.000000e+00>,
         <-2.000000e+00,-6.000000e+00,-4.000000e+00>>>>