text:
#include <Access/AccessEntityIO.h>
#include <Access/IAccessEntity.h>
#include <Access/IAccessStorage.h>
#include <Access/Quota.h>
#include <Access/Role.h>
#include <Access/RowPolicy.h>
#include <Access/SettingsProfile.h>
#include <Access/User.h>
#include <Core/Defines.h>
#include <IO/WriteHelpers.h>
#include <Interpreters/Access/InterpreterCreateQuotaQuery.h>
#include <Interpreters/Access/InterpreterCreateRoleQuery.h>
#include <Interpreters/Access/InterpreterCreateRowPolicyQuery.h>
#include <Interpreters/Access/InterpreterCreateSettingsProfileQuery.h>
#include <Interpreters/Access/InterpreterCreateUserQuery.h>
#include <Interpreters/Access/InterpreterGrantQuery.h>
#include <Interpreters/Access/InterpreterShowCreateAccessEntityQuery.h>
#include <Interpreters/Access/InterpreterShowGrantsQuery.h>
#include <Parsers/Access/ASTCreateQuotaQuery.h>
#include <Parsers/Access/ASTCreateRoleQuery.h>
#include <Parsers/Access/ASTCreateRowPolicyQuery.h>
#include <Parsers/Access/ASTCreateSettingsProfileQuery.h>
#include <Parsers/Access/ASTCreateUserQuery.h>
#include <Parsers/Access/ASTGrantQuery.h>
#include <Parsers/Access/ParserCreateQuotaQuery.h>
#include <Parsers/Access/ParserCreateRoleQuery.h>
#include <Parsers/Access/ParserCreateRowPolicyQuery.h>
#include <Parsers/Access/ParserCreateSettingsProfileQuery.h>
#include <Parsers/Access/ParserCreateUserQuery.h>
#include <Parsers/Access/ParserGrantQuery.h>
#include <Parsers/ASTSetQuery.h>
#include <Parsers/ParserSetQuery.h>
#include <Parsers/formatAST.h>
#include <Parsers/parseQuery.h>
#include <boost/range/algorithm/copy.hpp>
#include <boost/range/algorithm_ext/push_back.hpp>

namespace DB
{
namespace ErrorCodes
{
    extern const int INCORRECT_ACCESS_ENTITY_DEFINITION;
    extern const int RBAC_VERSION_IS_TOO_NEW;
}

extern const UInt64 RBAC_INITIAL_VERSION;
extern const UInt64 RBAC_LATEST_VERSION;

namespace
{
    constexpr const char RBAC_VERSION_SETTING_NAME[] = "rbac_version";

    /// Special parser for the 'ATTACH access entity' queries.
    class ParserAttachAccessEntity : public IParserBase
    {
    protected:
        const char * getName() const override { return "ATTACH access entity query"; }

        bool parseImpl(Pos & pos, ASTPtr & node, Expected & expected) override
        {
            ParserCreateUserQuery create_user_p;
            ParserCreateRoleQuery create_role_p;
            ParserCreateRowPolicyQuery create_policy_p;
            ParserCreateQuotaQuery create_quota_p;
            ParserCreateSettingsProfileQuery create_profile_p;
            ParserGrantQuery grant_p;
            ParserSetQuery set_p;

            create_user_p.useAttachMode();
            create_role_p.useAttachMode();
            create_policy_p.useAttachMode();
            create_quota_p.useAttachMode();
            create_profile_p.useAttachMode();
            grant_p.useAttachMode();

            return create_user_p.parse(pos, node, expected) || create_role_p.parse(pos, node, expected)
                || create_policy_p.parse(pos, node, expected) || create_quota_p.parse(pos, node, expected)
                || create_profile_p.parse(pos, node, expected) || grant_p.parse(pos, node, expected)
                || set_p.parse(pos, node, expected);
        }
    };

}

String serializeAccessEntity(const IAccessEntity & entity)
{
    ASTs queries;
    {
        /// Prepend the list with "SET rbac_version = ..." query.
        auto set_rbac_version_query = std::make_shared<ASTSetQuery>();
        set_rbac_version_query->changes.emplace_back(RBAC_VERSION_SETTING_NAME, RBAC_LATEST_VERSION);
        queries.push_back(set_rbac_version_query);
    }

    /// Build list of ATTACH queries.
    queries.push_back(InterpreterShowCreateAccessEntityQuery::getAttachQuery(entity));
    if ((entity.getType() == AccessEntityType::USER) || (entity.getType() == AccessEntityType::ROLE))
        boost::range::push_back(queries, InterpreterShowGrantsQuery::getAttachGrantQueries(entity));

    /// Serialize the list of ATTACH queries to a string.
    WriteBufferFromOwnString buf;
    for (const ASTPtr & query : queries)
    {
        formatAST(*query, buf, false, true);
        buf.write(";\n", 2);
    }
    return buf.str();
}

AccessEntityPtr deserializeAccessEntityImpl(const String & definition)
{
    ASTs queries;
    ParserAttachAccessEntity parser;
    const char * begin = definition.data(); /// begin of current query
    const char * pos = begin; /// parser moves pos from begin to the end of current query
    const char * end = begin + definition.size();
    while (pos < end)
    {
        queries.emplace_back(parseQueryAndMovePosition(parser, pos, end, "", true, 0, DBMS_DEFAULT_MAX_PARSER_DEPTH));
        while (isWhitespaceASCII(*pos) || *pos == ';')
            ++pos;
    }

    /// Number of queries interpreted.
    size_t query_index = 0;

    /// If there is no "SET rbac_version = ..." query we assume that it's the initial version.
    UInt64 rbac_version = RBAC_INITIAL_VERSION;

    /// Interpret the AST to build an access entity.
    std::shared_ptr<User> user;
    std::shared_ptr<Role> role;
    std::shared_ptr<RowPolicy> policy;
    std::shared_ptr<Quota> quota;
    std::shared_ptr<SettingsProfile> profile;
    AccessEntityPtr res;

    for (const auto & query : queries)
    {
        if (auto * set_query = query->as<ASTSetQuery>())
        {
            if ((set_query->changes.size() != 1) || (set_query->changes[0].name != RBAC_VERSION_SETTING_NAME))
                throw Exception(ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION, "SET query in this file is only allowed to set {}", RBAC_VERSION_SETTING_NAME);
            if (query_index != 0)
                throw Exception(ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION, "SET {} should be the first query in the file", RBAC_VERSION_SETTING_NAME);
            rbac_version = set_query->changes[0].value.safeGet<UInt64>();
            if (rbac_version < RBAC_INITIAL_VERSION)
                throw Exception(ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION, "{} must be >= {}", RBAC_VERSION_SETTING_NAME, RBAC_INITIAL_VERSION);
            if (rbac_version > RBAC_LATEST_VERSION)
                throw Exception(ErrorCodes::RBAC_VERSION_IS_TOO_NEW, "{} must be <= {}, {} {} is too new", RBAC_VERSION_SETTING_NAME, RBAC_LATEST_VERSION, RBAC_VERSION_SETTING_NAME, rbac_version);
        }
        else if (auto * create_user_query = query->as<ASTCreateUserQuery>())
        {
            if (res)
                throw Exception("Two access entities attached in the same file", ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION);
            res = user = std::make_unique<User>();
            InterpreterCreateUserQuery::updateUserFromQuery(*user, *create_user_query);
        }
        else if (auto * create_role_query = query->as<ASTCreateRoleQuery>())
        {
            if (res)
                throw Exception("Two access entities attached in the same file", ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION);
            res = role = std::make_unique<Role>();
            InterpreterCreateRoleQuery::updateRoleFromQuery(*role, *create_role_query);
        }
        else if (auto * create_policy_query = query->as<ASTCreateRowPolicyQuery>())
        {
            if (res)
                throw Exception("Two access entities attached in the same file", ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION);
            res = policy = std::make_unique<RowPolicy>();
            InterpreterCreateRowPolicyQuery::updateRowPolicyFromQuery(*policy, *create_policy_query, rbac_version);
        }
        else if (auto * create_quota_query = query->as<ASTCreateQuotaQuery>())
        {
            if (res)
                throw Exception("Two access entities attached in the same file", ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION);
            res = quota = std::make_unique<Quota>();
            InterpreterCreateQuotaQuery::updateQuotaFromQuery(*quota, *create_quota_query);
        }
        else if (auto * create_profile_query = query->as<ASTCreateSettingsProfileQuery>())
        {
            if (res)
                throw Exception("Two access entities attached in the same file", ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION);
            res = profile = std::make_unique<SettingsProfile>();
            InterpreterCreateSettingsProfileQuery::updateSettingsProfileFromQuery(*profile, *create_profile_query);
        }
        else if (auto * grant_query = query->as<ASTGrantQuery>())
        {
            if (!user && !role)
                throw Exception("A user or role should be attached before grant", ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION);
            if (user)
                InterpreterGrantQuery::updateUserFromQuery(*user, *grant_query);
            else
                InterpreterGrantQuery::updateRoleFromQuery(*role, *grant_query);
        }
        else
            throw Exception("No interpreter found for query " + query->getID(), ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION);

        ++query_index;
    }

    if (!res)
        throw Exception("No access entities attached", ErrorCodes::INCORRECT_ACCESS_ENTITY_DEFINITION);

    return res;
}

AccessEntityPtr deserializeAccessEntity(const String & definition, const String & file_path)
{
    if (file_path.empty())
        return deserializeAccessEntityImpl(definition);

    try
    {
        return deserializeAccessEntityImpl(definition);
    }
    catch (Exception & e)
    {
        e.addMessage("Could not parse " + file_path);
        e.rethrow();
        __builtin_unreachable();
    }
}

}

ext: cpp

max_stars_repo_name: hanqf-git/ClickHouse

max_issues_repo_name: hanqf-git/ClickHouse

max_forks_repo_name: hanqf-git/ClickHouse


text:
/*
* SVRTK : SVR reconstruction based on MIRTK
*
* Copyright 2008-2017 Imperial College London
* Copyright 2018-2021 King's College London
*
* 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.
*/

