text:
#include <mesh_array_io_obj.h>

#include <tiny.h>

#include <util_log.h>

#include <boost/algorithm/string.hpp>  //--- needed for split and is_any_of function
#include <boost/lexical_cast.hpp>      //--- needed for lexical cast function
#include <boost/filesystem.hpp>        //--- needed for boost::filesystem::path

#include <cassert>
#include <sstream>
#include <fstream>
#include <vector>
#include <stdexcept>

namespace mesh_array
{

  template<typename T>
  void read_obj(
                std::string const & filename
                , T3Mesh & mesh
                , VertexAttribute<T,T3Mesh> & X
                , VertexAttribute<T,T3Mesh> & Y
                , VertexAttribute<T,T3Mesh> & Z
                )
  {
    boost::filesystem::path  mypath(filename);

    std::string const ext = boost::filesystem::extension( mypath );

    assert( ext == ".obj" || !"read_obj(): illegal file name extension");

    if( ext != ".obj")
    {
      util::Log logging;

      logging << "mesh_array::read_obj(...): illegal file name extension" << util::Log::newline();

      throw std::logic_error("Filename extension must be obj");

      return;
    }

    bool const did_exist = boost::filesystem::exists( mypath );

    assert( did_exist || !"read_obj(): file did not exist");

    if (!did_exist)
    {
      util::Log logging;

      logging << "mesh_array::read_obj(...): file did not exist" << util::Log::newline();

      throw std::logic_error("Filename did not exist");

      return;
    }

    typedef std::vector<T>       fcontainer;
    typedef std::vector<size_t>  icontainer;

    std::ifstream file(filename.c_str());

    assert(file.is_open() || !"read_obj(): could not open file");

    fcontainer px;
    fcontainer py;
    fcontainer pz;
    fcontainer nx;
    fcontainer ny;
    fcontainer nz;
    fcontainer tx;
    fcontainer ty;
    fcontainer tz;
    icontainer I;
    icontainer J;
    icontainer K;

    while(!file.eof())
    {
      std::string line;

      if(! getline(file,line) )
        break;

      std::stringstream stream(line);

      std::string type;
      stream >> type;

      if(line[0] == '#')
      {
        continue;
      }
      else if(type == "g")
      {
        std::string group_name;
        stream >> group_name;

        // create new mesh with name=group name!!!

        continue;
      }
      else if(type == "usemtl")
      {
        std::string material_name;
        stream >> material_name;

        // material = materials[ material_name ];

        continue;
      }
      else if(type == "mtllib")
      {
        std::string material_filename;
        stream >> material_filename;

        // read materials from material_filename

        continue;
      }
      else if(type == "v")
      {
        T value;

        stream >> value;
        assert(is_number(value) || !"read_obj(): nan");
        assert(is_finite(value) || !"read_obj(): inf");
        px.push_back(value);

        stream >> value;
        assert(is_number(value) || !"read_obj(): nan");
        assert(is_finite(value) || !"read_obj(): inf");
        py.push_back(value);

        stream >> value;
        assert(is_number(value) || !"read_obj(): nan");
        assert(is_finite(value) || !"read_obj(): inf");
        pz.push_back(value);
      }
      else if(type == "vt")
      {
        T value;

        stream >> value;
        assert(is_number(value) || !"read_obj(): nan");
        assert(is_finite(value) || !"read_obj(): inf");
        tx.push_back(value);

        if(!stream.eof())
        {
          stream >> value;
          assert(is_number(value) || !"read_obj(): nan");
          assert(is_finite(value) || !"read_obj(): inf");
          ty.push_back(value);
        }

        if(!stream.eof())
        {
          stream >> value;
          assert(is_number(value) || !"read_obj(): nan");
          assert(is_finite(value) || !"read_obj(): inf");
          tz.push_back(value);
        }
      }
      else if(type == "vn")
      {
        T value;

        stream >> value;
        assert(is_number(value) || !"read_obj(): nan");
        assert(is_finite(value) || !"read_obj(): inf");
        nx.push_back(value);

        stream >> value;
        assert(is_number(value) || !"read_obj(): nan");
        assert(is_finite(value) || !"read_obj(): inf");
        ny.push_back(value);

        stream >> value;
        assert(is_number(value) || !"read_obj(): nan");
        assert(is_finite(value) || !"read_obj(): inf");
        nz.push_back(value);
      }
      else if(type == "f")
      {
        icontainer face;

        //--- parse:
        //---  f v0_idx//vn0_idx//vt0_idx  ...  vN_idx//vnN_idx//vtN_idx
        while(!stream.eof())
        {
          std::string raw;
          stream >> raw;

          // NOTE: hd 20060605 - hack to let linux version gracefully read windows-formatted files.
          if (raw.empty())
            break;

          std::vector<std::string> split_vector;
          boost::algorithm::split(
                                  split_vector
                                  , raw
                                  , boost::algorithm::is_any_of("/")
                                  , boost::algorithm::token_compress_on
                                  );

          assert(!split_vector.empty() || !"obj_read() : Internal error, no indices on face");

          //size_t const cnt = split_vector.size();  // if 1 then only coords, 2 coords+normal, 3 coords+normal+texture

          //std::vector<std::string>::iterator value = split_vector.begin();
          //std::vector<std::string>::iterator end   = split_vector.end();
          //for(;value != end; ++value)
          //  face.push_back(boost::lexical_cast<size_t>(*value));  // if we wanted to extract texture and normal indices too

          size_t const idx = boost::lexical_cast<size_t>( split_vector[0] );
          face.push_back( idx ); // we only care about vertex indices

        }

        // We assume that a polygonal face is convex so we perform ear clipping algorithm
        // if more than 3 vertices in polygon face
        size_t const triangles =  face.size() - 2u;
        for(size_t t = 0u;t<triangles;++t)
        {
          size_t const i = face[0u]     - 1u;   // obj files uses one-based indexing so we convert to zero-based indexing
          size_t const j = face[t+1u]   - 1u;
          size_t const k = face[t+2u]   - 1u;

          I.push_back(i);
          J.push_back(j);
          K.push_back(k);
        }

      }
      else
      {
        util::Log logging;

        logging << "read_obj_materials(): skipping "<< type << util::Log::newline();
      }
    }
    file.close();

    size_t const cntV = px.size();
    size_t const cntT = I.size();

    assert(cntV > 0u || !"read_obj(): no vertices");
    assert(cntT > 0u || !"read_obj(): no triangles");

    mesh.clear();

    mesh.set_capacity(cntV, cntT);

    X.bind(mesh);
    Y.bind(mesh);
    Z.bind(mesh);

    std::vector<Vertex> vertices (cntV);

    for(size_t idx=0u; idx<cntV; ++idx)
    {
      Vertex v = mesh.push_vertex();

      X(v) = px[idx];
      Y(v) = py[idx];
      Z(v) = pz[idx];

      vertices[idx] = v;
    }

    for(size_t idx=0u; idx<cntT; ++idx)
    {
      size_t const i = I[idx];
      size_t const j = J[idx];
      size_t const k = K[idx];

      mesh.push_triangle(vertices[i],vertices[j],vertices[k]);
    }
  }

  template
  void read_obj<float>(
                std::string const & filename
                , T3Mesh & mesh
                , VertexAttribute<float,T3Mesh> & X
                , VertexAttribute<float,T3Mesh> & Y
                , VertexAttribute<float,T3Mesh> & Z
                );

