mlib3/mlib3math.cpp
//v1 is moving, v2 is target
v1 is moving, v2 is target
v1 is moving v2 is target
bool MMath::MutualInfoNeib(Volume& v1, Volume& v2, Volume& vout, int neib_sz)
{
	if (v1.m_nElements != v2.m_nElements || neib_sz > v1.SX() / 8) return false;
	Matrix hist1, hist2, jhist;
	int nbins = 50;
	Range3 sub;
	Real res = 0, h1 = 0, h2 = 0;
	Real eps = 1e-14;
	int i, j;
	vout = v1; vout = 0;
	for (int z = neib_sz; z < v1.SZ() - neib_sz; z++) {
		cout << "slice " << z << "/" << v1.SZ() - neib_sz << "\r" << std::flush;
		for (int y = neib_sz; y < v1.SY() - neib_sz; y++) {
			for (int x = neib_sz; x < v1.SX() - neib_sz; x++) {
				sub.Set(x - neib_sz, x + neib_sz, y - neib_sz, y + neib_sz, z - neib_sz, z + neib_sz);
				v1.HistSubvol(hist1, nbins, sub); hist1 /= hist1.Sum();
				v2.HistSubvol(hist2, nbins, sub); hist2 /= hist2.Sum();
				v2.JointHistSubvol(jhist, v1, nbins, sub); jhist /= jhist.Sum();
				res = h1 = h2 = 0;
				for (i = 1; i <= nbins; i++)
				{
					h1 -= hist1(i, 1) * log(hist1(i, 1) + eps);
					h2 -= hist2(i, 1) * log(hist2(i, 1) + eps);
				}
				for (i = 1; i <= nbins; i++)
					for (j = 1; j <= nbins; j++)
						res -= jhist(i, j) * log(jhist(i, j) + eps);
				vout(x, y, z) = (h1 + h2) / res;
			}
		}
	}
	return true;
}
bool MMath :: MutualInfoNeib ( Volume & v1 , Volume & v2 , Volume & vout , int neib_sz ) { if ( v1 . m_nElements != v2 . m_nElements || neib_sz > v1 . SX ( ) / 8 ) return false ; Matrix hist1 , hist2 , jhist ; int nbins = 50 ; Range3 sub ; Real res = 0 , h1 = 0 , h2 = 0 ; Real eps = 1e-14 ; int i , j ; vout = v1 ; vout = 0 ; for ( int z = neib_sz ; z < v1 . SZ ( ) - neib_sz ; z ++ ) { cout << " " << z << " " << v1 . SZ ( ) - neib_sz << " \r " << std :: flush ; for ( int y = neib_sz ; y < v1 . SY ( ) - neib_sz ; y ++ ) { for ( int x = neib_sz ; x < v1 . SX ( ) - neib_sz ; x ++ ) { sub . Set ( x - neib_sz , x + neib_sz , y - neib_sz , y + neib_sz , z - neib_sz , z + neib_sz ) ; v1 . HistSubvol ( hist1 , nbins , sub ) ; hist1 /= hist1 . Sum ( ) ; v2 . HistSubvol ( hist2 , nbins , sub ) ; hist2 /= hist2 . Sum ( ) ; v2 . JointHistSubvol ( jhist , v1 , nbins , sub ) ; jhist /= jhist . Sum ( ) ; res = h1 = h2 = 0 ; for ( i = 1 ; i <= nbins ; i ++ ) { h1 -= hist1 ( i , 1 ) * log ( hist1 ( i , 1 ) + eps ) ; h2 -= hist2 ( i , 1 ) * log ( hist2 ( i , 1 ) + eps ) ; } for ( i = 1 ; i <= nbins ; i ++ ) for ( j = 1 ; j <= nbins ; j ++ ) res -= jhist ( i , j ) * log ( jhist ( i , j ) + eps ) ; vout ( x , y , z ) = ( h1 + h2 ) / res ; } } } return true ; }
v1 is moving, v2 is target
v1 is moving v2 is target