text:
/***************************************************************************
 *   Software License Agreement (BSD License)                              *
 *   Copyright (C) 2012 by Markus Bader <markus.bader@tuwien.ac.at>        *
 *                                                                         *
 *   Redistribution and use in source and binary forms, with or without    *
 *   modification, are permitted provided that the following conditions    *
 *   are met:                                                              *
 *                                                                         *
 *   1. Redistributions of source code must retain the above copyright     *
 *      notice, this list of conditions and the following disclaimer.      *
 *   2. Redistributions in binary form must reproduce the above copyright  *
 *      notice, this list of conditions and the following disclaimer in    *
 *      the documentation and/or other materials provided with the         *
 *      distribution.                                                      *
 *   3. Neither the name of the copyright holder nor the names of its      *
 *      contributors may be used to endorse or promote products derived    *
 *      from this software without specific prior written permission.      *
 *                                                                         *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS   *
 *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT     *
 *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS     *
 *   FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE        *
 *   COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,  *
 *   INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,  *
 *   BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;      *
 *   LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER      *
 *   CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT    *
 *   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY *
 *   WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE           *
 *   POSSIBILITY OF SUCH DAMAGE.                                           *
 ***************************************************************************/

#include <cstdio>
#include <stdlib.h>
#include <iostream>
#include <string>
#include <signal.h>
#include <boost/program_options.hpp>

#include <opencv/cv.h>
#include <opencv/highgui.h>
#include <shmfw/objects/image.h>
#include <shmfw/allocator.h>

#include "shmfw/log.h"

struct Prarmeters {
    bool clear;
    std::string shm_memory_name;
    unsigned int shm_memory_size;
    std::string variable_name;
    std::string configFile;
    std::string image_file;
    int device;
    int height;
    int width;
    bool imshow;
};

int program_options ( int argc, char* argv[], Prarmeters &params ) {
    boost::program_options::options_description desc ( "Allowed Parameters" );
    desc.add_options() //
    ( "help", "get this help message" ) //
    ( "cfg", boost::program_options::value<std::string> ( &params.configFile ), "Config file" ) //
    ( "imshow,s", "opens a highgui window and shows the images" ) //
    ( "device,d", boost::program_options::value<int> ( &params.device )->default_value ( 0 ), "device" )
    ( "image_file,i", boost::program_options::value<std::string> ( &params.image_file )->default_value ( "" ), "image used instead of an device" )
    ( "width,w", boost::program_options::value<int> ( &params.width )->default_value ( 0 ), "image width, zero means auto" )
    ( "height,h", boost::program_options::value<int> ( &params.height )->default_value ( 0 ), "image height, zero means auto" )
    ( "shm_memory_name", boost::program_options::value<std::string> ( &params.shm_memory_name )->default_value ( ShmFw::DEFAULT_SEGMENT_NAME() ), "shared memory segment name" )
    ( "shm_memory_size", boost::program_options::value<unsigned int> ( &params.shm_memory_size )->default_value ( ShmFw::DEFAULT_SEGMENT_SIZE() ), "shared memory segment size" )
    ( "variable_name,v", boost::program_options::value<std::string > ( &params.variable_name )->default_value ( "cvcam" ), "shared variable" );

    std::cout << std::endl;

    boost::program_options::variables_map vm;
    try {
        boost::program_options::store ( boost::program_options::parse_command_line ( argc, argv, desc ), vm );
    } catch ( const std::exception& ex ) {
        std::cout << desc << "\n";
        return 1;
    }
    boost::program_options::notify ( vm );

    if ( vm.count ( "help" ) ) {
        std::cout << desc << "\n";
        exit ( 0 );
    }

    if ( !params.configFile.empty() ) {
        std::ifstream file ( params.configFile.c_str(), std::ifstream::in );
        if ( file.is_open() ) {
            try {
                boost::program_options::store ( boost::program_options::parse_config_file ( file, desc ), vm );
                boost::program_options::notify ( vm );
            } catch ( const std::exception& ex ) {
                std::cout << "Error reading config file: " << ex.what() << std::endl;
                exit ( 0 );
            }
        } else {
            std::cout << "Error opening config file " << params.configFile << std::endl;
            exit ( 0 );
        }
    }
    params.imshow = ( vm.count ( "imshow" ) > 0 );
    return 0;
}

bool loop_program = true;

void terminate ( int param ) {
    std::cout << "Program was trying to abort or terminate." << std::endl;
    loop_program = false;
}

int main ( int argc, char *argv[] ) {

    signal ( SIGABRT,	terminate );
    signal ( SIGTERM,	terminate );
    Prarmeters params;
    program_options ( argc, argv, params );

    cv::Point x;

    cv::VideoCapture cap;
    cv::Mat imgSrc;
    if ( params.image_file.empty() ) {
        cap.open ( params.device );
        if ( !cap.isOpened() ) {
            std::cout << "***Could not initialize capturing...***\n";
            exit ( 0 );
        }
        if ( params.width > 0 ) cap.set ( CV_CAP_PROP_FRAME_WIDTH, params.width );
        if ( params.height > 0 ) cap.set ( CV_CAP_PROP_FRAME_HEIGHT, params.height );
        std::cout << "Video " <<
                  ": width=" << cap.get ( CV_CAP_PROP_FRAME_WIDTH ) <<
                  ", height=" << cap.get ( CV_CAP_PROP_FRAME_HEIGHT ) << std::endl;
        if ( cap.isOpened() ) {
            cap >> imgSrc;
        } else {
            std::cout << "***Could not get any image...***\n";
            exit ( 0 );
        }
    } else {
        std::cout << "using file: " << params.image_file << std::endl;
        imgSrc = cv::imread ( params.image_file );
    }
    if ( imgSrc.empty() ) {
        std::cout << "***image...***\n";
        return -1;
    }

    ShmFw::HandlerPtr shmHdl = ShmFw::Handler::create ( params.shm_memory_name, params.shm_memory_size );
    ShmFw::Alloc<ShmFw::ImageShm> shmImg ( params.variable_name, shmHdl );
    shmImg.unlock();
    cv::Mat shmMat;
    shmImg->copyFrom(imgSrc);
    shmImg->cvMat ( shmMat );
    std::cout << "You should be able to view the image with:  ";
    std::cout << "shmfw-image_view -v " << params.variable_name << std::endl;
    do {
        shmImg.lock();
        if ( cap.isOpened() ) {
            cap >> shmMat;
        }
        if ( shmMat.empty() ) {
            continue;
        }
        shmImg.unlock();
        shmImg.itHasChanged();
        if ( params.imshow ) {
            cv::imshow ( "image", shmMat );
        }
    } while ( ( cv::waitKey ( 30 ) < 0 ) && loop_program ) ;
    // the camera will be deinitialized automatically in VideoCapture destructor
    return 0;
}

max_stars_repo_name: ShmFw/shmfw

max_issues_repo_name: ShmFw/shmfw

max_forks_repo_name: ShmFw/shmfw


text:
#pragma once

#include <boost/polygon/polygon.hpp>
#include <boost/polygon/voronoi.hpp>
#include <cslibs_vectormaps/dxf/dxf_map.h>

namespace boost {
namespace polygon {
using PointType = cslibs_vectormaps::dxf::DXFMap::Point;
using SegmentType = cslibs_vectormaps::dxf::DXFMap::Vector;
template<>
struct geometry_concept<PointType> {
    typedef point_concept type;
};
template<>
struct point_traits<PointType> {
    typedef double coordinate_type;

    static inline coordinate_type get(const PointType &point, orientation_2d orient)
    {
        return (orient == HORIZONTAL) ? point.x() : point.y();
    }
};

template<>
struct geometry_concept<SegmentType> {
    typedef segment_concept type;
};
template<>
struct segment_traits<SegmentType> {
    typedef double coordinate_type;
    typedef PointType point_type;

    static inline point_type get(const SegmentType &segment, direction_1d dir)
    {
        return dir.to_int() ? segment.second : segment.first;
    }
};
}
}

typedef boost::polygon::voronoi_diagram<double> VoronoiType;

ext: hpp


text:
// Copyright (c) 2014, German Neuroinformatics Node (G-Node)
//
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted under the terms of the BSD License. See
// LICENSE source in the root of the Project.

#include <boost/python.hpp>
#include <boost/optional.hpp>

#include <nix.hpp>
#include <nix/util/util.hpp>

#include <accessors.hpp>
#include <transmorgify.hpp>
#include <docstrings.hpp>

#include <PyUtil.hpp>

using namespace nix;
using namespace boost::python;

namespace nixpy {

struct UnitWrap {};

bool isScalableSingleUnit(const std::string &unitA, const std::string &unitB) {
    return util::isScalable(unitA, unitB);
}

bool isScalableMultiUnits(const std::vector<std::string> &unitsA, const std::vector<std::string> &unitsB) {
    return util::isScalable(unitsA, unitsB);
}

std::vector<std::string> splitUnit(const std::string &unit) {
    std::string si, prefix, power;
    util::splitUnit(unit, prefix, si, power);
    return std::vector<std::string>{prefix, si, power};
}

std::vector<std::string> splitCompound(const std::string &unit) {
   std::vector<std::string> parts;
   util::splitCompoundUnit(unit, parts);
   return parts;
}

struct NameWrap {};

void PyUtil::do_export() {
  class_<UnitWrap> ("units")
  .def("sanitizer", util::unitSanitizer, doc::unit_sanitizer).staticmethod("sanitizer")
  .def("is_si", util::isSIUnit, doc::unit_is_si).staticmethod("is_si")
  .def("is_atomic", util::isAtomicSIUnit, doc::unit_is_atomic).staticmethod("is_atomic")
  .def("is_compound", util::isCompoundSIUnit, doc::unit_is_compound).staticmethod("is_compound")
  .def("scalable", &isScalableMultiUnits, doc::unit_scalable).staticmethod("scalable")
  .def("scaling", util::getSIScaling, doc::unit_scaling).staticmethod("scaling")
  .def("split", &splitUnit, doc::unit_split).staticmethod("split")
  .def("split_compound", &splitCompound, doc::unit_compound_split).staticmethod("split_compound")
  ;

  class_<NameWrap> ("names")
  .def("sanitizer", util::nameSanitizer, doc::name_sanitizer).staticmethod("sanitizer")
  .def("check", util::nameCheck, doc::name_check).staticmethod("check")
  .def("create_id", util::createId, doc::create_id).staticmethod("create_id")
  ;
}

}

ext: cpp

max_stars_repo_name: gicmo/nixpy

max_issues_repo_name: gicmo/nixpy

max_forks_repo_name: gicmo/nixpy


text:
// TRENTO: Reduced Thickness Event-by-event Nuclear Topology
// Copyright 2015 Jonah E. Bernhard, J. Scott Moreland
// MIT License

#include "nucleus.h"

#include <algorithm>
#include <cmath>
#include <memory>
#include <stdexcept>
#include <string>
#include <utility>

#include <boost/math/constants/constants.hpp>
#ifdef TRENTO_HDF5
// include multi_array for use with ManualNucleus
#ifdef NDEBUG
#define BOOST_DISABLE_ASSERTS
#endif
#include <boost/multi_array.hpp>
#endif

#include "hdf5_utils.h"
#include "random.h"