// MIRTK
#include "mirtk/Common.h"
#include "mirtk/Options.h"
#include "mirtk/NumericsConfig.h"
#include "mirtk/IOConfig.h"
#include "mirtk/TransformationConfig.h"
#include "mirtk/RegistrationConfig.h"
#include "mirtk/GenericImage.h"
#include "mirtk/GenericRegistrationFilter.h"
#include "mirtk/Transformation.h"
#include "mirtk/HomogeneousTransformation.h"
#include "mirtk/RigidTransformation.h"
#include "mirtk/ImageReader.h"

// SVRTK
#include "svrtk/Reconstruction.h"
#define SVRTK_TOOL
#include "svrtk/Profiling.h"

// Boost
#include <boost/program_options.hpp>
#include <boost/filesystem.hpp>
#include <boost/format.hpp>

// C++ Standard
#include <iostream>
#include <ctime>
#include <fstream>
#include <cmath>
#include <set>
#include <algorithm>
#include <thread>
#include <functional>
#include <vector>
#include <cstdlib>
#include <pthread.h>
#include <string>

using namespace std;
using namespace mirtk;
using namespace boost::program_options;

// =============================================================================
//
// =============================================================================

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

void PrintUsage(const options_description& opts) {
    // Print positional arguments
    cout << "Usage: reconstruct [reconstructed] [N] [stack_1] .. [stack_N] <options>\n" << endl;
    cout << "  [reconstructed]            Name for the reconstructed volume (Nifti or Analyze format)" << endl;
    cout << "  [N]                        Number of stacks" << endl;
    cout << "  [stack_1] .. [stack_N]     The input stacks (Nifti or Analyze format)" << endl << endl;
    // Print optional arguments
    cout << opts << endl;
}

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

