C++ VectorMatrix::getShape方法代码示例

VectorVectorConvolution_Simple::VectorVectorConvolution_Simple(const VectorMatrix &lhs, int dim_x, int dim_y, int dim_z)
	: lhs(lhs), dim_x(dim_x), dim_y(dim_y), dim_z(dim_z)
{
	assert(lhs.getShape().getDim(0) == 3);
	exp_x = lhs.getShape().getDim(1);
	exp_y = lhs.getShape().getDim(2);
	exp_z = lhs.getShape().getDim(3);
}

ConstVectorMatrixAccessor::ConstVectorMatrixAccessor(const VectorMatrix &mat) : mat(mat) 
{
	mat.readLock(0);
	data_x = static_cast<CPUArray*>(mat.getArray(0, 0))->ptr();
	data_y = static_cast<CPUArray*>(mat.getArray(0, 1))->ptr();
	data_z = static_cast<CPUArray*>(mat.getArray(0, 2))->ptr();

	// Precalculate strides
	const int rank = mat.getShape().getRank();
	strides[0] = 1;
	strides[1] = strides[0] * (rank > 0 ? mat.getShape().getDim(0) : 1);
	strides[2] = strides[1] * (rank > 1 ? mat.getShape().getDim(1) : 1);
	strides[3] = strides[2] * (rank > 2 ? mat.getShape().getDim(2) : 1);
}

Vector3d findExtremum_cpu(VectorMatrix &M, int z_slice, int component)
{
	if (M.getShape().getRank() != 3) {
		throw std::runtime_error("findExtremum: Fixme: Need matrix of rank 3");
	}

	if (component < 0 || component > 2) {
		throw std::runtime_error("findExtremum: Invalid 'component' value, must be 0, 1 or 2.");
	}

	const int dim_x = M.getShape().getDim(0);
	const int dim_y = M.getShape().getDim(1);

	VectorMatrix::const_accessor M_acc(M);

	// Find cell with maximum absolute value
	double max_val = -1.0;
	int max_x = -1, max_y = -1;
	for (int y=1; y<dim_y-1; ++y)
	for (int x=1; x<dim_x-1; ++x) {
		const int val = std::fabs(M_acc.get(x, y, z_slice)[component]);
		if (val > max_val) {
			max_val = val;
			max_x = x;
			max_y = y;
		}
	}
	assert(max_x > 0);
	assert(max_y > 0);
	
	// Refine maximum by fitting to sub-cell precision
	const double xdir_vals[3] = {
		M_acc.get(max_x-1, max_y+0, z_slice)[component],
		M_acc.get(max_x+0, max_y+0, z_slice)[component],
		M_acc.get(max_x+1, max_y+0, z_slice)[component]
	};

	const double ydir_vals[3] = {
		M_acc.get(max_x+0, max_y-1, z_slice)[component],
		M_acc.get(max_x+0, max_y+0, z_slice)[component],
		M_acc.get(max_x+0, max_y+1, z_slice)[component]
	};

	return Vector3d(
		fit(max_x-1, max_x+0, max_x+1, xdir_vals[0], xdir_vals[1], xdir_vals[2]),
		fit(max_y-1, max_y+0, max_y+1, ydir_vals[0], ydir_vals[1], ydir_vals[2]),
		static_cast<double>(z_slice)
	);
}

VectorMatrix linearInterpolate(const VectorMatrix &src, Shape dest_dim)
{
	Shape src_dim = src.getShape();

	if (src_dim.getRank() != dest_dim.getRank()) {
		throw std::runtime_error("linearInterpolate: Source and destination matrices need to have the same rank.");
	}

	if (src_dim.getRank() != 3) {
		throw std::runtime_error("linearInterpolate: Fixme: Need to have matrix of rank 3");
	}

	VectorMatrix dest(dest_dim);

	VectorMatrix::      accessor dest_acc(dest);
	VectorMatrix::const_accessor  src_acc(src);

	const bool sing_x = (src_dim.getDim(0) == 1);
	const bool sing_y = (src_dim.getDim(1) == 1);
	const bool sing_z = (src_dim.getDim(2) == 1);

	Vector3d scale(1.0, 1.0, 1.0);
	if (!sing_x) scale.x = double(dest_dim.getDim(0)-1) / double(src_dim.getDim(0)-1);
	if (!sing_y) scale.y = double(dest_dim.getDim(1)-1) / double(src_dim.getDim(1)-1);
	if (!sing_z) scale.z = double(dest_dim.getDim(2)-1) / double(src_dim.getDim(2)-1);

	for (int k=0; k<dest_dim.getDim(2); ++k)
	for (int j=0; j<dest_dim.getDim(1); ++j)
	for (int i=0; i<dest_dim.getDim(0); ++i) {
		// (x,y,z): coordinates of point with indices (i,j,k) in dst matrix
		const double x = i / scale.x;
		const double y = j / scale.y;
		const double z = k / scale.z;

		const double u = x - std::floor(x);
		const double v = y - std::floor(y);
		const double w = z - std::floor(z);

		const int I = std::floor(x);
		const int J = std::floor(y);
		const int K = std::floor(z);

		Vector3d tmp(0.0, 0.0, 0.0);
		if (true                         ) tmp = tmp + (1.0-u) * (1.0-v) * (1.0-w) * src_acc.get(I  , J  , K  );
		if (                      !sing_z) tmp = tmp + (1.0-u) * (1.0-v) *      w  * src_acc.get(I  , J  , K+1);
		if (           !sing_y           ) tmp = tmp + (1.0-u) *      v  * (1.0-w) * src_acc.get(I  , J+1, K  );
		if (           !sing_y && !sing_z) tmp = tmp + (1.0-u) *      v  *      w  * src_acc.get(I  , J+1, K+1);
		if (!sing_x                      ) tmp = tmp +      u  * (1.0-v) * (1.0-w) * src_acc.get(I+1, J  , K  );
		if (!sing_x            && !sing_z) tmp = tmp +      u  * (1.0-v) *      w  * src_acc.get(I+1, J  , K+1);
		if (!sing_x && !sing_y           ) tmp = tmp +      u  *      v  * (1.0-w) * src_acc.get(I+1, J+1, K  );
		if (!sing_x && !sing_y && !sing_z) tmp = tmp +      u  *      v  *      w  * src_acc.get(I+1, J+1, K+1);
		dest_acc.set(i, j, k, tmp);
	}

	return dest;
}