namespace trento {

NucleusPtr Nucleus::create(const std::string& species, double nucleon_dmin) {
  // W-S params ref. in header
  // XXX: remember to add new species to the help output in main() and the readme
  if (species == "p")
    return NucleusPtr{new Proton{}};
  else if (species == "d")
    return NucleusPtr{new Deuteron{}};
  else if (species == "Cu")
    return NucleusPtr{new WoodsSaxonNucleus{
       63, 4.20, 0.596, 0.0, nucleon_dmin
    }};
  else if (species == "Cu2")
    return NucleusPtr{new DeformedWoodsSaxonNucleus{
       63, 4.20, 0.596, 0.162, -0.006, nucleon_dmin
    }};
  else if (species == "O")
      return NucleusPtr{new WoodsSaxonNucleus{
              16, 2.608, 0.513, -0.051, nucleon_dmin
      }};
  else if (species == "Xe")
    return NucleusPtr{new WoodsSaxonNucleus{
      129, 5.36, 0.590, 0.0, nucleon_dmin
    }};
  else if (species == "Xe2")
    return NucleusPtr{new DeformedWoodsSaxonNucleus{
      129, 5.36, 0.590, 0.162, -0.003, nucleon_dmin
    }};
  else if (species == "XeTri")
    return NucleusPtr{new DeformedWoodsSaxonNucleus{
     129, 5.60, 0.49, 0.21, 0.0, nucleon_dmin
    }};
    else if (species == "XeALICE")
    return NucleusPtr{new DeformedWoodsSaxonNucleus{
     129, 5.36, 0.590, 0.18, 0.0, nucleon_dmin
    }};
  else if (species == "Au")
    return NucleusPtr{new WoodsSaxonNucleus{
      197, 6.38, 0.535, 0.0, nucleon_dmin
    }};
  else if (species == "Au2")
    return NucleusPtr{new DeformedWoodsSaxonNucleus{
      197, 6.38, 0.535, -0.131, -0.031, nucleon_dmin
    }};
  else if (species == "Pb")
    return NucleusPtr{new WoodsSaxonNucleus{
      208, 6.62, 0.546, 0.0, nucleon_dmin
    }};
  else if (species == "PbTri")
    return NucleusPtr{new WoodsSaxonNucleus{
      208, 6.65, 0.54, 0.0, nucleon_dmin
    }};
  else if (species == "U")
    return NucleusPtr{new DeformedWoodsSaxonNucleus{
      238, 6.81, 0.600, 0.280, 0.093, nucleon_dmin
    }};
  else if (species == "U2")
    return NucleusPtr{new DeformedWoodsSaxonNucleus{
      238, 6.86, 0.420, 0.265, 0.000, nucleon_dmin
    }};
  else if (species == "U3")
    return NucleusPtr{new DeformedWoodsSaxonNucleus{
      238, 6.67, 0.440, 0.280, 0.093, nucleon_dmin
    }};
  // Read nuclear configurations from HDF5.
  else if (hdf5::filename_is_hdf5(species)) {
#ifdef TRENTO_HDF5
    return ManualNucleus::create(species);
#else
    throw std::invalid_argument{"HDF5 output was not compiled"};
#endif  // TRENTO_HDF5
  }
  else
    throw std::invalid_argument{"unknown projectile species: " + species};
}

Nucleus::Nucleus(std::size_t A) : nucleons_(A), offset_(0) {}

void Nucleus::sample_nucleons(double offset) {
  offset_ = offset;
  sample_nucleons_impl();
}

void Nucleus::set_nucleon_position(
    NucleonData& nucleon, double x, double y, double z) {
  nucleon.set_position(x + offset_, y, z);
}

Proton::Proton() : Nucleus(1) {}

/// Always zero.
double Proton::radius() const {
  return 0.;
}

/// Always place the nucleon at the origin.
void Proton::sample_nucleons_impl() {
  set_nucleon_position(*begin(), 0., 0., 0.);
}

// Without loss of generality, let the internal a_ parameter be the minimum of
// the given (a, b) and the internal b_ be the maximum.
Deuteron::Deuteron(double a, double b)
    : Nucleus(2),
      a_(std::fmin(a, b)),
      b_(std::fmax(a, b))
{}

double Deuteron::radius() const {
  // The quantile function for the exponential distribution exp(-2*a*r) is
  // -log(1-q)/(2a).  Return the 99% quantile.
  return -std::log(.01)/(2*a_);
}

void Deuteron::sample_nucleons_impl() {
  // Sample the inter-nucleon radius using rejection sampling with an envelope
  // function.  The Hulthén wavefunction including the r^2 Jacobian expands to
  // three exponential terms:  exp(-2*a*r) + exp(-2*b*r) - 2*exp(-(a+b)*r).
  // This does not have a closed-form inverse CDF, however we can easily sample
  // exponential numbers from the term that falls off the slowest, i.e.
  // exp(-2*min(a,b)*r).  In the ctor initializer list the "a" parameter is
  // always set to the minimum, so we should sample from exp(-2*a*r).
  double r, prob;
  do {
    // Sample a uniform random number, u = exp(-2*a*r).
    auto u = random::canonical<double>();
    // Invert to find the actual radius.
    r = -std::log(u) / (2*a_);
    // The acceptance probability is now the radial wavefunction over the
    // envelope function, both evaluated at the proposal radius r.
    // Conveniently, the envelope evaluated at r is just the uniform random
    // number u.
    prob = std::pow(std::exp(-a_*r) - std::exp(-b_*r), 2) / u;
  } while (prob < random::canonical<double>());

  // Now sample spherical rotation angles.
  auto cos_theta = random::cos_theta<double>();
  auto phi = random::phi<double>();

  // And compute the Cartesian coordinates of one nucleon.
  auto r_sin_theta = r * std::sqrt(1. - cos_theta*cos_theta);
  auto x = r_sin_theta * std::cos(phi);
  auto y = r_sin_theta * std::sin(phi);
  auto z = r * cos_theta;

  // Place the first nucleon at the sampled coordinates.
  set_nucleon_position(*begin(), x, y, z);
  // Place the second nucleon opposite to the first.
  set_nucleon_position(*std::next(begin()), -x, -y, -z);
}

MinDistNucleus::MinDistNucleus(std::size_t A, double dmin)
    : Nucleus(A),
      dminsq_(dmin*dmin)
{}

bool MinDistNucleus::is_too_close(const_iterator nucleon) const {
  if (dminsq_ < 1e-10)
    return false;
  for (const_iterator nucleon2 = begin(); nucleon2 != nucleon; ++nucleon2) {
    auto dx = nucleon->x() - nucleon2->x();
    auto dy = nucleon->y() - nucleon2->y();
    auto dz = nucleon->z() - nucleon2->z();
    if (dx*dx + dy*dy + dz*dz < dminsq_)
      return true;
  }
  return false;
}

// Extend the W-S dist out to R + 10a; for typical values of (R, a), the
// probability of sampling a nucleon beyond this radius is O(10^-5).
WoodsSaxonNucleus::WoodsSaxonNucleus(
    std::size_t A, double R, double a, double w, double dmin)
    : MinDistNucleus(A, dmin),
      R_(R),
      a_(a),
      woods_saxon_dist_(1000, 0., R + 10.*a,
        [R, a, w](double r) { return r*r*(1+w*(r/R)*(r/R))/(1.+std::exp((r-R)/a)); })
{}

/// Return something a bit smaller than the true maximum radius.  The
/// Woods-Saxon distribution falls off very rapidly (exponentially), and since
/// this radius determines the impact parameter range, the true maximum radius
/// would cause far too many events with zero participants.
double WoodsSaxonNucleus::radius() const {
  return R_ + 3.*a_;
}

/// Sample Woods-Saxon nucleon positions.
void WoodsSaxonNucleus::sample_nucleons_impl() {
  // When placing nucleons with a minimum distance criterion, resample spherical
  // angles until the nucleon is not too close to a previously sampled nucleon,
  // but do not resample radius -- this could modify the Woods-Saxon dist.

  // Because of the r^2 Jacobian, there is less available space at smaller
  // radii.  Therefore, pre-sample all radii first, sort them, and then place
  // nucleons starting with the smallest radius and working outwards.  This
  // dramatically reduces the chance that a nucleon cannot be placed.
  std::vector<double> radii(size());
  for (auto&& r : radii)
    r = woods_saxon_dist_(random::engine);
  std::sort(radii.begin(), radii.end());

  // Place each nucleon at a pre-sampled radius.
  auto r_iter = radii.cbegin();
  for (iterator nucleon = begin(); nucleon != end(); ++nucleon) {
    // Get radius and advance iterator.
    auto& r = *r_iter++;

    // Sample angles until the minimum distance criterion is satisfied.
    auto ntries = 0;
    do {
      // Sample isotropic spherical angles.
      auto cos_theta = random::cos_theta<double>();
      auto phi = random::phi<double>();

      // Convert to Cartesian coordinates.
      auto r_sin_theta = r * std::sqrt(1. - cos_theta*cos_theta);
      auto x = r_sin_theta * std::cos(phi);
      auto y = r_sin_theta * std::sin(phi);
      auto z = r * cos_theta;

      set_nucleon_position(*nucleon, x, y, z);

      // Retry sampling a reasonable number of times.  If a nucleon cannot be
      // placed, give up and leave it at its last sampled position.  Some
      // approximate numbers for Pb nuclei:
      //
      //   dmin = 0.5 fm, < 0.001% of nucleons cannot be placed
      //          1.0 fm, ~0.005%
      //          1.5 fm, ~0.1%
      //          1.73 fm, ~1%
    } while (++ntries < 1000 && is_too_close(nucleon));
  }
  // XXX: re-center nucleon positions?
}

// Set rmax like the non-deformed case (R + 10a), but for the maximum
// "effective" radius.  The numerical coefficients for beta2 and beta4 are the
// approximate values of Y20 and Y40 at theta = 0.
DeformedWoodsSaxonNucleus::DeformedWoodsSaxonNucleus(
    std::size_t A, double R, double a, double beta2, double beta4, double dmin)
    : MinDistNucleus(A, dmin),
      R_(R),
      a_(a),
      beta2_(beta2),
      beta4_(beta4),
      rmax_(R*(1. + .63*std::fabs(beta2) + .85*std::fabs(beta4)) + 10.*a)
{}

/// Return something a bit smaller than the true maximum radius.  The
/// Woods-Saxon distribution falls off very rapidly (exponentially), and since
/// this radius determines the impact parameter range, the true maximum radius
/// would cause far too many events with zero participants.
double DeformedWoodsSaxonNucleus::radius() const {
  return rmax_ - 7.*a_;
}

double DeformedWoodsSaxonNucleus::deformed_woods_saxon_dist(
    double r, double cos_theta) const {
  auto cos_theta_sq = cos_theta*cos_theta;

  // spherical harmonics
  using math::double_constants::one_div_root_pi;
  auto Y20 = std::sqrt(5)/4. * one_div_root_pi * (3.*cos_theta_sq - 1.);
  auto Y40 = 3./16. * one_div_root_pi *
             (35.*cos_theta_sq*cos_theta_sq - 30.*cos_theta_sq + 3.);

  // "effective" radius
  auto Reff = R_ * (1. + beta2_*Y20 + beta4_*Y40);

  return 1. / (1. + std::exp((r - Reff) / a_));
}

/// Sample deformed Woods-Saxon nucleon positions.
void DeformedWoodsSaxonNucleus::sample_nucleons_impl() {
  // The deformed W-S distribution is defined so the symmetry axis is aligned
  // with the Z axis, so e.g. the long axis of uranium coincides with Z.
  //
  // After sampling positions, they must be randomly rotated.  In general this
  // requires three Euler rotations, but in this case we only need two
  // because there is no use in rotating about the nuclear symmetry axis.
  //
  // The two rotations are:
  //  - a polar "tilt", i.e. rotation about the X axis
  //  - an azimuthal "spin", i.e. rotation about the original Z axis

  // "tilt" angle
  const auto cos_a = random::cos_theta<double>();
  const auto sin_a = std::sqrt(1. - cos_a*cos_a);

  // "spin" angle
  const auto angle_b = random::phi<double>();
  const auto cos_b = std::cos(angle_b);
  const auto sin_b = std::sin(angle_b);

  // Pre-sample and sort (r, cos_theta) points from the deformed W-S dist.
  // See comments in WoodsSaxonNucleus (above) for rationale.
  struct Sample {
    double r, cos_theta;
  };

  std::vector<Sample> samples(size());

  for (auto&& sample : samples) {
    // Sample (r, cos_theta) using a standard rejection method.
    // Remember to include the phase-space factors.
    do {
      sample.r = rmax_ * std::cbrt(random::canonical<double>());
      sample.cos_theta = random::cos_theta<double>();
    } while (
      random::canonical<double>() >
      deformed_woods_saxon_dist(sample.r, sample.cos_theta)
    );
  }

  // Sort by radius.  Could also sort by e.g. the perpendicular distance from
  // the z-axis, or by descending W-S density.  Empirically, radius leads to the
  // smallest failure rate.
  std::sort(
    samples.begin(), samples.end(),
    [](const Sample& a, const Sample& b) {
      return a.r < b.r;
    }
  );

  // Place each nucleon at a pre-sampled (r, cos_theta).
  auto sample = samples.cbegin();
  for (iterator nucleon = begin(); nucleon != end(); ++nucleon, ++sample) {
    auto& r = sample->r;
    auto& cos_theta = sample->cos_theta;

    auto r_sin_theta = r * std::sqrt(1. - cos_theta*cos_theta);
    auto z = r * cos_theta;

    // Sample azimuthal angle until the minimum distance criterion is satisfied.
    auto ntries = 0;
    do {
      // Choose azimuthal angle.
      auto phi = random::phi<double>();

      // Convert to Cartesian coordinates.
      auto x = r_sin_theta * std::cos(phi);
      auto y = r_sin_theta * std::sin(phi);

      // Rotate.
      // The rotation formula was derived by composing the "tilt" and "spin"
      // rotations described above.
      auto x_rot = x*cos_b - y*cos_a*sin_b + z*sin_a*sin_b;
      auto y_rot = x*sin_b + y*cos_a*cos_b - z*sin_a*cos_b;
      auto z_rot =           y*sin_a       + z*cos_a;

      set_nucleon_position(*nucleon, x_rot, y_rot, z_rot);

      // In addition to resampling phi, flip the z-coordinate each time.  This
      // works because the deformed WS dist is symmetric in z.  Effectively
      // doubles the available space for the nucleon.
      z *= -1;

      // Retry a reasonable number of times.  Unfortunately the failure rate is
      // worse than non-deformed sampling because there is less freedom to place
      // each nucleon.  Some approximate numbers for U nuclei:
      //
      //   dmin = 0.5 fm, < 0.001% of nucleons cannot be placed
      //          1.0 fm, ~0.03%
      //          1.3 fm, ~0.3%
      //          1.5 fm, ~1.2%
    } while (++ntries < 1000 && is_too_close(nucleon));
  }
}

#ifdef TRENTO_HDF5

namespace {

// Read a slice of an HDF5 dataset into a boost::multi_array.
template <typename T, std::size_t FileDims, std::size_t MemDims>
boost::multi_array<T, MemDims>
read_dataset(
    const H5::DataSet& dataset,
    const std::array<hsize_t, FileDims>& count,
    const std::array<hsize_t, FileDims>& start,
    const std::array<hsize_t, MemDims>& shape) {
  boost::multi_array<T, MemDims> array{shape};
  auto filespace = dataset.getSpace();
  filespace.selectHyperslab(H5S_SELECT_SET, count.data(), start.data());
  auto memspace = hdf5::make_dataspace(shape);

  dataset.read(array.data(), hdf5::type<T>(), memspace, filespace);

  return array;
}

}  // unnamed namespace

std::unique_ptr<ManualNucleus> ManualNucleus::create(const std::string& path) {
  auto file = hdf5::try_open_file(path);

  // Check that there is a single dataset in the file.
  // Might relax this constraint in the future.
  if (file.getNumObjs() != 1)
    throw std::invalid_argument{
      "file '" + path + "' must contain exactly one object"
    };

  auto name = file.getObjnameByIdx(0);
#if H5_VERSION_GE(1, 8, 13)
  if (file.childObjType(name) != H5O_TYPE_DATASET)  // added v1.8.13
#else
  if (file.getObjTypeByIdx(0) != H5G_DATASET)  // deprecated fall back
#endif
    throw std::invalid_argument{
      "object '" + name + "' in file '" + path + "' is not a dataset"
    };

  // Make dataset object in a unique_ptr for eventual passing to ctor.
  auto dataset = std::unique_ptr<H5::DataSet>{
    new H5::DataSet{file.openDataSet(name)}
  };

  // Verify that the dataset has the correct dimensionality and shape.
  std::array<hsize_t, 3> shape;
  auto ndim = dataset->getSpace().getSimpleExtentDims(shape.data());

  if (ndim != 3)
    throw std::invalid_argument{
      "dataset '" + name + "' in file '" + path + "' has " +
      std::to_string(ndim) + " dimensions (need 3)"
    };

  if (shape[2] != 3)
    throw std::invalid_argument{
      "dataset '" + name + "' in file '" + path + "' has " +
      std::to_string(shape[2]) + " columns (need 3)"
    };

  // Deduce number of configs and number of nucleons (A) from the shape.
  const auto& nconfigs = shape[0];
  const auto& A = shape[1];

  // Estimate the max radius from at least 500 nucleon positions.
  auto n = std::min(500/A + 1, nconfigs);
  std::array<hsize_t, 3> count = {n, A, 3};
  std::array<hsize_t, 3> start = {0, 0, 0};
  std::array<hsize_t, 2> shape_n = {n*A, 3};
  auto positions = read_dataset<float>(*dataset, count, start, shape_n);

  auto rmax_sq = 0.;

  for (const auto& position : positions) {
    auto& x = position[0];
    auto& y = position[1];
    auto& z = position[2];
    auto r_sq = x*x + y*y + z*z;
    if (r_sq > rmax_sq)
      rmax_sq = r_sq;
  }

  auto rmax = std::sqrt(rmax_sq);

  return std::unique_ptr<ManualNucleus>{
    new ManualNucleus{std::move(dataset), nconfigs, A, rmax}
  };
}

ManualNucleus::ManualNucleus(std::unique_ptr<H5::DataSet> dataset,
                             std::size_t nconfigs, std::size_t A, double rmax)
    : Nucleus(A),
      dataset_(std::move(dataset)),
      rmax_(rmax),
      index_dist_(0, nconfigs - 1)
{}

ManualNucleus::~ManualNucleus() = default;

double ManualNucleus::radius() const {
  return rmax_;
}

void ManualNucleus::sample_nucleons_impl() {
  // Sample Euler rotation angles.
  // First is an azimuthal spin about the Z axis.
  const auto angle_1 = random::phi<double>();
  const auto c1 = std::cos(angle_1);
  const auto s1 = std::sin(angle_1);
  // Then a polar tilt about the original X axis, uniform in cos(theta).
  const auto c2 = random::cos_theta<double>();
  const auto s2 = std::sqrt(1. - c2*c2);
  // Finally another azimuthal spin about the original Z axis.
  const auto angle_3 = random::phi<double>();
  const auto c3 = std::cos(angle_3);
  const auto s3 = std::sin(angle_3);

  // Choose and read a random config from the dataset.
  std::array<hsize_t, 3> count = {1, size(), 3};
  std::array<hsize_t, 3> start = {index_dist_(random::engine), 0, 0};
  std::array<hsize_t, 2> shape = {size(), 3};
  const auto positions = read_dataset<float>(*dataset_, count, start, shape);

  // Loop over positions and nucleons.
  auto positions_iter = positions.begin();
  for (iterator nucleon = begin(); nucleon != end(); ++nucleon) {
    // Extract position vector and increment iterator.
    auto position = *positions_iter++;
    auto& x = position[0];
    auto& y = position[1];
    auto& z = position[2];

    // Rotate.
    auto x_rot = x*(c1*c3 - c2*s1*s3) - y*(c3*s1 + c1*c2*s3) + z*s2*s3;
    auto y_rot = x*(c1*s3 + c2*c3*s1) - y*(s1*s3 - c1*c2*c3) - z*c3*s2;
    auto z_rot = x*s1*s2              + y*c1*s2              + z*c2;

    set_nucleon_position(*nucleon, x_rot, y_rot, z_rot);
  }
}

#endif  // TRENTO_HDF5

}  // namespace trento