  template
  void read_obj<double>(
                std::string const & filename
                , T3Mesh & mesh
                , VertexAttribute<double,T3Mesh> & X
                , VertexAttribute<double,T3Mesh> & Y
                , VertexAttribute<double,T3Mesh> & Z
                );

  template<typename T>
  void write_obj(
                 std::string const & filename
                 , T3Mesh const & mesh
                 , VertexAttribute<T,T3Mesh> const & X
                 , VertexAttribute<T,T3Mesh> const & Y
                 , VertexAttribute<T,T3Mesh> const & Z
                 )
  {
    boost::filesystem::path  mypath(filename);

    std::string const ext = boost::filesystem::extension( mypath );

    assert( ext == ".obj" || !"write_obj(): illegal file name extension");

    std::ofstream file(filename.c_str());

    assert( file.is_open() || !"write_obj(): could not open file");

    file << "g group0"  << std::endl;

    file.precision(30);

    size_t const cntV = mesh.vertex_size();

    for(size_t idx = 0u; idx < cntV; ++idx)
    {
      Vertex const & vertex = mesh.vertex( idx );

      file << "v " << X(vertex) << ' ' << Y(vertex)  << ' ' << Z(vertex)  << std::endl;
      //file << "vt " << TX(vertex) << ' ' << TY(vertex)  << ' ' << TZ(vertex)  << std::endl;
      //file << "vn " << NX(vertex) << ' ' << NY(vertex)  << ' ' << NZ(vertex)  << std::endl;
    }

    size_t const cntT = mesh.triangle_size();
    for(size_t idx = 0u; idx <cntT; ++idx)
    {
      Triangle const & triangle = mesh.triangle(idx);

      file << "f";

      size_t const i = triangle.i() + 1u; //--- need one-based indices!!!
      size_t const j = triangle.j() + 1u; //--- need one-based indices!!!
      size_t const k = triangle.k() + 1u; //--- need one-based indices!!!

      file << " "
      << i   /*   << "/" << ni << "/" ti */
      << " "
      << j   /*   << "/" << nj << "/" tj */
      << " "
      << k   /*   << "/" << nk << "/" tk */
      << std::endl;
    }
    file.flush();
    file.close();
  }

  template
  void write_obj<float>(
                 std::string const & filename
                 , T3Mesh const & mesh
                 , VertexAttribute<float,T3Mesh> const & X
                 , VertexAttribute<float,T3Mesh> const & Y
                 , VertexAttribute<float,T3Mesh> const & Z
                 );

  template
  void write_obj<double>(
                 std::string const & filename
                 , T3Mesh const & mesh
                 , VertexAttribute<double,T3Mesh> const & X
                 , VertexAttribute<double,T3Mesh> const & Y
                 , VertexAttribute<double,T3Mesh> const & Z
                 );

} // namespace mesh_array

max_stars_repo_name: diku-dk/PROX

max_issues_repo_name: erleben/matchstick

max_forks_repo_name: erleben/matchstick


text:
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements.  See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership.  The ASF licenses this file
// to you under the Apache License, Version 2.0 (the
// "License"); you may not use this file except in compliance
// with the License.  You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License

#ifndef __WINDOWS__
#include <unistd.h>
#endif // __WINDOWS__

#include <sys/stat.h>

#include <algorithm>
#include <map>
#include <string>
#include <vector>

#include <boost/shared_array.hpp>

#include <process/defer.hpp>
#include <process/deferred.hpp> // TODO(benh): This is required by Clang.
#include <process/dispatch.hpp>
#include <process/future.hpp>
#include <process/help.hpp>
#include <process/http.hpp>
#include <process/io.hpp>
#include <process/mime.hpp>
#include <process/process.hpp>

#include <stout/error.hpp>
#include <stout/hashmap.hpp>
#include <stout/json.hpp>
#include <stout/lambda.hpp>
#include <stout/none.hpp>
#include <stout/numify.hpp>
#include <stout/option.hpp>
#include <stout/os.hpp>
#include <stout/path.hpp>
#include <stout/result.hpp>
#include <stout/stringify.hpp>
#include <stout/strings.hpp>
#include <stout/try.hpp>

#include <stout/os/constants.hpp>

#include "common/http.hpp"

#include "files/files.hpp"

#include "logging/logging.hpp"

namespace http = process::http;
namespace io = process::io;
namespace mime = process::mime;

using http::BadRequest;
using http::Forbidden;
using http::InternalServerError;
using http::NotFound;
using http::OK;

using mesos::Authorizer;

using process::AUTHENTICATION;
using process::AUTHORIZATION;
using process::defer;
using process::DESCRIPTION;
using process::Failure;
using process::Future;
using process::HELP;
using process::Process;
using process::TLDR;
using process::wait; // Necessary on some OS's to disambiguate.

using std::list;
using std::map;
using std::string;
using std::vector;

namespace mesos {
namespace internal {

class FilesProcess : public Process<FilesProcess>
{
public:
  FilesProcess(const Option<string>& _authenticationRealm,
               const Option<Authorizer*>& _authorizer);

  // Files implementation.
  Future<Nothing> attach(
      const string& path,
      const string& name,
      const Option<lambda::function<Future<bool>(const Option<string>&)>>&
          authorized);

  void detach(const string& name);

  Future<Try<list<FileInfo>, FilesError>> browse(
    const string& path,
    const Option<string>& principal);

protected:
  virtual void initialize();

private:
  // Resolves the virtual path to an actual path.
  // Returns the actual path if found.
  // Returns None if the file is not found.
  // Returns Error if we find the file but it cannot be resolved or it breaks
  // out of the chroot.
  Result<string> resolve(const string& path);

  Future<bool> authorize(string requestedPath, const Option<string>& principal);

  // HTTP endpoints.

  // Returns a file listing for a directory.
  // Requests have the following parameters:
  //   path: The directory to browse. Required.
  // The response will contain a list of JSON files and directories contained
  // in the path (see `FileInfo` model override for the format).
  Future<http::Response> _browse(
      const http::Request& request,
      const Option<string>& principal);

  // Reads data from a file at a given offset and for a given length.
  // See the jquery pailer for the expected behavior.
  Future<http::Response> read(
      const http::Request& request,
      const Option<string>& principal);

  Future<http::Response> _read(
      off_t offset,
      Option<size_t> length,
      const string& path,
      const Option<string>& jsonp);

  // Returns the raw file contents for a given path.
  // Requests have the following parameters:
  //   path: The directory to browse. Required.
  Future<http::Response> download(
      const http::Request& request,
      const Option<string>& principal);

  Future<http::Response> _download(const string& path);

  // Returns the internal virtual path mapping.
  Future<http::Response> debug(
      const http::Request& request,
      const Option<string>& principal);

  const static string BROWSE_HELP;
  const static string READ_HELP;
  const static string DOWNLOAD_HELP;
  const static string DEBUG_HELP;

  hashmap<string, string> paths;

  // Set of authorization functions. They will be called whenever
  // access to the path used as key is requested, and will pass
  // as parameter the principal returned by the HTTP authenticator.
  hashmap<string, lambda::function<Future<bool>(const Option<string>&)>>
      authorizations;

  // The authentication realm, if any, into which this process'
  // endpoints will be installed.
  Option<string> authenticationRealm;

