C++ MultidimArray类代码示例

void symmetriseMap(MultidimArray<DOUBLE> &img, FileName &fn_sym, bool do_wrap)
{

	if (img.getDim() != 3)
		REPORT_ERROR("symmetriseMap ERROR: symmetriseMap can only be run on 3D maps!");

	img.setXmippOrigin();

	SymList SL;
	SL.read_sym_file(fn_sym);

	Matrix2D<DOUBLE> L(4, 4), R(4, 4); // A matrix from the list
    MultidimArray<DOUBLE> sum, aux;
    sum = img;
    aux.resize(img);

	for (int isym = 0; isym < SL.SymsNo(); isym++)
    {
        SL.get_matrices(isym, L, R);
        applyGeometry(img, aux, R, IS_INV, do_wrap);
        sum += aux;
    }

	// Overwrite the input
	img = sum / (SL.SymsNo() + 1);

}

 /* Constructor */
 PyObject *
 FourierProjector_new(PyTypeObject *type, PyObject *args, PyObject *kwargs)
 {

     FourierProjectorObject *self = (FourierProjectorObject*) type->tp_alloc(type, 0);
     XMIPP_TRY

     if (self != NULL)
     {
         PyObject *image = NULL;
         double padding_factor, max_freq, spline_degree;
         if (PyArg_ParseTuple(args, "O|ddd", &image, &padding_factor, &max_freq, &spline_degree))
         {
        	 Image_Value(image).getDimensions(self->dims);
        	 MultidimArray<double> *pdata;
        	 Image_Value(image).data->getMultidimArrayPointer(pdata);
        	 pdata->setXmippOrigin();
        	 self->fourier_projector = new FourierProjector(*(pdata), padding_factor, max_freq, spline_degree);

         }
     }
     XMIPP_CATCH

     return (PyObject *) self;
 }

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 FourierProjector::produceSideInfo()
{
    // Zero padding
    MultidimArray<double> Vpadded;
    int paddedDim=(int)(paddingFactor*volumeSize);
    // JMRT: TODO: I think it is a very poor design to modify the volume passed
    // in the construct, it will be padded anyway, so new memory should be allocated
    volume->window(Vpadded,FIRST_XMIPP_INDEX(paddedDim),FIRST_XMIPP_INDEX(paddedDim),FIRST_XMIPP_INDEX(paddedDim),
                   LAST_XMIPP_INDEX(paddedDim),LAST_XMIPP_INDEX(paddedDim),LAST_XMIPP_INDEX(paddedDim));
    volume->clear();
    // Make Fourier transform, shift the volume origin to the volume center and center it
    MultidimArray< std::complex<double> > Vfourier;
    transformer3D.completeFourierTransform(Vpadded,Vfourier);
    ShiftFFT(Vfourier, FIRST_XMIPP_INDEX(XSIZE(Vpadded)), FIRST_XMIPP_INDEX(YSIZE(Vpadded)), FIRST_XMIPP_INDEX(ZSIZE(Vpadded)));
    CenterFFT(Vfourier,true);
    Vfourier.setXmippOrigin();

    // Compensate for the Fourier normalization factor
    double K=(double)(XSIZE(Vpadded)*XSIZE(Vpadded)*XSIZE(Vpadded))/(double)(volumeSize*volumeSize);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Vfourier)
    DIRECT_MULTIDIM_ELEM(Vfourier,n)*=K;
    Vpadded.clear();
    // Compute Bspline coefficients
    if (BSplineDeg==3)
    {
        MultidimArray< double > VfourierRealAux, VfourierImagAux;
        Complex2RealImag(Vfourier, VfourierRealAux, VfourierImagAux);
        Vfourier.clear();
        produceSplineCoefficients(BSPLINE3,VfourierRealCoefs,VfourierRealAux);
        produceSplineCoefficients(BSPLINE3,VfourierImagCoefs,VfourierImagAux);
        //VfourierRealAux.clear();
        //VfourierImagAux.clear();
    }
    else
        Complex2RealImag(Vfourier, VfourierRealCoefs, VfourierImagCoefs);

    // Allocate memory for the 2D Fourier transform
    projection().initZeros(volumeSize,volumeSize);
    projection().setXmippOrigin();
    transformer2D.FourierTransform(projection(),projectionFourier,false);

    // Calculate phase shift terms
    phaseShiftImgA.initZeros(projectionFourier);
    phaseShiftImgB.initZeros(projectionFourier);
    double shift=-FIRST_XMIPP_INDEX(volumeSize);
    double xxshift = -2 * PI * shift / volumeSize;
    for (size_t i=0; i<YSIZE(projectionFourier); ++i)
    {
        double phasey=(double)(i) * xxshift;
        for (size_t j=0; j<XSIZE(projectionFourier); ++j)
        {
            // Phase shift to move the origin of the image to the corner
            double dotp = (double)(j) * xxshift + phasey;
            sincos(dotp,&DIRECT_A2D_ELEM(phaseShiftImgB,i,j),&DIRECT_A2D_ELEM(phaseShiftImgA,i,j));
        }
    }
}