ext: cxx


text:
#include <vector>
#include <iterator>
#include <typeinfo>
#include <openssl/crypto.h>

#include <boost/algorithm/string.hpp>

#include "user.h"
#include "utils.h"

User::User(unsigned int keysize)
{
    this->keysize = keysize;
    this->numberG = BN_new();
    BIGNUM * number1 = BN_new();
    BN_one(number1);

    BN_add(this->numberG, number1, number1);

    BN_clear_free(number1);
};

User::~User(void)
{
    if (this->numberX_1)
        BN_clear_free(this->numberX_1);
    if (this->numberS)
        BN_clear_free(this->numberX_1);
    if (this->numberG)
        BN_clear_free(this->numberG);
    if (this->primeQ)
        BN_clear_free(this->primeQ);
}

void User::setUserKey(const BIGNUM * userKey)
{
    this->numberX_1 = BN_dup(userKey);
};

void User::setNumberG(const BIGNUM * serverKey)
{
    this->numberG = BN_dup(serverKey);
};

void User::setPrimeQ(const BIGNUM * Q)
{
    this->primeQ = BN_dup(Q);
};

void User::setSalt(const BIGNUM * saltKey)
{
    this->numberS = BN_dup(saltKey);
};

std::pair<BIGNUM *, BIGNUM *> User::UserTD(const unsigned char * component, unsigned int size_component)
{
    if(this->primeQ == NULL || this->numberG == NULL)
    {
        printf("Values not initialized\n");
        exit(-1);
    }
    /* Convert the user string into a number */
    BIGNUM * componentInBinary = BN_new();
    if (size_component > 160) // TODO: FIX THIS LIMIT ACCORDING TO THE KEYSIZE, SHOW A WARNING MESSAGE!
        size_component = 160;
    componentInBinary = BN_bin2bn(component, size_component, componentInBinary);

    BIGNUM * r = BN_new();
    BN_CTX * ctx = BN_CTX_new();
    BN_div(NULL, r, componentInBinary, this->primeQ, ctx);
    BN_CTX_free(ctx);

    BN_clear_free(componentInBinary);

    /* TODO: use a pseudorandom function */

    /* Calculate first part of the pair */
    BIGNUM * t_1 = BN_new();
    ctx = BN_CTX_new();
    BN_mod_exp(t_1, this->numberG, r, this->primeQ, ctx);
    BN_CTX_free(ctx);

    /* Calculate second part of the pair */
    BIGNUM * t_2 = BN_new();
    ctx = BN_CTX_new();
    BN_mod_exp(t_2, t_1, this->numberX_1, this->primeQ, ctx);
    BN_CTX_free(ctx);

    BN_clear_free(r);

    /* Return the key pair */

    return std::make_pair(t_1, t_2);
};

char * User::FullUserTD(const std::vector<std::pair<BIGNUM *, BIGNUM *> > * user_content_names)
{

    char * new_content_name;
    unsigned int offset = 0;
    new_content_name = (char *) calloc(sizeof(char), user_content_names->size() * ((SECURITY_KEYSIZE/4)+1 + 1) * 2); // 1 for the / 1 for the \0

    std::vector<std::pair<BIGNUM *, BIGNUM *> >::const_iterator it;

    /* Initialize name */
    memcpy((void *) (new_content_name + offset), (void *) "/", 1 + 1);
    offset +=1;

    for (it = user_content_names->begin() ; it != user_content_names->end(); it++)
    {
        std::pair<BIGNUM *, BIGNUM *> pair_name = *it;
        //res->push_back(encrypted_bignum);

        //TODO: me puedo estar olvidando un cero aca
        char * aux_component = BN_bn2hex(pair_name.first);
        unsigned int size_aux_component = strnlen(aux_component, SECURITY_KEYSIZE/4 + 1 + 1); // Potentil \0 and potential minus sign
        memcpy((void *) (new_content_name + offset), (void *) aux_component, size_aux_component);
        offset += size_aux_component;
        memcpy((void *) (new_content_name + offset), (void *) "/", 1 + 1);
        offset +=1;
        OPENSSL_free(aux_component);

        aux_component = BN_bn2hex(pair_name.second);
        size_aux_component = strnlen(aux_component, SECURITY_KEYSIZE/4 + 1 + 1); // Potentil \0 and potential minus sign
        memcpy((void *) (new_content_name + offset), (void *) aux_component, size_aux_component);
        offset +=(size_aux_component);
        memcpy((void *) (new_content_name + offset), (void *) "/", 1 + 1);
        offset +=1;

        OPENSSL_free(aux_component);
    }

    return new_content_name;
}

// TODO: untested
std::pair<BIGNUM *, BIGNUM *> User::ContentProviderEnc(const unsigned char * component, unsigned int size_component)
{
    if(this->primeQ == NULL || this->numberG == NULL)
    {
        printf("Values not initialized\n");
        exit(-1);
    }

    /* Convert the user string into a number */
    BIGNUM * componentInBinary = BN_new();
    if (size_component > 160)
        size_component = 160;
    componentInBinary = BN_bin2bn(component, size_component, componentInBinary);

	// Calculate random r
	BIGNUM * r = BN_new();
    BN_rand(r, this->keysize, -1, 0);

	// Calculate g^r
	BIGNUM * t_1 = BN_new();
	BN_CTX * ctx = BN_CTX_new();
    BN_mod_exp(t_1, this->numberG, r, this->primeQ, ctx);
    BN_CTX_free(ctx);

	// Calculate (g^r)^Xuser
	BIGNUM * aux = BN_new();
	ctx = BN_CTX_new();
	BN_mod_exp(aux, t_1, this->numberX_1, this->primeQ, ctx);
	BN_CTX_free(ctx);

	// Multiply the result
	BIGNUM * t_2 = BN_new();
	ctx = BN_CTX_new();
	BN_mod_mul(t_2, aux, componentInBinary, this->primeQ, ctx);
	BN_CTX_free(ctx);

	BN_clear_free(aux);
	BN_clear_free(r);
	BN_clear_free(componentInBinary);

	// Create pair (g^r, g^{r \times Xserver} \times d)
	return std::make_pair(t_1, t_2);
};

BIGNUM * User::UserDec(const std::pair<BIGNUM *, BIGNUM *> e)
{

	if(this->primeQ == NULL || this->numberG == NULL || this->numberX_1 == NULL)
    {
        printf("Values not initialized\n");
        exit(-1);
    }

	BIGNUM * aux = BN_new();
	BN_CTX * ctx = BN_CTX_new();
	BN_mod_exp(aux, e.first, this->numberX_1, this->primeQ, ctx);
	BN_CTX_free(ctx);

	ctx = BN_CTX_new();
	BIGNUM * aux2 = BN_mod_inverse(NULL, aux, this->primeQ, ctx);
	BN_CTX_free(ctx);

	BIGNUM * res = BN_new();
	ctx = BN_CTX_new();
	BN_mod_mul(res,  aux2, e.second, this->primeQ, ctx);
	BN_CTX_free(ctx);

	BN_clear_free(aux);
	BN_clear_free(aux2);

	return res;
}