// =============================================================================
// Main function
// =============================================================================

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

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

    // -----------------------------------------------------------------------------
    // INPUT VARIABLES, FLAG AND DEFAULT VALUES
    // -----------------------------------------------------------------------------

    // Initialisation of MIRTK image reader library
    InitializeIOLibrary();

    // Initialise profiling
    SVRTK_START_TIMING();

    // General utility variables
    int i, j, x, y, z;

    // Names for output info files
    string infoFilename = "slice_info.tsv";
    string logId;
    string strFlags;

    // Main variables:

    // Name for output volume
    string outputName;

    // Folder for slice-to-volume registrations, if given
    string folder;

    // Array of stacks and stack names
    Array<RealImage> stacks;
    vector<string> stackFiles;

    // Template stack
    RealImage templateStack;

    // Input mask
    unique_ptr<RealImage> mask;

    // Array of stack transformation to the template space
    Array<RigidTransformation> stackTransformations;

    // Default RigidTransformation variable for operations
    const unique_ptr<RigidTransformation> rigidTransformation(new RigidTransformation);

    // Array of stack slice thickness
    vector<double> thickness;

    // Number of input stacks
    int nStacks = 0;

    // Array of number of packages for each stack
    vector<int> packages;

    // Masked template stack
    RealImage maskedTemplate;

    // Output of stack average
    RealImage average;

    // Variables for forced exclusion of slices
    vector<int> forceExcluded;

    // Default values for reconstruction variables:

    int templateNumber = 0;
    int iterations = 3;
    int srIterations = 7;
    bool debug = false;
    double sigma = 20;
    double resolution = 0.75;
    double lambda = 0.02;
    double delta = 150;
    int levels = 3;
    double lastIterLambda = 0.01;
    int recIterations;
    double averageValue = 700;
    double smoothMask = 4;
    bool globalBiasCorrection = false;
    double lowIntensityCutoff = 0.01;

    // Flags for reconstruction options:

    // Flag whether the template (-template option)
    bool useTemplate = false;

    // Flag if the exact should be used
    bool flagNoOverlapThickness = false;

    // Flag to remove black background, e.g. when neonatal motion correction is performed
    bool removeBlackBackground = false;

    // Flag to switch the intensity matching on and off
    bool intensityMatching = true;

    bool rescaleStacks = false;

    // Flag whether registration should be used
    bool registrationFlag = true;

    // Flags to switch the robust statistics on and off
    bool robustStatistics = true;
    bool robustSlicesOnly = false;

    // No log flag
    bool noLog = false;

    // Flag for SVR registration to the template (skips 1st averaging SR iteration)
    bool svrOnly = false;

    // Flag for struture-based exclusion of slices
    bool structural = false;

    // Flag for switching off NMI and using NCC for SVR step
    bool nccRegFlag = false;

    // Flag for running registration step outside
    bool remoteFlag = false;

    // Flag for no global registration
    bool noGlobalFlag = false;

    // Flag that sets slice thickness to 1.5 of spacing (for testing purposes)
    bool thinFlag = false;

    // Flag for longer reconstruction
    bool fullRemoteRecon = false;

    // Flag for automated exclusion of low quality / similarity stacks
    bool excludeWrongStacks = false;
    int bestSelectedStacks = -1;

    // Paths of 'dofin' arguments
    vector<string> dofinPaths;

    // Create reconstruction object
    unique_ptr<Reconstruction> reconstruction(new Reconstruction());

    cout << "------------------------------------------------------" << endl;

    // -----------------------------------------------------------------------------
    // READ INPUT DATA AND OPTIONS
    // -----------------------------------------------------------------------------

    // Define required options
    options_description reqOpts;
    reqOpts.add_options()
        ("reconstructed", value<string>(&outputName)->required(), "Name for the reconstructed volume. Nifti or Analyze format.")
        ("N", value<int>(&nStacks)->required(), "Number of stacks")
        ("stack", value<vector<string>>(&stackFiles)->multitoken()->required(), "The input stacks. Nifti or Analyze format.");

    // Define positional options
    positional_options_description posOpts;
    posOpts.add("reconstructed", 1).add("N", 1).add("stack", -1);

    // Define optional options
    options_description opts("Options");
    opts.add_options()
        ("dofin", value<vector<string>>(&dofinPaths)->multitoken(), "The transformations of the input stack to template in \'dof\' format used in IRTK. Only rough alignment with correct orienation and some overlap is needed. Use \'id\' for an identity transformation for at leastone stack. The first stack with \'id\' transformationwill be resampled as template.")
        ("template", value<string>(), "Use template for initialisation of registration loop. [Default: average of stack registration]")
        ("thickness", value<vector<double>>(&thickness)->multitoken(), "Give slice thickness. [Default: twice voxel size in z direction]")
        ("default_thickness", value<double>(), "Default thickness for all stacks. [Default: twice voxel size in z direction]")
        ("mask", value<string>(), "Binary mask to define the region od interest. [Default: whole image]")
        ("packages", value<vector<int>>(&packages)->multitoken(), "Give number of packages used during acquisition for each stack. The stacks will be split into packages during registration iteration 1 and then into odd and even slices within each package during registration iteration 2. The method will then continue with slice to volume approach. [Default: slice to volume registration only]")
        ("template_number", value<int>(&templateNumber), "Number of the template stack [Default: 0]")
        ("iterations", value<int>(&iterations), "Number of registration-reconstruction iterations [Default: 3]")
        ("sr_iterations", value<int>(&srIterations), "Number of SR reconstruction iterations [Default: 7,...,7,7*3]")
        ("sigma", value<double>(&sigma), "Stdev for bias field [Default: 20mm]")
        ("resolution", value<double>(&resolution), "Isotropic resolution of the volume [Default: 0.75mm]")
        ("multires", value<int>(&levels), "Multiresolution smooting with given number of levels [Default: 3]")
        ("average", value<double>(&averageValue), "Average intensity value for stacks [Default: 700]")
        ("delta", value<double>(&delta), "Parameter to define what is an edge [Default: 150]")
        ("lambda", value<double>(&lambda), "Smoothing parameter [Default: 0.02]")
        ("lastIter", value<double>(&lastIterLambda), "Smoothing parameter for last iteration [Default: 0.01]")
        ("smooth_mask", value<double>(&smoothMask), "Smooth the mask to reduce artefacts of manual segmentation [Default: 4mm]")
        ("global_bias_correction", bool_switch(&globalBiasCorrection), "Correct the bias in reconstructed image against previous estimation")
        ("no_intensity_matching", "Switch off intensity matching")
        ("no_robust_statistics", "Switch off robust statistics")
        ("exclude_wrong_stacks", bool_switch(&excludeWrongStacks), "Automated exclusion of misregistered stacks")
        ("rescale_stacks", bool_switch(&rescaleStacks), "Rescale stacks to avoid nan pixel errors [Default: false]")
        ("svr_only", bool_switch(&svrOnly), "Only SVR registration to a template stack")
        ("no_global", bool_switch(&noGlobalFlag), "No global stack registration")
        ("exact_thickness", bool_switch(&flagNoOverlapThickness), "Exact slice thickness without negative gap [Default: false]")
        ("ncc", bool_switch(&nccRegFlag), "Use global NCC similarity for SVR steps [Default: NMI]")
        ("nmi_bins", value<int>(), "Number of NMI bins for registration [Default: 16]")
        ("structural", bool_switch(&structural), "Use structrural exclusion of slices at the last iteration")
        ("exclude_slices_only", bool_switch(&robustSlicesOnly), "Robust statistics for exclusion of slices only")
        ("remove_black_background", bool_switch(&removeBlackBackground), "Create mask from black background")
        ("transformations", value<string>(&folder), "Use existing slice-to-volume transformations to initialize the reconstruction")
        ("force_exclude", value<vector<int>>(&forceExcluded), "Force exclusion of slices with these indices")
        ("remote", bool_switch(&remoteFlag), "Run SVR registration as remote functions in case of memory issues [Default: false]")
        ("full_remote", bool_switch(&fullRemoteRecon), "Run SR+SVR steps as remote functions in case of memory issues (slower option) [Default: false]")
        ("no_registration", "Switch off registration")
        ("thin", bool_switch(&thinFlag), "Option for 1.5 x dz slice thickness (testing)")
        ("debug", bool_switch(&debug), "Debug mode - save intermediate results");

    // Combine all options
    options_description allOpts("Allowed options");
    allOpts.add(reqOpts).add(opts);

    // Parse arguments and catch errors
    variables_map vm;
    try {
        store(command_line_parser(argc, argv).options(allOpts).positional(posOpts)
            // Allow single dash (-) for long arguments
            .style(command_line_style::unix_style | command_line_style::allow_long_disguise).run(), vm);
        notify(vm);

        if (stackFiles.size() < nStacks)
            throw error("Count of input stacks should equal to stack count!");
        if (!dofinPaths.empty() && dofinPaths.size() < nStacks)
            throw error("Count of dof files should equal to stack count!");
        if (!thickness.empty() && thickness.size() < nStacks)
            throw error("Count of thickness values should equal to stack count!");
        if (!packages.empty() && packages.size() < nStacks)
            throw error("Count of package values should equal to stack count!");
    } catch (error& e) {
        // Delete -- from the argument name in the error message
        string err = e.what();
        size_t dashIndex = err.find("\'--");
        if (dashIndex != string::npos)
            err.erase(dashIndex + 1, 2);
        cerr << "Argument parsing error: " << err << "\n\n";
        PrintUsage(opts);
        return 1;
    }

    cout << "Reconstructed volume name : " << outputName << endl;
    cout << "Number of stacks : " << nStacks << endl;

    // Read input stacks
    for (i = 0; i < nStacks; i++) {
        cout << "Stack " << i << " : " << stackFiles[i];
        unique_ptr<RealImage> stack(new RealImage(stackFiles[i].c_str()));

        // Check if the intensity is not negative and correct if so
        double smin, smax;
        stack->GetMinMax(&smin, &smax);
        if (smin < 0 || smax < 0)
            stack->PutMinMaxAsDouble(0, 1000);

        // Print stack info
        double dx = stack->GetXSize(); double dy = stack->GetYSize();
        double dz = stack->GetZSize(); double dt = stack->GetTSize();
        double sx = stack->GetX(); double sy = stack->GetY();
        double sz = stack->GetZ(); double st = stack->GetT();

        cout << "  ;  size : " << sx << " - " << sy << " - " << sz << " - " << st << "  ;";
        cout << "  voxel : " << dx << " - " << dy << " - " << dz << " - " << dt << "  ;";
        cout << "  range : [" << smin << "; " << smax << "]" << endl;

        stacks.push_back(*stack);
    }

    // Default template stack
    templateStack = stacks[0];

    // Input stack transformations to the template space
    if (!dofinPaths.empty()) {
        for (i = 0; i < nStacks; i++) {
            cout << "Transformation " << i << " : " << dofinPaths[i] << endl;
            Transformation *t = Transformation::New(dofinPaths[i].c_str());
            RigidTransformation *rigidTransf = dynamic_cast<RigidTransformation*>(t);
            stackTransformations.push_back(*rigidTransf);
            delete rigidTransf;
        }
        reconstruction->InvertStackTransformations(stackTransformations);
    }

    // Slice thickness per stack
    if (!thickness.empty()) {
        cout << "Slice thickness : ";
        for (i = 0; i < nStacks; i++)
            cout << thickness[i] << " ";
        cout << endl;
    }

    // Number of packages for each stack
    if (!packages.empty()) {
        cout << "Package number : ";
        for (i = 0; i < nStacks; i++)
            cout << packages[i] << " ";
        cout << endl;
        reconstruction->SetNPackages(packages);
    }

    // Template for initialisation of registration
    if (vm.count("template")) {
        cout << "Template : " << vm["template"].as<string>() << endl;
        templateStack.Read(vm["template"].as<string>().c_str());

        // Check intensities and rescale if required
        double smin, smax;
        templateStack.GetMinMax(&smin, &smax);
        if (smin < 0 || smax < 0)
            templateStack.PutMinMaxAsDouble(0, 1000);

        // Extract 1st dynamic
        if (templateStack.GetT() > 1)
            templateStack = templateStack.GetRegion(0, 0, 0, 0, templateStack.GetX(), templateStack.GetY(), templateStack.GetZ(), 1);

        useTemplate = true;
        reconstruction->SetTemplateFlag(useTemplate);
    } else if (vm.count("template_number"))
        templateStack = stacks[templateNumber];

    // Binary mask for reconstruction / final volume
    if (vm.count("mask")) {
        cout << "Mask : " << vm["mask"].as<string>() << endl;
        mask = unique_ptr<RealImage>(new RealImage(vm["mask"].as<string>().c_str()));
    }

    // Number of registration-reconstruction iterations
    if (vm.count("iterations"))
        strFlags += " -iterations " + to_string(iterations);

    // The same thickness for all stacks
    if (vm.count("default_thickness")) {
        double defaultThickness = vm["default_thickness"].as<double>();
        cout << "Slice thickness (default for all stacks): ";
        for (i = 0; i < nStacks; i++) {
            thickness.push_back(defaultThickness);
            cout << thickness[i] << " ";
        }
        cout << endl;
    }

    // Number of SR iterations
    if (vm.count("sr_iterations"))
        strFlags += " -sr_iterations " + to_string(srIterations);

    // Variance of Gaussian kernel to smooth the bias field
    if (vm.count("sigma"))
        strFlags += " -sigma " + to_string(sigma);

    // SR smoothing parameter
    if (vm.count("lambda"))
        strFlags += " -lambda " + to_string(lambda);

    // Smoothing parameter for last iteration
    if (vm.count("lastIter"))
        strFlags += " -lastIter " + to_string(lastIterLambda);

    // SR parameter to define what is an edge
    if (vm.count("delta"))
        strFlags += " -delta " + to_string(delta);

    // Isotropic resolution for the reconstructed volume
    if (vm.count("resolution"))
        strFlags += " -resolution " + to_string(resolution);

    // Number of resolution levels
    if (vm.count("multires"))
        strFlags += " -multires " + to_string(levels);

    // Run registration as remote functions
    if (remoteFlag)
        strFlags += " -remote";

    // Long reconstruction (10x10x30 iterations)
    if (fullRemoteRecon) {
        remoteFlag = true;
        strFlags += " -remote";
    }

    // Use only SVR to the template (skip 1st SR averaging iteration)
    if (svrOnly)
        strFlags += " -svr_only";

    // Use NCC similarity metric for SVR
    if (nccRegFlag) {
        reconstruction->SetNCC(nccRegFlag);
        strFlags += " -ncc";
    }

    // Force removal of certain slices
    if (!forceExcluded.empty()) {
        cout << forceExcluded.size() << " force excluded slices: ";
        for (i = 0; i < forceExcluded.size(); i++)
            cout << forceExcluded[i] << " ";
        cout << endl;
    }

    // Switch off intensity matching
    if (vm.count("no_intensity_matching")) {
        intensityMatching = false;
        strFlags += " -no_intensity_matching";
    }

    // Switch off robust statistics
    if (vm.count("no_robust_statistics")) {
        robustStatistics = false;
        strFlags += " -no_robust_statistics";
    }

    // Use structural exclusion of slices (after 2nd iteration)
    if (structural)
        strFlags += " -structural";

    // Use robust statistics for slices only
    if (robustSlicesOnly)
        strFlags += " -exclude_slices_only";

    // Switch off registration
    if (vm.count("no_registration")) {
        registrationFlag = false;
        strFlags += " -no_registration";
    }

    // Perform bias correction of the reconstructed image agains the GW image in the same motion correction iteration
    if (globalBiasCorrection)
        strFlags += " -global_bias_correction";

    // Debug mode
    if (debug)
        strFlags += " -debug";

    // -----------------------------------------------------------------------------
    // SET RECONSTRUCTION OPTIONS AND PERFORM PREPROCESSING
    // -----------------------------------------------------------------------------

    // Set thickess to the exact dz value if specified
    if (flagNoOverlapThickness && thickness.size() < 1) {
        for (i = 0; i < stacks.size(); i++)
            thickness.push_back(stacks[i].GetZSize());
    }

    // Check if stacks have multiple dynamics and spit them if it is the case
    bool has4DStacks = false;
    for (i = 0; i < stacks.size(); i++) {
        if (stacks[i].GetT() > 1) {
            has4DStacks = true;
            break;
        }
    }

    if (has4DStacks) {
        cout << "Splitting stacks into dynamics ... ";

        Array<double> newThickness;
        Array<int> newPackages;
        Array<RigidTransformation> newStackTransformations;
        Array<RealImage> newStacks;

        for (i = 0; i < stacks.size(); i++) {
            if (stacks[i].GetT() == 1) {
                newStacks.push_back(stacks[i]);
                if (!stackTransformations.empty())
                    newStackTransformations.push_back(stackTransformations[i]);
                if (!packages.empty())
                    newPackages.push_back(packages[i]);
                if (!thickness.empty())
                    newThickness.push_back(thickness[i]);
            } else {
                for (int t = 0; t < stacks[i].GetT(); t++) {
                    RealImage stack = stacks[i].GetRegion(0, 0, 0, t, stacks[i].GetX(), stacks[i].GetY(), stacks[i].GetZ(), t + 1);
                    newStacks.push_back(move(stack));
                    if (!stackTransformations.empty())
                        newStackTransformations.push_back(stackTransformations[i]);
                    if (!packages.empty())
                        newPackages.push_back(packages[i]);
                    if (!thickness.empty())
                        newThickness.push_back(thickness[i]);
                }
            }
        }

        nStacks = newStacks.size();
        stacks = move(newStacks);
        thickness = move(newThickness);
        packages = move(newPackages);
        stackTransformations = move(newStackTransformations);

        cout << "New number of stacks : " << nStacks << endl;
    }

    // Read path to MIRTK executables for remote registration
    string strMirtkPath(argv[0]);
    strMirtkPath = strMirtkPath.substr(0, strMirtkPath.find_last_of("/"));
    const string strCurrentMainFilePath = boost::filesystem::current_path().string();

    // Create an empty file exchange directory
    const string strCurrentExchangeFilePath = strCurrentMainFilePath + "/tmp-file-exchange";
    boost::filesystem::remove_all(strCurrentExchangeFilePath.c_str());
    boost::filesystem::create_directory(strCurrentExchangeFilePath.c_str());

    // Rescale stack if specified
    if (rescaleStacks) {
        for (i = 0; i < nStacks; i++)
            reconstruction->Rescale(stacks[i], 1000);
    }

    // If transformations were not defined by user, set them to identity
    for (i = 0; i < nStacks; i++)
        stackTransformations.push_back(*rigidTransformation);

    // Initialise 2*slice thickness if not given by user
    if (thickness.size() == 0) {
        cout << "Slice thickness : ";
        for (i = 0; i < nStacks; i++) {
            double dx, dy, dz;
            stacks[i].GetPixelSize(&dx, &dy, &dz);
            thickness.push_back(thinFlag ? dz * 1.5 : dz * 2);
            cout << thickness[i] << " ";
        }
        cout << endl;
    }

    // Set debug mode option
    if (debug) reconstruction->DebugOn();
    else reconstruction->DebugOff();

    // Set verbose mode on with file
    reconstruction->VerboseOn("log-registration.txt");

    // Set force excluded slices option
    reconstruction->SetForceExcludedSlices(forceExcluded);

    // Set low intensity cutoff for bias estimation
    reconstruction->SetLowIntensityCutoff(lowIntensityCutoff);

    // If no mask was given - try to create mask from the template image in case it was padded
    if (mask == NULL) {
        mask = unique_ptr<RealImage>(new RealImage(stacks[templateNumber]));
        *mask = reconstruction->CreateMask(*mask);
        cout << "Warning : no mask was provided " << endl;
    }

    // Before creating the template we will crop template stack according to the given mask
    if (mask != NULL) {
        // First resample the mask to the space of the stack
        // For template stact the transformation is identity
        RealImage m = *mask;
        reconstruction->TransformMask(stacks[templateNumber], m, stackTransformations[templateNumber]);

        // Crop template stack and prepare template for global volumetric registration
        maskedTemplate = stacks[templateNumber] * m;
        reconstruction->CropImage(stacks[templateNumber], m);
        reconstruction->CropImage(maskedTemplate, m);

        if (debug) {
            m.Write("maskforTemplate.nii.gz");
            stacks[templateNumber].Write("croppedTemplate.nii.gz");
        }
    }

    // If the template was provided separately - crop and mask the template with the given mask
    if (useTemplate) {
        RealImage m = *mask;
        reconstruction->TransformMask(templateStack, m, *rigidTransformation);

        // Crop template stack and prepare template for global volumetric registration
        maskedTemplate = templateStack * m;
        reconstruction->CropImage(maskedTemplate, m);
        reconstruction->CropImage(templateStack, m);
    }

    // Create template volume with isotropic resolution
    // If resolution==0 it will be determined from in-plane resolution of the image
    resolution = reconstruction->CreateTemplate(maskedTemplate, resolution);

    // Set mask to reconstruction object
    reconstruction->SetMask(mask.get(), smoothMask);

    // If remove_black_background flag is set, create mask from black background of the stacks
    if (removeBlackBackground)
        reconstruction->CreateMaskFromBlackBackground(stacks, stackTransformations, smoothMask);

    // Set precision
    cout << setprecision(3);
    cerr << setprecision(3);

    // -----------------------------------------------------------------------------
    // RUN GLOBAL STACK REGISTRATION AND FURTHER PREPROCESSING
    // -----------------------------------------------------------------------------

    cout << "------------------------------------------------------" << endl;

    // Volumetric stack to template registration
    if (!noGlobalFlag)
        reconstruction->StackRegistrations(stacks, stackTransformations, templateNumber);

    // Create average volume
    average = reconstruction->CreateAverage(stacks, stackTransformations);
    if (debug)
        average.Write("average1.nii.gz");

    // Mask is transformed to the all other stacks and they are cropped
    for (i = 0; i < nStacks; i++) {
        // Template stack has been cropped already
        if ((i == templateNumber))
            continue;
        // Transform the mask
        RealImage m = reconstruction->GetMask();
        reconstruction->TransformMask(stacks[i], m, stackTransformations[i]);
        // Crop template stack
        reconstruction->CropImage(stacks[i], m);

        if (debug) {
            m.Write((boost::format("mask%1%.nii.gz") % i).str().c_str());
            stacks[i].Write((boost::format("cropped%1%.nii.gz") % i).str().c_str());
        }
    }

    // Remove small stacks (no intersection with ROI)
    Array<RealImage> selectedStacks;
    Array<RigidTransformation> selectedStackTransformations;
    int nNewStacks = 0;
    int newTemplateNumber = 0;
    for (i = 0; i < nStacks; i++) {
        if (stacks[i].GetX() == 1) {
            cerr << "stack " << i << " has no intersection with ROI" << endl;
            continue;
        }
        selectedStacks.push_back(stacks[i]);
        selectedStackTransformations.push_back(stackTransformations[i]);
        if (i == templateNumber)
            newTemplateNumber = nNewStacks;
        nNewStacks++;
    }
    nStacks = nNewStacks;
    templateNumber = newTemplateNumber;
    stacks = move(selectedStacks);
    stackTransformations = move(selectedStackTransformations);

    // Perform volumetric registration again
    if (!noGlobalFlag)
        reconstruction->StackRegistrations(stacks, stackTransformations, templateNumber);

    // Exclude low quality / similarity stacks (should be transferred to a separate function)
    if (excludeWrongStacks) {
        bestSelectedStacks = stacks.size();
        cout << "Selecting stacks : " << " " << endl;

        RealImage transformedTemplateMask = *mask;
        RealImage templateToCheck = templateStack;
        reconstruction->TransformMask(templateToCheck, transformedTemplateMask, *rigidTransformation);
        reconstruction->CropImage(templateToCheck, transformedTemplateMask);
        reconstruction->CropImage(transformedTemplateMask, transformedTemplateMask);

        Array<double> allNccArray;
        Array<double> allSliceNccArray;
        Array<int> allIndicesArray;
        Array<double> allCountArray;
        Array<double> sortedNccArray;
        Array<double> selectedNccArray;
        Array<int> selectedIndicesArray;
        double maxNcc = -1.0;
        int maxIndex = -1;
        int averageCountNcc = 0;
        double averageSliceNcc = 0;
        double averageVolumeNcc = 0;
        Array<RealImage> newStacks;
        Array<double> newThickness;
        Array<int> newPackages;
        Array<RigidTransformation> newStackTransformations;

        for (i = 0; i < stacks.size(); i++) {
            RealImage stackToCheck = stacks[i];
            Matrix m = stackTransformations[i].GetMatrix();
            stackToCheck.PutAffineMatrix(m, true);

            const double sliceNcc = reconstruction->VolumeNCC(stackToCheck, templateToCheck, transformedTemplateMask);
            double countNcc = -1;
            const double volumeNcc = reconstruction->ComputeNCC(stackToCheck, templateToCheck, 0.1, &countNcc);

            averageCountNcc += countNcc;
            averageSliceNcc += sliceNcc;
            averageVolumeNcc += volumeNcc;

            allNccArray.push_back(volumeNcc);
            allIndicesArray.push_back(i);
            allCountArray.push_back(countNcc);
            allSliceNccArray.push_back(sliceNcc);

            if (volumeNcc > maxNcc) {
                maxNcc = volumeNcc;
                maxIndex = i;
            }
        }

        averageCountNcc /= stacks.size();
        averageSliceNcc /= stacks.size();
        averageVolumeNcc /= stacks.size();

        double stdCountNcc = 0;
        double stdSliceNcc = 0;
        double stdVolumeNcc = 0;

        for (i = 0; i < stacks.size(); i++) {
            stdCountNcc += pow(allCountArray[i] - averageCountNcc, 2);
            stdSliceNcc += pow(allSliceNccArray[i] - averageSliceNcc, 2);
            stdVolumeNcc = averageVolumeNcc + pow(allNccArray[i] - averageVolumeNcc, 2);
        }

        stdSliceNcc /= stacks.size();
        stdVolumeNcc /= stacks.size();
        stdCountNcc /= stacks.size();

        cout << " - average values : volume ncc = " << averageVolumeNcc << " +/- " << stdVolumeNcc << " ; slice ncc = " << averageSliceNcc << " +/- " << stdSliceNcc << " ; volume [mm^3] = " << averageCountNcc << "+/- " << stdCountNcc << endl;

        sortedNccArray = allNccArray;
        sort(sortedNccArray.begin(), sortedNccArray.end(), greater<double>());

        cout << " - selected : " << endl;
        int totalSelected = 0;
        for (int j = 0; j < sortedNccArray.size(); j++) {
            for (i = 0; i < stacks.size(); i++) {
                if (totalSelected < bestSelectedStacks) {
                    if (sortedNccArray[j] == allNccArray[i] && allCountArray[i] > 0.65 * averageCountNcc && allNccArray[i] > (averageVolumeNcc - 2 * stdVolumeNcc)) {
                        selectedNccArray.push_back(allNccArray[i]);
                        selectedIndicesArray.push_back(allIndicesArray[i]);
                        totalSelected++;
                        newStacks.push_back(stacks[allIndicesArray[i]]);
                        cout << "" << allIndicesArray[i] << " : volume ncc = " << allNccArray[i] << " ; slice ncc = " << allSliceNccArray[i] << endl;

                        if (!stackTransformations.empty())
                            newStackTransformations.push_back(stackTransformations[i]);
                        if (!packages.empty())
                            newPackages.push_back(packages[i]);
                        if (!thickness.empty())
                            newThickness.push_back(thickness[i]);
                    }
                }
            }
        }

        nStacks = newStacks.size();
        stacks = move(newStacks);
        thickness = move(newThickness);
        packages = move(newPackages);
        stackTransformations = move(newStackTransformations);
    }

    cout << "------------------------------------------------------" << endl;

    // Rescale intensities of the stacks to have the same average
    if (intensityMatching)
        reconstruction->MatchStackIntensitiesWithMasking(stacks, stackTransformations, averageValue);

    // Create average of the registered stacks
    average = reconstruction->CreateAverage(stacks, stackTransformations);
    if (debug)
        average.Write("average2.nii.gz");

    // Create slices and slice-dependent transformations
    Array<RealImage> probabilityMaps;
    reconstruction->CreateSlicesAndTransformations(stacks, stackTransformations, thickness, probabilityMaps);

    // If full remote reconstruction option is used - save the model and all outputs
    if (fullRemoteRecon)
        reconstruction->SaveModelRemote(strCurrentExchangeFilePath, 1, 0);

    // Mask all the slices
    reconstruction->MaskSlices();

    // Set sigma for the bias field smoothing
    if (sigma > 0)
        reconstruction->SetSigma(sigma);
    else
        reconstruction->SetSigma(20);

    // Set global bias correction flag
    if (globalBiasCorrection)
        reconstruction->GlobalBiasCorrectionOn();
    else
        reconstruction->GlobalBiasCorrectionOff();

    // If given read slice-to-volume registrations
    if (!folder.empty())
        reconstruction->ReadTransformation((char*)folder.c_str());

    // Initialise data structures for EM
    reconstruction->InitializeEM();

    // If registration was switched off - only 1 iteration is required
    if (!registrationFlag)
        iterations = 1;

    // -----------------------------------------------------------------------------
    // RUN INTERLEAVED SVR-SR RECONSTRUCTION
    // -----------------------------------------------------------------------------

    int currentIteration = 0;

    if (fullRemoteRecon && packages.size() < 1) {
        // Run reconstruction remotely ("reconstruct-round" function) - can be removed
        for (int iter = 0; iter < iterations; iter++) {
            currentIteration = iter;
            const string reconstructCmd = strMirtkPath + "/reconstruct-round " + " " + strMirtkPath + " " + strCurrentMainFilePath + " " + strCurrentExchangeFilePath + " " + to_string(currentIteration) + " " + to_string(reconstruction->_number_of_slices_org) + " " + to_string(reconstruction->_average_thickness_org) + " " + strFlags;
            if (system(reconstructCmd.c_str()) == -1) {
                cerr << "The reconstruct command couldn't be executed!" << endl;
                return 1;
            }
        }
    } else {
        // Interleaved registration-reconstruction iterations
        for (int iter = 0; iter < iterations; iter++) {
            cout << "------------------------------------------------------" << endl;
            cout << "Iteration : " << iter << endl;

            reconstruction->MaskVolume();

            // If only SVR option is used - skip 1st SR only averaging
            if (svrOnly) {
                if (remoteFlag)
                    reconstruction->RemoteSliceToVolumeRegistration(iter, strMirtkPath, strCurrentMainFilePath, strCurrentExchangeFilePath);
                else
                    reconstruction->SliceToVolumeRegistration();
            } else if (iter > 0) {
                // Run package-based registartion is the number of packages was given
                if (!packages.empty() && iter < iterations - 1) {
                    reconstruction->PackageToVolume(stacks, packages, stackTransformations);
                } else {
                    // Run
                    if (remoteFlag)
                        reconstruction->RemoteSliceToVolumeRegistration(iter, strMirtkPath, strCurrentMainFilePath, strCurrentExchangeFilePath);
                    else
                        reconstruction->SliceToVolumeRegistration();
                }
            }

            // Run structure-based outlier rejection if specified
            if (structural && iter > 1)
                reconstruction->StructuralExclusion();

            // Set smoothing parameters
            // Amount of smoothing (given by lambda) is decreased with improving alignment
            // Delta (to determine edges) stays constant throughout
            if (iter == (iterations - 1))
                reconstruction->SetSmoothingParameters(delta, lastIterLambda);
            else {
                double l = lambda;
                for (i = 0; i < levels; i++) {
                    if (iter == iterations * (levels - i - 1) / levels)
                        reconstruction->SetSmoothingParameters(delta, l);
                    l *= 2;
                }
            }

            // Use faster reconstruction during iterations and slower for final reconstruction
            if (iter < iterations - 1)
                reconstruction->SpeedupOn();
            else
                reconstruction->SpeedupOff();

            if (robustSlicesOnly)
                reconstruction->ExcludeWholeSlicesOnly();

            reconstruction->InitializeEMValues();

            // Calculate matrix of transformation between voxels of slices and volume
            reconstruction->CoeffInit();

            // Initialise reconstructed image with Gaussian weighted reconstruction
            reconstruction->GaussianReconstruction();

            // Simulate slices (needs to be done after Gaussian reconstruction)
            reconstruction->SimulateSlices();

            // Initialize robust statistics parameters
            reconstruction->InitializeRobustStatistics();

            // EStep
            if (robustStatistics)
                reconstruction->EStep();

            // Set number of reconstruction iterations
            if (iter == iterations - 1)
                recIterations = srIterations * 3;
            else
                recIterations = srIterations;

            // SR reconstruction loop
            for (i = 0; i < recIterations; i++) {
                if (debug) {
                    cout << "------------------------------------------------------" << endl;
                    cout << "Reconstruction iteration : " << i << endl;
                }

                if (intensityMatching) {
                    // Calculate bias fields
                    if (sigma > 0)
                        reconstruction->Bias();

                    // Calculate scales
                    reconstruction->Scale();
                }

                // Update reconstructed volume - super-resolution reconstruction
                reconstruction->Superresolution(i + 1);

                // Run bias normalisation
                if (intensityMatching && sigma > 0 && !globalBiasCorrection)
                    reconstruction->NormaliseBias(i);

                // Simulate slices (needs to be done after the update of the reconstructed volume)
                reconstruction->SimulateSlices();

                // Run robust statistics for rejection of outliers
                if (robustStatistics) {
                    SVRTK_START_TIMING();
                    reconstruction->MStep(i + 1);
                    reconstruction->EStep();
                    SVRTK_END_TIMING("robust statistics");
                }

                if (debug) {
                    // Save intermediate reconstructed image
                    reconstruction->GetReconstructed().Write((boost::format("super%1%.nii.gz") % i).str().c_str());

                    // Evaluate reconstruction quality
                    double error = reconstruction->EvaluateReconQuality(1);
                    cout << "Total reconstruction error : " << error << endl;
                }

            } // End of SR reconstruction iterations

            // Mask reconstructed image to ROI given by the mask
            reconstruction->MaskVolume();

            // Save reconstructed image
            reconstruction->GetReconstructed().Write((boost::format("image%1%.nii.gz") % iter).str().c_str());

            // Compute and save quality metrics
            double outNcc = 0;
            double outNrmse = 0;
            double averageVolumeWeight = 0;
            double ratioExcluded = 0;
            reconstruction->ReconQualityReport(outNcc, outNrmse, averageVolumeWeight, ratioExcluded);
            cout << " - global metrics: ncc = " << outNcc << " ; nrmse = " << outNrmse << " ; average weight = " << averageVolumeWeight << " ; excluded slices = " << ratioExcluded << endl;

            ofstream ofsNcc("output-metric-ncc.txt");
            ofstream ofsNrmse("output-metric-nrmse.txt");
            ofstream ofsWeight("output-metric-average-weight.txt");
            ofstream ofsExcluded("output-metric-excluded-ratio.txt");

            ofsNcc << outNcc << endl;
            ofsNrmse << outNrmse << endl;
            ofsWeight << averageVolumeWeight << endl;
            ofsExcluded << ratioExcluded << endl;

        } // End of interleaved registration-reconstruction iterations

        // -----------------------------------------------------------------------------
        // SAVE RESULTS
        // -----------------------------------------------------------------------------

        cout << "------------------------------------------------------" << endl;

        reconstruction->RestoreSliceIntensities();

        if (debug) {
            reconstruction->SaveTransformations();
            reconstruction->SaveSlices();
            reconstruction->SaveWeights();
            reconstruction->SaveBiasFields();
            reconstruction->SimulateStacks(stacks);
            for (unsigned int i = 0; i < stacks.size(); i++)
                stacks[i].Write((boost::format("simulated%1%.nii.gz") % i).str().c_str());
        }

        reconstruction->ScaleVolume();
    }

    // Load results if remote reconstruction was used - this option can be removed
    if (fullRemoteRecon) {
        reconstruction->LoadResultsRemote(strCurrentExchangeFilePath, reconstruction->_number_of_slices_org, currentIteration);
        reconstruction->ScaleVolume();
    }

    // Remove the file exchange directory
    boost::filesystem::remove_all(strCurrentExchangeFilePath.c_str());

    reconstruction->GetReconstructed().Write(outputName.c_str());

    cout << "Output volume : " << outputName << endl;

    SVRTK_END_TIMING("all");

    cout << "------------------------------------------------------" << endl;

    return 0;
}

