text:
// Copyright (C) 2020 T. Zachary Laine
//
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_TEXT_CONFIG_HPP
#define BOOST_TEXT_CONFIG_HPP

#include <boost/config.hpp>

// Included for definition of __cpp_lib_concepts.
#include <iterator>

/** When you insert into a rope, the incoming sequence may be inserted as a
    new segment, or if it falls within an existing string-segment, it may be
    inserted into the string object used to represent that segment.  This only
    happens if the incoming sequence will fit within the existing segment's
    capacity, or if the segment is smaller than a certain limit.
    BOOST_TEXT_STRING_INSERT_MAX is that limit. */
#ifndef BOOST_TEXT_STRING_INSERT_MAX
#    define BOOST_TEXT_STRING_INSERT_MAX 4096
#endif

#if defined(__cpp_lib_concepts) && !defined(BOOST_TEXT_DISABLE_CONCEPTS)
#    define BOOST_TEXT_USE_CONCEPTS 1
#else
#    define BOOST_TEXT_USE_CONCEPTS 0
#endif

#if defined(__cpp_coroutines) && !defined(BOOST_TEXT_DISABLE_COROUTINES)
#    define BOOST_TEXT_USE_COROUTINES 1
#else
#    define BOOST_TEXT_USE_COROUTINES 0
#endif

#ifndef BOOST_TEXT_DOXYGEN

// The inline namespaces v1 and v2 represent pre- and post-C++20.  v1 is
// inline for standards before C++20, and v2 is inline for C++20 and later.
// Note that this only applies to code for which a v2 namespace alternative
// exists.  Some instances of the v1 namespace may still be inline, if there
// is no v2 version of its contents.
#if BOOST_TEXT_USE_CONCEPTS
#    define BOOST_TEXT_NAMESPACE_V1 namespace v1
#    define BOOST_TEXT_NAMESPACE_V2 inline namespace v2
#else
#    define BOOST_TEXT_NAMESPACE_V1 inline namespace v1
#    define BOOST_TEXT_NAMESPACE_V2 namespace v2
#endif

// Nothing before GCC 6 has proper C++14 constexpr support.
#if defined(__GNUC__) && __GNUC__ < 6 && !defined(__clang__)
#    define BOOST_TEXT_CXX14_CONSTEXPR
#    define BOOST_TEXT_NO_CXX14_CONSTEXPR
#elif defined(_MSC_VER) && _MSC_VER <= 1916
#    define BOOST_TEXT_CXX14_CONSTEXPR
#    define BOOST_TEXT_NO_CXX14_CONSTEXPR
#else
#    define BOOST_TEXT_CXX14_CONSTEXPR BOOST_CXX14_CONSTEXPR
#    if defined(BOOST_NO_CXX14_CONSTEXPR)
#        define BOOST_TEXT_NO_CXX14_CONSTEXPR
#    endif
#endif

// Implements separate compilation features as described in
// http://www.boost.org/more/separate_compilation.html

//  normalize macros

#if !defined(BOOST_TEXT_DYN_LINK) && !defined(BOOST_TEXT_STATIC_LINK) &&       \
    !defined(BOOST_ALL_DYN_LINK) && !defined(BOOST_ALL_STATIC_LINK)
#    define BOOST_TEXT_STATIC_LINK
#endif

#if defined(BOOST_ALL_DYN_LINK) && !defined(BOOST_TEXT_DYN_LINK)
#    define BOOST_TEXT_DYN_LINK
#elif defined(BOOST_ALL_STATIC_LINK) && !defined(BOOST_TEXT_STATIC_LINK)
#    define BOOST_TEXT_STATIC_LINK
#endif

#if defined(BOOST_TEXT_DYN_LINK) && defined(BOOST_TEXT_STATIC_LINK)
#    error Must not define both BOOST_TEXT_DYN_LINK and BOOST_TEXT_STATIC_LINK
#endif

//  enable dynamic or static linking as requested

