Deconvolution/2D+: Difference between revisions
→{{header|Phix}}
SqrtNegInf (talk | contribs) (→{{header|Perl 6}}: modernize with shaped arrays, more idiomatic code) |
|||
Line 788:
],
]</pre>
=={{header|Phix}}==
{{trans|Tcl}}
Quite frankly I'm fairly astonished that it actually works...<br>
(be warned this contains an exciting mix of 0- and 1- based indexes)
<lang Phix>function m_size(sequence m)
--
-- returns the size of a matrix as a list of lengths
--
sequence res = {}
object me = m
while sequence(me) do
res &= length(me)
me = me[1]
end while
return res
end function
function product(sequence s)
--
-- multiply all elements of s together
--
integer res = s[1]
for i=2 to length(s) do
res *= s[i]
end for
return res
end function
function make_coordset(sequence size)
--
-- returns all points in the matrix, in zero-based indexes,
-- eg {{0,0,0}..{3,3,5}} for a 4x4x6 matrix [96 in total]
--
sequence res = {}
integer count = product(size)
for i=0 to count-1 do
sequence coords = {}
integer j = i
for s=length(size) to 1 by -1 do
integer dimension = size[s]
coords &= mod(j,dimension)
j = floor(j/dimension)
end for
coords = reverse(coords)
res = append(res,coords)
end for
return res
end function
function row(sequence g, f, gs, gc, fs, hs)
--
--# 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.
--
sequence row = {},
coords = make_coordset(hs)
for i=1 to length(coords) do
sequence hc = coords[i]
object fn = f
for k=1 to length(gc) do
integer d = gc[k]-hc[k]
if d<0 or d>=fs[k] then
fn = 0
exit
end if
fn = fn[d+1]
end for
row = append(row,fn)
end for
object gn = g
for i=1 to length(gc) do
gn = gn[gc[i]+1]
end for
row = append(row,gn)
return row
end function
function toRREF(sequence m)
--
-- [renamed] copy of Reduced_row_echelon_form.htm#Phix
-- plus one small tweak, as noted below, exit->return,
-- not that said seems to make any actual difference.
--
integer lead = 1,
rows = length(m),
cols = length(m[1])
for r=1 to rows do
if lead>=cols then exit end if
integer i = r
while m[i][lead]=0 do
i += 1
if i=rows then
i = r
lead += 1
-- if lead=cols then exit end if
if lead=cols then return m end if
end if
end while
-- nb m[i] is assigned before m[r], which matters when i=r:
{m[r],m[i]} = {sq_div(m[i],m[i][lead]),m[r]}
for j=1 to rows do
if j!=r then
m[j] = sq_sub(m[j],sq_mul(m[j][lead],m[r]))
end if
end for
lead += 1
end for
return m
end function
function lset(sequence h, sequence idx, object v)
-- helper routine: store v somewhere deep inside h
integer i1 = idx[1]+1
if length(idx)=1 then
h[i1] = v
else
h[i1] = lset(h[i1],idx[2..$],v)
end if
return h
end function
function deconvolve(sequence g, f)
--
--# 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.
--
-- Compute the sizes of the various matrixes involved.
sequence gsize = m_size(g),
fsize = m_size(f),
hsize = sq_add(sq_sub(gsize,fsize),1)
-- Prepare the set of simultaneous equations to solve
sequence toSolve = {},
coords = make_coordset(gsize)
for i=1 to length(coords) do
toSolve = append(toSolve,row(g,f,gsize,coords[i],fsize,hsize))
end for
-- Solve the equations
sequence solved = toRREF(toSolve)
-- Create a result matrix of the right size
object h = 0
for i=length(hsize) to 1 by -1 do
h = repeat(h,hsize[i])
end for
-- Fill the results from the equations into the result matrix
coords = make_coordset(hsize)
for i=1 to length(coords) do
h = lset(h,coords[i],solved[i][$])
end for
return h
end function
constant 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}}}
ppOpt({pp_Nest,2,pp_IntFmt,"%3d"})
pp(deconvolve(g, f))</lang>
{{out}}
<pre>
{{{ -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|Tcl}}==
| |||