Deconvolution/2D+: Difference between revisions

Updated D entry
(Updated D entry)
Line 320:
 
=={{header|D}}==
<lang d>import std.stdio, std.conv, std.arrayalgorithm, std.algorithmnumeric, std.numeric,range;
std.range;
 
class M(T) {
Line 328 ⟶ 327:
private T[] d;
 
this(size_t[] dimension ...) pure nothrow {
setDimension(dimension);
d[] = 0; // init each entry to zero;
Line 334 ⟶ 333:
 
M!T dup() {
M!Tauto m = new M!T(dim);
return m.set1DArray(d);
}
 
M!T setDimension(size_t[] dimension ...) pure nothrow {
foreach (const e; dimension)
assert(e > 0, "no zero dimension");
dim = dimension.dup;
subsize = dim.dup;
foreach (immutable i; 0 .. dim.length)
subsize[i] = reduce!"q{a * b"}(1, dim[i + 1 .. $]);
autoimmutable dlength = dim[0] * subsize[0];
if (d.length != dlength)
d = new T[dlength];
Line 351 ⟶ 350:
}
 
M!T set1DArray(in T[] t ...) pure nothrow @nogc {
auto minLen = min(t.length, d.length);
d[] = 0;
Line 358 ⟶ 357:
}
 
size_t[] seq2idx(in size_t seq) const pure nothrow {
size_t acc = seq, tmp;
size_t[] idx;
foreach (immutable e; subsize) {
idx ~= tmp = acc / e;
acc = acc - tmp * e; // same as % (mod) e;.
}
return idx;
}
 
size_t size() const @propertypure {nothrow return@nogc d.length;@property }{
return d.length;
}
 
size_t rank() const @propertypure {nothrow return@nogc dim.length;@property }{
return dim.length;
}
 
size_t[] shape() const pure nothrow @property { return dim.dup; }
 
T[] raw() const pure nothrow @property { return d.dup; }
 
bool checkBound(size_t[] idx ...) const pure nothrow @nogc {
if (idx.length > dim.length)
return false;
foreach (immutable i, immutable dm; idx)
if (dm >= dim[i])
return false;
Line 385 ⟶ 388:
}
 
T opIndex(size_t[] idx ...) const pure nothrow @nogc {
assert(checkBound(idx), "OOPS");
return d[dotProduct(idx, subsize)];
}
 
T opIndexAssign(T v, size_t[] idx ...) pure nothrow @nogc {
assert(checkBound(idx), "OOPS");
d[dotProduct(idx, subsize)] = v;
Line 396 ⟶ 399:
}
 
override bool opEquals(Object o) const pure {
autoconst rhs = to!(M!T)(o);
return dim == rhs.dim && d == rhs.d;
}
 
int opApply(int delegate(ref size_t[]) dg) const {
size_t[] yieldIdx;
foreach (immutable i; 0 .. d.length) {
yieldIdx = seq2idx(i);
if (dg(yieldIdx))
Line 413 ⟶ 416:
int opApply(int delegate(ref size_t[], ref T) dg) {
size_t idx1d = 0;
foreach (size_t[] idx; this) {
if (dg(idx, d[idx1d++]))
break;
Line 420 ⟶ 423:
}
 
M!T convolute(M!T rhs) { // _this_ is h, rhs is f, output g.
M!T convolute(M!T rhs) const pure nothrow {
auto dm = dim.dup;
dm[] += rhs.dim[] - 1;
M!T m = new M!T(dm); // dm will be reused as m's idx.
auto bound = m.size;
foreach (immutable i; 0 .. d.length) {
auto thisIdx = seq2idx(i);
foreach (immutable j; 0 .. rhs.d.length) {
dm[] = thisIdx[] + rhs.seq2idx(j)[];
autoimmutable midx1d = dotProduct(dm, m.subsize);
if ( midx1d < bound)
m.d[midx1d] += d[i] * rhs.d[j];
else
break; // boundBound reach, okOK to break.
}
}
Line 439 ⟶ 443:
}
 
M!T deconvolute(M!T rhs) { // _this_ is g, rhs is f, output is h.
M!T deconvolute(M!T rhs) const pure nothrow {
auto dm = dim.dup;
foreach (i, e; dm)
Line 445 ⟶ 450:
"deconv : dimensions is zero or negative");
dm[] -= (rhs.dim[] - 1);
M!Tauto m = new M!T(dm); // dm will be reused as rhs' idx.
 
foreach (immutable i; 0 .. m.size) {
auto idx = m.seq2idx(i);
m.d[i] = this[idx];
foreach (immutable j; 0 .. i) {
autoimmutable jdx = m.seq2idx(j);
dm[] = idx[] - jdx[];
if (rhs.checkBound(dm))
m.d[i] -= m.d[j] * rhs[dm];
}
m.d[i] /= rhs.d[0];
}
return m;
}
 
override string toString() const pure { return to!string(d).text; }
}
 
auto fold(T)(T[] arr, ref size_t[] d) pure {
if (d.length == 0)
d ~= arr.length;
 
static if (is(T U : U[])) { // isIs arr an array of arrays?
assert(arr.length > 0, "no empty dimension");
d ~= arr[0].length;
foreach (e; arr)
assert(e.length == arr[0].length, "notNot rectangular");
return fold(arr.reduce!q{a ~ b}(arr), d);
} else {
assert(arr.length == d.reduce!q{a * b}(d), "notNot same size");
return arr;
}
}
 
auto arr2M(T)(T a) pure {
size_t[] dm;
auto d = fold(a, dm);
alias E = ElementType!(typeof(d)) E;
auto m = new M!E(dm);
m.set1DArray(d);
Line 489 ⟶ 495:
 
void main() {
alias Mi = M!int Mi;
auto hh = [[[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]],
[[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]];
Line 497 ⟶ 503:
auto f = arr2M(ff);
 
autoconst g = h.convolute(f);
 
writeln("g == f convolute h ? ", g == f.convolute(h));
Line 506 ⟶ 512:
}</lang>
''todo(may be not :): pretty print & convert to normal D array''
{{out}}
 
Output:
<pre>g == f convolute h ? true
h == g deconv f ? true