ext: cc

max_stars_repo_name: SVRTK/SVRTK

max_issues_repo_name: SVRTK/SVRTK

max_forks_repo_name: SVRTK/SVRTK


text:
              //
// Copyright (c) 2016-2018 CNRS
//

#include "pinocchio/multibody/model.hpp"
#include "pinocchio/multibody/data.hpp"
#include "pinocchio/algorithm/jacobian.hpp"
#include "pinocchio/algorithm/frames.hpp"
#include "pinocchio/algorithm/rnea.hpp"
#include "pinocchio/spatial/act-on-set.hpp"
#include "pinocchio/parsers/sample-models.hpp"
#include "pinocchio/utils/timer.hpp"
#include "pinocchio/algorithm/joint-configuration.hpp"

#include <iostream>

#include <boost/test/unit_test.hpp>
#include <boost/utility/binary.hpp>

template<typename Derived>
inline bool isFinite(const Eigen::MatrixBase<Derived> & x)
{
  return ((x - x).array() == (x - x).array()).all();
}

BOOST_AUTO_TEST_SUITE ( BOOST_TEST_MODULE )

BOOST_AUTO_TEST_CASE(frame_basic)
{
  using namespace pinocchio;
  Model model;
  buildModels::humanoidRandom(model);

  BOOST_CHECK(model.frames.size() >= size_t(model.njoints));
  for(Model::FrameVector::const_iterator it = model.frames.begin();
      it != model.frames.end(); ++it)
  {
    const Frame & frame = *it;
    BOOST_CHECK(frame == frame);
    Frame frame_copy(frame);
    BOOST_CHECK(frame_copy == frame);
  }

  std::ostringstream os;
  os << Frame() << std::endl;
  BOOST_CHECK(!os.str().empty());
}

