Commit ce6381d7125584b110ee8cd0a22ffb3f78d4d88a

Authored by David Mayerich
1 parent 35b5b79c

updating to TIRA

Showing 197 changed files with 108965 additions and 0 deletions   Show diff stats

Too many changes.

To preserve performance only 100 of 197 files are displayed.

tira/biology/fibernet.h 0 → 100644
tira/biomodels/cellset.h 0 → 100644
  1 +/*
  2 +Copyright <2017> <David Mayerich>
  3 +
  4 +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
  5 +
  6 +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
  7 +
  8 +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  9 +*/
  10 +
  11 +#ifndef STIM_CELLSET_H
  12 +#define STIM_CELLSET_H
  13 +
  14 +#include <stim/math/vec3.h>
  15 +#include <vector>
  16 +#include <list>
  17 +#include <unordered_map>
  18 +#include <fstream>
  19 +
  20 +namespace stim{
  21 +
  22 +class cellset{
  23 +private:
  24 + static const char delim = ' ';
  25 +protected:
  26 + std::vector<double*> cells; //vector storing field data for each cell
  27 + std::unordered_map<std::string, size_t> fields; //unordered map storing field->index information for each field
  28 + size_t ip[3]; //hard code to position indices (for speed)
  29 +
  30 + void init(){
  31 +
  32 + }
  33 +
  34 + /// Initialize fields for a standard cell set containing three points and a radius
  35 + void init_p3(){
  36 + fields.insert(std::pair<std::string, size_t>("x", 0));
  37 + ip[0] = 0;
  38 + fields.insert(std::pair<std::string, size_t>("y", 1));
  39 + ip[1] = 1;
  40 + fields.insert(std::pair<std::string, size_t>("z", 2));
  41 + ip[2] = 2;
  42 + }
  43 +
  44 + void init_p3r(){
  45 + init_p3();
  46 + fields.insert(std::pair<std::string, size_t>("radius", 3));
  47 + ip[3] = 3;
  48 + }
  49 +public:
  50 + /// Constructor - create an empty cell set
  51 + cellset(){
  52 + init();
  53 + }
  54 +
  55 + /// Constructor - load a cellset from a file
  56 + cellset(std::string filename){
  57 + init(); //initialize an empty cellset
  58 + load(filename); //load the cellset from an existing file
  59 + }
  60 +
  61 + /// Loads a cellset from a file
  62 + void load(std::string filename){
  63 + std::ifstream infile(filename);
  64 + std::string header; //allocate space for the file header
  65 + std::getline(infile, header); //get the file header
  66 +
  67 + // break the header into fields
  68 + std::stringstream ss(header); //create a string stream
  69 + std::string field; //store a single field name
  70 + size_t i = 0; //current field index
  71 + while (std::getline(ss, field, delim)) { //split the header into individual fields
  72 + std::pair<std::string, size_t> p(field, i); //create a pair associating the header name with the index
  73 + fields.insert(p); //insert the pair into the fields map
  74 + i++; //increment the data index
  75 + }
  76 + size_t nfields = fields.size(); //store the number of fields for each cell
  77 +
  78 + //load each cell and all associated fields
  79 + std::string cell_line; //string holds all information for a cell
  80 + std::list<std::string> cell_list; //list will be temporary storage for the cell fields
  81 + while(std::getline(infile, cell_line)){ //for each cell entry
  82 + cell_list.push_back(cell_line); //push the cell entry into the list
  83 + }
  84 +
  85 + //convert the list into actual data
  86 + size_t ncells = cell_list.size(); //count the number of cells
  87 + cells.resize(ncells); //allocate enough space in the array to store all cells
  88 + for(size_t c = 0; c < ncells; c++){ //for each cell entry in the list
  89 + cells[c] = (double*) malloc(sizeof(double) * nfields); //allocate enough space for each field
  90 + std::stringstream fss(cell_list.front()); //turn the string representing the cell list into a stringstream
  91 + for(size_t f = 0; f < nfields; f++){ //for each field
  92 + fss>>cells[c][f]; //load the field
  93 + }
  94 + cell_list.pop_front(); //pop the read string off of the front of the list
  95 + }
  96 + infile.close(); //close the input file
  97 +
  98 + ip[0] = fields["x"]; //hard code the position indices for speed
  99 + ip[1] = fields["y"]; // this assumes all cells have positions
  100 + ip[2] = fields["z"];
  101 + }
  102 +
  103 + /// Return the value a specified field for a cell
  104 + /// @param c is the cell index
  105 + /// @param f is the field
  106 + double value(size_t c, std::string f){
  107 + size_t idx = fields[f];
  108 + return cells[c][idx];
  109 + }
  110 +
  111 + /// returns an ID used to look up a field
  112 + bool exists(std::string f){
  113 + std::unordered_map<std::string, size_t>::iterator iter = fields.find(f);
  114 + if(iter == fields.end()) return false;
  115 + else return true;
  116 + }
  117 +
  118 + /// Return the position of cell [i]
  119 + stim::vec3<double> p(size_t i){
  120 + stim::vec3<double> pos(cells[i][ip[0]], cells[i][ip[1]], cells[i][ip[2]]);
  121 + return pos;
  122 + }
  123 +
  124 + /// Return the number of cells in the set
  125 + size_t size(){
  126 + return cells.size();
  127 + }
  128 +
  129 + /// Return the maximum value of a field in this cell set
  130 + double maximum(std::string field){
  131 + size_t idx = fields[field]; //get the field index
  132 + size_t ncells = cells.size(); //get the total number of cells
  133 + double maxval, val; //stores the current and maximum values
  134 + for(size_t c = 0; c < ncells; c++){ //for each cell
  135 + val = cells[c][idx]; //get the field value for this cell
  136 + if(c == 0) maxval = val; //if this is the first cell, just assign the maximum
  137 + else if(val > maxval) maxval = val; // otherwise text for the size of val and assign it as appropriate
  138 + }
  139 + return maxval;
  140 + }
  141 +
  142 + /// Return the maximum value of a field in this cell set
  143 + double minimum(std::string field){
  144 + size_t idx = fields[field]; //get the field index
  145 + size_t ncells = cells.size(); //get the total number of cells
  146 + double minval, val; //stores the current and maximum values
  147 + for(size_t c = 0; c < ncells; c++){ //for each cell
  148 + val = cells[c][idx]; //get the field value for this cell
  149 + if(c == 0) minval = val; //if this is the first cell, just assign the maximum
  150 + else if(val < minval) minval = val; // otherwise text for the size of val and assign it as appropriate
  151 + }
  152 + return minval;
  153 + }
  154 +
  155 + /// adds a cell to the cell set
  156 + void add(double x, double y, double z, double r = 0){
  157 +
  158 + if(cells.size() == 0){ //if the cell set is empty
  159 + if(r == 0) //if the radius is zero
  160 + init_p3(); //initialize without a radius
  161 + else
  162 + init_p3r(); //otherwise initialize with a radius
  163 + }
  164 +
  165 + size_t nf = fields.size(); //get the number of fields
  166 + size_t bytes = sizeof(double) * nf; //get the size of a cell field
  167 + double* newcell = (double*) malloc(bytes); //allocate memory for the new cell
  168 + memset(newcell, 0, bytes); //initialize all fields to zero
  169 + newcell[ip[0]] = x;
  170 + newcell[ip[1]] = y;
  171 + newcell[ip[2]] = z;
  172 +
  173 + if(r != 0){ //if the user specifies a radius
  174 + size_t ir = fields["radius"]; //get the index for the radius field
  175 + newcell[ir] = r; //add the radius to the field
  176 + }
  177 + cells.push_back(newcell); //push the new memory entry into the cell array
  178 + }
  179 +
  180 +void save(std::string filename){
  181 +
  182 +
  183 + size_t ncells = cells.size(); // get the number of cells
  184 + std::ofstream file(filename); //open a file to store the cell's coordinates
  185 + if (file.is_open()) {
  186 +
  187 + file << "x y z radius\n"; //write the file header
  188 + for (size_t c=0; c < ncells; c++){ //for each cell
  189 + if (cells[c][ip[3]] != NULL) //check if for the current cell, radius has been assigned
  190 + file << cells[c][ip[0]] << delim << cells[c][ip[1]] << delim << cells[c][ip[2]] << delim << cells[c][ip[3]] << '\n' ;
  191 + else //check if for the current cell, radius has not been assigned, set radius to zero
  192 + file << cells[c][ip[0]] << delim << cells[c][ip[1]] << delim << cells[c][ip[2]] << delim << 0 << '\n' ;
  193 + }
  194 +
  195 + }
  196 + file.close();
  197 +
  198 + }
  199 +
  200 +
  201 +}; //end class cellset
  202 +}; //end namespace stim
  203 +
  204 +#endif