/* WRAPPER FOR C */
extern "C" char * EncryptUserName(const char * original_content_name)
{
    User *u = new User(512); //TODO: change

    char * res_getUserK = retrieveKeyFromServer("userK:0");
    char * res_getPrimeQ = retrieveKeyFromServer("primeQ");

    BIGNUM * userKey = BN_new();
    BIGNUM * primeQ = BN_new();

    BN_hex2bn(&primeQ, res_getPrimeQ);
    BN_hex2bn(&userKey, res_getUserK);
    u->setUserKey(userKey);
    u->setPrimeQ(primeQ);

    BN_clear_free(primeQ);
    BN_clear_free(userKey);

    free(res_getPrimeQ);
    free(res_getUserK);

    /* split in keywords */
    std::vector<std::string> keywords;
    boost::split(keywords, original_content_name, boost::is_any_of("/"));

    std::vector<std::pair<BIGNUM *, BIGNUM *> > * vec = new std::vector<std::pair<BIGNUM *, BIGNUM *> >();
    for(std::vector<std::string>::iterator it = keywords.begin(); it != keywords.end(); it++)
    {
        std::string t = *it;
        if (t.compare("") == 0)
            continue;
        std::pair<unsigned char *, unsigned int> hashed_component = hash_message((const char*)(t.c_str()));

        std::pair<BIGNUM *, BIGNUM *> utrapdoor = u->UserTD(hashed_component.first, hashed_component.second);

        vec->push_back(utrapdoor);

        free(hashed_component.first);
    }

    char * new_user_content_names = u->FullUserTD(vec);
    //printf("%s\n", new_user_content_names);

    // Atach ccnx: at the begining
    int length = strnlen(new_user_content_names, 1 + ((SECURITY_KEYSIZE/4 + 1)* vec->size() * 2) + 1 + 1);
    char *str2 = (char*) calloc(sizeof(char), length + 5 + 1);
    strcpy(str2, "ccnx:");
    strcat(str2, new_user_content_names);

    for(std::vector<std::pair<BIGNUM *, BIGNUM *> >::iterator it2 = vec->begin(); it2 != vec->end(); it2++)
    {
        std::pair<BIGNUM *, BIGNUM *> p = *it2;

        BN_clear_free(p.first);
        BN_clear_free(p.second);

    }
    delete vec;

    delete u;

    free(new_user_content_names);

    return str2;
}

// Encode and encrypt the content
extern "C" char * EncryptUserContentNoNetwork(const unsigned char * original_content, const unsigned int content_len, const BIGNUM * primeQ, const BIGNUM * userKey)
{
	//std::pair<unsigned char *, unsigned int> hashed_component = encodeIntoBase64(original_content, content_len);
	User *u = new User(512); // TODO: CHANGE

    u->setPrimeQ(primeQ);
    u->setUserKey(userKey);
	std::pair<BIGNUM *, BIGNUM *> p1 = u->ContentProviderEnc(original_content, content_len);

	char * first_part_encr = BN_bn2hex(p1.first);
	char * second_part_encr = BN_bn2hex(p1.second);

	/* pack both in a string: append both parts with a pipe */
    char *encrypted_content = (char*) calloc(sizeof(char), (SECURITY_KEYSIZE/4 + 1) * 2 + 1 + 1);
    strcpy(encrypted_content, first_part_encr);
	strcat(encrypted_content, "|");
    strcat(encrypted_content, second_part_encr);

	/* Free all the elements */
	BN_clear_free(p1.first);
	BN_clear_free(p1.second);
	OPENSSL_free(first_part_encr);
	OPENSSL_free(second_part_encr);

	return encrypted_content;
}

extern "C" char * libprotector_EncryptUserContentWithKeys(const unsigned char * original_content, const unsigned int content_len, const char * user_key, const char * prime_q)
{
    BIGNUM * primeQ = BN_new();
    BIGNUM * userKey = BN_new();

    BN_hex2bn(&primeQ, prime_q);
    BN_hex2bn(&userKey, user_key);

    User *u = new User(512); // TODO: CHANGE

    u->setPrimeQ(primeQ);
    u->setUserKey(userKey);
	std::pair<BIGNUM *, BIGNUM *> p1 = u->ContentProviderEnc(original_content, content_len);
	delete u;

	char * first_part_encr = BN_bn2hex(p1.first);
	char * second_part_encr = BN_bn2hex(p1.second);

	/* pack both in a string: append both parts with a pipe */
    char *encrypted_content = (char*) calloc(sizeof(char), (SECURITY_KEYSIZE/4 + 1) * 2 + 1 + 1);
    strcpy(encrypted_content, first_part_encr);
	strcat(encrypted_content, "|");
    strcat(encrypted_content, second_part_encr);

	/* Free all the elements */
	BN_clear_free(p1.first);
	BN_clear_free(p1.second);
	OPENSSL_free(first_part_encr);
	OPENSSL_free(second_part_encr);

	BN_clear_free(primeQ);
	BN_clear_free(userKey);

	return encrypted_content;
}

// Encode and encrypt the content
extern "C" char * libprotector_EncryptUserContent(const unsigned char * original_content, const unsigned int content_len)
{
	//std::pair<unsigned char *, unsigned int> hashed_component = encodeIntoBase64(original_content, content_len);

	char * res_getUserK = retrieveKeyFromServer("userK:0");
    char * res_getPrimeQ = retrieveKeyFromServer("primeQ");

	char * encrypted_content = libprotector_EncryptUserContentWithKeys(original_content, content_len, res_getUserK, res_getPrimeQ);

	free(res_getPrimeQ);
	free(res_getUserK);

	return encrypted_content;
}

extern "C" unsigned char * libprotector_ReDecryptContentWithKeys(const char * encrypted_content, const char * res_getUserK, const char * res_getPrimeQ)
{
    BIGNUM * primeQ = BN_new();
    BIGNUM * userKey = BN_new();

    BN_hex2bn(&primeQ, res_getPrimeQ);
    BN_hex2bn(&userKey, res_getUserK);

    User *u = new User(512); // TODO: CHANGE

    u->setPrimeQ(primeQ);
    u->setUserKey(userKey);

	// Convert char into BIGNUM
	std::vector<std::string> keywords;
    boost::split(keywords, encrypted_content, boost::is_any_of("|"));

	BIGNUM *b1 = BN_new();
	BIGNUM *b2 = BN_new();

	const char * c1 = keywords[0].c_str();
	const char * c2 = keywords[1].c_str();
	BN_hex2bn(&b1, c1);
	BN_hex2bn(&b2, c2);

	std::pair<BIGNUM *, BIGNUM *> utrapdoor = std::make_pair(b1, b2);

	BIGNUM * res = u->UserDec(utrapdoor);
	char * res_char = (char *) calloc(sizeof(char), SECURITY_KEYSIZE/4 + 2);
	BN_bn2bin(res, (unsigned char *) res_char);

	delete u;

	BN_clear_free(res);

	BN_clear_free(b1);
	BN_clear_free(b2);

	BN_clear_free(userKey);
	BN_clear_free(primeQ);

	return res_char;
}

extern "C" unsigned char * libprotector_ReDecryptContent(const char * encrypted_content)
{
	char * res_getUserK = retrieveKeyFromServer("userK:0");
    char * res_getPrimeQ = retrieveKeyFromServer("primeQ");

	char * res_char = libprotector_ReDecryptContentWithKeys(encrypted_content, res_getUserK, res_getPrimeQ);

	free(res_getPrimeQ);
	free(res_getUserK);

	return (unsigned char *) res_char;//res_final;
}

extern "C" unsigned char * libprotector_ReDecryptAndSplitContent(const char * encrypted_content)
{
	char * res_getUserK = retrieveKeyFromServer("userK:0");
    char * res_getPrimeQ = retrieveKeyFromServer("primeQ");

	BIGNUM * primeQ = BN_new();
    BIGNUM * userKey = BN_new();

    BN_hex2bn(&primeQ, res_getPrimeQ);
    BN_hex2bn(&userKey, res_getUserK);

    User *u = new User(512); // TODO: CHANGE

    u->setPrimeQ(primeQ);
    u->setUserKey(userKey);

	// Convert char into BIGNUM
	std::vector<std::string> keywords;
    boost::split(keywords, encrypted_content, boost::is_any_of("|"));

    char * b64_message = (char *) calloc(sizeof(char), (SECURITY_KEYSIZE/4 + 2)*keywords.size()/2);
    b64_message[0] = '\0';
    int offset = 0;

    //printf("There are # %d keywords (%d)\n", keywords.size(), strlen(encrypted_content));

    for(unsigned int i=0; i<keywords.size(); i+=2)
    {
	    BIGNUM *b1 = BN_new();
	    BIGNUM *b2 = BN_new();

	    const char * c1 = keywords[i].c_str();
	    const char * c2 = keywords[i+1].c_str();
	    BN_hex2bn(&b1, c1);
	    BN_hex2bn(&b2, c2);

	    std::pair<BIGNUM *, BIGNUM *> utrapdoor = std::make_pair(b1, b2);

	    BIGNUM * res = u->UserDec(utrapdoor);
        //TODO: HERE IS THE PROBLEM!!!
	    char * aux = (char *) calloc(sizeof(char), SECURITY_KEYSIZE/4 + 2);
	    int copied_bytes = BN_bn2bin(res, (unsigned char *) aux);
	    BN_clear_free(res);
        memcpy(b64_message+offset, aux, copied_bytes);
        //printf("COPIED: %d\n", copied_bytes);
        offset += copied_bytes;

        free(aux);

        BN_clear_free(b1);
        BN_clear_free(b2);
    }
    b64_message[offset+1] = '\0';

    delete u;

    BN_clear_free(primeQ);
    BN_clear_free(userKey);

    free(res_getPrimeQ);
    free(res_getUserK);

    //printf("The base64 message is %s (%d vs %d)\n", b64_message, strlen(b64_message),  offset);

	return (unsigned char*) b64_message;
}

/* libprotector C interface */
extern "C" User * libprotector_User_new(unsigned int keysize)
{
    return new User(keysize);
}

extern "C" void libprotector_User_setUserKey(User * u, const char * user_key)
{
    BIGNUM * aux = BN_new();
    BN_hex2bn(&aux, user_key);
    reinterpret_cast<User *>(u)->setUserKey(aux);
    BN_clear_free(aux);
};

extern "C" void libprotector_User_setNumberG(User * u, const char * number_g)
{
    BIGNUM * aux = BN_new();
    BN_hex2bn(&aux, number_g);
    reinterpret_cast<User *>(u)->setNumberG(aux);
    BN_clear_free(aux);
};

extern "C" void libprotector_User_setPrimeQ(User * u, const char * prime_q)
{
    BIGNUM * aux = BN_new();
    BN_hex2bn(&aux, prime_q);
    reinterpret_cast<User *>(u)->setPrimeQ(aux);
    BN_clear_free(aux);
};

extern "C" void libprotector_User_setSalt(User * u, const char * salt_key)
{
    BIGNUM * aux = BN_new();
    BN_hex2bn(&aux, salt_key);
    reinterpret_cast<User *>(u)->setSalt(aux);
    BN_clear_free(aux);
};

extern "C" void libprotector_User_UserTD(User * u, const unsigned char * component, unsigned int size_component, char *& p1, char *& p2)
{
    std::pair<BIGNUM *, BIGNUM *> t = reinterpret_cast<User *>(u)->UserTD(component, size_component);
    unsigned int keysize = reinterpret_cast<User *>(u)->keysize;
    p1 = (char *) calloc(sizeof(char), ((keysize/4)) );
    char * aux_component = BN_bn2hex(t.first);
    memcpy((void *) (p1), (void *) aux_component, (keysize/4));

    p2 = (char *) calloc(sizeof(char), ((keysize/4)) );
    char * aux_component2 = BN_bn2hex(t.second);
    memcpy((void *) (p2), (void *) aux_component2, (keysize/4));

    free(aux_component);
    free(aux_component2);
}

extern "C" void libprotector_User_ContentTD(User * u, const unsigned char * block, const unsigned int block_size, char *& p1, char *& p2)
{
    std::pair<BIGNUM *, BIGNUM *> t = reinterpret_cast<User *>(u)->ContentProviderEnc(block, block_size);

    unsigned int keysize = reinterpret_cast<User *>(u)->keysize;
    p1 = (char *) calloc(sizeof(char), ((keysize/4)) );
    char * aux_component = BN_bn2hex(t.first);
    memcpy((void *) (p1), (void *) aux_component, (keysize/4));

    p2 = (char *) calloc(sizeof(char), ((keysize/4)) );
    char * aux_component2 = BN_bn2hex(t.second);
    memcpy((void *) (p2), (void *) aux_component2, (keysize/4));

    free(aux_component);
    free(aux_component2);
}

extern "C" void libprotector_User_ClientDec(User * u, const char * component1, const char * component2, char *& p1)
{
    unsigned int keysize = reinterpret_cast<User *>(u)->keysize;

    BIGNUM * aux = BN_new();
    BN_hex2bn(&aux, component1);

    BIGNUM * aux2 = BN_new();
    BN_hex2bn(&aux2, component2);

    std::pair<BIGNUM *, BIGNUM *> encrypted_payload = std::make_pair(aux, aux2);

    BIGNUM * decrypted_value = reinterpret_cast<User *>(u)->UserDec(encrypted_payload);

    p1 = (char *) calloc(sizeof(char), ((keysize/4)) );
    char * aux_component = BN_bn2hex(decrypted_value);

    BN_bn2bin(decrypted_value, (unsigned char *) p1);

    free(aux_component);
}

