C++ MultidimArray::resize方法代码示例

TEST_F( FftwTest, directFourierTransformComplex)
{

    MultidimArray< std::complex< double > > FFT1, complxDouble;
    FourierTransformer transformer1;
    typeCast(mulDouble, complxDouble);
    transformer1.FourierTransform(complxDouble, FFT1, false);
    transformer1.inverseFourierTransform();

    transformer1.inverseFourierTransform();
    MultidimArray<std::complex<double> > auxFFT;
    auxFFT.resize(3,3);
    DIRECT_A2D_ELEM(auxFFT,0,0) = std::complex<double>(2.77778,0);
    DIRECT_A2D_ELEM(auxFFT,0,1) = std::complex<double>(-0.0555556,0.096225);
    DIRECT_A2D_ELEM(auxFFT,0,2) = std::complex<double>(-0.0555556,-0.096225);

    DIRECT_A2D_ELEM(auxFFT,1,0) = std::complex<double>(-0.388889,0.673575) ;
    DIRECT_A2D_ELEM(auxFFT,1,1) = std::complex<double>(-0.388889,-0.096225);
    DIRECT_A2D_ELEM(auxFFT,1,2) = std::complex<double>(-0.0555556,-0.288675);

    DIRECT_A2D_ELEM(auxFFT,2,0) = std::complex<double>(-0.388889,-0.673575) ;
    DIRECT_A2D_ELEM(auxFFT,2,1) = std::complex<double>(-0.0555556,0.288675) ;
    DIRECT_A2D_ELEM(auxFFT,2,2) = std::complex<double>(-0.388889,0.096225) ;
    EXPECT_EQ(FFT1,auxFFT);
    transformer1.cleanup();
}

void resizeMap(MultidimArray<double >& img, int newsize)
{

	FourierTransformer transformer;
	MultidimArray<Complex > FT, FT2;
	transformer.FourierTransform(img, FT, false);
	windowFourierTransform(FT, FT2, newsize);
	if (img.getDim() == 2)
	{
		img.resize(newsize, newsize);
	}
	else if (img.getDim() == 3)
	{
		img.resize(newsize, newsize, newsize);
	}
	transformer.inverseFourierTransform(FT2, img);

}

// MORE INFO HERE: http://code.google.com/p/googletest/wiki/AdvancedGuide
// Modify this test so it uses Fixures as test_image and test_metadata
TEST( MultidimTest, Size)
{
    MultidimArray<int> md;
    md.resize(2,3);
    size_t x,y,z, n;
    md.getDimensions(x,y,z,n);
    EXPECT_EQ((size_t)1, n) << "MultidimArray: wrong n size";
    EXPECT_EQ((size_t)1, z) << "MultidimArray: wrong y size";
    EXPECT_EQ((size_t)2, y) << "MultidimArray: wrong z size";
    EXPECT_EQ((size_t)3, x) << "MultidimArray: wrong x size";
}

