text:
// ---------------------------------------------------------------------
// pion:  a Boost C++ framework for building lightweight HTTP interfaces
// ---------------------------------------------------------------------
// Copyright (C) 2007-2012 Cloudmeter, Inc.  (http://www.cloudmeter.com)
//
// Distributed under the Boost Software License, Version 1.0.
// See http://www.boost.org/LICENSE_1_0.txt
//

#include "EchoService.hpp"
#include <boost/bind.hpp>
#include <pion/algorithm.hpp>
#include <pion/http/response_writer.hpp>
#include <pion/user.hpp>

using namespace pion;

namespace pion {        // begin namespace pion
namespace plugins {     // begin namespace plugins

/// used by handle_request to write dictionary terms
void writeDictionaryTerm(http::response_writer_ptr& writer,
                         const ihash_multimap::value_type& val)
{
    // text is copied into writer text cache
    writer << val.first << http::types::HEADER_NAME_VALUE_DELIMITER
    << val.second
    << http::types::STRING_CRLF;
}

// EchoService member functions

/// handles requests for EchoService
void EchoService::operator()(http::request_ptr& http_request_ptr, tcp::connection_ptr& tcp_conn)
{
    // this web service uses static text to test the mixture of "copied" with
    // "static" (no-copy) text
    static const std::string REQUEST_ECHO_TEXT("[Request Echo]");
    static const std::string REQUEST_HEADERS_TEXT("[Request Headers]");
    static const std::string QUERY_PARAMS_TEXT("[Query Parameters]");
    static const std::string COOKIE_PARAMS_TEXT("[Cookie Parameters]");
    static const std::string POST_CONTENT_TEXT("[POST Content]");
    static const std::string USER_INFO_TEXT("[USER Info]");

    // Set Content-type to "text/plain" (plain ascii text)
    http::response_writer_ptr writer(http::response_writer::create(tcp_conn, *http_request_ptr,
                                                            boost::bind(&tcp::connection::finish, tcp_conn)));
    writer->get_response().set_content_type(http::types::CONTENT_TYPE_TEXT);

    // write request information
    writer->write_no_copy(REQUEST_ECHO_TEXT);
    writer->write_no_copy(http::types::STRING_CRLF);
    writer->write_no_copy(http::types::STRING_CRLF);
    writer
        << "Request method: "
        << http_request_ptr->get_method()
        << http::types::STRING_CRLF
        << "Resource originally requested: "
        << http_request_ptr->get_original_resource()
        << http::types::STRING_CRLF
        << "Resource delivered: "
        << http_request_ptr->get_resource()
        << http::types::STRING_CRLF
        << "Query string: "
        << http_request_ptr->get_query_string()
        << http::types::STRING_CRLF
        << "HTTP version: "
        << http_request_ptr->get_version_major() << '.' << http_request_ptr->get_version_minor()
        << http::types::STRING_CRLF
        << "Content length: "
        << (unsigned long)http_request_ptr->get_content_length()
        << http::types::STRING_CRLF
        << http::types::STRING_CRLF;

    // write request headers
    writer->write_no_copy(REQUEST_HEADERS_TEXT);
    writer->write_no_copy(http::types::STRING_CRLF);
    writer->write_no_copy(http::types::STRING_CRLF);
    std::for_each(http_request_ptr->get_headers().begin(), http_request_ptr->get_headers().end(),
                  boost::bind(&writeDictionaryTerm, writer, _1));
    writer->write_no_copy(http::types::STRING_CRLF);

    // write query parameters
    writer->write_no_copy(QUERY_PARAMS_TEXT);
    writer->write_no_copy(http::types::STRING_CRLF);
    writer->write_no_copy(http::types::STRING_CRLF);
    std::for_each(http_request_ptr->get_queries().begin(), http_request_ptr->get_queries().end(),
                  boost::bind(&writeDictionaryTerm, writer, _1));
    writer->write_no_copy(http::types::STRING_CRLF);

    // write cookie parameters
    writer->write_no_copy(COOKIE_PARAMS_TEXT);
    writer->write_no_copy(http::types::STRING_CRLF);
    writer->write_no_copy(http::types::STRING_CRLF);
    std::for_each(http_request_ptr->get_cookies().begin(), http_request_ptr->get_cookies().end(),
                  boost::bind(&writeDictionaryTerm, writer, _1));
    writer->write_no_copy(http::types::STRING_CRLF);

    // write POST content
    writer->write_no_copy(POST_CONTENT_TEXT);
    writer->write_no_copy(http::types::STRING_CRLF);
    writer->write_no_copy(http::types::STRING_CRLF);
    if (http_request_ptr->get_content_length() != 0) {
        writer->write(http_request_ptr->get_content(), http_request_ptr->get_content_length());
        writer->write_no_copy(http::types::STRING_CRLF);
        writer->write_no_copy(http::types::STRING_CRLF);
    }

    // if authenticated, write user info
    user_ptr user = http_request_ptr->get_user();
    if (user) {
        writer->write_no_copy(USER_INFO_TEXT);
        writer->write_no_copy(http::types::STRING_CRLF);
        writer->write_no_copy(http::types::STRING_CRLF);
        writer << "User authenticated, username: " << user->get_username();
        writer->write_no_copy(http::types::STRING_CRLF);
    }

    // send the writer
    writer->send();
}

}   // end namespace plugins
}   // end namespace pion

/// creates new EchoService objects
extern "C" PION_PLUGIN pion::plugins::EchoService *pion_create_EchoService(void)
{
    return new pion::plugins::EchoService();
}

/// destroys EchoService objects
extern "C" PION_PLUGIN void pion_destroy_EchoService(pion::plugins::EchoService *service_ptr)
{
    delete service_ptr;
}

max_stars_repo_name: paolovolpi/pion

max_issues_repo_name: pmed/pion

max_forks_repo_name: pmed/pion


text:
#include "RendermanPointsSchemaHandler.h"
#include "RendermanPolyMeshSchemaHandler.h"
#include <String2ArgcArgv.h>
#include <boost/program_options.hpp>
#include <boost/shared_ptr.hpp>
#include <stdio.h> // required when building as procedural "RunProgram"
#include <simba_logging.h>

namespace po = boost::program_options;

/*!
 * \todo Remove the i_level param as this is a procedural
 */
void locate_geometry_in_hierarchy(const Alembic::Abc::IObject& top,
								  const StringContainer&       i_hierachy_path,
								  size_t 					   i_requested_index,
								  float          			   i_relative_shutter_open,
								  float          			   i_relative_shutter_close,
								  Alembic::Abc::uint8_t        num_motion_samples,
								  size_t                       i_level = 0)
{
	size_t numChildren = top.getNumChildren();
	std::string unique_object_path;

	for (size_t i=0;i<numChildren;i++)
	{
		std::string child_name =top.getChildHeader(i).getName();
		for (size_t indent=0;indent<i_level;indent++)
			DLOG(INFO) << "  ";
		Alembic::Abc::IObject child(top,child_name);
		const Alembic::Abc::MetaData &child_md = child.getMetaData();
		std::string metadata_string = child_md.serialize();
		DLOG(INFO) << boost::format("metadata_string='%1%' child_name='%2%'") % metadata_string % child_name;
		if (Alembic::AbcGeom::IPolyMeshSchema::matches(child_md))
		{
			RendermanPolyMeshSchemaHandler apmsh;

			apmsh.ProcessPolyMesh(top, child_name, i_hierachy_path, i_requested_index, i_relative_shutter_open, i_relative_shutter_close, num_motion_samples);
		}
		else if (Alembic::AbcGeom::IPointsSchema::matches(child_md))
		{
			RendermanPointsSchemaHandler apsh;
			apsh.ProcessPoints(top, child_name, i_hierachy_path, i_requested_index, i_relative_shutter_open, i_relative_shutter_close, num_motion_samples);
		}
		DLOG(INFO) << std::endl;
		StringContainer       concatenated_hierachy_path = i_hierachy_path;
		concatenated_hierachy_path.push_back(child_name);
		locate_geometry_in_hierarchy(child,
								concatenated_hierachy_path,
								i_requested_index,
								i_relative_shutter_open,
								i_relative_shutter_close,
								num_motion_samples,
								i_level+1);
	}
}

int EmitGeometry(int argc, const char** argv)
{
	try {

	    size_t                frame_to_export = 1;
	    std::string           alembic_fileName;
	    std::string           velocity_fileName;
	    std::string           geometry_filter;
		float                 relative_shutter_open = -0.25f;
		float                 relative_shutter_close = 0.25f;
		size_t                num_motion_samples = 3;

		po::options_description desc("Allowed options");
		desc.add_options()
    		("version", "print version string")
    		("help", "produce help message")
    		("samples,s", po::value<size_t>(&num_motion_samples),
    				"number of samples to generate [default 3]")
    		("open,o", po::value<float>(&relative_shutter_open),
    				"relative shutter open [default -0.25]")
    		("close,c", po::value<float>(&relative_shutter_close),
    				"relative shutter close [default 0.25]")
    		("frame,f", po::value<size_t>(&frame_to_export),
    				"frame to export [default 1]")
    		("abc,a", po::value<std::string>(&alembic_fileName),
    				"alembic file")
			("velocity,v", po::value<std::string>(&velocity_fileName),
					"velocity file")
    				;

		po::variables_map vm;
		po::store(po::parse_command_line(argc, argv, desc), vm);
		po::notify(vm);

		if (vm.count("help") || alembic_fileName.empty()) {
			std::cout << desc << "\n";
			return 1;
		}

		Alembic::AbcCoreFactory::IFactory factory;
		Alembic::AbcCoreFactory::IFactory::CoreType oType;

	    Alembic::AbcGeom::IArchive alembic_archive(factory.getArchive(alembic_fileName, oType));

		StringContainer       hierachy_path;
		locate_geometry_in_hierarchy(alembic_archive.getTop(),hierachy_path,frame_to_export,relative_shutter_open,relative_shutter_close,num_motion_samples);
	}
	catch(std::exception& e) {
		LOG(ERROR) << "error: " << e.what() << "\n";
		return 1;
	}
	catch(...) {
		LOG(ERROR) << "Exception of unknown type!\n";
	}
	return 0;
}

#if defined(__cplusplus) || defined(c_plusplus)
extern "C" {
#endif

RtPointer ConvertParameters(RtString paramstr)
{
	std::string param_std_string = (boost::format("dummy %1%") % paramstr).str();
	DLOG(INFO) << boost::format("param_std_string = '%1%'") % param_std_string << std::endl;
	PI::String2ArgcArgv* s2aa = new PI::String2ArgcArgv(param_std_string);

	/* return the blind data pointer */
	return (RtPointer)s2aa;
}

RtVoid Subdivide(RtPointer data, RtFloat detail)
{
	PI::String2ArgcArgv* s2aa = (PI::String2ArgcArgv *)data;

	EmitGeometry(s2aa->argc(), s2aa->argv());

}

RtVoid Free(RtPointer data)
{
	PI::String2ArgcArgv* s2aa = (PI::String2ArgcArgv *)data;
	delete s2aa;
}

#if defined(__cplusplus) || defined(c_plusplus)
}
#endif

/*!
 * \note This should be harmless when compiled as a procedural "DynamicLoad" but will
 *       be the entry point for procedural "RunProgram"
 */
int main(int argc, char** argv)
{
	google::InitGoogleLogging(argv[0]);
	RtFloat normalize_pixel_coverage;
	DLOG(INFO) << boost::format("normalize_pixel_coverage = %1%") % normalize_pixel_coverage;
	char  args[BUFSIZ];

	while(fgets(args,BUFSIZ,stdin) != NULL) {
		std::string args_string(args);
		args_string.erase(std::remove(args_string.begin(), args_string.end(), '\n'), args_string.end());
		DLOG(INFO) << boost::format("args_string = '%1%'") % args_string << std::endl;
		// std::cerr << boost::format("args = '%1%'") % args << std::endl;
		PI::String2ArgcArgv s2aa(args_string);
		sscanf(args, "%g", &normalize_pixel_coverage);
		RiBegin("stdout");
		RiAttributeBegin();
		EmitGeometry(s2aa.argc(), s2aa.argv());
		RiAttributeEnd();
		RiArchiveRecord(RI_VERBATIM, "\377");
		RiEnd();
	}
	return 0;
}