ext: cc

max_stars_repo_name: mesarpe/libprotector

max_issues_repo_name: mesarpe/libprotector

max_forks_repo_name: mesarpe/libprotector


text:
// Author(s): Wieger Wesselink
// Copyright: see the accompanying file COPYING or copy at
// https://svn.win.tue.nl/trac/MCRL2/browser/trunk/COPYING
//
// 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)
//
/// \file action_test.cpp
/// \brief Add your file description here.

#include "mcrl2/core/print.h"
#include "mcrl2/data/function_symbol.h"
#include "mcrl2/process/process_expression.h"
#include <boost/test/minimal.hpp>
#include <iostream>
#include <set>
#include <string>

using namespace mcrl2;

int test_main(int argc, char** argv)
{
  data::sort_expression X(data::basic_sort("X"));
  data::sort_expression Y(data::basic_sort("Y"));

  process::action_label aX(core::identifier_string("a"),  { X });
  process::action_label aY(core::identifier_string("a"),  { Y });
  process::action_label aXY(core::identifier_string("a"), { X, Y });
  process::action_label bX(core::identifier_string("b"),  { X });

  data::data_expression x_X = data::variable("x", X);
  data::data_expression y_X = data::variable("y", X);
  data::data_expression y_Y = data::variable("y", Y);

  process::action a1(aX, { x_X });
  process::action a2(aX, { y_X });
  process::action a3(bX, { x_X });

  BOOST_CHECK(process::equal_signatures(a1, a1));
  BOOST_CHECK(process::equal_signatures(a1, a2));
  BOOST_CHECK(!process::equal_signatures(a1, a3));

  process::action_label_list l;
  l.push_front(aX);
  BOOST_CHECK(l.front()==aX);
  l.push_front(aY);
  l.push_front(aY);
  BOOST_CHECK(l.front()==aY);
  process::action_label_list l1=l;
  BOOST_CHECK(l1.front()==aY);
  BOOST_CHECK(l.front()==aY);
  process::action_label a=l.front();
  BOOST_CHECK(a==aY);

  return 0;
}


text:
/*
 *            Copyright 2009-2020 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      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.
 *
 */

// Third party includes
#include <boost/interprocess/sync/file_lock.hpp>

// VOTCA includes
#include <votca/tools/filesystem.h>

// Local VOTCA includes
#include "votca/xtp/statesaver.h"
#include "votca/xtp/topology.h"

namespace votca {
namespace xtp {

std::vector<Index> StateSaver::getFrames() const {
  CheckpointFile cpf(_hdf5file, CheckpointAccessLevel::READ);
  CheckpointReader r = cpf.getReader();
  std::vector<Index> frames = std::vector<Index>{};
  try {
    r(frames, "frames");
  } catch (std::runtime_error&) {
    ;
  }
  return frames;
}
void StateSaver::WriteFrame(const Topology& top) {
  if (!tools::filesystem::FileExists(_hdf5file)) {
    std::cout << "Creating statefile " << _hdf5file << std::endl;
    CheckpointFile cpf(_hdf5file, CheckpointAccessLevel::CREATE);
  }
  boost::interprocess::file_lock flock(_hdf5file.c_str());
  flock.lock();
  if (!TopStepisinFrames(top.getStep())) {
    std::vector<Index> frames = this->getFrames();
    frames.push_back(top.getStep());
    CheckpointFile cpf(_hdf5file, CheckpointAccessLevel::MODIFY);
    CheckpointWriter w = cpf.getWriter();
    w(frames, "frames");
    std::cout << "Frame with id " << top.getStep() << " was not in statefile "
              << _hdf5file << " ,adding it now." << std::endl;
  }

  CheckpointFile cpf(_hdf5file, CheckpointAccessLevel::MODIFY);
  CheckpointWriter w = cpf.getWriter("/frame_" + std::to_string(top.getStep()));
  top.WriteToCpt(w);

  flock.unlock();

  std::cout << "Wrote MD topology (step = " << top.getStep()
            << ", time = " << top.getTime() << ") to " << _hdf5file
            << std::endl;
  std::cout << "... ";

  std::cout << ". " << std::endl;
  return;
}

Topology StateSaver::ReadFrame(Index frameid) const {
  if (!tools::filesystem::FileExists(_hdf5file)) {
    throw std::runtime_error("Statefile " + _hdf5file + " does not exist.");
  }
  std::cout << "Import MD Topology (i.e. frame " << frameid << ")"
            << " from " << _hdf5file << std::endl;
  std::cout << "...";
  boost::interprocess::file_lock flock(_hdf5file.c_str());
  flock.lock();
  if (!TopStepisinFrames(frameid)) {
    throw std::runtime_error("Frame with id " + std::to_string(frameid) +
                             " is not in statefile.");
  }

  CheckpointFile cpf(_hdf5file, CheckpointAccessLevel::READ);
  CheckpointReader r = cpf.getReader("/frame_" + std::to_string(frameid));
  Topology top;
  top.ReadFromCpt(r);
  flock.unlock();
  std::cout << ". " << std::endl;
  return top;
}

bool StateSaver::TopStepisinFrames(Index frameid) const {
  std::vector<Index> frames = this->getFrames();
  return std::find(frames.begin(), frames.end(), frameid) != frames.end();
}

}  // namespace xtp
}  // namespace votca

max_stars_repo_name: fossabot/xtp

max_issues_repo_name: fossabot/xtp

max_forks_repo_name: fossabot/xtp


text:
// Copyright 2011-2021 Google LLC
//
// 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.

#include "third_party/zynamics/binexport/reader/graph_utility.h"

#include <utility>
#include <vector>

#include <boost/graph/adjacency_list.hpp>  // NOLINT

#include "gmock/gmock.h"
#include "gtest/gtest.h"

namespace security::binexport {

class GraphUtilityTest : public testing::Test {
 protected:
  const std::vector<Address> addresses_{0,  1,  2,  3,  4,  5,  6,  7,  8,  9,
                                        10, 11, 12, 13, 14, 15, 16, 16, 18, 19};

  using edge = std::pair<uint32_t, uint32_t>;
  const std::vector<edge> edges_{
      {0, 1},   {1, 2},   {2, 3},   {3, 4},   {4, 5},   {5, 6},   {5, 7},
      {5, 8},   {5, 9},   {6, 10},  {7, 10},  {8, 10},  {9, 10},  {10, 11},
      {11, 12}, {12, 13}, {13, 14}, {14, 15}, {14, 15}, {15, 16}, {16, 17},
      {17, 18}, {18, 19}, {1, 18},  {17, 2}};

  //    [0]->[1]->[2]->[3]->[4]->[5]->[6]->[10]->[11]->[12]->[13]->[14]->[15]
  //          ^    ^              +-->[7]----^                       +----^|
  //          |    |              +-->[8]----^                             |
  //          |    |              +-->[9]----^                             |
  //          |    +---------------------------------------------+         |
  //          +-------------------------------------------+      |         |
  //                                               [19]<-[18]<-[17]<-[16]<-+
  // [20] <- alone to test if the graph hashing works if a vertex is not
  // connected by any edge.

  struct VertexInfo {
    Address address;
  };

 public:
  struct EdgeProperty {};

  using Graph =
      boost::compressed_sparse_row_graph<boost::bidirectionalS, VertexInfo,
                                         EdgeProperty, boost::no_property,
                                         uint32_t, uint32_t>;

  Graph graph_;

  using Vertex = boost::graph_traits<Graph>::vertex_descriptor;
  using Edge = boost::graph_traits<Graph>::edge_descriptor;
  using EdgeIterator = boost::graph_traits<Graph>::edge_iterator;
  using OutEdgeIterator = boost::graph_traits<Graph>::out_edge_iterator;
  using InEdgeIterator = boost::graph_traits<Graph>::in_edge_iterator;
  using AdjacencyIterator = boost::graph_traits<Graph>::adjacency_iterator;

  using VertexIndexMap =
      boost::property_map<Graph, boost::vertex_index_t>::type;

  using UndirectedGraph =
      boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
                            VertexInfo, EdgeProperty>;

 protected:
  void SetUp() override {
    graph_ = Graph(boost::edges_are_unsorted_multi_pass, edges_.begin(),
                   edges_.end(), addresses_.size());
    using VertexIterator = boost::graph_traits<Graph>::vertex_iterator;
    VertexIterator vertices_it, vertices_end;
    int j = 0;
    for (boost::tie(vertices_it, vertices_end) = boost::vertices(graph_);
         vertices_it != vertices_end; ++vertices_it, ++j) {
      graph_[j].address = addresses_.at(j);
    }
  }
};

TEST_F(GraphUtilityTest, GetEdgeVector) {
  const std::pair<Edge, bool> edge_pair(
      boost::edge(Vertex(5), Vertex(6), graph_));
  const EdgeDegrees edge_vector =
      GetEdgeDegrees<GraphUtilityTest>(graph_, edge_pair.first);
  EXPECT_EQ(edge_vector.source_in_degree, 1);
  EXPECT_EQ(edge_vector.source_out_degree, 4);
  EXPECT_EQ(edge_vector.target_in_degree, 1);
  EXPECT_EQ(edge_vector.target_out_degree, 1);
}

}  // namespace security::binexport

max_stars_repo_name: clayne/binexport

max_issues_repo_name: gamozolabs/binexport

max_forks_repo_name: gamozolabs/binexport


text:
#include <pcl/io/pcd_io.h>
#include <pcl/point_cloud.h>
#include <pcl/correspondence.h>

#include <pcl/features/normal_3d_omp.h>   //法线
#include <pcl/features/3dsc.h>     //3dsc
#include <pcl/features/board.h>

//#include <pcl/filters/uniform_sampling.h>
#include <pcl/keypoints/uniform_sampling.h>
#include <pcl/keypoints/iss_3d.h>

#include <pcl/recognition/cg/hough_3d.h>
#include <pcl/recognition/cg/geometric_consistency.h>

#include <pcl/visualization/pcl_visualizer.h>

#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/kdtree/impl/kdtree_flann.hpp>

#include <pcl/common/transforms.h>
#include <pcl/console/parse.h>

// ensaio de acerto de resolução
#include <iostream>

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <boost/lexical_cast.hpp>

#include <pcl/point_types.h>
#include <pcl/filters/voxel_grid.h>

#include <sstream>

// tempo
#include <ctime>
//conversão integer to string
#include <boost/lexical_cast.hpp>

typedef pcl::PointXYZ PointType;
typedef pcl::Normal NormalType;
typedef pcl::ReferenceFrame RFType;

typedef pcl::PointWithViewpoint PointTypeViewPoint;

typedef pcl::ShapeContext1980 DescriptorType;

std::string model_filename_;
std::string scene_filename_;

//Algorithm params
bool show_keypoints_ (false);
bool show_correspondences_ (false);
bool show_normals_(false);
bool use_cloud_resolution_ (false);
bool use_hough_ (true);

//Keypoints parameters
bool use_uniformsampling(true);// default
bool use_iss(false);
bool use_harris3d(false);

//Features parameters
bool use_shot(true);//default
bool use_pfh(false);
bool use_usd(false);

float model_ss_ (0.01f); //amostragem do modelo em milimetros
float scene_ss_ (0.01f); //amostragem da cena em milimetros

float norm_rad_ (0.1f); //raio da esfera para o cálculo das normmais

float rf_rad_ (0.05f);
float descr_rad_ (0.04f); //raio da esfera para o cálculo dos descritores

// Resolution of the consensus set used to cluster correspondences together,
	// in metric units. Default is 1.0.
float cg_size_ (0.08f); //tamanho minimo do cluster dos objectos
// Minimum cluster size. Default is 3 (as at least 3 correspondences
	// are needed to compute the 6 DoF pose).
float cg_thresh_ (3000.0f); //

float Lim_dist_ (200000.0f); //Distancia minima entre descritores para que se considere válida a correspondencia

// Cronometro
std::clock_t start;
std::clock_t end;

double diffclock(std::clock_t clock1,std::clock_t clock2)
{
  double diffticks = clock2 - clock1;
  double diffms = diffticks / (CLOCKS_PER_SEC / 1000 );
  return diffms;
}