/** Kullback-Leibner divergence */
double getKullbackLeibnerDivergence(MultidimArray<Complex >& Fimg,
                                    MultidimArray<Complex >& Fref, MultidimArray<double>& sigma2,
                                    MultidimArray<double>& p_i, MultidimArray<double>& q_i, int highshell, int lowshell)
{
	// First check dimensions are OK
	if (!Fimg.sameShape(Fref))
	{
		REPORT_ERROR("getKullbackLeibnerDivergence ERROR: Fimg and Fref are not of the same shape.");
	}

	if (highshell < 0)
	{
		highshell = XSIZE(Fimg) - 1;
	}
	if (lowshell < 0)
	{
		lowshell = 0;
	}

	if (highshell > XSIZE(sigma2))
	{
		REPORT_ERROR("getKullbackLeibnerDivergence ERROR: highshell is larger than size of sigma2 array.");
	}

	if (highshell < lowshell)
	{
		REPORT_ERROR("getKullbackLeibnerDivergence ERROR: highshell is smaller than lowshell.");
	}

	// Initialize the histogram
	MultidimArray<int> histogram;
	int histogram_size = 101;
	int histogram_origin = histogram_size / 2;
	double sigma_max = 10.;
	double histogram_factor = histogram_origin / sigma_max;
	histogram.initZeros(histogram_size);

	// This way this will work in both 2D and 3D
	FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM(Fimg)
	{
		int ires = ROUND(sqrt(kp * kp + ip * ip + jp * jp));
		if (ires >= lowshell && ires <= highshell)
		{
			// Use FT of masked image for noise estimation!
			double diff_real = (DIRECT_A3D_ELEM(Fref, k, i, j)).real - (DIRECT_A3D_ELEM(Fimg, k, i, j)).real;
			double diff_imag = (DIRECT_A3D_ELEM(Fref, k, i, j)).imag - (DIRECT_A3D_ELEM(Fimg, k, i, j)).imag;
			double sigma = sqrt(DIRECT_A1D_ELEM(sigma2, ires));

			// Divide by standard deviation to normalise all the difference
			diff_real /= sigma;
			diff_imag /= sigma;

			// Histogram runs from -10 sigma to +10 sigma
			diff_real += sigma_max;
			diff_imag += sigma_max;

			// Make histogram on-the-fly;
			// Real part
			int ihis = ROUND(diff_real * histogram_factor);
			if (ihis < 0)
			{
				ihis = 0;
			}
			else if (ihis >= histogram_size)
			{
				ihis = histogram_size - 1;
			}
			histogram(ihis)++;
			// Imaginary part
			ihis = ROUND(diff_imag * histogram_factor);
			if (ihis < 0)
			{
				ihis = 0;
			}
			else if (ihis > histogram_size)
			{
				ihis = histogram_size;
			}
			histogram(ihis)++;

		}
	}

	// Normalise the histogram and the discretised analytical Gaussian
	double norm = (double)histogram.sum();
	double gaussnorm = 0.;
	for (int i = 0; i < histogram_size; i++)
	{
		double x = (double)i / histogram_factor;
		gaussnorm += gaussian1D(x - sigma_max, 1. , 0.);
	}

	// Now calculate the actual Kullback-Leibner divergence
	double kl_divergence = 0.;
	p_i.resize(histogram_size);
	q_i.resize(histogram_size);
	for (int i = 0; i < histogram_size; i++)
	{
		// Data distribution
//.........这里部分代码省略.........

// Shift an image through phase-shifts in its Fourier Transform (without pretabulated sine and cosine)
void shiftImageInFourierTransform(MultidimArray<Complex >& in,
                                  MultidimArray<Complex >& out,
                                  double oridim, Matrix1D<double> shift)
{
	out.resize(in);
	shift /= -oridim;
	double dotp, a, b, c, d, ac, bd, ab_cd, x, y, z, xshift, yshift, zshift;
	switch (in.getDim())
	{
	case 1:
		xshift = XX(shift);
		if (ABS(xshift) < XMIPP_EQUAL_ACCURACY)
		{
			out = in;
			return;
		}
		for (long int j = 0; j < XSIZE(in); j++)
		{
			x = j;
			dotp = 2 * PI * (x * xshift);
			a = cos(dotp);
			b = sin(dotp);
			c = DIRECT_A1D_ELEM(in, j).real;
			d = DIRECT_A1D_ELEM(in, j).imag;
			ac = a * c;
			bd = b * d;
			ab_cd = (a + b) * (c + d); // (ab_cd-ac-bd = ad+bc : but needs 4 multiplications)
			DIRECT_A1D_ELEM(out, j) = Complex(ac - bd, ab_cd - ac - bd);
		}
		break;
	case 2:
		xshift = XX(shift);
		yshift = YY(shift);
		if (ABS(xshift) < XMIPP_EQUAL_ACCURACY && ABS(yshift) < XMIPP_EQUAL_ACCURACY)
		{
			out = in;
			return;
		}
		for (long int i = 0; i < XSIZE(in); i++)
			for (long int j = 0; j < XSIZE(in); j++)
			{
				x = j;
				y = i;
				dotp = 2 * PI * (x * xshift + y * yshift);
				a = cos(dotp);
				b = sin(dotp);
				c = DIRECT_A2D_ELEM(in, i, j).real;
				d = DIRECT_A2D_ELEM(in, i, j).imag;
				ac = a * c;
				bd = b * d;
				ab_cd = (a + b) * (c + d);
				DIRECT_A2D_ELEM(out, i, j) =  Complex(ac - bd, ab_cd - ac - bd);
			}
		for (long int i = YSIZE(in) - 1; i >= XSIZE(in); i--)
		{
			y = i - YSIZE(in);
			for (long int j = 0; j < XSIZE(in); j++)
			{
				x = j;
				dotp = 2 * PI * (x * xshift + y * yshift);
				a = cos(dotp);
				b = sin(dotp);
				c = DIRECT_A2D_ELEM(in, i, j).real;
				d = DIRECT_A2D_ELEM(in, i, j).imag;
				ac = a * c;
				bd = b * d;
				ab_cd = (a + b) * (c + d);
				DIRECT_A2D_ELEM(out, i, j) = Complex(ac - bd, ab_cd - ac - bd);
			}
		}
		break;
	case 3:
		xshift = XX(shift);
		yshift = YY(shift);
		zshift = ZZ(shift);
		if (ABS(xshift) < XMIPP_EQUAL_ACCURACY && ABS(yshift) < XMIPP_EQUAL_ACCURACY && ABS(zshift) < XMIPP_EQUAL_ACCURACY)
		{
			out = in;
			return;
		}
		for (long int k = 0; k < ZSIZE(in); k++)
		{
			z = (k < XSIZE(in)) ? k : k - ZSIZE(in);
			for (long int i = 0; i < YSIZE(in); i++)
			{
				y = (i < XSIZE(in)) ? i : i - YSIZE(in);
				for (long int j = 0; j < XSIZE(in); j++)
				{
					x = j;
					dotp = 2 * PI * (x * xshift + y * yshift + z * zshift);
					a = cos(dotp);
					b = sin(dotp);
					c = DIRECT_A3D_ELEM(in, k, i, j).real;
					d = DIRECT_A3D_ELEM(in, k, i, j).imag;
					ac = a * c;
					bd = b * d;
					ab_cd = (a + b) * (c + d);
					DIRECT_A3D_ELEM(out, k, i, j) = Complex(ac - bd, ab_cd - ac - bd);
				}
//.........这里部分代码省略.........