#if defined(BOOST_ALL_DYN_LINK) || defined(BOOST_TEXT_DYN_LINK)
#    if defined(BOOST_TEXT_SOURCE)
#        define BOOST_TEXT_DECL BOOST_SYMBOL_EXPORT
#    else
#        define BOOST_TEXT_DECL BOOST_SYMBOL_IMPORT
#    endif
#else
#    define BOOST_TEXT_DECL
#endif

#if 0 // TODO: Disabled for now.
//  enable automatic library variant selection

#if !defined(BOOST_TEXT_SOURCE) && !defined(BOOST_ALL_NO_LIB) &&               \
    !defined(BOOST_TEXT_NO_LIB)
//
// Set the name of our library, this will get undef'ed by auto_link.hpp
// once it's done with it:
//
#define BOOST_LIB_NAME boost_text
//
// If we're importing code from a dll, then tell auto_link.hpp about it:
//
#if defined(BOOST_ALL_DYN_LINK) || defined(BOOST_TEXT_DYN_LINK)
#    define BOOST_DYN_LINK
#endif
//
// And include the header that does the work:
//
#include <boost/config/auto_link.hpp>
#endif  // auto-linking disabled
#endif

#endif // doxygen

#endif

max_forks_repo_name: tzlaine/yaml


text:
#pragma once

#include <El.hpp>

#include <boost/filesystem.hpp>

void read_text_block(El::DistMatrix<El::BigFloat> &block,
                     const boost::filesystem::path &block_path);
void read_text_block(El::DistMatrix<El::BigFloat> &block,
                     const boost::filesystem::path &checkpoint_directory,
                     const std::string &prefix, const size_t &block_index);

ext: hxx

max_stars_repo_name: maneandrea/sdpb

max_issues_repo_name: maneandrea/sdpb

max_forks_repo_name: maneandrea/sdpb


text:
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/calib3d/calib3d.hpp>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <Eigen/SVD>
#include <g2o/core/auto_differentiation.h>
#include <g2o/core/base_vertex.h>
#include <g2o/core/base_unary_edge.h>
#include <g2o/core/block_solver.h>
#include <g2o/core/optimization_algorithm_gauss_newton.h>
#include <g2o/core/optimization_algorithm_levenberg.h>
#include <g2o/solvers/dense/linear_solver_dense.h>
#include <chrono>
#include <sophus/se3.hpp>

using namespace std;
using namespace cv;

void find_feature_matches(const Mat &img_1, const Mat &img_2,
                          std::vector<KeyPoint> &keypoints_1,
                          std::vector<KeyPoint> &keypoints_2,
                          std::vector<DMatch> &matches) {
  Mat descriptors_1, descriptors_2;
  // used in OpenCV3
  Ptr<FeatureDetector> detector = ORB::create();
  Ptr<DescriptorExtractor> descriptor = ORB::create();
  // use this if you are in OpenCV2
  // Ptr<FeatureDetector> detector = FeatureDetector::create ( "ORB" );
  // Ptr<DescriptorExtractor> descriptor = DescriptorExtractor::create ( "ORB" );
  Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce-Hamming");
  detector->detect(img_1, keypoints_1);
  detector->detect(img_2, keypoints_2);

  descriptor->compute(img_1, keypoints_1, descriptors_1);
  descriptor->compute(img_2, keypoints_2, descriptors_2);

  vector<DMatch> match;
  // BFMatcher matcher ( NORM_HAMMING );
  matcher->match(descriptors_1, descriptors_2, match);

  double min_dist = 10000, max_dist = 0;

  for (int i = 0; i < descriptors_1.rows; i++) {
    double dist = match[i].distance;
    if (dist < min_dist) min_dist = dist;
    if (dist > max_dist) max_dist = dist;
  }

  printf("-- Max dist : %f \n", max_dist);
  printf("-- Min dist : %f \n", min_dist);

  for (int i = 0; i < descriptors_1.rows; i++) {
    if (match[i].distance <= max(2 * min_dist, 30.0)) {
      matches.push_back(match[i]);
    }
  }
}