void Steerable::generate3DFilter(MultidimArray<double>& h3D,
    std::vector< MultidimArray<double> > &hx1,
    std::vector< MultidimArray<double> > &hy1,
    std::vector< MultidimArray<double> > &hz1)
{
    h3D.initZeros(XSIZE(hz1[0]),XSIZE(hy1[0]),XSIZE(hx1[0]));
    h3D.setXmippOrigin();
    FOR_ALL_ELEMENTS_IN_ARRAY3D(h3D)
        for (int n=0; n<6; n++)
            h3D(k,i,j)+=(hz1[n](k)*hy1[n](i)*hx1[n](j));    
}

// 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";
}

void getSpectrum(MultidimArray<double>& Min,
                 MultidimArray<double>& spectrum,
                 int spectrum_type)
{

	MultidimArray<Complex > Faux;
	int xsize = XSIZE(Min);
	Matrix1D<double> f(3);
	MultidimArray<double> count(xsize);
	FourierTransformer transformer;

	spectrum.initZeros(xsize);
	count.initZeros();
	transformer.FourierTransform(Min, Faux, false);
	FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM(Faux)
	{
		long int idx = ROUND(sqrt(kp * kp + ip * ip + jp * jp));
		if (spectrum_type == AMPLITUDE_SPECTRUM)
		{
			spectrum(idx) += abs(dAkij(Faux, k, i, j));
		}
		else
		{
			spectrum(idx) += norm(dAkij(Faux, k, i, j));
		}
		count(idx) += 1.;
	}

	for (long int i = 0; i < xsize; i++)
		if (count(i) > 0.)
		{
			spectrum(i) /= count(i);
		}

}

    // Computes the average of a number of frames in movies
    void computeAvg(const FileName &movieFile, int begin, int end, MultidimArray<double> &avgImg)
    {
        ImageGeneric movieStack;
        MultidimArray<double> frameImage, shiftedFrame;
        Matrix1D<double> shiftMatrix(2);
        int N=end-begin+1;

        for (size_t i=begin;i<=end;i++)
        {
            movieStack.readMapped(movieFile,i);
            movieStack().getImage(frameImage);
            if (i==begin)
                avgImg.initZeros(YSIZE(frameImage), XSIZE(frameImage));
            if (darkImageCorr)
                frameImage-=darkImage;
            if (gainImageCorr)
                frameImage/=gainImage;
            if (globalShiftCorr)
            {
                XX(shiftMatrix)=XX(shiftVector[i-1]);
                YY(shiftMatrix)=YY(shiftVector[i-1]);
                translate(LINEAR, shiftedFrame, frameImage, shiftMatrix, WRAP);
                avgImg+=shiftedFrame;
            }
            else
                avgImg+=frameImage;
        }
        avgImg/=double(N);
        frameImage.clear();
        movieStack.clear();
    }

void Projector::griddingCorrect(MultidimArray<DOUBLE> &vol_in)
{
	// Correct real-space map by dividing it by the Fourier transform of the interpolator(s)
	vol_in.setXmippOrigin();
	FOR_ALL_ELEMENTS_IN_ARRAY3D(vol_in)
	{
		DOUBLE r = sqrt((DOUBLE)(k*k+i*i+j*j));
		// if r==0: do nothing (i.e. divide by 1)
		if (r > 0.)
		{
			DOUBLE rval = r / (ori_size * padding_factor);
			DOUBLE sinc = sin(PI * rval) / ( PI * rval);
			//DOUBLE ftblob = blob_Fourier_val(rval, blob) / blob_Fourier_val(0., blob);
			// Interpolation (goes with "interpolator") to go from arbitrary to fine grid
			if (interpolator==NEAREST_NEIGHBOUR && r_min_nn == 0)
			{
				// NN interpolation is convolution with a rectangular pulse, which FT is a sinc function
            	A3D_ELEM(vol_in, k, i, j) /= sinc;
			}
			else if (interpolator==TRILINEAR || (interpolator==NEAREST_NEIGHBOUR && r_min_nn > 0) )
			{
				// trilinear interpolation is convolution with a triangular pulse, which FT is a sinc^2 function
				A3D_ELEM(vol_in, k, i, j) /= sinc * sinc;
			}
			else
				REPORT_ERROR("BUG Projector::griddingCorrect: unrecognised interpolator scheme.");
//#define DEBUG_GRIDDING_CORRECT
#ifdef DEBUG_GRIDDING_CORRECT
			if (k==0 && i==0 && j > 0)
				std::cerr << " j= " << j << " sinc= " << sinc << std::endl;
#endif
		}
	}
}

 // Converts an OpenCV matrix to XMIPP MultidimArray
 void opencv2Xmipp(const cv::Mat &opencvMat, MultidimArray<double> &xmippArray)
 {
     int h = opencvMat.rows;
     int w = opencvMat.cols;
     xmippArray.initZeros(h, w);
     FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xmippArray)
     DIRECT_A2D_ELEM(xmippArray,i,j) = opencvMat.at<float>(i,j);
 }

