mlib3/mlib3math.cpp
//return conjugate gradient computed in fixed size neighborhood, for the entire volume
return conjugate gradient computed in fixed size neighborhood, for the entire volume
return conjugate gradient computed in fixed size neighborhood for the entire volume
bool MMath::EtaNeib(Volume& v1, Volume& v2, Volume& vout, int neib_sz)
{
	if (v1.m_nElements != v1.m_nElements || neib_sz > v1.SX() / 4) return false;
	struct EX {
		static void EtaLocalSum(Volume& u11, Volume& u12, Volume& u22, int x0, int y0, int z0, int neib_sz, Real& s11, Real& s12, Real& s22, int& n)
		{
			int x, y, z; n = 0; s11 = s12 = s22 = 0;
			for (z = z0 - neib_sz; z < z0 + neib_sz; z++) {
				for (y = y0 - neib_sz; y < y0 + neib_sz; y++) {
					for (x = x0 - neib_sz; x < x0 + neib_sz; x++) {
						if (u11(x, y, z) * u22(x, y, z) <= 0) continue;
						s11 += u11(x, y, z); s12 += u12(x, y, z); s22 += u22(x, y, z); n++;
					}
				}
			}
		}
	};
	Volume g1x, g1y, g1z, g2x, g2y, g2z, u11, u22, u12;
	Gradient1D(v1, g1x, 1); Gradient1D(v1, g1y, 2); Gradient1D(v1, g1z, 3);
	Gradient1D(v2, g2x, 1); Gradient1D(v2, g2y, 2); Gradient1D(v2, g2z, 3);
	int x, y, z;
	u11 = v1; u11 = 0; u22 = v1; u22 = 0; u12 = v1; u12 = 0;
	Real a11, a22, a12, q;
	for (z = 0; z < v1.SZ(); z++) {
		for (y = 0; y < v1.SY(); y++) {
			for (x = 0; x < v1.SX(); x++) {
				a11 = g1x(x, y, z) * g1x(x, y, z) + g1y(x, y, z) * g1y(x, y, z) + g1z(x, y, z) * g1z(x, y, z);
				a22 = g2x(x, y, z) * g2x(x, y, z) + g2y(x, y, z) * g2y(x, y, z) + g2z(x, y, z) * g2z(x, y, z);
				q = a11 * a22;
				if (q <= 0) continue;
				a12 = g1x(x, y, z) * g2x(x, y, z) + g1y(x, y, z) * g2y(x, y, z) + g1z(x, y, z) * g2z(x, y, z);
				u11(x, y, z) = a11; u22(x, y, z) = a22; u12(x, y, z) = (a12 * a12) / sqrt(q);
			}
		}
	}
	Real s11 = 0, s12 = 0, s22 = 0;
	int n = 0;
	vout = v1; vout = 0;
	for (z = neib_sz; z < v1.SZ() - neib_sz; z++) {
		cout << "slice " << z << "/" << v1.SZ() - neib_sz << "\r" << std::flush;
		for (y = neib_sz; y < v1.SY() - neib_sz; y++) {
			for (x = neib_sz; x < v1.SX() - neib_sz; x++) {
				EX::EtaLocalSum(u11, u12, u22, x, y, z, neib_sz, s11, s12, s22, n);
				vout(x, y, z) = (n > 0) ? s12 / sqrt(s11 * s22) : 0;
			}
		}
	}
	return true;
}
bool MMath :: EtaNeib ( Volume & v1 , Volume & v2 , Volume & vout , int neib_sz ) { if ( v1 . m_nElements != v1 . m_nElements || neib_sz > v1 . SX ( ) / 4 ) return false ; struct EX { static void EtaLocalSum ( Volume & u11 , Volume & u12 , Volume & u22 , int x0 , int y0 , int z0 , int neib_sz , Real & s11 , Real & s12 , Real & s22 , int & n ) { int x , y , z ; n = 0 ; s11 = s12 = s22 = 0 ; for ( z = z0 - neib_sz ; z < z0 + neib_sz ; z ++ ) { for ( y = y0 - neib_sz ; y < y0 + neib_sz ; y ++ ) { for ( x = x0 - neib_sz ; x < x0 + neib_sz ; x ++ ) { if ( u11 ( x , y , z ) * u22 ( x , y , z ) <= 0 ) continue ; s11 += u11 ( x , y , z ) ; s12 += u12 ( x , y , z ) ; s22 += u22 ( x , y , z ) ; n ++ ; } } } } } ; Volume g1x , g1y , g1z , g2x , g2y , g2z , u11 , u22 , u12 ; Gradient1D ( v1 , g1x , 1 ) ; Gradient1D ( v1 , g1y , 2 ) ; Gradient1D ( v1 , g1z , 3 ) ; Gradient1D ( v2 , g2x , 1 ) ; Gradient1D ( v2 , g2y , 2 ) ; Gradient1D ( v2 , g2z , 3 ) ; int x , y , z ; u11 = v1 ; u11 = 0 ; u22 = v1 ; u22 = 0 ; u12 = v1 ; u12 = 0 ; Real a11 , a22 , a12 , q ; for ( z = 0 ; z < v1 . SZ ( ) ; z ++ ) { for ( y = 0 ; y < v1 . SY ( ) ; y ++ ) { for ( x = 0 ; x < v1 . SX ( ) ; x ++ ) { a11 = g1x ( x , y , z ) * g1x ( x , y , z ) + g1y ( x , y , z ) * g1y ( x , y , z ) + g1z ( x , y , z ) * g1z ( x , y , z ) ; a22 = g2x ( x , y , z ) * g2x ( x , y , z ) + g2y ( x , y , z ) * g2y ( x , y , z ) + g2z ( x , y , z ) * g2z ( x , y , z ) ; q = a11 * a22 ; if ( q <= 0 ) continue ; a12 = g1x ( x , y , z ) * g2x ( x , y , z ) + g1y ( x , y , z ) * g2y ( x , y , z ) + g1z ( x , y , z ) * g2z ( x , y , z ) ; u11 ( x , y , z ) = a11 ; u22 ( x , y , z ) = a22 ; u12 ( x , y , z ) = ( a12 * a12 ) / sqrt ( q ) ; } } } Real s11 = 0 , s12 = 0 , s22 = 0 ; int n = 0 ; vout = v1 ; vout = 0 ; for ( z = neib_sz ; z < v1 . SZ ( ) - neib_sz ; z ++ ) { cout << " " << z << " " << v1 . SZ ( ) - neib_sz << " \r " << std :: flush ; for ( y = neib_sz ; y < v1 . SY ( ) - neib_sz ; y ++ ) { for ( x = neib_sz ; x < v1 . SX ( ) - neib_sz ; x ++ ) { EX :: EtaLocalSum ( u11 , u12 , u22 , x , y , z , neib_sz , s11 , s12 , s22 , n ) ; vout ( x , y , z ) = ( n > 0 ) ? s12 / sqrt ( s11 * s22 ) : 0 ; } } } return true ; }
return conjugate gradient computed in fixed size neighborhood, for the entire volume
return conjugate gradient computed in fixed size neighborhood for the entire volume
//helper functions; potential further optimization to compute partial sums. //phase 1. pre-compute conjugate gradients. /*
		//phase 2. compute value for the first neighborhood.
		for (z = 0; z < neib_sz; z++) {
			for (y = 0; y < neib_sz; y++) {
				for (x = 0; x < neib_sz; x++) {

					if (u11(x,y,z)*u22(x,y,z)<=0) continue;
					s11+=u11(x,y,z); s12+=u12(x,y,z); s22+=u22(x,y,z); n++;
				}
			}
		}
	*/ //phase 3. compute value for the remaining neighborhoods.