1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
#include "multutils.h"
#include <iostream>
#include <rumba/matrixio.h>
using namespace RUMBA;
// out had better be type double !!!
void multiply ( const ManifoldMatrix & matrix, ManifoldFile* DataFile, ManifoldFile* out )
{
writeManifoldMatrix(matrix,std::cout);
BaseManifold* x;
intPoint ext = DataFile->extent();
ext.z() = 1;
Manifold<double> Tmp (ext);
// ManifoldMatrix mTmp = makeMatrix( Tmp, true );
ManifoldMatrix mTmp = makeMatrix( Tmp, false );
double* it;
for ( int i = 0; i < DataFile->depth(); ++i )
{
x = DataFile->get ( intPoint(0,0,i,0), ext );
it = mTmp.begin();
for ( int j = 0; j < x->size(); ++j )
*it++ = x->getElementDouble(j);
delete x;
ManifoldMatrix result = matrix * mTmp.transpose();
out->put ( &result.M, intPoint(0,0,i,0) );
}
}
double product_i_j ( int i, int j, BaseManifold* r )
{
int skip = r->pixels();
double sum = 0;
int c1 = i * skip;
int c2 = j * skip;
while ( skip-- )
{
sum += r->getElementDouble(c1) * r->getElementDouble(c2);
++c1;
++c2;
}
return sum;
}
ManifoldMatrix r_t_r ( ManifoldFile* r)
{
int sz = r->timepoints();
ManifoldMatrix M = makeMatrix( sz, sz );
fill (M.begin(), M.end(), 0);
BaseManifold* v;
BaseManifold* w;
for ( int k = 0; k < r->depth(); ++k )
{
w = r->get (
intPoint(0,0,k,0),
intPoint(r->width(),r->height(),1,r->timepoints())
);
for ( int i = 0; i < sz; ++i )
for ( int j = 0; j <= i; ++j )
M.element(i,j) += product_i_j (i,j,w);
}
for ( int i = 0; i < r->timepoints(); ++i )
{
for ( int j = i+1; j < r->timepoints(); ++j )
M.element(i,j) = M.element(j,i);
}
return M;
}
|