Point2d pixel2cam(const Point2d &p, const Mat &K) {
  return Point2d(
    (p.x - K.at<double>(0, 2)) / K.at<double>(0, 0),
    (p.y - K.at<double>(1, 2)) / K.at<double>(1, 1)
  );
}

// Solve ICP with linear algebra (SVD solution)
void pose_estimation_3d3d(vector<Eigen::Vector3d> pts1,
                          vector<Eigen::Vector3d> pts2,
                          Eigen::Matrix3d &R, Eigen::Vector3d &t) {
  Eigen::Vector3d c1(0.0, 0.0, 0.0), c2(0.0, 0.0, 0.0);
  for (auto& p : pts1)
    c1 += p;
  c1 /= pts1.size();
  for (auto& p : pts2)
    c2 += p;
  c2 /= pts2.size();

  for (auto& p : pts1)
    p -= c1;
  for (auto& p : pts2)
    p -= c2;

  Eigen::Matrix3d W = Eigen::Matrix3d::Zero();
  for (int i = 0; i < pts1.size(); ++i)
  {
    W += pts1[i] * pts2[i].transpose();
  }
  Eigen::JacobiSVD<Eigen::Matrix3d> svd(W, Eigen::ComputeFullU | Eigen::ComputeFullV);
  R = svd.matrixU() * svd.matrixV().transpose();
  if (R.determinant() < 0)
    R = -R;

  t = c1 - R * c2;
}

class VertexPose: public g2o::BaseVertex<6, Sophus::SE3d>
{
  public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;

    virtual void setToOriginImpl() override {
      _estimate = Sophus::SE3d();
    }

    virtual void oplusImpl(const double* update) override {
      Eigen::Matrix<double, 6, 1> u;
      u << update[0], update[1], update[2], update[3], update[4], update[5];
      _estimate = Sophus::SE3d::exp(u) * _estimate;
    }

  virtual bool read(istream &in) override {}
  virtual bool write(ostream &out) const override {}
};

class EdgeICP: public g2o::BaseUnaryEdge<3, Eigen::Vector3d, VertexPose>
{
  public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;

    EdgeICP(const Eigen::Vector3d& X2): _X2(X2) {}

    virtual void computeError() override {
      const VertexPose* v = static_cast<VertexPose*>(_vertices[0]);
      Sophus::SE3d Rt = v->estimate();
      _error = measurement() - Rt * _X2;
    }

    virtual void linearizeOplus() override {
      // fill _jacobianOPlusXi;
      const VertexPose* v = static_cast<VertexPose*>(_vertices[0]);
      Sophus::SE3d Rt = v->estimate();
      Eigen::Vector3d X2_transf = Rt * _X2;
      _jacobianOplusXi.block<3, 3>(0, 0) = -Eigen::Matrix3d::Identity();
      _jacobianOplusXi.block<3, 3>(0, 3) = Sophus::SO3d::hat(X2_transf);
    }

    // Auto-diff is not working well. No idea why
    // template <class T>
    // bool operator() (const T* params, T* errors) const {
    //   const Eigen::Map<const Sophus::SE3<T>> Rt_sophus(params);
    //   Eigen::Matrix<T, 3, 1> X2(T(_X2.x()), T(_X2.y()), T(_X2.z()));
    //   Eigen::Matrix<T, 3, 1> X2_transf = Rt_sophus * X2;
    //   errors[0] = T(_measurement[0]) - X2_transf.x();
    //   errors[1] = T(_measurement[1]) - X2_transf.y();
    //   errors[2] = T(_measurement[2]) - X2_transf.z();
    //   return true;
    // }
    // G2O_MAKE_AUTO_AD_FUNCTIONS;

    virtual bool read(std::istream&) override {}
    virtual bool write(std::ostream&) const override {}

  private:
    Eigen::Vector3d _X2;
};

