text:
/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2012 by Wenzel Jakob and others.

    Mitsuba is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License Version 3
    as published by the Free Software Foundation.

    Mitsuba is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#include <mitsuba/render/scene.h>
#include <mitsuba/render/renderproc.h>
#include <mitsuba/core/autodiff.h>
#include <mitsuba/core/statistics.h>
#include <boost/algorithm/string.hpp>

DECLARE_DIFFSCALAR_BASE();

MTS_NAMESPACE_BEGIN

static StatsCounter statsConverged(
		" manifold", "Converged manifold walks", EPercentage);

/*!\plugin{motion}{Motion and specular motion vector integrator}
 * \parameters{
 *     \parameter{time}{\Float}{
 *       Denotes the time stamp of the target frame of the motion vectors.
 *       The current frame is specified via the sensor's \code{shutterOpen}
 *       and \code{shutterClose} parameters, which should both be set to
 *       the same value. \default{0}
 *     }
 *     \parameter{time}{\String}{
 *        Path configuration of the desired motion vectors:
 *        \begin{enumerate}[(i)]
 *           \item \textbf{d}: Primary (non-specular) hit points\vspace{-1mm}
 *           \item \textbf{rd}: A non-specular surface seen through a specular reflection\vspace{-1mm}
 *           \item \textbf{ttd}: A non-specular surface seen through a pair of specular refractions\vspace{-1mm}
 *           \item \textbf{trtd}: A non-specular surface seen through a sequence of refraction, reflection, and refraction events.\vspace{-1mm}
 *        \end{enumerate}
 *        etc.
 *     }
 *     \parameter{derivativesOnly}{\Boolean}{
 *       By default, the Manifold Exploration technique is used to accurately
 *       solve the underlying specular flow problem. When this parameter is set to \code{true},
 *       the nonlinear solver is deactivated, and only first-order extrapolations are provided.
 *       \default{\code{false}}
 *     }
 *     \parameter{glossyThreshold}{\Float}{
 *        Threshold on the roughness parameter of reflectance models
 *        to be classified as specular.
 *        \default{\texttt{0}, i.e.~only perfectly specular materials are classified as specular}
 *     }
 *     \parameter{maxTimeSteps}{\Integer}{
 *        Maximum number of temporal sub-steps \default{5}
 *     }
 *     \parameter{maxSpaceSteps}{\Integer}{
 *        Maximum number of spatial sub-steps \default{10}
 *     }
 * }
 * This integrator extracts motion vectors for animated input scenes, alternatively
 * at primary hit points or at hit points observed through sequences of reflective
 * and refractive objects. The first two color components (R and G) of the resulting
 * rendering specify the screen-space motion in 2D pixel coordinates, and the last
 * component (B) denotes the change in distance of the observed 3D point to the
 * camera position. Sometimes a specular path cannot be tracked from one frame to the other,
 * e.g. because it does not exist, or because the solver did not converge. In this case,
 * the pixel color is set to infinity. The images on the following page show
 * motion vectors obtained for a sphere that is moving from the left to the right.
 *
 * \renderings{
 *     \rendering{Input scene at time $t=0$}{integrator_motion_sphere_1}
 *     \rendering{Input scene at time $t=1$}{integrator_motion_sphere_2}
 * }
 * \renderings{
 *     \medrendering{\code{config="d"}}{integrator_motion_path_d}
 *     \medrendering{\code{config="rd"}}{integrator_motion_path_rd}
 *     \medrendering{\code{config="ttd"}}{integrator_motion_path_ttd}
 * }
 * \renderings{
 *     \rendering{\code{config="trtd"}}{integrator_motion_path_trtd}
 *     \rendering{\code{config="trrtd"}}{integrator_motion_path_trrtd}
 * }
 * \clearpage
 *
 * \begin{xml}[caption={Exemplary scene configuration for computing specular motion vectors}]
 * <scene>
 *     <integrator type="motion">
 *         <string name="config" value="ttd"/>
 *         <float name="time" value="1"/>
 *     </integrator>
 *
 *     <shape type="serialized">
 *         <string name="filename" value="..."/>
 *         <animation name="toWorld">
 *             <transform time="0">
 *                 <!-- Transformation at time zero -->
 *             </transform>
 *             <transform time="1">
 *                 <!-- Transformation at time one -->
 *             </transform>
 *         </animation>
 *          <bsdf type="dielectric"/>
 *     </shape>
 *
 *     <sensor type="perspective">
 *         <float name="shutterOpen" value="0"/>
 *         <float name="shutterClose" value="0"/>
 *
 *         <sampler type="ldsampler">
 *             <integer name="sampleCount" value="1"/>
 *             <boolean name="pixelCenters" value="true"/>
 *         </sampler>
 *
 *         <film type="hdrfilm" id="film">
 *             <string name="pixelFormat" value="rgb"/>
 *             <boolean name="banner" value="false"/>
 *             <rfilter type="box"/>
 *         </film>
 *     </sensor>
 * </scene>
 * \end{xml}
 */

class MotionIntegrator : public SamplingIntegrator {
public:
	typedef Eigen::Matrix<Float, Eigen::Dynamic, Eigen::Dynamic> EMatrix;
	typedef Eigen::Matrix<Float, Eigen::Dynamic, 1> EVector;
	typedef Eigen::Matrix<Float, 1, 7> Gradient;
	typedef DScalar1<Float, Gradient> DScalar;
	typedef DScalar::DVector3 DVector;

	MotionIntegrator(const Properties &props) : SamplingIntegrator(props) {
		m_time = props.getFloat("time");
		m_config = boost::to_lower_copy(props.getString("config", "d"));
		if (m_config.length() == 0)
			Log(EError, "Path configuration string must have at least one entry!");
		if (m_config[m_config.length()-1] != 'd')
			Log(EError, "Configuration string must end with 'd'!");

		m_derivativesOnly = props.getBoolean("derivativesOnly", false);
		m_maxTimeSteps = props.getInteger("maxTimeSteps", 5);
		m_maxSpaceSteps = props.getInteger("maxSpaceSteps", 10);
		m_glossyThreshold = props.getFloat("glossyThreshold", 0);
		m_subSteps = props.getInteger("subSteps", 1);
	}

	MotionIntegrator(Stream *stream, InstanceManager *manager)
	 : SamplingIntegrator(stream, manager) {
		 m_time = stream->readFloat();
		 m_config = stream->readString();
		 m_derivativesOnly = stream->readBool();
		 m_maxTimeSteps = stream->readInt();
		 m_maxSpaceSteps = stream->readInt();
		 m_glossyThreshold = stream->readFloat();
		 m_subSteps = stream->readInt();
	}

	void serialize(Stream *stream, InstanceManager *manager) const {
		SamplingIntegrator::serialize(stream, manager);
		stream->writeFloat(m_time);
		stream->writeString(m_config);
		stream->writeBool(m_derivativesOnly);
		stream->writeInt(m_maxTimeSteps);
		stream->writeInt(m_maxSpaceSteps);
		stream->writeFloat(m_glossyThreshold);
		stream->writeInt(m_subSteps);
	}

	Spectrum Li(const RayDifferential &r, RadianceQueryRecord &rRec) const {
		const Point2 apertureSample(0.5f);
		Point p0, p1;

		if (m_config.length() == 1) {
			Intersection &its = rRec.its;

			// compute motion of environment
			if (!rRec.rayIntersect(r)) {
				BSphere sphere = rRec.scene->getBSphere();
				Float nearT, farT;
				sphere.radius = 1e6;
				sphere.rayIntersect(r, nearT, farT);
				p0 = r(nearT);
			} else {
				p0 = its.p;
			}

			p1 = p0;

			// reproject point according to the motion of the intersected shape
			if (its.isValid() && its.instance != NULL)  {
				Intersection its2(its);
				its2.adjustTime(m_time);
				p1 = its2.p;
			}
		} else {
			std::vector<Intersection> source, target, temp, temp2;
			Ray ray(r);

			/* Trace an initial light path with the given configuration*/
			if (!tracePath(rRec, ray, source))
				return Spectrum(0.0f);

			m_scene = rRec.scene;
			p0 = source[1].p;

			int timeIteration = 0;

			Float stepSizeReduction = 1.0f;
			while (true) {
				if (++timeIteration > m_maxTimeSteps)
					return Spectrum(std::numeric_limits<Float>::infinity());

				Float maxStepSize = (m_time-r.time) / m_subSteps;
				Float timeStepSize = std::min((Float) 1.0f, maxStepSize / (m_time-source[0].time));
				timeStepSize *= stepSizeReduction;

				/* Compute updated intersection records for time 'm_time' */
				adjustTime(rRec, apertureSample, source, target, timeStepSize);

				if (timeIteration == 1) {
					bool moved = false;
					for (size_t i=0; i<source.size(); ++i)
						if ((source[i].p-target[i].p).length() > 1e-4f)
							moved = true;

					if (!moved)
						return Spectrum(0.0f);

					if (m_derivativesOnly) {
						p1 = extrapolateTimePoint(source, target);
						break;
					} else {
						statsConverged.incrementBase();
					}
				}

				if (!timeStep(rRec, source, target, temp, temp2)) {
					stepSizeReduction *= 0.5f;
				} else {
					p1 = source[1].p;
					stepSizeReduction = std::min((Float) 1.0f, stepSizeReduction * 2);

					if (std::abs(source[0].time-m_time) < 1e-5f) {
						++statsConverged;
						break;
					}
				}
			}
		}

		const Sensor *sensor = rRec.scene->getSensor();
		DirectSamplingRecord dRec0(p0, r.time), dRec1(p1, m_time);
		sensor->sampleDirect(dRec0, apertureSample);
		sensor->sampleDirect(dRec1, apertureSample);

		/* Step 4: Compute depth difference */
		Float dDelta = dRec1.dist - dRec0.dist;
		Spectrum result(0.0f);
		result.fromLinearRGB(dRec1.uv.x-dRec0.uv.x,
				dRec1.uv.y-dRec0.uv.y,
				std::isfinite(dDelta) ? dDelta : (Float) 0);
		return result;
	}

	bool timeStep(RadianceQueryRecord &rRec, std::vector<Intersection> &source, const std::vector<Intersection> &target, std::vector<Intersection> &temp, std::vector<Intersection> &temp2) const {
		Ray ray = extrapolateTimeRay(source, target);

		if (!tracePath(rRec, ray, temp))
			return false;

		Float error = computeError(temp, target);

		Float spaceStepSize = 1.0f;
		int spaceIteration = 0;
		while (error > 1e-5f) {
			++spaceIteration;

			if (spaceIteration > m_maxSpaceSteps) {
				return false;
			}

			Ray candidateRay = extrapolateSpaceRay(temp, target, spaceStepSize);

			Float candidateError = 0;
			if (!tracePath(rRec, candidateRay, temp2))
				candidateError = std::numeric_limits<Float>::infinity();
			else
				candidateError = computeError(temp2, target);

			if (candidateError < error) {
				temp = temp2;
				error = candidateError;
				spaceStepSize = std::min((Float) 1.0f, spaceStepSize * 2);
			} else {
				spaceStepSize *= 0.5f;
			}
		}

		source = temp;
		return true;
	}

