text:
#include <pybind11/pybind11.h>
#include <boost/function_output_iterator.hpp>
#include <boost/iterator/transform_iterator.hpp>
#include <boost/range/iterator_range_core.hpp>

#include "modules/io/utils.h"
#include "python/biograph/internal/assemble.h"
#include "python/biograph/internal/anchor.h"

variant_wrapper::variant_wrapper(const sv_out& var,
                                 std::shared_ptr<assembly> src,
                                 std::shared_ptr<reference> ref)
    : m_var(var),
      m_assembly(src),
      m_ref(ref),
      m_min_depth(255),
      m_max_depth(0),
      m_avg_depth(0.0) {
  for (int i = m_var.seq_begin - 1; i < m_var.seq_end + 1; i++) {
    uint8_t depth = m_assembly->depth[i];
    m_min_depth = std::min(depth, m_min_depth);
    m_max_depth = std::max(depth, m_max_depth);
    m_avg_depth += depth;
  }
  m_avg_depth /= 2 + m_var.seq_end - m_var.seq_begin;

  // canonicalize positions
  if (!left_forward() and !right_forward()) {
    m_var = m_var.flip(m_assembly->assembly.size());
  }
  // It can be reversed, just want left less than right
  if (left_position() > right_position()) {
    m_var = m_var.flip(m_assembly->assembly.size());
  }

}

std::string variant_wrapper::str() const {
  if (m_var.is_structural) {
    return printstring("%s:%lu%c_%d_%s:%lu%c", left_contig().c_str(),
                       left_position(), left_forward() ? '+' : '-',
                       m_var.seq_end - m_var.seq_begin, right_contig().c_str(),
                       right_position(), right_forward() ? '+' : '-');
  } else {
    return printstring(
        "%s:%lu_%lu:%d", left_contig().c_str(), left_position(),
        m_var.right_ref - m_var.left_ref, m_var.seq_end - m_var.seq_begin);
  }
}

std::string variant_wrapper::repr() const {
  return std::string("<biograph.Variant ") + str() + std::string(">");
}

std::string variant_wrapper::left_contig() const {
  size_t flat = m_var.left_ref.get_offset();
  return m_ref->get_assembly()
      .scaffold_order[m_ref->get_seq_position(flat).scaffold_id];
}

std::string variant_wrapper::right_contig() const {
  size_t flat = m_var.right_ref.get_offset();
  return m_ref->get_assembly()
      .scaffold_order[m_ref->get_seq_position(flat).scaffold_id];
}

unsigned long variant_wrapper::left_position() const {
  size_t flat = m_var.left_ref.get_offset();
  return m_ref->get_seq_position(flat).position;
}

unsigned long variant_wrapper::right_position() const {
  size_t flat = m_var.right_ref.get_offset();
  return m_ref->get_seq_position(flat).position;
}

dna_sequence variant_wrapper::variant_sequence() const {
  return m_assembly->assembly.subseq(m_var.seq_begin,
                                     m_var.seq_end - m_var.seq_begin);
}

dna_sequence variant_wrapper::assembly_sequence() const {
  return m_assembly->assembly;
}

reference_range variant_wrapper::range() const {
  size_t lflat = m_var.left_ref.get_offset();
  size_t rflat = m_var.right_ref.get_offset();
  return reference_range(m_ref, std::min(lflat, rflat),
                         std::max(lflat, rflat));
}

py_list variant_wrapper::depths() const {
  py_list ret;
  for (int i = m_var.seq_begin - 1; i < m_var.seq_end + 1; i++) {
    ret.append(m_assembly->depth[i]);
  }
  return ret;
}

py_list variant_wrapper::assembly_depths() const {
  py_list ret;
  for (uint8_t depth : m_assembly->depth) {
    ret.append(depth);
  }
  return ret;
}

variant_wrapper variant_wrapper::flip() const {
  variant_wrapper out;
  out.m_var = m_var.flip(m_assembly->assembly.size());
  out.m_assembly = std::make_shared<assembly>(m_assembly->flip());
  out.m_ref = m_ref;
  out.m_min_depth = m_min_depth;
  out.m_max_depth = m_max_depth;
  out.m_avg_depth = m_avg_depth;
  return out;
}

static void add_reference_depths(py_list out, dna_const_iterator a,
                                 dna_const_iterator b, const uint8_t* depths,
                                 const reference& ref) {
  bool reverse = false;
  if (a.is_rev_comp()) {
    dna_const_iterator t = a;
    a = b.rev_comp();
    b = t.rev_comp();
    reverse = true;
  }

  if (a.is_rev_comp() || b.is_rev_comp()) {
    throw io_exception("Impossible case #1 in add_reference_depths");
  }

  seq_position sp_a = ref.get_seq_position(a.get_offset());
  seq_position sp_b = ref.get_seq_position(b.get_offset());
  if (sp_a.scaffold_id != sp_b.scaffold_id) {
    throw io_exception("Impossible case #2 in add_reference_depths");
  }

  py_list py_dp;
  for (size_t i = 0; i < sp_b.position - sp_a.position + 1; i++) {
    if (reverse) {
      py_dp.append(depths[sp_b.position - sp_a.position - i]);
    } else {
      py_dp.append(depths[i]);
    }
  }

  std::string scaf = ref.get_assembly().scaffold_order[sp_a.scaffold_id];
  // printf("REF: %s: [%lu, %lu]\n", scaf.c_str(), sp_a.position,
  // sp_b.position);

  py_list all;
  all.append(scaf);
  all.append(sp_a.position + 1);
  all.append(py_dp);
  out.append(all);
}

py_list py_assemble(py_list src, py_list dest, uint8_t min_overlap,
                    uint32_t max_steps, bool skip_ambig,
                    const std::shared_ptr<readmap>& wrap) {
  std::shared_ptr<seqset_bitmap_base> bm = wrap;

  py_list ret;
  py_list ref_seq;

  if (pybind11::len(src) == 0 || pybind11::len(dest) == 0) {
    return ret;
  }

  anchor_wrapper x = pybind11::cast<anchor_wrapper>(src[0]);
  anchor_wrapper y = pybind11::cast<anchor_wrapper>(dest[0]);

  if (x.m_reference != y.m_reference || x.m_seqset != y.m_seqset) {
    throw io_exception("Using anchors with mismatched References or Seqsets");
  }

  const reference& ref = *x.m_reference;

  auto out_it =
      boost::make_function_output_iterator([&](const assembly& assem) {

        std::shared_ptr<assembly> pa = std::make_shared<assembly>(assem);
        dna_slice left_range = ref.get_supercontig(assem.left.get_offset());
        dna_slice right_range = ref.get_supercontig(assem.right.get_offset());

        if (assem.left.is_rev_comp()) {
          left_range = left_range.rev_comp();
        }

        if (assem.right.is_rev_comp()) {
          right_range = right_range.rev_comp();
        }

        left_range = dna_slice(assem.left, left_range.end());
        right_range = dna_slice(right_range.begin(), assem.right + 1);

        std::vector<sv_out> bits =
            call_structural(assem.assembly, left_range, right_range, 200);
        if (bits.size() == 0) {
          return;
        }

        std::shared_ptr<assembly> par;
        dna_const_iterator it_ref = assem.left;
        const uint8_t* depth_start = &assem.depth[0];

        for (size_t i = 0; i < bits.size(); i++) {
          if (bits[i].anchor_drop) {
            // printf("Drop in assembly: %d\n", assem.id);
            if (bits[i].seq_begin == 0) {
              it_ref = bits[i].right_ref;
              pa->left = it_ref;
            } else {
              pa->right = bits[i].left_ref;
            }
            continue;
          }

          add_reference_depths(ref_seq, it_ref, bits[i].left_ref, depth_start,
                               ref);
          it_ref = bits[i].right_ref;
          depth_start = &assem.depth[bits[i].seq_end];

          // Get Ref between it_ref + left_ref
          // Update it_ref to right_ref
          if (!bits[i].is_structural && bits[i].left_ref.is_rev_comp()) {
            // TODO: This code is UNTESTED!!!
            sv_out rb = bits[i];
            rb.left_ref = bits[i].right_ref.rev_comp();
            rb.right_ref = bits[i].left_ref.rev_comp();
            rb.seq_begin = assem.assembly.size() - bits[i].seq_end;
            rb.seq_end = assem.assembly.size() - bits[i].seq_begin;
            if (!par) {
              par = std::make_shared<assembly>(pa->flip());
            }
            ret.append(variant_wrapper(rb, par, x.m_reference));
          } else {
            ret.append(variant_wrapper(bits[i], pa, x.m_reference));
          }
        }

        add_reference_depths(ref_seq, it_ref, pa->right, depth_start, ref);
      });

  auto transform_func = [&](pybind11::handle h) {
    return pybind11::cast<const anchor_wrapper&>(h).m_anchor;
  };

  auto left = boost::make_iterator_range(
      boost::make_transform_iterator(src.begin(), transform_func),
      boost::make_transform_iterator(src.end(), transform_func));

  auto right = boost::make_iterator_range(
      boost::make_transform_iterator(dest.begin(), transform_func),
      boost::make_transform_iterator(dest.end(), transform_func));

  seqset_assemble(out_it, x.m_seqset->get_seqset(), left, right, min_overlap,
                  max_steps, skip_ambig, *wrap);

  py_list both;
  both.append(ret);
  both.append(ref_seq);
  return both;
}

bool variant_wrapper::py_lessthan(variant_wrapper other) {
  unsigned long mStart =
      std::min(this->left_position(), this->right_position());
  unsigned long oStart =
      std::min(other.left_position(), other.right_position());
  return mStart < oStart;
}

int variant_wrapper::py_cmp(variant_wrapper other) {
  if (this->left_contig() == other.left_contig() &&
      this->right_contig() == other.right_contig() &&
      this->left_position() == other.left_position() &&
      this->right_position() == other.right_position() &&
      this->variant_sequence() == other.variant_sequence()) {
    return 0;
  }
  if (this->py_lessthan(other)) {
    return -1;
  }

  return 1;
}