... ...
tira/biomodels/centerline.h 0 → 100644
  1 +#ifndef JACK_CENTERLINE_H
  2 +#define JACK_CENTERLINE_H
  3 +
  4 +#include <stdlib.h>
  5 +#include <vector>
  6 +#include <stim/math/vec3.h>
  7 +#include <stim/structures/kdtree.cuh>
  8 +
  9 +namespace stim {
  10 +
  11 + /// we always assume that one centerline has a flow direction even it actually does not have. Also, we allow loop centerline
  12 + /// NOTE: centerline is not derived from std::vector<stim::vec3<T>> class!!!
  13 + template<typename T>
  14 + class centerline {
  15 +
  16 + private:
  17 + size_t n; // number of points on the centerline, can be zero if NULL
  18 +
  19 + // update length information at each point (distance from starting point) starting from index "start"
  20 + void update_L(size_t start = 0) {
  21 + L.resize(n);
  22 +
  23 + if (start == 0) {
  24 + L[0] = 0.0f;
  25 + start++;
  26 + }
  27 +
  28 + stim::vec3<T> dir; // temp direction vector for calculating length between two points
  29 + for (size_t i = start; i < n; i++) {
  30 + dir = C[i] - C[i - 1]; // calculate the certerline extending direction
  31 + L[i] = L[i - 1] + dir.len(); // addition
  32 + }
  33 + }
  34 +
  35 + protected:
  36 + std::vector<stim::vec3<T> > C; // points on the centerline
  37 + std::vector<T> L; // stores the integrated length along the fiber
  38 +
  39 + public:
  40 + /// constructors
  41 + // empty constructor
  42 + centerline() {
  43 + n = 0;
  44 + }
  45 +
  46 + // constructor that allocate memory
  47 + centerline(size_t s) {
  48 + n = s;
  49 + C.resize(s); // allocate memory for points
  50 + L.resize(s); // allocate memory for lengths
  51 +
  52 + update_L();
  53 + }
  54 +
  55 + // constructor that constructs a centerline based on a list of points
  56 + centerline(std::vector<stim::vec3<T> > rhs) {
  57 + n = rhs.size(); // get the number of points
  58 + C.resize(n);
  59 + for (size_t i = 0; i < n; i++)
  60 + C[i] = rhs[i]; // copy data
  61 + update_L();
  62 + }
  63 +
  64 +
  65 + /// vector operations
  66 + // add a new point to current centerline
  67 + void push_back(stim::vec3<T> p) {
  68 + C.push_back(p);
  69 + n++; // increase the number of points
  70 + update_L(n - 1);
  71 + }
  72 +
  73 + // insert a new point at specific location to current centerline
  74 + void insert(typename std::vector<stim::vec3<T> >::iterator pos, stim::vec3<T> p) {
  75 + C.insert(pos, p); // insert a new point
  76 + n++;
  77 + size_t d = std::distance(C.begin(), pos); // get the index
  78 + update_L(d);
  79 + }
  80 + // insert a new point at C[idx]
  81 + void insert(size_t idx, stim::vec3<T> p) {
  82 + n++;
  83 + C.resize(n);
  84 + for (size_t i = n - 1; i > idx; i--) // move point location
  85 + C[i] = C[i - 1];
  86 + C[idx] = p;
  87 + update_L(idx);
  88 + }
  89 +
  90 + // assign a point at specific location to current centerline
  91 + void assign(size_t idx, stim::vec3<T> p) {
  92 + C[idx] = p;
  93 + update_L(idx);
  94 + }
  95 +
  96 + // erase a point at specific location on current centerline
  97 + void erase(typename std::vector<stim::vec3<T> >::iterator pos) {
  98 + C.erase(pos); // erase a point
  99 + n--;
  100 + size_t d = std::distance(C.begin(), pos); // get the index
  101 + update_L(d);
  102 + }
  103 + // erase a point at C[idx]
  104 + void erase(size_t idx) {
  105 + n--;
  106 + for (size_t i = idx; i < n; i++)
  107 + C[i] = C[i + 1];
  108 + C.resize(n);
  109 + update_L(idx);
  110 + }
  111 +
  112 + // clear up all the points
  113 + void clear() {
  114 + C.clear(); // clear list
  115 + L.clear(); // clear length information
  116 + n = 0; // set number to zero
  117 + }
  118 +
  119 + // reverse current centerline in terms of points order
  120 + centerline<T> reverse() {
  121 + centerline<T> result = *this;
  122 +
  123 + std::reverse(result.C.begin(), result.C.end());
  124 + result.update_L();
  125 +
  126 + return result;
  127 + }
  128 +
  129 + /// functions for reading centerline information
  130 + // return the number of points on current centerline
  131 + size_t size() {
  132 + return n;
  133 + }
  134 +
  135 + // return the length
  136 + T length() {
  137 + return L.back();
  138 + }
  139 +
  140 + // finds the index of the point closest to the length "l" on the lower bound
  141 + size_t findIdx(T l) {
  142 + for (size_t i = 1; i < L.size(); i++) {
  143 + if (L[i] > l)
  144 + return i - 1;
  145 + }
  146 +
  147 + return L.size() - 1;
  148 + }
  149 +
  150 + // get a position vector at the given length into the fiber (based on the pvalue), interpolate
  151 + stim::vec3<T> p(T l, size_t idx) {
  152 + T rate = (l - L[idx]) / (L[idx + 1] - L[idx]);
  153 + stim::vec3<T> v1 = C[idx];
  154 + stim::vec3<T> v2 = C[idx + 1];
  155 +
  156 + return (v1 + (v2 - v1) * rate);
  157 + }
  158 + // get a position vector at the given pvalue(pvalue[0.0f, 1.0f])
  159 + stim::vec3<T> p(T pvalue) {
  160 + // degenerated cases
  161 + if (pvalue <= 0.0f) return C[0];
  162 + if (pvalue >= 1.0f) return C.back();
  163 +
  164 + T l = pvalue * L.back(); // get the length based on the given pvalue
  165 + size_t idx = findIdx(l);
  166 +
  167 + return p(l, idx); // interpolation
  168 + }
  169 +
  170 + // get the normalized direction vector at point idx (average of the incoming and outgoing directions)
  171 + stim::vec3<T> d(size_t idx) {
  172 + if (n <= 1) return stim::vec3<T>(0.0f, 0.0f, 0.0f); // if there is insufficient information to calculate the direction, return null
  173 + if (n == 2) return (C[1] - C[0]).norm(); // if there are only two points, the direction vector at both is the direction of the line segment
  174 +
  175 + // degenerate cases at two ends
  176 + if (idx == 0) return (C[1] - C[0]).norm(); // the first direction vector is oriented towards the first line segment
  177 + if (idx == n - 1) return (C[n - 1] - C[n - 2]).norm(); // the last direction vector is oriented towards the last line segment
  178 +
  179 + // all other direction vectors are the average direction of the two joined line segments
  180 + stim::vec3<T> a = C[idx] - C[idx - 1];
  181 + stim::vec3<T> b = C[idx + 1] - C[idx];
  182 + stim::vec3<T> ab = a.norm() + b.norm();
  183 + return ab.norm();
  184 + }
  185 +
  186 +
  187 + /// arithmetic operations
  188 + // '=' operation
  189 + centerline<T> & operator=(centerline<T> rhs) {
  190 + n = rhs.n;
  191 + C = rhs.C;
  192 + L = rhs.L;
  193 +
  194 + return *this;
  195 + }
  196 +
  197 + // "[]" operation
  198 + stim::vec3<T> & operator[](size_t idx) {
  199 + return C[idx];
  200 + }
  201 +
  202 + // "==" operation
  203 + bool operator==(centerline<T> rhs) const {
  204 +
  205 + if (n != rhs.size())
  206 + return false;
  207 + else {
  208 + size_t num = rhs.size();
  209 + stim::vec3<T> tmp; // weird situation that I can only use tmp instead of C itself in comparison
  210 + for (size_t i = 0; i < num; i++) {
  211 + stim::vec3<T> tmp = C[i];
  212 + if (tmp[0] != rhs[i][0] || tmp[1] != rhs[i][1] || tmp[2] != rhs[i][2])
  213 + return false;
  214 + }
  215 + return true;
  216 + }
  217 + }
  218 +
  219 + // "+" operation
  220 + centerline<T> operator+(stim::vec3<T> p) const {
  221 + centerline<T> result(*this);
  222 + result.C.push_back(p);
  223 + result.n++;
  224 + result.update_L(n - 1);
  225 +
  226 + return result;
  227 + }
  228 + centerline<T> operator+(centerline<T> c) const {
  229 + centerline<T> result(*this);
  230 + size_t num1 = result.size();
  231 + size_t num2 = c.size();
  232 + for (size_t i = 0; i < num2; i++)
  233 + result.push_back(c[i]);
  234 + result.update_L(num1);
  235 +
  236 + return result;
  237 + }
  238 +
  239 +
  240 + /// advanced operation
  241 + // stitch two centerlines if possible (mutual-stitch and self-stitch)
  242 + static std::vector<centerline<T> > stitch(centerline<T> c1, centerline<T> c2, size_t end = 0) {
  243 + std::vector<centerline<T> > result;
  244 + centerline<T> new_centerline;
  245 + stim::vec3<T> new_vertex;
  246 + // ********** for Pavel **********
  247 + // ********** JACK thinks that ultimately we want it AUTOMATEDLY! **********
  248 +
  249 + // check stitch case
  250 + if (c1 == c2) { // self-stitch case
  251 + // ***** don't know how it works *****
  252 + }
  253 + else { // mutual-stitch case
  254 + size_t num1 = c1.size(); // get the numbers of two centerlines
  255 + size_t num2 = c2.size();
  256 +
  257 + T* reference = (T*)malloc(sizeof(T) * num1 * 3); // c1 as reference set
  258 + T* query = (T*)malloc(sizeof(T) * num2 * 3); // c2 as query set
  259 + for (size_t p = 0; p < num1; p++) // read points
  260 + for (size_t d = 0; d < 3; d++)
  261 + reference[p * 3 + d] = c1[p][d]; // KDTREE is stilla close code, it has its own structure
  262 + for (size_t p = 0; p < num2; p++) // read points
  263 + for (size_t d = 0; d < 3; d++)
  264 + query[p * 3 + d] = c2[p][d];
  265 +
  266 + stim::kdtree<T, 3> kdt;
  267 + kdt.create(reference, num1, 5); // create a tree based on reference set, max_tree_level is set to 5
  268 +
  269 + std::vector<size_t> index(num2); // stores index results
  270 + std::vector<float> dist(num2);
  271 +
  272 + kdt.search(query, num2, &index[0], &dist[0]); // find the nearest neighbor in c1 for c2
  273 +
  274 + // clear up
  275 + free(reference);
  276 + free(query);
  277 +
  278 + std::vector<float>::iterator sm = std::min_element(dist.begin(), dist.end()); // find smallest distance
  279 + size_t id = std::distance(dist.begin(), sm); // find the target index
  280 +
  281 + size_t id1 = index[id];
  282 + size_t id2 = id;
  283 +
  284 + // for centerline c1
  285 + bool flag = false; // flag indicates whether it has already added the bifurcation point to corresponding new centerline
  286 + if (id1 == 0 || id1 == num1 - 1) { // splitting bifurcation is on the end
  287 + new_centerline = c1;
  288 + new_vertex = c2[id2];
  289 + if (id1 == 0) {
  290 + new_centerline.insert(0, new_vertex);
  291 + flag = true;
  292 + }
  293 + if (id1 == num1 - 1) {
  294 + new_centerline.push_back(new_vertex);
  295 + flag = true;
  296 + }
  297 +
  298 + result.push_back(new_centerline);
  299 + }
  300 + else { // splitting bifurcation is on the centerline
  301 + std::vector<centerline<T> > tmp_centerline = c1.split(id1);
  302 + result = tmp_centerline;
  303 + }
  304 +
  305 + // for centerline c2
  306 + if (id2 == 0 || id2 == num2 - 1) { // splitting bidurcation is on the end
  307 + new_centerline = c2;
  308 + if (flag)
  309 + result.push_back(new_centerline);
  310 + else { // add the bifurcation point to this centerline
  311 + new_vertex = c1[id1];
  312 + if (id2 == 0) {
  313 + new_centerline.insert(0, new_vertex);
  314 + flag = true;
  315 + }
  316 + if (id2 == num2 - 1) {
  317 + new_centerline.push_back(new_vertex);
  318 + flag = true;
  319 + }
  320 +
  321 + result.push_back(new_centerline);
  322 + }
  323 + }
  324 + else { // splitting bifurcation is on the centerline
  325 + std::vector<centerline<T> > tmp_centerline = c2.split(id2);
  326 + result.push_back(tmp_centerline[0]);
  327 + result.push_back(tmp_centerline[1]);
  328 + }
  329 + }
  330 +
  331 + return result;
  332 + }
  333 +
  334 + // split current centerline at specific position
  335 + std::vector<centerline<T> > split(size_t idx) {
  336 + std::vector<centerline<T> > result;
  337 +
  338 + // won't split
  339 + if (idx <= 0 || idx >= n - 1) {
  340 + result.resize(1);
  341 + result[0] = *this; // return current centerline
  342 + }
  343 + // do split
  344 + else {
  345 + size_t n1 = idx + 1; // vertex idx would appear twice
  346 + size_t n2 = n - idx; // in total n + 1 points
  347 +
  348 + centerline<T> tmp; // temp centerline
  349 +
  350 + result.resize(2);
  351 +
  352 + for (size_t i = 0; i < n1; i++) // first half
  353 + tmp.push_back(C[i]);
  354 + tmp.update_L();
  355 + result[0] = tmp;
  356 + tmp.clear(); // clear up for next computation
  357 +
  358 + for (size_t i = 0; i < n2; i++) // second half
  359 + tmp.push_back(C[i + idx]);
  360 + tmp.update_L();
  361 + result[1] = tmp;
  362 + }
  363 +
  364 + return result;
  365 + }
  366 +
  367 + // resample current centerline
  368 + centerline<T> resample(T spacing) {
  369 +
  370 + stim::vec3<T> dir; // direction vector
  371 + stim::vec3<T> tmp; // intermiate point to be added
  372 + stim::vec3<T> p1; // starting point
  373 + stim::vec3<T> p2; // ending point
  374 +
  375 + centerline<T> result;
  376 +
  377 + for (size_t i = 0; i < n - 1; i++) {
  378 + p1 = C[i];
  379 + p2 = C[i + 1];
  380 +
  381 + dir = p2 - p1; // compute the direction of current segment
  382 + T seg_len = dir.len();
  383 +
  384 + if (seg_len > spacing) { // current segment can be sampled
  385 + for (T step = 0.0f; step < seg_len; step += spacing) {
  386 + tmp = p1 + dir * (step / seg_len); // add new point
  387 + result.push_back(tmp);
  388 + }
  389 + }
  390 + else
  391 + result.push_back(p1); // push back starting point
  392 + }
  393 + result.push_back(p2); // push back ending point
  394 +
  395 + return result;
  396 + }
  397 + };
  398 +}
  399 +
  400 +#endif
