mlib3/mlib3math.cpp
//return correlation coefficient computed in fixed size neighborhood for entire volume
return correlation coefficient computed in fixed size neighborhood for entire volume
return correlation coefficient computed in fixed size neighborhood for entire volume
bool MMath::CorrCoeffNeib(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 Real CCLocal(Volume& v1, Volume& v2, int x0, int y0, int z0, int neib_sz)
		{
			int x, y, z;
			Real n = 8 * neib_sz * neib_sz * neib_sz, s1 = 0, s2 = 0, t1, t2, acc1 = 0, acc2 = 0, acc3 = 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++) {
						s1 += v1(x, y, z); s2 += v2(x, y, z);
					}
				}
			}
			s1 /= n; s2 /= n;
			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++) {
						t1 = v1(x, y, z) - s1; t2 = v2(x, y, z) - s2; acc1 += t1 * t2; acc2 += t1 * t1; acc3 += t2 * t2;
					}
				}
			}
			if (acc2 * acc3 == 0) return 0;
			return fabs(acc1 / sqrt(acc2 * acc3));
		}
	};
	vout = v1; vout = 0; int x, y, z;
	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++) {
				vout(x, y, z) = EX::CCLocal(v1, v2, x, y, z, neib_sz);
			}
		}
	}
	return true;
}
bool MMath :: CorrCoeffNeib ( 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 Real CCLocal ( Volume & v1 , Volume & v2 , int x0 , int y0 , int z0 , int neib_sz ) { int x , y , z ; Real n = 8 * neib_sz * neib_sz * neib_sz , s1 = 0 , s2 = 0 , t1 , t2 , acc1 = 0 , acc2 = 0 , acc3 = 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 ++ ) { s1 += v1 ( x , y , z ) ; s2 += v2 ( x , y , z ) ; } } } s1 /= n ; s2 /= n ; 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 ++ ) { t1 = v1 ( x , y , z ) - s1 ; t2 = v2 ( x , y , z ) - s2 ; acc1 += t1 * t2 ; acc2 += t1 * t1 ; acc3 += t2 * t2 ; } } } if ( acc2 * acc3 == 0 ) return 0 ; return fabs ( acc1 / sqrt ( acc2 * acc3 ) ) ; } } ; vout = v1 ; vout = 0 ; int x , y , z ; 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 ++ ) { vout ( x , y , z ) = EX :: CCLocal ( v1 , v2 , x , y , z , neib_sz ) ; } } } return true ; }
return correlation coefficient computed in fixed size neighborhood for entire volume
return correlation coefficient computed in fixed size neighborhood for entire volume
//this structure is for potential further optimization, to compute partial sums. //pass 1, sums. //pass 2, cc accumulators //phase 3. compute value for the remaining neighborhoods.