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

/* Filter generation ------------------------------------------------------- */
void Steerable::generate1DFilters(double sigma,
    const MultidimArray<double> &Vtomograph,
    std::vector< MultidimArray<double> > &hx1,
    std::vector< MultidimArray<double> > &hy1,
    std::vector< MultidimArray<double> > &hz1){

    // Initialization 
    MultidimArray<double> aux;
    aux.initZeros(XSIZE(Vtomograph));
    aux.setXmippOrigin();
    for (int i=0; i<6; i++) hx1.push_back(aux);
    
    aux.initZeros(YSIZE(Vtomograph));
    aux.setXmippOrigin();
    for (int i=0; i<6; i++) hy1.push_back(aux);

    aux.initZeros(ZSIZE(Vtomograph));
    aux.setXmippOrigin();
    for (int i=0; i<6; i++) hz1.push_back(aux);

    double sigma2=sigma*sigma;       
    double k1 =  1.0/pow((2.0*PI*sigma),(3.0/2.0));
    double k2 = -1.0/(sigma2);
    
    FOR_ALL_ELEMENTS_IN_ARRAY1D(hx1[0])
    {        
        double i2=i*i;
        double g = -exp(-i2/(2.0*sigma2));
	hx1[0](i) = k1*k2*g*(1.0-(i2/sigma2));
	hx1[1](i) = k1*k2*g;
	hx1[2](i) = k1*k2*g;
	hx1[3](i) = k1*k2*k2*g*i;
	hx1[4](i) = k1*k2*k2*g*i;
	hx1[5](i) = k1*k2*k2*g;
    }    
    FOR_ALL_ELEMENTS_IN_ARRAY1D(hy1[0])
    {
        double i2=i*i;
        double g = -exp(-i2/(2.0*sigma2));
        hy1[0](i) = g;
        hy1[1](i) = g*(1.0-(i2/sigma2));
        hy1[2](i) = g;
        hy1[3](i) = g*i;
        hy1[4](i) = g;
        hy1[5](i) = g*i;
    }
    FOR_ALL_ELEMENTS_IN_ARRAY1D(hz1[0])
    {
        double i2=i*i;
        double g = -exp(-i2/(2.0*sigma2));
	hz1[0](i) = g;
	hz1[1](i) = g;
	hz1[2](i) = g*(1.0-(i2/sigma2));
	hz1[3](i) = g;
	hz1[4](i) = g*i;
	hz1[5](i) = g*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 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);
		}

}

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