  // FilesProcess needs an authorizer object to add authorization in
  // `/files/debug` endpoint.
  Option<Authorizer*> authorizer;
};

FilesProcess::FilesProcess(
    const Option<string>& _authenticationRealm,
    const Option<Authorizer*>& _authorizer)
  : ProcessBase("files"),
    authenticationRealm(_authenticationRealm),
    authorizer(_authorizer) {}

void FilesProcess::initialize()
{
  if (authenticationRealm.isSome()) {
    // TODO(ijimenez): Remove these endpoints at the end of the
    // deprecation cycle on 0.26.
    route("/browse.json",
          authenticationRealm.get(),
          FilesProcess::BROWSE_HELP,
          &FilesProcess::_browse);
    route("/read.json",
          authenticationRealm.get(),
          FilesProcess::READ_HELP,
          &FilesProcess::read);
    route("/download.json",
          authenticationRealm.get(),
          FilesProcess::DOWNLOAD_HELP,
          &FilesProcess::download);
    route("/debug.json",
          authenticationRealm.get(),
          FilesProcess::DEBUG_HELP,
          &FilesProcess::debug);

    route("/browse",
          authenticationRealm.get(),
          FilesProcess::BROWSE_HELP,
          &FilesProcess::_browse);
    route("/read",
          authenticationRealm.get(),
          FilesProcess::READ_HELP,
          &FilesProcess::read);
    route("/download",
          authenticationRealm.get(),
          FilesProcess::DOWNLOAD_HELP,
          &FilesProcess::download);
    route("/debug",
          authenticationRealm.get(),
          FilesProcess::DEBUG_HELP,
          &FilesProcess::debug);
  } else {
    // TODO(ijimenez): Remove these endpoints at the end of the
    // deprecation cycle on 0.26.
    route("/browse.json",
          FilesProcess::BROWSE_HELP,
          lambda::bind(&FilesProcess::_browse, this, lambda::_1, None()));
    route("/read.json",
          FilesProcess::READ_HELP,
          lambda::bind(&FilesProcess::read, this, lambda::_1, None()));
    route("/download.json",
          FilesProcess::DOWNLOAD_HELP,
          lambda::bind(&FilesProcess::download, this, lambda::_1, None()));
    route("/debug.json",
          FilesProcess::DEBUG_HELP,
          lambda::bind(&FilesProcess::debug, this, lambda::_1, None()));

    route("/browse",
          FilesProcess::BROWSE_HELP,
          lambda::bind(&FilesProcess::_browse, this, lambda::_1, None()));
    route("/read",
          FilesProcess::READ_HELP,
          lambda::bind(&FilesProcess::read, this, lambda::_1, None()));
    route("/download",
          FilesProcess::DOWNLOAD_HELP,
          lambda::bind(&FilesProcess::download, this, lambda::_1, None()));
    route("/debug",
          FilesProcess::DEBUG_HELP,
          lambda::bind(&FilesProcess::debug, this, lambda::_1, None()));
  }
}

Future<Nothing> FilesProcess::attach(
    const string& path,
    const string& name,
    const Option<lambda::function<Future<bool>(const Option<string>&)>>&
        authorized)
{
  Result<string> result = os::realpath(path);

  if (!result.isSome()) {
    return Failure(
        "Failed to get realpath of '" + path + "': " +
        (result.isError()
         ? result.error()
         : "No such file or directory"));
  }

  // Make sure we have permissions to read the file/dir.
  Try<bool> access = os::access(result.get(), R_OK);

  if (access.isError() || !access.get()) {
    return Failure("Failed to access '" + path + "': " +
                   (access.isError() ? access.error() : "Access denied"));
  }

  // To simplify the read/browse logic, strip any trailing / from the name.
  string cleanedName = strings::remove(name, "/", strings::SUFFIX);

  // TODO(bmahler): Do we want to always wipe out the previous path?
  paths[cleanedName] = result.get();

  if (authorized.isSome()) {
    authorizations[cleanedName] = authorized.get();
  }

  return Nothing();
}

void FilesProcess::detach(const string& name)
{
  paths.erase(name);
  authorizations.erase(name);
}

const string FilesProcess::BROWSE_HELP = HELP(
    TLDR(
        "Returns a file listing for a directory."),
    DESCRIPTION(
        "Lists files and directories contained in the path as",
        "a JSON object.",
        "",
        "Query parameters:",
        "",
        ">        path=VALUE          The path of directory to browse."),
    AUTHENTICATION(true),
    AUTHORIZATION(
        "Browsing files requires that the request principal is",
        "authorized to do so for the target virtual file path.",
        "",
        "Authorizers may categorize different virtual paths into",
        "different ACLs, e.g. logs in one and task sandboxes in",
        "another.",
        "",
        "See authorization documentation for details."));

Future<bool> FilesProcess::authorize(
    string requestedPath,
    const Option<string>& principal)
{
  // The path may contain a trailing forward slash. Since we store the
  // authorization callbacks without the trailing slash, we must remove it here,
  // if present.
  string trimmedPath = strings::remove(requestedPath, "/", strings::SUFFIX);

  if (authorizations.count(trimmedPath) > 0) {
    return authorizations[trimmedPath](principal);
  }

  do {
    requestedPath = Path(requestedPath).dirname();

    if (authorizations.count(requestedPath) > 0) {
      return authorizations[requestedPath](principal);
    }
  } while (Path(requestedPath).dirname() != requestedPath);

  // requestedPath does not require authorization.
  return true;
}

Future<http::Response> FilesProcess::_browse(
    const http::Request& request,
    const Option<string>& principal)
{
  Option<string> path = request.url.query.get("path");

  if (!path.isSome() || path.get().empty()) {
    return BadRequest("Expecting 'path=value' in query.\n");
  }

  string requestedPath = path.get();
  Option<string> jsonp = request.url.query.get("jsonp");

  return browse(requestedPath, principal)
  .then([jsonp](const Try<list<FileInfo>, FilesError>& result)
    -> Future<http::Response> {
    if (result.isError()) {
      const FilesError& error = result.error();

      switch (error.type) {
        case FilesError::Type::INVALID:
          return BadRequest(error.message);

        case FilesError::Type::NOT_FOUND:
          return NotFound(error.message);

        case FilesError::Type::UNAUTHORIZED:
          return Forbidden(error.message);

        case FilesError::Type::UNKNOWN:
          return InternalServerError(error.message);
      }

      UNREACHABLE();
    }

    JSON::Array listing;
    foreach (const FileInfo& fileInfo, result.get()) {
      listing.values.push_back(model(fileInfo));
    }

    return OK(listing, jsonp);
  });
}

Future<Try<list<FileInfo>, FilesError>> FilesProcess::browse(
    const string& path,
    const Option<string>& principal)
{
  return authorize(path, principal)
    .then(defer(self(),
        [this, path](bool authorized)
          -> Future<Try<list<FileInfo>, FilesError>> {
      if (!authorized) {
        return FilesError(FilesError::Type::UNAUTHORIZED);
      }

      Result<string> resolvedPath = resolve(path);

      if (resolvedPath.isError()) {
        return FilesError(
            FilesError::Type::INVALID,
            resolvedPath.error() + ".\n");
      } else if (resolvedPath.isNone()) {
        return FilesError(FilesError::Type::NOT_FOUND);
      }

      // The result will be a sorted (on path) list of files and dirs.
      map<string, FileInfo> files;
      Try<list<string> > entries = os::ls(resolvedPath.get());
      if (entries.isSome()) {
        foreach (const string& entry, entries.get()) {
          struct stat s;
          string fullPath = path::join(resolvedPath.get(), entry);

          if (stat(fullPath.c_str(), &s) < 0) {
            PLOG(WARNING) << "Found " << fullPath << " in ls but stat failed";
            continue;
          }

          files[fullPath] =
            protobuf::createFileInfo(path::join(path, entry), s);
        }
      }

      list<FileInfo> listing;
      foreachvalue (const FileInfo& fileInfo, files) {
        listing.push_back(fileInfo);
      }

      return listing;
    }));
}

// TODO(benh): Remove 'const &' from size after fixing libprocess.
Future<http::Response> __read(int fd,
                       const size_t& size,
                       off_t offset,
                       const boost::shared_array<char>& data,
                       const Option<string>& jsonp) {
  JSON::Object object;

  object.values["offset"] = offset;
  object.values["data"] = string(data.get(), size);

  return OK(object, jsonp);
}

const string FilesProcess::READ_HELP = HELP(
    TLDR(
        "Reads data from a file."),
    DESCRIPTION(
        "This endpoint reads data from a file at a given offset and for",
        "a given length."
        "",
        "Query parameters:",
        "",
        ">        path=VALUE          The path of directory to browse.",
        ">        offset=VALUE        Value added to base address to obtain "
        "a second address",
        ">        length=VALUE        Length of file to read."),
    AUTHENTICATION(true),
    AUTHORIZATION(
        "Reading files requires that the request principal is",
        "authorized to do so for the target virtual file path.",
        "",
        "Authorizers may categorize different virtual paths into",
        "different ACLs, e.g. logs in one and task sandboxes in",
        "another.",
        "",
        "See authorization documentation for details."));

Future<http::Response> FilesProcess::read(
    const http::Request& request,
    const Option<string>& principal)
{
  Option<string> path = request.url.query.get("path");

  if (!path.isSome() || path.get().empty()) {
    return BadRequest("Expecting 'path=value' in query.\n");
  }

  off_t offset = -1;

  if (request.url.query.get("offset").isSome()) {
    Try<off_t> result = numify<off_t>(request.url.query.get("offset").get());

    if (result.isError()) {
      return BadRequest("Failed to parse offset: " + result.error() + ".\n");
    }

    offset = result.get();
  }

  Option<size_t> length;

  if (request.url.query.get("length").isSome()) {
    Try<ssize_t> result = numify<ssize_t>(
        request.url.query.get("length").get());

    if (result.isError()) {
      return BadRequest("Failed to parse length: " + result.error() + ".\n");
    }

    // TODO(tomxing): The pailer in the webui sends `length=-1` at first to
    // determine the length of the file, so we allow a length of -1.
    // Setting `length=-1` has the same effect as not providing a length: we
    // read to the end of the file, up to the maximum read length.
    // Will change this logic in MESOS-5334.
    if (result.get() < -1) {
      return BadRequest(
        strings::format("Negative length provided: %d.\n", result.get()).get());
    }

    if (result.get() > -1){
      length = result.get();
    }
  }

  string requestedPath = path.get();
  Option<string> jsonp = request.url.query.get("jsonp");

  return authorize(requestedPath, principal)
    .then(defer(self(),
        [this, offset, length, path, jsonp](bool authorized)
          -> Future<http::Response> {
      if (authorized) {
        return _read(offset, length, path.get(), jsonp);
      }

      return Forbidden();
    }));
}

Future<http::Response> FilesProcess::_read(
    off_t offset,
    Option<size_t> length,
    const string& path,
    const Option<string>& jsonp)
{
  Result<string> resolvedPath = resolve(path);

  if (resolvedPath.isError()) {
    return BadRequest(resolvedPath.error() + ".\n");
  } else if (!resolvedPath.isSome()) {
    return NotFound();
  }

  // Don't read directories.
  if (os::stat::isdir(resolvedPath.get())) {
    return BadRequest("Cannot read a directory.\n");
  }

  // TODO(benh): Cache file descriptors so we aren't constantly
  // opening them and paging the data in from disk.
  Try<int> fd = os::open(resolvedPath.get(), O_RDONLY | O_CLOEXEC);

  if (fd.isError()) {
    string error = strings::format(
        "Failed to open file at '%s': %s",
        resolvedPath.get(),
        fd.error()).get();
    LOG(WARNING) << error;
    return InternalServerError(error + ".\n");
  }

  off_t size = lseek(fd.get(), 0, SEEK_END);

  if (size == -1) {
    string error = strings::format(
        "Failed to open file at '%s': %s",
        resolvedPath.get(),
        os::strerror(errno)).get();

    LOG(WARNING) << error;
    os::close(fd.get());
    return InternalServerError(error + ".\n");
  }

  if (offset == -1) {
    offset = size;
  }

  if (length.isNone()) {
    length = size - offset;
  }

  // Cap the read length at 16 pages.
  length = std::min<size_t>(length.get(), os::pagesize() * 16);

  if (offset >= size) {
    os::close(fd.get());

    JSON::Object object;
    object.values["offset"] = size;
    object.values["data"] = "";
    return OK(object, jsonp);
  }

  // Seek to the offset we want to read from.
  if (lseek(fd.get(), offset, SEEK_SET) == -1) {
    string error = strings::format(
        "Failed to seek file at '%s': %s",
        resolvedPath.get(),
        os::strerror(errno)).get();

    LOG(WARNING) << error;
    os::close(fd.get());
    return InternalServerError(error);
  }

  Try<Nothing> nonblock = os::nonblock(fd.get());
  if (nonblock.isError()) {
    string error =
        "Failed to set file descriptor nonblocking: " + nonblock.error();
    LOG(WARNING) << error;
    os::close(fd.get());
    return InternalServerError(error);
  }

  // Read 'length' bytes (or to EOF).
  boost::shared_array<char> data(new char[length.get()]);

  return io::read(fd.get(), data.get(), length.get())
    .then(lambda::bind(
        __read,
        fd.get(),
        lambda::_1,
        offset,
        data,
        jsonp))
    .onAny([fd]() { os::close(fd.get()); });
}

const string FilesProcess::DOWNLOAD_HELP = HELP(
    TLDR(
        "Returns the raw file contents for a given path."),
    DESCRIPTION(
        "This endpoint will return the raw file contents for the",
        "given path.",
        "",
        "Query parameters:",
        "",
        ">        path=VALUE          The path of directory to browse."),
    AUTHENTICATION(true),
    AUTHORIZATION(
        "Downloading files requires that the request principal is",
        "authorized to do so for the target virtual file path.",
        "",
        "Authorizers may categorize different virtual paths into",
        "different ACLs, e.g. logs in one and task sandboxes in",
        "another.",
        "",
        "See authorization documentation for details."));

Future<http::Response> FilesProcess::download(
    const http::Request& request,
    const Option<string>& principal)
{
  Option<string> path = request.url.query.get("path");

  if (!path.isSome() || path.get().empty()) {
    return BadRequest("Expecting 'path=value' in query.\n");
  }

  string requestedPath = path.get();

  return authorize(requestedPath, principal)
    .then(defer(self(),
        [this, path](bool authorized) -> Future<http::Response> {
      if (authorized) {
        return _download(path.get());
      }

      return Forbidden();
    }));
}

Future<http::Response> FilesProcess::_download(const string& path)
{
  Result<string> resolvedPath = resolve(path);

  if (resolvedPath.isError()) {
    return BadRequest(resolvedPath.error() + ".\n");
  } else if (!resolvedPath.isSome()) {
    return NotFound();
  }

  // Don't download directories.
  if (os::stat::isdir(resolvedPath.get())) {
    return BadRequest("Cannot download a directory.\n");
  }

  string basename = Path(resolvedPath.get()).basename();

  OK response;
  response.type = response.PATH;
  response.path = resolvedPath.get();
  response.headers["Content-Type"] = "application/octet-stream";
  response.headers["Content-Disposition"] =
    strings::format("attachment; filename=%s", basename).get();

  // Attempt to detect the mime type.
  Option<string> extension = Path(resolvedPath.get()).extension();

  if (extension.isSome() && mime::types.count(extension.get()) > 0) {
    response.headers["Content-Type"] = mime::types[extension.get()];
  }

  return response;
}

const string FilesProcess::DEBUG_HELP = HELP(
    TLDR(
        "Returns the internal virtual path mapping."),
    DESCRIPTION(
        "This endpoint shows the internal virtual path map as a",
        "JSON object."),
    AUTHENTICATION(true),
    AUTHORIZATION(
        "The request principal should be authorized to query this endpoint.",
        "See the authorization documentation for details."));

Future<http::Response> FilesProcess::debug(
    const http::Request& request,
    const Option<string>&  principal )
{
  JSON::Object object;
  foreachpair (const string& name, const string& path, paths) {
    object.values[name] = path;
  }

  const Option<string>& jsonp = request.url.query.get("jsonp");

  return authorizeEndpoint(
      request.url.path,
      request.method,
      authorizer,
      principal)
    .then(defer([object, jsonp](bool authorized)
                  -> Future<http::Response> {
          if (!authorized) {
            return Forbidden();
          }

          return OK(object, jsonp);
        }));
}

Result<string> FilesProcess::resolve(const string& path)
{
  // Suppose we have: /1/2/hello_world.txt
  // And we attach: /1/2 as /sandbox
  // Then this function would resolve the following virtual path
  // into the actual path:
  // input: /sandbox/hello_world.txt
  // output: /1/2/hello_world.txt
  //
  // Try and see if this path has been attached. We check for the
  // longest possible prefix match and if found append any suffix to
  // the attached path (provided the path is to a directory).
  vector<string> tokens = strings::split(
      strings::remove(path, stringify(os::PATH_SEPARATOR), strings::SUFFIX),
      stringify(os::PATH_SEPARATOR));

  string suffix;
  while (!tokens.empty()) {
    string prefix = path::join(tokens);

    if (!paths.contains(prefix)) {
      if (suffix.empty()) {
        suffix = tokens.back();
      } else {
        suffix = path::join(tokens.back(), suffix);
      }

      tokens.pop_back();
      continue;
    }

    // Determine the final path: if it's a directory, append the
    // suffix, if it's not a directory and there is a suffix, return
    // 'Not Found'.
    string path = paths[prefix];
    if (os::stat::isdir(path)) {
      path = path::join(path, suffix, os::PATH_SEPARATOR);

      // Canonicalize the absolute path.
      Result<string> realpath = os::realpath(path);
      if (realpath.isError()) {
        return Error(
            "Failed to determine canonical path of '" + path +
            "': " + realpath.error());
      } else if (realpath.isNone()) {
        return None();
      }

      // Make sure the canonicalized absolute path is accessible
      // (i.e., not outside the "chroot").
      if (!strings::startsWith(realpath.get(), paths[prefix])) {
        return Error("'" + path + "' is inaccessible");
      }

      path = realpath.get();
    } else if (suffix != "") {
      // Request is assuming attached path is a directory, but it is
      // not! Rather than 'Bad Request', treat this as 'Not Found'.
      return None();
    }

    return path;
  }

  return None();
}

Files::Files(const Option<string>& authenticationRealm,
             const Option<Authorizer*>& authorizer)
{
  process = new FilesProcess(authenticationRealm, authorizer);
  spawn(process);
}

Files::~Files()
{
  terminate(process);
  wait(process);
  delete process;
}

Future<Nothing> Files::attach(
    const string& path,
    const string& name,
    const Option<lambda::function<Future<bool>(const Option<string>&)>>&
        authorized)
{
  return dispatch(process, &FilesProcess::attach, path, name, authorized);
}

void Files::detach(const string& name)
{
  dispatch(process, &FilesProcess::detach, name);
}

Future<Try<list<FileInfo>, FilesError>> Files::browse(
    const string& path,
    const Option<string>& principal)
{
  return dispatch(process, &FilesProcess::browse, path, principal);
}

} // namespace internal {
} // namespace mesos {