// Solve ICP with non-linear optimization
void bundleAdjustment(
  const vector<Eigen::Vector3d> &pts1,
  const vector<Eigen::Vector3d> &pts2,
  Eigen::Matrix3d &R, Eigen::Vector3d &t) {

    Sophus::SE3d pose_init(R, t);

    typedef g2o::BlockSolverX BlockSolverType;
    typedef g2o::LinearSolverDense<BlockSolverType::PoseMatrixType> LinearSolverType;
    auto solver = new g2o::OptimizationAlgorithmLevenberg(
      g2o::make_unique<BlockSolverType>(g2o::make_unique<LinearSolverType>()));

    g2o::SparseOptimizer optimizer;
    optimizer.setAlgorithm(solver);
    optimizer.setVerbose(true);

    VertexPose *vertex_pose = new VertexPose();
    vertex_pose->setId(0);
    vertex_pose->setEstimate(pose_init);
    optimizer.addVertex(vertex_pose);

    for (int i = 0; i < pts1.size(); ++i)
    {
      EdgeICP *edge = new EdgeICP(pts2[i]);
      edge->setVertex(0, vertex_pose);
      edge->setMeasurement(pts1[i]);
      edge->setInformation(Eigen::Matrix3d::Identity());
      optimizer.addEdge(edge);
    }

    chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
    optimizer.initializeOptimization();
    optimizer.optimize(50);
    chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
    chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
    cout << "optimization ICP costs time: " << time_used.count() << " seconds." << endl;

    Sophus::SE3d optim_pose = vertex_pose->estimate();
    R = optim_pose.so3().unit_quaternion().toRotationMatrix();
    t = optim_pose.translation();
  }

int main(int argc, char **argv) {
  // if (argc != 5) {
  //   cout << "usage: pose_estimation_3d3d img1 img2 depth1 depth2" << endl;
  //   return 1;
  // }
  string f1 = "../1.png"; //argv[1];
  string f2 = "../2.png"; //argv[2];
  string f3 = "../1_depth.png"; //argv[3];
  string f4 = "../2_depth.png"; //argv[3];

  Mat img_1 = imread(f1, CV_LOAD_IMAGE_COLOR);
  Mat img_2 = imread(f2, CV_LOAD_IMAGE_COLOR);

  vector<KeyPoint> keypoints_1, keypoints_2;
  vector<DMatch> matches;
  find_feature_matches(img_1, img_2, keypoints_1, keypoints_2, matches);

  Mat depth1 = imread(f3, CV_LOAD_IMAGE_UNCHANGED);
  Mat depth2 = imread(f4, CV_LOAD_IMAGE_UNCHANGED);
  Mat K = (Mat_<double>(3, 3) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1);
  Eigen::Matrix3d K_eigen;
  K_eigen << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1;
  std::vector<Eigen::Vector3d> pts1, pts2;

  for (DMatch m:matches) {
    ushort d1 = depth1.ptr<unsigned short>(int(keypoints_1[m.queryIdx].pt.y))[int(keypoints_1[m.queryIdx].pt.x)];
    ushort d2 = depth2.ptr<unsigned short>(int(keypoints_2[m.trainIdx].pt.y))[int(keypoints_2[m.trainIdx].pt.x)];
    if (d1 == 0 || d2 == 0)   // bad depth
      continue;
    Point2d p1 = pixel2cam(keypoints_1[m.queryIdx].pt, K);
    Point2d p2 = pixel2cam(keypoints_2[m.trainIdx].pt, K);
    float dd1 = float(d1) / 5000.0;
    float dd2 = float(d2) / 5000.0;
    pts1.push_back(Eigen::Vector3d(p1.x * dd1, p1.y * dd1, dd1));
    pts2.push_back(Eigen::Vector3d(p2.x * dd2, p2.y * dd2, dd2));
  }

  cout << "3d-3d pairs: " << pts1.size() << endl;
  Eigen::Matrix3d R = Eigen::Matrix3d::Identity();
  Eigen::Vector3d t = Eigen::Vector3d::Zero();
  pose_estimation_3d3d(pts1, pts2, R, t);
  cout << "ICP via SVD results: " << endl;
  cout << "R = " << R << endl;
  cout << "t = " << t.transpose() << endl;
  // cout << "R_inv = " << R.t() << endl;
  // cout << "t_inv = " << -R.t() * t << endl;

  // verify p1 = R * p2 + t
  double total_error = 0.0;
  for (int i = 0; i < pts1.size(); i++) {
    total_error += (pts1[i] - (R * pts2[i] + t)).norm();
  }
  std::cout << "Mean error (SVD): " << total_error / pts1.size() << "\n";

  cout << "calling bundle adjustment" << endl;
  R = Eigen::Matrix3d::Identity();
  t = Eigen::Vector3d::Zero();
  bundleAdjustment(pts1, pts2, R, t);

  // verify p1 = R * p2 + t
  total_error = 0.0;
  for (int i = 0; i < pts1.size(); i++) {
    total_error += (pts1[i] - (R * pts2[i] + t)).norm();
  }
  std::cout << "Mean error (BA): " << total_error / pts1.size() << "\n";
}