// == Emacs ================
// -------------------------
// Local variables:
// tab-width: 4
// indent-tabs-mode: t
// c-basic-offset: 4
// end:
//
// == vi ===================
// -------------------------
// Format block
// ex:ts=4:sw=4:expandtab
// -------------------------

ext: cpp


text:
#pragma once

#include <iostream>
#include <boost/thread/thread.hpp>
#include <boost/asio.hpp>
#include <boost/bind.hpp>
#include <vector>

//-------------------------------------------------------------------------------------------------------------------//

/* reference: https://kezunlin.me/post/f241bd30/ */

class ThreadPool {
public:
    ThreadPool(unsigned int num_threads) : io_work(ios) {
        for (int i = 0; i < num_threads; i++) {
            threads.create_thread(boost::bind(&boost::asio::io_service::run, &ios));
        }
    }

    ~ThreadPool() {
        ios.stop();
        try { // suppress all exceptions
            threads.join_all(); // wait for all threads to terminate
        }
        catch ( const std::exception& ) {}
    }

    template<class Function>
    void execute(Function func) {
        ios.post(func); // add the function to the io_service queue
                        // to be run in the threads created in the constructor
        // non-blocking, return immediately

        /* if there aren't available threads in the pool, i.e. every
           thread in the pool is already busy executing some other
           function, the new coming thread has to wait until one of
           the running thread is finished. */

        /* No counter provided to track the number of free threads
           due to implementation overheads with wrapping func with
           a counter addition when function finishes, and at the same
           time supporting boost::bind externally for passing func.
           Refer to this post for more details:
                https://stackoverflow.com/questions/12215395/thread-pool-using-boost-asio/12267138#12267138
        */
       num_tasks++;
    }

    /* total CUMULATIVE number of tasks ever being posted by ThreadPool::execute,
     * which also include those tasks that are already finished
     * */
    unsigned int num_posted_funcs() {
        return num_tasks;
    }

private:
    boost::thread_group threads;
    boost::asio::io_service ios;
    boost::asio::io_service::work io_work;
    unsigned int num_tasks = 0; // this includes those tasks that finished early and got dequeued,
                            // this measurement can't deduce anything about the
                            // the current available free threads

};

//-------------------------------------------------------------------------------------------------------------------//

ext: hpp


text:
/*
  [auto_generated]
  boost/numeric/odeint/iterator/detail/adaptive_time_iterator_controlled_impl.hpp

  [begin_description]
  tba.
  [end_description]

  Copyright 2009-2012 Karsten Ahnert
  Copyright 2009-2012 Mario Mulansky

  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or
  copy at http://www.boost.org/LICENSE_1_0.txt)
*/

#ifndef BOOST_NUMERIC_ODEINT_ITERATOR_DETAIL_ADAPTIVE_TIME_ITERATOR_CONTROLLED_IMPL_HPP_DEFINED
#define BOOST_NUMERIC_ODEINT_ITERATOR_DETAIL_ADAPTIVE_TIME_ITERATOR_CONTROLLED_IMPL_HPP_DEFINED

#include <boost/numeric/odeint/util/unit_helper.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/iterator/detail/ode_time_iterator_base.hpp>

namespace boost {
namespace numeric {
namespace odeint {

    /*
     * Specilization for controlled steppers
     */
    /**
     * \brief ODE Iterator with adaptive step size control. The value type of this iterator is a pair of state type and time type of the stepper.
     *
     * Implements an ODE iterator with adaptive step size control. Uses controlled steppers. adaptive_iterator is a model
     * of single-pass iterator.
     *
     * The value type of this iterator is a pair of state type and time type of the stepper.
     *
     * \tparam Stepper The stepper type which should be used during the iteration.
     * \tparam System The type of the system function (ODE) which should be solved.
     */
    template< class Stepper , class System >
    class adaptive_time_iterator< Stepper , System , controlled_stepper_tag > : public detail::ode_time_iterator_base
    <
        adaptive_time_iterator< Stepper , System , controlled_stepper_tag > ,
        Stepper , System
    >
    {
    private:

        typedef Stepper stepper_type;
        typedef System system_type;
        typedef typename boost::numeric::odeint::unwrap_reference< stepper_type >::type unwrapped_stepper_type;
        typedef typename unwrapped_stepper_type::state_type state_type;
        typedef typename unwrapped_stepper_type::time_type time_type;
        typedef typename unwrapped_stepper_type::value_type ode_value_type;
        typedef detail::ode_time_iterator_base<
            adaptive_time_iterator< Stepper , System , controlled_stepper_tag > ,
            Stepper , System > base_type;

    public:

        /**
         * \brief Constructs an adaptive_time_iterator. This constructor should be used to construct the begin iterator.
         *
         * \param stepper The stepper to use during the iteration.
         * \param sys The system function (ODE) to solve.
         * \param s The initial state. adaptive_time_iterator stores a reference of s and changes its value during the iteration.
         * \param t The initial time.
         * \param t_end The end time, at which the iteration should stop.
         * \param dt The initial time step.
         */
        adaptive_time_iterator( stepper_type stepper , system_type sys , state_type &s , time_type t , time_type t_end , time_type dt )
            : base_type( stepper , sys , s , t , t_end , dt ) {}

        /**
         * \brief Constructs an adaptive_time_iterator. This constructor should be used to construct the end iterator.
         *
         * \param stepper The stepper to use during the iteration.
         * \param sys The system function (ODE) to solve.
         * \param s The initial state. adaptive_time_iterator stores a reference of s and changes its value during the iteration.
         */
        adaptive_time_iterator( stepper_type stepper , system_type sys , state_type &s )
            : base_type( stepper , sys , s ) {}

    private:

        friend class boost::iterator_core_access;

        void increment()
        {
            unwrapped_stepper_type &stepper = this->m_stepper;
            const size_t max_attempts = 1000;
            size_t trials = 0;
            controlled_step_result res = success;
            do
            {
                res = stepper.try_step( this->m_system , *this->m_state , this->m_t , this->m_dt );
                ++trials;
            }
            while( ( res == fail ) && ( trials < max_attempts ) );
            if( trials == max_attempts )
            {
                throw std::overflow_error( "Adaptive iterator : Maximal number of iterations reached. A step size could not be found." );
            }
            this->check_end();
        }
    };

} // namespace odeint
} // namespace numeric
} // namespace boost

#endif // BOOST_NUMERIC_ODEINT_ITERATOR_DETAIL_ADAPTIVE_TIME_ITERATOR_CONTROLLED_IMPL_HPP_DEFINED

max_stars_repo_name: MINATILO/packing-generation

max_issues_repo_name: MINATILO/packing-generation

max_forks_repo_name: MINATILO/packing-generation


text:
#pragma once

#include "conf_make_install.hpp"

#include <bunsan/stream_enum.hpp>
#include <bunsan/utility/maker.hpp>

#include <boost/optional.hpp>
#include <boost/serialization/access.hpp>
#include <boost/serialization/nvp.hpp>

#include <string>
#include <unordered_map>
#include <vector>

namespace bunsan::utility::builders {

class meson : public conf_make_install {
 public:
  struct config {
    template <typename Archive>
    void serialize(Archive &ar, const unsigned int) {
      ar & BOOST_SERIALIZATION_NVP(meson);
      ar & BOOST_SERIALIZATION_NVP(maker);
    }

    struct {
      template <typename Archive>
      void serialize(Archive &ar, const unsigned int) {
        ar & BOOST_SERIALIZATION_NVP(options);
        ar & BOOST_SERIALIZATION_NVP(flags);
      }

      /// for example {"prefix": "/usr"}
      /// translated into -Dprefix=/usr
      std::unordered_map<std::string, std::string> options;
      /// for example {"buildtype": "release"}
      /// translated into --buildtype=release
      std::unordered_map<std::string, std::string> flags;
    } meson;

    boost::property_tree::ptree maker;
  };

 public:
  meson(const utility_config &ptree, resolver &resolver_);

 protected:
  void configure_(const boost::filesystem::path &src,
                  const boost::filesystem::path &bin) override;

  void make_(const boost::filesystem::path &src,
             const boost::filesystem::path &bin) override;

  void install_(const boost::filesystem::path &src,
                const boost::filesystem::path &bin,
                const boost::filesystem::path &root) override;

 private:
  std::vector<std::string> arguments_(const boost::filesystem::path &src) const;

 private:
  const config m_config;
  const boost::filesystem::path m_meson_exe;
  const maker_ptr m_maker;
};

struct meson_error : virtual error {};

}  // namespace bunsan::utility::builders

max_stars_repo_name: bacsorg/bacs

max_issues_repo_name: bacsorg/bacs

max_forks_repo_name: bacsorg/bacs


text:
#include <libbruce/be/disk.h>
#include <libbruce/util/be_registry.h>
#include <boost/make_shared.hpp>
#include <cstdio>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
#include <errno.h>

void sha1_compress(uint32_t state[5], const uint8_t block[64]) asm ("sha1_compress");

void sha1_hash(const uint8_t *message, uint32_t len, uint32_t hash[5])
{
	hash[0] = UINT32_C(0x67452301);
	hash[1] = UINT32_C(0xEFCDAB89);
	hash[2] = UINT32_C(0x98BADCFE);
	hash[3] = UINT32_C(0x10325476);
	hash[4] = UINT32_C(0xC3D2E1F0);

	uint32_t i;
	for (i = 0; len - i >= 64; i += 64)
		sha1_compress(hash, message + i);

	uint8_t block[64];
	uint32_t rem = len - i;
	memcpy(block, message + i, rem);

	block[rem] = 0x80;
	rem++;
	if (64 - rem >= 8)
		memset(block + rem, 0, 56 - rem);
	else {
		memset(block + rem, 0, 64 - rem);
		sha1_compress(hash, block);
		memset(block, 0, 56);
	}

	uint64_t longLen = ((uint64_t)len) << 3;
	for (i = 0; i < 8; i++)
		block[64 - 1 - i] = (uint8_t)(longLen >> (i * 8));
	sha1_compress(hash, block);
}

namespace libbruce { namespace be {

disk::disk(std::string pathPrefix, uint32_t maxBlockSize, uint32_t editQueueSize)
    : m_pathPrefix(pathPrefix), m_maxBlockSize(maxBlockSize), m_editQueueSize(editQueueSize)
{
}

mempage disk::get(const nodeid_t &id)
{
    int f = open((m_pathPrefix + boost::lexical_cast<std::string>(id)).c_str(), O_RDONLY);
    if (f == -1) throw std::runtime_error("Error opening file for reading");

    struct stat stat_info;
    if (fstat(f, &stat_info) != 0)
        throw std::runtime_error("Error stat'ing file");

    mempage ret(stat_info.st_size);
    read(f, ret.ptr(), stat_info.st_size); // FIXME: Check return code
    close(f);
    return ret;
}

getblockresult_t disk::get_all(const blockidlist_t &ids)
{
    getblockresult_t ret;
    for (blockidlist_t::const_iterator it = ids.begin(); it != ids.end(); ++it)
        ret[*it] = get(*it);
    return ret;
}

nodeid_t disk::id(const mempage &block)
{
    uint32_t hash[5];
    sha1_hash(block.ptr(), block.size(), hash);

    BOOST_STATIC_ASSERT(sizeof(nodeid_t) == sizeof(hash));

    nodeid_t ret;
    memcpy(ret.data(), hash, sizeof(hash));
    return ret;
}

void disk::put_all(putblocklist_t &blocklist)
{
    if (!m_maxBlockSize)
        throw be_error("Can't put; engine is read-only");

    for (int i = 0; i < blocklist.size(); i++)
    {
        put_one(blocklist[i]);
    }
}

void disk::put_one(putblock_t &block)
{
    int f = open((m_pathPrefix + boost::lexical_cast<std::string>(block.id)).c_str(), O_WRONLY|O_CREAT, 0644);
    // Used to use fstream here but it's too slow
    if (f == -1) throw std::runtime_error(strerror(errno));
    write(f, block.mem.ptr(), block.mem.size()); // FIXME: Check return code
    close(f);

    block.success = true;
}

void disk::del_all(delblocklist_t &ids)
{
    if (!m_maxBlockSize)
        throw be_error("Can't delete; engine is read-only");

    for (int i = 0; i < ids.size(); i++)
    {
        del_one(ids[i]);
    }
}

void disk::del_one(delblock_t &block)
{
    std::remove((m_pathPrefix + boost::lexical_cast<std::string>(block.id)).c_str());
    block.success = true;
}

uint32_t disk::maxBlockSize()
{
    return m_maxBlockSize;
}

uint32_t disk::editQueueSize()
{
    return m_editQueueSize;
}

be_ptr create_disk_engine(const std::string &location, size_t block_size, size_t queue_size, const util::options_t &options)
{
    return boost::make_shared<disk>(location, block_size, queue_size);
}

void register_disk_engine()
{
    util::register_be_factory("file", &create_disk_engine);
}

}}