ext: cpp


text:
#include "fatwater.h"

#include "hoMatrix.h"
#include "hoNDArray_linalg.h"
#include "hoNDArray_utils.h"
#include "hoNDArray_elemwise.h"
#include "hoNDArray_reductions.h"
#include "hoArmadillo.h"

#include <boost/config.hpp>

#include <boost/graph/graph_utility.hpp>
#include <boost/graph/boykov_kolmogorov_max_flow.hpp>

#include <boost/timer/timer.hpp>
#include <boost/iterator/function_input_iterator.hpp>
#include <iterator>
#include "graph_cut.h"

#include "hoNDArray_math.h"

#include <boost/random.hpp>
#include <boost/math/constants/constants.hpp>
#include <armadillo>
#include <random>
#include <cpu/hoNDArray_fileio.h>
#include <GadgetronTimer.h>
#include <complex>
#include <cpu/math/hoNDImage_util.h>

#include <algorithm>
#include "bounded_field_map.h"

using namespace boost;

namespace Gadgetron {
    namespace FatWater {
        using namespace std::complex_literals;
        static constexpr float PI = boost::math::constants::pi<float>();
        static std::mt19937 rng_state(4242);

        hoNDArray<std::vector<uint16_t>> find_local_minima(const hoNDArray<float> &residuals, float threshold = 0.00f) {

            auto threshold_signal = sum(residuals, 0);
            threshold_signal /= max(threshold_signal);
            sqrt_inplace(&threshold_signal);

            auto min_residuals = min(residuals, 0);
            auto max_residuals = max(residuals, 0);

            const auto Y = residuals.get_size(2);
            const auto X = residuals.get_size(1);
            hoNDArray<std::vector<uint16_t>> result(X, Y);
            const auto steps = residuals.get_size(0);
            size_t count = 0;
            for (size_t k2 = 0; k2 < Y; k2++) {
                for (size_t k1 = 0; k1 < X; k1++) {

                    std::vector<uint16_t> minima;
//                    if (threshold_signal(k1, k2) > threshold) {
                    for (size_t k0 = 1; k0 < steps - 1; k0++) {
                        if ((residuals(k0, k1, k2) < residuals(k0 - 1, k1, k2)) &&
                            (residuals(k0 + 1, k1, k2) >= residuals(k0, k1, k2))) {
                            minima.push_back(k0);
                        }

                    }
//                    }
                    if (minima.size() < 2) count++;
                    result(k1, k2) = std::move(minima);
                }
            }

            GDEBUG("COUNT DRACULA %i\n", count);
            return result;
        }