ext: cpp


text:
//
// $Id: peakpickerqtof.hpp 6337 2014-06-06 20:04:10Z witek96 $
//
//
// Original author: Witold Wolski <wewolski@gmail.com>
//
// Copyright : ETH Zurich
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//

#ifndef PEAKPICKER_H
#define PEAKPICKER_H

#include <boost/math/special_functions.hpp>
#include "pwiz/utility/findmf/base/resample/convert2dense.hpp"
#include "pwiz/utility/findmf/base/filter/filter.hpp"
#include "pwiz/utility/findmf/base/ms/simplepicker.hpp"
#include "pwiz/utility/findmf/base/filter/gaussfilter.hpp"
#include "pwiz/utility/findmf/base/base/interpolate.hpp"
#include "pwiz/utility/findmf/base/resample/determinebinwidth.hpp"
#include "pwiz/utility/findmf/base/base/copyif.hpp"

namespace ralab{
  namespace base{
    namespace ms{

      /// resamples spectrum, apply smoothing,
      /// determines zero crossings,
      /// integrates peaks.

      template<typename TReal>
      struct SimplePeakArea{
        TReal integwith_;

        SimplePeakArea(TReal integwith):integwith_(integwith){}

        /// intagrates the peak intesnities
        template<typename Tzerocross, typename Tintensity, typename Tout>
        void operator()( Tzerocross beginZ,
                         Tzerocross endZ,
                         Tintensity intensity,
                         Tintensity resmpled,
                         Tout area)const
        {
          typedef typename std::iterator_traits<Tout>::value_type AreaType;
          for( ; beginZ != endZ ; ++beginZ , ++area )
          {
            size_t idx = static_cast<size_t>( *beginZ );
            size_t start = static_cast<size_t>( boost::math::round( idx - integwith_ ) );
            size_t end = static_cast<size_t>( boost::math::round( idx + integwith_ + 2.) );
            AreaType aread = 0.;
            for( ; start != end ; ++start )
            {
              aread += *(resmpled + start);
            }
            *area = aread;
          }
        }
      };

      /// extends peak to the left and to the right to the next local minimum or a predefined threshol
      /// or a maximum allowed extension.
      template<typename TReal>
      struct LocalMinPeakArea{
        typedef TReal value_type;
        TReal integwith_;
        TReal threshold_;

        LocalMinPeakArea(TReal integwith,//!<maximal allowed peak width +- in pixel
                         TReal threshold = .1// minimum intensity
            ):integwith_(integwith),threshold_(threshold){}