	bool tracePath(RadianceQueryRecord &rRec, Ray ray, std::vector<Intersection> &intersections) const {
		int depth = 0;

		Intersection its;
		memset(&its, 0, sizeof(Intersection));
		its.p = ray.o;
		its.time = ray.time;

		intersections.clear();
		intersections.push_back(its);

		while (depth < (int) m_config.size()) {
			rRec.scene->rayIntersect(ray, its);

			char interactionType = m_config[depth++];
			if (interactionType == 'd') {
				if (!its.isValid()) {
					BSphere sphere = rRec.scene->getBSphere();
					Float nearT, farT;
					sphere.radius *= 1000;
					bool success = sphere.rayIntersect(ray, nearT, farT);
					Assert(success && nearT < 0 && farT > 0);
					its.p = ray(farT);
					its.geoFrame = its.shFrame = Frame(-ray.d);
					its.dpdu = its.geoFrame.s;
					its.dpdv = its.geoFrame.t;
					its.time = ray.time;
					its.uv = Point2(0.0f);
					its.shape = its.instance = NULL;
				}

				if (its.isValid() && !(its.shape->getBSDF()->getType() & BSDF::EDiffuseReflection) && !its.isEmitter())
					return false;

				intersections.push_back(its);

				break;
			}

			if (!its.isValid())
				return false;

			intersections.push_back(its);

			/* Sample BSDF * cos(theta) */
			const BSDF *bsdf = its.shape->getBSDF();

			BSDFSamplingRecord bRec(its, rRec.sampler, ERadiance);

			if (interactionType == 'r') {
				if (bsdf->getType() & BSDF::EDeltaReflection) {
					bRec.typeMask = BSDF::EDeltaReflection;
				} else if (bsdf->getType() & BSDF::EGlossyReflection) {
					for (int i=0; i<bsdf->getComponentCount(); ++i)
						if ((bsdf->getType(i) & BSDF::EGlossyReflection) && bsdf->getRoughness(its, i) < m_glossyThreshold)
							bRec.typeMask = BSDF::EGlossyReflection;
				}
			} else if (interactionType == 't') {
				if (bsdf->getType() & BSDF::EDeltaTransmission) {
					bRec.typeMask = BSDF::EDeltaTransmission;
				} else if (bsdf->getType() & BSDF::EGlossyTransmission) {
					for (int i=0; i<bsdf->getComponentCount(); ++i)
						if ((bsdf->getType(i) & BSDF::EGlossyTransmission) && bsdf->getRoughness(its, i) < m_glossyThreshold)
							bRec.typeMask = BSDF::EGlossyTransmission;
				}
			}

			if (bRec.typeMask == BSDF::EAll)
				return false;

			Spectrum bsdfWeight = bsdf->sample(bRec, Point2(0.0f));
			if (bsdfWeight.isZero())
				return false;

			const Vector wo = its.toWorld(bRec.wo);

			ray = Ray(its.p, wo, ray.time);
		}

		return true;
	}

	void adjustTime(const RadianceQueryRecord &rRec, const Point2 &apertureSample, const std::vector<Intersection> &source, std::vector<Intersection> &target, Float timeStepSize) const {
		target = source;

		Float targetTime = (1-timeStepSize) * source[0].time + timeStepSize * m_time;

		const Sensor *sensor = rRec.scene->getSensor();
		DirectSamplingRecord dRec0(Point(0.0f), source[0].time),
							 dRec1(Point(0.0f), targetTime);
		sensor->sampleDirect(dRec0, apertureSample);
		sensor->sampleDirect(dRec1, apertureSample);
		Assert((source[0].p-dRec0.p).lengthSquared() < 1e-6);

		target[0].p = dRec1.p;
		target[0].time = targetTime;

		for (size_t i=1; i<target.size(); ++i)
			target[i].adjustTime(targetTime);
	}

	DVector getVertexPosition(const std::vector<Intersection> &source, const std::vector<Intersection> &target, int i, int rel) const {
		DScalar u(rel*2, 0), v(rel*2+1, 0), time(6, 0);

		return DScalar::vector(source[i].p)
		     + DScalar::vector(source[i].dpdu) * u
		     + DScalar::vector(source[i].dpdv) * v
		     + DScalar::vector(target[i].p-source[i].p) * time;
	}

	void getVertexFrame(const std::vector<Intersection> &source, const std::vector<Intersection> &target, int i, DVector &s, DVector &t, DVector &n) const {
		DScalar u(2, 0), v(3, 0), time(6, 0);
		const BSDF *bsdf = source[i].shape->getBSDF();

		Frame frame, du, dv;
		frame = bsdf->getFrame(source[i]);
		bsdf->getFrameDerivative(source[i], du, dv);

		DVector n0 = DScalar::vector(frame.n)
			+ DScalar::vector(du.n) * u
			+ DScalar::vector(dv.n) * v;
		DVector s0 = DScalar::vector(frame.s)
			+ DScalar::vector(du.s) * u
			+ DScalar::vector(dv.s) * v;
		DVector t0 = DScalar::vector(frame.t)
			+ DScalar::vector(du.t) * u
			+ DScalar::vector(dv.t) * v;

		frame = bsdf->getFrame(target[i]);
		bsdf->getFrameDerivative(target[i], du, dv);

		DVector n1 = DScalar::vector(frame.n)
			+ DScalar::vector(du.n) * u
			+ DScalar::vector(dv.n) * v;
		DVector s1 = DScalar::vector(frame.s)
			+ DScalar::vector(du.s) * u
			+ DScalar::vector(dv.s) * v;
		DVector t1 = DScalar::vector(frame.t)
			+ DScalar::vector(du.t) * u
			+ DScalar::vector(dv.t) * v;

		n = (1-time)*n0 + time*n1;
		s = (1-time)*s0 + time*s1;
		t = (1-time)*t0 + time*t1;
	}

	void assembleMatrix(const std::vector<Intersection> &source, const std::vector<Intersection> &target, EMatrix &M) const {
		DScalar::setVariableCount(7);
		M.resize(2*(source.size()-2), 2*source.size()+1);
		M.setZero();

		for (int i=1; i < (int) source.size()-1; ++i) {
			DVector p_pred = getVertexPosition(source, target, i-1, 0);
			DVector p_cur  = getVertexPosition(source, target, i, 1);
			DVector p_succ = getVertexPosition(source, target, i+1, 2);

			DVector s, t, n;
			getVertexFrame(source, target, i, s, t, n);

			DVector wi = p_pred - p_cur, wo = p_succ - p_cur;
			wi *= inverse(sqrt(wi.dot(wi)));
			wo *= inverse(sqrt(wo.dot(wo)));

			Float eta = source[i].shape->getBSDF()->getEta();

			if (m_config[i-1] == 'r')
				eta = 1;
			else if (wi.dot(n).getValue() < 0)
				eta = 1/eta;

			DVector H = wi + wo * DScalar(eta);
			H *= inverse(sqrt(H.dot(H)));

			Gradient H_s = H.dot(s).getGradient(),
			         H_t = H.dot(t).getGradient();

			M.block<1, 6>((i-1)*2+0, (i-1)*2) = H_s.segment<6>(0);
			M.block<1, 6>((i-1)*2+1, (i-1)*2) = H_t.segment<6>(0);
			M((i-1)*2+0, M.cols()-1) = H_s(6);
			M((i-1)*2+1, M.cols()-1) = H_t(6);
		}
	}

	Point extrapolateTimePoint(const std::vector<Intersection> &source, const std::vector<Intersection> &target) const {
		EMatrix M;
		assembleMatrix(source, target, M);
		EVector b = -M.block(0, 2, (source.size()-2)*2, (source.size()-2)*2).lu().solve(M.col(M.cols()-1));
		return target[1].p + target[1].dpdu * b[0] + target[1].dpdv * b[1];
	}

	Ray extrapolateTimeRay(const std::vector<Intersection> &source, const std::vector<Intersection> &target) const {
		EMatrix M;
		assembleMatrix(source, target, M);

		EVector b = -M.block(0, 2, (source.size()-2)*2, (source.size()-2)*2).lu().solve(M.col(M.cols()-1));
		Point rayOrigin = target[0].p;
		Point rayTarget = target[1].p + target[1].dpdu * b[0] + target[1].dpdv * b[1];
		return Ray(rayOrigin, normalize(rayTarget-rayOrigin), target[1].time);
	}

	Float computeError(const std::vector<Intersection> &source, const std::vector<Intersection> &target) const {
		int last = source.size()-1;
		Float scale = std::max(Epsilon,std::max(std::abs(target[last].p.x), std::max(std::abs(target[last].p.y), std::abs(target[last].p.z))));
		return (target[last].p-source[last].p).length() / scale;
	}

	Ray extrapolateSpaceRay(const std::vector<Intersection> &source, const std::vector<Intersection> &target, Float stepSize) const {
		EMatrix M;
		assembleMatrix(source, target, M);

		Eigen::PartialPivLU<EMatrix> lu = M.block(0, 2, (source.size()-2)*2, (source.size()-2)*2).lu();

		int last = source.size()-1;
		Vector rel = target[last].p-source[last].p,
		       dpdu = source[last].dpdu,
			   dpdv = source[last].dpdv;

		Float b1 = dot(rel, dpdu),
			  b2 = dot(rel, dpdv),
			  a11 = dot(dpdu, dpdu), a12 = dot(dpdu, dpdv),
			  a22 = dot(dpdv, dpdv),
			  det = a11 * a22 - a12 * a12;

		Float invDet = 1.0f / det,
		      du = ( a22 * b1 - a12 * b2) * invDet,
		      dv = (-a12 * b1 + a11 * b2) * invDet;

		EVector b = -(du*lu.solve(M.col(M.cols()-3)) + dv*lu.solve(M.col(M.cols()-2)));

		Point rayTarget = source[1].p + stepSize * (source[1].dpdu * b[0] + source[1].dpdv * b[1]);

		return Ray(source[0].p, normalize(rayTarget-source[0].p), source[1].time);
	}

	std::string toString() const {
		return "MotionIntegrator[]";
	}

	MTS_DECLARE_CLASS()
private:
	Float m_time;
	std::string m_config;
	bool m_derivativesOnly;
	int m_maxSpaceSteps;
	int m_maxTimeSteps;
	int m_subSteps;
	Float m_glossyThreshold;
	mutable const Scene *m_scene;
};

MTS_IMPLEMENT_CLASS_S(MotionIntegrator, false, SamplingIntegrator)
MTS_EXPORT_PLUGIN(MotionIntegrator, " motion vector integrator");
MTS_NAMESPACE_END

max_stars_repo_name: NTForked-ML/pbrs

max_issues_repo_name: NTForked-ML/pbrs

max_forks_repo_name: NTForked-ML/pbrs


text:
/**
 *  Software License Agreement (BSD License)
 *
 *  Point Cloud Library (PCL) - www.pointclouds.org
 *  Copyright (c) 2010-2012, Willow Garage, Inc.
 *
 *  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 Willow Garage, Inc. nor the names of its
 *     contributors may be used to endorse or promote products derived
 *     from this software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 *  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 *  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 *  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 *  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 *  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 *  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 *  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 * $Id: $
 * @brief This file is the execution node of the Human Tracking
 * @copyright Copyright (2011) Willow Garage
 * @authors Koen Buys, Anatoly Baksheev
 **/
#include <pcl/pcl_base.h>
#include <pcl/point_types.h>
#include <pcl/point_cloud.h>
#include <pcl/common/time.h>
#include <pcl/exceptions.h>
#include <pcl/console/parse.h>
#include <pcl/console/print.h>
#include <pcl/gpu/containers/initialization.h>
#include <pcl/gpu/people/people_detector.h>
#include <pcl/gpu/people/colormap.h>
#include <pcl/visualization/image_viewer.h>
#include <pcl/io/openni_grabber.h>
#include <pcl/io/oni_grabber.h>
#include <pcl/io/pcd_grabber.h>
#include <pcl/io/pcd_io.h>
#include <pcl/io/png_io.h>
#include <boost/filesystem.hpp>

#include <iostream>

namespace pc = pcl::console;
using namespace pcl::visualization;
using namespace pcl::gpu;
using namespace pcl;
using namespace std;

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////