using namespace pybind11;
void bind_assemble(module& m) {
  class_<variant_wrapper>(m,
      "Variant",
      "An object representing a variation from reference.\n"
      "\n"
      "Left and Right correspond to where assemblies that contain variation from the reference align\n"
      "For example, a deletion's assembly would look like the following:\n"
      "  ATCAG----CAGAT\n"
      "  ATCAGAGATCAGAT\n"
      "Where:\n"
      "  variant.left_position = 4\n"
      "  variant.right_position = 9\n"
      "  variant.sequence = ''\n"
      "  variant.range.sequence = 'GAGAT'\n"
      "For indels, left_forward and right_forward will not match (one True,\n"
      "one False).\n"
      "\n"
      "Variants with both left_forward and right_forward True or False are\n"
      "inversions.\n"
      "")
      .def("__str__", &variant_wrapper::str,
           "A unique identifier for this variant.")
      .def("__repr__", &variant_wrapper::repr,
           "The python representation of this variant object.\n")
      .def("__lt__", &variant_wrapper::py_lessthan,
           "Is this variant's chrom:start less than another's\n")
      .def("__cmp__", &variant_wrapper::py_cmp,
           "Compare these variants' position and sequence")
      .def_property_readonly("is_structural", &variant_wrapper::is_structural,
                    "True if this is a structural variant.")
      .def_property_readonly(
           "left_contig", &variant_wrapper::left_contig,
           "The left side contig for this variant. Structural variants may\n"
           "cross supercontig boundaries.")
      .def_property_readonly(
           "right_contig", &variant_wrapper::right_contig,
           "The right side contig for this variant. Structural variants may\n"
           "cross supercontig boundaries.")
      .def_property_readonly("left_position", &variant_wrapper::left_position,
                    "The left position for this variant (pointing to the "
                    "reference base\n"
                    "immediately before the variant).")
      .def_property_readonly("right_position", &variant_wrapper::right_position,
                    "The right position for this variant (pointing to the "
                    "reference base\n"
                    "immediately after the variant).")
      .def_property_readonly("left_forward", &variant_wrapper::left_forward,
                    "True if the left side is in the forward direction. Always "
                    "true for\n"
                    "non-structural variants.")
      .def_property_readonly("right_forward", &variant_wrapper::right_forward,
                    "True if the right side is in the forward direction. "
                    "Always true for\n"
                    "non-structural variants.")
      .def_property_readonly("sequence", &variant_wrapper::variant_sequence,
                    "A Sequence object representing the variant nucleotide\n"
                    "sequence (if any).")
      .def_property_readonly("range", &variant_wrapper::range,
                    "A ReferenceRange over the [self.left_position:self.right_position) in\n"
                    "the reference where this variant occurs.")
      .def_property_readonly(
           "depths", &variant_wrapper::depths,
           "A list of read depth counts for each base in this variant,\n"
           "including the bases immediately preceding and following it in\n"
           "the reference.\n"
           "\n"
           "For example, a SNP would have depth list containing three values:\n"
           "the count of reads that match the last base in the reference "
           "prior\n"
           "to the SNP, the SNP itself, and the next base in the reference.")
      .def_property_readonly(
           "min_depth", &variant_wrapper::min_depth,
           "The minimum depth of coverage of any base in this variant.")
      .def_property_readonly(
           "max_depth", &variant_wrapper::max_depth,
           "The maximum depth of coverage of any base in this variant.")
      .def_property_readonly(
           "avg_depth", &variant_wrapper::avg_depth,
           "The average depth of coverage of all bases in this variant (float)")
      .def_property_readonly("min_overlap", &variant_wrapper::min_overlap,
                    "The lowest overlap found in the entire assembly.")
      .def_property_readonly(
           "assembly_id", &variant_wrapper::assembly_id,
           "A unique identifier for the assembly where this variant was "
           "found.\n"
           "Note that several variants may be found within a single assembly.")
      .def_property_readonly("assembly_sequence", &variant_wrapper::assembly_sequence,
                    "A Sequence object representing the entire assembly\n"
                    "including all variants.")
      .def_property_readonly("assembly_depths", &variant_wrapper::assembly_depths,
                    "A list of coverage depths for every base in the assembled "
                    "sequence.")
      .def_property_readonly("assembly_begin", &variant_wrapper::assembly_begin,
                    "The offset within this assembly where the variant begins "
                    "(inclusive).")
      .def_property_readonly("assembly_end", &variant_wrapper::assembly_end,
                    "The offset within this assembly where the variant ends "
                    "(exclusive).")
      .def("flip", &variant_wrapper::flip,
           "A Sequence object for the reverse complement of this variant.");

  m.def("assemble", py_assemble,
      "assemble([src], [dest], min_overlap, max_steps, skip_ambiguous, "
      "readmap)\n"
      "\n"
      "Find all possible assemblies between the lists of anchors in src\n"
      "and dest, for all reads matching the given Readmap.\n"
      "\n"
      "Stop when all assemblies are found or max_steps have\n"
      "been taken. If skip_ambiguous is True, don't assemble across\n"
      "ambiguous regions of the reference.\n"
      "\n"
      "Returns a list containing two sublists: \n"
      "\n"
      "  * assemblies[0] contains the variants called from all assemblies\n"
      "    found, with one entry for each variant. Since there can be\n"
      "    many variants for each assembly, this list can be much larger\n"
      "    than the total number of assemblies.\n"
      "\n"
      "  * assemblies[1] contains a list of each assembly's coverage\n"
      "    depth. There is one entry for each assembly found.\n"
      "\n"
      "NOTE: Some of the coverage of each assembly contributes to the\n"
      "reference coverage (ie. the parts of the anchor that match\n"
      "reference). This coverage is not included in reference_range\n"
      "coverage, since that can only include reads that precisely\n"
      "match reference.\n"
      "\n"
      "To calculate the true reference coverage for the anchors,\n"
      "you will need to add the assembly coverage for bases that\n"
      "lie in the reference back to the ReferenceRange coverage.\n"
      "\n"
      "Each coverage list is of the format:\n"
      "\n"
      "    ['supercontig', position, [25,25,26,26,26...]]\n"
      "\n"
      "Example:\n"
      "\n"
      "  from biograph import BioGraph, Reference, find_anchors, assemble\n"
      "  my_bg = BioGraph('NA12878.bg')\n"
      "  ref = Reference('/reference/human_g1k_v37/')\n"
      "  rr = ref.make_range('2', 1, 200000, False)\n"
      "  fwd = find_anchors(my_bg, rr, True, 70, 10000)\n"
      "  rev = find_anchors(my_bg, rr, False, 70, 10000)\n"
      "  both = assemble(fwd, rev, 70, 100000, True)\n"
      "  print both[0][0].sequence # the first variant\n");
}

max_stars_repo_name: spiralgenetics/biograph

max_issues_repo_name: spiralgenetics/biograph

max_forks_repo_name: spiralgenetics/biograph


text:
/* */
/**
 * Simple stand-alone ROS node that takes data from NetFT sensor and
 * Publishes it ROS topic
 */

#include "ros/ros.h"
#include "marvin_common/WrenchData.h"
#include "marvin_common/TimeRequest.h"
#include "marvin_common/MarvinUtils.hpp"
#include "marvin_common_rw/RwRos.hpp"
#include <rw/common/PropertyMap.hpp>
#include <rw/common/TimerUtil.hpp>
#include <rw/loaders/xml/XMLPropertyLoader.hpp>

#include <rwhw/netft/NetFTLogging.hpp>
#include <rwhw/netft/FTCompensation.hpp>
#include <boost/foreach.hpp>
#include <iostream>
#include <rw/models/WorkCell.hpp>

using namespace rw::common;
using namespace rw::loaders;
using namespace rwhw;
using namespace rw::models;
using namespace rw::math;

rwhw::FTCompensation *_ftCompensation;
std::pair<rw::math::Vector3D<>, rw::math::Vector3D<> > _wrench;

Wrench3D _ft_current;
ros::Time _lastTime;
void handleRobotState(const marvin_common::RobotState& state){
	//
	Q q = RwRos::toRw(state.q);
	Q dq = RwRos::toRw(state.dq);
	ros::Time _currentTime = ros::Time::now();
    _ftCompensation->update(_wrench, q, dq, (_currentTime-_lastTime).toSec() );
    _lastTime = _currentTime;
}

// initialize the ros stuff for publishing raw wrench data
void initWrenchRaw();
// initialize the ros stuff for publishing unbiased wrench data
void initWrenchUnbiased();