void
showHelp (char *filename)
{
  std::cout << std::endl;
  std::cout << "***************************************************************************" << std::endl;
  std::cout << "*                            Usage Guide                                  *" << std::endl;
  std::cout << "***************************************************************************" << std::endl << std::endl;
  std::cout << "Usage: " << filename << " model_filename.pcd scene_filename.pcd [Options]" << std::endl << std::endl;
  std::cout << "Options:" << std::endl;
  std::cout << "     -h:                     Show this help." << std::endl;
  std::cout << "     -k:                     Show used keypoints." << std::endl;
  std::cout << "     -c:                     Show used correspondences." << std::endl;
  std::cout << "     -n:                     Show used Normals." << std::endl;
  std::cout << "     -r:                     Compute the model cloud resolution and multiply" << std::endl;
  std::cout << "                             each radius given by that value." << std::endl;
  std::cout << "     --Keypoints: keypoints downsample algorithm used (default UniformSampling)." << std::endl;
  std::cout << "                  ISS (Intrinsic Shape Signatures)" << std::endl;
  //std::cout << "                  Harris3D                        " << std::endl;
  std::cout << "     --Features:  descritor features algorithm used (default SHOT)." << std::endl;
  //std::cout << "                  PFH (Point Features Histogram)" << std::endl;
  //std::cout << "                  ISS (Intrinsic Shape Signatures)" << std::endl;
  //std::cout << "                  USD (Unique Shape Context)" << std::endl;
  std::cout << "     --algorithm (Hough|GC): Clustering algorithm used (default Hough)." << std::endl;
  std::cout << "     --model_ss val:         Model uniform sampling radius (default 0.01)" << std::endl;
  std::cout << "     --scene_ss val:         Scene uniform sampling radius (default 0.03)" << std::endl;
  std::cout << "     --rf_rad val:           Reference frame radius (default 0.015)" << std::endl;
  std::cout << "     --descr_rad val:        Descriptor radius (default 0.02)" << std::endl;
  std::cout << "     --cg_size val:          Cluster size (default 0.01)" << std::endl;
  std::cout << "     --cg_thresh val:        Clustering threshold (default 5)" << std::endl << std::endl;
}

void
parseCommandLine (int argc, char *argv[])
{
  //Show help
  if (pcl::console::find_switch (argc, argv, "-h"))
  {
    showHelp (argv[0]);
    exit (0);
  }

  //Model & scene filenames
  std::vector<int> filenames;
  filenames = pcl::console::parse_file_extension_argument (argc, argv, ".pcd");
  if (filenames.size () != 2)
  {
    std::cout << "Filenames missing.\n";
    showHelp (argv[0]);
    exit (-1);
  }

  model_filename_ = argv[filenames[0]];
  scene_filename_ = argv[filenames[1]];

  //Program visualization behavior
  if (pcl::console::find_switch (argc, argv, "-k"))
  {
    show_keypoints_ = true;
  }
  if (pcl::console::find_switch (argc, argv, "-c"))
  {
    show_correspondences_ = true;
  }
  if (pcl::console::find_switch (argc, argv, "-r"))
  {
    use_cloud_resolution_ = true;
  }

  if (pcl::console::find_switch (argc, argv, "-n"))
  {
    show_normals_ = true;
  }

  // selection of algorithm for clustering
  std::string used_algorithm;
  if (pcl::console::parse_argument (argc, argv, "--algorithm", used_algorithm) != -1)
  {
    if (used_algorithm.compare ("Hough") == 0)  //霍夫变化
    {
      use_hough_ = true;
    }
    else if (used_algorithm.compare ("GC") == 0)   //geometry consietency几何一致性
    {
      use_hough_ = false;
    }
    else
    {
      std::cout << "Wrong algorithm name.\n";
      showHelp (argv[0]);
      exit (-1);
    }
  }

  // Selection of algorithm for Keypoint extration
  std::string used_keypoints;
  if (pcl::console::parse_argument (argc, argv, "--keypoint", used_keypoints) != -1)
  {
    if (used_keypoints.compare ("UniformSample") == 0)
    {

    }
    else if (used_keypoints.compare ("ISS") == 0)
    {
      use_uniformsampling = false;
      use_iss = true;
    }
    else if (used_keypoints.compare ("Harris3D") == 0)
    {
      use_uniformsampling = false;
      use_harris3d = true;
    }
    else
    {
      std::cout << "Wrong Keypoint name.\n";
      showHelp (argv[0]);
      exit (-1);
    }
  }

  //General parameters
  pcl::console::parse_argument (argc, argv, "--model_ss", model_ss_);
  pcl::console::parse_argument (argc, argv, "--scene_ss", scene_ss_);
  pcl::console::parse_argument (argc, argv, "--rf_rad", rf_rad_);
  pcl::console::parse_argument (argc, argv, "--descr_rad", descr_rad_);
  pcl::console::parse_argument (argc, argv, "--cg_size", cg_size_);
  pcl::console::parse_argument (argc, argv, "--cg_thresh", cg_thresh_);
}

//***********************************************************************************************
//                                  Função Cálculo resolução nuvem
//***********************************************************************************************

//计算点云的分辨率
double
computeCloudResolution (const pcl::PointCloud<PointType>::ConstPtr &cloud)
{
  double res = 0.0;
  int n_points = 0;
  int nres;

  std::vector<int> indices (2);
  std::vector<float> sqr_distances (2);
  pcl::search::KdTree<PointType> tree;
  tree.setInputCloud (cloud);

  for (size_t i = 0; i < cloud->size (); ++i)
  {
    if (! pcl_isfinite ((*cloud)[i].x))//判断点的x的值是否位无效点也就是NAN点   如果是无效点就不会计算
    {
      continue;
    }
    //Considering the second neighbor since the first is the point itself.
  //int pcl::search::Search<PointT>::nearSearch();是计算给定点的k个邻域的点
   //i 代表索引点，2是寻找周围的点树（其中第一个点是它本身）indices寻找到的周围点的索引sqr_distances与相邻点的平方距离
    nres = tree.nearestKSearch (i, 2, indices, sqr_distances);  //返回寻找到的点数
    if (nres == 2)
    {
      res += sqrt (sqr_distances[1]);
      ++n_points;
    }
  }
  if (n_points != 0)
  {
    res /= n_points;
  }
  return res;
}