/* 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 MlModel::initialise()
{

	// Auxiliary vector with relevant size in Fourier space
	MultidimArray<double > aux;
	aux.initZeros(ori_size / 2 + 1);

	// Now resize all relevant vectors
	Iref.resize(nr_classes);
	pdf_class.resize(nr_classes, 1. / (double)nr_classes);
	pdf_direction.resize(nr_classes);
	group_names.resize(nr_groups, "");
	sigma2_noise.resize(nr_groups, aux);
	nr_particles_group.resize(nr_groups);
	tau2_class.resize(nr_classes, aux);
	fsc_halves_class.resize(nr_classes, aux);
	sigma2_class.resize(nr_classes, aux);
	data_vs_prior_class.resize(nr_classes, aux);
	// TODO handle these two correctly.
	bfactor_correction.resize(nr_groups, 0.);
	scale_correction.resize(nr_groups, 1.);

	acc_rot.resize(nr_classes, 0);
	acc_trans.resize(nr_classes, 0);

	if (ref_dim == 2)
	{
		Matrix1D<double> empty(2);
		prior_offset_class.resize(nr_classes, empty);
	}
	// These arrays will be resized when they are filled
	orientability_contrib.resize(nr_classes);

	Projector ref(ori_size, interpolator, padding_factor, r_min_nn);
	PPref.clear();
	// Now fill the entire vector with instances of "ref"
	PPref.resize(nr_classes, ref);

}

// Fourier ring correlation -----------------------------------------------
// from precalculated Fourier Transforms, and without sampling rate etc.
void getFSC(MultidimArray< Complex >& FT1,
            MultidimArray< Complex >& FT2,
            MultidimArray< double >& fsc)
{
	if (!FT1.sameShape(FT2))
	{
		REPORT_ERROR("fourierShellCorrelation ERROR: MultidimArrays have different shapes!");
	}

	MultidimArray< int > radial_count(XSIZE(FT1));
	MultidimArray<double> num, den1, den2;
	Matrix1D<double> f(3);
	num.initZeros(radial_count);
	den1.initZeros(radial_count);
	den2.initZeros(radial_count);
	fsc.initZeros(radial_count);
	FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM(FT1)
	{
		int idx = ROUND(sqrt(kp * kp + ip * ip + jp * jp));
		if (idx >= XSIZE(FT1))
		{
			continue;
		}
		Complex z1 = DIRECT_A3D_ELEM(FT1, k, i, j);
		Complex z2 = DIRECT_A3D_ELEM(FT2, k, i, j);
		double absz1 = abs(z1);
		double absz2 = abs(z2);
		num(idx) += (conj(z1) * z2).real;
		den1(idx) += absz1 * absz1;
		den2(idx) += absz2 * absz2;
		radial_count(idx)++;
	}

	FOR_ALL_ELEMENTS_IN_ARRAY1D(fsc)
	{
		fsc(i) = num(i) / sqrt(den1(i) * den2(i));
	}

}

    void run()
    {
        // Get angles ==============================================================
        MetaData angles;
        angles.read(fnIn);
        size_t AngleNo = angles.size();
        if (AngleNo == 0 || !angles.containsLabel(MDL_ANGLE_ROT))
            REPORT_ERROR(ERR_MD_BADLABEL, "Input file doesn't contain angular information");

        double maxWeight = -99.e99;
        MultidimArray<double> weight;
        weight.initZeros(AngleNo);
        if (angles.containsLabel(MDL_WEIGHT))
        {
            // Find maximum weight
            int i=0;
            FOR_ALL_OBJECTS_IN_METADATA(angles)
            {
                double w;
                angles.getValue(MDL_WEIGHT,w,__iter.objId);
                DIRECT_A1D_ELEM(weight,i++)=w;
                maxWeight=XMIPP_MAX(w,maxWeight);
            }
        }

void ProgVolumePCA::run()
{
    show();
    produce_side_info();

    const MultidimArray<int> &imask=mask.imask;
    size_t Nvoxels=imask.sum();
    MultidimArray<float> v;
    v.initZeros(Nvoxels);

    // Add all volumes to the analyzer
    FileName fnVol;
    FOR_ALL_OBJECTS_IN_METADATA(mdVols)
    {
        mdVols.getValue(MDL_IMAGE,fnVol,__iter.objId);
        V.read(fnVol);

        // Construct vector
        const MultidimArray<double> &mV=V();
        size_t idx=0;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mV)
        {
            if (DIRECT_MULTIDIM_ELEM(imask,n))
                DIRECT_MULTIDIM_ELEM(v,idx++)=DIRECT_MULTIDIM_ELEM(mV,n);
        }

        analyzer.addVector(v);
    }

    // Construct PCA basis
    analyzer.subtractAvg();
    analyzer.learnPCABasis(NPCA,100);

    // Project onto the PCA basis
    Matrix2D<double> proj;
    analyzer.projectOnPCABasis(proj);
    std::vector<double> dimredProj;
    dimredProj.resize(NPCA);
    int i=0;
    FOR_ALL_OBJECTS_IN_METADATA(mdVols)
    {
        memcpy(&dimredProj[0],&MAT_ELEM(proj,i,0),NPCA*sizeof(double));
        mdVols.setValue(MDL_DIMRED,dimredProj,__iter.objId);
        i++;
    }
    if (fnVolsOut!="")
        mdVols.write(fnVolsOut);
    else
        mdVols.write(fnVols);

    // Save the basis
    const MultidimArray<double> &mV=V();
    for (int i=NPCA-1; i>=0; --i)
    {
        V().initZeros();
        size_t idx=0;
        const MultidimArray<double> &mPCA=analyzer.PCAbasis[i];
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mV)
        {
            if (DIRECT_MULTIDIM_ELEM(imask,n))
                DIRECT_MULTIDIM_ELEM(mV,n)=DIRECT_MULTIDIM_ELEM(mPCA,idx++);
        }
        if (fnBasis!="")
            V.write(fnBasis,i+1,true,WRITE_OVERWRITE);
    }

    // Generate the PCA volumes
    if (listOfPercentiles.size()>0 && fnOutStack!="" && fnAvgVol!="")
    {
        Image<double> Vavg;
        if (fnAvgVol!="")
            Vavg.read(fnAvgVol);
        else
            Vavg().initZeros(V());

        Matrix1D<double> p;
        proj.toVector(p);
        Matrix1D<double> psorted=p.sort();

        Image<double> Vpca;
        Vpca()=Vavg();
        createEmptyFile(fnOutStack,(int)XSIZE(Vavg()),(int)YSIZE(Vavg()),(int)ZSIZE(Vavg()),listOfPercentiles.size());
        std::cout << "listOfPercentiles.size()=" << listOfPercentiles.size() << std::endl;
        for (size_t i=0; i<listOfPercentiles.size(); i++)
        {
            int idx=(int)round(textToFloat(listOfPercentiles[i].c_str())/100.0*VEC_XSIZE(p));
            std::cout << "Percentile " << listOfPercentiles[i] << " -> idx=" << idx << " p(idx)=" << psorted(idx) << std::endl;
            Vpca()+=psorted(idx)*V();
            Vpca.write(fnOutStack,i+1,true,WRITE_REPLACE);
        }
    }
}

EnsembleNaiveBayes::EnsembleNaiveBayes(
    const std::vector < MultidimArray<double> >  &features,
    const Matrix1D<double> &priorProbs,
    int discreteLevels, int numberOfClassifiers,
    double samplingFeatures, double samplingIndividuals,
    const std::string &newJudgeCombination)
{
    int NFeatures=XSIZE(features[0]);
    int NsubFeatures=CEIL(NFeatures*samplingFeatures);
    K=features.size();
    judgeCombination=newJudgeCombination;

#ifdef WEIGHTED_SAMPLING
    // Measure the classification power of each variable
    NaiveBayes *nb_weights=new NaiveBayes(features, priorProbs, discreteLevels);
    MultidimArray<double> weights=nb_weights->__weights;
    delete nb_weights;
    double sumWeights=weights.sum();
#endif

    for (int n=0; n<numberOfClassifiers; n++)
    {
        // Produce the set of features for this subclassifier
        MultidimArray<int> subFeatures(NsubFeatures);
        FOR_ALL_ELEMENTS_IN_ARRAY1D(subFeatures)
        {
#ifdef WEIGHTED_SAMPLING
            double random_sum_weight=rnd_unif(0,sumWeights);
            int j=0;
            do
            {
                double wj=DIRECT_A1D_ELEM(weights,j);
                if (wj<random_sum_weight)
                {
                    random_sum_weight-=wj;
                    j++;
                    if (j==NFeatures)
                    {
                        j=NFeatures-1;
                        break;
                    }
                }
                else
                    break;
            }
            while (true);
            DIRECT_A1D_ELEM(subFeatures,i)=j;
#else

            DIRECT_A1D_ELEM(subFeatures,i)=round(rnd_unif(0,NFeatures-1));
#endif

        }

        // Container for the new training sample
        std::vector< MultidimArray<double> >  newFeatures;

        // Produce the data set for each class
        for (int k=0; k<K; k++)
        {
            int NIndividuals=YSIZE(features[k]);
            int NsubIndividuals=CEIL(NIndividuals*samplingIndividuals);
            MultidimArray<int> subIndividuals(NsubIndividuals);
            FOR_ALL_ELEMENTS_IN_ARRAY1D(subIndividuals)
            subIndividuals(i)=ROUND(rnd_unif(0,NsubIndividuals-1));

            MultidimArray<double> newFeaturesK;
            newFeaturesK.initZeros(NsubIndividuals,NsubFeatures);
            const MultidimArray<double>& features_k=features[k];
            FOR_ALL_ELEMENTS_IN_ARRAY2D(newFeaturesK)
            DIRECT_A2D_ELEM(newFeaturesK,i,j)=DIRECT_A2D_ELEM(features_k,
                                              DIRECT_A1D_ELEM(subIndividuals,i),
                                              DIRECT_A1D_ELEM(subFeatures,j));

            newFeatures.push_back(newFeaturesK);
        }

        // Create a Naive Bayes classifier with this data
        NaiveBayes *nb=new NaiveBayes(newFeatures, priorProbs, discreteLevels);
        ensemble.push_back(nb);
        ensembleFeatures.push_back(subFeatures);
    }
}

//#define DEBUG
void normalize_tomography(MultidimArray<double> &I, double tilt, double &mui,
                          double &sigmai, bool tiltMask, bool tomography0,
                          double mu0, double sigma0)
{
    // Tilt series are always 2D
    I.checkDimension(2);

    const int L=2;

    // Build a mask using the tilt angle
    I.setXmippOrigin();
    MultidimArray<int> mask;
    mask.initZeros(I);
    int Xdimtilt=(int)XMIPP_MIN(FLOOR(0.5*(XSIZE(I)*cos(DEG2RAD(tilt)))),
                           0.5*(XSIZE(I)-(2*L+1)));
    int N=0;
    for (int i=STARTINGY(I); i<=FINISHINGY(I); i++)
        for (int j=-Xdimtilt; j<=Xdimtilt;j++)
        {
            A2D_ELEM(mask,i,j)=1;
            N++;
        }

    // Estimate the local variance
    MultidimArray<double> localVariance;
    localVariance.initZeros(I);
    double meanVariance=0;
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    {
        // Center a mask of size 5x5 and estimate the variance within the mask
        double meanPiece=0, variancePiece=0;
        double Npiece=0;
        for (int ii=i-L; ii<=i+L; ii++)
        {
            if (INSIDEY(I,ii))
                for (int jj=j-L; jj<=j+L; jj++)
                {
                    if (INSIDEX(I,jj))
                    {
                        double pixval=A2D_ELEM(I,ii,jj);
                        meanPiece+=pixval;
                        variancePiece+=pixval*pixval;
                        ++Npiece;
                    }
                }
        }
        meanPiece/=Npiece;
        variancePiece=variancePiece/(Npiece-1)-
                      Npiece/(Npiece-1)*meanPiece*meanPiece;
        A2D_ELEM(localVariance,i,j)=variancePiece;
        meanVariance+=variancePiece;
    }
    meanVariance*=1.0/N;

    // Test the hypothesis that the variance in this piece is
    // the same as the variance in the whole image
    double iFu=1/icdf_FSnedecor(4*L*L+4*L,N-1,0.975);
    double iFl=1/icdf_FSnedecor(4*L*L+4*L,N-1,0.025);
    double iMeanVariance=1.0/meanVariance;
    FOR_ALL_ELEMENTS_IN_ARRAY2D(localVariance)
    {
        if (A2D_ELEM(localVariance,i,j)==0)
            A2D_ELEM(mask,i,j)=-2;
        if (A2D_ELEM(mask,i,j)==1)
        {
            double ratio=A2D_ELEM(localVariance,i,j)*iMeanVariance;
            double thl=ratio*iFu;
            double thu=ratio*iFl;
            if (thl>1 || thu<1)
                A2D_ELEM(mask,i,j)=-1;
        }
    }
#ifdef DEBUG
    Image<double> save;
    save()=I;
    save.write("PPP.xmp");
    save()=localVariance;
    save.write("PPPLocalVariance.xmp");
    Image<int> savemask;
    savemask()=mask;
    savemask.write("PPPmask.xmp");
#endif

    // Compute the statistics again in the reduced mask
    double avg, stddev, min, max;
    computeStats_within_binary_mask(mask, I, min, max, avg, stddev);
    double cosTilt=cos(DEG2RAD(tilt));
    double iCosTilt=1.0/cosTilt;
    if (tomography0)
    {
        double iadjustedStddev=1.0/(sigma0*cosTilt);
        FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
        switch (A2D_ELEM(mask,i,j))
        {
        case -2:
            A2D_ELEM(I,i,j)=0;
            break;
        case 0:
            if (tiltMask)
//.........这里部分代码省略.........

    int main2()
    {

        MultidimArray<double> preImg, avgCurr, mappedImg;
        MultidimArray<double> outputMovie;
        Matrix1D<double> meanStdev;
        ImageGeneric movieStack;
        Image<double> II;
        MetaData MD; // To save plot information
        FileName motionInfFile, flowFileName, flowXFileName, flowYFileName;
        ArrayDim aDim;

        // For measuring times (both for whole process and for each level of the pyramid)
        clock_t tStart, tStart2;

#ifdef GPU
        // Matrix that we required in GPU part
        GpuMat d_flowx, d_flowy, d_dest;
        GpuMat d_avgcurr, d_preimg;
#endif

        // Matrix required by Opencv
        cv::Mat flow, dest, flowx, flowy;
        cv::Mat flowxPre, flowyPre;
        cv::Mat avgcurr, avgstep, preimg, preimg8, avgcurr8;
        cv::Mat planes[]={flowxPre, flowyPre};

        int imagenum, cnt=2, div=0, flowCounter;
        int h, w, levelNum, levelCounter=1;

        motionInfFile=foname.replaceExtension("xmd");
        std::string extension=fname.getExtension();
        if (extension=="mrc")
            fname+=":mrcs";
        movieStack.read(fname,HEADER);
        movieStack.getDimensions(aDim);
        imagenum = aDim.ndim;
        h = aDim.ydim;
        w = aDim.xdim;
        if (darkImageCorr)
        {
            II.read(darkRefFilename);
            darkImage=II();
        }
        if (gainImageCorr)
        {
            II.read(gianRefFilename);
            gainImage=II();
        }
        meanStdev.initZeros(4);
        //avgcurr=cv::Mat::zeros(h, w,CV_32FC1);
        avgCurr.initZeros(h, w);
        flowxPre=cv::Mat::zeros(h, w,CV_32FC1);
        flowyPre=cv::Mat::zeros(h, w,CV_32FC1);
#ifdef GPU

        // Object for optical flow
        FarnebackOpticalFlow d_calc;
        setDevice(gpuDevice);

        // Initialize the parameters for optical flow structure
        d_calc.numLevels=6;
        d_calc.pyrScale=0.5;
        d_calc.fastPyramids=true;
        d_calc.winSize=winSize;
        d_calc.numIters=1;
        d_calc.polyN=5;
        d_calc.polySigma=1.1;
        d_calc.flags=0;
#endif
        // Initialize the stack for the output movie
        if (saveCorrMovie)
            outputMovie.initZeros(imagenum, 1, h, w);
        // Correct for global motion from a cross-correlation based algorithms
        if (globalShiftCorr)
        {
            Matrix1D<double> shiftMatrix(2);
            shiftVector.reserve(imagenum);
            shiftMD.read(globalShiftFilename);
            FOR_ALL_OBJECTS_IN_METADATA(shiftMD)
            {
                shiftMD.getValue(MDL_SHIFT_X, XX(shiftMatrix), __iter.objId);
                shiftMD.getValue(MDL_SHIFT_Y, YY(shiftMatrix), __iter.objId);
                shiftVector.push_back(shiftMatrix);
            }
        }
        tStart2=clock();
        // Compute the average of the whole stack
        fstFrame++; // Just to adapt to Li algorithm
        lstFrame++; // Just to adapt to Li algorithm
        if (lstFrame>=imagenum || lstFrame==1)
            lstFrame=imagenum;
        imagenum=lstFrame-fstFrame+1;
        levelNum=sqrt(double(imagenum));
        computeAvg(fname, fstFrame, lstFrame, avgCurr);
        // if the user want to save the PSD
        if (doAverage)
        {
            II()=avgCurr;
            II.write(foname);
//.........这里部分代码省略.........

void ProgSortByStatistics::processInputPrepare(MetaData &SF)
{
    PCAMahalanobisAnalyzer tempPcaAnalyzer;
    tempPcaAnalyzer.clear();

    Image<double> img;
    MultidimArray<double> img2;
    MultidimArray<int> radial_count;
    MultidimArray<double> radial_avg;
    Matrix1D<int> center(2);
    center.initZeros();

    if (verbose>0)
        std::cout << " Processing training set ..." << std::endl;

    int nr_imgs = SF.size();
    if (verbose>0)
        init_progress_bar(nr_imgs);
    int c = XMIPP_MAX(1, nr_imgs / 60);
    int imgno = 0, imgnoPCA=0;
    MultidimArray<float> v;
    MultidimArray<int> distance;
    int dim;

    bool thereIsEnable=SF.containsLabel(MDL_ENABLED);
    bool first=true;
    FOR_ALL_OBJECTS_IN_METADATA(SF)
    {
        if (thereIsEnable)
        {
            int enabled;
            SF.getValue(MDL_ENABLED,enabled,__iter.objId);
            if (enabled==-1)
                continue;
        }
        img.readApplyGeo(SF,__iter.objId);
        if (targetXdim!=-1 && targetXdim!=XSIZE(img()))
        	selfScaleToSize(LINEAR,img(),targetXdim,targetXdim,1);
        MultidimArray<double> &mI=img();
        mI.setXmippOrigin();
        mI.statisticsAdjust(0,1);

        // Overall statistics
        Histogram1D hist;
        compute_hist(mI,hist,-4,4,31);

        // Radial profile
        img2.resizeNoCopy(mI);
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(img2)
        {
            double val=DIRECT_MULTIDIM_ELEM(mI,n);
            DIRECT_MULTIDIM_ELEM(img2,n)=val*val;
        }
        if (first)
        {
            radialAveragePrecomputeDistance(img2, center, distance, dim);
            first=false;
        }
        fastRadialAverage(img2, distance, dim, radial_avg, radial_count);

        // Build vector
        v.initZeros(XSIZE(hist)+XSIZE(img2)/2);
        int idx=0;
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(hist)
        v(idx++)=(float)DIRECT_A1D_ELEM(hist,i);
        for (size_t i=0; i<XSIZE(img2)/2; i++)
            v(idx++)=(float)DIRECT_A1D_ELEM(radial_avg,i);

        tempPcaAnalyzer.addVector(v);

        if (imgno % c == 0 && verbose>0)
            progress_bar(imgno);
        imgno++;
        imgnoPCA++;
    }
    if (verbose>0)
        progress_bar(nr_imgs);

    MultidimArray<double> vavg,vstddev;
    tempPcaAnalyzer.computeStatistics(vavg,vstddev);
    tempPcaAnalyzer.evaluateZScore(2,20,false);
    pcaAnalyzer.insert(pcaAnalyzer.begin(), tempPcaAnalyzer);
}

void Steerable::singleFilter(const MultidimArray<double>& Vin,
    MultidimArray<double> &hx1, MultidimArray<double> &hy1, MultidimArray<double> &hz1,
    MultidimArray<double> &Vout){

    MultidimArray< std::complex<double> > H, Aux;
    Vout.initZeros(Vin);

    // Filter in X
    #define MINUS_ONE_POWER(n) (((n)%2==0)? 1:-1)
    FourierTransformer transformer;
    transformer.FourierTransform(hx1,H);
    
    FOR_ALL_ELEMENTS_IN_ARRAY1D(H)
          H(i)*= MINUS_ONE_POWER(i);

    FourierTransformer transformer2;
    
    MultidimArray<double> aux(XSIZE(Vin));
        	   
    transformer2.setReal(aux);		   
		   
    for (size_t k=0; k<ZSIZE(Vin); k++)
        for (size_t i=0; i<YSIZE(Vin); i++)
        {
            for (size_t j=0; j<XSIZE(Vin); j++)
                DIRECT_A1D_ELEM(aux,j)=DIRECT_A3D_ELEM(Vin,k,i,j);
			    
	    transformer2.FourierTransform( );	    
	    transformer2.getFourierAlias( Aux );
	    Aux*=H;
	    transformer2.inverseFourierTransform( );
            	    
	    for (size_t j=0; j<XSIZE(Vin); j++)
                DIRECT_A3D_ELEM(Vout,k,i,j)=XSIZE(aux)*DIRECT_A1D_ELEM(aux,j);
        }

    // Filter in Y
    transformer.FourierTransform(hy1,H);
    
    FOR_ALL_ELEMENTS_IN_ARRAY1D(H)
          H(i)*= MINUS_ONE_POWER(i);

    aux.initZeros(YSIZE(Vin));
    transformer2.setReal(aux);		   
    
    for (size_t k=0; k<ZSIZE(Vin); k++)
        for (size_t j=0; j<XSIZE(Vin); j++)
        {
            for (size_t i=0; i<YSIZE(Vin); i++)
                DIRECT_A1D_ELEM(aux,i)=DIRECT_A3D_ELEM(Vout,k,i,j);

	    transformer2.FourierTransform( );	    
	    transformer2.getFourierAlias( Aux );
	    Aux*=H;
	    transformer2.inverseFourierTransform( );
            
	    for (size_t i=0; i<YSIZE(Vin); i++)
                DIRECT_A3D_ELEM(Vout,k,i,j)=XSIZE(aux)*DIRECT_A1D_ELEM(aux,i);
        }

    // Filter in Z

    transformer.FourierTransform(hz1,H);

    FOR_ALL_ELEMENTS_IN_ARRAY1D(H)
          H(i)*= MINUS_ONE_POWER(i);

    aux.initZeros(ZSIZE(Vin));    
    transformer2.setReal(aux);		   

    for (size_t i=0; i<YSIZE(Vin); i++)
        for (size_t j=0; j<XSIZE(Vin); j++)
        {
            for (size_t k=0; k<ZSIZE(Vin); k++)
                DIRECT_A1D_ELEM(aux,k)=DIRECT_A3D_ELEM(Vout,k,i,j);

	    transformer2.FourierTransform( );	    
	    transformer2.getFourierAlias( Aux );
	    Aux*=H;
	    transformer2.inverseFourierTransform( );

            for (size_t k=0; k<ZSIZE(Vin); k++)
                DIRECT_A3D_ELEM(Vout,k,i,j)=XSIZE(aux)*DIRECT_A1D_ELEM(aux,k);
        }
    
    // If Missing wedge
    if (MW!=NULL)
        MW->removeWedge(Vout);
}

/** 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
//.........这里部分代码省略.........