text:
// File generated on Fri Aug 22, 2014 10:40:26 AM by xcpp.
/*Armadillo Test

Collection of sample functions showing how to use the
[Armadillo](http://arma.sourceforge.net/) library.*/
#define XC_PYTHON
#include <armadillo>
#include <xc>

extern "C" {
    void inverse_det_py(size_t, char*, size_t&, char*&);
    void inverse_det_py_clear();
    void qr_py(size_t, char*, size_t&, char*&);
    void qr_py_clear();
    void solve_py(size_t, char*, size_t&, char*&);
    void solve_py_clear();

}

std::string inverse_det_py_str;
void inverse_det_py(size_t ncin, char* pcin, size_t& ncout, char*& pcout) {
try {
    excentury::STextInterface<excentury::load_mode> XC_LI_(pcin, ncin);
    arma::mat A(2, 2); XC_LI_.load(A, A(0,0));
    XC_LI_.close();

    double d = det(A);
    arma::mat Ainv;
    try {
        Ainv = inv(A);
    } catch (std::runtime_error& run_error) {
        excentury::error(run_error.what());
    }

    excentury::STextInterface<excentury::dump_mode> XC_DI_;
    XC_DI_.dump(Ainv, "Ainv", Ainv(0, 0));
    XC_DI_.dump(d, "detA");
    XC_DI_.close();
    inverse_det_py_str = XC_DI_.str();
    ncout = inverse_det_py_str.size();
    pcout = (char*)inverse_det_py_str.data();
} catch (excentury::RuntimeError& run_error) {
    inverse_det_py_str = run_error.msg;
    ncout = inverse_det_py_str.size();
    pcout = (char*)inverse_det_py_str.data();
}
}
void inverse_det_py_clear() {
    inverse_det_py_str.clear();
}

std::string qr_py_str;
void qr_py(size_t ncin, char* pcin, size_t& ncout, char*& pcout) {
try {
    excentury::STextInterface<excentury::load_mode> XC_LI_(pcin, ncin);
    arma::mat X(2, 2); XC_LI_.load(X, X(0,0));
    XC_LI_.close();

    arma::mat Q, R;
    int status = arma::qr(Q, R, X);

    excentury::STextInterface<excentury::dump_mode> XC_DI_;
    XC_DI_.dump(Q, "Q", Q(0, 0));
    XC_DI_.dump(R, "R", R(0, 0));
    XC_DI_.dump(status, "status");
    XC_DI_.close();
    qr_py_str = XC_DI_.str();
    ncout = qr_py_str.size();
    pcout = (char*)qr_py_str.data();
} catch (excentury::RuntimeError& run_error) {
    qr_py_str = run_error.msg;
    ncout = qr_py_str.size();
    pcout = (char*)qr_py_str.data();
}
}
void qr_py_clear() {
    qr_py_str.clear();
}

std::string solve_py_str;
void solve_py(size_t ncin, char* pcin, size_t& ncout, char*& pcout) {
try {
    excentury::STextInterface<excentury::load_mode> XC_LI_(pcin, ncin);
    arma::mat A(2, 2); XC_LI_.load(A, A(0,0));
    arma::mat B(2, 2); XC_LI_.load(B, B(0,0));
    XC_LI_.close();

    arma::mat X;
    try {
        X = arma::solve(A,B);
    } catch (std::exception& run_error) {
        excentury::error(run_error.what());
    }

    excentury::STextInterface<excentury::dump_mode> XC_DI_;
    XC_DI_.dump(X, "X", X(0, 0));
    XC_DI_.close();
    solve_py_str = XC_DI_.str();
    ncout = solve_py_str.size();
    pcout = (char*)solve_py_str.data();
} catch (excentury::RuntimeError& run_error) {
    solve_py_str = run_error.msg;
    ncout = solve_py_str.size();
    pcout = (char*)solve_py_str.data();
}
}
void solve_py_clear() {
    solve_py_str.clear();
}


text:
// Copyright (c) 2015-2016, ETH Zurich, Wyss Zurich, Zurich Eye
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//     * 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.
//     * Neither the name of the ETH Zurich, Wyss Zurich, Zurich Eye 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 ETH Zurich, Wyss Zurich, Zurich Eye 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 <ze/common/matrix.hpp>

#include <Eigen/SVD>
#include <Eigen/LU>

namespace ze {

std::tuple<int, real_t, VectorX> directLinearTransform(const MatrixX& A, real_t rank_tol)
{
  int n = A.rows();
  int p = A.cols();
  int m = std::min(n,p);
  Eigen::JacobiSVD<MatrixX> svd(A, Eigen::ComputeFullV);
  VectorX s = svd.singularValues();
  MatrixX V = svd.matrixV();

  // Find rank
  int rank = 0;
  for (int i = 0; i < m; ++i)
  {
    if (s(i) > rank_tol)
    {
      ++rank;
    }
  }

  // Return rank, error, and corresponding column of V
  real_t error = (m < p) ? 0.0 : s(m-1);
  return std::make_tuple(rank, error, V.col(p-1));
}

} // namespace ze


text:
/**************************************************************
 *
 * Licensed to the Apache Software Foundation (ASF) under one
 * or more contributor license agreements.  See the NOTICE file
 * distributed with this work for additional information
 * regarding copyright ownership.  The ASF licenses this file
 * to you 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.
 *
 *************************************************************/

// MARKER(update_precomp.py): autogen include statement, do not remove
#include "precompiled_canvas.hxx"

#include <ctype.h> // don't ask. msdev breaks otherwise...
#include <basegfx/numeric/ftools.hxx>

#include <canvas/debug.hxx>
#include <canvas/verbosetrace.hxx>
#include <tools/diagnose_ex.h>

#include <com/sun/star/lang/XServiceInfo.hpp>
#include <com/sun/star/lang/XUnoTunnel.hpp>
#include <com/sun/star/geometry/RealPoint2D.hpp>
#include <com/sun/star/geometry/IntegerRectangle2D.hpp>

#include <basegfx/matrix/b2dhommatrix.hxx>
#include <basegfx/range/b2irectangle.hxx>
#include <basegfx/range/b2drectangle.hxx>
#include <basegfx/polygon/b2dpolygon.hxx>
#include <basegfx/polygon/b2dpolypolygon.hxx>
#include <basegfx/tools/canvastools.hxx>

#include <canvas/canvastools.hxx>
#include <canvas/verifyinput.hxx>

#include "dx_impltools.hxx"
#include "dx_vcltools.hxx"
#include "dx_linepolypolygon.hxx"
#include "dx_canvasbitmap.hxx"
#include "dx_canvasfont.hxx"
#include "dx_canvas.hxx"
#include "dx_spritecanvas.hxx"

#include <boost/scoped_array.hpp>

#include <vector>
#include <algorithm>

using namespace ::com::sun::star;

namespace dxcanvas
{
    namespace tools
    {
        ::basegfx::B2DPolyPolygon polyPolygonFromXPolyPolygon2D( const uno::Reference< rendering::XPolyPolygon2D >& xPoly )
        {
            LinePolyPolygon* pPolyImpl = dynamic_cast< LinePolyPolygon* >( xPoly.get() );

            if( pPolyImpl )
            {
                return pPolyImpl->getPolyPolygon();
            }
            else
            {
                const sal_Int32 nPolys( xPoly->getNumberOfPolygons() );

                // not a known implementation object - try data source
                // interfaces
                uno::Reference< rendering::XBezierPolyPolygon2D > xBezierPoly(
                    xPoly,
                    uno::UNO_QUERY );

                if( xBezierPoly.is() )
                {
                    return ::basegfx::unotools::polyPolygonFromBezier2DSequenceSequence(
                        xBezierPoly->getBezierSegments( 0,
                                                        nPolys,
                                                        0,
                                                        -1 ) );
                }
                else
                {
                    uno::Reference< rendering::XLinePolyPolygon2D > xLinePoly(
                        xPoly,
                        uno::UNO_QUERY );

                    // no implementation class and no data provider
                    // found - contract violation.
                    ENSURE_ARG_OR_THROW( xLinePoly.is(),
                                     "VCLCanvas::polyPolygonFromXPolyPolygon2D(): Invalid input "
                                     "poly-polygon, cannot retrieve vertex data" );

                    return ::basegfx::unotools::polyPolygonFromPoint2DSequenceSequence(
                        xLinePoly->getPoints( 0,
                                              nPolys,
                                              0,
                                              -1 ) );
                }
            }
        }

        void setupGraphics( Gdiplus::Graphics& rGraphics )
        {
            // setup graphics with (somewhat arbitrary) defaults
            //rGraphics.SetCompositingQuality( Gdiplus::CompositingQualityHighQuality );
            rGraphics.SetCompositingQuality( Gdiplus::CompositingQualityHighSpeed );
            //rGraphics.SetInterpolationMode( Gdiplus::InterpolationModeHighQualityBilinear ); // with prefiltering for shrinks
            rGraphics.SetInterpolationMode( Gdiplus::InterpolationModeBilinear );

            // #122683# Switched precedence of pixel offset
            // mode. Seemingly, polygon stroking needs
            // PixelOffsetModeNone to achieve visually pleasing
            // results, whereas all other operations (e.g. polygon
            // fills, bitmaps) look better with PixelOffsetModeHalf.
            rGraphics.SetPixelOffsetMode( Gdiplus::PixelOffsetModeHalf ); // Pixel center at (0.5, 0.5) etc.
            //rGraphics.SetPixelOffsetMode( Gdiplus::PixelOffsetModeNone );

            //rGraphics.SetSmoothingMode( Gdiplus::SmoothingModeHighSpeed ); // no line/curve antialiasing
            //rGraphics.SetSmoothingMode( Gdiplus::SmoothingModeHighQuality );
            rGraphics.SetSmoothingMode( Gdiplus::SmoothingModeAntiAlias );
            //rGraphics.SetTextRenderingHint( Gdiplus::TextRenderingHintAntiAlias );
            rGraphics.SetTextRenderingHint( Gdiplus::TextRenderingHintSystemDefault );
			rGraphics.SetPageUnit(Gdiplus::UnitPixel);
        }

        Gdiplus::Graphics* createGraphicsFromHDC(HDC aHDC)
        {
            Gdiplus::Graphics* pRet = new Gdiplus::Graphics(aHDC);
			if( pRet )
                setupGraphics( *pRet );
            return pRet;
        }

        Gdiplus::Graphics* createGraphicsFromBitmap(const BitmapSharedPtr& rBitmap)
        {
            Gdiplus::Graphics* pRet = Gdiplus::Graphics::FromImage(rBitmap.get());
            if( pRet )
                setupGraphics( *pRet );
            return pRet;
        }

        void gdiPlusMatrixFromB2DHomMatrix( Gdiplus::Matrix& rGdiplusMatrix, const ::basegfx::B2DHomMatrix& rMatrix )
        {
            rGdiplusMatrix.SetElements( static_cast<Gdiplus::REAL>(rMatrix.get(0,0)),
                                        static_cast<Gdiplus::REAL>(rMatrix.get(1,0)),
                                        static_cast<Gdiplus::REAL>(rMatrix.get(0,1)),
                                        static_cast<Gdiplus::REAL>(rMatrix.get(1,1)),
                                        static_cast<Gdiplus::REAL>(rMatrix.get(0,2)),
                                        static_cast<Gdiplus::REAL>(rMatrix.get(1,2)) );
        }

        void gdiPlusMatrixFromAffineMatrix2D( Gdiplus::Matrix& 					rGdiplusMatrix,
                                              const geometry::AffineMatrix2D&	rMatrix )
        {
            rGdiplusMatrix.SetElements( static_cast<Gdiplus::REAL>(rMatrix.m00),
                                        static_cast<Gdiplus::REAL>(rMatrix.m10),
                                        static_cast<Gdiplus::REAL>(rMatrix.m01),
                                        static_cast<Gdiplus::REAL>(rMatrix.m11),
                                        static_cast<Gdiplus::REAL>(rMatrix.m02),
                                        static_cast<Gdiplus::REAL>(rMatrix.m12) );
        }