int
main (int argc, char *argv[])
{
  parseCommandLine (argc, argv);

  pcl::PointCloud<PointTypeViewPoint>::Ptr modelviewpoint (new pcl::PointCloud<PointTypeViewPoint> ());

  pcl::PointCloud<PointType>::Ptr model (new pcl::PointCloud<PointType> ());
  pcl::PointCloud<PointType>::Ptr scene (new pcl::PointCloud<PointType> ());

  pcl::PointCloud<PointType>::Ptr model_keypoints (new pcl::PointCloud<PointType> ());
  pcl::PointCloud<PointType>::Ptr scene_keypoints (new pcl::PointCloud<PointType> ());

  pcl::PointCloud<NormalType>::Ptr model_normals (new pcl::PointCloud<NormalType> ());
  pcl::PointCloud<NormalType>::Ptr scene_normals (new pcl::PointCloud<NormalType> ());

  pcl::PointCloud<NormalType>::Ptr model_keypoints_normals (new pcl::PointCloud<NormalType> ());
  pcl::PointCloud<NormalType>::Ptr scene_keypoints_normals (new pcl::PointCloud<NormalType> ());

  pcl::PointCloud<DescriptorType>::Ptr model_descriptors (new pcl::PointCloud<DescriptorType> ());
  pcl::PointCloud<DescriptorType>::Ptr scene_descriptors (new pcl::PointCloud<DescriptorType> ());

  //
  //  Load clouds
  //
  if (pcl::io::loadPCDFile (model_filename_, *modelviewpoint) < 0)
  {
    std::cout << "Error loading model cloud." << std::endl;
    showHelp (argv[0]);
    return (-1);
  }
  if (pcl::io::loadPCDFile (scene_filename_, *scene) < 0)
  {
    std::cout << "Error loading scene cloud." << std::endl;
    showHelp (argv[1]);
    return (-1);
  }

  // Recolha dos view points para calculo das normais   //设置视点
  float vpx = modelviewpoint->points[0].vp_x;
  float vpy = modelviewpoint->points[0].vp_y;
  float vpz = modelviewpoint->points[0].vp_z;

  // Copia da nuvem para tipologia sem view points
  pcl::copyPointCloud<PointTypeViewPoint,PointType>(*modelviewpoint, *model);

  // see the resolution of the clouds
  start = std::clock();                                    // Debug time

  double res_scene = computeCloudResolution(scene);
  double res_model = computeCloudResolution(model);

  end = std::clock();                                       // Debug time
  std::cout << "Time computing clouds resolution:  " << diffclock(start, end) << " ms" << endl << endl;

  std::cout << "Scene resolution:  " << res_scene << std::endl;
  std::cout << "Model resolution:  " << res_model << std::endl;

  //
  //  Set up resolution invariance
  //
  if (use_cloud_resolution_)
  {
    float resolution = static_cast<float> (computeCloudResolution (model));
    if (resolution != 0.0f)
    {
      model_ss_   *= resolution;
      scene_ss_   *= resolution;
      rf_rad_     *= resolution;
      descr_rad_  *= resolution;
      cg_size_    *= resolution;
    }

    std::cout << "Model resolution:       " << resolution << std::endl;
    std::cout << "Model sampling size:    " << model_ss_ << std::endl;
    std::cout << "Scene sampling size:    " << scene_ss_ << std::endl;
    std::cout << "LRF support radius:     " << rf_rad_ << std::endl;
    std::cout << "SHOT descriptor radius: " << descr_rad_ << std::endl;
    std::cout << "Clustering bin size:    " << cg_size_ << std::endl << std::endl;
  }

  //
  //  Compute Normals
  //
  /*pcl::NormalEstimationOMP<PointType, NormalType> norm_est;
  norm_est.setKSearch (100);

  norm_est.setInputCloud (scene);
  norm_est.compute (*scene_normals);

  norm_est.setInputCloud (model);
  norm_est.setViewPoint(vpx,vpy,vpz);
  norm_est.compute (*model_normals);*/

//
//  Downsample Clouds to Extract keypoints
//
        if(use_uniformsampling)
        {
        pcl::UniformSampling<PointType> uniform_sampling;
        uniform_sampling.setInputCloud (model);
        uniform_sampling.setRadiusSearch (model_ss_);
        //uniform_sampling.filter (*model_keypoints);
        pcl::PointCloud<int> keypointIndices1;
        uniform_sampling.compute(keypointIndices1);

        pcl::copyPointCloud(*model, keypointIndices1.points, *model_keypoints);
      //  pcl::copyPointCloud(*model_normals, keypointIndices1.points, *model_keypoints_normals);

        std::cout << "Model total points: " << model->size () << "; Selected Keypoints: " << model_keypoints->size () << std::endl;

        uniform_sampling.setInputCloud (scene);
        uniform_sampling.setRadiusSearch (scene_ss_);
        //uniform_sampling.filter (*scene_keypoints);
        pcl::PointCloud<int> keypointIndices2;
        uniform_sampling.compute(keypointIndices2);

        pcl::copyPointCloud(*scene, keypointIndices2.points, *scene_keypoints);
      //  pcl::copyPointCloud(*scene_normals, keypointIndices2.points, *scene_keypoints_normals);

        std::cout << "Scene total points: " << scene->size () << "; Selected Keypoints: " << scene_keypoints->size () << std::endl;
        }

        //
        //  Compute Normals
        //
        pcl::NormalEstimation<PointType, NormalType> norm_est;
        norm_est.setRadiusSearch(norm_rad_);
//norm_est.setKSearch (100);
        norm_est.setInputCloud (scene_keypoints);
        norm_est.compute (*scene_keypoints_normals);
std::cout << scene_keypoints_normals->size() << std::endl;//Debug
        norm_est.setInputCloud (model_keypoints);
        norm_est.setViewPoint(vpx,vpy,vpz);
        norm_est.compute (*model_keypoints_normals);
std::cout << model_keypoints_normals->size() << std::endl;//Debug

std::cout << "PASSOU" << std::endl;
  //
  //  Downsample Clouds to Extract keypoints ////////////////////////////////////////////////
  //
  /*
        if(use_iss)
        {
        // ISS keypoint detector model.
        pcl::ISSKeypoint3D<pcl::PointXYZ, pcl::PointXYZ> detector;
        detector.setInputCloud(model);
        pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree_model(new pcl::search::KdTree<pcl::PointXYZ>);
        detector.setSearchMethod(kdtree_model);
        double resolution_model = computeCloudResolution(model);
        // Set the radius of the spherical neighborhood used to compute the scatter matrix.
        detector.setSalientRadius(6 * resolution_model);
        // Set the radius for the application of the non maxima supression algorithm.
        detector.setNonMaxRadius(4 * resolution_model);
        // Set the minimum number of neighbors that has to be found while applying the non maxima suppression algorithm.
        detector.setMinNeighbors(5);
        // Set the upper bound on the ratio between the second and the first eigenvalue.
        detector.setThreshold21(0.975);
        // Set the upper bound on the ratio between the third and the second eigenvalue.
        detector.setThreshold32(0.975);
        // Set the number of prpcessing threads to use. 0 sets it to automatic.
        detector.setNumberOfThreads(4);

        detector.setNormalRadius (4 * resolution_model);
        detector.setBorderRadius (1 * resolution_model);

        detector.compute(*model_keypoints);
        std::cout << "model_keypoints: " << model_keypoints->size () << std::endl;//debug

        // ISS keypoint detector scene.
        //pcl::ISSKeypoint3D<pcl::PointXYZ, pcl::PointXYZ> detector;
        detector.setInputCloud(scene);
        pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree_scene(new pcl::search::KdTree<pcl::PointXYZ>);
        detector.setSearchMethod(kdtree_scene);
        double resolution_scene = computeCloudResolution(scene);
        // Set the radius of the spherical neighborhood used to compute the scatter matrix.
        detector.setSalientRadius(6 * resolution_scene);
        // Set the radius for the application of the non maxima supression algorithm.
        detector.setNonMaxRadius(4 * resolution_scene);
        // Set the minimum number of neighbors that has to be found while applying the non maxima suppression algorithm.
        detector.setMinNeighbors(5);
        // Set the upper bound on the ratio between the second and the first eigenvalue.
        detector.setThreshold21(0.975);
        // Set the upper bound on the ratio between the third and the second eigenvalue.
        detector.setThreshold32(0.975);
        // Set the number of prpcessing threads to use. 0 sets it to automatic.
        detector.setNumberOfThreads(4);

        detector.setNormalRadius (4 * resolution_scene);
        detector.setBorderRadius (1 * resolution_scene);

        detector.compute(*scene_keypoints);
        std::cout << "scene_keypoints: " << scene_keypoints->size () << std::endl;//debug
        }
        */
  //
  //  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //

  //
  //  Compute Descriptor for keypoints SHOT
  //
/*
        pcl::SHOTEstimationOMP<PointType, NormalType, DescriptorType> descr_est;
        descr_est.setRadiusSearch (descr_rad_);

        descr_est.setInputCloud (model_keypoints);
        descr_est.setInputNormals (model_keypoints_normals);
        descr_est.setSearchSurface (model_keypoints);
        descr_est.compute (*model_descriptors);

        descr_est.setInputCloud (scene_keypoints);
        descr_est.setInputNormals (scene_keypoints_normals);
        descr_est.setSearchSurface (scene_keypoints);
        descr_est.compute (*scene_descriptors);
*/

  //
  //  Compute Descriptor for keypoints 3DSC
  //
        pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree(new pcl::search::KdTree<pcl::PointXYZ>);

        pcl::ShapeContext3DEstimation<PointType, NormalType, DescriptorType> descr_est;
        // Search radius, to look for neighbors. It will also be the radius of the support sphere.
        descr_est.setRadiusSearch (descr_rad_);

        descr_est.setInputCloud (model_keypoints);
        descr_est.setInputNormals (model_keypoints_normals);
        descr_est.setSearchMethod(kdtree);
        // The minimal radius value for the search sphere, to avoid being too sensitive
	      // in bins close to the center of the sphere.
	      descr_est.setMinimalRadius(0.05 / 10.0);
        // Radius used to compute the local point density for the neighbors
        // (the density is the number of points within that radius).
        descr_est.setPointDensityRadius(0.05 / 5.0);
        descr_est.compute (*model_descriptors);

        descr_est.setInputCloud (scene_keypoints);
        descr_est.setInputNormals (scene_keypoints_normals);
        descr_est.setSearchMethod(kdtree);
        // The minimal radius value for the search sphere, to avoid being too sensitive
	      // in bins close to the center of the sphere.
	      descr_est.setMinimalRadius(0.05 / 10.0);
        // Radius used to compute the local point density for the neighbors
        // (the density is the number of points within that radius).
        descr_est.setPointDensityRadius(0.05 / 5.0);
        descr_est.compute (*scene_descriptors);

  //
  //  Compute Descriptor for keypoints FPFH
  //
        //if(use_pfh)
        //{
  /*      pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree(new pcl::search::KdTree<pcl::PointXYZ>);

        pcl::FPFHEstimation<PointType, NormalType, DescriptorType> descr_est;
        descr_est.setRadiusSearch (descr_rad_);

        descr_est.setInputCloud (model_keypoints);
        descr_est.setInputNormals (model_keypoints_normals);
        descr_est.setSearchMethod(kdtree);
        // Search radius, to look for neighbors. Note: the value given here has to be
        // larger than the radius used to estimate the normals.
        descr_est.compute (*model_descriptors);

        descr_est.setInputCloud (scene_keypoints);
        descr_est.setInputNormals (scene_keypoints_normals);
        descr_est.setSearchMethod(kdtree);
        // Search radius, to look for neighbors. Note: the value given here has to be
        // larger than the radius used to estimate the normals.
        descr_est.compute (*scene_descriptors);
*/
        //Debug----------------------------------------------------------------------------------------
      /*  for (size_t i = 0; i < model_descriptors->size (); ++i)
        {

              std::cout << " desccriptors " << i << ":" << model_descriptors->at (i).histogram[0] <<
                                                   " " << model_descriptors->at (i).histogram[1] <<
                                                   " " << model_descriptors->at (i).histogram[2] <<
                                                   " " << model_descriptors->at (i).histogram[3] <<
                                                   " " << model_descriptors->at (i).histogram[4] <<
                                                   " " << model_descriptors->at (i).histogram[5] <<
                                                   " " << model_descriptors->at (i).histogram[6] <<
                                                   " " << model_descriptors->at (i).histogram[7] <<
                                                   " " << model_descriptors->at (i).histogram[8] <<
                                                   " " << model_descriptors->at (i).histogram[9] <<
                                                   " " << model_descriptors->at (i).histogram[10] <<
                                                   " " << model_descriptors->at (i).histogram[11] <<
                                                   " " << model_descriptors->at (i).histogram[12] <<
                                                   " " << model_descriptors->at (i).histogram[13] <<
                                                   " " << model_descriptors->at (i).histogram[14] <<
                                                   " " << model_descriptors->at (i).histogram[15] <<
                                                   " " << model_descriptors->at (i).histogram[16] <<
                                                   " " << model_descriptors->at (i).histogram[17] <<
                                                   " " << model_descriptors->at (i).histogram[18] <<
                                                   " " << model_descriptors->at (i).histogram[19] <<
                                                   " " << model_descriptors->at (i).histogram[20] <<
                                                   " " << model_descriptors->at (i).histogram[21] <<
                                                   " " << model_descriptors->at (i).histogram[22] <<
                                                   " " << model_descriptors->at (i).histogram[23] <<
                                                   " " << model_descriptors->at (i).histogram[24] <<
                                                   " " << model_descriptors->at (i).histogram[25] <<
                                                   " " << model_descriptors->at (i).histogram[26] <<
                                                   " " << model_descriptors->at (i).histogram[27] <<
                                                   " " << model_descriptors->at (i).histogram[28] <<
                                                   " " << model_descriptors->at (i).histogram[29] <<
                                                   " " << model_descriptors->at (i).histogram[30] <<
                                                   " " << model_descriptors->at (i).histogram[31] <<
                                                   " " << model_descriptors->at (i).histogram[32] <<
                                                   " " << model_descriptors->at (i).histogram[33] <<  std::endl;
        }*/
        //Debug----------------------------------------------------------------------------------------

        //}

/*
  //
  //  Compute Descriptor for keypoints Spin
  //
  pcl::SpinImageEstimation<PointType, NormalType, DescriptorType> descr_est;
  // Radius of the support cylinder.
  descr_est.setRadiusSearch(descr_rad_);
        // Set the resolution of the spin image (the number of bins along one dimension).
        // Note: you must change the output histogram size to reflect this.
  descr_est.setImageWidth(8);

  descr_est.setInputCloud (model_keypoints);
  descr_est.setInputNormals (model_normals);
  descr_est.compute (*model_descriptors);

  descr_est.setInputCloud (scene_keypoints);
  descr_est.setInputNormals (scene_normals);
  descr_est.compute (*scene_descriptors);
*/

  //
  //  Find Model-Scene Correspondences with KdTree for SHOT
  //
if(use_shot)
{
  pcl::CorrespondencesPtr model_scene_corrs (new pcl::Correspondences ());

  pcl::KdTreeFLANN<DescriptorType> match_search;
  match_search.setInputCloud (model_descriptors);

  //  For each scene keypoint descriptor, find nearest neighbor into the model keypoints descriptor cloud and add it to the correspondences vector.
  //  Para cada descritor de ponto chave da cena, encontra o mais proximo vizinho dentro dos descritores de pontos chave do modelo e junta os ao vector de correspondecias.
  for (size_t i = 0; i < scene_descriptors->size (); ++i)
  {
    std::vector<int> neigh_indices (1);
    std::vector<float> neigh_sqr_dists (1);

    if (!pcl_isfinite (scene_descriptors->at (i).descriptor[0])) //skipping NaNs
    {
      continue;
    }

    int found_neighs = match_search.nearestKSearch (scene_descriptors->at (i), 1, neigh_indices, neigh_sqr_dists);
    if(found_neighs == 1 && neigh_sqr_dists[0] < 0.30f) //  add match only if the squared descriptor distance is less than 0.25 (SHOT descriptor distances are between 0 and 1 by design)
    {
      pcl::Correspondence corr (neigh_indices[0], static_cast<int> (i), neigh_sqr_dists[0]);
      model_scene_corrs->push_back (corr);
    }
  }
  std::cout << "Correspondences found: " << model_scene_corrs->size () << std::endl;
}

  //
  //  Find Model-Scene Correspondences with KdTree for PFH
  //

//if(use_pfh)
//{
  pcl::CorrespondencesPtr model_scene_corrs (new pcl::Correspondences ());

  pcl::KdTreeFLANN<DescriptorType> match_search;
  match_search.setInputCloud (model_descriptors);

  //  For each scene keypoint descriptor, find nearest neighbor into the model keypoints descriptor cloud and add it to the correspondences vector.
  //  Para cada descritor de ponto chave da cena, encontra o mais proximo vizinho dentro dos descritores de pontos chave do modelo e junta os ao vector de correspondecias.
  for (size_t i = 0; i < scene_descriptors->size (); ++i)
  {
    std::vector<int> neigh_indices (1);
    std::vector<float> neigh_sqr_dists (1);

            /*  std::cout << " desccriptors " << i << ":" << scene_descriptors->at (i).histogram[0] <<
                                                   " " << scene_descriptors->at (i).histogram[1] <<
                                                   " " << scene_descriptors->at (i).histogram[2] <<
                                                   " " << scene_descriptors->at (i).histogram[3] <<
                                                   " " << scene_descriptors->at (i).histogram[4] <<
                                                   " " << scene_descriptors->at (i).histogram[5] <<
                                                   " " << scene_descriptors->at (i).histogram[6] <<
                                                   " " << scene_descriptors->at (i).histogram[7] <<
                                                   " " << scene_descriptors->at (i).histogram[8] <<
                                                   " " << scene_descriptors->at (i).histogram[9] <<
                                                   " " << scene_descriptors->at (i).histogram[10] <<
                                                   " " << scene_descriptors->at (i).histogram[11] <<
                                                   " " << scene_descriptors->at (i).histogram[12] <<
                                                   " " << scene_descriptors->at (i).histogram[13] <<
                                                   " " << scene_descriptors->at (i).histogram[14] <<
                                                   " " << scene_descriptors->at (i).histogram[15] <<
                                                   " " << scene_descriptors->at (i).histogram[16] <<
                                                   " " << scene_descriptors->at (i).histogram[17] <<
                                                   " " << scene_descriptors->at (i).histogram[18] <<
                                                   " " << scene_descriptors->at (i).histogram[19] <<
                                                   " " << scene_descriptors->at (i).histogram[20] <<
                                                   " " << scene_descriptors->at (i).histogram[21] <<
                                                   " " << scene_descriptors->at (i).histogram[22] <<
                                                   " " << scene_descriptors->at (i).histogram[23] <<
                                                   " " << scene_descriptors->at (i).histogram[24] <<
                                                   " " << scene_descriptors->at (i).histogram[25] <<
                                                   " " << scene_descriptors->at (i).histogram[26] <<
                                                   " " << scene_descriptors->at (i).histogram[27] <<
                                                   " " << scene_descriptors->at (i).histogram[28] <<
                                                   " " << scene_descriptors->at (i).histogram[29] <<
                                                   " " << scene_descriptors->at (i).histogram[30] <<
                                                   " " << scene_descriptors->at (i).histogram[31] <<
                                                   " " << scene_descriptors->at (i).histogram[32] <<
                                                   " " << scene_descriptors->at (i).histogram[33] <<  std::endl;*/

    if (!pcl_isfinite (scene_descriptors->at (i).descriptor[0])) //skipping NaNs
    {
      continue;
    }

    int found_neighs = match_search.nearestKSearch (scene_descriptors->at (i), 1, neigh_indices, neigh_sqr_dists);

          std::cout << "neigh_sqr_dists: " << neigh_sqr_dists[0] << std::endl; //debugg

    if(found_neighs == 1 && neigh_sqr_dists[0] < Lim_dist_) //  add match only if the squared descriptor distance is less than 0.25 (SHOT descriptor distances are between 0 and 1 by design)
    {
      pcl::Correspondence corr (neigh_indices[0], static_cast<int> (i), neigh_sqr_dists[0]);
      model_scene_corrs->push_back (corr);
    }
  }
  std::cout << "Correspondences found: " << model_scene_corrs->size () << std::endl;
//}

  //
  //  Actual Clustering
  //
  std::vector<Eigen::Matrix4f, Eigen::aligned_allocator<Eigen::Matrix4f> > rototranslations;
  std::vector<pcl::Correspondences> clustered_corrs;

  //  Using Hough3D
  if (use_hough_)
  {
    //
    //  Compute (Keypoints) Reference Frames only for Hough
    //
    pcl::PointCloud<RFType>::Ptr model_rf (new pcl::PointCloud<RFType> ());
    pcl::PointCloud<RFType>::Ptr scene_rf (new pcl::PointCloud<RFType> ());

    pcl::BOARDLocalReferenceFrameEstimation<PointType, NormalType, RFType> rf_est;
    rf_est.setFindHoles (false);
    rf_est.setRadiusSearch (rf_rad_);

    rf_est.setInputCloud (model_keypoints);
    rf_est.setInputNormals (model_keypoints_normals);
    rf_est.setSearchSurface (model_keypoints);
    rf_est.compute (*model_rf);

    rf_est.setInputCloud (scene_keypoints);
    rf_est.setInputNormals (scene_keypoints_normals);
    rf_est.setSearchSurface (scene_keypoints);
    rf_est.compute (*scene_rf);

    //  Clustering
    pcl::Hough3DGrouping<PointType, PointType, RFType, RFType> clusterer;
    clusterer.setHoughBinSize (cg_size_);
    clusterer.setHoughThreshold (cg_thresh_);
    clusterer.setUseInterpolation (true);
    clusterer.setUseDistanceWeight (false);

    clusterer.setInputCloud (model_keypoints);
    clusterer.setInputRf (model_rf);
    clusterer.setSceneCloud (scene_keypoints);
    clusterer.setSceneRf (scene_rf);
    clusterer.setModelSceneCorrespondences (model_scene_corrs);

    //clusterer.cluster (clustered_corrs);
    clusterer.recognize (rototranslations, clustered_corrs);
  }
  else // Using GeometricConsistency
  {
    pcl::GeometricConsistencyGrouping<PointType, PointType> gc_clusterer;
    gc_clusterer.setGCSize (cg_size_);
    gc_clusterer.setGCThreshold (cg_thresh_);

    gc_clusterer.setInputCloud (model_keypoints);
    gc_clusterer.setSceneCloud (scene_keypoints);
    gc_clusterer.setModelSceneCorrespondences (model_scene_corrs);

    //gc_clusterer.cluster (clustered_corrs);
    gc_clusterer.recognize (rototranslations, clustered_corrs);
    std::cout << "H" << std::endl; //debug
  }

  //
  //  Output results
  //
  std::cout << "Model instances found: " << rototranslations.size () << std::endl;
  for (size_t i = 0; i < rototranslations.size (); ++i)
  {
    std::cout << "\n    Instance " << i + 1 << ":" << std::endl;
    std::cout << "        Correspondences belonging to this instance: " << clustered_corrs[i].size () << std::endl;

    // Print the rotation matrix and translation vector
    Eigen::Matrix3f rotation = rototranslations[i].block<3,3>(0, 0);
    Eigen::Vector3f translation = rototranslations[i].block<3,1>(0, 3);

    printf ("\n");
    printf ("            | %6.3f %6.3f %6.3f | \n", rotation (0,0), rotation (0,1), rotation (0,2));
    printf ("        R = | %6.3f %6.3f %6.3f | \n", rotation (1,0), rotation (1,1), rotation (1,2));
    printf ("            | %6.3f %6.3f %6.3f | \n", rotation (2,0), rotation (2,1), rotation (2,2));
    printf ("\n");
    printf ("        t = < %0.3f, %0.3f, %0.3f >\n", translation (0), translation (1), translation (2));
  }

  //
  //  Visualization
  //
  pcl::visualization::PCLVisualizer viewer ("Correspondence Grouping");

  pcl::visualization::PointCloudColorHandlerCustom<PointType> scene_color_handler (scene, 255, 0, 255);
  viewer.addPointCloud (scene,scene_color_handler, "scene_cloud");

  std::string       id = "reference";
  viewer.addCoordinateSystem(0.1, id, 0);

  pcl::PointCloud<PointType>::Ptr off_scene_model (new pcl::PointCloud<PointType> ());
  pcl::PointCloud<PointType>::Ptr off_scene_model_keypoints (new pcl::PointCloud<PointType> ());

  if (show_correspondences_ || show_keypoints_)
  {
    //  We are translating the model so that it doesn't end in the middle of the scene representation
    pcl::transformPointCloud (*model, *off_scene_model, Eigen::Vector3f (-1,0,0), Eigen::Quaternionf (1, 0, 0, 0));
    pcl::transformPointCloud (*model_keypoints, *off_scene_model_keypoints, Eigen::Vector3f (-1,0,0), Eigen::Quaternionf (1, 0, 0, 0));

    pcl::visualization::PointCloudColorHandlerCustom<PointType> off_scene_model_color_handler (off_scene_model, 255, 255, 0);
    viewer.addPointCloud (off_scene_model, off_scene_model_color_handler, "off_scene_model");
  }

//std::cout << off_scene_model_keypoints->size() << std::endl; //debug

  if (show_keypoints_)
  {
    pcl::visualization::PointCloudColorHandlerCustom<PointType> scene_keypoints_color_handler (scene_keypoints, 0, 128, 255);
    viewer.addPointCloud (scene_keypoints, scene_keypoints_color_handler, "scene_keypoints");
    viewer.setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 5, "scene_keypoints");

    pcl::visualization::PointCloudColorHandlerCustom<PointType> off_scene_model_keypoints_color_handler (off_scene_model_keypoints, 0, 128, 255);
    viewer.addPointCloud (off_scene_model_keypoints, off_scene_model_keypoints_color_handler, "off_scene_model_keypoints");
    viewer.setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 5, "off_scene_model_keypoints");
  }

  if (show_normals_)
  {
  viewer.addPointCloudNormals<pcl::PointXYZ, pcl::Normal> (scene_keypoints, scene_keypoints_normals, 1, 0.005, "normals");
  viewer.addPointCloudNormals<pcl::PointXYZ, pcl::Normal> (off_scene_model_keypoints, model_keypoints_normals, 1, 0.005, "normal");
  }

