Deconvolution/2D+: Difference between revisions
Content added Content deleted
m (→{{header|J}}: whitespace) |
(→Tcl: Added implementation) |
||
Line 95: | Line 95: | ||
1</lang> |
1</lang> |
||
=={{header|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: |
|||
<pre> |
|||
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]] |
|||
</pre> |
|||
=={{header|Ursala}}== |
=={{header|Ursala}}== |
||
This is done mostly with list operations that are either primitive or standard library functions in the language (e.g., <code>zipp</code>, <code>zipt</code>, and <code>pad</code>). The equations are solved by |
This is done mostly with list operations that are either primitive or standard library functions in the language (e.g., <code>zipp</code>, <code>zipt</code>, and <code>pad</code>). The equations are solved by |