        namespace
        {
            // TODO(P2): Check whether this gets inlined. If not, make functor
            // out of it
            inline Gdiplus::PointF implGdiPlusPointFromRealPoint2D( const ::com::sun::star::geometry::RealPoint2D& rPoint )
            {
                return Gdiplus::PointF( static_cast<Gdiplus::REAL>(rPoint.X),
                                        static_cast<Gdiplus::REAL>(rPoint.Y) );
            }

            void graphicsPathFromB2DPolygon( GraphicsPathSharedPtr&             rOutput,
                                             ::std::vector< Gdiplus::PointF >&  rPoints,
                                             const ::basegfx::B2DPolygon&       rPoly,
                                             bool bNoLineJoin)
            {
                const sal_uInt32 nPoints( rPoly.count() );

                if( nPoints < 2 )
                    return;

                rOutput->StartFigure();

                const bool bClosedPolygon( rPoly.isClosed() );

                if( rPoly.areControlPointsUsed() )
                {
                    // control points used -> for now, add all
                    // segments as curves to GraphicsPath

                    // If the polygon is closed, we need to add the
                    // first point, thus, one more (can't simply
                    // GraphicsPath::CloseFigure() it, since the last
                    // point cannot have any control points for GDI+)
                    rPoints.resize( 3*nPoints + bClosedPolygon );

                    sal_uInt32 nCurrOutput=0;
                    for( sal_uInt32 nCurrPoint=0; nCurrPoint<nPoints; ++nCurrPoint )
                    {
                        const ::basegfx::B2DPoint& rPoint( rPoly.getB2DPoint( nCurrPoint ) );
                        rPoints[nCurrOutput++] = Gdiplus::PointF( static_cast<Gdiplus::REAL>(rPoint.getX()),
                                                                  static_cast<Gdiplus::REAL>(rPoint.getY()) );

                        const ::basegfx::B2DPoint& rControlPointA( rPoly.getNextControlPoint( nCurrPoint ) );
                        rPoints[nCurrOutput++] = Gdiplus::PointF( static_cast<Gdiplus::REAL>(rControlPointA.getX()),
                                                                  static_cast<Gdiplus::REAL>(rControlPointA.getY()) );

                        const ::basegfx::B2DPoint& rControlPointB( rPoly.getPrevControlPoint( (nCurrPoint + 1) % nPoints) );
                        rPoints[nCurrOutput++] = Gdiplus::PointF( static_cast<Gdiplus::REAL>(rControlPointB.getX()),
                                                                  static_cast<Gdiplus::REAL>(rControlPointB.getY()) );
                    }

                    if( bClosedPolygon )
                    {
                        // add first point again (to be able to pass
                        // control points for the last point, see
                        // above)
                        const ::basegfx::B2DPoint& rPoint( rPoly.getB2DPoint(0) );
                        rPoints[nCurrOutput++] = Gdiplus::PointF( static_cast<Gdiplus::REAL>(rPoint.getX()),
                                                                  static_cast<Gdiplus::REAL>(rPoint.getY()) );

                        if(bNoLineJoin && nCurrOutput > 7)
                        {
                            for(sal_uInt32 a(3); a < nCurrOutput; a+=3)
                            {
        				        rOutput->StartFigure();
                                rOutput->AddBezier(rPoints[a - 3], rPoints[a - 2], rPoints[a - 1], rPoints[a]);
                            }
                        }
                        else
                        {
                            rOutput->AddBeziers( &rPoints[0], nCurrOutput );
                        }
                    }
                    else
                    {
                        // GraphicsPath expects 3(n-1)+1 points (i.e. the
                        // last point must not have any trailing control
                        // points after it).
                        // Therefore, simply don't pass the last two
                        // points here.
                        if( nCurrOutput > 3 )
                        {
                            if(bNoLineJoin && nCurrOutput > 7)
                            {
                                for(sal_uInt32 a(3); a < nCurrOutput; a+=3)
                                {
        				            rOutput->StartFigure();
                                    rOutput->AddBezier(rPoints[a - 3], rPoints[a - 2], rPoints[a - 1], rPoints[a]);
                                }
                            }
                            else
                            {
                                rOutput->AddBeziers( &rPoints[0], nCurrOutput-2 );
                            }
                        }
                    }
                }
                else
                {
                    // no control points -> no curves, simply add
                    // straigt lines to GraphicsPath
                    rPoints.resize( nPoints );

                    for( sal_uInt32 nCurrPoint=0; nCurrPoint<nPoints; ++nCurrPoint )
                    {
                        const ::basegfx::B2DPoint& rPoint( rPoly.getB2DPoint( nCurrPoint ) );
                        rPoints[nCurrPoint] = Gdiplus::PointF( static_cast<Gdiplus::REAL>(rPoint.getX()),
                                                               static_cast<Gdiplus::REAL>(rPoint.getY()) );
                    }

                    if(bNoLineJoin && nPoints > 2)
                    {
                        for(sal_uInt32 a(1); a < nPoints; a++)
                        {
        				    rOutput->StartFigure();
                            rOutput->AddLine(rPoints[a - 1], rPoints[a]);
                        }

                        if(bClosedPolygon)
                        {
        				    rOutput->StartFigure();
                            rOutput->AddLine(rPoints[nPoints - 1], rPoints[0]);
                        }
                    }
                    else
                    {
                        rOutput->AddLines( &rPoints[0], nPoints );
                    }
                }

                if( bClosedPolygon && !bNoLineJoin )
                    rOutput->CloseFigure();
            }
        }

        Gdiplus::PointF gdiPlusPointFromRealPoint2D( const ::com::sun::star::geometry::RealPoint2D& rPoint )
        {
            return implGdiPlusPointFromRealPoint2D( rPoint );
        }

        Gdiplus::Rect gdiPlusRectFromIntegerRectangle2D( const geometry::IntegerRectangle2D& rRect )
        {
            return Gdiplus::Rect( rRect.X1,
                                  rRect.Y1,
                                  rRect.X2 - rRect.X1,
                                  rRect.Y2 - rRect.Y1 );
        }

        Gdiplus::RectF gdiPlusRectFFromRectangle2D( const geometry::RealRectangle2D& rRect )
        {
            return Gdiplus::RectF( static_cast<Gdiplus::REAL>(rRect.X1),
                                   static_cast<Gdiplus::REAL>(rRect.Y1),
                                   static_cast<Gdiplus::REAL>(rRect.X2 - rRect.X1),
                                   static_cast<Gdiplus::REAL>(rRect.Y2 - rRect.Y1) );
        }

        RECT gdiRectFromB2IRect( const ::basegfx::B2IRange& rRect )
        {
            RECT aRect = {rRect.getMinX(),
                          rRect.getMinY(),
                          rRect.getMaxX(),
                          rRect.getMaxY()};

            return aRect;
        }

        geometry::RealPoint2D realPoint2DFromGdiPlusPointF( const Gdiplus::PointF& rPoint )
        {
            return geometry::RealPoint2D( rPoint.X, rPoint.Y );
        }

        geometry::RealRectangle2D realRectangle2DFromGdiPlusRectF( const Gdiplus::RectF& rRect )
        {
            return geometry::RealRectangle2D( rRect.X, rRect.Y,
                                              rRect.X + rRect.Width,
                                              rRect.Y + rRect.Height );
        }

        ::basegfx::B2DPoint	b2dPointFromGdiPlusPointF( const Gdiplus::PointF& rPoint )
        {
            return ::basegfx::B2DPoint( rPoint.X, rPoint.Y );
        }

        ::basegfx::B2DRange	b2dRangeFromGdiPlusRectF( const Gdiplus::RectF& rRect )
        {
            return ::basegfx::B2DRange( rRect.X, rRect.Y,
                                        rRect.X + rRect.Width,
                                        rRect.Y + rRect.Height );
        }

        uno::Sequence< double > argbToDoubleSequence( const Gdiplus::ARGB& rColor )
        {
            // TODO(F1): handle color space conversions, when defined on canvas/graphicDevice
            uno::Sequence< double > aRet(4);

            aRet[0] = ((rColor >> 16) & 0xFF) / 255.0; 	// red
            aRet[1] = ((rColor >> 8) & 0xFF) / 255.0;	// green
            aRet[2] = (rColor & 0xFF) / 255.0;			// blue
            aRet[3] = ((rColor >> 24) & 0xFF) / 255.0;	// alpha

            return aRet;
        }

        uno::Sequence< sal_Int8 > argbToIntSequence( const Gdiplus::ARGB& rColor )
        {
            // TODO(F1): handle color space conversions, when defined on canvas/graphicDevice
            uno::Sequence< sal_Int8 > aRet(4);

            aRet[0] = static_cast<sal_Int8>((rColor >> 16) & 0xFF);	// red
            aRet[1] = static_cast<sal_Int8>((rColor >> 8) & 0xFF);	// green
            aRet[2] = static_cast<sal_Int8>(rColor & 0xFF);			// blue
            aRet[3] = static_cast<sal_Int8>((rColor >> 24) & 0xFF);	// alpha

            return aRet;
        }

        Gdiplus::ARGB sequenceToArgb( const uno::Sequence< sal_Int8 >& rColor )
        {
            ENSURE_OR_THROW( rColor.getLength() > 2,
                              "sequenceToArgb: need at least three channels" );

            // TODO(F1): handle color space conversions, when defined on canvas/graphicDevice
            Gdiplus::ARGB aColor;

            aColor = (static_cast<sal_uInt8>(rColor[0]) << 16) | (static_cast<sal_uInt8>(rColor[1]) << 8) | static_cast<sal_uInt8>(rColor[2]);

            if( rColor.getLength() > 3 )
                aColor |= static_cast<sal_uInt8>(rColor[3]) << 24;

            return aColor;
        }

        Gdiplus::ARGB sequenceToArgb( const uno::Sequence< double >& rColor )
        {
            ENSURE_OR_THROW( rColor.getLength() > 2,
                              "sequenceToColor: need at least three channels" );

            // TODO(F1): handle color space conversions, when defined on canvas/graphicDevice
            Gdiplus::ARGB aColor;

            ::canvas::tools::verifyRange(rColor[0],0.0,1.0);
            ::canvas::tools::verifyRange(rColor[1],0.0,1.0);
            ::canvas::tools::verifyRange(rColor[2],0.0,1.0);

            aColor =
                (static_cast<sal_uInt8>( ::basegfx::fround( 255*rColor[0] ) ) << 16) |
                (static_cast<sal_uInt8>( ::basegfx::fround( 255*rColor[1] ) ) << 8) |
                static_cast<sal_uInt8>( ::basegfx::fround( 255*rColor[2] ) );

            if( rColor.getLength() > 3 )
            {
                ::canvas::tools::verifyRange(rColor[3],0.0,1.0);
                aColor |= static_cast<sal_uInt8>( ::basegfx::fround( 255*rColor[3] ) ) << 24;
            }

            return aColor;
        }

        GraphicsPathSharedPtr graphicsPathFromRealPoint2DSequence( const uno::Sequence< uno::Sequence< geometry::RealPoint2D > >& points )
        {
            GraphicsPathSharedPtr pRes( new Gdiplus::GraphicsPath() );
            ::std::vector< Gdiplus::PointF > aPoints;

            sal_Int32 nCurrPoly;
            for( nCurrPoly=0; nCurrPoly<points.getLength(); ++nCurrPoly )
            {
                const sal_Int32 nCurrSize( points[nCurrPoly].getLength() );
                if( nCurrSize )
                {
                    aPoints.resize( nCurrSize );

                    // TODO(F1): Closed/open polygons

                    // convert from RealPoint2D array to Gdiplus::PointF array
                    ::std::transform( const_cast< uno::Sequence< geometry::RealPoint2D >& >(points[nCurrPoly]).getArray(),
                                      const_cast< uno::Sequence< geometry::RealPoint2D >& >(points[nCurrPoly]).getArray()+nCurrSize,
                                      aPoints.begin(),
                                      implGdiPlusPointFromRealPoint2D );

                    pRes->AddLines( &aPoints[0], nCurrSize );
                }
            }

            return pRes;
        }

