Deconvolution/2D+: Difference between revisions

Content added Content deleted
(Formatted D entry, not improved)
Line 320: Line 320:


=={{header|D}}==
=={{header|D}}==
<lang d>import std.stdio, std.conv ;
<lang d>import std.stdio, std.conv, std.array, std.algorithm, std.numeric,
import std.array, std.algorithm, std.numeric, std.range ;
std.range;


class M(T) {
class M(T) {
private size_t[] dim ;
private size_t[] dim;
private size_t[] subsize ;
private size_t[] subsize;
private T[] d ;
private T[] d;


this(size_t[] dimension ...) {
this(size_t[] dimension ...) {
setDimension(dimension) ;
setDimension(dimension);
d[] = 0 ; // init each entry to zero ;
d[] = 0; // init each entry to zero;
}
}


M!T dup() {
M!T dup() {
M!T m = new M!T(dim) ;
M!T m = new M!T(dim);
return m.set1DArray(d) ;
return m.set1DArray(d);
}
}

M!T setDimension(size_t[] dimension ...) {
M!T setDimension(size_t[] dimension ...) {
foreach(e ; dimension)
foreach (e; dimension)
assert(e > 0, "no zero dimension") ;
assert(e > 0, "no zero dimension");
dim = dimension.dup ;
dim = dimension.dup;
subsize = dim.dup ;
subsize = dim.dup;
foreach(i;0..dim.length)
foreach (i; 0 .. dim.length)
subsize[i] = reduce!"a*b"(1,dim[i+1..$]) ;
subsize[i] = reduce!"a*b"(1, dim[i + 1 .. $]);
auto dlength = dim[0] * subsize[0] ;
auto dlength = dim[0] * subsize[0];
if(d.length != dlength)
if (d.length != dlength)
d = new T[](dlength) ;
d = new T[dlength];
return this ;
return this;
}
}

M!T set1DArray(T[] t ...) {
M!T set1DArray(T[] t ...) {
auto minLen = min(t.length, d.length) ;
auto minLen = min(t.length, d.length);
d[] = 0 ;
d[] = 0;
d[0..minLen] = t[0..minLen] ;
d[0 .. minLen] = t[0 .. minLen];
return this ;
return this;
}
}


size_t[] seq2idx(size_t seq) {
size_t[] seq2idx(size_t seq) {
size_t acc = seq, tmp ;
size_t acc = seq, tmp;
size_t[] idx ;
size_t[] idx;
foreach(e ; subsize) {
foreach (e; subsize) {
idx ~= tmp = acc / e ;
idx ~= tmp = acc / e;
acc = acc - tmp*e ; // same as % (mod) e ;
acc = acc - tmp * e; // same as % (mod) e;
}
}
return idx ;
return idx;
}
}


size_t size() const @property { return d.length ; }
size_t size() const @property { return d.length; }

⚫
size_t rank() const @property { return dim.length ; }
size_t[] shape() const @property { return dim.dup ; }
size_t rank() const @property { return dim.length; }

⚫
T[] raw() const @property { return d.dup ; }
⚫
size_t[] shape() const @property { return dim.dup; }

⚫
T[] raw() const @property { return d.dup; }


bool checkBound(size_t[] idx ...) const {
bool checkBound(size_t[] idx ...) const {
if(idx.length > dim.length) return false ;
if (idx.length > dim.length)
foreach(i, dm ; idx)
return false;
if(dm >= dim[i]) return false ;
foreach (i, dm; idx)
return true ;
if (dm >= dim[i])
return false;
return true;
}
}

T opIndex(size_t[] idx ...) {
T opIndex(size_t[] idx ...) {
assert(checkBound(idx), "OOPS") ;
assert(checkBound(idx), "OOPS");
return d[dotProduct(idx,subsize)] ;
return d[dotProduct(idx, subsize)];
}
}

T opIndexAssign(T v, size_t[] idx ...) {
T opIndexAssign(T v, size_t[] idx ...) {
assert(checkBound(idx), "OOPS") ;
assert(checkBound(idx), "OOPS");
d[dotProduct(idx,subsize)] = v ;
d[dotProduct(idx, subsize)] = v;
return v ;
return v;
}
}

bool opEquals(Object o) {
override bool opEquals(Object o) {
auto rhs = to!(M!T)(o) ;
return dim == rhs.dim && d == rhs.d ;
auto rhs = to!(M!T)(o);
⚫
return dim == rhs.dim && d == rhs.d;
}
}


int opApply(int delegate(ref size_t[]) dg) {
int opApply(int delegate(ref size_t[]) dg) {
size_t[] yieldIdx ;
size_t[] yieldIdx;
foreach(i;0..d.length) {
foreach (i; 0 .. d.length) {
yieldIdx = seq2idx(i) ;
yieldIdx = seq2idx(i);
if(dg(yieldIdx))
if (dg(yieldIdx))
break ;
break;
}
}
return 0 ;
return 0;
}
}

int opApply(int delegate(ref size_t[], ref T) dg) {
int opApply(int delegate(ref size_t[], ref T) dg) {
size_t idx1d = 0 ;
size_t idx1d = 0;
foreach(size_t[] idx ; this) {
foreach (size_t[] idx; this) {
if(dg(idx, d[idx1d++]))
if (dg(idx, d[idx1d++]))
break ;
break;
}
}
return 0 ;
return 0;
}
}


M!T convolute(M!T rhs) { // _this_ is h, rhs is f, output g
M!T convolute(M!T rhs) { // _this_ is h, rhs is f, output g
auto dm = dim.dup ;
auto dm = dim.dup;
dm[] += (rhs.dim[] - 1) ;
dm[] += rhs.dim[] - 1;
M!T m = new M!T(dm) ; // dm will be reused as m's idx
M!T m = new M!T(dm); // dm will be reused as m's idx
auto bound = m.size ;
auto bound = m.size;
foreach(i;0..d.length) {
foreach (i; 0 .. d.length) {
auto thisIdx = seq2idx(i) ;
auto thisIdx = seq2idx(i);
foreach(j;0..rhs.d.length) {
foreach (j; 0 .. rhs.d.length) {
dm[] = thisIdx[] + rhs.seq2idx(j)[] ;
dm[] = thisIdx[] + rhs.seq2idx(j)[];
auto midx1d = dotProduct(dm, m.subsize) ;
auto midx1d = dotProduct(dm, m.subsize);
if( midx1d < bound)
if ( midx1d < bound)
m.d[midx1d] += d[i]*rhs.d[j] ;
m.d[midx1d] += d[i] * rhs.d[j];
else
else
break ; // bound reach, ok to break
break; // bound reach, ok to break
}
}
}
}
return m ;
return m;
}
}


M!T deconvolute(M!T rhs) { // _this_ is g, rhs is f, output is h
M!T deconvolute(M!T rhs) { // _this_ is g, rhs is f, output is h
auto dm = dim.dup ;
auto dm = dim.dup;
foreach(i, e ; dm)
foreach (i, e; dm)
assert(e + 1 > rhs.dim[i], "deconv : dimensions is zero or negative") ;
assert(e + 1 > rhs.dim[i],
"deconv : dimensions is zero or negative");
⚫
dm[] -= (rhs.dim[] - 1) ;
dm[] -= (rhs.dim[] - 1);
M!T m = new M!T(dm) ; // dm will be reused as rhs' idx
M!T m = new M!T(dm); // dm will be reused as rhs' idx


foreach(i;0..m.size) {
foreach (i; 0 .. m.size) {
auto idx = m.seq2idx(i) ;
auto idx = m.seq2idx(i);
m.d[i] = this[idx] ;
m.d[i] = this[idx];
foreach(j;0..i) {
foreach (j; 0 .. i) {
auto jdx = m.seq2idx(j) ;
auto jdx = m.seq2idx(j);
dm[] = idx[] - jdx[];
dm[] = idx[] - jdx[];
if(rhs.checkBound(dm))
if (rhs.checkBound(dm))
m.d[i] -= m.d[j]*rhs[dm] ;
m.d[i] -= m.d[j] * rhs[dm];
}
}
m.d[i] /= rhs.d[0] ;
m.d[i] /= rhs.d[0];
}
}
return m ;
return m;
}
}


string toString() { return to!string(d) ; }
override string toString() { return to!string(d); }
}
}