        hoNDArray<float> approx_second_derivative(const hoNDArray<float> &residuals,
                                                  const hoNDArray<std::vector<uint16_t>> &local_min_indices,
                                                  float step_size) {
            hoNDArray<float> second_deriv(local_min_indices.dimensions());

            const auto Y = second_deriv.get_size(1);
            const auto X = second_deriv.get_size(0);
            const auto nfields = residuals.get_size(0);

            for (uint16_t k2 = 0; k2 < Y; k2++) {
                for (uint16_t k1 = 0; k1 < X; k1++) {

                    auto minimum = std::min_element(&residuals(1, k1, k2), &residuals(nfields - 1, k1, k2)) -
                                   &residuals(0, k1, k2);

                    auto sd =
                            (residuals(minimum - 1, k1, k2) + residuals(minimum + 1, k1, k2) -
                             2 * residuals(minimum, k1, k2)) / (step_size * step_size);

                    second_deriv(k1, k2) = sd;

                }
            }

            return second_deriv;

        }

        hoNDArray<float>
        create_field_map(const hoNDArray<uint16_t> &field_map_index, const std::vector<float> &field_map_strengths) {
            const uint16_t max_val = field_map_strengths.size() - 1;
            hoNDArray<float> field_map(field_map_index.dimensions());
            std::transform(field_map_index.begin(), field_map_index.end(), field_map.begin(),
                           [&](uint16_t i) { return field_map_strengths[std::min(i, max_val)]; });
            return field_map;
        }