        GraphicsPathSharedPtr graphicsPathFromB2DPolygon( const ::basegfx::B2DPolygon& rPoly, bool bNoLineJoin )
        {
            GraphicsPathSharedPtr 				pRes( new Gdiplus::GraphicsPath() );
            ::std::vector< Gdiplus::PointF > 	aPoints;

            graphicsPathFromB2DPolygon( pRes, aPoints, rPoly, bNoLineJoin );

            return pRes;
        }

        GraphicsPathSharedPtr graphicsPathFromB2DPolyPolygon( const ::basegfx::B2DPolyPolygon& rPoly, bool bNoLineJoin )
        {
            GraphicsPathSharedPtr 				pRes( new Gdiplus::GraphicsPath() );
            ::std::vector< Gdiplus::PointF > 	aPoints;

            const sal_uInt32 nPolies( rPoly.count() );
            for( sal_uInt32 nCurrPoly=0; nCurrPoly<nPolies; ++nCurrPoly )
            {
                graphicsPathFromB2DPolygon( pRes,
                                            aPoints,
                                            rPoly.getB2DPolygon( nCurrPoly ),
                                            bNoLineJoin);
            }

            return pRes;
        }

        GraphicsPathSharedPtr graphicsPathFromXPolyPolygon2D( const uno::Reference< rendering::XPolyPolygon2D >& xPoly, bool bNoLineJoin )
        {
            LinePolyPolygon* pPolyImpl = dynamic_cast< LinePolyPolygon* >( xPoly.get() );

            if( pPolyImpl )
            {
                return pPolyImpl->getGraphicsPath( bNoLineJoin );
            }
            else
            {
                return tools::graphicsPathFromB2DPolyPolygon(
                    polyPolygonFromXPolyPolygon2D( xPoly ), bNoLineJoin );
            }
        }

        bool drawGdiPlusBitmap( const GraphicsSharedPtr& rGraphics,
                                const BitmapSharedPtr&	 rBitmap )
        {
            Gdiplus::PointF	aPoint;
            return (Gdiplus::Ok == rGraphics->DrawImage( rBitmap.get(),
                                                         aPoint ) );
        }

        bool drawDIBits( const GraphicsSharedPtr& rGraphics,
                         const BITMAPINFO&		  rBI,
                         const void*			  pBits )
        {
            BitmapSharedPtr pBitmap(
                Gdiplus::Bitmap::FromBITMAPINFO( &rBI,
                                                 (void*)pBits ) );

            return drawGdiPlusBitmap( rGraphics,
                                      pBitmap );
        }

        bool drawRGBABits( const GraphicsSharedPtr& rGraphics,
                           const RawRGBABitmap&		rRawRGBAData )
        {
            BitmapSharedPtr pBitmap( new Gdiplus::Bitmap( rRawRGBAData.mnWidth,
                                                          rRawRGBAData.mnHeight,
                                                          PixelFormat32bppARGB ) );

            Gdiplus::BitmapData aBmpData;
            aBmpData.Width		 = rRawRGBAData.mnWidth;
            aBmpData.Height		 = rRawRGBAData.mnHeight;
            aBmpData.Stride 	 = 4*aBmpData.Width; // bottom-up format
            aBmpData.PixelFormat = PixelFormat32bppARGB;
            aBmpData.Scan0		 = rRawRGBAData.mpBitmapData.get();

            const Gdiplus::Rect aRect( 0,0,aBmpData.Width,aBmpData.Height );
            if( Gdiplus::Ok != pBitmap->LockBits( &aRect,
                                                  Gdiplus::ImageLockModeWrite | Gdiplus::ImageLockModeUserInputBuf,
                                                  PixelFormat32bppARGB,
                                                  &aBmpData ) )
            {
                return false;
            }

            // commit data to bitmap
            pBitmap->UnlockBits( &aBmpData );

            return drawGdiPlusBitmap( rGraphics,
                                      pBitmap );
        }

        BitmapSharedPtr bitmapFromXBitmap( const uno::Reference< rendering::XBitmap >& xBitmap )
        {
            BitmapProvider* pBitmapProvider = dynamic_cast< BitmapProvider* >(xBitmap.get());

            if( pBitmapProvider )
            {
				IBitmapSharedPtr pBitmap( pBitmapProvider->getBitmap() );
				return pBitmap->getBitmap();
            }
            else
            {
                // not a native CanvasBitmap, extract VCL bitmap and
                // render into GDI+ bitmap of similar size
                // =================================================

                const geometry::IntegerSize2D aBmpSize( xBitmap->getSize() );
                BitmapSharedPtr 			  pBitmap;

                if( xBitmap->hasAlpha() )
                {
                    // TODO(P2): At least for the alpha bitmap case, it
                    // would be possible to generate the corresponding
                    // bitmap directly
                    pBitmap.reset( new Gdiplus::Bitmap( aBmpSize.Width,
                                                        aBmpSize.Height,
                                                        PixelFormat32bppARGB ) );
                }
                else
                {
                    // TODO(F2): Might be wise to create bitmap compatible
                    // to the VCL bitmap. Also, check whether the VCL
                    // bitmap's system handles can be used to create the
                    // GDI+ bitmap (currently, it does not seem so).
                    pBitmap.reset( new Gdiplus::Bitmap( aBmpSize.Width,
                                                        aBmpSize.Height,
                                                        PixelFormat24bppRGB ) );
                }

				GraphicsSharedPtr pGraphics(createGraphicsFromBitmap(pBitmap));
		        tools::setupGraphics(*pGraphics);
                if( !drawVCLBitmapFromXBitmap(
                        pGraphics,
                        xBitmap) )
                {
                    pBitmap.reset();
                }

                return pBitmap;
            }
        }

        CanvasFont::ImplRef	canvasFontFromXFont( const uno::Reference< rendering::XCanvasFont >& xFont )
        {
            CanvasFont* pCanvasFont = dynamic_cast< CanvasFont* >(xFont.get());

            ENSURE_ARG_OR_THROW( pCanvasFont,
                             "canvasFontFromXFont(): Invalid XFont (or incompatible font for this XCanvas)" );

            return CanvasFont::ImplRef( pCanvasFont );
        }

        void setModulateImageAttributes( Gdiplus::ImageAttributes& o_rAttr,
                                         double					   nRedModulation,
                                         double					   nGreenModulation,
                                         double					   nBlueModulation,
                                         double					   nAlphaModulation )
        {
            // This gets rather verbose, but we have to setup a color
            // transformation matrix, in order to incorporate the global
            // alpha value mfAlpha into the bitmap rendering.
            Gdiplus::ColorMatrix	 aColorMatrix;

            aColorMatrix.m[0][0] = static_cast<Gdiplus::REAL>(nRedModulation);
            aColorMatrix.m[0][1] = 0.0;
            aColorMatrix.m[0][2] = 0.0;
            aColorMatrix.m[0][3] = 0.0;
            aColorMatrix.m[0][4] = 0.0;

            aColorMatrix.m[1][0] = 0.0;
            aColorMatrix.m[1][1] = static_cast<Gdiplus::REAL>(nGreenModulation);
            aColorMatrix.m[1][2] = 0.0;
            aColorMatrix.m[1][3] = 0.0;
            aColorMatrix.m[1][4] = 0.0;

            aColorMatrix.m[2][0] = 0.0;
            aColorMatrix.m[2][1] = 0.0;
            aColorMatrix.m[2][2] = static_cast<Gdiplus::REAL>(nBlueModulation);
            aColorMatrix.m[2][3] = 0.0;
            aColorMatrix.m[2][4] = 0.0;

            aColorMatrix.m[3][0] = 0.0;
            aColorMatrix.m[3][1] = 0.0;
            aColorMatrix.m[3][2] = 0.0;
            aColorMatrix.m[3][3] = static_cast<Gdiplus::REAL>(nAlphaModulation);
            aColorMatrix.m[3][4] = 0.0;

            aColorMatrix.m[4][0] = 0.0;
            aColorMatrix.m[4][1] = 0.0;
            aColorMatrix.m[4][2] = 0.0;
            aColorMatrix.m[4][3] = 0.0;
            aColorMatrix.m[4][4] = 1.0;

            o_rAttr.SetColorMatrix( &aColorMatrix,
                                    Gdiplus::ColorMatrixFlagsDefault,
                                    Gdiplus::ColorAdjustTypeDefault );
        }

    } // namespace tools
} // namespace dxcanvas

ext: cxx

max_stars_repo_name: Grosskopf/openoffice

max_issues_repo_name: Grosskopf/openoffice

max_forks_repo_name: Grosskopf/openoffice


text:
#pragma once

// system includes ------------------------------------------------------------
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <memory>
#include <iostream>
#include <vector>

#include "base/exceptions.hpp"
#include "spectral/basis/spectral_basis_factory_ks.hpp"

namespace boltzmann {
class CollisionTensorGalerkin
{
 private:
  typedef Eigen::SparseMatrix<double> sparse_matrix_t;
  typedef std::shared_ptr<sparse_matrix_t> ptr_t;
  typedef Eigen::SparseLU<sparse_matrix_t> lu_t;
  typedef typename SpectralBasisFactoryKS::basis_type basis_t;

  typedef Eigen::VectorXd vec_t;
  typedef Eigen::MatrixXd matrix_t;
  typedef Eigen::Matrix4d m4_t;
  typedef Eigen::DiagonalMatrix<double, -1> diag_t;

 public:
  CollisionTensorGalerkin(const basis_t &basis);
  CollisionTensorGalerkin()
      : N_(0)
  {
  }

  void add(ptr_t &slice, unsigned int j);
  void apply(double *out, const double *in) const;
  void apply(double *out, const double *in, const unsigned int L) const;

  void apply_adaptive(double *out, const double *in, int nmax) const;

  void project(double *out, const double *in) const;

  const sparse_matrix_t &get(int j) const;

  /**
   * @brief read tensor from file
   *
   * @param fname Filename
   */
  void read_hdf5(const char *fname);

 private:
  /// basis size
  unsigned int N_;
  /// tensor entries
  std::vector<ptr_t> slices_;

  matrix_t Ht_;
  m4_t HtHinv_;
  diag_t Sinv_;
};

// ------------------------------------------------------------
inline void
CollisionTensorGalerkin::apply(double *out, const double *in) const
{
  typedef Eigen::Map<const vec_t> cvec_t;  // constant vector
  typedef Eigen::Map<vec_t> mvec_t;        // mutable vector

  mvec_t vout(out, N_);
  cvec_t vin(in, N_);

#pragma omp parallel for
  for (unsigned int i = 0; i < N_; ++i) {
    vout(i) = vin.dot((*slices_[i]) * vin);
  }

  vout = Sinv_ * vout;
}

// ------------------------------------------------------------
inline void
CollisionTensorGalerkin::project(double *out, const double *in) const
{
  typedef Eigen::Map<const vec_t> cvec_t;  // constant vector
  typedef Eigen::Map<vec_t> mvec_t;        // mutable vector

  cvec_t vin(in, N_);
  mvec_t vout(out, N_);

  Eigen::Vector4d lambda = HtHinv_ * Ht_ * (vout - vin);
  vout -= Sinv_ * Ht_.transpose() * lambda;
}

// ------------------------------------------------------------
inline void
CollisionTensorGalerkin::apply_adaptive(double *out, const double *in, int nmax) const
{
  typedef Eigen::Map<const vec_t> cvec_t;  // constant vector
  typedef Eigen::Map<vec_t> mvec_t;        // mutable vector

  mvec_t vout(out, N_);
  cvec_t vin(in, nmax);

#pragma omp parallel for
  for (unsigned int i = 0; i < N_; ++i) {
    vout(i) = vin.dot((*slices_[i]).topLeftCorner(nmax, nmax) * vin);
  }
  vout = Sinv_ * vout;
}

// ----------------------------------------------------------------------
inline void
CollisionTensorGalerkin::apply(double *out, const double *in, const unsigned int L) const
{
  typedef Eigen::Map<Eigen::VectorXd> vec_t;
  typedef Eigen::Map<const Eigen::VectorXd> const_vec_t;

  for (unsigned int j = 0; j < N_; j++) {
    for (unsigned int l = 0; l < L; ++l) {
      const_vec_t vin(in + N_ * l, N_);
      out[N_ * l + j] = vin.dot(*slices_[j] * vin);
    }
  }

  for (unsigned int l = 0; l < L; ++l) {
    vec_t vout(out + N_ * l, N_);
    vout = Sinv_ * vout;
  }
}

}  // end boltzmann