        /// intagrates the peak intesnities
        template< typename Tzerocross, typename Tintensity, typename Tout >
        void operator()( Tzerocross beginZ,
                         Tzerocross endZ,
                         Tintensity intensity,
                         Tintensity resampled,
                         Tout area) const
        {
          typedef typename std::iterator_traits<Tout>::value_type AreaType;
          for( ; beginZ != endZ ; ++beginZ , ++area )
          {
            size_t idx = static_cast<size_t>( *beginZ );
            size_t start = static_cast<size_t>( boost::math::round( idx - integwith_ ) );
            size_t end = static_cast<size_t>( boost::math::round( idx + integwith_ + 2) );

            Tintensity st = intensity + start;
            Tintensity en = intensity + end;
            Tintensity center = intensity + idx;
            std::ptrdiff_t x1 = std::distance(st, center);
            std::ptrdiff_t y1 = std::distance(center,en);
            mextend(st , en , center);
            std::ptrdiff_t x2 = std::distance(intensity,st);
            std::ptrdiff_t y2 = std::distance(intensity,en);
            std::ptrdiff_t pp = std::distance(st,en);
            AreaType areav = std::accumulate(resampled+x2,resampled+y2,0.);
            *area = areav;
          }
        }

      private:
        ///exend peak to left and rigth
        template<typename TInt >
        void mextend( TInt &start, TInt &end, TInt idx) const
        {
          typedef typename std::iterator_traits<TInt>::value_type Intensitytype;
          //
          for(TInt intens = idx ; intens >= start;  --intens){
            Intensitytype val1 = *intens;
            Intensitytype val2 = *(intens-1);
            if(val1 > threshold_){
              if(val1 < val2 ){
                start = intens;
                break;
              }
            }
            else{
              start = intens;
              break;
            }
          }

          for(TInt intens = idx ; intens <= end;  ++intens){
            Intensitytype val1 = *intens;
            Intensitytype val2 = *(intens+1);
            if(val1 > threshold_){
              if(val1 < val2 ){
                end = intens;
                break;
              }
            }
            else{
              end = intens;
              break;
            }
          }
        }
      };

      /// resamples spectrum, apply smoothing,
      /// determines zero crossings,
      /// integrates peaks.
      template<typename TReal, template <typename B> class TIntegrator >
      struct PeakPicker{
        typedef TReal value_type;
        typedef TIntegrator<value_type> PeakIntegrator;

        TReal resolution_;
        ralab::base::resample::Convert2Dense c2d_; // resamples spectrum
        std::vector<TReal> resampledmz_, resampledintensity_; // keeps result of convert to dense
        std::vector<TReal> filter_, zerocross_, smoothedintensity_; // working variables
        std::vector<TReal> peakmass_, peakarea_; //results
        TReal smoothwith_;
        TReal integrationWidth_;
        ralab::base::ms::SimplePicker<TReal> simplepicker_;
        ralab::base::resample::SamplingWith sw_;
        PeakIntegrator integrator_;
        TReal intensitythreshold_;
        bool area_;
        uint32_t maxnumbersofpeaks_;

        PeakPicker(TReal resolution, //!< instrument resolution
                   std::pair<TReal, TReal> & massrange, //!< mass range of spectrum
                   TReal width = 2., //!< smooth width
                   TReal intwidth = 2., //!< integration width used for area compuation
                   TReal intensitythreshold = 10., // intensity threshold
                   bool area = true,//!< compute area or height? default - height.
                   uint32_t maxnumberofpeaks = 0, //!< maximum of peaks returned by picker
                   double c2d = 1e-5  //!< instrument resampling with small default dissables automatic determination
            ): resolution_(resolution),c2d_( c2d ) ,smoothwith_(width),
          integrationWidth_(intwidth),sw_(),integrator_(integrationWidth_),
          intensitythreshold_(intensitythreshold),area_(area),maxnumbersofpeaks_(maxnumberofpeaks)
        {
          c2d_.defBreak(massrange,ralab::base::resample::resolution2ppm(resolution));
          c2d_.getMids(resampledmz_);
          ralab::base::filter::getGaussianFilterQuantile(filter_,width);
        }