BOOST_AUTO_TEST_CASE(cast)
{
  using namespace pinocchio;
  Frame frame("toto",0,0,SE3::Random(),OP_FRAME);

  BOOST_CHECK(frame.cast<double>() == frame);
  BOOST_CHECK(frame.cast<long double>().cast<double>() == frame);
}

BOOST_AUTO_TEST_CASE ( test_kinematics )
{
  using namespace Eigen;
  using namespace pinocchio;

  pinocchio::Model model;
  pinocchio::buildModels::humanoidRandom(model);
  Model::Index parent_idx = model.existJointName("rarm2_joint")?model.getJointId("rarm2_joint"):(Model::Index)(model.njoints-1);
  const std::string & frame_name = std::string( model.names[parent_idx]+ "_frame");
  const SE3 & framePlacement = SE3::Random();
  model.addFrame(Frame (frame_name, parent_idx, 0, framePlacement, OP_FRAME));
  pinocchio::Data data(model);

  VectorXd q = VectorXd::Ones(model.nq);
  q.middleRows<4> (3).normalize();
  framesForwardKinematics(model, data, q);

  BOOST_CHECK(data.oMf[model.getFrameId(frame_name)].isApprox(data.oMi[parent_idx]*framePlacement));

}

BOOST_AUTO_TEST_CASE ( test_update_placements )
{
  using namespace Eigen;
  using namespace pinocchio;

  pinocchio::Model model;
  pinocchio::buildModels::humanoidRandom(model);
  Model::Index parent_idx = model.existJointName("rarm2_joint")?model.getJointId("rarm2_joint"):(Model::Index)(model.njoints-1);
  const std::string & frame_name = std::string( model.names[parent_idx]+ "_frame");
  const SE3 & framePlacement = SE3::Random();
  model.addFrame(Frame (frame_name, parent_idx, 0, framePlacement, OP_FRAME));
  Model::FrameIndex frame_idx = model.getFrameId(frame_name);
  pinocchio::Data data(model);
  pinocchio::Data data_ref(model);

  VectorXd q = VectorXd::Ones(model.nq);
  q.middleRows<4> (3).normalize();

  forwardKinematics(model, data, q);
  updateFramePlacements(model, data);

  framesForwardKinematics(model, data_ref, q);

  BOOST_CHECK(data.oMf[frame_idx].isApprox(data_ref.oMf[frame_idx]));
}