std::cout << rototranslations.size () << std::endl; //debug
  for (size_t i = 0; i < rototranslations.size (); i++)
  {
    pcl::PointCloud<PointType>::Ptr rotated_model (new pcl::PointCloud<PointType> ());
    pcl::transformPointCloud (*model, *rotated_model, rototranslations[i]);

    int inst = i;// + ((v + 1) * 10);

    std::stringstream ss_cloud;
    ss_cloud << "instance" << inst;

    pcl::visualization::PointCloudColorHandlerCustom<PointType> rotated_model_color_handler (rotated_model, 255, 0, 0);
    viewer.addPointCloud (rotated_model, rotated_model_color_handler, ss_cloud.str ());

          // geração da matriz de rotação para os referenciais dos objectos detectados
          Eigen::Affine3f t;

          Eigen::Matrix3f rotation = rototranslations[i].block<3,3>(0, 0);
          Eigen::Vector3f translation = rototranslations[i].block<3,1>(0, 3);

          t.linear() = rotation;
          t.translation() = translation;

          //std::string car = i;

          std::string number;
          std::stringstream strstream;

          //strstream << inst;
          //strstream >> number;
          //std::string &id = i;
number = boost::lexical_cast<std::string>(i);

          viewer.addCoordinateSystem(0.05, t, number, 0);

    if (show_correspondences_)
    {
      for (size_t j = 0; j < clustered_corrs[i].size (); ++j)
      {
        std::stringstream ss_line;
        ss_line << "correspondence_line" << i << "_" << j;
        PointType& model_point = off_scene_model_keypoints->at (clustered_corrs[i][j].index_query);
        PointType& scene_point = scene_keypoints->at (clustered_corrs[i][j].index_match);

        //  We are drawing a line for each pair of clustered correspondences found between the model and the scene
        viewer.addLine<PointType, PointType> (model_point, scene_point, 0, 255, 0, ss_line.str ());
      }
    }
  }

//
//Visualização dos descritores da cena
//
std::cout << "scene_descriptors size:" << scene_descriptors->size() << std::endl;
  for (size_t i = 0; i < scene_descriptors->size (); i++) //i+=10
  {
    std::string number;
    std::string ValDes;

    std::ostringstream ss1;

        ss1 << std::fixed << std::setprecision(4);
        ss1 << i;

          number = boost::lexical_cast<std::string>(i);

          double textScale = 0.001;
          double         r = 1.0;
          double         g = 1.0;
          double         b = 1.0;
          const std::string & id = number;
          int            viewport = 1;

          for(size_t j = 0; j < 125; j++)
          {
            std::ostringstream ss;

                ss << std::fixed << std::setprecision(4);
                ss << scene_descriptors->at(i).descriptor[j];

                if (j==0) {ValDes.append(ss1.str() + ": ");}

          ValDes.append(ss.str() + ", ");

          //ValDes.push_back();

          }

        //  std::cout << ValDes << std::endl;

          //std::cout << " desccriptor " << i << ":" << scene_descriptors->at (i).histogram[5] << std::endl;

          viewer.addText3D   (  number, scene_keypoints->points[i], textScale, r, g, b, id, viewport);
  }

  //
  //Visualização dos descritores do modelo
  //
  std::cout << "model_descriptors size:" << model_descriptors->size() << std::endl;
    for (size_t i = 0; i < model_descriptors->size (); i++) //i+=10
    {
      std::string number;
      std::string ValDes;

      std::ostringstream ss2;

          ss2 << std::fixed << std::setprecision(4);
          ss2 << i;

            number = boost::lexical_cast<std::string>(i);

            double textScale = 0.001;
            double         r = 1.0;
            double         g = 1.0;
            double         b = 1.0;
            const std::string & id = number;
            int            viewport = 1;

            for(size_t j = 0; j < 125; j++)
            {
              std::ostringstream ss;

                  ss << std::fixed << std::setprecision(2);
                  ss << model_descriptors->at(i).descriptor[j];

                  if (j==0) {ValDes.append(ss2.str() + ": ");}

            ValDes.append(ss.str() + ", ");

            //ValDes.push_back();

            }

          //  std::cout << ValDes << std::endl;

            //std::cout << " desccriptor " << i << ":" << scene_descriptors->at (i).histogram[5] << std::endl;

            viewer.addText3D   (  number, off_scene_model_keypoints->points[i], textScale, r, g, b, id, viewport);
    }

  while (!viewer.wasStopped ())
  {
    viewer.spinOnce ();
  }

  return (0);
}


text:
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
// This file is part of the Matrix Template Library
//
// See also license.mtl.txt in the distribution.

#include <iostream>
#include <boost/numeric/mtl/mtl.hpp>

using namespace std;

int main(int, char**)
{
    typedef mtl::dense2D<double> Matrix;

    Matrix A(3,4);
    A[0][0]=1;A[0][1]=2;A[0][2]=3;A[0][3]=4;
    A[1][0]=5;A[1][1]=6;A[1][2]=7;A[1][3]=8;
    A[2][0]=9;A[2][1]=10;A[2][2]=11;A[2][3]=12;

    std::cout << "First Matrix A is: \n" << A << "\n";

    Matrix B = mtl::clone(sub_matrix(A,0,1,0,1));
    std::cout << "Matrix B is: \n" << B << "\n";

    Matrix C = mtl::clone(sub_matrix(A,0,1,0,2));
    std::cout << "Matrix C is: \n" << C << "\n";

    Matrix D = mtl::clone(sub_matrix(A,0,1,0,3));
    std::cout << "Matrix D is: \n" << D << "\n";

    Matrix E = mtl::clone(sub_matrix(A,0,1,0,4));
    std::cout << "Matrix E is: \n" << E << "\n";

    Matrix f = mtl::clone(sub_matrix(A,0,2,0,1));
    std::cout << "Matrix f is: \n" << f << "\n";

    MTL_THROW_IF(f[1][0] != 5.0, mtl::runtime_error("Clown has wrong value."));
    f[1][0]= 6.0;
    std::cout << "Matrix f after f[1][0]= 6.0 is: \n" << f << "\n";
    std::cout << "Matrix A after f[1][0]= 6.0 is: \n" << A << "\n";

    MTL_THROW_IF(f[1][0] != 6.0, mtl::runtime_error("Clown has wrong value after change."));
    MTL_THROW_IF(A[1][0] != 5.0, mtl::runtime_error("Original matrix was changed."));

    Matrix g = mtl::clone(sub_matrix(A,0,3,0,1));
    std::cout << "Matrix g is: \n" << g << "\n";

    Matrix h = mtl::clone(sub_matrix(A,0,3,0,1));
    std::cout << "Matrix h is: \n" << h << "\n";

    Matrix i = mtl::clone(sub_matrix(A,0,2,0,4));
    std::cout << "Matrix i is: \n" << i << "\n";

    Matrix k = mtl::clone(sub_matrix(A,0,2,0,2));
    std::cout << "Matrix k is: \n" << k << "\n";

    Matrix l = mtl::clone(sub_matrix(A,1,3,0,2));
    std::cout << "Matrix l is: \n" << l << "\n";

    return 0;
}