text:
/*
Copyright 2009-2021 Nicolas Colombe

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/

#pragma once

#include <engine/common/world.hpp>
#include <BulletCollision/CollisionDispatch/btGhostObject.h>
#include <boost/container/flat_set.hpp>

#include <math/mathtools.hpp>
#include "physicsdef.hpp"

class	btTriggerGhostObject : public btGhostObject
{
public:

  btTriggerGhostObject();

  virtual ~btTriggerGhostObject();

  void SetChangeCB(std::function<void()> iChangeCallback);

  void addOverlappingObjectInternal(btBroadphaseProxy* otherProxy, btBroadphaseProxy* thisProxy = 0);
  void removeOverlappingObjectInternal(btBroadphaseProxy* otherProxy, btDispatcher* dispatcher, btBroadphaseProxy* thisProxy = 0);

  eXl::SmallVector<eXl::ObjectHandle, 1> StealAdded();
  eXl::SmallVector<eXl::ObjectHandle, 1> StealRemoved();

protected:
  std::function<void()> m_ChangeCallback;
  eXl::UnorderedSet<eXl::ObjectHandle> m_Inside;
  eXl::SmallVector<eXl::ObjectHandle, 1> m_Added;
  eXl::SmallVector<eXl::ObjectHandle, 1> m_Removed;
};

namespace eXl
{
  struct PhysicsSystem_Impl;

  class TriggerManager
  {

  public:

    TriggerManager(PhysicsSystem_Impl& iImpl);

    void AddTrigger(ObjectHandle iObj, TriggerDef const& iDef, TriggerCallbackHandle iCallback, ContactFilterCallback iFilter = ContactFilterCallback());
    void DeleteComponent(ObjectHandle iObj);

    inline void UpdateTransform(ObjectHandle iObj, Matrix4f const& iTrans)
    {
      auto iter = m_ObjectToGhost.find(iObj);
      if (iter != m_ObjectToGhost.end())
      {
        Vector3f const& pos = MathTools::GetPosition(iTrans);
        btTransform newTransform;
        newTransform.getBasis().setValue(
          iTrans.m_Data[0], iTrans.m_Data[1], iTrans.m_Data[2],
          iTrans.m_Data[4], iTrans.m_Data[5], iTrans.m_Data[6],
          iTrans.m_Data[8], iTrans.m_Data[9], iTrans.m_Data[10]);
        newTransform.setOrigin(btVector3(pos.X(), pos.Y(), pos.Z()));
        iter->second->setWorldTransform(newTransform);
      }
    }

    void Tick(World& iWorld, float);

    TriggerCallbackHandle AddCallback(std::unique_ptr<TriggerCallback> iCallback);

  protected:

    struct TriggerEntry
    {
      btTriggerGhostObject triggerObj;
      ContactFilterCallback contactCb;
      ObjectHandle object;
    };

    struct CallbackEntry
    {
      std::unique_ptr<TriggerCallback> m_Callback;
      ObjectTable<TriggerEntry> m_Entries;
      Vector<TriggerCallback::ObjectPair> m_NewEnter;
      Vector<TriggerCallback::ObjectPair> m_NewLeave;

    };

    using TriggerKey = std::pair<ObjectTable<CallbackEntry>::Handle, ObjectTable<TriggerEntry>::Handle>;
    boost::container::flat_set<TriggerKey, std::less<TriggerKey>, Allocator<TriggerKey>> m_Dirty;

    ObjectTable<CallbackEntry> m_Callbacks;
    UnorderedMap<ObjectHandle, TriggerKey> m_ObjectToEntry;
    UnorderedMap<ObjectHandle, btTriggerGhostObject*> m_ObjectToGhost;

    PhysicsSystem_Impl& m_Impl;
  };
}

max_stars_repo_name: eXl-Nic/eXl

max_issues_repo_name: eXl-Nic/eXl

max_forks_repo_name: eXl-Nic/eXl


text:
// Boost.Geometry (aka GGL, Generic Geometry Library)

// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland.

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_GEOMETRY_EXTENSIONS_ALGEBRA_ALGORITHMS_TRANSLATION_HPP
#define BOOST_GEOMETRY_EXTENSIONS_ALGEBRA_ALGORITHMS_TRANSLATION_HPP

#include <boost/geometry/extensions/algebra/geometries/concepts/vector_concept.hpp>

//#include <boost/geometry/geometries/concepts/point_concept.hpp>
//#include <boost/geometry/util/for_each_coordinate.hpp>
#include <boost/geometry/arithmetic/arithmetic.hpp>

namespace boost { namespace geometry
{

// TODO - implement point_operation2 taking two arguments

template <typename Point1, typename Point2, typename Vector>
inline void translation(Point1 const& p1, Point2 const& p2, Vector & v)
{
    concepts::check_concepts_and_equal_dimensions<Point1 const, Point2 const>();
    // TODO - replace the following by check_equal_dimensions
    concepts::check_concepts_and_equal_dimensions<Point1 const, Vector>();

    for_each_coordinate(v, detail::point_assignment<Point2>(p2));
    for_each_coordinate(v, detail::point_operation<Point1, std::minus>(p1));
}

template <typename Vector, typename Point1, typename Point2>
inline Vector return_translation(Point1 const& p1, Point2 const& p2)
{
    Vector v;
    translation(p1, p2, v);
    return v;
}

}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_EXTENSIONS_ALGEBRA_ALGORITHMS_ASSIGN_HPP

max_stars_repo_name: jonasdmentia/geometry

max_issues_repo_name: jonasdmentia/geometry

max_forks_repo_name: jonasdmentia/geometry


text:
/*
 *  SimpleStaircase.cpp
 *  MWorksCore
 *
 *  Created by David Cox on 2/28/07.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "SimpleStaircase.h"
#include "ComponentRegistry.h"
#include <boost/lexical_cast.hpp>

BEGIN_NAMESPACE_MW

SimpleStaircaseOptimizer::SimpleStaircaseOptimizer(shared_ptr<Variable> _watch,
													 shared_ptr<Variable> _output,
													 shared_ptr<Variable> _lower_limit,
													 shared_ptr<Variable> _upper_limit,
													 shared_ptr<Variable> _step_size,
													 shared_ptr<Variable> _up_criterion,
													 shared_ptr<Variable> _down_criterion,
													 int _grace_period,
													 int _history){

	watch = _watch;
	output = _output;
	lower_limit = _lower_limit;
	upper_limit = _upper_limit;
	step_size = _step_size;
	history = _history;
	up_criterion = _up_criterion;
	down_criterion = _down_criterion;
	clear_history_after_step = true;

	if(_grace_period > history) {
		grace_period = _grace_period;
	} else {
		grace_period = _history;
	}

    count = 0;

	registerInput(watch);
	registerOutput(output);

    averager = new BoxcarFilter1D(history, false);

    // Enforce the limits to start
    if((double)(*output) > (double)(*upper_limit)){
        postResults(0, (Datum)((double)*upper_limit));
    }

    if((double)(*output) < (double)(*lower_limit)){
        postResults(0, (Datum)((double)*lower_limit));
    }

}

void SimpleStaircaseOptimizer::stepUp(){
	double stepped_up = (double)(*output) + (double)(*step_size);
	if((double)(*upper_limit) >= stepped_up && (double)(*lower_limit) <= stepped_up){
		postResults(0, (Datum)stepped_up);
	} else {
		if((double)(*upper_limit) < stepped_up) {
			postResults(0, (Datum)((double)*upper_limit));
		} else if((double)(*lower_limit) > stepped_up) {
			postResults(0, (Datum)((double)*lower_limit));
		}
	}

	if(clear_history_after_step){
		averager->reset();
		count = 0;
	}
}

void SimpleStaircaseOptimizer::stepDown(){
	double stepped_down = (double)(*output) - (double)(*step_size);
	if((double)(*upper_limit) >= stepped_down && (double)(*lower_limit) <= stepped_down){
		postResults(0, (Datum)stepped_down);
	} else {
		if((double)(*upper_limit) < stepped_down) {
			postResults(0, (Datum)((double)*upper_limit));
		} else if((double)(*lower_limit) > stepped_down) {
			postResults(0, (Datum)((double)*lower_limit));
		}
	}

	if(clear_history_after_step){
		averager->reset();
		count = 0;
	}
}

void SimpleStaircaseOptimizer::newDataReceived(int inputIndex, const Datum& data,
                                                MWTime timeUS){

	if(inputIndex != 0){
		// error
		return;
	}

	count++;

	averager->input(data);

    if(count < grace_period){
		return;
	}

	double val;
	averager->output(&val);

	if(val >= (double)*up_criterion){
		stepUp();
	} else if (val <= (double)*down_criterion){
		stepDown();
	}
}

void SimpleStaircaseOptimizer::reset(){
	averager->reset();
	count = 0;
}

shared_ptr<mw::Component> SimpleStaircaseOptimizerFactory::createObject(std::map<std::string, std::string> parameters,
															ComponentRegistry *reg) {
	REQUIRE_ATTRIBUTES(parameters,
					   "watch",
					   "output",
					   "lower_limit",
					   "upper_limit",
					   "up_criterion",
					   "down_criterion",
					   "history",
					   "grace_period",
					   "step_size");

	shared_ptr <Variable> watch = reg->getVariable(parameters.find("watch")->second);
	shared_ptr <Variable> output = reg->getVariable(parameters.find("output")->second);
	shared_ptr <Variable> lower_limit = reg->getVariable(parameters.find("lower_limit")->second);
	shared_ptr <Variable> upper_limit = reg->getVariable(parameters.find("upper_limit")->second);
	shared_ptr <Variable> up_criterion = reg->getVariable(parameters.find("up_criterion")->second);
	shared_ptr <Variable> down_criterion = reg->getVariable(parameters.find("down_criterion")->second);
	shared_ptr <Variable> step_size = reg->getVariable(parameters.find("step_size")->second);

	unsigned int grace_period = 0;
	try {
		grace_period = boost::lexical_cast< unsigned int >(parameters.find("grace_period")->second);
	} catch(boost::bad_lexical_cast &) {
		throw InvalidReferenceException(parameters["reference_id"], "grace_period", parameters.find("grace_period")->second);
	}

	unsigned int history = 0;
	try {
		history = boost::lexical_cast< unsigned int >(parameters.find("history")->second);
	} catch(boost::bad_lexical_cast &) {
		throw InvalidReferenceException(parameters["reference_id"], "history", parameters.find("history")->second);
	}

	checkAttribute(watch, parameters["reference_id"], "watch", parameters.find("watch")->second);
	checkAttribute(output, parameters["reference_id"], "output", parameters.find("output")->second);
	checkAttribute(lower_limit, parameters["reference_id"], "lower_limit", parameters.find("lower_limit")->second);
	checkAttribute(upper_limit, parameters["reference_id"], "upper_limit", parameters.find("upper_limit")->second);
	checkAttribute(up_criterion, parameters["reference_id"], "up_criterion", parameters.find("up_criterion")->second);
	checkAttribute(down_criterion, parameters["reference_id"], "down_criterion", parameters.find("down_criterion")->second);
	checkAttribute(step_size, parameters["reference_id"], "step_size", parameters.find("step_size")->second);

	shared_ptr <mw::Component> newSimpleStaircaseOptimizer = shared_ptr<mw::Component>(new SimpleStaircaseOptimizer(watch,
																											   output,
																											   lower_limit,
																											   upper_limit,
																											   step_size,
																											   up_criterion,
																											   down_criterion,
																											   grace_period,
																											   history));
	return newSimpleStaircaseOptimizer;
}

END_NAMESPACE_MW

max_stars_repo_name: esayui/mworks

max_issues_repo_name: esayui/mworks

max_forks_repo_name: esayui/mworks


text:
#ifndef BLUB_LOG_GLOBALLOGGER_HPP
#define BLUB_LOG_GLOBALLOGGER_HPP

#include <boost/log/sources/global_logger_storage.hpp>

#define BLUB_LOG_GLOBAL_LOGGER(name, logger_type) BOOST_LOG_GLOBAL_LOGGER(name, logger_type)
#define BLUB_LOG_GLOBAL_LOGGER_INIT(name, logger_type) BOOST_LOG_GLOBAL_LOGGER_INIT(name, logger_type)

#endif // BLUB_LOG_GLOBALLOGGER_HPP


text:
//  Copyright (c) 2001-2008 Hartmut Kaiser
//
//  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)

#if !defined(SPIRIT_KARMA_CORE_MARCH_06_2007_0833PM)
#define SPIRIT_KARMA_CORE_MARCH_06_2007_0833PM

#if defined(_MSC_VER) && (_MSC_VER >= 1020)
#pragma once      // MS compatible compilers support #pragma once
#endif

#include <boost/spirit/home/karma/char.hpp>
#include <boost/spirit/home/karma/string.hpp>
#include <boost/spirit/home/karma/numeric.hpp>
#include <boost/spirit/home/karma/operator.hpp>
#include <boost/spirit/home/karma/nonterminal.hpp>
#include <boost/spirit/home/karma/action.hpp>
#include <boost/spirit/home/karma/directive.hpp>
#include <boost/spirit/home/karma/auxiliary.hpp>
#include <boost/spirit/home/karma/generate.hpp>
#include <boost/spirit/home/karma/what.hpp>

#endif

max_stars_repo_name: jaredhoberock/gotham

max_issues_repo_name: jaredhoberock/gotham

max_forks_repo_name: jaredhoberock/gotham


text:
#include "mcubes.h"

#include <climits>
#include <algorithm>
#include <unordered_map>

// {{
#include <Eigen/Dense>
// }}

#include "common/array3d.h"
#include "common/progress.h"
#include "mcubes_utils.h"

double getThresholdOtsu(const Volume &volume) {
    // {{ NOT_IMPL_ERROR();
    const double total = volume.size(0) * volume.size(1) * volume.size(2);
    double hist[USHRT_MAX];
    std::memset(hist, 0, sizeof(hist));
    for (uint64_t z = 0; z < volume.size(2); z++) {
        for (uint64_t y = 0; y < volume.size(1); y++) {
            for (uint64_t x = 0; x < volume.size(0); x++) {
                const uint16_t val = volume(x, y, z);
                if (val != 0) {
                    hist[val] += 1;
                }
            }
        }
    }

    double sum1 = 0.0;
    double c1 = 0.0;
    for (int i = 0; i < USHRT_MAX; i++) {
        hist[i] /= total;
        sum1 += hist[i] * i / (double)USHRT_MAX;
        c1 += hist[i];
    }

    double sum2 = 0.0;
    double c2 = 0.0;
    double maxVar = 0.0;
    double threshold = 0.0;
    for (int i = 0; i < USHRT_MAX; i++) {
        const double mu1 = c1 != 0 ? sum1 / c1 : 0.0;
        const double mu2 = c2 != 0 ? sum2 / c2 : 0.0;

        const double diff = mu1 - mu2;
        const double var = c1 * c2 * diff * diff;
        if (maxVar < var) {
            maxVar = var;
            threshold = i / (double)USHRT_MAX;
        }

        sum1 -= hist[i] * i / (double)USHRT_MAX;
        c1 -= hist[i];
        sum2 += hist[i] * i / (double)USHRT_MAX;
        c2 += hist[i];
    }

    return threshold;
    // }}
}

void marchCubes(const Volume &volume, std::vector<Vec3> *vertices, std::vector<uint32_t> *indices, double threshold, bool flipFaces) {
    // Clear arrays
    vertices->clear();
    indices->clear();

    // Compute threshold with Otsu's method, if threshold is not specified.
    if (threshold < 0.0) {
        threshold = getThresholdOtsu(volume);
    }
    printf("Threshold: %.5f\n", threshold);

    // Marching cubes
    GRIDCELL cell;
    TRIANGLE tris[16];
    const Vec3 resolution = Vec3(1.0, 1.0, 1.0);

    // See: "http://paulbourke.net/geometry/polygonise/"
    // Vertex ordering in the cube is rather wired, and
    // the array below re-orders the vertices for easy bit masking.
    static int indexTable[8] = { 0, 1, 4, 5, 3, 2, 7, 6 };

    ProgressBar pbar((volume.size(1) - 1) * (volume.size(2) - 1));
    std::unordered_map<Vec3, uint32_t> uniqueVertices;
    for (uint64_t z = 0; z < volume.size(2) - 1; z++) {
        for (uint64_t y = 0; y < volume.size(1) - 1; y++) {
            for (uint64_t x = 0; x < volume.size(0) - 1; x++) {
                // {{ NOT_IMPL_ERROR();
                for (int i = 0; i < 8; i++) {
                    const int dx = (i >> 0) & 0x01;
                    const int dy = (i >> 1) & 0x01;
                    const int dz = (i >> 2) & 0x01;
                    cell.p[indexTable[i]] = Vec3(x + dx, y + dy, z + dz) * resolution;
                    cell.val[indexTable[i]] = volume(x + dx, y + dy, z + dz) / (double)USHRT_MAX;
                }

                std::fill(tris, tris + 16, TRIANGLE());
                const int ntris = Polygonise(cell, threshold, tris);

                for (int i = 0; i < ntris; i++) {
                    uint32_t tri[3];
                    for (int j = 0; j < 3; j++) {
                        const int k = flipFaces ? 2 - j : j;
                        const Vec3 &v = tris[i].p[k];
                        if (uniqueVertices.count(v) == 0) {
                            uniqueVertices[v] = static_cast<uint32_t>(vertices->size());
                            vertices->push_back(v);
                        }
                        tri[j] = uniqueVertices[v];
                    }

                    if (tri[0] != tri[1] && tri[0] != tri[2] && tri[1] != tri[2]) {
                        indices->push_back(tri[0]);
                        indices->push_back(tri[1]);
                        indices->push_back(tri[2]);
                    }
                }
                // }}
            }
            pbar.step();
        }
    }

    printf("#vert: %d\n", (int)vertices->size());
    printf("#face: %d\n", (int)indices->size() / 3);
}

// {{

void marchTets(const Volume &volume, std::vector<Vec3> *vertices, std::vector<uint32_t> *indices, double threshold, bool flipFaces) {
    // Clear arrays
    vertices->clear();
    indices->clear();

    // Compute threshold with Otsu's method, if threshold is not specified.
    if (threshold < 0.0) {
        threshold = getThresholdOtsu(volume);
    }
    printf("Threshold: %.5f\n", threshold);

    // Marching tetrahedra
    GRIDCELL cell;
    TETRAHEDRON tet;
    TRIANGLE tris[2];
    const Vec3 resolution = Vec3(1.0, 1.0, 1.0);

    // See: "http://paulbourke.net/geometry/polygonise/"
    // Vertex ordering in the cube is rather wired, and
    // the array below re-orders the vertices for easy bit masking.
    static int indexTable[8] = { 0, 1, 4, 5, 3, 2, 7, 6 };
    static const int tetsTable[6][4] = {
            { 6, 0, 5, 1 }, { 6, 0, 4, 5 },
            { 6, 2, 0, 1 }, { 6, 0, 7, 4 },
            { 6, 2, 3, 0 }, { 6, 0, 3, 7 }
    };

    ProgressBar pbar((volume.size(1) - 1) * (volume.size(2) - 1));
    std::unordered_map<Vec3, uint32_t> uniqueVertices;
    for (uint64_t z = 0; z < volume.size(2) - 1; z++) {
        for (uint64_t y = 0; y < volume.size(1) - 1; y++) {
            for (uint64_t x = 0; x < volume.size(0) - 1; x++) {
                for (int i = 0; i < 8; i++) {
                    const int dx = (i >> 0) & 0x01;
                    const int dy = (i >> 1) & 0x01;
                    const int dz = (i >> 2) & 0x01;
                    cell.p[indexTable[i]] = Vec3(x + dx, y + dy, z + dz) * resolution;
                    cell.val[indexTable[i]] = volume(x + dx, y + dy, z + dz) / (double)USHRT_MAX;
                }

                for (int t = 0; t < 6; t++) {
                    for (int j = 0; j < 4; j++) {
                        tet.p[j] = cell.p[tetsTable[t][j]];
                        tet.val[j] = cell.val[tetsTable[t][j]];
                    }

                    std::fill(tris, tris + 16, TRIANGLE());
                    const int ntris = PolygonizeTet(tet, threshold, tris);

                    for (int i = 0; i < ntris; i++) {
                        uint32_t tri[3];
                        for (int j = 0; j < 3; j++) {
                            const int k = flipFaces ? 2 - j : j;
                            const Vec3 &v = tris[i].p[k];
                            if (uniqueVertices.count(v) == 0) {
                                uniqueVertices[v] = static_cast<uint32_t>(vertices->size());
                                vertices->push_back(v);
                            }
                            tri[j] = uniqueVertices[v];
                        }

                        if (tri[0] != tri[1] && tri[0] != tri[2] && tri[1] != tri[2]) {
                            indices->push_back(tri[0]);
                            indices->push_back(tri[1]);
                            indices->push_back(tri[2]);
                        }
                    }
                }
            }
            pbar.step();
        }
    }

    printf("#vert: %d\n", (int)vertices->size());
    printf("#face: %d\n", (int)indices->size() / 3);
}

void dualContour(const Volume &volume, std::vector<Vec3> *vertices, std::vector<uint32_t> *indices, double threshold, bool flipFaces) {
    // Clear arrays
    vertices->clear();
    indices->clear();

    // Compute threshold with Otsu's method, if threshold is not specified.
    if (threshold < 0.0) {
        threshold = getThresholdOtsu(volume);
    }
    printf("Threshold: %.5f\n", threshold);

    // Compute normals
    Array3D<Vec3> normals(volume.size(0), volume.size(1), volume.size(2));
    for (int64_t z = 0; z < (int64_t)volume.size(2); z++) {
        for (int64_t y = 0; y < (int64_t)volume.size(1); y++) {
            for (int64_t x = 0; x < (int64_t)volume.size(0); x++) {
                const int64_t x0 = std::max(0ll, x - 1);
                const int64_t x1 = std::min(x + 1, (int64_t)volume.size(0) - 1);
                const int64_t y0 = std::max(0ll, y - 1);
                const int64_t y1 = std::min(y + 1, (int64_t)volume.size(1) - 1);
                const int64_t z0 = std::max(0ll, z - 1);
                const int64_t z1 = std::min(z + 1, (int64_t)volume.size(2) - 1);
                const double dx = ((volume(x1, y, z) / (double)USHRT_MAX) - (volume(x0, y, z) / (double)USHRT_MAX)) / (double)(x1 - x0);
                const double dy = ((volume(x, y1, z) / (double)USHRT_MAX) - (volume(x, y0, z) / (double)USHRT_MAX)) / (double)(y1 - y0);
                const double dz = ((volume(x, y, z1) / (double)USHRT_MAX) - (volume(x, y, z0) / (double)USHRT_MAX)) / (double)(z1 - z0);
                if (dx != 0.0 || dy != 0.0 || dz != 0.0) {
                    normals(x, y, z) = normalize(Vec3(dx, dy, dz)) * (flipFaces ? -1.0 : 1.0);
                } else {
                    normals(x, y, z) = Vec3(0.0, 0.0, 0.0);
                }
            }
        }
    }

    // Check intersection between cube edges and iso-contours
    Array3D<std::tuple<Vec3, Vec3>> edges[3];

    // X axis
    edges[0] = Array3D<std::tuple<Vec3, Vec3>>(volume.size(0) - 1, volume.size(1), volume.size(2));
    for (int64_t z = 0; z < (int64_t)volume.size(2); z++) {
        for (int64_t y = 0; y < (int64_t)volume.size(1); y++) {
            for (int64_t x = 0; x < (int64_t)volume.size(0) - 1; x++) {
                const double v0 = volume(x, y, z) / (double)USHRT_MAX;
                const double v1 = volume(x + 1, y, z) / (double)USHRT_MAX;
                if ((v0 < threshold) != (v1 < threshold)) {
                    const Vec3 p0 = Vec3(x, y, z);
                    const Vec3 p1 = Vec3(x + 1, y, z);
                    const Vec3 p = VertexInterp(threshold, p0, p1, v0, v1);
                    const Vec3 n0 = normals(x, y, z);
                    const Vec3 n1 = normals(x + 1, y, z);
                    const Vec3 n = normalize(std::abs(v1 - threshold) * n0 + std::abs(v0 - threshold) * n1);
                    edges[0](x, y, z) = std::make_tuple(p, n);
                }
            }
        }
    }

    // Y axis
    edges[1] = Array3D<std::tuple<Vec3, Vec3>>(volume.size(0), volume.size(1) - 1, volume.size(2));
    for (int64_t z = 0; z < (int64_t)volume.size(2); z++) {
        for (int64_t y = 0; y < (int64_t)volume.size(1) - 1; y++) {
            for (int64_t x = 0; x < (int64_t)volume.size(0); x++) {
                const double v0 = volume(x, y, z) / (double)USHRT_MAX;
                const double v1 = volume(x, y + 1, z) / (double)USHRT_MAX;
                if ((v0 < threshold) != (v1 < threshold)) {
                    const Vec3 p0 = Vec3(x, y, z);
                    const Vec3 p1 = Vec3(x, y + 1, z);
                    const Vec3 p = VertexInterp(threshold, p0, p1, v0, v1);
                    const Vec3 n0 = normals(x, y, z);
                    const Vec3 n1 = normals(x, y + 1, z);
                    const Vec3 n = normalize(std::abs(v1 - threshold) * n0 + std::abs(v0 - threshold) * n1);
                    edges[1](x, y, z) = std::make_tuple(p, n);
                }
            }
        }
    }

    // Z axis
    edges[2] = Array3D<std::tuple<Vec3, Vec3>>(volume.size(0), volume.size(1), volume.size(2) - 1);
    for (int64_t z = 0; z < (int64_t)volume.size(2) - 1; z++) {
        for (int64_t y = 0; y < (int64_t)volume.size(1); y++) {
            for (int64_t x = 0; x < (int64_t)volume.size(0); x++) {
                const double v0 = volume(x, y, z) / (double)USHRT_MAX;
                const double v1 = volume(x, y, z + 1) / (double)USHRT_MAX;
                if ((v0 < threshold) != (v1 < threshold)) {
                    const Vec3 p0 = Vec3(x, y, z);
                    const Vec3 p1 = Vec3(x, y, z + 1);
                    const Vec3 p = VertexInterp(threshold, p0, p1, v0, v1);
                    const Vec3 n0 = normals(x, y, z);
                    const Vec3 n1 = normals(x, y, z + 1);
                    const Vec3 n = normalize(std::abs(v1 - threshold) * n0 + std::abs(v0 - threshold) * n1);
                    edges[2](x, y, z) = std::make_tuple(p, n);
                }
            }
        }
    }

    // Define vertex positions of iso-surface
    using EigenMatrix3 = Eigen::Matrix<double, 3, 3>;
    using EigenVector3 = Eigen::Matrix<double, 3, 1>;
    static const int offset0[4] = { 0, 1, 0, 1 };
    static const int offset1[4] = { 0, 0, 1, 1 };
    Array3D<Vec3> cubes(volume.size(0) - 1, volume.size(1) - 1, volume.size(2) - 1);
    for (int64_t z = 0; z < (int64_t)volume.size(2) - 1; z++) {
        for (int64_t y = 0; y < (int64_t)volume.size(1) - 1; y++) {
            for (int64_t x = 0; x < (int64_t)volume.size(0) - 1; x++) {
                EigenMatrix3 AA;
                AA.setZero();
                EigenVector3 bb;
                bb.setZero();
                int count = 0;
                Vec3 avg(0.0);

                // X axis
                for (int k = 0; k < 4; k++) {
                    const int64_t nx = x;
                    const int64_t ny = y + offset0[k];
                    const int64_t nz = z + offset1[k];
                    const Vec3 &p = std::get<0>(edges[0](nx, ny, nz)) - Vec3(x, y, z);
                    const Vec3 &n = std::get<1>(edges[0](nx, ny, nz));
                    if (length(p) != 0.0 && length(n) != 0.0) {
                        EigenVector3 nn;
                        nn << n.x, n.y, n.z;
                        EigenVector3 pp;
                        pp << p.x, p.y, p.z;

                        const EigenMatrix3 M = nn * nn.transpose();
                        AA += M;
                        bb += M * pp;
                        avg += p;
                        count += 1;
                    }
                }

                // Y axis
                for (int k = 0; k < 4; k++) {
                    const int64_t nx = x + offset1[k];
                    const int64_t ny = y;
                    const int64_t nz = z + offset0[k];
                    const Vec3 &p = std::get<0>(edges[1](nx, ny, nz)) - Vec3(x, y, z);
                    const Vec3 &n = std::get<1>(edges[1](nx, ny, nz));
                    if (length(p) != 0.0 && length(n) != 0.0) {
                        EigenVector3 nn;
                        nn << n.x, n.y, n.z;
                        EigenVector3 pp;
                        pp << p.x, p.y, p.z;

                        const EigenMatrix3 M = nn * nn.transpose();
                        AA += M;
                        bb += M * pp;
                        avg += p;
                        count += 1;
                    }
                }

                // Z axis
                for (int k = 0; k < 4; k++) {
                    const int64_t nx = x + offset0[k];
                    const int64_t ny = y + offset1[k];
                    const int64_t nz = z;
                    const Vec3 &p = std::get<0>(edges[2](nx, ny, nz)) - Vec3(x, y, z);
                    const Vec3 &n = std::get<1>(edges[2](nx, ny, nz));
                    if (length(p) != 0.0 && length(n) != 0.0) {
                        EigenVector3 nn;
                        nn << n.x, n.y, n.z;
                        EigenVector3 pp;
                        pp << p.x, p.y, p.z;

                        const EigenMatrix3 M = nn * nn.transpose();
                        AA += M;
                        bb += M * pp;
                        avg += p;
                        count += 1;
                    }
                }

                if (count != 0) {
                    Eigen::JacobiSVD<EigenMatrix3> svd;
                    svd.compute(AA, Eigen::ComputeFullU | Eigen::ComputeFullV);
                    double det = 1.0;
                    auto values = svd.singularValues();
                    for (int dim = 0; dim < values.size(); dim++) {
                        const double v = std::abs(values(dim)) < 0.1 ? 0.0 : values(dim);
                        det *= v;
                        values(dim) = v == 0.0 ? 0.0 : 1.0 / v;
                    }
                    const EigenMatrix3 AAinv = svd.matrixU() * values.asDiagonal() * svd.matrixV().transpose();

                    EigenVector3 xx;
                    if (det != 0.0) {
                        xx = AAinv * bb;
                    } else {
                        avg = avg / count;
                        xx << avg.x, avg.y, avg.z;
                    }

//                    EigenVector3 cc;
//                    cc << avg.x, avg.y, avg.z;

//                    const double beta = 1.0e-2;
//                    const EigenVector3 xx = (AA + beta * II).ldlt().solve(bb + beta * cc);

                    const double cx = x + xx(0); //std::max(0.0, std::min(xx(0), 1.0));
                    const double cy = y + xx(1); //std::max(0.0, std::min(xx(1), 1.0));
                    const double cz = z + xx(2); //std::max(0.0, std::min(xx(2), 1.0));
                    cubes(x, y, z) = Vec3(cx, cy, cz);
                } else {
                    cubes(x, y, z) =  Vec3(x + 0.5, y + 0.5, z + 0.5);
                }
            }
        }
    }

    // Dual contouring
    Vec3 rectangle[4];
    int triindex[2][3] = { {0, 1, 3}, {0, 3, 2} };
    std::unordered_map<Vec3, uint32_t> uniqueVertices;

    const int64_t totalSteps = 3 * (volume.size(1) - 1) * (volume.size(2) - 1);
    ProgressBar pbar(totalSteps);

    // X axis
    for (int64_t z = 1; z < (int64_t)volume.size(2); z++) {
        for (int64_t y = 1; y < (int64_t)volume.size(1); y++) {
            for (int64_t x = 0; x < (int64_t)volume.size(0) - 1; x++) {
                const double v0 = volume(x, y, z) / (double)USHRT_MAX;
                const double v1 = volume(x + 1, y, z) / (double)USHRT_MAX;
                if ((v0 < threshold) != (v1 < threshold)) {
                    rectangle[0] = cubes(x, y - 1, z - 1);
                    rectangle[1] = cubes(x, y, z - 1);
                    rectangle[2] = cubes(x, y - 1, z);
                    rectangle[3] = cubes(x, y, z);

                    for (int t = 0; t < 2; t++) {
                        uint32_t tri[3];
                        for (int k = 0; k < 3; k++) {
                            const Vec3 &v = (v0 > v1) ? rectangle[triindex[t][k]] : rectangle[triindex[t][3 - k - 1]];
                            if (uniqueVertices.count(v) == 0) {
                                uniqueVertices[v] = static_cast<uint32_t>(vertices->size());
                                vertices->push_back(v);
                            }
                            tri[k] = uniqueVertices[v];
                        }

                        if (tri[0] != tri[1] && tri[0] != tri[2] && tri[1] != tri[2]) {
                            indices->push_back(tri[0]);
                            indices->push_back(tri[1]);
                            indices->push_back(tri[2]);
                        }
                    }
                }
            }
            pbar.step();
        }
    }

    // Y axis
    for (int64_t z = 1; z < (int64_t)volume.size(2); z++) {
        for (int64_t y = 0; y < (int64_t)volume.size(1) - 1; y++) {
            for (int64_t x = 1; x < (int64_t)volume.size(0); x++) {
                const double v0 = volume(x, y, z) / (double)USHRT_MAX;
                const double v1 = volume(x, y + 1, z) / (double)USHRT_MAX;
                if ((v0 < threshold) != (v1 < threshold)) {
                    rectangle[0] = cubes(x - 1, y, z - 1);
                    rectangle[1] = cubes(x - 1, y, z);
                    rectangle[2] = cubes(x, y,  z - 1);
                    rectangle[3] = cubes(x, y, z);

                    for (int t = 0; t < 2; t++) {
                        uint32_t tri[3];
                        for (int k = 0; k < 3; k++) {
                            const Vec3 &v = (v0 > v1) ? rectangle[triindex[t][k]] : rectangle[triindex[t][3 - k - 1]];
                            if (uniqueVertices.count(v) == 0) {
                                uniqueVertices[v] = static_cast<uint32_t>(vertices->size());
                                vertices->push_back(v);
                            }
                            tri[k] = uniqueVertices[v];
                        }

                        if (tri[0] != tri[1] && tri[0] != tri[2] && tri[1] != tri[2]) {
                            indices->push_back(tri[0]);
                            indices->push_back(tri[1]);
                            indices->push_back(tri[2]);
                        }
                    }
                }
            }
            pbar.step();
        }
    }

    // Z axis
    for (int64_t z = 0; z < (int64_t)volume.size(2) - 1; z++) {
        for (int64_t y = 1; y < (int64_t)volume.size(1); y++) {
            for (int64_t x = 1; x < (int64_t)volume.size(0); x++) {
                const double v0 = volume(x, y, z) / (double)USHRT_MAX;
                const double v1 = volume(x, y, z + 1) / (double)USHRT_MAX;
                if ((v0 < threshold) != (v1 < threshold)) {
                    rectangle[0] = cubes(x - 1, y - 1, z);
                    rectangle[1] = cubes(x, y - 1, z);
                    rectangle[2] = cubes(x - 1, y, z);
                    rectangle[3] = cubes(x, y, z);

                    for (int t = 0; t < 2; t++) {
                        uint32_t tri[3];
                        for (int k = 0; k < 3; k++) {
                            const Vec3 &v = (v0 > v1) ? rectangle[triindex[t][k]] : rectangle[triindex[t][3 - k - 1]];
                            if (uniqueVertices.count(v) == 0) {
                                uniqueVertices[v] = static_cast<uint32_t>(vertices->size());
                                vertices->push_back(v);
                            }
                            tri[k] = uniqueVertices[v];
                        }

                        if (tri[0] != tri[1] && tri[0] != tri[2] && tri[1] != tri[2]) {
                            indices->push_back(tri[0]);
                            indices->push_back(tri[1]);
                            indices->push_back(tri[2]);
                        }
                    }
                }
            }
            pbar.step();
        }
    }

    printf("#vert: %d\n", (int)vertices->size());
    printf("#face: %d\n", (int)indices->size() / 3);
}

// }}

max_stars_repo_name: tatsy/cpp-python-beginners

max_issues_repo_name: tatsy/cpp-python-beginners

max_forks_repo_name: tatsy/cpp-python-beginners