Deconvolution/2D+: Difference between revisions
Content added Content deleted
SqrtNegInf (talk | contribs) (→{{header|Perl 6}}: modernize with shaped arrays, more idiomatic code) |
|||
Line 788: | Line 788: | ||
], |
], |
||
]</pre> |
]</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}}== |
=={{header|Tcl}}== |