int main(int argc, char **argv)
{
	// TODO: add boost options here if we want to parse command line options
	// for now we use parameter server

	ros::init(argc, argv, "netft");
	ros::NodeHandle nh("~");
	std::string name = nh.getNamespace();
	std::string propertyFile;
<<<<<<< .mine
	std::string host, robotName, robotState;
=======
	std::string host, devName, calibfile;

>>>>>>> .r1128
	int publishRate = 400;
	nh.getParam("properties", propertyFile);
	nh.param("host", host, std::string(""));
	nh.param("rate", publishRate, 400);
<<<<<<< .mine
	nh.param("robot", robotName, "");
	nh.param("robotState", robotState, "");

	bool unbiasedWrenchEnabled = false;

	// try and load the workcell, this will determine if the calibrated wrench is published or not
    WorkCell::Ptr workcell;
    Device::Ptr dev;
    try {
        workcell = RwRos::getWorkCell();
        dev = workcell->findDevice(robotName);
        if (dev == NULL) {
            ROS_INFO_STREAM("Unable to find device "<<robotName<<" in work cell");
        }

        // todo allso
        unbiasedWrenchEnabled = true;
=======
	nh.getParam("robot", devName);
	nh.getParam("calibfile",calibfile);
>>>>>>> .r1128

<<<<<<< .mine
    } catch (const Exception& exp) {
        ROS_WARN_STREAM("Unable to open file " << exp.what());
        ROS_WARN("Unbiased wrench data will not be available!");
    }

    ROS_INFO("WorkCell loaded!");
    ROS_INFO("WorkCell loaded! Starting ur component.");

=======
    // get the workcell
	WorkCell::Ptr _pWorkCell = RwRos::getWorkCell();
	if(_pWorkCell==NULL){
		ROS_ERROR("No workcell added to the parameter server!");
		ROS_BREAK();
	}

	// device
	Device::Ptr _dev = _pWorkCell->findDevice(devName);
	rw::kinematics::State tmpState = _pWorkCell->getDefaultState();
	if (_dev == NULL) {
		ROS_ERROR_STREAM("No device by name " << devName << " Possible devices are:");
		BOOST_FOREACH(rw::models::Device::Ptr dev, _pWorkCell->getDevices()){
			ROS_ERROR_STREAM("dev: " << dev->getName() );
		}
		ROS_BREAK();
	}

	_ftCompensation = new rwhw::FTCompensation(_dev, _pWorkCell->getDefaultState(), calibfile);
	_lastTime = ros::Time::now();
	ros::Subscriber robotState = nh.subscribe("/" + devName + "/RobotState", 1, &handleRobotState);
>>>>>>> .r1128

	try {

		if(!propertyFile.empty()){
			PropertyMap properties = XMLPropertyLoader::load(propertyFile);
			host = properties.get<std::string>("Host","");
			publishRate = properties.get<int>("PublishRate", 400);
		}

		// if there is no host then the sensor cannot start
		if(host.empty()){
			ROS_ERROR("No valid host defined for netft sensor: ");
			ROS_BREAK();
		}

		NetFTLogging netft(host);
		netft.start();

		ros::Publisher pub = nh.advertise<marvin_common::WrenchData>(name + "/WrenchData_raw", 10);
		ros::Publisher pub_calib = nh.advertise<marvin_common::WrenchData>(name + "/WrenchData", 10);

		ros::Rate pub_rate(publishRate);

//        marvinTime.initialize(nh);

        //ros::Duration driverTimeOffset = ros::Time::now() - ros::Time(netft.driverTime());
       // marvinTime.setDriverTime(netft.driverTime());
		int cnt = 0;
		std::vector<Wrench3D> biases;
        while (ros::ok())
        {
        	NetFTLogging::NetFTData data = netft.getAllData();
        	marvin_common::WrenchData wdata, wdata_calib;
        	wdata.header.frame_id = Marvin::FT_ID;
        	wdata.header.stamp = ros::Time::now();

        	wdata_calib.header.frame_id = Marvin::FT_ID;
        	wdata_calib.header.stamp = ros::Time::now();

        //	wdata.timestamp = marvinTime.convertDriverToMarvinTime(data.timestamp);//ros::Time(data.timestamp + driverTimeOffset.toSec()) + syncTimeOffset;

        	Wrench3D ft = data.data;
        	_ft_current = ft;
        	wdata.wrench.force.x = ft.first[0];
        	wdata.wrench.force.y = ft.first[1];
        	wdata.wrench.force.z = ft.first[2];

        	wdata.wrench.torque.x = ft.second[0];
        	wdata.wrench.torque.y = ft.second[1];
        	wdata.wrench.torque.z = ft.second[2];

        	ft = _ftCompensation->getFT();
        	wdata_calib.wrench.force.x = ft.first[0];
        	wdata_calib.wrench.force.y = ft.first[1];
        	wdata_calib.wrench.force.z = ft.first[2];

        	wdata_calib.wrench.torque.x = ft.second[0];
        	wdata_calib.wrench.torque.y = ft.second[1];
        	wdata_calib.wrench.torque.z = ft.second[2];
        	cnt++;
        	if( cnt > 2000  && cnt < 2100){
        		biases.push_back(ft);
        	}
        	if(cnt==2101){
        		ROS_INFO("BIASING");
        		_ftCompensation->updateBias(biases);
        	}

        	pub.publish(wdata);
        	pub_calib.publish( wdata_calib );
        	ros::spinOnce();
        	pub_rate.sleep();
        }

	} catch (const Exception& exp) {
		std::cout<<"Exception = "<<exp.what()<<std::endl;
	}

	return 0;
}

ext: cpp


text:
/**
 * Copyright (c) 2017 Melown Technologies SE
 *
 * 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.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */
#include <atomic>
#include <limits>
#include <sstream>
#include <iomanip>

#include <boost/algorithm/string/split.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/io/ios_state.hpp>

#include "dbglog/dbglog.hpp"
#include "utility/streams.hpp"
#include "utility/uri.hpp"

#include "tileop.hpp"

#include "urltemplate.hpp"