... ...
tira/biomodels/centerline_dep.h 0 → 100644
  1 +/*
  2 +Copyright <2017> <David Mayerich>
  3 +
  4 +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
  5 +
  6 +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
  7 +
  8 +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  9 +*/
  10 +
  11 +#ifndef STIM_CENTERLINE_H
  12 +#define STIM_CENTERLINE_H
  13 +
  14 +#include <vector>
  15 +#include <stim/math/vec3.h>
  16 +
  17 +namespace stim{
  18 +
  19 +/** This class stores information about a single fiber represented as a set of geometric points
  20 + * between two branch or end points. This class is used as a fundamental component of the stim::network
  21 + * class to describe an interconnected (often biological) network.
  22 + */
  23 +template<typename T>
  24 +class centerline{
  25 +
  26 +protected:
  27 + unsigned int N; //number of points in the fiber
  28 + double **c; //centerline (array of double pointers)
  29 +
  30 + /// Initialize an empty fiber
  31 + void init(){
  32 + N=0;
  33 + c=NULL;
  34 + }
  35 +
  36 + /// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
  37 + void init(unsigned int n){
  38 +
  39 + N = n; //set the number of points
  40 + c = (double**) malloc(sizeof(double*) * N); //allocate the array pointer
  41 +
  42 + for(unsigned int i = 0; i < N; i++) //allocate space for each point
  43 + c[i] = (double*) malloc(sizeof(double) * 3);
  44 + }
  45 +
  46 + /// Copies an existing fiber to the current fiber
  47 +
  48 + /// @param cpy stores the new copy of the fiber
  49 + void copy( const stim::centerline<T>& cpy, bool kd = 0){
  50 +
  51 + ///allocate space for the new fiber
  52 + init(cpy.N);
  53 +
  54 + ///copy the points
  55 + for(unsigned int i = 0; i < N; i++){
  56 + for(unsigned int d = 0; d < 3; d++) //for each dimension
  57 + c[i][d] = cpy.c[i][d]; //copy the coordinate
  58 + }
  59 + }
  60 +
  61 + /// find distance between two points
  62 + double dist(double* p0, double* p1){
  63 +
  64 + double sum = 0; // initialize variables
  65 + float v;
  66 + for(unsigned int d = 0; d < 3; d++)
  67 + {
  68 + v = p1[d] - p0[d];
  69 + sum +=v * v;
  70 + }
  71 + return sqrt(sum);
  72 + }
  73 +
  74 + /// Returns a stim::vec representing the point at index i
  75 +
  76 + /// @param i is an index of the desired centerline point
  77 + stim::vec<T> get_vec(unsigned i){
  78 + stim::vec3<T> r;
  79 + r.resize(3);
  80 + r[0] = c[i][0];
  81 + r[1] = c[i][1];
  82 + r[2] = c[i][2];
  83 +
  84 + return r;
  85 + }
  86 +
  87 +
  88 +public:
  89 +
  90 + centerline(){
  91 + init();
  92 + }
  93 +
  94 + /// Copy constructor
  95 + centerline(const stim::centerline<T> &obj){
  96 + copy(obj);
  97 + }
  98 +
  99 + //initialize a centerline with n points
  100 + centerline(int n){
  101 + init(n);
  102 + }
  103 +
  104 + /// Constructor takes a list of stim::vec points, the radius at each point is set to zero
  105 + centerline(std::vector< stim::vec<T> > p, bool kd = 0){
  106 + init(p.size()); //initialize the fiber
  107 +
  108 + //for each point, set the centerline position and radius
  109 + for(unsigned int i = 0; i < N; i++){
  110 +
  111 + //set the centerline position
  112 + for(unsigned int d = 0; d < 3; d++)
  113 + c[i][d] = (double) p[i][d];
  114 + }
  115 + }
  116 +
  117 + /// constructor takes a list of points
  118 + centerline(std::vector< stim::vec3< T > > pos, bool kd = 0){
  119 + init(pos.size()); //initialize the fiber
  120 +
  121 + //for each point, set the centerline position and radius
  122 + for(unsigned int i = 0; i < N; i++){
  123 + //set the centerline position
  124 + for(unsigned int d = 0; d < 3; d++)
  125 + c[i][d] = (double) pos[i][d];
  126 + }
  127 + }
  128 +
  129 + /// Assignment operation
  130 + centerline& operator=(const centerline &rhs){
  131 + if(this == &rhs) return *this; //test for and handle self-assignment
  132 + copy(rhs);
  133 + return *this;
  134 + }
  135 +
  136 +
  137 + /// Returns the fiber centerline as an array of stim::vec points
  138 + std::vector< stim::vec<T> > get_centerline(){
  139 +
  140 + //create an array of stim vectors
  141 + std::vector< stim::vec3<T> > pts(N);
  142 +
  143 + //cast each point to a stim::vec, keeping only the position information
  144 + for(unsigned int i = 0; i < N; i++)
  145 + pts[i] = stim::vec3<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]);
  146 +
  147 + //return the centerline array
  148 + return pts;
  149 + }
  150 +
  151 + /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
  152 + std::vector< stim::centerline<T> > split(unsigned int idx){
  153 +
  154 + std::vector< stim::centerline<T> > fl; //create an array to store up to two fibers
  155 +
  156 + //if the index is an end point, only the existing fiber is returned
  157 + if(idx == 0 || idx == N-1){
  158 + fl.resize(1); //set the size of the fiber to 1
  159 + fl[0] = *this; //copy the current fiber
  160 + }
  161 +
  162 + //if the index is not an end point
  163 + else{
  164 +
  165 + unsigned int N1 = idx + 1; //calculate the size of both fibers
  166 + unsigned int N2 = N - idx;
  167 +
  168 + fl.resize(2); //set the array size to 2
  169 +
  170 + fl[0].init(N1); //set the size of each fiber
  171 + fl[1].init(N2);
  172 +
  173 + //copy both halves of the fiber
  174 + unsigned int i, d;
  175 +
  176 + //first half
  177 + for(i = 0; i < N1; i++){ //for each centerline point
  178 + for(d = 0; d < 3; d++)
  179 + fl[0].c[i][d] = c[i][d]; //copy each coordinate
  180 + }
  181 +
  182 + //second half
  183 + for(i = 0; i < N2; i++){
  184 + for(d = 0; d < 3; d++)
  185 + fl[1].c[i][d] = c[idx + i][d];
  186 +
  187 + }
  188 + }
  189 +
  190 + return fl; //return the array
  191 +
  192 + }
  193 +
  194 + /// Outputs the fiber as a string
  195 + std::string str(){
  196 + std::stringstream ss;
  197 +
  198 + //create an iterator for the point list
  199 + //typename std::list< point<T> >::iterator i;
  200 + for(unsigned int i = 0; i < N; i++){
  201 + ss<<" [ ";
  202 + for(unsigned int d = 0; d < 3; d++){
  203 + ss<<c[i][d]<<" ";
  204 + }
  205 + }
  206 +
  207 + return ss.str();
  208 + }
  209 + /// Returns the number of centerline points in the fiber
  210 + unsigned int size(){
  211 + return N;
  212 + }
  213 +
  214 +
  215 + /// Bracket operator returns the element at index i
  216 +
  217 + /// @param i is the index of the element to be returned as a stim::vec
  218 + stim::vec<T> operator[](unsigned i){
  219 + return get_vec(i);
  220 + }
  221 +
  222 + /// Back method returns the last point in the fiber
  223 + stim::vec<T> back(){
  224 + return get_vec(N-1);
  225 + }
  226 + ////resample a fiber in the network
  227 + stim::centerline<T> resample(T spacing)
  228 + {
  229 + std::cout<<"fiber::resample()"<<std::endl;
  230 +
  231 + std::vector<T> v(3); //v-direction vector of the segment
  232 + stim::vec<T> p(3); //- intermediate point to be added
  233 + stim::vec<T> p1(3); // p1 - starting point of an segment on the fiber,
  234 + stim::vec<T> p2(3); // p2 - ending point,
  235 + double sum=0; //distance summation
  236 + std::vector<stim::vec<T> > fiberPositions = centerline();
  237 + std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
  238 + // for each point on the centerline (skip if it is the last point on centerline)
  239 + for(unsigned int f=0; f< N-1; f++)
  240 + {
  241 +
  242 + p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1;
  243 + for(unsigned int d = 0; d < 3; d++){
  244 + sum +=v[d] * v[d];} //length of segment-distance between starting and ending point
  245 +
  246 + T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment
  247 +
  248 + if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample
  249 + {
  250 + // repeat resampling until accumulated stepsize is equsl to length of the segment
  251 + for(T step=0.0; step<lengthSegment; step+=spacing)
  252 + {
  253 + // calculate the resampled point by travelling step size in the direction of normalized gradient vector
  254 + for(unsigned int i=0; i<3;i++)
  255 + {
  256 + p[i] = p1[i] + v[i]*(step/lengthSegment);
  257 + } //for each dimension
  258 + // add this resampled points to the new fiber list
  259 + newPointList.push_back(p);
  260 + }
  261 + }
  262 + else // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
  263 + newPointList.push_back(fiberPositions[f+1]);
  264 + }
  265 + newPointList.push_back(fiberPositions[N-1]); //add the last point on the fiber to the new fiber list
  266 + centerline newFiber(newPointList);
  267 + return newFiber;
  268 + }
  269 +
  270 +};
  271 +
  272 +
  273 +
  274 +} //end namespace stim
  275 +
  276 +
  277 +
  278 +#endif