        hoNDArray<uint16_t> create_field_map_proposal1(const hoNDArray<uint16_t> &field_map_index,
                                                       const hoNDArray<std::vector<uint16_t>> &minima,
                                                       const hoNDArray<float> &residuals,
                                                       const std::vector<float> &field_map_strengths, float fat_freq,
                                                       float dF, float dTE) {

            const size_t elements = field_map_index.get_number_of_elements();
            hoNDArray<uint16_t> proposed_field_map_index(field_map_index.dimensions());
            const size_t field_maps = field_map_strengths.size();
            std::uniform_int_distribution<int> coinflip(0, 1);
            int jump;
            if (coinflip(rng_state)) {
                jump = round(std::abs(fat_freq / dF));

            } else {
                jump = round((1.0 / dTE - std::abs(fat_freq)) / dF);
            }

            for (size_t i = 0; i < elements; i++) {
                auto &mins = minima[i];
                auto fbi = field_map_index[i];
                auto fqmi = std::find_if(mins.begin(), mins.end(),
                                         [&](auto fqi) { return fqi > fbi + 20; }); //Find smallest
                proposed_field_map_index[i] = (fqmi == mins.end()) ? std::min<int>(fbi + jump, field_maps - 1) : *fqmi;
            }

            return proposed_field_map_index;

        }

        hoNDArray<uint16_t> create_field_map_proposal2(const hoNDArray<uint16_t> &field_map_index,
                                                       const hoNDArray<std::vector<uint16_t>> &minima,
                                                       const hoNDArray<float> &residuals,
                                                       const std::vector<float> &field_map_strengths, float fat_freq,
                                                       float dF, float dTE) {

            const size_t elements = field_map_index.get_number_of_elements();
            hoNDArray<uint16_t> proposed_field_map_index(field_map_index.dimensions());
            std::uniform_int_distribution<int> coinflip(0, 1);
            int jump;
            if (coinflip(rng_state)) {
                jump = round(std::abs(fat_freq / dF));
            } else {
                jump = round((1.0 / dTE - std::abs(fat_freq)) / dF);
            }

            for (size_t i = 0; i < elements; i++) {
                auto &mins = minima[i];
                int fbi = field_map_index[i];
                auto fqmi = std::find_if(mins.rbegin(), mins.rend(),
                                         [&](auto fqi) { return fqi < (fbi - 20); }); //Find smallest
                proposed_field_map_index[i] = (fqmi == mins.rend()) ? std::max<int>(fbi - jump, 0) : *fqmi;
            }

            return proposed_field_map_index;

        }

        hoNDArray<uint16_t>
        create_field_map_proposal_standard(const hoNDArray<uint16_t> &field_map_index, int sign,
                                           uint16_t max_field_value) {

            std::uniform_int_distribution<int> rng(1, 3);

            int step_size = sign * rng(rng_state);
//        int step_size = sign;
            hoNDArray<uint16_t> proposed_field_map_index(field_map_index.dimensions());
            std::transform(field_map_index.begin(), field_map_index.end(), proposed_field_map_index.begin(),
                           [&](uint16_t j) {
                               return uint16_t(std::min(std::max(j + step_size, 0), int(max_field_value)));
                           });

            return proposed_field_map_index;
        }

        hoNDArray<uint16_t>
        solve_MRF(const Config &config, const std::vector<float> &field_map_strengths,
                  const hoNDArray<float> &residual, const hoNDArray<std::vector<uint16_t>> &local_min_indices,
                  const hoNDArray<float> &second_deriv, float fat_freq, float dF, float dTE) {

            hoNDArray<uint16_t> fmIndex(local_min_indices.dimensions());

            std::uniform_int_distribution<int> coinflip(0, 2);
            fmIndex.fill(field_map_strengths.size() / 2);

            hoNDArray<uint16_t> fmIndex_update;
            for (int i = 0; i < config.number_of_iterations; i++) {
                if (coinflip(rng_state) == 0 || i < 15) {
                    if (!(i % 2)) {
                        fmIndex_update = create_field_map_proposal1(fmIndex, local_min_indices, residual,
                                                                    field_map_strengths, fat_freq, dF, dTE);
                    } else {
                        fmIndex_update = create_field_map_proposal2(fmIndex, local_min_indices, residual,
                                                                    field_map_strengths, fat_freq, dF, dTE);
                    }
                } else {
                    fmIndex_update = create_field_map_proposal_standard(fmIndex, std::pow(-1, i),
                                                                        field_map_strengths.size() - 1);
                }

                fmIndex = update_field_map(fmIndex, fmIndex_update, residual, second_deriv);
            }

            return fmIndex;
        }

        arma::Mat<std::complex<float>>
        calculate_psi_matrix(const std::vector<float> &echoTimes,
                             const arma::Mat<std::complex<float>> &phiMatrix,  float r2star) {
            arma::Mat<std::complex<float>> psiMatrix(phiMatrix.n_rows, phiMatrix.n_cols);
            for (int k1 = 0; k1 < phiMatrix.n_rows; k1++) {
                auto curModulation = exp(-r2star * (echoTimes[k1] - echoTimes[0]));
                for (int k2 = 0; k2 < phiMatrix.n_cols; k2++) {
                    psiMatrix(k1, k2) = phiMatrix(k1, k2) * curModulation;
                }
            }

            return psiMatrix;
        }

        hoNDArray<arma::Mat<std::complex<float>>>
        calculate_projection_matrices(const std::vector<float> &echo_times,
                                      const arma::Mat<std::complex<float>> &phiMatrix,
                                      const std::vector<float> &field_map_strengths,
                                      const std::vector<float> &r2stars) {

            auto num_fm = field_map_strengths.size();
            auto num_r2star = r2stars.size();
            hoNDArray<arma::Mat<std::complex<float>>> Ps(num_fm, num_r2star);
            size_t nte = echo_times.size();

#ifdef WIN32
#pragma omp parallel for
#else
#pragma omp parallel for
#endif
            for (int k3 = 0; k3 < num_fm; k3++) {

                float fm = field_map_strengths[k3];
                arma::Col<std::complex<float>> b_shifts(nte);
                for (int kt = 0; kt < nte; kt++)
                    b_shifts[kt] = std::exp(2if * PI * (echo_times[kt] - echo_times[0]) * fm);
                for (int k4 = 0; k4 < num_r2star; k4++) {
                    float r2star = r2stars[k4];

                    arma::Mat<std::complex<float>> psiMatrix = calculate_psi_matrix(echo_times, phiMatrix,
                                                                                r2star);
                    Ps(k3, k4) = arma::diagmat(b_shifts) * (arma::eye<arma::Mat<std::complex<float>>>(nte, nte) -
                                                            psiMatrix * arma::pinv(psiMatrix)) *
                                 arma::diagmat(arma::conj(b_shifts));

                }
            }
            return Ps;
        }