namespace vtslibs { namespace vts {

namespace {

namespace ba = boost::algorithm;

inline unsigned int variableValue(UrlTemplate::Variable what
                                  , const UrlTemplate::Vars &vars)
{
    typedef UrlTemplate::Variable Variable;

    // expand variable
    switch (what) {
    case Variable::lod: return vars.tileId.lod;
    case Variable::x: return vars.tileId.x;
    case Variable::y: return vars.tileId.y;
    case Variable::loclod: return vars.localId.lod;
    case Variable::locx: return vars.localId.x;
    case Variable::locy: return vars.localId.y;
    case Variable::sub: return vars.subMesh;

    case Variable::rf:
    case Variable::srs:
        LOG(warn1) << "Using variable <" << what << "> in numeric context.";
        return std::numeric_limits<unsigned int>::max();
    }
    return 0;
}

struct Expand {
    UrlTemplate::Variable what;
    const UrlTemplate::Vars &vars;
};

inline std::ostream& operator<<(std::ostream &os, const Expand &e)
{
    typedef UrlTemplate::Variable Variable;

    // expand variable
    switch (e.what) {
    case Variable::lod: os << e.vars.tileId.lod; break;
    case Variable::x: os << e.vars.tileId.x; break;
    case Variable::y: os << e.vars.tileId.y; break;
    case Variable::loclod: os << e.vars.localId.lod; break;
    case Variable::locx: os << e.vars.localId.x; break;
    case Variable::locy: os << e.vars.localId.y; break;
    case Variable::sub: os << e.vars.subMesh; break;
    case Variable::srs: os << e.vars.srs; break;
    case Variable::rf: os << e.vars.rf; break;
    }
    return os;
}

inline std::string variableStringValue(UrlTemplate::Variable what
                                       , const UrlTemplate::Vars &vars)
{
    return boost::lexical_cast<std::string>(Expand{what, vars});
}

unsigned int hash(unsigned int x) {
    x = ((x >> 16) ^ x) * 0x45d9f3b;
    x = ((x >> 16) ^ x) * 0x45d9f3b;
    x = (x >> 16) ^ x;
    return x;
}

UrlTemplate::Expander makeAltExpander(const std::string &str
                                   , const std::string::size_type &start
                                   , const std::string::size_type &end)
{
    std::vector<std::string> alternatives;
    auto range(std::make_pair(str.data() + start, str.data() + end));
    ba::split(alternatives, range, ba::is_any_of(", ")
              , ba::token_compress_on);

    if (alternatives.empty()) { return {}; }

    return [alternatives](std::ostream &os,
            const UrlTemplate::Vars &vars)
    {
        unsigned int sum = vars.subMesh
                + vars.tileId.lod
                + vars.tileId.x
                + vars.tileId.y
                + vars.localId.lod
                + vars.localId.x
                + vars.localId.y;
        os << alternatives[hash(sum) % alternatives.size()];
    };
}

UrlTemplate::Expander makePpExpander(const std::string &str
                                  , const std::string::size_type &start
                                  , const std::string::size_type &end
                                  , bool xAxis)
{
    // split args
    std::vector<std::string> args;
    auto range(std::make_pair(str.data() + start, str.data() + end));
    ba::split(args, range, ba::is_any_of(", ")
              , ba::token_compress_on);

    if (args.size() != 2) { return {}; }

    try {
        auto lodVar(boost::lexical_cast<UrlTemplate::Variable>(args[0]));
        auto posVar(boost::lexical_cast<UrlTemplate::Variable>(args[1]));

        if (xAxis) {
            // ppx
            return [lodVar, posVar](std::ostream &os
                                           , const UrlTemplate::Vars &vars)
            {
                vts::Lod lod(variableValue(lodVar, vars));
                unsigned int x(variableValue(posVar, vars));

                boost::io::ios_base_all_saver ps(os);
                os << std::setw(7) << std::hex << std::setfill('0')
                   << (x << (28 - lod));
            };
        }

        // ppy
        return [lodVar, posVar](std::ostream &os
                                       , const UrlTemplate::Vars &vars)
        {
            vts::Lod lod(variableValue(lodVar, vars));
            unsigned int y(variableValue(posVar, vars));

            boost::io::ios_base_all_saver ps(os);
            os << std::setw(7) << std::hex << std::setfill('0')
               << ((1 << 28) - ((y + 1) << (28 - lod)));
        };

    } catch (const boost::bad_lexical_cast&) {}

    return {};
}

UrlTemplate::Expander makeMsDigitExpander(const std::string &str
                                       , const std::string::size_type &start
                                       , const std::string::size_type &end)
{
    // split args
    std::vector<std::string> args;
    auto range(std::make_pair(str.data() + start, str.data() + end));
    ba::split(args, range, ba::is_any_of(", ")
              , ba::token_compress_on);

    if (args.size() != 2) { return {}; }

    try {
        auto xVar(boost::lexical_cast<UrlTemplate::Variable>(args[0]));
        auto yVar(boost::lexical_cast<UrlTemplate::Variable>(args[1]));

        return [xVar, yVar](std::ostream &os, const UrlTemplate::Vars &vars)
        {
            unsigned int x(variableValue(xVar, vars));
            unsigned int y(variableValue(yVar, vars));
            os << (((y & 3) << 1) + (x & 1));
        };

    } catch (const boost::bad_lexical_cast&) {}

    return {};
}

UrlTemplate::Expander makeQuadExpander(const std::string &str
                                    , const std::string::size_type &start
                                    , const std::string::size_type &end)
{
    // split args
    std::vector<std::string> args;
    auto range(std::make_pair(str.data() + start, str.data() + end));
    ba::split(args, range, ba::is_any_of(", ")
              , ba::token_compress_on);

    if (args.size() != 3) { return {}; }

    try {
        auto lodVar(boost::lexical_cast<UrlTemplate::Variable>(args[0]));
        auto xVar(boost::lexical_cast<UrlTemplate::Variable>(args[1]));
        auto yVar(boost::lexical_cast<UrlTemplate::Variable>(args[2]));

        return [lodVar, xVar, yVar](std::ostream &os
                                    , const UrlTemplate::Vars &vars)
        {
            int lod(variableValue(lodVar, vars));
            unsigned int x(variableValue(xVar, vars));
            unsigned int y(variableValue(yVar, vars));

            for (int i(lod - 1); i >= 0; --i) {
                os << (((((y >> i) & 1) << 1) + ((x >> i) & 1)));
            }
        };

    } catch (const boost::bad_lexical_cast&) {}

    return {};
}

UrlTemplate::Expander makeSwitchExpander(const std::string &str
                                         , const std::string::size_type &start
                                         , const std::string::size_type &end)
{
    // split args
    std::vector<std::string> args;
    auto range(std::make_pair(str.data() + start, str.data() + end));
    ba::split(args, range, ba::is_any_of(", ")
              , ba::token_compress_on);

    if (args.size() < 1) { return {}; }

    typedef std::map<std::string, std::string> Mapping;
    Mapping mapping;
    std::string default_;

    for (std::size_t i(1), e(args.size()); i != e; ++i) {
        const auto &arg(args[i]);

        auto colon(arg.find(':'));
        if (colon == std::string::npos) {
            // as-is
            mapping.insert(Mapping::value_type(arg, arg));
            continue;
        }

        Mapping::value_type pair(arg.substr(0, colon), arg.substr(colon + 1));
        if (pair.first == "*") {
            // default value
            default_ = std::move(pair.second);
            continue;
        }

        // store mapping
        mapping.insert(std::move(pair));
    }

    try {
        auto var(boost::lexical_cast<UrlTemplate::Variable>(args[0]));

        return [var, mapping, default_](std::ostream &os
                                        , const UrlTemplate::Vars &vars)
        {
            const auto value(variableStringValue(var, vars));

            const auto fmapping(mapping.find(value));
            if (fmapping == mapping.end()) {
                // not found, use default
                os << ((default_ == "*") ? value : default_);
            } else {
                os << ((fmapping->second == "*") ? value : fmapping->second);
            }
        };

    } catch (const boost::bad_lexical_cast&) {}

    return {};
}

UrlTemplate::Expander makeParamExpander(const std::string &str
                                        , const std::string::size_type &start
                                        , const std::string::size_type &end)
{
    // split args
    std::vector<std::string> args;
    auto range(std::make_pair(str.data() + start, str.data() + end));
    ba::split(args, range, ba::is_any_of(", ")
              , ba::token_compress_on);

    if (args.size() != 1) { return {}; }

    try {
        auto index(boost::lexical_cast<unsigned int>(args[0]));

        return [index](std::ostream &os
                       , const UrlTemplate::Vars &vars)
        {
            if (index >= vars.params.size()) { return; }
            os << utility::urlEncode(vars.params[index]);
        };

    } catch (const boost::bad_lexical_cast&) {}

    return {};
}

UrlTemplate::Expander parseFunction(const std::string &str
                                 , const std::string::size_type &start
                                 , const std::string::size_type &end)
{
    auto open(str.find('(', start));
    if (open >= end) {
        return {};
    }

    auto close(str.find(')', open));
    if (close >= end) {
        return {};
    }

    auto name(str.substr(start, open - start));

    if (name == "alt") {
        return makeAltExpander(str, open + 1, close);
    }

    if (name == "ppx") {
        return makePpExpander(str, open + 1, close, true);
    }

    if (name == "ppy") {
        return makePpExpander(str, open + 1, close, false);
    }

    if (name == "ms_digit") {
        return makeMsDigitExpander(str, open + 1, close);
    }

    if (name == "quad") {
        return makeQuadExpander(str, open + 1, close);
    }

    if (name == "switch") {
        return makeSwitchExpander(str, open + 1, close);
    }

    if (name == "param") {
        return makeParamExpander(str, open + 1, close);
    }

    return {};
}

} // namespace

void UrlTemplate::parse(const std::string &str)
{
    tokens_.clear();

    std::string::size_type index(0);
    std::string::size_type end(str.size());
    while (index < end) {
        // find opening brace
        auto open(str.find('{', index));
        if (open != index) {
            // some verbatim text
            tokens_.emplace_back(str.substr(index, (open - index)));
        }

        // eof?
        if (open == std::string::npos) { break; }

        // find close brace
        auto close(str.find('}', open + 1));
        if (close == std::string::npos) {
            // unclosed, consume everything
            LOG(warn2)
                << "Invalid expansion string <" << str
                << ">: unmatched '{'.";
            tokens_.emplace_back(str.substr(open));
            break;
        }

        // we have something inside { braces }

        // try a variable
        Variable result = Variable();
        if (boost::conversion::try_lexical_convert(
            str.substr(open + 1, close - open - 1), result)) {
            tokens_.emplace_back(result);
        } else {
            // not a variable; try a function
            if (auto func = parseFunction(str, open + 1, close)) {
                tokens_.emplace_back(str.substr(open + 1, close - open - 1)
                                     , func);
            } else {
                // emplace the name of the substituent, excluding the braces
                tokens_.emplace_back(str.substr(open + 1, close - open - 1));

                // ignore viewspec, otherwise issue a warning
                if (tokens_.back().value != "viewspec") {
                    LOG(warn2)
                        << "Unprocessed directive <" << tokens_.back().value
                        << ">.";
                }
            }
        }

        // next round
        index = close + 1;
    }
}

void UrlTemplate::expand(std::ostream &os, const Vars &vars) const
{
    for (const auto &token : tokens_) {
        token.expand(os, vars);
    }
}

void UrlTemplate::Token::expand(std::ostream &os, const Vars &vars)
    const
{
    if (expander) {
        // expander callback
        expander(os, vars);
        return;
    }

    if (variable) {
        // expand variable
        os << Expand{*variable, vars};
        return;
    }

    // plain value
    os << value;
}

void UrlTemplate::dump(std::ostream &os) const
{
    for (const auto &token : tokens_) {
        token.dump(os);
    }
}

void UrlTemplate::Token::dump(std::ostream &os)
    const
{
    if (expander) {
        os << "[expander:" << value << "]";
        return;
    }

    if (variable) {
        os << "[var:" << *variable << "]";
        return;
    }

    // plain value
    os << "[text:" << value << "]";
}

} } // namespace vtslibs::vts

max_stars_repo_name: melowntech/vts-libs

max_issues_repo_name: HanochZhu/vts-browser-unity-plugin

max_forks_repo_name: HanochZhu/vts-browser-unity-plugin


text:
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012-2015, all rights reserved,
//      Christian Schulte zu Berge <christian.szb@in.tum.de>
//      Chair for Computer Aided Medical Procedures
//      Technische Universitaet Muenchen
//      Boltzmannstr. 3, 85748 Garching b. Muenchen, Germany
//
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
// except in compliance with the License. You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the
// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
// either express or implied. See the License for the specific language governing permissions
// and limitations under the License.
//
// ================================================================================================

#include "tensoranalyzer.h"

#include "cgt/cgt_math.h"
#include "cgt/logmanager.h"

#include <Eigen/Eigenvalues>
#include <tbb/tbb.h>
#include <tbb/atomic.h>

#include "core/datastructures/imagedata.h"
#include "core/tools/stringutils.h"

namespace campvis {

    static const GenericOption<TensorAnalyzer::DegeneratedEvHandling> handlingModes[4] = {
        GenericOption<TensorAnalyzer::DegeneratedEvHandling>("none", "None", TensorAnalyzer::NONE),
        GenericOption<TensorAnalyzer::DegeneratedEvHandling>("mask", "Mask", TensorAnalyzer::MASK),
        GenericOption<TensorAnalyzer::DegeneratedEvHandling>("invert", "Invert", TensorAnalyzer::INVERT),
        GenericOption<TensorAnalyzer::DegeneratedEvHandling>("shift", "Shift", TensorAnalyzer::SHIFT)
    };

    static const GenericOption<std::string> measurementOptions[15] = {
        GenericOption<std::string>("Disabled", "Disabled"),
        GenericOption<std::string>("EigenValue1", "Eigenvalue 1"),
        GenericOption<std::string>("EigenValue2", "Eigenvalue 2"),
        GenericOption<std::string>("EigenValue3", "Eigenvalue 3"),
        GenericOption<std::string>("MainEigenvector", "Main Eigenvector"),
        GenericOption<std::string>("VolumeRatio", "Volume Ratio"),
        GenericOption<std::string>("FractionalAnisotropy", "Fractional Anisotropy"),
        GenericOption<std::string>("RelativeAnisotropy", "Relative Anisotropy"),
        GenericOption<std::string>("MeanDiffusivity", "Mean Diffusivity"),
        GenericOption<std::string>("Trace", "Trace"),
        GenericOption<std::string>("AxialDiffusivity", "Axial Diffusivity"),
        GenericOption<std::string>("RadialDiffusivity", "Radial Diffusivity"),
        GenericOption<std::string>("LinearAnisotropy", "Linear Anisotropy"),
        GenericOption<std::string>("PlanarAnisotropy", "Planar Anisotropy"),
        GenericOption<std::string>("Isotropy", "Isotropy")
    };

    const std::string TensorAnalyzer::loggerCat_ = "CAMPVis.modules.classification.TensorAnalyzer";

    TensorAnalyzer::OutputPropertyPair::OutputPropertyPair(size_t index)
        : _imageId("OutputId" + StringUtils::toString(index), "Output " + StringUtils::toString(index) + " Image", "TensorAnalyzer.output" + StringUtils::toString(index), DataNameProperty::WRITE)
        , _imageType("OutputType" + StringUtils::toString(index), "Output " + StringUtils::toString(index) + " Image Type", measurementOptions, 15)
    {}