... ...
tira/biomodels/flow.h 0 → 100644
  1 +/*
  2 +Copyright <2017> <David Mayerich>
  3 +
  4 +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
  5 +
  6 +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
  7 +
  8 +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  9 +*/
  10 +#ifndef STIM_FLOW_H
  11 +#define STIM_FLOW_H
  12 +
  13 +#include <vector>
  14 +#include <algorithm>
  15 +
  16 +//STIM include
  17 +#include <stim/math/vec3.h>
  18 +#include <stim/parser/arguments.h>
  19 +#include <stim/biomodels/network.h>
  20 +
  21 +#ifdef __CUDACC__
  22 +#include <cublas_v2.h>
  23 +#include <stim/cuda/cudatools/error.h>
  24 +#endif
  25 +
  26 +namespace stim {
  27 + template <typename A, typename B, typename C>
  28 + struct triple {
  29 + A first;
  30 + B second;
  31 + C third;
  32 + };
  33 +
  34 + template <typename T>
  35 + struct bridge {
  36 + std::vector<unsigned> v; // vertices' indices
  37 + std::vector<typename stim::vec3<T> > V; // vertices' coordinates
  38 + T l; // length
  39 + T r; // radii
  40 + T deltaP; // pressure drop
  41 + T Q; // volume flow rate
  42 + };
  43 +
  44 + template <typename T>
  45 + class flow {
  46 +
  47 + private:
  48 +
  49 + // calculate the cofactor of elemen[row][col]
  50 + void get_minor(T** src, T** dest, int row, int col, int order) {
  51 +
  52 + // index of element to be copied
  53 + int rowCount = 0;
  54 + int colCount = 0;
  55 +
  56 + for (int i = 0; i < order; i++) {
  57 + if (i != row) {
  58 + colCount = 0;
  59 + for (int j = 0; j < order; j++) {
  60 + // when j is not the element
  61 + if (j != col) {
  62 + dest[rowCount][colCount] = src[i][j];
  63 + colCount++;
  64 + }
  65 + }
  66 + rowCount++;
  67 + }
  68 + }
  69 + }
  70 +
  71 + // calculate the det()
  72 + T determinant(T** mat, int order) {
  73 +
  74 + // degenate case when n = 1
  75 + if (order == 1)
  76 + return mat[0][0];
  77 +
  78 + T det = 0.0; // determinant value
  79 +
  80 + // allocate the cofactor matrix
  81 + T** minor = (T**)malloc((order - 1) * sizeof(T*));
  82 + for (int i = 0; i < order - 1; i++)
  83 + minor[i] = (T*)malloc((order - 1) * sizeof(T));
  84 +
  85 +
  86 + for (int i = 0; i < order; i++) {
  87 +
  88 + // get minor of element(0, i)
  89 + get_minor(mat, minor, 0, i, order);
  90 +
  91 + // recursion
  92 + det += (i % 2 == 1 ? -1.0 : 1.0) * mat[0][i] * determinant(minor, order - 1);
  93 + }
  94 +
  95 + // release memory
  96 + for (int i = 0; i < order - 1; i++)
  97 + free(minor[i]);
  98 + free(minor);
  99 +
  100 + return det;
  101 + }
  102 +
  103 + public:
  104 + T** C; // Conductance
  105 + std::vector<typename stim::triple<unsigned, unsigned, float> > Q; // volume flow rate
  106 + std::vector<T> QQ; // Q' vector
  107 + std::vector<T> P; // initial pressure
  108 + std::vector<T> pressure; // final pressure
  109 +
  110 + //std::vector<typename stim::triple<unsigned, unsigned, T> > V; // velocity
  111 + //std::vector<typename stim::triple<unsigned, unsigned, T> > Q; // volume flow rate
  112 + //std::vector<typename stim::triple<unsigned, unsigned, T> > deltaP; // pressure drop
  113 +
  114 + flow() {} // default constructor
  115 +
  116 + void init(unsigned n_e, unsigned n_v) {
  117 +
  118 + C = new T*[n_v]();
  119 + for (unsigned i = 0; i < n_v; i++) {
  120 + C[i] = new T[n_v]();
  121 + }
  122 +
  123 + QQ.resize(n_v);
  124 + P.resize(n_v);
  125 + pressure.resize(n_v);
  126 +
  127 + Q.resize(n_e);
  128 + }
  129 +
  130 + void reset(unsigned n_v) {
  131 +
  132 + for (unsigned i = 0; i < n_v; i++) {
  133 + for (unsigned j = 0; j < n_v; j++) {
  134 + C[i][j] = 0;
  135 + }
  136 + }
  137 + }
  138 +
  139 + void clear(unsigned n_v) {
  140 +
  141 + for (unsigned i = 0; i < n_v; i++)
  142 + delete[] C[i];
  143 + delete[] C;
  144 + }
  145 +
  146 + /// Calculate the inverse of A and store the result in C
  147 + void inversion(T** A, int order, T* C) {
  148 +
  149 +#ifdef __CUDACC__
  150 +
  151 + // convert from double pointer to single pointer, make it flat
  152 + T* Aflat = (T*)malloc(order * order * sizeof(T));
  153 + for (unsigned i = 0; i < order; i++)
  154 + for (unsigned j = 0; j < order; j++)
  155 + Aflat[i * order + j] = A[i][j];
  156 +
  157 + // create device pointer
  158 + T* d_Aflat; // flat original matrix
  159 + T* d_Cflat; // flat inverse matrix
  160 + T** d_A; // put the flat original matrix into another array of pointer
  161 + T** d_C;
  162 + int *d_P;
  163 + int *d_INFO;
  164 +
  165 + // allocate memory on device
  166 + HANDLE_ERROR(cudaMalloc((void**)&d_Aflat, order * order * sizeof(T)));
  167 + HANDLE_ERROR(cudaMalloc((void**)&d_Cflat, order * order * sizeof(T)));
  168 + HANDLE_ERROR(cudaMalloc((void**)&d_A, sizeof(T*)));
  169 + HANDLE_ERROR(cudaMalloc((void**)&d_C, sizeof(T*)));
  170 + HANDLE_ERROR(cudaMalloc((void**)&d_P, order * 1 * sizeof(int)));
  171 + HANDLE_ERROR(cudaMalloc((void**)&d_INFO, 1 * sizeof(int)));
  172 +
  173 + // copy matrix from host to device
  174 + HANDLE_ERROR(cudaMemcpy(d_Aflat, Aflat, order * order * sizeof(T), cudaMemcpyHostToDevice));
  175 +
  176 + // copy matrix from device to device
  177 + HANDLE_ERROR(cudaMemcpy(d_A, &d_Aflat, sizeof(T*), cudaMemcpyHostToDevice));
  178 + HANDLE_ERROR(cudaMemcpy(d_C, &d_Cflat, sizeof(T*), cudaMemcpyHostToDevice));
  179 +
  180 + // calculate the inverse of matrix based on cuBLAS
  181 + cublasHandle_t handle;
  182 + CUBLAS_HANDLE_ERROR(cublasCreate_v2(&handle)); // create cuBLAS handle object
  183 +
  184 + CUBLAS_HANDLE_ERROR(cublasSgetrfBatched(handle, order, d_A, order, d_P, d_INFO, 1));
  185 +
  186 + int INFO = 0;
  187 + HANDLE_ERROR(cudaMemcpy(&INFO, d_INFO, sizeof(int), cudaMemcpyDeviceToHost));
  188 + if (INFO == order)
  189 + {
  190 + std::cout << "Factorization Failed : Matrix is singular." << std::endl;
  191 + cudaDeviceReset();
  192 + exit(1);
  193 + }
  194 +
  195 + CUBLAS_HANDLE_ERROR(cublasSgetriBatched(handle, order, (const T **)d_A, order, d_P, d_C, order, d_INFO, 1));
  196 +
  197 + CUBLAS_HANDLE_ERROR(cublasDestroy_v2(handle));
  198 +
  199 + // copy inverse matrix from device to device
  200 + HANDLE_ERROR(cudaMemcpy(&d_Cflat, d_C, sizeof(T*), cudaMemcpyDeviceToHost));
  201 +
  202 + // copy inverse matrix from device to host
  203 + HANDLE_ERROR(cudaMemcpy(C, d_Cflat, order * order * sizeof(T), cudaMemcpyDeviceToHost));
  204 +
  205 + // clear up
  206 + free(Aflat);
  207 + HANDLE_ERROR(cudaFree(d_Aflat));
  208 + HANDLE_ERROR(cudaFree(d_Cflat));
  209 + HANDLE_ERROR(cudaFree(d_A));
  210 + HANDLE_ERROR(cudaFree(d_C));
  211 + HANDLE_ERROR(cudaFree(d_P));
  212 + HANDLE_ERROR(cudaFree(d_INFO));
  213 +
  214 +#else
  215 + // get the determinant of a
  216 + double det = 1.0 / determinant(A, order);
  217 +
  218 + // memory allocation
  219 + T* tmp = (T*)malloc((order - 1)*(order - 1) * sizeof(T));
  220 + T** minor = (T**)malloc((order - 1) * sizeof(T*));
  221 + for (int i = 0; i < order - 1; i++)
  222 + minor[i] = tmp + (i * (order - 1));
  223 +
  224 + for (int j = 0; j < order; j++) {
  225 + for (int i = 0; i < order; i++) {
  226 + // get the co-factor (matrix) of A(j,i)
  227 + get_minor(A, minor, j, i, order);
  228 + C[i][j] = det * determinant(minor, order - 1);
  229 + if ((i + j) % 2 == 1)
  230 + C[i][j] = -C[i][j];
  231 + }
  232 + }
  233 +
  234 + // release memory
  235 + free(tmp);
  236 + free(minor);
  237 +#endif
  238 + }
  239 + };
  240 +}
  241 +
  242 +#endif
  243 +
  244 +
  245 +
  246 +//// calculate the flow rate of 3D model(circle cross section)
  247 +//void calculate_flow_rate(unsigned e, T r) {
  248 +// stim::triple<unsigned, unsigned, T> tmp_Q;
  249 +// tmp_Q.first = V[e].first; // copy the vertices information
  250 +// tmp_Q.second = V[e].second;
  251 +// tmp_Q.third = V[e].third * stim::PI * pow(r, 2); // UNITS: uL/s
  252 +// Q.push_back(tmp_Q); // push back the volume flow rate information for every edge
  253 +//}
  254 +
  255 +//// calculate the flow rate of 2D model(rectangular cross section)
  256 +//void calculate_flow_rate(unsigned e, T r, T h) {
  257 +// stim::triple<unsigned, unsigned, T> tmp_Q;
  258 +// tmp_Q.first = V[e].first; // copy the vertices information
  259 +// tmp_Q.second = V[e].second;
  260 +// tmp_Q.third = V[e].third * h * r; // UNITS: uL/s = mm^3/s
  261 +// Q.push_back(tmp_Q); // push back the volume flow rate information for every edge
  262 +//}
  263 +
  264 +//// calculate the pressure drop of 3D model(circle cross section)
  265 +//void calculate_deltaP(unsigned e, T u, T l, T r) {
  266 +// stim::triple<unsigned, unsigned, T> tmp_deltaP;
  267 +// tmp_deltaP.first = V[e].first; // copy the vertices information
  268 +// tmp_deltaP.second = V[e].second;
  269 +// tmp_deltaP.third = (8 * u * l * Q[e].third) / (stim::PI * pow(r, 4)); // UNITS: g/mm/s^2 = Pa
  270 +// deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge
  271 +//}
  272 +
  273 +//// calculate the pressure drop of 2D model(rectangular cross section)
  274 +//void calculate_deltaP(unsigned e, T u, T l, T r, T h) {
  275 +// stim::triple<unsigned, unsigned, T> tmp_deltaP;
  276 +// tmp_deltaP.first = V[e].first; // copy the vertices information
  277 +// tmp_deltaP.second = V[e].second;
  278 +// tmp_deltaP.third = (12 * u * l * Q[e].third) / (h * pow(r, 3)); // UNITS: g/mm/s^2 = Pa
  279 +// deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge
  280 +//}
  281 +
  282 +//// better way to do this???
  283 +//// find the maximum and minimum pressure positions
  284 +//void find_max_min_pressure(size_t n_e, size_t n_v, unsigned& max, unsigned& min) {
  285 +// std::vector<T> P(n_v, FLT_MAX);
  286 +// // set one to reference
  287 +// P[Q[0].first] = 0.0;
  288 +// unsigned first = 0;
  289 +// unsigned second = 0;
  290 +// // calculate all the relative pressure in brute force manner
  291 +// for (unsigned e = 0; e < n_e; e++) {
  292 +// // assuming the obj file stores in a straight order, in other words, like swc file
  293 +// first = Q[e].first;
  294 +// second = Q[e].second;
  295 +// if (P[first] != FLT_MAX) // if pressure at start vertex is known
  296 +// P[second] = P[first] - deltaP[e].third;
  297 +// else if (P[second] != FLT_MAX) // if pressure at end vertex is known
  298 +// P[first] = P[second] + deltaP[e].third;
  299 +// }
  300 +
  301 +// // find the maximum and minimum pressure position
  302 +// auto m1 = std::max_element(P.begin(), P.end()); // temporarily max number
  303 +// auto m2 = std::min_element(P.begin(), P.end()); // temporarily min number
  304 +
  305 +// max = std::distance(P.begin(), m1);
  306 +// min = std::distance(P.begin(), m2);
  307 +
  308 +// T tmp_m = *m2;
  309 +// // Now set the lowest pressure port to reference pressure(0.0 Pa)
  310 +// for (unsigned i = 0; i < n_v; i++)
  311 +// P[i] -= tmp_m;
  312 +
  313 +// for (unsigned i = 0; i < n_v; i++)
  314 +// pressure.push_back(P[i]);
  315 +//}
