C++ VectorXf::mean方法代码示例

float math_util::get_variance(std::vector<float> distribution_)
{
    // Initialize an Eigen vector representing the distribution
    // In the following the formula for variance is implemented:
    // Var(x) = 1/N * SUM((x-mean(x))^2)
    Eigen::VectorXf distribution;
    distribution = Eigen::Map<Eigen::VectorXf>(&distribution_[0],distribution_.size(),1);
    // Retrieve mean value, mean(x)
    float distribution_mean = distribution.mean();
    // Get deviation vector. (x-mean(x))^2
    distribution = (distribution.array() - distribution_mean) * (distribution.array() - distribution_mean);
    // Get variance from deviation vector. Var(x)
    return distribution.mean();

}

void
mutualproximity::set_normfacts(
        musly_trackid trackid,
        Eigen::VectorXf& sim)
{
    double mu = sim.mean();
    Eigen::VectorXd sim_mu = sim.cast<double>().array() - mu;
    double std = (sim_mu.transpose() * sim_mu);
    std /= (static_cast<double>(sim.size()) - 1.0);

    // allocate space
    if (trackid >= (int)norm_facts.size()) {
        norm_facts.resize(trackid+1);
        norm_facts[trackid].mu = mu;
        norm_facts[trackid].std = sqrt(std);
    }
}

//.........这里部分代码省略.........
    const InputImageType * inputy = readery->GetOutput();
    const InputImageType * inputz = readerz->GetOutput();


    //create the output image
    OutputImageType::Pointer image = OutputImageType::New();
    image->SetRegions(inputx->GetLargestPossibleRegion());
    image->SetSpacing(inputx->GetSpacing());
    image->SetOrigin(inputx->GetOrigin());
    image->Allocate();
    image->FillBuffer(0);

    
    InputImageType::SizeType radius;
    radius[0] = neighborhood_radius;
    radius[1] = neighborhood_radius;
    radius[2] = neighborhood_radius;
    
    itk::ConstNeighborhoodIterator<InputImageType> 
            iteratorx(radius, inputx, image->GetRequestedRegion());
    itk::ConstNeighborhoodIterator<InputImageType> 
            iteratory(radius, inputy, image->GetRequestedRegion());
    itk::ConstNeighborhoodIterator<InputImageType> 
            iteratorz(radius, inputz, image->GetRequestedRegion());

    itk::ImageRegionIterator<OutputImageType>  it( image, image->GetRequestedRegion() );
    
    iteratorx.GoToBegin();
    iteratory.GoToBegin();
    iteratorz.GoToBegin();
    it.GoToBegin();
  
    const unsigned int nneighbors = (neighborhood_radius*2+1)*(neighborhood_radius*2+1)*(neighborhood_radius*2+1);
    Eigen::MatrixXf neighbors( nneighbors-1, 3 ); //to store the vectors
    Eigen::Vector3f central_vector;
    
    int c_id = (int) (nneighbors / 2); // get offset of center pixel
    
    
    //check image size for the progress indicator
    OutputImageType::RegionType region = image->GetLargestPossibleRegion();
    OutputImageType::SizeType size = region.GetSize();
    const unsigned long int npixels = size[0]*size[1]*size[2];
    
    unsigned long int k = 0;
    while ( ! iteratorx.IsAtEnd() )
    {
        if(k%10000 == 0)
        {
            std::cout<<"Processed "<<k<<"/"<<npixels<<"\r"<<std::flush;
        }
            
        k++;
            
        unsigned int j = 0;
        for (unsigned int i = 0; i < nneighbors; ++i)
        {    
            if(i!=c_id) //exclude the center
            {
                neighbors(j, 0) = iteratorx.GetPixel(i);  
                neighbors(j, 1) = iteratory.GetPixel(i);  
                neighbors(j, 2) = iteratorz.GetPixel(i);  
                j++;
            }
            else
            {
                central_vector(0) = iteratorx.GetPixel(i);  
                central_vector(1) = iteratory.GetPixel(i);  
                central_vector(2) = iteratorz.GetPixel(i);                  
            }
        }
        
        // calculate the angles 
        Eigen::VectorXf dp = neighbors*central_vector;
        
        const float variance = ((dp.array()-dp.mean()).pow(2)).mean();
        it.Set(variance);        
        
        ++iteratorx;
        ++iteratory;
        ++iteratorz;
        ++it;
    }    
    std::cout<<std::endl;
    
    //  The output of the resampling filter is connected to a writer and the
    //  execution of the pipeline is triggered by a writer update.
    try {
        writer->SetInput(image);
        writer->Update();
    } catch (itk::ExceptionObject & excep) {
        std::cerr << "Exception caught !" << std::endl;
        std::cerr << excep << std::endl;
    }

    

    
    return EXIT_SUCCESS;
}