    TensorAnalyzer::TensorAnalyzer()
        : AbstractProcessor()
        , p_inputImage("InputImage", "Input Tensor Image", "tensors", DataNameProperty::READ)
        , p_evalsImage("EvalsImage", "Output Eigenvalues Image", "TensorAnalyzer.eigenvalues", DataNameProperty::WRITE)
        , p_evecsImage("EvecsImage", "Output Eigenvectors Image", "TensorAnalyzer.eigenvectors", DataNameProperty::WRITE)
        , p_degeneratedHandling("DegeneratedHandling", "Handling of Degenerated Tensors", handlingModes, 4)
        , p_maskMixedTensors("MaskMixedTensors", "Mask Mixed Tensors", true)
        , p_addOutputButton("AddOutputButton", "Add Output")
        , _eigenvalues(0)
        , _eigenvectors(0)
    {
        addProperty(p_inputImage, INVALID_RESULT | EIGENSYSTEM_INVALID);
        addProperty(p_evalsImage);
        addProperty(p_evecsImage);
        addProperty(p_degeneratedHandling);
        addProperty(p_maskMixedTensors);
        addProperty(p_addOutputButton, VALID);
        addOutput();

        p_addOutputButton.s_clicked.connect(this, &TensorAnalyzer::addOutput);
    }

    TensorAnalyzer::~TensorAnalyzer() {

    }

    void TensorAnalyzer::updateResult(DataContainer& data) {
        if (getInvalidationLevel() & EIGENSYSTEM_INVALID) {
            computeEigensystem(data);
        }

        if (_eigenvalues.getData() != 0 && _eigenvectors.getData() != 0) {
            for (size_t i = 0; i < p_outputProperties.size(); ++i) {
                computeOutput(data, i);
            }
        }
        else {
            LERROR("Could not compute Eigensystem");
        }
    }

    void TensorAnalyzer::computeEigensystem(DataContainer& data) {
        // get input tensor image
        GenericImageRepresentationLocal<float, 6>::ScopedRepresentation input(data, p_inputImage.getValue());

        if (input != 0) {
            // create output images for eigenvalues (stored as vec3) and eigenvectors (stored as mat3)
            ImageData* evals = new ImageData(input->getDimensionality(), input->getSize(), 3);
            GenericImageRepresentationLocal<float, 3>* evalRep = GenericImageRepresentationLocal<float, 3>::create(evals, 0);

            ImageData* evecs = new ImageData(input->getDimensionality(), input->getSize(), 9);
            GenericImageRepresentationLocal<float, 9>* evecRep = GenericImageRepresentationLocal<float, 9>::create(evecs, 0);

            tbb::atomic<size_t> countDiscarded;
            countDiscarded = 0;
            DegeneratedEvHandling evh = p_degeneratedHandling.getOptionValue();

            // perform eigen decomposition in parallel
            tbb::parallel_for(tbb::blocked_range<size_t>(0, input->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const Tensor2<float>& t = input->getElement(i);

                    if (t.Dxx == 0 && t.Dxy == 0 && t.Dxz == 0 && t.Dyy == 0 && t.Dyz == 0 && t.Dzz == 0) {
                        evalRep->setElement(i, cgt::vec3::zero);
                        evecRep->setElement(i, cgt::mat3::zero);
                        continue;
                    }

                    // use Eigen library for eigensystem calculation
                    Eigen::Matrix3f tensorMatrix;
                    tensorMatrix << t.Dxx, t.Dxy, t.Dxz, t.Dxy, t.Dyy, t.Dyz, t.Dxz, t.Dyz, t.Dzz;

                    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> solver(tensorMatrix);

                    // gather data from eigen library
                    const Eigen::Matrix<float, 3, 1>& values = solver.eigenvalues();
                    const Eigen::Matrix3f& vectors = solver.eigenvectors();

                    // kill NaN values
                    if (cgt::isNaN(values(0)) || cgt::isNaN(values(1)) || cgt::isNaN(values(2))) {
                        evalRep->setElement(i, cgt::vec3::zero);
                        evecRep->setElement(i, cgt::mat3::zero);
                        continue;
                    }

                    // determine order of eigenvalues (i.e. sort'em)
                    int high, mid, low;
                    bool negative = false;
                    if (values(0) > values(1)) {
                        if (values(1) > values(2)) {
                            high = 0; mid = 1; low = 2;
                        }
                        else if (values(0) > values(2)) {
                            high = 0; mid = 2; low = 1;
                        }
                        else {
                            high = 2; mid = 0; low = 1;
                        }
                    }
                    else {
                        if (values(0) > values(2)) {
                            high = 1; mid = 0; low = 2;
                        }
                        else if (values(1) > values(2)) {
                            high = 1; mid = 2; low = 0;
                        }
                        else {
                            high = 2; mid = 1; low = 0;
                        }
                    }

                    // perform handling of degenerated tensors (i.e. negative/mixed eigenvalues):
                    if (values(0) < 0 && values(1) < 0 && values(2) < 0) {
                        negative = true;
                        if (evh == INVERT)
                            std::swap(high, low);
                    }
                    else if (p_maskMixedTensors.getValue() && (values(0) < 0 || values(1) < 0 || values(2) < 0)) {
                        // We assume that either all eigenvalues are positive or all eigenvalues are negative.
                        // If we encounter both positive and negative eigenvalues this must be due to severe noise
                        // (e.g. area outside brain) so it is reasonable to discard these voxels.
                        ++countDiscarded;
                        evalRep->setElement(i, cgt::vec3::zero);
                        evecRep->setElement(i, cgt::mat3::zero);
                        continue;
                    }

                    if (negative && evh == MASK) {
                        evalRep->setElement(i, cgt::vec3::zero);
                        evecRep->setElement(i, cgt::mat3::zero);
                    }
                    else {
                        evalRep->setElement(i, cgt::vec3(values(high), values(mid), values(low)));
                        evecRep->setElement(i, cgt::mat3(vectors(0, high), vectors(1, high), vectors(2, high),
                            vectors(0, mid),  vectors(1, mid),  vectors(2, mid),
                            vectors(0, low),  vectors(1, low),  vectors(2, low)));

                        if (negative && evh == SHIFT) {
                            evalRep->getElement(i) -= (evalRep->getElement(i).x + evalRep->getElement(i).z);
                        }
                    }
                }
            });

            // write results to DataContainer and also cache them in local members
            _eigenvalues = data.addData(p_evalsImage.getValue(), evals);
            _eigenvectors = data.addData(p_evecsImage.getValue(), evecs);
        }
        else {
            LDEBUG("No suitable input image found.");
        }

        validate(EIGENSYSTEM_INVALID);
    }

    void TensorAnalyzer::computeOutput(DataContainer& data, size_t index) {
        // sanity check
        if (index >= p_outputProperties.size()) {
            LERROR("Index out of bounds while computing output #" << index);
            return;
        }

        // gather eigensystem
        const GenericImageRepresentationLocal<float, 3>* evalRep = 0;
        const GenericImageRepresentationLocal<float, 9>* evecRep = 0;
        if (_eigenvalues.getData() != 0 && _eigenvectors.getData() != 0) {
            evalRep = static_cast<const ImageData*>(_eigenvalues.getData())->getRepresentation< GenericImageRepresentationLocal<float, 3> >(false);
            evecRep = static_cast<const ImageData*>(_eigenvectors.getData())->getRepresentation< GenericImageRepresentationLocal<float, 9> >(false);
        }
        if (evalRep == 0 || evecRep == 0) {
            LERROR("Could not compute output, no eigensystem present.");
            return;
        }

        OutputPropertyPair* opp = p_outputProperties[index];
        const std::string& type = opp->_imageType.getOptionValue();

        if (type == "Disabled") {
            return;
        }
        else if (type == "EigenValue1") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    float ev = vals.x;
                    if (ev == 0.0f || cgt::isNaN(ev))
                        output->setElement(i, 0.0f);
                    else
                        output->setElement(i, ev);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "EigenValue2") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    float ev = vals.y;
                    if (ev == 0.0f || cgt::isNaN(ev))
                        output->setElement(i, 0.0f);
                    else
                        output->setElement(i, ev);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "EigenValue3") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    float ev = vals.z;
                    if (ev == 0.0f || cgt::isNaN(ev))
                        output->setElement(i, 0.0f);
                    else
                        output->setElement(i, ev);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }

        else if (type == "MainEigenvector") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 3);
            GenericImageRepresentationLocal<float, 3>* output = GenericImageRepresentationLocal<float, 3>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    output->setElement(i, evecRep->getElement(i)[0]);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }

        // = Anisotropy Measures: =========================================================================
        else if (type == "VolumeRatio") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = (vals.x * vals.y * vals.z) / pow((vals.x + vals.y + vals.z)/3.f, 3);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "FractionalAnisotropy") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);
                    const float root = sqrt(.5f);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = root * sqrt((vals.x-vals.y)*(vals.x-vals.y) + (vals.y-vals.z)*(vals.y-vals.z) + (vals.z-vals.x)*(vals.z-vals.x)) / sqrt(vals.x*vals.x + vals.y*vals.y + vals.z*vals.z);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "RelativeAnisotropy") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);
                    const float root = sqrt(.5f);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = root * sqrt((vals.x-vals.y)*(vals.x-vals.y) + (vals.y-vals.z)*(vals.y-vals.z) + (vals.z-vals.x)*(vals.z-vals.x)) / (vals.x + vals.y + vals.z);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "MeanDiffusivity") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = (vals.x + vals.y + vals.z) / 3;
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "Trace") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = vals.x + vals.y + vals.z;
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "AxialDiffusivity") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = vals.x;
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "RadialDiffusivity") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = (vals.y + vals.z) / 2;
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "LinearAnisotropy") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = (vals.x - vals.y) / (vals.x + vals.y + vals.z);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "PlanarAnisotropy") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = 2.f*(vals.y - vals.z) / (vals.x + vals.y + vals.z);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }
        else if (type == "Isotropy") {
            ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
            GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
            tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
                for (size_t i = range.begin(); i != range.end(); ++i) {
                    const cgt::vec3& vals = evalRep->getElement(i);

                    if (vals == cgt::vec3::zero || cgt::isNaN(vals))
                        output->getElement(i) = 0;
                    else
                        output->getElement(i) = (3.f*vals.z) / (vals.x + vals.y + vals.z);
                }
            });
            data.addData(opp->_imageId.getValue(), id);
        }

    }

    void TensorAnalyzer::addOutput() {
        OutputPropertyPair* opp = new OutputPropertyPair(p_outputProperties.size() + 1);
        p_outputProperties.push_back(opp);
        addProperty(opp->_imageId);
        addProperty(opp->_imageType);
    }

    void TensorAnalyzer::deinit() {
        _eigenvectors = DataHandle(0);
        _eigenvalues = DataHandle(0);
        AbstractProcessor::deinit();
    }

}