/* Normalizations ---------------------------------------------------------- */
void normalize_OldXmipp(MultidimArray<double> &I)
{
    double mean,std;
    I.computeAvgStdev(mean,std);
    double istd=1.0/std;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I)
    DIRECT_MULTIDIM_ELEM(I,n)=(DIRECT_MULTIDIM_ELEM(I,n)-mean)*istd;
}

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);

}

void normalize_Near_OldXmipp(MultidimArray<double> &I, const MultidimArray<int> &bg_mask)
{
    double avg=0., stddev, min, max;
    double avgbg, stddevbg, minbg, maxbg;
    I.computeStats(avg, stddev, min, max);
    computeStats_within_binary_mask(bg_mask, I, minbg, maxbg, avgbg,
                                    stddevbg);
    I -= avg;
    I /= stddevbg;
}

/* Do inference ------------------------------------------------------------ */
int EnsembleNaiveBayes::doInference(const Matrix1D<double> &newFeatures,
                                    double &cost, MultidimArray<int> &votes,
                                    Matrix1D<double> &classesProbs, Matrix1D<double> &allCosts)
{
    int nmax=ensemble.size();
    MultidimArray<double> minCost, maxCost;
    votes.initZeros(K);
    minCost.initZeros(K);
    minCost.initConstant(1);
    maxCost.initZeros(K);
    maxCost.initConstant(1);
    double bestMinCost=0;
    int bestClass=0;
    MultidimArray<double> newFeaturesn;
    for (int n=0; n<nmax; n++)
    {
        double costn;
        newFeaturesn.initZeros(XSIZE(ensembleFeatures[n]));
        FOR_ALL_ELEMENTS_IN_ARRAY1D(newFeaturesn)
        newFeaturesn(i)=newFeatures(ensembleFeatures[n](i));
        int k=ensemble[n]->doInference(newFeaturesn, costn, classesProbs, allCosts);
        votes(k)++;
        if (minCost(k)>0 || minCost(k)>costn)
            minCost(k)=costn;
        if (maxCost(k)>0 || maxCost(k)<costn)
            maxCost(k)=costn;
        if (minCost(k)<bestMinCost)
        {
            bestMinCost=minCost(k);
            bestClass=k;
        }
    }
    if      (judgeCombination[bestClass]=='m')
        cost=minCost(bestClass);
    else if (judgeCombination[bestClass]=='M')
        cost=maxCost(bestClass);
    else
        cost=minCost(bestClass);
    return bestClass;
}

void FourierProjector::produceSideInfo()
{
    // Zero padding
    MultidimArray<double> Vpadded;
    int paddedDim=(int)(paddingFactor*volumeSize);
    volume->window(Vpadded,FIRST_XMIPP_INDEX(paddedDim),FIRST_XMIPP_INDEX(paddedDim),FIRST_XMIPP_INDEX(paddedDim),
                   LAST_XMIPP_INDEX(paddedDim),LAST_XMIPP_INDEX(paddedDim),LAST_XMIPP_INDEX(paddedDim));
    volume->clear();
    // Make Fourier transform, shift the volume origin to the volume center and center it
    MultidimArray< std::complex<double> > Vfourier;
    transformer3D.completeFourierTransform(Vpadded,Vfourier);
    ShiftFFT(Vfourier, FIRST_XMIPP_INDEX(XSIZE(Vpadded)), FIRST_XMIPP_INDEX(YSIZE(Vpadded)), FIRST_XMIPP_INDEX(ZSIZE(Vpadded)));
    CenterFFT(Vfourier,true);
    Vfourier.setXmippOrigin();

    // Compensate for the Fourier normalization factor
    double K=(double)(XSIZE(Vpadded)*XSIZE(Vpadded)*XSIZE(Vpadded))/(double)(volumeSize*volumeSize);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Vfourier)
    DIRECT_MULTIDIM_ELEM(Vfourier,n)*=K;
    Vpadded.clear();
    // Compute Bspline coefficients
    if (BSplineDeg==3)
    {
        MultidimArray< double > VfourierRealAux, VfourierImagAux;
        Complex2RealImag(Vfourier, VfourierRealAux, VfourierImagAux);
        Vfourier.clear();
        produceSplineCoefficients(BSPLINE3,VfourierRealCoefs,VfourierRealAux);
        produceSplineCoefficients(BSPLINE3,VfourierImagCoefs,VfourierImagAux);
        //VfourierRealAux.clear();
        //VfourierImagAux.clear();
    }
    else
        Complex2RealImag(Vfourier, VfourierRealCoefs, VfourierImagCoefs);

    // Allocate memory for the 2D Fourier transform
    projection().initZeros(volumeSize,volumeSize);
    projection().setXmippOrigin();
    transformer2D.FourierTransform(projection(),projectionFourier,false);
}