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

/* Do inference for class ------------------------------------------------- */
int EnsembleNaiveBayes::doInferenceForClass(int classNumber, const Matrix1D<double> &newFeatures, double &cost,
        Matrix1D<double> &classesProbs, Matrix1D<double> &allCosts)
{
    int nmax=ensemble.size();
    double minCost=1, maxCost=1;
    int votes=0;
    MultidimArray<double> newFeaturesn;
    for (int n=0; n<nmax; n++)
    {
        double costn;
        const MultidimArray<int> &ensembleFeatures_n=ensembleFeatures[n];
        newFeaturesn.resizeNoCopy(XSIZE(ensembleFeatures_n));
        FOR_ALL_ELEMENTS_IN_ARRAY1D(newFeaturesn)
        {
            int idx=A1D_ELEM(ensembleFeatures_n,i);
            A1D_ELEM(newFeaturesn,i)=VEC_ELEM(newFeatures,idx);
        }

        int k=ensemble[n]->doInference(newFeaturesn, costn, classesProbs, allCosts);
        if (k==classNumber)
        {
            votes++;
            if (minCost>0 || minCost>costn)
                minCost=costn;
            if (maxCost>0 || maxCost<costn)
                maxCost=costn;
        }
    }
    if      (judgeCombination[classNumber]=='m')
        cost=minCost;
    else if (judgeCombination[classNumber]=='M')
        cost=maxCost;
    else
        cost=minCost;
    return votes;
}