vector<string> getPcdFilesInDir(const string &directory) {
    namespace fs = boost::filesystem;
    fs::path dir(directory);

    if (!fs::exists(dir) || !fs::is_directory(dir))
        PCL_THROW_EXCEPTION(pcl::IOException, "Wrong PCD directory");

    vector<string> result;
    fs::directory_iterator pos(dir);
    fs::directory_iterator end;

    for (; pos != end; ++pos)
        if (fs::is_regular_file(pos->status()))
            if (fs::extension(*pos) == ".pcd")
                result.push_back(pos->path().string());

    return result;
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////

struct SampledScopeTime : public StopWatch {
    enum { EACH = 33 };
    SampledScopeTime(int &time_ms) : time_ms_(time_ms) {}
    ~SampledScopeTime() {
        static int i_ = 0;
        time_ms_ += getTime();
        if (i_ % EACH == 0 && i_) {
            cout << "Average frame time = " << time_ms_ / EACH << "ms ( "
                 << 1000.f * EACH / time_ms_ << "fps )" << endl;
            time_ms_ = 0;
        }
        ++i_;
    }

  private:
    int &time_ms_;
};

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////

string make_name(int counter, const char *suffix) {
    char buf[4096];
    sprintf(buf, "./people%04d_%s.png", counter, suffix);
    return buf;
}

template <typename T>
void savePNGFile(const std::string &filename,
                 const pcl::gpu::DeviceArray2D<T> &arr) {
    int c;
    pcl::PointCloud<T> cloud(arr.cols(), arr.rows());
    arr.download(cloud.points, c);
    pcl::io::savePNGFile(filename, cloud);
}

template <typename T>
void savePNGFile(const std::string &filename, const pcl::PointCloud<T> &cloud) {
    pcl::io::savePNGFile(filename, cloud);
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
class PeoplePCDApp {
  public:
    typedef pcl::gpu::people::PeopleDetector PeopleDetector;

    enum { COLS = 640, ROWS = 480 };

    PeoplePCDApp(pcl::Grabber &capture)
        : capture_(capture), exit_(false), time_ms_(0), cloud_cb_(true),
          counter_(0), final_view_("Final labeling"), depth_view_("Depth") {
        final_view_.setSize(COLS, ROWS);
        depth_view_.setSize(COLS, ROWS);

        final_view_.setPosition(0, 0);
        depth_view_.setPosition(650, 0);

        cmap_device_.create(ROWS, COLS);
        cmap_host_.points.resize(COLS * ROWS);
        depth_device_.create(ROWS, COLS);
        image_device_.create(ROWS, COLS);

        depth_host_.points.resize(COLS * ROWS);

        rgba_host_.points.resize(COLS * ROWS);
        rgb_host_.resize(COLS * ROWS * 3);

        people::uploadColorMap(color_map_);
    }

    void visualizeAndWrite(bool write = false) {
        const PeopleDetector::Labels &labels =
            people_detector_.rdf_detector_->getLabels();
        people::colorizeLabels(color_map_, labels, cmap_device_);
        // people::colorizeMixedLabels(

        int c;
        cmap_host_.width = cmap_device_.cols();
        cmap_host_.height = cmap_device_.rows();
        cmap_host_.points.resize(cmap_host_.width * cmap_host_.height);
        cmap_device_.download(cmap_host_.points, c);

        final_view_.showRGBImage<pcl::RGB>(cmap_host_);
        final_view_.spinOnce(1, true);

        if (cloud_cb_) {
            depth_host_.width = people_detector_.depth_device1_.cols();
            depth_host_.height = people_detector_.depth_device1_.rows();
            depth_host_.points.resize(depth_host_.width * depth_host_.height);
            people_detector_.depth_device1_.download(depth_host_.points, c);
        }

        depth_view_.showShortImage(&depth_host_.points[0], depth_host_.width,
                                   depth_host_.height, 0, 5000, true);
        depth_view_.spinOnce(1, true);

        if (write) {
            if (cloud_cb_)
                savePNGFile(make_name(counter_, "ii"), cloud_host_);
            else
                savePNGFile(make_name(counter_, "ii"), rgba_host_);
            savePNGFile(make_name(counter_, "c2"), cmap_host_);
            savePNGFile(make_name(counter_, "s2"), labels);
            savePNGFile(make_name(counter_, "d1"),
                        people_detector_.depth_device1_);
            savePNGFile(make_name(counter_, "d2"),
                        people_detector_.depth_device2_);
        }
    }

    void
    source_cb1(const boost::shared_ptr<const PointCloud<PointXYZRGBA>> &cloud) {
        {
            boost::mutex::scoped_lock lock(data_ready_mutex_);
            if (exit_)
                return;

            pcl::copyPointCloud(*cloud, cloud_host_);
        }
        data_ready_cond_.notify_one();
    }

    void source_cb2(
        const boost::shared_ptr<openni_wrapper::Image> &image_wrapper,
        const boost::shared_ptr<openni_wrapper::DepthImage> &depth_wrapper,
        float) {
        {
            boost::mutex::scoped_try_lock lock(data_ready_mutex_);

            if (exit_ || !lock)
                return;

            // getting depth
            int w = depth_wrapper->getWidth();
            int h = depth_wrapper->getHeight();
            int s = w * PeopleDetector::Depth::elem_size;
            const unsigned short *data =
                depth_wrapper->getDepthMetaData().Data();
            depth_device_.upload(data, s, h, w);

            depth_host_.points.resize(w * h);
            depth_host_.width = w;
            depth_host_.height = h;
            std::copy(data, data + w * h, &depth_host_.points[0]);

            // getting image
            w = image_wrapper->getWidth();
            h = image_wrapper->getHeight();
            s = w * PeopleDetector::Image::elem_size;

            // fill rgb array
            rgb_host_.resize(w * h * 3);
            image_wrapper->fillRGB(w, h, (unsigned char *)&rgb_host_[0]);

            // convert to rgba, TODO image_wrapper should be updated to support
            // rgba directly
            rgba_host_.points.resize(w * h);
            rgba_host_.width = w;
            rgba_host_.height = h;
            for (int i = 0; i < rgba_host_.size(); ++i) {
                const unsigned char *pixel = &rgb_host_[i * 3];
                RGB &rgba = rgba_host_.points[i];
                rgba.r = pixel[0];
                rgba.g = pixel[1];
                rgba.b = pixel[2];
            }
            image_device_.upload(&rgba_host_.points[0], s, h, w);
        }
        data_ready_cond_.notify_one();
    }

    void startMainLoop() {
        cloud_cb_ = false;

        PCDGrabberBase *ispcd = dynamic_cast<pcl::PCDGrabberBase *>(&capture_);
        if (ispcd)
            cloud_cb_ = true;

        typedef boost::shared_ptr<openni_wrapper::DepthImage> DepthImagePtr;
        typedef boost::shared_ptr<openni_wrapper::Image> ImagePtr;

        boost::function<void(
            const boost::shared_ptr<const PointCloud<PointXYZRGBA>> &)>
            func1 = boost::bind(&PeoplePCDApp::source_cb1, this, _1);
        boost::function<void(const ImagePtr &, const DepthImagePtr &,
                             float constant)>
            func2 = boost::bind(&PeoplePCDApp::source_cb2, this, _1, _2, _3);
        boost::signals2::connection c = cloud_cb_
                                            ? capture_.registerCallback(func1)
                                            : capture_.registerCallback(func2);

        {
            boost::unique_lock<boost::mutex> lock(data_ready_mutex_);

            try {
                capture_.start();
                while (!exit_ && !final_view_.wasStopped()) {
                    bool has_data = data_ready_cond_.timed_wait(
                        lock, boost::posix_time::millisec(100));
                    if (has_data) {
                        SampledScopeTime fps(time_ms_);

                        if (cloud_cb_)
                            process_return_ = people_detector_.process(
                                cloud_host_.makeShared());
                        else
                            process_return_ = people_detector_.process(
                                depth_device_, image_device_);

                        ++counter_;
                    }

                    if (has_data && (process_return_ == 2))
                        visualizeAndWrite();
                }
                final_view_.spinOnce(3);
            } catch (const std::bad_alloc & /*e*/) {
                cout << "Bad alloc" << endl;
            } catch (const std::exception & /*e*/) {
                cout << "Exception" << endl;
            }

            capture_.stop();
        }
        c.disconnect();
    }

    boost::mutex data_ready_mutex_;
    boost::condition_variable data_ready_cond_;

    pcl::Grabber &capture_;

    bool cloud_cb_;
    bool exit_;
    int time_ms_;
    int counter_;
    int process_return_;
    PeopleDetector people_detector_;
    PeopleDetector::Image cmap_device_;
    pcl::PointCloud<pcl::RGB> cmap_host_;

    PeopleDetector::Depth depth_device_;
    PeopleDetector::Image image_device_;

    pcl::PointCloud<unsigned short> depth_host_;
    pcl::PointCloud<pcl::RGB> rgba_host_;
    std::vector<unsigned char> rgb_host_;

    PointCloud<PointXYZRGB> cloud_host_;

    ImageViewer final_view_;
    ImageViewer depth_view_;

    DeviceArray<pcl::RGB> color_map_;
};

void print_help() {
    cout << "\nPeople tracking app options (help):" << endl;
    cout << "\t -numTrees    \t<int> \tnumber of trees to load" << endl;
    cout << "\t -tree0       \t<path_to_tree_file>" << endl;
    cout << "\t -tree1       \t<path_to_tree_file>" << endl;
    cout << "\t -tree2       \t<path_to_tree_file>" << endl;
    cout << "\t -tree3       \t<path_to_tree_file>" << endl;
    cout << "\t -gpu         \t<GPU_device_id>" << endl;
    cout << "\t -dev         \t<Kinect_device_id>" << endl;
    cout << "\t -pcd         \t<path_to_pcd_file>" << endl;
    cout << "\t -oni         \t<path_to_oni_file>" << endl;
    cout << "\t -pcd_folder  \t<path_to_folder_with_pcf_files>" << endl;
}

int main(int argc, char **argv) {
    // answering for help
    PCL_INFO("People tracking App version 0.1\n");
    if (pc::find_switch(argc, argv, "--help") ||
        pc::find_switch(argc, argv, "-h"))
        return print_help(), 0;

    // selecting GPU and prining info
    int device = 0;
    pc::parse_argument(argc, argv, "-gpu", device);
    pcl::gpu::setDevice(device);
    pcl::gpu::printShortCudaDeviceInfo(device);

    // selecting data source
    boost::shared_ptr<pcl::Grabber> capture;
    string openni_device, oni_file, pcd_file, pcd_folder;

    try {
        if (pc::parse_argument(argc, argv, "-dev", openni_device) > 0) {
            capture.reset(new pcl::OpenNIGrabber(openni_device));
        } else if (pc::parse_argument(argc, argv, "-oni", oni_file) > 0) {
            capture.reset(new pcl::ONIGrabber(oni_file, true, true));
        } else if (pc::parse_argument(argc, argv, "-pcd", pcd_file) > 0) {
            capture.reset(new pcl::PCDGrabber<PointXYZRGBA>(
                vector<string>(31, pcd_file), 30, true));
        } else if (pc::parse_argument(argc, argv, "-pcd_folder", pcd_folder) >
                   0) {
            vector<string> pcd_files = getPcdFilesInDir(pcd_folder);
            capture.reset(
                new pcl::PCDGrabber<PointXYZRGBA>(pcd_files, 30, true));
        } else {
            capture.reset(new pcl::OpenNIGrabber());
            // capture.reset( new pcl::ONIGrabber("d:/onis/20111013-224932.oni",
            // true, true) );

            // vector<string> pcd_files(31, "d:/3/0008.pcd");
            // vector<string> pcd_files(31,
            // "d:/git/pcl/gpu/people/tools/test.pcd"); vector<string> pcd_files
            // = getPcdFilesInDir("d:/3/"); capture.reset( new
            // pcl::PCDGrabber<PointXYZRGBA>(pcd_files, 30, true) );
        }
    } catch (const pcl::PCLException & /*e*/) {
        return cout << "Can't open depth source" << endl, -1;
    }

    // selecting tree files
    vector<string> tree_files;
    tree_files.push_back("Data/forest1/tree_20.txt");
    tree_files.push_back("Data/forest2/tree_20.txt");
    tree_files.push_back("Data/forest3/tree_20.txt");
    tree_files.push_back("Data/forest4/tree_20.txt");

    pc::parse_argument(argc, argv, "-tree0", tree_files[0]);
    pc::parse_argument(argc, argv, "-tree1", tree_files[1]);
    pc::parse_argument(argc, argv, "-tree2", tree_files[2]);
    pc::parse_argument(argc, argv, "-tree3", tree_files[3]);

    int num_trees = (int)tree_files.size();
    pc::parse_argument(argc, argv, "-numTrees", num_trees);

    tree_files.resize(num_trees);
    if (num_trees == 0 || num_trees > 4)
        return cout << "Invalid number of trees" << endl, -1;

    try {
        // loading trees
        typedef pcl::gpu::people::RDFBodyPartsDetector RDFBodyPartsDetector;
        RDFBodyPartsDetector::Ptr rdf(new RDFBodyPartsDetector(tree_files));
        PCL_INFO("Loaded files into rdf");

        // Create the app
        PeoplePCDApp app(*capture);
        app.people_detector_.rdf_detector_ = rdf;

        // executing
        app.startMainLoop();
    } catch (const pcl::PCLException &e) {
        cout << "PCLException: " << e.detailedMessage() << endl;
    } catch (const std::runtime_error &e) {
        cout << e.what() << endl;
    } catch (const std::bad_alloc & /*e*/) {
        cout << "Bad alloc" << endl;
    } catch (const std::exception & /*e*/) {
        cout << "Exception" << endl;
    }

    return 0;
}

ext: cpp

max_stars_repo_name: yxlao/StanfordPCL

max_issues_repo_name: yxlao/StanfordPCL

max_forks_repo_name: yxlao/StanfordPCL


text:
#pragma once
#include <stdint.h>
#include <stdlib.h>
#include <boost/noncopyable.hpp>
#include <glog/logging.h>

namespace petuum {

/**
 * A buffer that allows appending records to it.
 */
class RecordBuff : boost::noncopyable {
public:
  RecordBuff() {
    LOG(FATAL) << "Default constructor called";
  }
  RecordBuff(void *mem, size_t size):
      mem_(reinterpret_cast<uint8_t*>(mem)),
      mem_size_(size),
      offset_(0) { }
  ~RecordBuff() { }

  RecordBuff(RecordBuff && other):
      mem_(other.mem_),
      mem_size_(other.mem_size_),
      offset_(other.offset_) {
    //VLOG(0) << "mem_size_ = " << mem_size_;
  };

  /**
   * does not take ownership of the memory
   */
  void *ResetMem(void *mem, size_t size) {
    void *old_mem = mem_;
    mem_ = reinterpret_cast<uint8_t*>(mem);
    mem_size_ = size;
    offset_ = 0;
    return old_mem;
  }

  void ResetOffset() {
    offset_ = 0;
  }

  bool Append(int32_t record_id, const void *record, size_t record_size) {
    if (offset_ + sizeof(int32_t) + record_size + sizeof(size_t) > mem_size_) {
      return false;
    }
    *(reinterpret_cast<int32_t*>(mem_ + offset_)) = record_id;
    offset_ += sizeof(int32_t);
    *(reinterpret_cast<size_t*>(mem_ + offset_)) = record_size;
    offset_ += sizeof(size_t);
    memcpy(mem_ + offset_, record, record_size);
    offset_ += record_size;
    //VLOG(0) << "Append() end offset = " << offset_;
    return true;
  }

  size_t GetMemUsedSize() {
    return offset_;
  }

  int32_t *GetMemPtrInt32() {
    if (offset_ + sizeof(int32_t) > mem_size_) {
      return 0;
    }
    int32_t *ret_ptr = reinterpret_cast<int32_t*>(mem_ + offset_);
    offset_ += sizeof(int32_t);

    return ret_ptr;
  }

  void PrintInfo() const {
    VLOG(0) << "mem_ = " << mem_
            << " mem_size_ = " << mem_size_
            << " offset = " << offset_;
  }

private:
  uint8_t *mem_;
  size_t mem_size_;
  size_t offset_;
};

}  // namespace petuum

ext: hpp

max_stars_repo_name: zuowang/bosen

max_issues_repo_name: zuowang/bosen

max_forks_repo_name: zuowang/bosen


text:
/***********************************************************************************************************************
 *  OpenStudio(R), Copyright (c) 2008-2017, Alliance for Sustainable Energy, LLC. All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the
 *  following conditions are met:
 *
 *  (1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following
 *  disclaimer.
 *
 *  (2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the
 *  following disclaimer in the documentation and/or other materials provided with the distribution.
 *
 *  (3) Neither the name of the copyright holder nor the names of any contributors may be used to endorse or promote
 *  products derived from this software without specific prior written permission from the respective party.
 *
 *  (4) Other than as required in clauses (1) and (2), distributions in any form of modifications or other derivative
 *  works may not use the "OpenStudio" trademark, "OS", "os", or any other confusingly similar designation without
 *  specific prior written permission from Alliance for Sustainable Energy, LLC.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 *  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 *  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER, THE UNITED STATES GOVERNMENT, OR ANY CONTRIBUTORS BE LIABLE FOR
 *  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 *  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 *  AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 **********************************************************************************************************************/

#include "IddFileFactoryData.hpp"
#include "WriteEnums.hpp"

#include "../utilities/idd/IddRegex.hpp"

#include <boost/regex.hpp>
#include <boost/algorithm/string.hpp>

#include <iostream>
#include <sstream>
#include <exception>

namespace openstudio {

IddFileFactoryData::IddFileFactoryData(const std::string& fileNameAndPathPair) {
  std::cout << "Creating new IddFileFactoryData object from input argument '"
            << fileNameAndPathPair << "'." << std::endl;
  std::stringstream ss;

  // decompose based on comma
  boost::regex re("([^,]+),?([^,]*)");
  boost::smatch match;
  bool ok = boost::regex_match(fileNameAndPathPair,match,re);
  if (!ok) {
    ss << "Unable to decompose " << fileNameAndPathPair
       << " into name and path. Type 'GenerateIddFactory --help' if usage is unclear.";
    throw std::runtime_error(ss.str().c_str());
  }
  std::string firstStr = std::string(match[1].first,match[1].second);
  std::string secondStr = std::string(match[2].first,match[2].second);
  boost::trim(firstStr);
  boost::trim(secondStr);
  if (secondStr.empty()) {
    secondStr = firstStr;
    firstStr.clear();
  }
  m_fileName = firstStr;
  m_filePath = path(secondStr);

  // validate filePath
  m_filePath = openstudio::filesystem::system_complete(m_filePath);
  if (!openstudio::filesystem::is_regular_file(m_filePath)) {
    ss << "Unable to locate intended idd file from input argument '"
       << fileNameAndPathPair << "'. User-supplied path resolved to '"
       << m_filePath.string() << "'.";
    throw std::runtime_error(ss.str().c_str());
  }

  // validate fileName
  if (m_fileName.empty()) {
    m_fileName = m_filePath.stem().string();
  }
  re = boost::regex("[a-zA-Z]\\w*");
  ok = boost::regex_match(m_fileName,match,re);
  if (!ok) {
    ss << "'" << m_fileName << "' is an invalid Idd file name. "
       << "Please use Perl regex word characters, and start the name with a letter.";
    throw std::runtime_error(ss.str().c_str());
  }

  std::cout << "IddFileFactoryData object created: " << std::endl
            << "  fileName = " << m_fileName << std::endl
            << "  filePath = " << m_filePath.string() << std::endl << std::endl;
}

void IddFileFactoryData::parseFile(const path& outPath,
                                   const std::string& outFileHeader,
                                   GenerateIddFactoryOutFiles& outFiles,
                                   int iddFileIndex)
{
  std::stringstream ss;

  openstudio::filesystem::ifstream iddFile(m_filePath);
  if (!iddFile) {
    ss << "Unable to open Idd file " << m_fileName << " located at " << m_filePath.string() << ".";
    throw std::runtime_error(ss.str().c_str());
  }

  std::stringstream header;
  std::string line,trimLine;
  boost::smatch matches;
  int lineNum = 1;
  std::shared_ptr<IddFactoryOutFile>& cxxFile = outFiles.iddFactoryIddFileCxxs[iddFileIndex];

  // get version
  std::getline(iddFile,line);
  trimLine = line; boost::trim(trimLine);
  bool ok = boost::regex_search(trimLine,matches,iddRegex::version());
  if (!ok) {
    ss << "Idd file " << m_fileName << " located at " << m_filePath.string()
       << " does not list its version on the first line, which is: " << std::endl
       << line;
    throw std::runtime_error(ss.str().c_str());
  }
  m_version = std::string(matches[1].first,matches[1].second);
  boost::trim(m_version);
  header << m_readyLineForOutput(line) << std::endl;

  // get rest of header
  while (std::getline(iddFile,line)) {
    ++lineNum; trimLine = line; boost::trim(trimLine);
    if (trimLine.empty()) { break; }
    if (!boost::regex_match(trimLine,iddRegex::commentOnlyLine())) {
      ss << "Idd file " << m_fileName << "'s header is not separated from other data by a blank line.";
      throw std::runtime_error(ss.str().c_str());
    }
    header << m_readyLineForOutput(line) << std::endl;
  }
  m_header = header.str();

  // parse body of file
  std::string group;
  std::string includeFile;
  while (std::getline(iddFile,line)) {
    ++lineNum; trimLine = line; boost::trim(trimLine);

    if (trimLine.empty()) { continue; }
    if (boost::regex_match(trimLine,iddRegex::commentOnlyLine())) { continue; }

    // look for group identifier
    if (boost::regex_search(trimLine,matches,iddRegex::group())) {
      group = std::string(matches[1].first,matches[1].second);
      boost::trim(group);
      continue;
    }

    // look for include file
    if (boost::regex_search(trimLine,matches,iddRegex::includeFile())) {
      includeFile = std::string(matches[1].first,matches[1].second);
      boost::trim(includeFile);
      m_includedFiles.push_back(FileNameRemovedObjectsPair(includeFile,std::vector<std::string>()));
      continue;
    }

    // look for object to remove from include file
    if (boost::regex_search(trimLine,matches,iddRegex::removeObject())) {
      if (includeFile.empty()) {
        ss << "Cannot process \\remove-object identifier: " << std::endl
           << line << std::endl
           << "found on line " << lineNum << " of Idd file '" << m_fileName << "', because it " << std::endl
           << "is not preceded by a non-empty \\include-file identifier.";
        throw std::runtime_error(ss.str().c_str());
      }
      std::string objectToRemove = std::string(matches[1].first,matches[1].second);
      objectToRemove = m_convertName(objectToRemove);
      m_includedFiles.back().second.push_back(objectToRemove);
      continue;
    }

    // look for IddObject
    StringPair objectName;
    ok = boost::regex_search(trimLine,matches,iddRegex::line());
    if (!ok) {
      ss << "Unable to extract object name from line " << lineNum << "," << std::endl
         << line << std::endl
         << ", of Idd file '" << m_fileName << "'.";
      throw std::runtime_error(ss.str().c_str());
    }
    objectName.second = std::string(matches[1].first,matches[1].second);
    boost::trim(objectName.second);
    objectName.first = m_convertName(objectName.second);
    m_objectNames.push_back(objectName);

    // start writing create function
    cxxFile->tempFile
      << std::endl
      << "IddObject create" << objectName.first << "IddObject() {" << std::endl
      << std::endl
      << "  static IddObject object;" << std::endl
      << std::endl
      << "  if (object.type() == IddObjectType::Catchall) {" << std::endl
      << "    std::stringstream ss;" << std::endl
      << "    ss << \"" << m_readyLineForOutput(line) << "\\n\";";

    // start collecting field names
    // (requires \field tag, which is expected to occur one per line)
    std::vector<std::string> fieldNames;
    std::vector<std::string> extensibleFieldNames;

    // get rest of object
    while (std::getline(iddFile,line)) {
      ++lineNum; trimLine = line; boost::trim(trimLine);
      if (trimLine.empty()) {
        // finish writing create function
        cxxFile->tempFile
          << std::endl
          << std::endl
          << "    IddObjectType objType(IddObjectType::" << objectName.first << ");" << std::endl
          << "    OptionalIddObject oObj = IddObject::load(\"" << objectName.second << "\"," << std::endl
          << "                                             \"" << group << "\"," << std::endl
          << "                                             ss.str()," << std::endl
          << "                                             objType);" << std::endl
          << "    OS_ASSERT(oObj);" << std::endl
          << "    object = *oObj;" << std::endl
          << "  }" << std::endl
          << std::endl
          << "  OS_ASSERT(object.type() == IddObjectType::" << objectName.first << ");" << std::endl
          << "  return object;" << std::endl
          << "}" << std::endl;

        // write field enums
        if (!fieldNames.empty() || !extensibleFieldNames.empty()) {
          // file for this object
          ss << objectName.first << "_FieldEnums.hxx";
          std::string filename = ss.str(); ss.str("");

          IddFactoryOutFile fieldEnumHxx(filename,outPath,outFileHeader);
          std::string upperObjectName(objectName.first);
          boost::algorithm::to_upper(upperObjectName);

          std::stringstream tempSS;

          fieldEnumHxx.tempFile
            << "#ifndef UTILITIES_IDD_" << upperObjectName << "_FIELDENUMS_HXX" << std::endl
            << "#define UTILITIES_IDD_" << upperObjectName << "_FIELDENUMS_HXX" << std::endl
            << std::endl
            << "#include <utilities/UtilitiesAPI.hpp>" << std::endl
            << "#include <utilities/core/Enum.hpp>" << std::endl
            << std::endl
            << "namespace openstudio {";

          if (!fieldNames.empty()) {
            std::list<std::pair<std::string, std::string>> fieldNamePairs;
            for (const std::string& name : fieldNames) {
              fieldNamePairs.emplace_back(m_convertName(name), name);
            }
            writeEnumFast(tempSS, objectName.first + "Fields", fieldNamePairs);

            fieldEnumHxx.tempFile
              << std::endl
              << "/** \\class " << objectName.first << "Fields" << std::endl
              << " *  \\brief Enumeration of " << objectName.second << "'s Idd and Idf non-extensible fields. " << std::endl
              << " *  \\details This enumeration provides readable, and semi-robust tags to use instead" << std::endl
              << " *  of bare unsigned integer field indices. (The enumeration value names are more likely" << std::endl
              << " *  to be valid across EnergyPlus and OpenStudio releases, and will break if they become" << std::endl
              << " *  invalid.) See the OPENSTUDIO_ENUM documentation in utilities/core/Enum.hpp. The actual macro" << std::endl
              << " *  call is:" << std::endl
              << " *" << std::endl
              << " *  \\code" << std::endl
//              << tempSS.str() << std::endl
              << " *  \\endcode */" << std::endl
              << tempSS.str() << std::endl;

            tempSS.str("");
          }

          if (!extensibleFieldNames.empty()) {
            std::list<std::pair<std::string, std::string>> extensibleFieldNamePairs;
            for (const std::string& name : extensibleFieldNames) {
              extensibleFieldNamePairs.emplace_back(m_convertName(name), name);
            }
            writeEnumFast(tempSS, objectName.first + "ExtensibleFields", extensibleFieldNamePairs);

            fieldEnumHxx.tempFile
              << std::endl
              << "/** \\class " << objectName.first << "ExtensibleFields" << std::endl
              << " *  \\brief Enumeration of " << objectName.second << "'s Idd and Idf extensible fields" << std::endl
              << " *  to be used in conjunction with ExtensibleIndex and ExtensibleGroup." << std::endl
              << " *  \\details This enumeration provides readable, and semi-robust tags to use instead" << std::endl
              << " *  of bare unsigned integer ExtensibleIndex.field indices. (The enumeration value names" << std::endl
              << " *  are more likely to be valid across EnergyPlus and OpenStudio releases, and will break" << std::endl
              << " *  if they become invalid.)See the OPENSTUDIO_ENUM documentation in utilities/core/Enum.hpp. The actual macro" << std::endl
              << " *  call is:" << std::endl
              << " *" << std::endl
              << " *  \\code" << std::endl
//              << tempSS.str() << std::endl
              << " *  \\endcode */" << std::endl
              << tempSS.str() << std::endl;

            tempSS.str("");
          }

          fieldEnumHxx.tempFile
            << std::endl
            << "} // openstudio" << std::endl
            << std::endl
            << "#endif // UTILITIES_IDD_" << upperObjectName << "_FIELDENUMS_HXX" << std::endl;
          outFiles.finalizeIddFactoryOutFile(fieldEnumHxx);

          // add to aggregate headers
          outFiles.iddFieldEnumsHxx.tempFile
            << "#include <utilities/idd/" << filename << ">" << std::endl;
          outFiles.iddFieldEnumsIxx.tempFile
            << "%include <utilities/idd/" << filename << ">" << std::endl;
        }

        break;
      }

      // continue writing create function
      cxxFile->tempFile
        << std::endl
        << "    ss << \"" << m_readyLineForOutput(line) << "\\n\";";

      // look for field name
      std::string fieldName;
      if (boost::regex_search(trimLine,matches,iddRegex::name())) {
        fieldName = std::string(matches[1].first,matches[1].second);
        fieldName = boost::regex_replace(fieldName,boost::regex(","),"");
        boost::trim(fieldName);
        if (extensibleFieldNames.empty()) {
          fieldNames.push_back(fieldName);
        }
        else {
          // strip out numbers
          fieldName = boost::regex_replace(fieldName,boost::regex("\\s?[0-9]+"),"");
          boost::trim(fieldName);
          extensibleFieldNames.push_back(fieldName);
        }
      }

      // look for extensible identifier
      if (boost::regex_search(trimLine,iddRegex::beginExtensible())) {
        if (fieldNames.empty()) {
          std::cerr << "Extensible processing unable to proceed because " << objectName.second
                    << " is extensible but is missing \\begin-extensible flag." << std::endl << std::endl;
        }
        assert(!fieldNames.empty());
        fieldName = fieldNames.back();
        // strip out numbers and transfer from fieldNames to extensibleFieldNames
        fieldName = boost::regex_replace(fieldName,boost::regex("\\s?[0-9]+"),"");
        boost::trim(fieldName);
        extensibleFieldNames.push_back(fieldName);
        fieldNames.pop_back();
      }

    } // while -- IddObject

  } // while -- IddFile

  iddFile.close();
  std::cout << "Parsed Idd file " << m_fileName << " located at " << m_filePath.string() << "," << std::endl
            << "which contains " << m_objectNames.size() << " objects." << std::endl << std::endl;
  if (!m_includedFiles.empty()) {
    for (const FileNameRemovedObjectsPair& p : m_includedFiles) {
      std::cout << "Idd file '" << m_fileName << "' includes all but " << p.second.size()
                << " objects of Idd file '" << p.first << "'." << std::endl << std::endl;
    }
  }

  // register objects with CallbackMap
  cxxFile->tempFile
    << "void IddFactorySingleton::register" << fileName() << "ObjectsInCallbackMap() {" << std::endl;
  for (const StringPair& objectName : objectNames()) {
    cxxFile->tempFile
      << "  m_callbackMap.insert(IddObjectCallbackMap::value_type(IddObjectType::"
      << objectName.first << ",create" << objectName.first << "IddObject));" << std::endl;
  }
  cxxFile->tempFile
    << "}" << std::endl
    << std::endl;
}

std::string IddFileFactoryData::fileName() const {
  return m_fileName;
}

std::string IddFileFactoryData::version() const {
  return m_version;
}

std::string IddFileFactoryData::header() const {
  return m_header;
}

std::vector<std::pair<std::string,std::string> > IddFileFactoryData::objectNames() const {
  return m_objectNames;
}

unsigned IddFileFactoryData::numIncludedFiles() const {
  return m_includedFiles.size();
}

IddFileFactoryData::FileNameRemovedObjectsPair IddFileFactoryData::includedFile(unsigned index) const {
  assert(index < numIncludedFiles());
  return m_includedFiles[index];
}

std::string IddFileFactoryData::m_convertName(const std::string& originalName) const {
  std::string result(originalName);
  boost::trim(result);
  result = boost::regex_replace(result,boost::regex("^100 ?%"),"All");
  result = boost::regex_replace(result,boost::regex("[ \\-]+"),"\\u");
  result = boost::regex_replace(result,boost::regex("[:/,\\(\\)%]+"),"_");
  result = boost::regex_replace(result,boost::regex("[\\.\\?\\]\\[]"),"");
  result = boost::regex_replace(result,boost::regex("="),"_EQUAL_");
  result = boost::regex_replace(result,boost::regex("\\*\\*"),"_POW_");
  result = boost::regex_replace(result,boost::regex("\\+"),"_PLUS_");
  result = boost::regex_replace(result,boost::regex("\\*"),"_TIMES_");
  return result;
}

std::string IddFileFactoryData::m_readyLineForOutput(const std::string& line) const {
  std::string result(line);
  result = boost::regex_replace(result,boost::regex("\\\\"),"\\\\\\\\");
  result = boost::regex_replace(result,boost::regex("\""),"\\\\\"");
  boost::trim(result);
  return result;
}

} // openstudio

max_stars_repo_name: OpenStudioThailand/OpenStudio

max_issues_repo_name: OpenStudioThailand/OpenStudio

max_forks_repo_name: OpenStudioThailand/OpenStudio


text:
/*=============================================================================
    Copyright (c) 2003 Joel de Guzman

    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)
==============================================================================*/
#if !defined(FUSION_ALGORITHM_PUSH_BACK_HPP)
#define FUSION_ALGORITHM_PUSH_BACK_HPP

#include <boost/spirit/fusion/sequence/single_view.hpp>
#include <boost/spirit/fusion/sequence/joint_view.hpp>

namespace boost { namespace fusion
{
    namespace meta
    {
        template <typename Sequence, typename T>
        struct push_back
        {
            typedef joint_view<Sequence, single_view<T> > type;
        };
    }

    namespace function
    {
        struct push_back
        {
            template <typename Sequence, typename T>
            struct apply : meta::push_back<Sequence, T> {};

            template <typename Sequence, typename T>
            inline typename apply<Sequence const, T>::type
            operator()(Sequence const& seq, T const& x) const
            {
                typedef joint_view<Sequence const, single_view<T> > result;
                single_view<T> val(x);
                return result(seq, val);
            }

            template <typename Sequence, typename T>
            inline typename apply<Sequence, T>::type
            operator()(Sequence& seq, T const& x) const
            {
                typedef joint_view<Sequence, single_view<T> > result;
                single_view<T> val(x);
                return result(seq, val);
            }
        };
    }

    function::push_back const push_back = function::push_back();
}}

#endif

max_stars_repo_name: Imperator-Knoedel/Sunset

max_issues_repo_name: Imperator-Knoedel/Sunset

max_forks_repo_name: Imperator-Knoedel/Sunset


text:
// -*- mode: c++; coding: utf-8; tab-width: 4 -*-

// Copyright (C) 2014, Oscar Acena <oscaracena@gmail.com>
// This software is under the terms of Apache License v2 or later.

#include <bluetooth/bluetooth.h>
#include <bluetooth/hci.h>
#include <bluetooth/hci_lib.h>

#include <exception>

#include "gattlib.h"
#include "gattservices.h"

#include <boost/algorithm/string.hpp>

DiscoveryService::DiscoveryService(const std::string device) :
	_device(device),
	_device_desc(-1) {

	int dev_id = hci_devid(device.c_str());
	if (dev_id < 0)
		throw std::runtime_error("Invalid device!");

	_device_desc = hci_open_dev(dev_id);
	if (_device_desc < 0)
		throw std::runtime_error("Could not open device!");
	}

DiscoveryService::~DiscoveryService() {
	if (_device_desc != -1)
		hci_close_dev(_device_desc);
}

void
DiscoveryService::enable_scan_mode() {
	int result;
	uint8_t scan_type = 0x01;
	uint16_t interval = htobs(0x0010);
	uint16_t window = htobs(0x0010);
	uint8_t own_type = 0x00;
	uint8_t filter_policy = 0x00;

	result = hci_le_set_scan_parameters
		(_device_desc, scan_type, interval, window,
		 own_type, filter_policy, 10000);

	if (result < 0)
		throw std::runtime_error
			("Set scan parameters failed (are you root?)");

	result = hci_le_set_scan_enable(_device_desc, 0x01, 1, 10000);
	if (result < 0)
		throw std::runtime_error("Enable scan failed");
}

void
DiscoveryService::get_advertisements(int timeout, boost::python::dict & ret) {
	struct hci_filter old_options;
	socklen_t slen = sizeof(old_options);
	if (getsockopt(_device_desc, SOL_HCI, HCI_FILTER,
				   &old_options, &slen) < 0)
		throw std::runtime_error("Could not get socket options");

	struct hci_filter new_options;
	hci_filter_clear(&new_options);
	hci_filter_set_ptype(HCI_EVENT_PKT, &new_options);
	hci_filter_set_event(EVT_LE_META_EVENT, &new_options);

	if (setsockopt(_device_desc, SOL_HCI, HCI_FILTER,
				   &new_options, sizeof(new_options)) < 0)
		throw std::runtime_error("Could not set socket options\n");

	int len;
	unsigned char buffer[HCI_MAX_EVENT_SIZE];
	struct timeval wait;
	fd_set read_set;
	wait.tv_sec = timeout;
	int ts = time(NULL);

	while(1) {
		FD_ZERO(&read_set);
		FD_SET(_device_desc, &read_set);

		int err = select(FD_SETSIZE, &read_set, NULL, NULL, &wait);
		if (err <= 0)
			break;

		len = read(_device_desc, buffer, sizeof(buffer));
        process_input(buffer, len, ret);

		int elapsed = time(NULL) - ts;
		if (elapsed >= timeout)
			break;

		wait.tv_sec = timeout - elapsed;
	}

	setsockopt(_device_desc, SOL_HCI, HCI_FILTER,
			   &old_options, sizeof(old_options));
}

void
DiscoveryService::process_input(unsigned char* buffer, int size, boost::python::dict & ret) {

    for (int i=0; i < size; i++) {
        unsigned char* ptr = buffer + HCI_EVENT_HDR_SIZE + 1;
        evt_le_meta_event* meta = (evt_le_meta_event*) ptr;

        if (meta->subevent != EVT_LE_ADVERTISING_REPORT) {
        //|| (uint8_t)buffer[BLE_EVENT_TYPE] != BLE_SCAN_RESPONSE) {
            //std::cout << "not an advertising report" << std::endl;
            continue;
        }

        le_advertising_info* info;
        info = (le_advertising_info*) (meta->data + 1);

        char addr[18];
        ba2str(&info->bdaddr, addr);

        int flags = parse_flags(info->data, info->length);
        int appearance = parse_appearance(info->data, info->length);
        std::string name = parse_name(info->data, info->length);
        boost::python::list uuids = parse_uuids(info->data, info->length);

        if (!ret.contains(addr)) {
        	boost::python::dict new_object;
        	ret[addr] = new_object;
    	}

    	if (!ret[addr].contains("flags") || flags != 0)
    		ret[addr]["flags"] = flags;

    	if (!ret[addr].contains("appearance") || appearance != 0)
    		ret[addr]["appearance"] = appearance;

    	// it's safe to cast the length to an unsigned int, since length > 0.
    	if (!ret[addr].contains("name") || name.length() > (unsigned int)boost::python::len(ret[addr]["name"]))
    		ret[addr]["name"] = name;

    	if (!ret[addr].contains("uuids") || boost::python::len(uuids) > boost::python::len(ret[addr]["uuids"]))
    		ret[addr]["uuids"] = uuids;

        ptr += sizeof(evt_le_meta_event);
    }
}

int
DiscoveryService::parse_appearance(uint8_t* data, size_t size) {
	size_t offset = 0;

	while (offset < size) {
		uint8_t field_len = data[0];

		if (field_len == 0 || offset + field_len > size) {
			// this probably ought to be an exception
			return 0;
		}

		if (data[1] == EIR_APPEARANCE) {
			uint16_t appearance = ((uint16_t *)(data + 2))[0];
			return appearance;
		}

		offset += field_len + 1;
		data += field_len + 1;
	}

	return 0;
}

boost::python::list
DiscoveryService::parse_uuids(uint8_t * data, size_t size) {
	size_t offset = 0;
	boost::python::list uuids;

	while (offset < size) {
		uint8_t field_len = data[0];

		if (field_len == 0 || offset + field_len > size) {
			// this probably ought to be an exception
			return uuids;
		}

		if (data[1] == EIR_16BIT_UUIDS_COMPLETE || data[1] == EIR_16BIT_UUIDS_INCOMPLETE) {
			int number_of_uuids = (field_len - 1) / 2;
			uint16_t *big_uuids = (uint16_t *)(data + 2);

			for (int i = 0; i < number_of_uuids; i++) {
				// add this uuid to the list
				uuids.append(big_uuids[i]);
			}
		} else if (data[1] == EIR_32BIT_UUIDS_COMPLETE || data[1] == EIR_32BIT_UUIDS_INCOMPLETE) {
			int number_of_uuids = (field_len - 1) / 4;
			uint32_t *big_uuids = (uint32_t *)(data + 2);

			for (int i = 0; i < number_of_uuids; i++) {
				// add this uuid to the list
				uuids.append(big_uuids[i]);
			}
		} else if (data[1] == EIR_128BIT_UUIDS_COMPLETE || data[1] == EIR_128BIT_UUIDS_INCOMPLETE) {
			int number_of_uuids = (field_len - 1) / 16;
			uint32_t *big_uuids = (uint32_t *)(data + 2);

			for (int i = 0; i < number_of_uuids * 4; i+=4) {
				// add this uuid to the list
				uuids.append(boost::python::make_tuple(big_uuids[i], big_uuids[i+1], big_uuids[i+2], big_uuids[i+3]));
			}
		}

		offset += field_len + 1;
		data += field_len + 1;
	}

	return uuids;
}

int
DiscoveryService::parse_flags(uint8_t* data, size_t size) {
	size_t offset = 0;

	while (offset < size) {
		uint8_t field_len = data[0];

		if (field_len == 0 || offset + field_len > size) {
			// this probably ought to be an exception
			return 0;
		}

		if (data[1] == EIR_FLAGS)
			return data[2];

		offset += field_len + 1;
		data += field_len + 1;
	}

	return 0;
}

std::string
DiscoveryService::parse_name(uint8_t* data, size_t size) {
	size_t offset = 0;
	std::string unknown = "";

	while (offset < size) {
		uint8_t field_len = data[0];
		size_t name_len;

		if (field_len == 0 || offset + field_len > size)
			return unknown;

		switch (data[1]) {
		case EIR_NAME_SHORT:
		case EIR_NAME_COMPLETE:
			name_len = field_len - 1;

            // remove any trailing \x00s
			if (data[1 + name_len] == '\0')
				name_len--;

			if (name_len > size)
				return unknown;

			return std::string((const char*)(data + 2), name_len);
		}

		offset += field_len + 1;
		data += field_len + 1;
	}

	return unknown;
}

void
DiscoveryService::disable_scan_mode() {
	if (_device_desc == -1)
		throw std::runtime_error("Could not disable scan, not enabled yet");

	int result = hci_le_set_scan_enable(_device_desc, 0x00, 1, 10000);
	if (result < 0)
		throw std::runtime_error("Disable scan failed");
}

boost::python::dict
DiscoveryService::discover(int timeout) {
	boost::python::dict retval;
	enable_scan_mode();
	get_advertisements(timeout, retval);
	disable_scan_mode();

	return retval;
}

max_stars_repo_name: mobilinkd/pygattlib

max_issues_repo_name: mobilinkd/pygattlib

max_forks_repo_name: mobilinkd/pygattlib


text:
#include <iostream>
#include <vector>
#include <unordered_map>

#include <boost/range/adaptors.hpp>
namespace ba = boost::adaptors;

#include <dionysus/simplex.h>
#include <dionysus/fields/zp.h>
#include <dionysus/fields/z2.h>
#include <dionysus/distances.h>
#include <dionysus/rips.h>
#include <dionysus/zigzag-persistence.h>
namespace d = dionysus;

#include <dionysus/dlog/progress.h>

#include <opts/opts.h>

#include <common.h>     // read_points()

typedef         std::vector<float>                                      Point;
typedef         std::vector<Point>                                      PointContainer;

typedef         d::PairwiseDistances<PointContainer,
                                     d::L2Distance<Point>>              PairDistances;
typedef         PairDistances::DistanceType                             DistanceType;
typedef         PairDistances::IndexType                                Vertex;

typedef         d::Rips<PairDistances>                                  Generator;
typedef         Generator::Simplex                                      Simplex;
typedef         std::set<Simplex>                                       SimplexSet;

typedef         std::vector<Vertex>                                     VertexVector;
typedef         std::vector<DistanceType>                               EpsilonVector;
typedef         std::tuple<Vertex,Vertex>                               Edge;
typedef         std::vector<Edge>                                       EdgeVector;

//typedef         d::Z2Field                                              K;
typedef         d::ZpField<>                                            K;
typedef         d::Simplex<>                                            Simplex;
typedef         d::ZigzagPersistence<K>                                 Persistence;
typedef         typename Persistence::Index                             Index;

typedef         std::unordered_map<Simplex, Index>                      Complex;
typedef         d::ChainEntry<K, Simplex>                               SimplexChainEntry;
typedef         d::ChainEntry<K, Index>                                 ChainEntry;

// debug
typedef         std::unordered_map<Index, Simplex>                      RComplex;

// Information we need to know when a class dies
struct      BirthInfo
{
    typedef         short unsigned                                      Dimension;

                    BirthInfo(DistanceType dist = DistanceType(), Dimension dim = Dimension()):
                        distance(dist), dimension(dim)              {}
    DistanceType    distance;
    Dimension       dimension;
};

typedef         std::unordered_map<Index, BirthInfo>                    BirthMap;

int main(int argc, char** argv)
{
    using opts::Options;
    using opts::Option;
    using opts::PosOption;

    short unsigned          skeleton = 2;
    DistanceType            multiplier = 6;
    short unsigned          p = 11;
    std::string             infilename, diagram_name;
    bool                    help;

    Options ops;
    ops
        >> Option('s', "skeleton",      skeleton,           "dimension of the Rips complex we want to compute")
        >> Option('m', "multiplier",    multiplier,         "multiplier for epsilon (distance to the next maxmin point)")
        >> Option('p', "prime",         p,                  "prime for arithmetic")
        >> Option('h', "help",          help,               "show help message")
    ;

    if (!ops.parse(argc,argv) || !(ops >> PosOption(infilename)) || !(ops >> PosOption(diagram_name)))
    {
        std::cout << "Usage: " << argv[0] << " input-points diagram.out" << std::endl;
        std::cout << ops;
        return 1;
    }

    PointContainer          points;
    read_points(infilename, points);

    std::ofstream   dgm_out(diagram_name);
    std::ostream&   out = dgm_out;

    // Construct distances and Rips generator
    PairDistances           distances(points);
    Generator               rips(distances);
    Generator::Evaluator    size(distances);

    // Order vertices and epsilons (in maxmin fashion)
    VertexVector        vertices;
    EpsilonVector       epsilons;
    EdgeVector          edges;
    DistanceType        inf     = std::numeric_limits<DistanceType>::infinity();

    {
        EpsilonVector   dist(distances.size(), inf);

        vertices.push_back(distances.begin());
        //epsilons.push_back(inf);
        while (vertices.size() < distances.size())
        {
            for (Vertex v = distances.begin(); v != distances.end(); ++v)
                dist[v] = std::min(dist[v], distances(v, vertices.back()));
            auto max = std::max_element(dist.begin(), dist.end());
            vertices.push_back(max - dist.begin());
            epsilons.push_back(*max);
        }
        epsilons.push_back(0);
    }

    // Generate and sort all the edges
    for (unsigned i = 0; i != vertices.size(); ++i)
        for (unsigned j = i+1; j != vertices.size(); ++j)
        {
            Vertex u = vertices[i];
            Vertex v = vertices[j];
            if (distances(u,v) <= multiplier*epsilons[j-1])
                edges.emplace_back(u,v);
        }
    std::sort(edges.begin(), edges.end(),
              [&distances](const Edge& e1, const Edge& e2)
              { return distances(std::get<0>(e1), std::get<1>(e1)) < distances(std::get<0>(e2), std::get<1>(e2)); });

    // Construct zigzag
    //K               k;
    K               k(p);
    Persistence     persistence(k);
    Complex         simplices;
#ifdef DIONYSUS_ZIGZAG_DEBUG
    RComplex        rsimplices;
#endif

    // Insert vertices
    Index       op   = 0;
    Index       cell = 0;
    BirthMap    births;
    for (auto v : vertices)
    {
        // Add a vertex
        Simplex s = {v};

        // We don't actually need to transform the boundary here,
        // since it's empty anyway, but we keep it for the sake of completeness
        Index pair = persistence.add(s.boundary(persistence.field()) |
                                                ba::transformed([&simplices](const SimplexChainEntry& e)
                                                { return ChainEntry(e.element(), simplices.find(e.index())->second); }));

#ifdef DIONYSUS_ZIGZAG_DEBUG
        rsimplices.emplace(cell, s);
        persistence.check_boundaries([&simplices](const Simplex& s) { return simplices[s]; },
                                     [&rsimplices](Index i)         { return rsimplices.find(i)->second; });
#endif

        births[op++] = BirthInfo(0,0);                  // record the birth
        simplices.emplace(std::move(s), cell++);        // record the cell id
    }

    // Process vertices
    dlog::progress progress(vertices.size());
    unsigned    ce = 0;         // index of the current one past last edge in the complex
    SimplexSet  cofaces;        // record the cofaces of all the simplices that need to be removed and reinserted
    for (unsigned stage = 0; stage != vertices.size() - 1; ++stage)
    {
        unsigned i = vertices.size() - 1 - stage;

        /* Increase epsilon */
        cofaces.clear();

        // Add anything else that needs to be inserted into the complex
        while (ce < edges.size())
        {
            Vertex u,v;
            std::tie(u,v) = edges[ce];
            if (distances(u,v) <= multiplier*epsilons[i-1])
                ++ce;
            else
                break;
            //std::cout << "Adding cofaces of " << u << ' ' << v << std::endl;
            rips.edge_cofaces(u, v,
                              skeleton,
                              multiplier*epsilons[i-1],
                              [&cofaces](Simplex&& s) { cofaces.insert(s); },
                              vertices.begin(),
                              vertices.begin() + i + 1);
        }

        // Insert all the cofaces
        for (auto& s : cofaces)
        {
            //std::cout << "Inserting: " << s << std::endl;

            Index pair = persistence.add(s.boundary(persistence.field()) |
                                                    ba::transformed([&simplices](const SimplexChainEntry& e)
                                                    { return ChainEntry(e.element(), simplices.find(e.index())->second); }));
            simplices.emplace(std::move(s), cell);      // record the cell id
#ifdef DIONYSUS_ZIGZAG_DEBUG
            rsimplices.emplace(cell, s);
            persistence.check_boundaries([&simplices](const Simplex& s) { return simplices[s]; },
                                         [&rsimplices](Index i)         { return rsimplices.find(i)->second; });
#endif
            ++cell;

            if (pair == Persistence::unpaired())
                births[op++] = BirthInfo(epsilons[i-1],s.dimension());              // record the birth
            else
            {
                const BirthInfo& birth = births[pair];
                if ((birth.distance - epsilons[i-1]) != 0 && birth.dimension < skeleton)
                    out << birth.dimension << " " << birth.distance << " " << epsilons[i-1] << std::endl;
                births.erase(pair);
                ++op;
            }
        }

        /* Remove the vertex */
        //std::cout << "Removing vertex: " << vertices[i] << std::endl;
        cofaces.clear();
        rips.vertex_cofaces(vertices[i],
                            skeleton,
                            multiplier*epsilons[i-1],
                            [&cofaces](Simplex&& s) { cofaces.insert(s); },
                            vertices.begin(),
                            vertices.begin() + i + 1);
        //std::cout << "Total cofaces: " << cofaces.size() << std::endl;

        for (auto& s : cofaces | ba::reversed)
        {
            //std::cout << "Removing: " << s << std::endl;
            Complex::const_iterator  it    = simplices.find(s);
            Index                    c     = it->second;
            simplices.erase(it);

            Index pair  = persistence.remove(c);
#ifdef DIONYSUS_ZIGZAG_DEBUG
            rsimplices.erase(c);
            persistence.check_boundaries([&simplices](const Simplex& s) { return simplices[s]; },
                                         [&rsimplices](Index i)         { return rsimplices.find(i)->second; });
#endif

            if (pair == Persistence::unpaired())
                births[op++] = BirthInfo(epsilons[i-1],s.dimension() - 1);          // record the birth
            else
            {
                const BirthInfo& birth = births[pair];
                if ((birth.distance - epsilons[i-1]) != 0 && birth.dimension < skeleton)
                    out << birth.dimension << " " << birth.distance << " " << epsilons[i-1] << std::endl;
                births.erase(pair);
                ++op;
            }
        }

        ++progress;
    }

    // Remove the last vertex
    Index pair = persistence.remove(0);
    simplices.erase((Complex::const_iterator) simplices.begin());     // TODO: add an assertion that the complex has only 1 simplex
#ifdef DIONYSUS_ZIGZAG_DEBUG
    rsimplices.erase(0);
    persistence.check_boundaries([&simplices](const Simplex& s) { return simplices[s]; },
                                 [&rsimplices](Index i)         { return rsimplices.find(i)->second; });
#endif

    const BirthInfo& birth = births[pair];
    out << birth.dimension << " " << birth.distance << " " << epsilons[0] << std::endl;
    ++progress;

    std::cout << "Finished" << std::endl;
}

max_stars_repo_name: dlm/dionysus

max_issues_repo_name: dlm/dionysus

max_forks_repo_name: dlm/dionysus


text:
//  (C) Copyright Gennadiy Rozental 2001.
//  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)

//  See http://www.boost.org/libs/test for the library home page.
//
//!@file
//!@brief Defines Unit Test Framework mono-state interfaces.
//! The framework interfaces are based on Monostate design pattern.
// ***************************************************************************

#ifndef BOOST_TEST_FRAMEWORK_HPP_020805GER
#define BOOST_TEST_FRAMEWORK_HPP_020805GER

// Boost.Test
#include <boost/test/detail/fwd_decl.hpp>
#include <boost/test/detail/global_typedef.hpp>
#include <boost/test/detail/suppress_warnings.hpp>
#include <boost/test/detail/throw_exception.hpp>

// STL
#include <stdexcept>

//____________________________________________________________________________//

namespace boost
{

/// Main namespace for the Unit Test Framework interfaces and implementation
namespace unit_test
{

// ************************************************************************** //
// **************              init_unit_test_func             ************** //
// ************************************************************************** //

/// Test module initialization routine signature

/// Different depending on whether BOOST_TEST_ALTERNATIVE_INIT_API is defined or
/// not
#ifdef BOOST_TEST_ALTERNATIVE_INIT_API
typedef bool (*init_unit_test_func)();
#else
typedef test_suite* (*init_unit_test_func)(int, char*[]);
#endif

// ************************************************************************** //
// **************                   framework                  ************** //
// ************************************************************************** //

/// Namespace of the Unit Test Framework mono-state
namespace framework
{

/// @name Unit Test Framework initialization and shutdown
/// @{

/// @brief This function performs initialization of the framework mono-state.
///
/// It needs to be called every time before the test is started.
/// @param[in] init_func test module initialization routine
/// @param[in] argc command line arguments collection
/// @param[in] argv command line arguments collection
BOOST_TEST_DECL void init(init_unit_test_func init_func, int argc,
                          char* argv[]);

/// This function applies all the decorators and figures out default run status.
/// This argument facilitates an ability of the test cases to prepare some other
/// test units (primarily used internally for self testing).
/// @param[in] tu Optional id of the test unit representing root of test tree.
/// If absent, master test suite is used
BOOST_TEST_DECL void finalize_setup_phase(test_unit_id tu = INV_TEST_UNIT_ID);

/// This function returns true when testing is in progress (setup is finished).
BOOST_TEST_DECL bool test_in_progress();

/// This function shuts down the framework and clears up its mono-state.
///
/// It needs to be at the very end of test module execution
BOOST_TEST_DECL void shutdown();
/// @}

/// @name Test unit registration
/// @{

/// Provides both read and write access to current "leaf" auto test suite during
/// the test unit registration phase.
///
/// During auto-registration phase the framework maintain a FIFO queue of test
/// units being registered. New test units become children of the current "leaf"
/// test suite and if this is test suite it is pushed back into queue and
/// becomes a new leaf. When test suite registration is completed, a test suite
/// is popped from the back of the queue. Only automatically registered test
/// suites should be added to this queue. Master test suite is always a zero
/// element in this queue, so if no other test suites are registered all test
/// cases are added to master test suite.

/// This function facilitates all three possible actions:
///    - if no argument are provided it returns the current queue leaf test
///    suite
///    - if test suite is provided and no second argument are set, test suite is
///    added to the queue
///    - if no test suite are provided and last argument is false, the semantic
///    of this function is similar to queue pop: last element is popped from the
///    queue
/// @param[in] ts test suite to push back to the queue
/// @param[in] push_or_pop should we push ts to the queue or pop leaf test suite
/// instead
/// @returns a reference to the currently active/"leaf" test suite
BOOST_TEST_DECL test_suite& current_auto_test_suite(test_suite* ts = 0,
                                                    bool push_or_pop = true);

/// This function add new test case into the global collection of test units the
/// framework aware of.

/// This function also assignes unique test unit id for every test case. Later
/// on one can use this id to locate the test case if necessary. This is the way
/// for the framework to maintain weak references between test units.
/// @param[in]  tc  test case to register
BOOST_TEST_DECL void register_test_unit(test_case* tc);

/// This function add new test suite into the global collection of test units
/// the framework aware of.

/// This function also assignes unique test unit id for every test suite. Later
/// on one can use this id to locate the test case if necessary. This is the way
/// for the framework to maintain weak references between test units.
/// @param[in]  ts  test suite to register
BOOST_TEST_DECL void register_test_unit(test_suite* ts);

/// This function removes the test unit from the collection of known test units
/// and destroys the test unit object.

/// This function also assigns unique test unit id for every test case. Later on
/// one can use this id to located the test case if necessary. This is the way
/// for the framework to maintain weak references between test units.
/// @param[in]  tu  test unit to deregister
BOOST_TEST_DECL void deregister_test_unit(test_unit* tu);

// This function clears up the framework mono-state.

/// After this call the framework can be reinitialized to perform a second test
/// run during the same program lifetime.
BOOST_TEST_DECL void clear();
/// @}

/// @name Test observer registration
/// @{
/// Adds new test execution observer object into the framework's list of test
/// observers.

/// Observer lifetime should exceed the the testing execution timeframe
/// @param[in]  to  test observer object to add
BOOST_TEST_DECL void register_observer(test_observer& to);

/// Excludes the observer object form the framework's list of test observers
/// @param[in]  to  test observer object to exclude
BOOST_TEST_DECL void deregister_observer(test_observer& to);

/// @}

/// @name Global fixtures registration
/// @{

/// Adds a new global fixture to be setup before any other tests starts and tore
/// down after any other tests finished. Test unit fixture lifetime should
/// exceed the testing execution timeframe
/// @param[in]  tuf  fixture to add
BOOST_TEST_DECL void register_global_fixture(test_unit_fixture& tuf);

/// Removes a test global fixture from the framework
///
/// Test unit fixture lifetime should exceed the testing execution timeframe
/// @param[in]  tuf  fixture to remove
BOOST_TEST_DECL void deregister_global_fixture(test_unit_fixture& tuf);
/// @}

/// @name Assertion/uncaught exception context support
/// @{
/// Context accessor
struct BOOST_TEST_DECL context_generator
{
    context_generator() : m_curr_frame(0)
    {
    }

    /// Is there any context?
    bool is_empty() const;

    /// Give me next frame; empty - last frame
    const_string next() const;

  private:
    // Data members
    mutable unsigned m_curr_frame;
};

/// Records context frame message.

/// Some context frames are sticky - they can only explicitly cleared by
/// specifying context id. Other (non sticky) context frames cleared after every
/// assertion.
/// @param[in] context_descr context frame message
/// @param[in] sticky is this sticky frame or not
/// @returns id of the newly created frame
BOOST_TEST_DECL int add_context(lazy_ostream const& context_descr, bool sticky);
/// Erases context frame (when test exits context scope)

/// If context_id is passed clears that specific context frame identified by
/// this id, otherwise clears all non sticky contexts.
BOOST_TEST_DECL void clear_context(int context_id = -1);
/// Produces an instance of small "delegate" object, which facilitates access to
/// collected context.
BOOST_TEST_DECL context_generator get_context();
/// @}

/// @name Access to registered test units.
/// @{
/// This function provides access to the master test suite.

/// There is only only master test suite per test module.
/// @returns a reference the master test suite instance
BOOST_TEST_DECL master_test_suite_t& master_test_suite();

/// This function provides an access to the test unit currently being executed.

/// The difference with current_test_case is about the time between a test-suite
/// is being set up or torn down (fixtures) and when the test-cases of that
/// suite start.

/// This function is only valid during test execution phase.
/// @see current_test_case_id, current_test_case
BOOST_TEST_DECL test_unit const& current_test_unit();

/// This function provides an access to the test case currently being executed.

/// This function is only valid during test execution phase.
/// @see current_test_case_id
BOOST_TEST_DECL test_case const& current_test_case();

/// This function provides an access to an id of the test case currently being
/// executed.

/// This function safer than current_test_case, cause if wont throw if no test
/// case is being executed.
/// @see current_test_case
BOOST_TEST_DECL test_unit_id current_test_case_id(); /* safe version of above */

/// This function provides access to a test unit by id and type combination. It
/// will throw if no test unit located.
/// @param[in]  tu_id    id of a test unit to locate
/// @param[in]  tu_type  type of a test unit to locate
/// @returns located test unit
BOOST_TEST_DECL test_unit& get(test_unit_id tu_id, test_unit_type tu_type);

/// This function template provides access to a typed test unit by id

/// It will throw if you specify incorrect test unit type
/// @tparam UnitType compile time type of test unit to get (test_suite or
/// test_case)
/// @param  id id of test unit to get
template <typename UnitType>
inline UnitType& get(test_unit_id id)
{
    return static_cast<UnitType&>(
        get(id, static_cast<test_unit_type>(UnitType::type)));
}
///@}

/// @name Test initiation interface
/// @{

/// Initiates test execution

/// This function is used to start the test execution from a specific "root"
/// test unit. If no root provided, test is started from master test suite. This
/// second argument facilitates an ability of the test cases to start some other
/// test units (primarily used internally for self testing).
/// @param[in] tu Optional id of the test unit or test unit itself from which
/// the test is started. If absent, master test suite is used
/// @param[in] continue_test true == continue test if it was already started,
/// false == restart the test from scratch regardless
BOOST_TEST_DECL void run(test_unit_id tu = INV_TEST_UNIT_ID,
                         bool continue_test = true);
/// Initiates test execution. Same as other overload
BOOST_TEST_DECL void run(test_unit const* tu, bool continue_test = true);
/// @}

/// @name Test events dispatchers
/// @{
/// Reports results of assertion to all test observers
BOOST_TEST_DECL void assertion_result(unit_test::assertion_result ar);
/// Reports uncaught exception to all test observers
BOOST_TEST_DECL void exception_caught(execution_exception const&);
/// Reports aborted test unit to all test observers
BOOST_TEST_DECL void test_unit_aborted(test_unit const&);
/// Reports aborted test module to all test observers
BOOST_TEST_DECL void test_aborted();
/// @}

namespace impl
{
// exclusively for self test
BOOST_TEST_DECL void setup_for_execution(test_unit const&);
BOOST_TEST_DECL void setup_loggers();
} // namespace impl

// ************************************************************************** //
// **************                framework errors              ************** //
// ************************************************************************** //

/// This exception type is used to report internal Boost.Test framework errors.
struct BOOST_TEST_DECL internal_error : public std::runtime_error
{
    internal_error(const_string m) :
        std::runtime_error(std::string(m.begin(), m.size()))
    {
    }
};

//____________________________________________________________________________//

/// This exception type is used to report test module setup errors.
struct BOOST_TEST_DECL setup_error : public std::runtime_error
{
    setup_error(const_string m) :
        std::runtime_error(std::string(m.begin(), m.size()))
    {
    }
};

#define BOOST_TEST_SETUP_ASSERT(cond, msg)                                     \
    BOOST_TEST_I_ASSRT(cond, unit_test::framework::setup_error(msg))

//____________________________________________________________________________//

struct nothing_to_test
{
    explicit nothing_to_test(int rc) : m_result_code(rc)
    {
    }

    int m_result_code;
};

//____________________________________________________________________________//

} // namespace framework
} // namespace unit_test
} // namespace boost

#include <boost/test/detail/enable_warnings.hpp>

#endif // BOOST_TEST_FRAMEWORK_HPP_020805GER

max_stars_repo_name: sotaoverride/backup

max_issues_repo_name: sotaoverride/backup

max_forks_repo_name: sotaoverride/backup


text:
#include <cctbx/boost_python/flex_fwd.h>

#include <scitbx/array_family/boost_python/flex_wrapper.h>
#include <scitbx/array_family/boost_python/flex_pickle_double_buffered.h>
#include <scitbx/array_family/boost_python/ref_pickle_double_buffered.h>
#include <boost/python/args.hpp>
#include <boost/python/return_internal_reference.hpp>

#include <spotfinder/core_toolbox/distl.h>

namespace scitbx { namespace af { namespace boost_python {

namespace {

  struct to_string_ir : pickle_double_buffered::to_string
  {
    using pickle_double_buffered::to_string::operator<<;

    to_string_ir()
    {
      int version = 1;
      *this << version;
    }

    to_string_ir& operator<<(Distl::icering const& val)
    {
      *this << val.lowerr2
            << val.upperr2
            << val.lowerresol
            << val.upperresol
            << val.strength;
      return *this;
    }
  };

  struct from_string_ir : pickle_double_buffered::from_string
  {
    from_string_ir(const char* str_ptr)
    : pickle_double_buffered::from_string(str_ptr)
    {
          *this >> version;
    }

    using pickle_double_buffered::from_string::operator>>;

    from_string_ir& operator>>(Distl::icering& val)
    {
      Distl::point lpeak;

      *this >> val.lowerr2
            >> val.upperr2
            >> val.lowerresol
            >> val.upperresol
            >> val.strength;
      return *this;
    }

    int version;
  };

  struct to_string : pickle_double_buffered::to_string
  {
    using pickle_double_buffered::to_string::operator<<;

    to_string()
    {
      int version = 2; //Current version for flex_distl_spot
      //printf("Pickle flex_distl_spot; version: %d\n", version);
      *this << version;
    }

    to_string& operator<<(scitbx::vec2<double> const& vec)
    {
      //special time saver for 2D profiles:  ignore z-direction
      *this << (float)vec[0] << (float)vec[1];
      return *this;
    }

    to_string& operator<<(spotfinder::distltbx::w_spot const& val)
    {
      *this << val.max_pxl_x()
            << val.max_pxl_y()
            << val.peakintensity
            << val.bodypixels.size();

      for (int i=0; i<val.bodypixels.size(); ++i){
      *this << val.bodypixels[i].x
            << val.bodypixels[i].y;
      }

      *this << val.total_mass;
      to_string::operator<<(val.model_m->center_of_mass());
      *this << (float)val.model_m->eigenvalue(0)
            << (float)val.model_m->eigenvalue(1);

      to_string::operator<<(val.model_m->eigenvector(0));
      to_string::operator<<(val.model_m->eigenvector(1));

      return *this;
    }
  };

  struct from_string : pickle_double_buffered::from_string
  {
    from_string(const char* str_ptr)
    : pickle_double_buffered::from_string(str_ptr)
    {
      *this >> version;
      //printf("Unpickle flex_distl_spot; version: %d\n", version);
    }

    from_string& operator>>(scitbx::vec2<double>& vec)
    {
      //special time saver for 2D profiles:  ignore z-direction
      float a,b;
      *this >> a >> b;
      vec[0] = (double)a; vec[1] = (double)b;
      return *this;
    }

    using pickle_double_buffered::from_string::operator>>;

    from_string& operator>>(spotfinder::distltbx::w_spot& val)
    {
      int iptx,ipty,bodylen;
      Distl::point lpeak;

      *this >> iptx
            >> ipty;
      lpeak = Distl::point(iptx,ipty);
      *this >> val.peakintensity
            >> bodylen;

      val.bodypixels = scitbx::af::shared<Distl::point>();
      for (int i=0; i<bodylen; ++i){
      *this >> iptx
            >> ipty;
        val.bodypixels.push_back(Distl::point(iptx,ipty));
      }

      val.setstate(lpeak);
      if (version == 2) {
        double mass;
        *this >> mass;
        scitbx::vec2<double> com,value,axis0,axis1;
        float a,b;
        from_string::operator>>(com);
        *this >> a >> b;
        value[0] = (double)a; value[1] = (double)b;
        from_string::operator>>(axis0);
        from_string::operator>>(axis1);
        val.ss_setstate(com,value,axis0,axis1,mass);
      }
      return *this;
    }

    int version;
  };

}}}} // namespace scitbx::af::boost_python::<anonymous>

namespace scitbx { namespace af { namespace boost_python {

  void wrap_flex_icering()
  {
    flex_wrapper<Distl::icering,
                 boost::python::return_internal_reference<>
                 >::plain("distl_icering")
    .def_pickle(flex_pickle_double_buffered<
                 Distl::icering, to_string_ir, from_string_ir>())
    ;
  }

  af::shared<double>
  ctr_mass_distances_from_direct_beam(
    af::const_ref<spotfinder::distltbx::w_spot> const& spots,
    scitbx::vec2<double> detector_size,
    scitbx::vec2<int> detector_pixels,
    scitbx::vec2<double> xy_beam)
  {
    af::shared<double> result(
      spots.size(),
      af::init_functor_null<double>());
    scitbx::vec2<double> sop;
    for(std::size_t i=0;i<2;i++) {
      sop[i] = detector_size[i] / detector_pixels[i];
    }
    double const* b = xy_beam.begin();
    for(std::size_t i=0;i<spots.size();i++) {
      spotfinder::distltbx::w_spot const& spot = spots[i];
      double dx = spot.ctr_mass_x() * sop[0] - b[0];
      double dy = spot.ctr_mass_y() * sop[1] - b[1];
      result[i] = std::sqrt(dx*dx + dy*dy);
    }
    return result;
  }

  void wrap_flex_w_spot()
  {
    using boost::python::arg;
    flex_wrapper<spotfinder::distltbx::w_spot,
                 boost::python::return_internal_reference<>
                 >::plain("distl_spot")
    .def_pickle(flex_pickle_double_buffered<
                 spotfinder::distltbx::w_spot, to_string, from_string>())
    .def("ctr_mass_distances_from_direct_beam",
      ctr_mass_distances_from_direct_beam, (
        arg("detector_size"),
        arg("detector_pixels"),
        arg("xy_beam")))
    ;
  }

  void wrap_flex_point()
  {
    flex_wrapper<Distl::point>::plain("distl_point");
  }

}}} // namespace scitbx::af::boost_python