BOOST_AUTO_TEST_CASE ( test_update_single_placement )
{
  using namespace Eigen;
  using namespace pinocchio;

  pinocchio::Model model;
  pinocchio::buildModels::humanoidRandom(model);
  Model::Index parent_idx = model.existJointName("rarm2_joint")?model.getJointId("rarm2_joint"):(Model::Index)(model.njoints-1);
  const std::string & frame_name = std::string( model.names[parent_idx]+ "_frame");
  const SE3 & framePlacement = SE3::Random();
  model.addFrame(Frame (frame_name, parent_idx, 0, framePlacement, OP_FRAME));
  Model::FrameIndex frame_idx = model.getFrameId(frame_name);
  pinocchio::Data data(model);
  pinocchio::Data data_ref(model);

  VectorXd q = VectorXd::Ones(model.nq);
  q.middleRows<4> (3).normalize();

  forwardKinematics(model, data, q);
  updateFramePlacement(model, data, frame_idx);

  framesForwardKinematics(model, data_ref, q);

  BOOST_CHECK(data.oMf[frame_idx].isApprox(data_ref.oMf[frame_idx]));
}

BOOST_AUTO_TEST_CASE ( test_velocity )
{
  using namespace Eigen;
  using namespace pinocchio;

  pinocchio::Model model;
  pinocchio::buildModels::humanoidRandom(model);
  Model::Index parent_idx = model.existJointName("rarm2_joint")?model.getJointId("rarm2_joint"):(Model::Index)(model.njoints-1);
  const std::string & frame_name = std::string( model.names[parent_idx]+ "_frame");
  const SE3 & framePlacement = SE3::Random();
  model.addFrame(Frame (frame_name, parent_idx, 0, framePlacement, OP_FRAME));
  Model::FrameIndex frame_idx = model.getFrameId(frame_name);
  pinocchio::Data data(model);

  VectorXd q = VectorXd::Ones(model.nq);
  q.middleRows<4> (3).normalize();
  VectorXd v = VectorXd::Ones(model.nv);
  forwardKinematics(model, data, q, v);

  Motion vf = getFrameVelocity(model, data, frame_idx);

  BOOST_CHECK(vf.isApprox(framePlacement.actInv(data.v[parent_idx])));
}

BOOST_AUTO_TEST_CASE ( test_acceleration )
{
  using namespace Eigen;
  using namespace pinocchio;

  pinocchio::Model model;
  pinocchio::buildModels::humanoidRandom(model);
  Model::Index parent_idx = model.existJointName("rarm2_joint")?model.getJointId("rarm2_joint"):(Model::Index)(model.njoints-1);
  const std::string & frame_name = std::string( model.names[parent_idx]+ "_frame");
  const SE3 & framePlacement = SE3::Random();
  model.addFrame(Frame (frame_name, parent_idx, 0, framePlacement, OP_FRAME));
  Model::FrameIndex frame_idx = model.getFrameId(frame_name);
  pinocchio::Data data(model);

  VectorXd q = VectorXd::Ones(model.nq);
  q.middleRows<4> (3).normalize();
  VectorXd v = VectorXd::Ones(model.nv);
  VectorXd a = VectorXd::Ones(model.nv);
  forwardKinematics(model, data, q, v, a);

  Motion af = getFrameAcceleration(model, data, frame_idx);

  BOOST_CHECK(af.isApprox(framePlacement.actInv(data.a[parent_idx])));
}

BOOST_AUTO_TEST_CASE ( test_get_frame_jacobian )
{
  using namespace Eigen;
  using namespace pinocchio;

  Model model;
  buildModels::humanoidRandom(model);
  Model::Index parent_idx = model.existJointName("rarm2_joint")?model.getJointId("rarm2_joint"):(Model::Index)(model.njoints-1);
  const std::string & frame_name = std::string( model.names[parent_idx]+ "_frame");
  const SE3 & framePlacement = SE3::Random();
  model.addFrame(Frame (frame_name, parent_idx, 0, framePlacement, OP_FRAME));
  BOOST_CHECK(model.existFrame(frame_name));

  pinocchio::Data data(model);
  pinocchio::Data data_ref(model);

  model.lowerPositionLimit.head<7>().fill(-1.);
  model.upperPositionLimit.head<7>().fill( 1.);
  VectorXd q = randomConfiguration(model);
  VectorXd v = VectorXd::Ones(model.nv);

  /// In local frame
  Model::Index idx = model.getFrameId(frame_name);
  const Frame & frame = model.frames[idx];
  BOOST_CHECK(frame.placement.isApprox_impl(framePlacement));
  Data::Matrix6x Jjj(6,model.nv); Jjj.fill(0);
  Data::Matrix6x Jff(6,model.nv); Jff.fill(0);

  computeJointJacobians(model,data,q);
  updateFramePlacement(model, data, idx);
  getFrameJacobian(model,     data,        idx, LOCAL, Jff);
  computeJointJacobians(model,data_ref,q);
  getJointJacobian(model, data_ref, parent_idx, LOCAL, Jjj);

  Motion nu_frame = Motion(Jff*v);
  Motion nu_joint = Motion(Jjj*v);

  const SE3::ActionMatrixType jXf = frame.placement.toActionMatrix();
  Data::Matrix6x Jjj_from_frame(jXf * Jff);
  BOOST_CHECK(Jjj_from_frame.isApprox(Jjj));

  BOOST_CHECK(nu_frame.isApprox(frame.placement.actInv(nu_joint), 1e-12));

  // In world frame
  getFrameJacobian(model,data,idx,WORLD,Jff);
  getJointJacobian(model, data_ref, parent_idx,WORLD, Jjj);
  BOOST_CHECK(Jff.isApprox(Jjj));
}

BOOST_AUTO_TEST_CASE ( test_frame_jacobian )
{
  using namespace Eigen;
  using namespace pinocchio;

  Model model;
  buildModels::humanoidRandom(model);
  Model::Index parent_idx = model.existJointName("rarm2_joint")?model.getJointId("rarm2_joint"):(Model::Index)(model.njoints-1);
  const std::string & frame_name = std::string( model.names[parent_idx]+ "_frame");
  const SE3 & framePlacement = SE3::Random();
  model.addFrame(Frame (frame_name, parent_idx, 0, framePlacement, OP_FRAME));
  BOOST_CHECK(model.existFrame(frame_name));

  pinocchio::Data data(model);
  pinocchio::Data data_ref(model);

  model.lowerPositionLimit.head<7>().fill(-1.);
  model.upperPositionLimit.head<7>().fill( 1.);
  VectorXd q = randomConfiguration(model);
  VectorXd v = VectorXd::Ones(model.nv);

  Model::Index idx = model.getFrameId(frame_name);
  const Frame & frame = model.frames[idx];
  BOOST_CHECK(frame.placement.isApprox_impl(framePlacement));
  Data::Matrix6x Jf(6,model.nv); Jf.fill(0);
  Data::Matrix6x Jf_ref(6,model.nv); Jf_ref.fill(0);

  frameJacobian(model, data_ref, q, idx, Jf);

  computeJointJacobians(model, data_ref, q);
  updateFramePlacement(model,  data_ref, idx);
  getFrameJacobian(model,      data_ref, idx, LOCAL, Jf_ref);

  BOOST_CHECK(Jf.isApprox(Jf_ref));
}