//majorAxis and minorAxis is the estimated particle size in px
void ProgSortByStatistics::processInprocessInputPrepareSPTH(MetaData &SF, bool trained)
{
    //#define DEBUG
    PCAMahalanobisAnalyzer tempPcaAnalyzer0;
    PCAMahalanobisAnalyzer tempPcaAnalyzer1;
    PCAMahalanobisAnalyzer tempPcaAnalyzer2;
    PCAMahalanobisAnalyzer tempPcaAnalyzer3;
    PCAMahalanobisAnalyzer tempPcaAnalyzer4;

    //Morphology
    tempPcaAnalyzer0.clear();
    //Signal to noise ratio
    tempPcaAnalyzer1.clear();
    tempPcaAnalyzer2.clear();
    tempPcaAnalyzer3.clear();
    //Histogram analysis, to detect black points and saturated parts
    tempPcaAnalyzer4.clear();

    double sign = 1;//;-1;
    int numNorm = 3;
    int numDescriptors0=numNorm;
    int numDescriptors2=4;
    int numDescriptors3=11;
    int numDescriptors4 = 10;

    MultidimArray<float> v0(numDescriptors0);
    MultidimArray<float> v2(numDescriptors2);
    MultidimArray<float> v3(numDescriptors3);
    MultidimArray<float> v4(numDescriptors4);

    if (verbose>0)
    {
        std::cout << " Sorting particle set by new xmipp method..." << 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;

    bool thereIsEnable=SF.containsLabel(MDL_ENABLED);
    bool first=true;

    // We assume that at least there is one particle
    size_t Xdim, Ydim, Zdim, Ndim;
    getImageSize(SF,Xdim,Ydim,Zdim,Ndim);

    //Initialization:
    MultidimArray<double> nI, modI, tempI, tempM, ROI;
    MultidimArray<bool> mask;
    nI.resizeNoCopy(Ydim,Xdim);
    modI.resizeNoCopy(Ydim,Xdim);
    tempI.resizeNoCopy(Ydim,Xdim);
    tempM.resizeNoCopy(Ydim,Xdim);
    mask.resizeNoCopy(Ydim,Xdim);
    mask.initConstant(true);

    MultidimArray<double> autoCorr(2*Ydim,2*Xdim);
    MultidimArray<double> smallAutoCorr;

    Histogram1D hist;
    Matrix2D<double> U,V,temp;
    Matrix1D<double> D;

    MultidimArray<int> radial_count;
    MultidimArray<double> radial_avg;
    Matrix1D<int> center(2);
    MultidimArray<int> distance;
    int dim;
    center.initZeros();

    v0.initZeros(numDescriptors0);
    v2.initZeros(numDescriptors2);
    v3.initZeros(numDescriptors3);
    v4.initZeros(numDescriptors4);

    ROI.resizeNoCopy(Ydim,Xdim);
    ROI.setXmippOrigin();
    FOR_ALL_ELEMENTS_IN_ARRAY2D(ROI)
    {
        double temp = std::sqrt(i*i+j*j);
        if ( temp < (Xdim/2))
            A2D_ELEM(ROI,i,j)= 1;
        else
            A2D_ELEM(ROI,i,j)= 0;
    }

    Image<double> img;
    FourierTransformer transformer(FFTW_BACKWARD);

    FOR_ALL_OBJECTS_IN_METADATA(SF)
    {
        if (thereIsEnable)
        {
            int enabled;
            SF.getValue(MDL_ENABLED,enabled,__iter.objId);
            if ( (enabled==-1)  )
//.........这里部分代码省略.........

void runThread(ThreadArgument &thArg)
{
    int thread_id = thArg.thread_id;
    ProgXrayImport * ptrProg= (ProgXrayImport *)thArg.workClass;

    MetaData localMD;
    Image<double> Iaux;
    FileName fnImgIn, fnImgOut;
    size_t first = 0, last = 0;
    MultidimArray<char> mask;

    while (ptrProg->td->getTasks(first, last))
    {
        for (size_t i=first; i<=last; i++)
        {
            ptrProg->inMD.getValue(MDL_IMAGE, fnImgIn, ptrProg->objIds[i]);


            MDRow rowGeo;
            ptrProg->readGeoInfo(fnImgIn, rowGeo);
//            ptrProg->readAndCrop(fnImgIn, Iaux, ptrProg->cropSizeX, ptrProg->cropSizeY);

            Iaux.read(fnImgIn);
            Iaux().selfWindow(ptrProg->cropSizeYi,ptrProg->cropSizeXi,
            		(int)(YSIZE(Iaux())-ptrProg->cropSizeYe-1),(int)(XSIZE(Iaux())-ptrProg->cropSizeXe-1));

            Iaux().resetOrigin();




            if (XSIZE(ptrProg->IavgDark())!=0)
            {
                Iaux()-=ptrProg->IavgDark();
                forcePositive(Iaux());
            }


            double currentBeam = 1;
            double expTime = 1;
            double slitWidth = 1;

            if ( ptrProg->dSource == ptrProg->MISTRAL )
            {
                size_t idx = fnImgIn.getPrefixNumber();
                currentBeam = dMi(ptrProg->cBeamArray, idx-1);
                expTime = dMi(ptrProg->expTimeArray, idx-1);
                slitWidth = dMi(ptrProg->slitWidthArray, idx-1);
            }
            else
                ptrProg->readCorrectionInfo(fnImgIn, currentBeam, expTime, slitWidth);

            Iaux() *= 1.0/(currentBeam*expTime*slitWidth);
            if (XSIZE(ptrProg->IavgFlat())!=0)
                Iaux()/=ptrProg->IavgFlat();

            // Assign median filter to zero valued pixels to avoid -inf when applying log10
            Iaux().equal(0,mask);
            mask.resizeNoCopy(Iaux());

            if (XSIZE(ptrProg->bpMask()) != 0)
                mask += ptrProg->bpMask();

            boundMedianFilter(Iaux(), mask);

            if (ptrProg->logFix)
            {
                Iaux().selfLog();
                if (ptrProg->selfAttFix)
                    Iaux() *= -1.;
            }

            fnImgOut.compose(i+1, ptrProg->fnOut);


            size_t objId = localMD.addObject();
            localMD.setValue(MDL_IMAGE,fnImgOut,objId);
            localMD.setRow(rowGeo, objId); //
            //            localMD.setValue(MDL_ANGLE_TILT,Iaux.tilt(),objId);
            Iaux.write(fnImgOut);
            if (thread_id==0)
                progress_bar(i);
        }
    }
    //Lock for update the total counter
    ptrProg->mutex.lock();
    ptrProg->outMD.unionAll(localMD);
    ptrProg->mutex.unlock();
}

void PolyZernikes::fit(const Matrix1D<int> & coef, MultidimArray<double> & im, MultidimArray<double> &weight,
                       MultidimArray<bool> & ROI, int verbose)
{
    this->create(coef);

    size_t xdim = XSIZE(im);
    size_t ydim = YSIZE(im);
    //int numZer = (size_t)coef.sum();
    int numZer = (size_t)coef.sum();

    //Actually polOrder corresponds to the polynomial order +1
    int polOrder=(int)ZERNIKE_ORDER(coef.size());

    im.setXmippOrigin();

    Matrix2D<double> polValue(polOrder,polOrder);

    //First argument means number of images
    //Second argument means number of pixels
    WeightedLeastSquaresHelper weightedLeastSquaresHelper;
    Matrix2D<double>& zerMat=weightedLeastSquaresHelper.A;

    zerMat.resizeNoCopy((size_t)ROI.sum(), numZer);
    double iMaxDim2 = 2./std::max(xdim,ydim);

    size_t pixel_idx=0;

    weightedLeastSquaresHelper.b.resizeNoCopy((size_t)ROI.sum());
    weightedLeastSquaresHelper.w.resizeNoCopy(weightedLeastSquaresHelper.b);

    FOR_ALL_ELEMENTS_IN_ARRAY2D(im)
    {
        if ( (A2D_ELEM(ROI,i,j)))
        {
            //For one i we swap the different j
            double y=i*iMaxDim2;
            double x=j*iMaxDim2;

            //polValue = [ 0    y   y2    y3   ...
            //             x   xy  xy2    xy3  ...
            //             x2  x2y x2y2   x2y3 ]
            //dMij(polValue,py,px) py es fila, px es columna

            for (int py = 0; py < polOrder; ++py)
            {
                double ypy=std::pow(y,py);
                for (int px = 0; px < polOrder; ++px)
                    dMij(polValue,px,py) = ypy*std::pow(x,px);
            }

            Matrix2D<int> *fMat;

            //We generate the representation of the Zernike polynomials
            for (int k=0; k < numZer; ++k)
            {
                fMat = &fMatV[k];

                if (fMat == NULL)
                    continue;

                double temp = 0;
                for (size_t px = 0; px < (*fMat).Xdim(); ++px)
                    for (size_t py = 0; py < (*fMat).Ydim(); ++py)
                        temp += dMij(*fMat,py,px)*dMij(polValue,py,px);

                dMij(zerMat,pixel_idx,k) = temp;
            }

            VEC_ELEM(weightedLeastSquaresHelper.b,pixel_idx)=A2D_ELEM(im,i,j);
            VEC_ELEM(weightedLeastSquaresHelper.w,pixel_idx)=std::abs(A2D_ELEM(weight,i,j));
            ++pixel_idx;
        }
    }

    Matrix1D<double> zernikeCoefficients;
    weightedLeastSquares(weightedLeastSquaresHelper, zernikeCoefficients);
    fittedCoeffs = zernikeCoefficients;

    // Pointer to the image to be fitted
    MultidimArray<double> reconstructed;

    reconstructed.resizeNoCopy(im);
    pixel_idx=0;

    FOR_ALL_ELEMENTS_IN_ARRAY2D(im)
    if (A2D_ELEM(ROI,i,j))
    {
        double temp=0;
        for (int k=0; k < numZer; ++k)
            temp+=dMij(zerMat,pixel_idx,k)*VEC_ELEM(fittedCoeffs,k);

        A2D_ELEM(reconstructed,i,j)=temp;

        if ( fabs(A2D_ELEM(reconstructed,i,j)-A2D_ELEM(im,i,j)) > PI)
            A2D_ELEM(ROI,i,j) = false;

        ++pixel_idx;
    }

    pixel_idx=0;
//.........这里部分代码省略.........