        hoNDArray<float>
        calculate_r2star_map(const hoNDArray<std::complex<float> > &data, const hoNDArray<float> &field_map,
                             const std::vector<float> &r2star_values,
                             const arma::Mat<std::complex<float>> &phiMatrix, const std::vector<float> &echoTimes) {
            using cMat = arma::Mat<std::complex<float>>;
            uint16_t X = data.get_size(0);
            uint16_t Y = data.get_size(1);
            uint16_t Z = data.get_size(2);
            uint16_t CHA = data.get_size(3);
            uint16_t N = data.get_size(4);
            uint16_t S = data.get_size(5);
            uint16_t LOC = data.get_size(6);
            auto nte = phiMatrix.n_rows;

            auto data_corrected = data;

            for (int kS = 0; kS < S; kS++) {
                for (int kN = 0; kN < N; kN++) {
                    for (int kcha = 0; kcha < CHA; kcha++) {
                        for (int ky = 0; ky < Y; ky++) {
                            for (int kx = 0; kx < X; kx++) {
                                data_corrected(kx, ky, 0, kcha, kN, kS) *= std::exp(
                                        -2if * PI * field_map(kx, ky) * echoTimes[kS]);
                            }
                        }
                    }
                }
            }

            std::vector<arma::Mat<std::complex<float>>> Ps;
            for (auto r2star : r2star_values) {

                auto psiMatrix = calculate_psi_matrix(echoTimes, phiMatrix,
                                                      r2star);
                Ps.emplace_back(
                        arma::eye<arma::Mat<std::complex<float>>>(nte, nte) - psiMatrix * arma::pinv(psiMatrix)
                );

            }

            hoNDArray<float> r2star_map(field_map.dimensions());
#ifdef WIN32
#pragma omp parallel for
#else
#pragma omp parallel for collapse(2)
#endif
            for (int k2 = 0; k2 < Y; k2++) {
                for (int k1 = 0; k1 < X; k1++) {
                    // Get current signal
                    std::vector<cMat> signals(CHA, cMat(S, N));
                    for (int cha = 0; cha < CHA; cha++) {
                        auto &tempSignal = signals[cha];
                        for (int k4 = 0; k4 < N; k4++) {
                            for (int k5 = 0; k5 < S; k5++) {
                                tempSignal(k5, k4) = data_corrected(k1, k2, 0, cha, k4, k5, 0);

                            }
                        }
                    }

                    float minResidual = std::numeric_limits<float>::max();
                    std::vector<float> field_map_strengths = {field_map(k1, k2)};

                    for (int kr2 = 0; kr2 < r2star_values.size(); kr2++) {

                        float curResidual = 0;

                        for (int cha = 0; cha < CHA; cha++) {
                            // Apply projector
                            arma::Mat<std::complex<float>> projected = Ps[kr2] * signals[cha];
                            curResidual += std::accumulate(projected.begin(), projected.end(), 0.0f,
                                                           [](auto v1, auto v2) {
                                                               return v1 +
                                                                      std::norm(v2);
                                                           });
                        }
                        if (curResidual < minResidual) {
                            minResidual = curResidual;
                            r2star_map(k1, k2) = r2star_values[kr2];
                        }

                    }

                }
            }

            return r2star_map;
        }

        std::tuple<hoNDArray<float>, hoNDArray<uint16_t>>
        calculate_residual_and_r2star(const hoNDArray<std::complex<float>> &data, const Parameters &parameters,
                                      const arma::Mat<std::complex<float>> &phi,
                                      const std::vector<float> &field_strengths,
                                      const std::vector<float> &r2star_values) {
            using cMat = arma::Mat<std::complex<float>>;
            uint16_t X = data.get_size(0);
            uint16_t Y = data.get_size(1);
            uint16_t Z = data.get_size(2);
            uint16_t CHA = data.get_size(3);
            uint16_t N = data.get_size(4);
            uint16_t S = data.get_size(5);
            uint16_t LOC = data.get_size(6);

            auto projection_matrices = calculate_projection_matrices(parameters.echo_times_s, phi, field_strengths,
                                                                     r2star_values);

            auto result = std::make_tuple(hoNDArray<float>(field_strengths.size(), X, Y, Z),
                                          hoNDArray<uint16_t>(X, Y, Z, field_strengths.size()));

            auto &residual = std::get<0>(result);
            auto &r2starIndex = std::get<1>(result);

#ifdef WIN32
#pragma omp parallel for
#else
#pragma omp parallel for collapse(3)
#endif
            for (int kz = 0; kz < Z; kz++) {
                for (int ky = 0; ky < Y; ky++) {
                    for (int kx = 0; kx < X; kx++) {

                        std::vector<cMat> signals(CHA, cMat(S, N));
                        for (int cha = 0; cha < CHA; cha++) {
                            auto &tempSignal = signals[cha];
                            for (int kn = 0; kn < N; kn++) {
                                for (int ks = 0; ks < S; ks++) {
                                    tempSignal(ks, kn) = data(kx, ky, kz, cha, kn, ks, 0);
                                }
                            }
                        }

                        for (int kf = 0; kf < field_strengths.size(); kf++) {

                            float minResidual = std::numeric_limits<float>::max();

                            for (int kr = 0; kr < r2star_values.size(); kr++) {
                                // Apply projector
                                float curResidual = 0;
                                for (int cha = 0; cha < CHA; cha++) {
                                    arma::Mat<std::complex<float>> projected =
                                            projection_matrices(kf, kr) * signals[cha];
                                    curResidual += std::accumulate(projected.begin(), projected.end(), 0.0f,
                                                                   [](auto v1, auto v2) {
                                                                       return v1 +
                                                                              std::norm(v2);
                                                                   });
                                }
                                if (curResidual < minResidual) {
                                    minResidual = curResidual;
                                    r2starIndex(kx, ky, kz, kf) = kr;
                                }
                            }
                            residual(kf, kx, ky, kz) = minResidual;

                        }
                    }
                }
            }

            return result;
        }

        hoNDArray<std::complex<float>>
        separate_species(const hoNDArray<std::complex<float>> &data,
                         const arma::Mat<std::complex<float>> &phiMatrix, hoNDArray<float> &r2star_map,
                         hoNDArray<float> &field_map, const Parameters &parameters) {

            using cMat = arma::Mat<std::complex<float>>;
            uint16_t X = data.get_size(0);
            uint16_t Y = data.get_size(1);
            uint16_t Z = data.get_size(2);
            uint16_t CHA = data.get_size(3);
            uint16_t N = data.get_size(4);
            uint16_t S = data.get_size(5);
            uint16_t LOC = data.get_size(6);
            hoNDArray<std::complex<float> > out(X, Y, Z, CHA, N, parameters.species.size(),
                                                LOC); // S dimension gets replaced by water/fat stuff

#ifdef WIN32
#pragma omp parallel for
#else
#pragma omp parallel for collapse(3)
#endif
            for (int kz = 0; kz < Z; kz++) {
                for (int ky = 0; ky < Y; ky++) {
                    for (int kx = 0; kx < X; kx++) {
                        std::vector<cMat> signals(CHA, cMat(S, N));

                        // Get current signal
                        for (int cha = 0; cha < CHA; cha++) {
                            auto &tempSignal = signals[cha];
                            for (int kn = 0; kn < N; kn++) {
                                for (int ks = 0; ks < S; ks++) {
                                    tempSignal(ks, kn) = data(kx, ky, kz, cha, kn, ks, 0);
                                }
                            }
                        }

                        auto fm = field_map(kx, ky);
                        auto r2star = r2star_map(kx, ky);

                        arma::Mat<std::complex<float>> psiMatrix = calculate_psi_matrix(parameters.echo_times_s,
                                                                                phiMatrix,r2star);

                        auto nte = parameters.echo_times_s.size();
                        arma::Col<std::complex<float>> b_shifts(nte);
                        for (int kt = 0; kt < nte; kt++)
                            b_shifts[kt] = std::exp(
                                    2if * PI * (parameters.echo_times_s[kt] - parameters.echo_times_s[0]) * fm);
                        psiMatrix = arma::diagmat(b_shifts) * psiMatrix;
                        // Solve for water and fat

                        for (int cha = 0; cha < CHA; cha++) {
                            arma::Mat<std::complex<float>> curWaterFat = arma::solve(psiMatrix, signals[cha]);
                            for (int kn = 0; kn < N; kn++) {
                                for (int kspecies = 0; kspecies <
                                                       parameters.species.size(); kspecies++) { // 2 elements for water and fat currently
                                    out(kx, ky, kz, cha, kn, kspecies, 0) = curWaterFat(kspecies, kn);
                                }
                            }
                        }

                    }
                }
            }
            return out;
        }