BOOST_AUTO_TEST_CASE ( test_frame_jacobian_time_variation )
{
  using namespace Eigen;
  using namespace pinocchio;

  pinocchio::Model model;
  pinocchio::buildModels::humanoidRandom(model);
  Model::Index parent_idx = model.existJointName("rarm2_joint")?model.getJointId("rarm2_joint"):(Model::Index)(model.njoints-1);
  const std::string & frame_name = std::string( model.names[parent_idx]+ "_frame");
  const SE3 & framePlacement = SE3::Random();
  model.addFrame(Frame (frame_name, parent_idx, 0, framePlacement, OP_FRAME));
  pinocchio::Data data(model);
  pinocchio::Data data_ref(model);

  VectorXd q = randomConfiguration(model, -1 * Eigen::VectorXd::Ones(model.nq), Eigen::VectorXd::Ones(model.nq) );
  VectorXd v = VectorXd::Random(model.nv);
  VectorXd a = VectorXd::Random(model.nv);

  computeJointJacobiansTimeVariation(model,data,q,v);
  updateFramePlacements(model,data);

  forwardKinematics(model,data_ref,q,v,a);
  updateFramePlacements(model,data_ref);

  BOOST_CHECK(isFinite(data.dJ));

  Model::Index idx = model.getFrameId(frame_name);
  const Frame & frame = model.frames[idx];
  BOOST_CHECK(frame.placement.isApprox_impl(framePlacement));

  Data::Matrix6x J(6,model.nv); J.fill(0.);
  Data::Matrix6x dJ(6,model.nv); dJ.fill(0.);

  // Regarding to the world origin
  getFrameJacobian(model,data,idx,WORLD,J);
  getFrameJacobianTimeVariation(model,data,idx,WORLD,dJ);

  Motion v_idx(J*v);
  const Motion & v_ref_local = frame.placement.actInv(data_ref.v[parent_idx]);
  const Motion & v_ref = data_ref.oMf[idx].act(v_ref_local);
  BOOST_CHECK(v_idx.isApprox(v_ref));

  Motion a_idx(J*a + dJ*v);
  const Motion & a_ref_local = frame.placement.actInv(data_ref.a[parent_idx]);
  const Motion & a_ref = data_ref.oMf[idx].act(a_ref_local);
  BOOST_CHECK(a_idx.isApprox(a_ref));

  J.fill(0.);  dJ.fill(0.);
  // Regarding to the local frame
  getFrameJacobian(model,data,idx,LOCAL,J);
  getFrameJacobianTimeVariation(model,data,idx,LOCAL,dJ);

  v_idx = (Motion::Vector6)(J*v);
  BOOST_CHECK(v_idx.isApprox(v_ref_local));

  a_idx = (Motion::Vector6)(J*a + dJ*v);
  BOOST_CHECK(a_idx.isApprox(a_ref_local));

  // compare to finite differencies
  {
    Data data_ref(model), data_ref_plus(model);

    const double alpha = 1e-8;
    Eigen::VectorXd q_plus(model.nq);
    q_plus = integrate(model,q,alpha*v);

    //data_ref
    Data::Matrix6x J_ref_world(6,model.nv), J_ref_local(6,model.nv);
    J_ref_world.fill(0.);     J_ref_local.fill(0.);
    computeJointJacobians(model,data_ref,q);
    updateFramePlacements(model,data_ref);
    const SE3 & oMf_q = data_ref.oMf[idx];
    getFrameJacobian(model,data_ref,idx,WORLD,J_ref_world);
    getFrameJacobian(model,data_ref,idx,LOCAL,J_ref_local);

    //data_ref_plus
    Data::Matrix6x J_ref_plus_world(6,model.nv), J_ref_plus_local(6,model.nv);
    J_ref_plus_world.fill(0.);    J_ref_plus_local.fill(0.);
    computeJointJacobians(model,data_ref_plus,q_plus);
    updateFramePlacements(model,data_ref_plus);
    const SE3 & oMf_qplus = data_ref_plus.oMf[idx];
    getFrameJacobian(model,data_ref_plus,idx,WORLD,J_ref_plus_world);
    getFrameJacobian(model,data_ref_plus,idx,LOCAL,J_ref_plus_local);

    //Move J_ref_plus_local to reference frame
    J_ref_plus_local = (oMf_q.inverse()*oMf_qplus).toActionMatrix()*(J_ref_plus_local);

    Data::Matrix6x dJ_ref_world(6,model.nv), dJ_ref_local(6,model.nv);
    dJ_ref_world.fill(0.);        dJ_ref_local.fill(0.);
    dJ_ref_world = (J_ref_plus_world - J_ref_world)/alpha;
    dJ_ref_local = (J_ref_plus_local - J_ref_local)/alpha;

    //data
    computeJointJacobiansTimeVariation(model,data,q,v);
    forwardKinematics(model,data,q,v);
    updateFramePlacements(model,data);
    Data::Matrix6x dJ_world(6,model.nv), dJ_local(6,model.nv);
    dJ_world.fill(0.);    dJ_local.fill(0.);
    getFrameJacobianTimeVariation(model,data,idx,WORLD,dJ_world);
    getFrameJacobianTimeVariation(model,data,idx,LOCAL,dJ_local);

    BOOST_CHECK(dJ_world.isApprox(dJ_ref_world,sqrt(alpha)));
    BOOST_CHECK(dJ_local.isApprox(dJ_ref_local,sqrt(alpha)));
  }
}

BOOST_AUTO_TEST_SUITE_END ()

max_stars_repo_name: mkatliar/pinocchio

max_issues_repo_name: mkatliar/pinocchio

max_forks_repo_name: mkatliar/pinocchio


text:
#include <boost/foreach.hpp>
#include "stdhdrs.h"

#include "Log.h"
#include "Server.h"
#include "Friend.h"
#include "CmdMsg.h"
#include "doFunc.h"

CFriendMember::CFriendMember()
	: m_strFriendName(MAX_CHAR_NAME_LENGTH + 1)
{
	m_nChaIndex = -1;
//	m_nLastConnectDate = 0; //최근접속일자
	m_nJobType = -1; //종족
	m_nCondition = 	MSG_FRIEND_CONDITION_OFFLINE; //상태
	m_groupIndex = 0;
	m_pPc = NULL;
}

CFriendMember::~CFriendMember()
{
}
/*
int CFriendMember::GetChaIndex()
{
	return m_nChaIndex;
}
*/
void CFriendMember::SetChaIndex(int chaindex)
{
	m_nChaIndex = chaindex;
}

void CFriendMember::SetName(const char* name)
{
	if (name)
	{
		m_strFriendName = name;
	}
	else
	{
		m_strFriendName = "";
	}
}

/*
void CFriendMember::SetFriendListIndex(int listindex)
{
	m_nListIndex = listindex;
}
*/

void CFriendMember::SetCondition(int condition)
{
	m_nCondition = condition;
}

void CFriendMember::SetJob(int job)
{
	m_nJobType = job;
}

CFriend::CFriend()
{
	memset(m_pFriendMember, 0, sizeof(CFriendMember*) * FRIEND_MAX_MEMBER);
//	memset(m_aFriendMeChaIndex, -1, sizeof(int) * FRIEND_MAX_MEMBER);
	m_nFriendCnt = 0;
}

CFriend::~CFriend()
{
	for (int i = 0; i < FRIEND_MAX_MEMBER; i++)
	{
		if (m_pFriendMember[i])
			delete m_pFriendMember[i];
	}
	memset(m_pFriendMember, 0, sizeof(CFriendMember*) * FRIEND_MAX_MEMBER);
	m_nFriendCnt=0;
}

CFriendMember* CFriend::FindFriendMember(int chaindex)
{
	for (int i = 0; i < FRIEND_MAX_MEMBER; i++)
	{
		if (m_pFriendMember[i])
		{
			if (m_pFriendMember[i]->GetChaIndex() == chaindex)
				return m_pFriendMember[i];
		}
	}
	return NULL;
}

CFriendMember* CFriend::FindFriendMember(const char* name)
{
	for (int i = 0; i < FRIEND_MAX_MEMBER; i++)
	{
		if (m_pFriendMember[i])
		{
			if ( strcmp(m_pFriendMember[i]->GetName(), name) == 0 )
				return m_pFriendMember[i];
		}
	}
	return NULL;
}

//친구등록
int CFriend::AddFriend(int chaindex, const char* name, int job, int nCon, int groupIndex)
{
	if (m_nFriendCnt >= FRIEND_MAX_MEMBER)
	{
		return -1;
	}

	CFriendMember* pFriendMember;
	for (int i = 0; i < FRIEND_MAX_MEMBER; i++)
	{
		if (m_pFriendMember[i] == NULL)
		{
			pFriendMember = new CFriendMember();
			pFriendMember->SetChaIndex(chaindex);
			pFriendMember->SetName(name);
			pFriendMember->SetCondition(nCon);
			pFriendMember->SetJob(job);
			map_t::iterator it = m_groupList.find(groupIndex);
			if (it != m_groupList.end())
			{
				pFriendMember->SetGroup(groupIndex);
			}
			else
			{
				pFriendMember->SetGroup(0);
			}

			m_pFriendMember[i] = pFriendMember;
			m_nFriendCnt++;
			return i;	//리스트의 인덱스가 되겠지..
		}
	}

	return -1;
}

//친구삭제
void CFriend::RemoveFriend(int chaindex)
{
	for (int i = 0; i < FRIEND_MAX_MEMBER; i++)
	{
		if (m_pFriendMember[i])
		{
			if (m_pFriendMember[i]->GetChaIndex() == chaindex)
			{
				m_pFriendMember[i] = NULL;
				m_nFriendCnt--;
				return;
			}
		}
	}
}

/*
void CFriend::SetFriendList(int index, int friendindex, const char* name, int job, int condition)
{
	if(m_pFriendMember)
	{
		SetFriendCount(index+1);
//		m_pFriendMember[index]->SetFriendListIndex(index);
		m_pFriendMember[index]->SetChaIndex(friendindex);
		m_pFriendMember[index]->SetName(name);
		m_pFriendMember[index]->SetJob(job);
		m_pFriendMember[index]->SetCondition(condition);
	}
}
*/