auto fold(T)(T[] arr, ref size_t[] d) {
auto fold(T)(T[] arr, ref size_t[] d) {
if(d.length == 0) d ~= arr.length ;
if (d.length == 0)
static if(is(T U : U[])) {// is arr an array of arrays?
d ~= arr.length;
static if (is(T U : U[])) { // is arr an array of arrays?
assert(arr.length > 0, "no empty dimension") ;
d ~= arr[0].length ;
assert(arr.length > 0, "no empty dimension");
foreach(e ; arr)
d ~= arr[0].length;
assert(e.length == arr[0].length, "not rectangular") ;
foreach (e; arr)
return fold(reduce!"a~b"(arr), d) ;
assert(e.length == arr[0].length, "not rectangular");
return fold(reduce!q{a ~ b}(arr), d);
} else {
} else {
assert(arr.length == reduce!"a*b"(d), "not same size") ;
assert(arr.length == reduce!q{a * b}(d), "not same size");
return arr ;
return arr;
}
}
}
}


auto arr2M(T)(T a) {
auto arr2M(T)(T a) {
size_t[] dm ;
size_t[] dm;
auto d = fold(a, dm) ;
auto d = fold(a, dm);
alias ElementType!(typeof(d)) E ;
alias ElementType!(typeof(d)) E;
auto m = new M!E(dm) ;
auto m = new M!E(dm);
m.set1DArray(d) ;
m.set1DArray(d);
return m ;
return m;
}
}


void main() {
void main() {
alias M!int Mi ;
alias M!int Mi;
auto hh = [[[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]],
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]]] ;
[[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]];
auto ff = [[[-9, 5, -8], [3, 5, 1]],[[-1, -7, 2], [-5, -6, 6]],
auto ff = [[[-9, 5, -8], [3, 5, 1]],[[-1, -7, 2], [-5, -6, 6]],
[[8, 5, 8],[-2, -6, -4]]] ;
[[8, 5, 8],[-2, -6, -4]]];
auto h = arr2M(hh) ;
auto h = arr2M(hh);
auto f = arr2M(ff) ;
auto f = arr2M(ff);


auto g = h.convolute(f) ;
auto g = h.convolute(f);


writeln("g == f convolute h ? ", g == f.convolute(h)) ;
writeln("g == f convolute h ? ", g == f.convolute(h));
writeln("h == g deconv f ? ", h == g.deconvolute(f)) ;
writeln("h == g deconv f ? ", h == g.deconvolute(f));
writeln("f == g deconv h ? ", f == g.deconvolute(h)) ;
writeln("f == g deconv h ? ", f == g.deconvolute(h));
writeln(" f = ", f) ;
writeln(" f = ", f);
writeln("g deconv h = ", g.deconvolute(h)) ;
writeln("g deconv h = ", g.deconvolute(h));
}</lang>
}</lang>
''todo(may be not :): pretty print & convert to normal D array''
''todo(may be not :): pretty print & convert to normal D array''
Retrieved from "https://rosettacode.org/wiki/Deconvolution/2D%2B"