        arma::Mat<std::complex<float>>
        calculatePhiMatrix(const Parameters parameters) {

            auto echoTimes = parameters.echo_times_s;

            auto nte = echoTimes.size();
            auto nspecies = parameters.species.size();
            typedef arma::Mat<std::complex<float>> Cmat;
            Cmat phiMatrix = arma::zeros<Cmat>(nte, nspecies);
            for (int k1 = 0; k1 < nte; k1++) {
                for (int k2 = 0; k2 < nspecies; k2++) {
                    auto &species = parameters.species[k2];
                    auto npeaks = species.amplitude_frequency_pairs.size();
                    for (int k3 = 0; k3 < npeaks; k3++) {
                        auto relAmp = species.amplitude_frequency_pairs[k3].first;
                        auto freq_hz = parameters.field_strength_T * parameters.gyromagnetic_ratio_Mhz *
                                       species.amplitude_frequency_pairs[k3].second;
                        phiMatrix(k1, k2) += relAmp * exp(2if * PI * echoTimes[k1] * freq_hz);

                    }

                }
            }
            return phiMatrix;
        }

        std::vector<float> linspace(const std::pair<float, float> &range, unsigned int count) {
            std::vector<float> result(count);
            result[0] = range.first;
            for (int i = 1; i < count; i++) {
                result[i] = range.first + i * (range.second - range.first) / (count - 1);
            }
            return result;
        }

        std::complex<float> mean_frequency(const Parameters &parameters) {
            ChemicalSpecies fat = parameters.species[1];
            auto average_fat_freq =
                    accumulate(fat.amplitude_frequency_pairs.begin(), fat.amplitude_frequency_pairs.end(), std::complex<float>(0.0),
                               [](auto val, auto tup) {
                                   return val + std::get<0>(tup) * std::get<1>(tup);
                               }) /
                    accumulate(fat.amplitude_frequency_pairs.begin(), fat.amplitude_frequency_pairs.end(), std::complex<float>(0.0),
                               [](auto val, auto tup) { return val + std::get<0>(tup); });
            average_fat_freq *= parameters.field_strength_T * parameters.gyromagnetic_ratio_Mhz;
            return average_fat_freq;
        }

        std::tuple<hoNDArray<float>, hoNDArray<float>>
        calculate_field_map(const hoNDArray<std::complex<float>> &data, const Parameters &parameters,
                            const Config &config,
                            const std::complex<float> &average_fat_freq, const arma::Mat<std::complex<float>> &phi,
                            const std::vector<float> &field_map_strengths);

        FatWater::Output
        fatwater_separation(const hoNDArray<std::complex<float> > &data, Parameters parameters,
                            Config config) {

            uint16_t X = data.get_size(0);
            uint16_t Y = data.get_size(1);
            uint16_t Z = data.get_size(2);
            uint16_t CHA = data.get_size(3);
            uint16_t N = data.get_size(4);
            uint16_t S = data.get_size(5);
            uint16_t LOC = data.get_size(6);

            auto echo_times = parameters.echo_times_s;
            GDEBUG("In toolbox - Field Strength: %f T \n", parameters.field_strength_T);
            for (auto &te: echo_times) {
                GDEBUG("In toolbox - Echo time: %f seconds \n", te);
            }
            GDEBUG("In toolbox - PrecessionIsClockwise: %d \n", parameters.precession_is_clockwise);

            //Get or set some algorithm parameters

            std::complex<float> average_fat_freq = mean_frequency(parameters);

            arma::Mat<std::complex<float>> phi = calculatePhiMatrix(parameters);

            std::vector<float> field_map_strengths = linspace(config.frequency_range,
                                                              config.number_of_frequency_samples);

            hoNDArray<float> field_map, r2star_map;
            std::tie(field_map, r2star_map) = calculate_field_map(data, parameters, config, average_fat_freq, phi,
                                                                  field_map_strengths);

            auto species = separate_species(data, phi, r2star_map, field_map, parameters);

            if (config.do_gradient_descent) {
                bounded_field_map(field_map, data, parameters, (field_map_strengths[1] - field_map_strengths[0])*2);
                species = separate_species(data, phi, r2star_map, field_map, parameters);
            }

            return Output{std::move(species), std::move(field_map), std::move(r2star_map)};
        }

        std::tuple<hoNDArray<float>, hoNDArray<float>>
        calculate_field_map(const hoNDArray<std::complex<float>> &data, const Parameters &parameters,
                            const Config &config,
                            const std::complex<float> &average_fat_freq, const arma::Mat<std::complex<float>> &phi,
                            const std::vector<float> &field_map_strengths) {

            std::vector<float> r2stars_fine = linspace(config.r2_range, config.number_of_r2_fine_samples);
            std::vector<float> r2stars = linspace(config.r2_range, config.number_of_r2_samples);

            auto data_scaled = data;
            for (int downsamples = 0; downsamples < config.downsamples; downsamples++)
                data_scaled = downsample<std::complex<float>, 2>(data_scaled);

            hoNDArray<float> residual;
            hoNDArray<uint16_t> r2starIndex;
            std::tie(residual, r2starIndex) = calculate_residual_and_r2star(data_scaled, parameters, phi,
                                                                            field_map_strengths,
                                                                            r2stars);

            auto dF = field_map_strengths[1] - field_map_strengths[0];

            hoNDArray<std::vector<uint16_t>> local_min_indices = find_local_minima(residual);

            hoNDArray<float> lambda_map = approx_second_derivative(residual, local_min_indices, dF);
            lambda_map += mean(&lambda_map) * config.lambda_extra;
            lambda_map *= config.lambda * dF * dF;

            hoNDArray<uint16_t> fmIndex = solve_MRF(config, field_map_strengths, residual, local_min_indices,
                                                    lambda_map, abs(average_fat_freq), dF,
                                                    parameters.echo_times_s[1] - parameters.echo_times_s[0]);

            hoNDArray<float> field_map = create_field_map(fmIndex, field_map_strengths);

            if (config.downsamples) {
                field_map = upsample_spline<float, 2>(field_map,std::pow(2,config.downsamples));
            }
//            fmIndex = upsample<uint16_t,2>(&fmIndex);

            auto r2star_map = calculate_r2star_map(data, field_map, r2stars_fine, phi,
                                                   parameters.echo_times_s);

//            r2star_map = upsample<float,2>(&r2star_map);
            return std::make_tuple(std::move(field_map), std::move(r2star_map));
        }

    }
}

max_stars_repo_name: roopchansinghv/gadgetron