Deconvolution/2D+: Difference between revisions

=={{header|C}}==

Very tedious code: unpacks 2D or 3D matrix into a vector with padding, do 1D FFT, then pack result back into matrix.

#include <stdlib.h>

#include <math.h>

printf("\n");

}

-6 -8 -5 9

-7 9 -6 -8

 

8 5 8

 

=={{header|D}}==

 

class M(T) {

writeln(" f = ", f);

writeln("g deconv h = ", g.deconvolute(h));

''todo(may be not :): pretty print & convert to normal D array''

{{out}}

=={{header|Go}}==

{{trans|C}}

 

import (

}

}

 

{{out}}

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

 

sz =. x >:@-&$ y NB. shape of z

poi =. ,<"1 ($y) ,"0/&(,@i.) sz NB. pair of indexes

sz $ 1e_8 round ({:"1 %. }:"1) T1

)

 

f1=: 6 _9 _7 _5

g1=: _48 84 _16 95 125 _70 7 29 54 10

_42 _31 _103 _30 _23 _8

6 4 _26 _10 26 12

 

1

h2 -: g2 deconv3 f2

1

h3 -: g3 deconv3 f3 NB. -: checks for matching structure and data

 

=={{header|Julia}}==

Julia has a deconv() function that works on Julia's builtin multidimensional arrays, but not on the nested type 2D and 3D arrays used in the task. So, the solution function, deconvn(), sets up repackaging for 1D fft. The actual solving work is done on one line of ifft/fft, and the rest of the code is merely to repackage the nested arrays.

 

const h1 = [-8, 2, -9, -2, 9, -8, -2]

println(deconvn(f3nested, g3nested,

(length(h3nested), length(h3nested[1]), length(h3nested[1][1])))) # 3D

<pre>

Array{Array{Int64,1},1}[[[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]], [[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]]

 

=={{header|Mathematica}} / {{header|Wolfram Language}}==

 

Round[ListDeconvolve[{{-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}},

{-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},

 

{{out}}

=={{header|Nim}}==

{{trans|D}}

 

type Size = uint64

echo "f == g deconv h ? ", f == g.deconvolute(h)

echo " f = ", f

 

{{out}}

{{libheader|ntheory}}

{{trans|Raku}}

use ntheory qw/forsetproduct/;

 

pretty_print(0,@h);

print "\nff =\n";

{{out}}

<pre>3D arrays:

Quite frankly I'm fairly astonished that it actually works...<br>

(be warned this contains an exciting mix of 0- and 1- based indexes)

{{out}}

<pre>

https://math.stackexchange.com/questions/380720/is-deconvolution-simply-division-in-frequency-domain

 

"""

 

 

pprint.pprint(deconv(g,f))

 

Output:

=={{header|Raku}}==

(formerly Perl 6)

 

Translation of Tcl.

sub deconvolve-N ( @g, @f ) {

my @hsize = @g.shape »-« @f.shape »+» 1;

@h.AT-POS(|$_) = $v;

}

}

 

for ^@hc -> $i {

my $window = @gcoord[$i] - @hc[$i];

last;

}

}

@row.push: @g.AT-POS(|@gcoord);

}

 

# Constructs an AoA of coordinates to all elements of N dimensional array

sub coords ( @dim ) {

}

 

# Can handle over-specified systems (N unknowns in N + M equations)

sub rref (@m) {

for ^$rows -> $r {

my $i = $r;

until @m[$i;$lead] {

$i = $r;

}

@m[$i, $r] = @m[$r, $i] if $r != $i;

for ^$rows -> $n {

next if $n == $r;

}

++$lead;

}

 

}

}

}

 

# Pretty printer for N dimensional arrays

# Assumes if first element in level is an array, then all are

if @array[0] ~~ Array {

say ' ' x $indent,"[";

say ' ' x $indent, "]{$indent??','!!''}";

} else {

my @h = deconvolve-N( @g, @f );

say "h =";

my @h-shaped[2;3;4] = @(deconvolve-N( @g, @f ));

my @ff = deconvolve-N( @g, @h-shaped );

say "\nff =";

 

Output:

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.

 

namespace path {::tcl::mathfunc ::tcl::mathop}

 

 

return $h

Demonstrating how to use for the 3-D case:

proc pretty matrix {

set size [rank $matrix]

 

# Now do the deconvolution and print it out

Output:

<pre>

the [http://www.netlib.org/lapack/lug/node27.html <code>dgelsd</code>] function from the Lapack library.

The <code>break</code> function breaks a long list into a sequence of sublists according to a given template, and the <code>band</code> function is taken from the [[Deconvolution/1D]] solution.

#import nat

 

lapack..dgelsd^^(

(~&||0.!**+ ~&B^?a\~&Y@a ^lriFhNSS2iDrlYSK7LS2SL2rQ/~&alt band@alh2faltPrDPMX)^|\~&+ gang,

The equations tend to become increasingly sparse in higher dimensions,

so the following alternative implementation uses the sparse matrix

instead of Lapack, which is also callable in Ursala, adjusted as shown for the different

[http://www.basis.uklinux.net/avram/refman/umf-input-parameters.html calling convention].

 

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

~&rFS+ num*rSS+ zipt^*D/~&r ^lrrSPT/~&ltK33tx zipt^/~&r ~&lSNyCK33+ zipp0,

^/~&rx ~&B->NlNSPC ~&bt+-+-+-,

UMFPACK doesn't solve systems with more equations than unknowns, so

the system is pruned to a square matrix by first selecting an equation

However, some improvement may be possible by averaging the results over several runs.

Here is a test program.

f = <<<-9.,5.,-8.>,<3.,5.,1.>>,<<-1.,-7.,2.>,<-5.,-6.,6.>>,<<8.,5.,8.>,<-2.,-6.,-4.>>>

 

'h': deconv3(g,f),

'f': deconv3(g,h)>

output:

<pre style="height:25ex;overflow:scroll">

{{libheader|Wren-complex}}

{{libheader|Wren-fmt}}

 

var fft2 // recursive

}

if (i < kx-1) System.print()

 

{{out}}