max_stars_repo_name: hmsgit/campvis

max_issues_repo_name: hmsgit/campvis

max_forks_repo_name: hmsgit/campvis


text:
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/program_options.hpp>
#include <boost/regex.hpp>
#include <boost/stacktrace.hpp>
#include <boost/system/error_code.hpp>
#include <fstream>
#include <iostream>
#include <memory>
#include <sstream>

#include "QStruct.h"
#include "compiler.hpp"
#include "config.h"  // Add an automatically generated configuration file

using namespace boost;

using boost::lexical_cast;

extern std::string parser(std::string sInputFile);
extern std::string storeParseResult(std::string sOutputFile);

extern qTree coreInstance;

int main(int argc, char *argv[]) {
  // Clarification: When gcc has been upgraded to 9.x version some tests fails.
  // Bug appear when data are passing to program via script .sh
  // additional 13 (\r) character was append - this code normalize argv list.
  // C99: The parameters argc and argv and the strings pointed to by the argv
  // array shall be modifiable by the program, and retain their last-stored
  // values between program startup and program termination.
  for (int i = 0; i < argc; i++) {
    auto len = strlen(argv[i]);
    if (len > 0)
      if (argv[i][len - 1] == 13) argv[i][len - 1] = 0;
  }
  try {
    std::unique_ptr<std::ostream> p_ofs;
    namespace po = boost::program_options;
    std::string sOutputFile;
    std::string sInputFile;
    po::options_description desc("Avaiable options");
    desc.add_options()("help,h", "Show program options")(
        "queryfile,q", po::value<std::string>(&sInputFile), "query set file")(
        "outfile,o",
        po::value<std::string>(&sOutputFile)->default_value("query.qry"),
        "output file")("dumpcross,d",
                       "dump diagnostic cross compilation forms");
    po::positional_options_description
        p;  // Assume that infile is the first option
    p.add("queryfile", -1);
    po::variables_map vm;
    po::store(
        po::command_line_parser(argc, argv).options(desc).positional(p).run(),
        vm);
    po::notify(vm);
    if (vm.count("help")) {
      std::cout << desc;
      std::cout << CONFIG_LINE;
      return system::errc::success;
    }
    if (vm.count("queryfile") == 0) {
      std::cout << argv[0] << ": fatal error: no input file" << std::endl;
      std::cout << "compilation terminated." << std::endl;
      return EPERM;  // ERROR defined in errno-base.h
    }
    std::string sQueryFile = vm["queryfile"].as<std::string>();
    std::string sOutFile = vm["outfile"].as<std::string>();
    std::string sSubsetFile = "query.sub";
    // override defaults if filename is matching regexp
    {
      boost::regex filenameRe("(\\w*).(\\w*)");  // filename.extension
      boost::cmatch what;
      if (regex_match(sQueryFile.c_str(), what, filenameRe)) {
        std::string filenameNoExt = std::string(what[1]);
        sSubsetFile = filenameNoExt + ".sub";
      }
    }
    std::cout << "Input file:" << sQueryFile << std::endl;
    std::cout << "Compile result:" << parser(sQueryFile) << std::endl;
    //
    // Main Algorithm body
    //
    std::istringstream iss(dumpInstance(""), std::ios::binary);
    boost::archive::text_iarchive ia(iss);
    ia >> coreInstance;
    if (coreInstance.size() == 0)
      throw std::out_of_range("No queries to process found");
    if (vm.count("dumpcross")) {
      std::string sOutFileLog1 = vm["outfile"].as<std::string>() + ".lg1";
      dumpInstance(sOutFileLog1);
    }
    int dAfterParserSize((int)coreInstance.size());
    std::string response;
    response = simplifyLProgram();
    if (vm.count("dumpcross")) {
      std::string sOutFileLog2 = vm["outfile"].as<std::string>() + ".lg2";
      dumpInstance(sOutFileLog2);
    }
    response = prepareFields();
    assert(response == "OK");
    if (vm.count("dumpcross")) {
      std::string sOutFileLog3 = vm["outfile"].as<std::string>() + ".lg3";
      dumpInstance(sOutFileLog3);
    }
    response = intervalCounter();
    assert(response == "OK");
    response = convertReferences();
    assert(response == "OK");
    response = replicateIDX();
    assert(response == "OK");
    int dAfterCompilingSize((int)coreInstance.size());
    dumpInstance(sOutFile);
  } catch (std::exception &e) {
    std::cerr << e.what() << "\n";
    // cerr << boost::stacktrace::stacktrace() << endl ;
    return system::errc::interrupted;
  }
  return system::errc::success;
}

max_stars_repo_name: michalwidera/abracadabradb

max_issues_repo_name: michalwidera/abracadabradb

max_forks_repo_name: michalwidera/abracadabradb


text:
/**
*  @file    mdm_RunTools_madym_MakeXtr.cxx
*  @brief   Implementation of mdm_RunTools_madym_MakeXtr class
*
*  Original author MA Berks 24 Oct 2018
*  (c) Copyright QBI, University of Manchester 2020
*/

#ifndef MDM_API_EXPORTS
#define MDM_API_EXPORTS
#endif // !MDM_API_EXPORTS

#include "mdm_RunTools_madym_MakeXtr.h"

#include <madym/image_io/xtr/mdm_XtrFormat.h>
#include <madym/image_io/mdm_ImageIO.h>
#include <madym/mdm_Image3D.h>
#include <madym/mdm_SequenceNames.h>

#include <madym/mdm_exception.h>

#include <string>
#include <vector>

#include <boost/filesystem.hpp>
#include <boost/format.hpp>
#include <boost/algorithm/string.hpp>

#include <algorithm>
#include <numeric>

namespace fs = boost::filesystem;

//
MDM_API mdm_RunTools_madym_MakeXtr::mdm_RunTools_madym_MakeXtr()
{
}

MDM_API mdm_RunTools_madym_MakeXtr::~mdm_RunTools_madym_MakeXtr()
{
}

//
MDM_API void mdm_RunTools_madym_MakeXtr::run()
{
  //Set curent working dir
  set_up_cwd();

  if (options_.makeT1Inputs())
    makeT1InputXtr();

  if (options_.makeDyn())
    makeDynamicXtr();

  mdm_ProgramLogger::logProgramMessage("Finished processing!");
}

//
MDM_API int mdm_RunTools_madym_MakeXtr::parseInputs(int argc, const char *argv[])
{
  po::options_description cmdline_options("madym_MakeXtr options");
  po::options_description config_options("madym_MakeXtr config options");

  options_parser_.add_option(cmdline_options, options_.configFile);
  options_parser_.add_option(cmdline_options, options_.dataDir);

	//General output options_
  options_parser_.add_option(config_options, options_.dynDir);
  options_parser_.add_option(config_options, options_.dynName);
  options_parser_.add_option(config_options, options_.sequenceFormat);
  options_parser_.add_option(config_options, options_.sequenceStart);
  options_parser_.add_option(config_options, options_.sequenceStep);
  options_parser_.add_option(config_options, options_.T1inputNames);
  options_parser_.add_option(config_options, options_.T1Dir);
  options_parser_.add_option(config_options, options_.nDyns);
  options_parser_.add_option(config_options, options_.T1method);

  // XTR options
  options_parser_.add_option(config_options, options_.makeDyn);
  options_parser_.add_option(config_options, options_.makeT1Inputs);
  options_parser_.add_option(config_options, options_.dynTimesFile);
  options_parser_.add_option(config_options, options_.temporalResolution);
  options_parser_.add_option(config_options, options_.TR);
  options_parser_.add_option(config_options, options_.FA);
  options_parser_.add_option(config_options, options_.VFAs);

  return options_parser_.parseInputs(
    cmdline_options,
    config_options,
    options_.configFile(),
    who(),
    argc, argv);
}

MDM_API std::string mdm_RunTools_madym_MakeXtr::who() const
{
	return "madym_DicomConvert";
}

//*******************************************************************************
// Private:
//*******************************************************************************

//---------------------------------------------------------------------
void mdm_RunTools_madym_MakeXtr::readDynamicTimes()
{
  //Try and open the file and read in the times
  std::ifstream dynTimesStream(
    fs::absolute(options_.dynTimesFile()).string(), std::ios::in);

  if (!dynTimesStream.is_open())
  {
    throw mdm_exception(__func__, "error opening dynamic times file, Check it exists");
  }

  dynamicTimes_.clear();
  for (int i = 0; i < options_.nDyns(); i++)
  {
    double t;
    dynTimesStream >> t;
    dynamicTimes_.push_back(t);
  }
}

//---------------------------------------------------------------------
double mdm_RunTools_madym_MakeXtr::getDynamicTime(int dynNum)
{
  double secs;
  if (dynamicTimes_.size() == options_.nDyns())
    secs = 60 * dynamicTimes_[dynNum];
  else
    secs = dynNum * options_.temporalResolution();

  return mdm_Image3D::secsToTimestamp(secs);

}