void CFriend::SendToAll(CNetMsg::SP& msg)
{
	int i;
	for (i = 0; i < FRIEND_MAX_MEMBER; i++)
	{
		if (m_pFriendMember[i])
		{
			if (m_pFriendMember[i]->GetPC())
			{
				SEND_Q(msg, m_pFriendMember[i]->GetPC()->m_desc);
			}
		}
	}
}

int CFriend::addGroup(const char* gName, int gIndex)
{
	// 문자열 검사
	if (strinc(gName, "'") || strinc(gName, "%") || strinc(gName, " "))
	{
		return -1;
	}

	// 최대 갯수가 넘으면
	if (this->m_groupList.size() >= FRIEND_GROUP_MAX_NUMBER)
	{
		return -1;
	}

	// 동일 그룹명이 있는지 검색
	BOOST_FOREACH(map_t::value_type& p, this->m_groupList)
	{
		CFriendGroup* pGroup = p.second;
		if( strcmp( (const char*) pGroup->GetGroupName(), gName) == 0 )
			return -1;
	}

	// 디비에서 읽은 데이터이면 바로 입력 처리
	CFriendGroup* pGroup = NULL;
	if (gIndex != -1)
	{
		pGroup = new CFriendGroup(gIndex, gName);
	}
	else
	{
		// 최종 인덱스를 구함
		if (this->m_groupList.empty())
		{
			gIndex = 0;
		}
		else
		{
			// std::map은 자동으로 오름차순으로 정리되므로 맨마지막 key값이 가장 큼
			map_t::iterator it = --m_groupList.end();
			gIndex = (it->first) + 1;
		}

		pGroup = new CFriendGroup(gIndex, gName);
	}

	m_groupList.insert(map_t::value_type(gIndex, pGroup));

	return gIndex;
}

//void	delGroup(int gIndex);
//int		changeGroupName(int gIndex, const char* newName);
//int		moveGroupFriend(int gIndex, const char* charName);

int CFriend::delGroup(int gIndex)
{
	map_t::iterator it = m_groupList.find(gIndex);
	if (it != m_groupList.end())
	{
		m_groupList.erase(it);
		return gIndex;
	}

	return -1;
}

int CFriend::changeGroupName(int gIndex, const char* newName)
{
	if( strinc(newName, "'") || strinc(newName, "%") || strinc(newName, " ") )
	{
		return -1;
	}

	map_t::iterator it = m_groupList.find(gIndex);
	if (it != m_groupList.end())
	{
		CFriendGroup* fgroup = it->second;
		fgroup->SetName(newName);
		return gIndex;
	}

	return -1;
}

int CFriend::moveGroupFriend(int gIndex, int charIndex)
{
	CFriendMember* fmember = this->FindFriendMember(charIndex);

	if( fmember )
	{
		fmember->SetGroup(gIndex);
		return gIndex;
	}

	return -1;
}

void CFriend::GetGroupIndexString(CLCString& gIndexList, CLCString& gNameList)
{
	CFriendGroup* fgroup = NULL;
	map_t::iterator it = m_groupList.begin();
	map_t::iterator endit = m_groupList.end();
	char tmpbuf[1024] = {0,};
	for(; it != endit; ++it)
	{
		CFriendGroup* fgroup = it->second;

		sprintf(tmpbuf, "%d", fgroup->GetGroupIndex() );
		gIndexList += tmpbuf;
		gIndexList += " ";
		gNameList += fgroup->GetGroupName();
		gNameList += " ";
	}
}

CFriendGroup::~CFriendGroup()
{
}
//

max_stars_repo_name: openlastchaos/lastchaos-source-server

max_issues_repo_name: openlastchaos/lastchaos-source-server

max_forks_repo_name: openlastchaos/lastchaos-source-server


text: #include <boost/geometry/formulas/elliptic_arc_length.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/*
 [auto_generated]
 libs/numeric/odeint/examples/bind_member_functions.hpp

 [begin_description]
 tba.
 [end_description]

 Copyright 2012 Karsten Ahnert
 Copyright 2012 Mario Mulansky

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

#include <iostream>

#include <boost/numeric/odeint.hpp>

namespace odeint = boost::numeric::odeint;

typedef boost::array< double , 3 > state_type;

//[ ode_wrapper
template< class Obj , class Mem >
class ode_wrapper
{
    Obj m_obj;
    Mem m_mem;

public:

    ode_wrapper( Obj obj , Mem mem ) : m_obj( obj ) , m_mem( mem ) { }

    template< class State , class Deriv , class Time >
    void operator()( const State &x , Deriv &dxdt , Time t )
    {
        (m_obj.*m_mem)( x , dxdt , t );
    }
};

template< class Obj , class Mem >
ode_wrapper< Obj , Mem > make_ode_wrapper( Obj obj , Mem mem )
{
    return ode_wrapper< Obj , Mem >( obj , mem );
}
//]

template< class Obj , class Mem >
class observer_wrapper
{
    Obj m_obj;
    Mem m_mem;

public:

    observer_wrapper( Obj obj , Mem mem ) : m_obj( obj ) , m_mem( mem ) { }

    template< class State , class Time >
    void operator()( const State &x , Time t )
    {
        (m_obj.*m_mem)( x , t );
    }
};

template< class Obj , class Mem >
observer_wrapper< Obj , Mem > make_observer_wrapper( Obj obj , Mem mem )
{
    return observer_wrapper< Obj , Mem >( obj , mem );
}

//[ bind_member_function
struct lorenz
{
    void ode( const state_type &x , state_type &dxdt , double t ) const
    {
        dxdt[0] = 10.0 * ( x[1] - x[0] );
        dxdt[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -8.0 / 3.0 * x[2] + x[0] * x[1];
    }
};

int main( int argc , char *argv[] )
{
    using namespace boost::numeric::odeint;
    state_type x = {{ 10.0 , 10.0 , 10.0 }};
    integrate_const( runge_kutta4< state_type >() , make_ode_wrapper( lorenz() , &lorenz::ode ) ,
                     x , 0.0 , 10.0 , 0.01 );
    return 0;
}
//]

/*
struct lorenz
{
    void ode( const state_type &x , state_type &dxdt , double t ) const
    {
        dxdt[0] = 10.0 * ( x[1] - x[0] );
        dxdt[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -8.0 / 3.0 * x[2] + x[0] * x[1];
    }

    void obs( const state_type &x , double t ) const
    {
        std::cout << t << " " << x[0] << " " << x[1] << " " << x[2] << "\n";
    }
};

int main( int argc , char *argv[] )
{
    using namespace boost::numeric::odeint;

    state_type x = {{ 10.0 , 10.0 , 10.0 }};
    integrate_const( runge_kutta4< state_type >() ,
                     make_ode_wrapper( lorenz() , &lorenz::ode ) ,
                     x , 0.0 , 10.0 , 0.01 ,
                     make_observer_wrapper( lorenz() , &lorenz::obs ) );

    return 0;
}
*/

max_stars_repo_name: cpp-pm/boost

max_issues_repo_name: cpp-pm/boost

max_forks_repo_name: cpp-pm/boost


text:
// Copyright (C) 2019-2020 Zilliz. All rights reserved.
//
// 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 "common/Schema.h"
#include <google/protobuf/text_format.h>
#include <boost/algorithm/string/case_conv.hpp>
#include <boost/lexical_cast.hpp>
#include "common/SystemProperty.h"
#include <optional>

namespace milvus {

using std::string;
static std::map<string, string>
RepeatedKeyValToMap(const google::protobuf::RepeatedPtrField<proto::common::KeyValuePair>& kvs) {
    std::map<string, string> mapping;
    for (auto& kv : kvs) {
        AssertInfo(!mapping.count(kv.key()), "repeat key(" + kv.key() + ") in protobuf");
        mapping.emplace(kv.key(), kv.value());
    }
    return mapping;
}

std::shared_ptr<Schema>
Schema::ParseFrom(const milvus::proto::schema::CollectionSchema& schema_proto) {
    auto schema = std::make_shared<Schema>();
    // schema->set_auto_id(schema_proto.autoid());

    // NOTE: only two system

    for (const milvus::proto::schema::FieldSchema& child : schema_proto.fields()) {
        auto field_offset = FieldOffset(schema->size());
        auto field_id = FieldId(child.fieldid());
        auto name = FieldName(child.name());

        if (field_id.get() < 100) {
            // system field id
            auto is_system = SystemProperty::Instance().SystemFieldVerify(name, field_id);
            AssertInfo(is_system,
                       "invalid system type: name(" + name.get() + "), id(" + std::to_string(field_id.get()) + ")");
            continue;
        }

        auto data_type = DataType(child.data_type());

        if (datatype_is_vector(data_type)) {
            auto type_map = RepeatedKeyValToMap(child.type_params());
            auto index_map = RepeatedKeyValToMap(child.index_params());

            AssertInfo(type_map.count("dim"), "dim not found");
            auto dim = boost::lexical_cast<int64_t>(type_map.at("dim"));
            if (!index_map.count("metric_type")) {
                schema->AddField(name, field_id, data_type, dim, std::nullopt);
            } else {
                auto metric_type = GetMetricType(index_map.at("metric_type"));
                schema->AddField(name, field_id, data_type, dim, metric_type);
            }
        } else {
            schema->AddField(name, field_id, data_type);
        }

        if (child.is_primary_key()) {
            AssertInfo(!schema->get_primary_key_offset().has_value(), "repetitive primary key");
            Assert(!schema_proto.autoid());
            schema->set_primary_key(field_offset);
        }
    }
    if (schema->get_is_auto_id()) {
        AssertInfo(!schema->get_primary_key_offset().has_value(), "auto id mode: shouldn't have primary key");
    } else {
        AssertInfo(schema->get_primary_key_offset().has_value(), "primary key should be specified when autoId is off");
    }

    return schema;
}

const FieldMeta FieldMeta::RowIdMeta(FieldName("RowID"), FieldId(0), DataType::INT64);

}  // namespace milvus

max_stars_repo_name: mfkiwl/milvus

max_issues_repo_name: mfkiwl/milvus

max_forks_repo_name: mfkiwl/milvus