0 316 \ No newline at end of file
... ...
tira/biomodels/flow_dep.h 0 → 100644
  1 +/*
  2 +Copyright <2017> <David Mayerich>
  3 +
  4 +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
  5 +
  6 +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
  7 +
  8 +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  9 +*/
  10 +
  11 +#pragma once
  12 +#include <fstream> // Required for ofstream, etc.
  13 +#include <iomanip> // Required for setw
  14 +#include <iostream> // Required for cout, cin, etc.
  15 +#include <tuple> // Required for returning multiple values from a function
  16 +
  17 +using namespace std;
  18 +
  19 +
  20 +class flow
  21 +{
  22 +public:
  23 + void backupToTxt(unsigned int nL, double **D, char filename[]);
  24 + tuple<int, int> copySrcDesRadLen(char filename[]);
  25 + void copyToArray(int *src, int *dest, double *radii, double *len);
  26 + int getDangleNodes(int datarow, int numNodes, int *row, int *column, int *dangleNodes);
  27 + void inversion(double **a, int n, double **b);
  28 +
  29 +protected:
  30 + float determinant(double **a, int n);
  31 + int minor(double **src, double **dest, int row, int col, int order);
  32 +};
  33 +
  34 +/* Function to find the dangle nodes in a network */
  35 +// Created by Cherub P. Harder (8/10/2015), U of Houston
  36 +// Modified by Cherub P. Harder on 8/12/2015
  37 +int flow::getDangleNodes(int datarow, int numNodes, int *column1, int *column2, int *dangleNodes)
  38 +{
  39 + int count = datarow, diff1 = 0, diff2 = 0, numPress = 0, st = 0;
  40 +
  41 + // Find matching nodes within column2
  42 + for( int i = 0; i < count; i++ )
  43 + {
  44 + for( int y = i+1; y < datarow; y++ )
  45 + {
  46 + if( column2[i] == column2[y] ) // Is there a match?
  47 + {
  48 + st = column2[i]; // Save the matching node
  49 +// cout << endl << column2[i] << " = " << column2[y] << endl; // Display the matching nodes
  50 + memmove(column2+i, column2+i+1, (datarow-(i+1)) * sizeof(column2[0])); // Move up the rows
  51 + // taking the places of the rows before them starting
  52 + // with where the matching node is located
  53 + column2[datarow-1] = st; // Place the matching node at the very end of the array--
  54 + // this is for comparison purpose so that the other match-
  55 + // ing node will be moved as well and then omitted later.
  56 + diff1++; // Count the matching node
  57 +
  58 + // Display the updated array (with the matching node moved to the bottommost row)
  59 +/* cout << "Updated array:" << endl;
  60 + for( int k = 0; k < datarow; k++ )
  61 + cout << column2[k] << endl;
  62 +*/
  63 + // Decrement the counters
  64 + // NOTE: The counters need to be decremented because the rows had been moved up, so the same
  65 + // locations need to be read again because they contain different values now after the move.
  66 + i--; // Decrement i to read the node that took over the place
  67 + // of the matching node. Otherwise, it will be skipped.
  68 + y--; // Decrement y to read the following node for comparison
  69 + count--; // The maximum count need to be decremented so that the
  70 + // matching nodes that had been moved will not be read again.
  71 + // However, the maximum count (datarow) for finding a match
  72 + // will not be decremented because the remaining matching
  73 + // node that has not been moved yet needs to be moved and
  74 + // the only way to do that is to match it with its duplicate.
  75 + }
  76 + }
  77 + }
  78 +
  79 + // Store the nodes that have no duplicates
  80 + // NOTE: This will only save the nodes that have not been moved to the bottom.
  81 +// cout << "\ndangleNodes array:" << endl;
  82 + for( int j = 0; j < datarow-diff1; j++ )
  83 + {
  84 + dangleNodes[numPress] = column2[j];
  85 +// cout << dangleNodes[j] << endl; // DELETE!!!
  86 + numPress++; // Count the non-duplicated node
  87 + }
  88 +
  89 + // Find if the non-duplicated nodes have a match from column1
  90 + count = datarow-diff1; // Reinitialize the counter
  91 +
  92 + for( int i = 0; i < count; i++ )
  93 + {
  94 + for( int j = 0; j < datarow; j++ )
  95 + {
  96 + if( dangleNodes[i] == column1[j] ) // Is there a match?
  97 + {
  98 + st = column1[j]; // Save the matching node
  99 +// cout << endl << dangleNodes[i] << " = " << column1[j] << endl; // Display the matching nodes
  100 + memmove(dangleNodes+i, dangleNodes+i+1, (datarow-diff1-(i+1)) * sizeof(dangleNodes[0]));
  101 + dangleNodes[count-1] = st; // Move the matching node to the bottom of the array
  102 + diff2++; // Count the matching node
  103 +
  104 + // Display the updated array
  105 +/* cout << "Updated dangleNodes array:" << endl;
  106 + for( int k = 0; k < count-1; k++ )
  107 + {
  108 + cout << dangleNodes[k] << endl;
  109 + }
  110 +*/
  111 + // Decrement the counters
  112 + i--;
  113 + j--;
  114 + count--;
  115 + numPress--; // Decrement to the exact number of dangle nodes
  116 + }
  117 + }
  118 + }
  119 +
  120 + return numPress; // Number of dangle nodes
  121 +}
  122 +
  123 +
  124 +// Function to make a backup copy of the contents of a matrix to a .txt file
  125 +// Created by Cherub P. Harder (8/10/2015), U of Houston
  126 +void flow::backupToTxt(unsigned int nL, double **D, char filename[])
  127 +{
  128 + ofstream output_file(filename);
  129 +
  130 + for( unsigned int i = 0; i < nL; i++ )
  131 + {
  132 + for( int j = 0; j < 4; j++ )
  133 + {
  134 + if( j < 3 )
  135 + output_file << D[i][j] << "\t";
  136 +
  137 + else
  138 + output_file << D[i][j];
  139 + }
  140 +
  141 + output_file << "\n";
  142 + }
  143 +
  144 + output_file.close( );
  145 +}
  146 +
  147 +
  148 +// Function to make separate copies of the source nodes, destination nodes, radii, and lengths
  149 +// Created by Cherub P. Harder (8/10/2015), U of Houston
  150 +tuple<int, int> flow::copySrcDesRadLen(char filename[])
  151 +{
  152 + int cnt = 0, numElements = 0, numNodes = 0;
  153 + float number = 0.0;
  154 + ofstream srcData("srcCol.txt"); // A .txt file to store the source nodes
  155 + ofstream desData("destCol.txt"); // A .txt file to store the destination nodes
  156 + ofstream radiiData("radii.txt"); // A .txt file to store the radii
  157 + ofstream lenData("lengths.txt"); // A .txt file to store the lengths
  158 + FILE *fp = fopen(filename, "r"); // Create a variable of type FILE* and open the file using
  159 + // the fopen function and assign the file to the variable
  160 + // Check if the file exists
  161 + if(fp == NULL) // Alternative: if(!fp)
  162 + {
  163 + printf("Error! File does not exist.\n");
  164 + getchar( ); // Pause
  165 + exit(-1); // NOTE: Must include stdlib.h.
  166 + }
  167 +
  168 + // Store data to their respective .txt files
  169 + while(fscanf(fp, "%f", &number) == 1)
  170 + {
  171 + cnt++; // Increment counter
  172 +
  173 + // Store to srcCol.txt
  174 + if(cnt == 1)
  175 + srcData << number << endl;
  176 +
  177 + // Store to destCol.txt
  178 + if(cnt == 2)
  179 + desData << number << endl;
  180 +
  181 + // Save the current number of nodes
  182 + if(cnt < 3)
  183 + {
  184 + if(number > numNodes)
  185 + numNodes = (int)number;
  186 + }
  187 +
  188 + // Store to radii.txt
  189 + if(cnt == 3)
  190 + radiiData << number << endl;
  191 +
  192 + // Store to lengths.txt
  193 + if(cnt == 4)
  194 + {
  195 + lenData << number << endl;
  196 +
  197 + numElements++; // Count the elements
  198 + cnt = 0; // Reset counter
  199 + }
  200 + }
  201 +
  202 + srcData.close( );
  203 + desData.close( );
  204 + radiiData.close( );
  205 + lenData.close( );
  206 +
  207 + return make_tuple(numNodes, numElements); // Return two values
  208 +}
  209 +
  210 +
  211 +// Function to copy data for .txt files to their respective arrays
  212 +// Created by Cherub P. Harder (8/11/2015), U of Houston
  213 +void flow::copyToArray(int *src, int *dest, double *radii, double *len)
  214 +{
  215 + int v = 0;
  216 + double tmp = 0, R = 0, L = 0;
  217 +
  218 + // Store source node values to the array src
  219 + ifstream readSrc("srcCol.txt");
  220 +
  221 + while( readSrc >> tmp )
  222 + {
  223 + src[v] = (int)tmp;
  224 + v++;
  225 + }
  226 +
  227 + readSrc.close( );
  228 +
  229 + // Store destination node values to the array dest
  230 + v = 0; // Reset counter
  231 + ifstream readDest("destCol.txt");
  232 +
  233 + while( readDest >> tmp )
  234 + {
  235 + dest[v] = (int)tmp;
  236 + v++;
  237 + }
  238 +
  239 + readDest.close( );
  240 +
  241 + // Store radius values to the array radii
  242 + v = 0; // Reset counter
  243 + ifstream readRad("radii.txt");
  244 +
  245 + while( readRad >> tmp )
  246 + {
  247 + radii[v] = tmp;
  248 + v++;
  249 + }
  250 +
  251 + readRad.close( );
  252 +
  253 + // Store length values to the array len
  254 + v = 0; // Reset counter
  255 + ifstream readLen("lengths.txt");
  256 +
  257 + while( readLen >> tmp )
  258 + {
  259 + len[v] = tmp;
  260 + v++;
  261 + }
  262 +
  263 + readLen.close( );
  264 +}
  265 +
  266 +
  267 +// Function to find the inverse of a square matrix
  268 +void flow::inversion(double **a, int n, double **b)
  269 +{
  270 + // Get 1 over the determinant of A
  271 + double det = (double)(1.0/determinant(a, n));
  272 +// cerr << "\n1/det(C) = " << det << endl; // DELETE!!!
  273 +
  274 + // Memory allocation
  275 + double *tmp = new double[(n-1) * (n-1)];
  276 + double **m = new double * [n-1];
  277 + for( int i = 0; i < n-1; i++ )
  278 + m[i] = tmp + ( i * (n-1) );
  279 +
  280 + for( int j = 0; j < n; j++)
  281 + {
  282 + for( int i = 0; i < n; i++ )
  283 + {
  284 + // Get the cofactor (matrix) of a(j,i)
  285 + minor(a, m, j, i, n);
  286 + b[i][j] = det * determinant( m, n-1 );
  287 + if( (i+j)%2 == 1 )
  288 + b[i][j] = -b[i][j];
  289 + }
  290 + }
  291 +
  292 + // Release memory
  293 + // Delete [] minor[0];
  294 + delete [] tmp;
  295 + delete [] m;
  296 +}
  297 +
  298 +
  299 +// Function to find the determinant of a matrix using recursion
  300 +// Contribution by Edward Popko
  301 +// Modified by Cherub P. Harder (7/15/2015), U of Houston
  302 +// Arguments: a(double **) - pointer to a pointer of an arbitrary square matrix
  303 +// n(int) - dimension of the square matrix
  304 +float flow::determinant(double **a, int n)
  305 +{
  306 + int i, j, j1, j2; // General loop and matrix subscripts
  307 + double det = 0; // Initialize determinant
  308 + double **m = NULL; // Pointer to pointer to implement 2D square array
  309 +
  310 + // Display contents of matrix C (DELETE!!!)
  311 +/* std::cout << "\nThe updated matrix C:\n";
  312 + for( int j = 0; j < n; ++j )
  313 + {
  314 + std::cerr << "\t";
  315 +
  316 + for( int k = 0; k < n; ++k )
  317 + std::cerr << left << setw(15) << a[j][k];
  318 +
  319 + std::cerr << endl;
  320 + }
  321 +
  322 + getchar(); // DELETE!!!*/
  323 +
  324 + if(n < 1) { } // Error condition - should never get here
  325 +
  326 + else if (n == 1) // Should never get here
  327 + {
  328 + det = a[0][0];
  329 + }
  330 +
  331 + else if(n == 2) // Basic 2x2 sub-matrix determinate definition
  332 + { // When n == 2, this ends the recursion series
  333 + det = a[0][0] * a[1][1] - a[1][0] * a[0][1];
  334 + }
  335 + // Recursion continues, solve next sub-matrix
  336 + else // Solve the next minor by building a sub-matrix
  337 + {
  338 + det = 0; // Initialize determinant of sub-matrix
  339 +
  340 + for (j1 = 0; j1 < n; j1++) // For each column in sub-matrix get space for the
  341 + { // pointer list
  342 + m = (double **) malloc((n-1) * sizeof(double *));
  343 +
  344 + for (i = 0; i < n-1; i++)
  345 + m[i] = (double *) malloc((n-1)* sizeof(double));
  346 + // i[0][1][2][3] first malloc
  347 + // m -> + + + + space for 4 pointers
  348 + // | | | | j second malloc
  349 + // | | | +-> _ _ _ [0] pointers to
  350 + // | | +----> _ _ _ [1] and memory for
  351 + // | +-------> _ a _ [2] 4 doubles
  352 + // +----------> _ _ _ [3]
  353 + //
  354 + // a[1][2]
  355 + // Build sub-matrix with minor elements excluded
  356 +
  357 + for (i = 1; i < n; i++)
  358 + {
  359 + j2 = 0 ; // Start at first sum-matrix column position
  360 + // Loop to copy source matrix less one column
  361 + for (j = 0; j < n; j++)
  362 + {
  363 + if (j == j1) continue; // Do not copy the minor column element
  364 +
  365 + m[i-1][j2] = a[i][j]; // Copy source element into new sub-matrix
  366 + // i-1 because new sub-matrix is one row
  367 + // (and column) smaller with excluded minors
  368 + j2++; // Move to next sub-matrix column position
  369 + }
  370 + }
  371 +
  372 + det += (double)pow(-1.0, 1.0 + j1 + 1.0) * a[0][j1] * determinant(m, n-1);
  373 + // Sum x raised to y power
  374 + // recursively get determinant of next
  375 + // sub-matrix which is now one
  376 + // row & column smaller
  377 +
  378 + for (i = 0; i < n-1; i++) free(m[i]); // Free the storage allocated to
  379 + // this minor's set of pointers
  380 + free(m); // Free the storage for the original
  381 + // pointer to pointer
  382 + }
  383 + }
  384 +
  385 + return(det);
  386 +}
  387 +
  388 +
  389 +// Function to calculate the cofactor of element (row, col)
  390 +int flow::minor(double **src, double **dest, int row, int col, int order)
  391 +{
  392 + // Indicate which col and row is being copied to dest
  393 + int colCount=0,rowCount=0;
  394 +
  395 + for(int i = 0; i < order; i++)
  396 + {
  397 + if(i != row)
  398 + {
  399 + colCount = 0;
  400 + for(int j = 0; j < order; j++)
  401 + {
  402 + // When j is not the element
  403 + if( j != col )
  404 + {
  405 + dest[rowCount][colCount] = src[i][j];
  406 + colCount++;
  407 + }
  408 + }
  409 +
  410 + rowCount++;
  411 + }
  412 + }
  413 +
  414 + return 1;
  415 +}