//---------------------------------------------------------------------
void mdm_RunTools_madym_MakeXtr::makeT1InputXtr()
{
  //Parse T1 method from string, will abort if method type not recognised
  auto methodType = mdm_T1MethodGenerator::parseMethodName(options_.T1method(), false);

  switch (methodType)
  {
  case mdm_T1MethodGenerator::VFA: makeVFAXtr(); break;
  default:
    throw mdm_exception(__func__, "T1 method not recognised");
  }
}

//---------------------------------------------------------------------
void mdm_RunTools_madym_MakeXtr::makeVFAXtr()
{
  //Check required inputs are set
  auto nInputs = options_.T1inputNames().size();

  if (!nInputs)
    throw mdm_exception(__func__, "T1 input names must be set");

  if (nInputs != options_.VFAs().size())
    throw mdm_exception(__func__,
      boost::format("Number of elements in VFAs (%1%) does not match number input names (%2%)")
      % options_.VFAs().size() % nInputs);

  auto TR = options_.TR();

  if (!TR)
    throw mdm_exception(__func__, "TR must not be zero for VFA T1 mapping");

  for (size_t i_t1 = 0; i_t1 < nInputs; i_t1++)
  {
    auto T1path = fs::absolute(
      fs::path(options_.T1Dir()) / options_.T1inputNames()[i_t1]);
    T1path.replace_extension();

    auto FA = options_.VFAs()[i_t1];

    //Set up mean image
    mdm_Image3D img;
    img.info().flipAngle.setValue(FA);
    img.info().TR.setValue(TR);

    mdm_XtrFormat::writeAnalyzeXtr(T1path.string(), img, mdm_XtrFormat::NEW_XTR);
    mdm_ProgramLogger::logProgramMessage("Created T1 input XTR file " + T1path.string() + ".xtr");
  }
}

//---------------------------------------------------------------------
void mdm_RunTools_madym_MakeXtr::makeDynamicXtr()
{
  //Check required inputs are set
  if (!options_.nDyns())
    throw mdm_exception(__func__, "nDyns must not be zero for dynamic series");

  if (!options_.dynTimesFile().empty())
    readDynamicTimes();

  else if (!options_.temporalResolution())
    throw mdm_exception(__func__, "Either temporalResolution or a dynamic times file must be set");

  //Get basename of dynamic volumes from input options
  fs::path dynPath = fs::absolute(fs::path(options_.dynDir()) / options_.dynName());
  std::string dynPrefix = dynPath.filename().string();
  std::string dynBasePath = dynPath.remove_filename().string();

  auto TR = options_.TR();
  auto FA = options_.FA();

  if (!TR)
    throw mdm_exception(__func__, "TR must not be zero for dynamic series");

  if (!FA)
    throw mdm_exception(__func__, "FA must not be zero for dynamic series");

  //Get write format from options
  for (int i_dyn = 0; i_dyn < options_.nDyns(); i_dyn++)
  {
    mdm_Image3D img;
    img.info().flipAngle.setValue(FA);
    img.info().TR.setValue(TR);

    double acquisitionTime = getDynamicTime(i_dyn);
    img.setTimeStampFromDoubleStr(acquisitionTime);

    auto outputName = mdm_SequenceNames::makeSequenceFilename(
      dynBasePath, dynPrefix, i_dyn + 1, options_.sequenceFormat(),
      options_.sequenceStart(), options_.sequenceStep());

    mdm_XtrFormat::writeAnalyzeXtr(outputName, img, mdm_XtrFormat::NEW_XTR);
    mdm_ProgramLogger::logProgramMessage("Created dynamic XTR file " + outputName + ".xtr");
  }
}

ext: cxx


text:
#include "ReadVideo.h"

#include <iostream>
#include <boost/exception/all.hpp>
#include <boost/log/trivial.hpp>
#include <opencv2/imgproc.hpp>

ReadVideo::ReadVideo() {

}

ReadVideo::~ReadVideo() {
    this->result.reset();
}

void ReadVideo::Start() {
    this->status = TaskStatus::Running;

    this->result = std::make_shared<ReadVideoResult>();

    if(this->readVideo()) {
        this->status = TaskStatus::Finished;
        return;
    } else {
        this->status = TaskStatus::Errored;
    }
}

bool ReadVideo::readVideo() {

    // open input file context
    AVFormatContext* formatContext = nullptr;

    int ret = avformat_open_input(
        &formatContext,
        this->videoFilePath.c_str(),
        nullptr,
        nullptr);

    if (ret < 0) {
        BOOST_LOG_TRIVIAL(error) << "Failed to avformat_open_input: \"" << this->videoFilePath.c_str() << "\" Code:" << ret;
        return false;
    }

    ret = avformat_find_stream_info(formatContext, nullptr);
    if (ret < 0) {
        BOOST_LOG_TRIVIAL(error) << "Failed to avformat_find_stream_info: Code:" << ret;
        return false;
    }

    AVCodec* vcodec = nullptr;
    int bestStream = av_find_best_stream(formatContext, AVMEDIA_TYPE_VIDEO, -1, -1, &vcodec, 0);

    if (bestStream < 0) {
        BOOST_LOG_TRIVIAL(error) << "Fail to av_find_best_stream: Code:" << bestStream;
        return false;
    }

    // Select best stream
    AVStream* vstrm  = formatContext->streams[bestStream];

    AVCodecContext *codecContext = avcodec_alloc_context3(vcodec);
    avcodec_parameters_to_context(codecContext, vstrm->codecpar);

    ret = avcodec_open2(codecContext, vcodec, nullptr);
    if (ret < 0) {
        BOOST_LOG_TRIVIAL(error) << "Failed to avcodec_open2: Code:" << ret;
        return false;
    }

    BOOST_LOG_TRIVIAL(info)
        << "File       : " << this->videoFilePath << "\n"
        << "Format     : " << formatContext->iformat->name << "\n"
        << "Codec      : " << vcodec->name << "\n"
        << "Size       : " << codecContext->width << 'x' << codecContext->height << "\n"
        << "FPS        : " << av_q2d(codecContext->framerate) << " [fps]\n"
        << "Length     : " << av_rescale_q(vstrm->duration, vstrm->time_base, {1,1000}) / 1000. << " [sec]\n"
        << "Pix format : " << av_get_pix_fmt_name(codecContext->pix_fmt) << "\n"
        << "Framecount : " << vstrm->nb_frames << "\n";

    SwsContext*  sws_ctx =
        sws_getContext(   // [13]
            codecContext->width,
            codecContext->height,
            codecContext->pix_fmt,
            codecContext->width,
            codecContext->height,
            //AVPixelFormat::AV_PIX_FMT_BGR24,   // sws_scale destination color scheme
            AVPixelFormat::AV_PIX_FMT_RGB24,
            SWS_BILINEAR,
            nullptr,
            nullptr,
            nullptr
            );

    // Convert frames to opencv Mat
    AVPacket pkt;
    AVFrame* decFrame = av_frame_alloc();

    int idx = 0;
    while(av_read_frame(formatContext, &pkt)  >= 0) {

        if (pkt.stream_index != bestStream)
        {
            av_packet_unref(&pkt);
            continue;
        }

        int ret = avcodec_send_packet(codecContext, &pkt);
        if (ret < 0)
        {

            // could not send packet for decoding
            BOOST_LOG_TRIVIAL(error) << "Error sending packet for decoding";

            return false;
        }

        while(ret >= 0) {
            ret = avcodec_receive_frame(codecContext, decFrame);
            if (ret == AVERROR(EAGAIN) || ret == AVERROR_EOF)
            {
                // EOF exit loop
                break;
            }
            else if (ret < 0)
            {
                // could not decode packet
                BOOST_LOG_TRIVIAL(error) << "Error while decoding.\n";

                // exit with error
                return false;
            }

            if(idx >= this->startFrame && idx <= this->endFrame) { // 0 - 99 resulted frame count should be 100 / [99; 99] -> 1 frame / [0; 1] -> 2 frames

                cv::Mat image = cv::Mat(decFrame->height, decFrame->width, CV_8UC3);
                int cvLinesizes[1];
                cvLinesizes[0] = image.step1();

                sws_scale(sws_ctx, decFrame->data, decFrame->linesize, 0, decFrame->height, &image.data, cvLinesizes);

                this->result->VideoFrameCollection.push_back(image.clone());
                idx++;

            } else {
                idx++;
                break;
            }

        }

        av_packet_unref(&pkt);
    }

    av_frame_free(&decFrame);
    sws_freeContext(sws_ctx);
    avcodec_close(codecContext);
    avformat_close_input(&formatContext);

    this->result->Successful = true;

    return true;
}

void ReadVideo::Stop() {
    this->status = TaskStatus::Stopped;
}

void ReadVideo::Pause() {
    this->status = TaskStatus::Paused;
}

TaskStatus ReadVideo::GetStatus() {
    return this->status;
}

std::shared_ptr<void> ReadVideo::GetResult() {
    return this->result;
}

void ReadVideo::SetVideoFilePath(std::string videoFilePath) {
    this->videoFilePath = videoFilePath;
}

void ReadVideo::SetReadRange(int startFrame, int endFrame) {

    this->startFrame = startFrame;
    this->endFrame   = endFrame;

}

max_stars_repo_name: CoreTrackProject/Prototype

max_issues_repo_name: CoreTrackProject/Prototype

max_forks_repo_name: CoreTrackProject/Prototype


text:
//
// Created by enrico on 19.03.21.
//

#ifndef TAGPROCESSORLIST_HPP
#define TAGPROCESSORLIST_HPP

#include "OrthancPluginCppWrapper.h"

#include <dcmtk/dcmdata/dctagkey.h>
#include <dcmtk/dcmdata/dcdatset.h>
#include <dcmtk/dcmdata/dcdict.h>
#include <dcmtk/dcmdata/dcdicent.h>
#include <boost/optional.hpp>
#include <boost/filesystem.hpp>
#include <json/json.h>
#include <regex>
#include <list>

struct TagProcessingGenerator {
	DcmTagKey tag;
	boost::optional<std::regex> regex;
	std::string replacement;
	std::string generate(DcmDataset& dset)const;
};
struct TagProcessor{
	DcmTagKey tag;
	boost::optional<std::regex> mask;
	std::list<TagProcessingGenerator> generators;
};