        template<typename Tmass, typename Tintensity>
        void operator()(Tmass begmz, Tmass endmz, Tintensity begint )
        {
          typename std::iterator_traits<Tintensity>::value_type minint = *std::upper_bound(begint,begint+std::distance(begmz,endmz),0.1);

          //determine sampling with
          double a = sw_(begmz,endmz);
          //resmpale the spectrum
          c2d_.am_ = a;
          c2d_.convert2dense(begmz,endmz, begint, resampledintensity_);

          //smooth the resampled spectrum
          ralab::base::filter::filter(resampledintensity_ , filter_ , smoothedintensity_ , true);
          //determine zero crossings
          zerocross_.resize( smoothedintensity_.size()/2 );
          size_t nrzerocross = simplepicker_( smoothedintensity_.begin( ) , smoothedintensity_.end() , zerocross_.begin(), zerocross_.size());

          peakmass_.resize(nrzerocross);
          //determine mass of zerocrossing
          ralab::base::base::interpolate_linear( resampledmz_.begin() , resampledmz_.end() ,
                                                 zerocross_.begin(),  zerocross_.begin()+nrzerocross ,
                                                 peakmass_.begin());

          //determine peak area
          if(area_){
            peakarea_.resize(nrzerocross);
            integrator_( zerocross_.begin(), zerocross_.begin() + nrzerocross ,
                         smoothedintensity_.begin(),resampledintensity_.begin(), peakarea_.begin() );
          }else{
            //determine intensity
            peakarea_.resize(nrzerocross);
            ralab::base::base::interpolate_cubic( smoothedintensity_.begin() , smoothedintensity_.end() ,
                                                  zerocross_.begin(),  zerocross_.begin()+nrzerocross ,
                                                  peakarea_.begin());
          }

          TReal threshold = static_cast<TReal>(minint) * intensitythreshold_;

          if(maxnumbersofpeaks_ > 0){
            double threshmax = getNToppeaks();
            if(threshmax > threshold)
              threshold = threshmax;
          }

          if(threshold > 0.01){
            filter(threshold);
          }
        }

        /// get min instensity of peak to qualify for max-intensity;
        TReal getNToppeaks(){
          TReal intthres  = 0.;
          if(maxnumbersofpeaks_ < peakarea_.size())
          {
            std::vector<TReal> tmparea( peakarea_.begin() , peakarea_.end() );
            std::nth_element(tmparea.begin(),tmparea.end() - maxnumbersofpeaks_ , tmparea.end());
            intthres = *(tmparea.end() - maxnumbersofpeaks_);
          }
          return intthres;
        }

        /// clean the masses using the threshold
        void filter(TReal threshold){
          typename std::vector<TReal>::iterator a = ralab::base::utils::copy_if(peakarea_.begin(),peakarea_.end(),peakmass_.begin(),
                                                                                peakmass_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
          peakmass_.resize(std::distance(peakmass_.begin(),a));
          typename std::vector<TReal>::iterator b = ralab::base::utils::copy_if(peakarea_.begin(),peakarea_.end(),
                                                                                peakarea_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
          peakarea_.resize(std::distance(peakarea_.begin(),b));
          //int x = 1;
        }

        const std::vector<TReal> & getPeakMass(){
          return peakmass_;
        }

        const std::vector<TReal> & getPeakArea(){
          return peakarea_;
        }

        const std::vector<TReal> & getResampledMZ(){
          return resampledmz_;
        }

        const std::vector<TReal> & getResampledIntensity(){
          return resampledintensity_;
        }

        const std::vector<TReal> & getSmoothedIntensity(){
          return smoothedintensity_;
        }
      };
    }//ms
  }//base
}//ralab

#endif // PEAKPICKER_H

ext: hpp

max_stars_repo_name: edyp-lab/pwiz-mzdb

max_issues_repo_name: shze/pwizard-deb

max_forks_repo_name: shze/pwizard-deb