0 416 \ No newline at end of file
... ...
tira/biomodels/network.h 0 → 100644
  1 +#ifndef JACK_NETWORK_H
  2 +#define JACK_NETWORK_H
  3 +
  4 +#include "centerline.h"
  5 +#include<stim/visualization/obj.h>
  6 +#include<stim/visualization/swc.h>
  7 +#include <stim/math/circle.h>
  8 +#include <stim/structures/kdtree.cuh>
  9 +
  10 +
  11 +// *****help function*****
  12 +
  13 +template<typename T>
  14 +CUDA_CALLABLE T gaussian(T x, T std = 25.0f) {
  15 + return exp(-x / (2.0f * std * std));
  16 +}
  17 +
  18 +#ifdef __CUDACC__
  19 +template<typename T>
  20 +__global__ void find_metric(T* M, size_t n, T* D, T sigma) {
  21 + size_t x = blockDim.x * blockIdx.x + threadIdx.x;
  22 + if (x >= n) return; // segfault
  23 + M[x] = 1.0f - gaussian<T>(D[x], sigma);
  24 +}
  25 +#endif
  26 +
  27 +namespace stim{
  28 +
  29 + template<typename T>
  30 + class network {
  31 + public:
  32 + // define edge class that extends centerline class with radius information
  33 + class edge : public stim::centerline<T> {
  34 + protected:
  35 + std::vector<T> R; // radius at each point on current edge
  36 +
  37 + public:
  38 + unsigned int v[2]; // unique idx for starting and ending point
  39 + using stim::centerline<T>::d;
  40 + using stim::centerline<T>::C;
  41 +
  42 +
  43 + /// constructors
  44 + // empty constructor
  45 + edge() : stim::centerline<T>() {
  46 + v[0] = UINT_MAX; // set to default value, risky!
  47 + v[1] = UINT_MAX;
  48 + }
  49 +
  50 + // constructor that contructs an edge based on a centerline
  51 + edge(stim::centerline<T> c) : stim::centerline<T>(c) {
  52 + size_t num = c.size();
  53 + R.resize(num);
  54 + }
  55 +
  56 + // constructor that constructs an edge based on a centerline and a list of radii
  57 + edge(stim::centerline<T> c, std::vector<T> r) : stim::centerline<T>(c) {
  58 + R = r; // copy radii
  59 + }
  60 +
  61 + // constructor that constructs an edge based on a list of points and a list of radii
  62 + edge(std::vector<stim::vec3<T> > c, std::vector<T> r) : stim::centerline<T>(c) {
  63 + R = r; // copy radii
  64 + }
  65 +
  66 +
  67 + /// basic operations
  68 + // get radius
  69 + T r(size_t idx) {
  70 + return R[idx];
  71 + }
  72 +
  73 + // set radius
  74 + void set_r(T value) {
  75 + size_t num = R.size();
  76 + for (size_t i = 0; i < num; i++)
  77 + R[i] = value;
  78 + }
  79 + void set_r(size_t idx, T value) {
  80 + R[idx] = value;
  81 + }
  82 + void set_r(std::vector<T> value) {
  83 + size_t num = value.size();
  84 + for (size_t i = 0; i < num; i++)
  85 + R[i] = value[i];
  86 + }
  87 +
  88 + // push back a new radius
  89 + void push_back_r(T value) {
  90 + R.push_back(value);
  91 + }
  92 +
  93 +
  94 + /// vector operation
  95 + // insert a new point and radius at specific location
  96 + void insert(size_t idx, stim::vec3<T> p, T r) {
  97 + centerline<T>::insert(idx, p); // insert a new point on current centerline
  98 +
  99 + R.insert(R.begin() + idx, r); // insert a new radius for that point
  100 + }
  101 +