class TagProcessorList: public std::list<TagProcessor>
{
public:
	TagProcessorList(const OrthancPlugins::OrthancConfiguration &configuration);
	void FixTags(DcmDataset *dset)
	{
		for(const auto &proc : *this){
			OFString original;
			auto erg=dset->findAndGetOFString(proc.tag,original);
			if(proc.mask && !std::regex_match(original.begin(), original.end(), *proc.mask)){
				continue;//existing mask did not match, don't process this Tag
			}

			if(proc.generators.empty()){ // is a simple delete
				dset->findAndDeleteElement(proc.tag,true,true);
			} else {
				std::string replacement;
				for (const auto &generator:proc.generators) {
					replacement+=generator.generate(*dset);
				}
				if(!dset->putAndInsertOFStringArray(proc.tag,replacement.c_str()).good())
					OrthancPlugins::LogError(std::string("Failed set ") + DcmTag(proc.tag).getTagName() + " to " + replacement);
			}
		}
	}

};

#endif //TAGPROCESSORLIST_HPP

ext: hpp

max_stars_repo_name: DerOrfa/Orthanc-plugins-cpp

max_issues_repo_name: DerOrfa/Orthanc-plugins-cpp

max_forks_repo_name: DerOrfa/Orthanc-plugins-cpp


text:
// Copyright (C) 2015-2016 Internet Systems Consortium, Inc. ("ISC")
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "config.h"

#include <dhcpsrv/dhcpsrv_log.h>
#include <dhcpsrv/host_data_source_factory.h>
#include <dhcpsrv/hosts_log.h>

#ifdef HAVE_MYSQL
#include <dhcpsrv/mysql_host_data_source.h>
#endif

#ifdef HAVE_PGSQL
#include <dhcpsrv/pgsql_host_data_source.h>
#endif

#include <boost/algorithm/string.hpp>
#include <boost/foreach.hpp>
#include <boost/scoped_ptr.hpp>

#include <algorithm>
#include <iostream>
#include <iterator>
#include <map>
#include <sstream>
#include <utility>

using namespace std;

namespace isc {
namespace dhcp {

HostDataSourcePtr&
HostDataSourceFactory::getHostDataSourcePtr() {
    static HostDataSourcePtr hostDataSourcePtr;
    return (hostDataSourcePtr);
}

void
HostDataSourceFactory::create(const std::string& dbaccess) {
    // Parse the access string and create a redacted string for logging.
    DatabaseConnection::ParameterMap parameters =
            DatabaseConnection::parse(dbaccess);

    // Get the databaase type and open the corresponding database
    DatabaseConnection::ParameterMap::iterator it = parameters.find("type");
    if (it == parameters.end()) {
        isc_throw(InvalidParameter, "Host database configuration does not "
                  "contain the 'type' keyword");
    }

    std::string db_type = it->second;

#ifdef HAVE_MYSQL
    if (db_type == "mysql") {
        LOG_INFO(dhcpsrv_logger, DHCPSRV_MYSQL_HOST_DB)
            .arg(DatabaseConnection::redactedAccessString(parameters));
        getHostDataSourcePtr().reset(new MySqlHostDataSource(parameters));
        return;
    }
#endif

#ifdef HAVE_PGSQL
    if (db_type == "postgresql") {
        LOG_INFO(dhcpsrv_logger, DHCPSRV_PGSQL_HOST_DB)
            .arg(DatabaseConnection::redactedAccessString(parameters));
        getHostDataSourcePtr().reset(new PgSqlHostDataSource(parameters));
        return;
    }
#endif

#ifdef HAVE_CQL
    if (db_type == "cql") {
        isc_throw(NotImplemented, "Sorry, CQL backend for host reservations "
                  "is not implemented yet.");
    }
#endif

    // Get here on no match.
    isc_throw(InvalidType, "Hosts database access parameter 'type': " <<
                           db_type << " is invalid");
}

void
HostDataSourceFactory::destroy() {
    // Destroy current host data source instance.  This is a no-op if no host
    // data source is available.
    if (getHostDataSourcePtr()) {
        LOG_DEBUG(dhcpsrv_logger, DHCPSRV_DBG_TRACE, HOSTS_CFG_CLOSE_HOST_DATA_SOURCE)
            .arg(getHostDataSourcePtr()->getType());
    }
    getHostDataSourcePtr().reset();
}

#if 0
BaseHostDataSource&
HostDataSourceFactory::instance() {
    BaseHostDataSource* hdsptr = getHostDataSourcePtr().get();
    if (hdsptr == NULL) {
        isc_throw(NoHostDataSourceManager,
                "no current host data source instance is available");
    }
    return (*hdsptr);
}
#endif

}; // namespace dhcp
}; // namespace isc

ext: cc

max_stars_repo_name: sebschrader/debian-pkg-isc-kea

max_issues_repo_name: sebschrader/debian-pkg-isc-kea

max_forks_repo_name: sebschrader/debian-pkg-isc-kea


text:
/// Test the Config class
///
///    tconfig  --log_level=all

/*   Part of the rsked package.
 *
 *   Copyright 2020 Steven A. Harp
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 *
 */

/// Dynamically link boost test framework
#define BOOST_TEST_MODULE config_test
#ifndef BOOST_TEST_DYN_LINK
#define BOOST_TEST_DYN_LINK 1
#endif
#include <boost/test/unit_test.hpp>
#include <boost/test/data/test_case.hpp>
#include <boost/test/data/monomorphic.hpp>

#include <iostream>
#include "logging.hpp"
#include "config.hpp"

/// Simple test fixture that just handles logging setup/teardown.
///
struct LogFixture {
    LogFixture() {
        init_logging("tconfig","tconfig_%5N.log",LF_FILE|LF_DEBUG);
    }
    ~LogFixture() {
        finish_logging();
    }
};

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

/// Do some tests on test file. All of these should return normally,
/// none should throw.
///
BOOST_AUTO_TEST_CASE( Valid_params_test )
{
    namespace fs = boost::filesystem;
    LogFixture lf;

    LOG_INFO(Lgr) << "Unit test: Valid_params_test";

    const char *confname = "../test/tconfig.json";
    Config cfg(confname);
    cfg.read_config();

    // Bool
    bool sdr_player_enabled {false};
    BOOST_TEST(cfg.get_bool("Sdr_player","enabled", sdr_player_enabled));
    BOOST_TEST(sdr_player_enabled);

    bool mp3_player_enabled {true};
    BOOST_TEST(cfg.get_bool("Mp3_player","enabled",mp3_player_enabled));
    BOOST_TEST(mp3_player_enabled == false);

    // Int
    int mpd_port=6600;
    BOOST_TEST(cfg.get_int("Mpd_player","mpd_port",mpd_port));
    BOOST_TEST(mpd_port == 6666);

    // Unsigned
    unsigned umpd_port=6600;
    BOOST_TEST(cfg.get_unsigned("Mpd_player","mpd_port",umpd_port));
    BOOST_CHECK_EQUAL(umpd_port,6666);

    unsigned unexpected_port = 1313;
    BOOST_TEST(not cfg.get_unsigned("Sdr_player","gqrx_port",unexpected_port));
    BOOST_CHECK_EQUAL(unexpected_port,1313);

    // Double
    double low_s=13.0;
    BOOST_TEST(cfg.get_double("Sdr_player","gqrx_low_s",low_s));
    BOOST_CHECK_EQUAL(low_s, -20.0);

    // Pathnames
    fs::path path1_expected  { "/usr/bin/whoami" };
    BOOST_TEST( fs::exists(path1_expected) );
    fs::path path1 {"/foo/bar"};
    BOOST_TEST(cfg.get_pathname("Ogg_player","ogg_bin_path",
                                FileCond::MustExist, path1));
    BOOST_CHECK_EQUAL(path1,path1_expected);

    fs::path path2_default("/bin/bash"); // this file should exist
    BOOST_TEST( fs::exists(path2_default) );
    fs::path path2 {path2_default};
    BOOST_TEST(not cfg.get_pathname("NA_player","well_trodden_path",
                                    FileCond::MustExist, path2));
    BOOST_CHECK_EQUAL(path2,path2_default);

    // this path exists but must not; the key is not present in the JSON.
    fs::path validpath("tconfig");
    BOOST_CHECK_THROW(cfg.get_pathname("NA_player","well_trodden_path",
                                       FileCond::MustNotExist, validpath),
                      Config_path_error);

    //  the path specified in the JSON does not exist, should throw
    fs::path invalidpath {"/tmp"};
    BOOST_CHECK_THROW(cfg.get_pathname("NA_player","secret_path",
                                       FileCond::MustExist, invalidpath),
                      Config_path_error );
}

/// Parse the player_preference nested object in tconfig.json
///
BOOST_AUTO_TEST_CASE( Json_parse_test )
{
    LogFixture lf;
    LOG_INFO(Lgr) << "Unit test: Json_parse_test";
    const char *confname = "../test/tconfig.json";
    Config cfg(confname);
    cfg.read_config();

    Json::Value jppref = cfg.get_root()["player_preference"];
    BOOST_TEST(not jppref.isNull());
    BOOST_TEST( jppref.isObject() );
    for (auto medname : jppref.getMemberNames()) {
        LOG_INFO(Lgr) << medname << " :"; // * validate medium here
        auto jmed = jppref[medname];
        BOOST_TEST( jmed.isObject() );
        for (auto encname : jmed.getMemberNames()) {
            LOG_INFO(Lgr) << "    " << encname << ":"; // * validate enc here
            auto jenc = jmed[encname];
            BOOST_TEST( jenc.isArray() );
            unsigned i=0;
            for (auto j : jenc) {
                LOG_INFO(Lgr) << "       (" << i++ << ") " << j.asString();
            }
        }
    }
}

max_stars_repo_name: farlies/rsked

max_issues_repo_name: farlies/rsked

max_forks_repo_name: farlies/rsked