diff --git a/stim/visualization/network.h b/stim/visualization/network.h index 9711fd6..528b91f 100644 --- a/stim/visualization/network.h +++ b/stim/visualization/network.h @@ -1,13 +1,20 @@ #ifndef STIM_NETWORK_H #define STIM_NETWORK_H +#include +#include +#include +#include +#include +#include +#include #include #include -#include +#include #include -#include "fiber.h" -#include +#include +#define multfactor 2.51 namespace stim{ /** This is the a class that interfaces with gl_spider in order to store the currently @@ -77,10 +84,7 @@ class network{ setEndpoints(int id1, int id2) { Nodes[1] = id1; Nodes[2] = id2; - } - - - + } }; @@ -480,40 +484,74 @@ class network{ ///@param pos: std::vector of stim vecs with the positions of the point in the fiber. ///@param radii: std::vector of floats with the radii of the fiber at positions in pos. - ///adds an edge from two std::vectors with positions and radii. - void - addEdge(std::vector< stim::vec > pos, std::vector radii) + ///returns the forest as a std::string. For testing only. + std::string + edges_to_str() { - edge *a = new edge(pos,radii); - E.push_back(a); + std::stringstream ss; + for(unsigned int i = 0; i < E.size(); i++) + { + ss << i << ": " << E[i]->str() << std::endl; + } + return(ss.str()); } - void - addNode(stim::vec nodes) + // total number of points on all edges!=fibers in the network + int + totalEdges() { - V.push_back(nodes); + int totalEdgeVertices=0;int N=0; + for (int i=0; i < sizeE(); i ++) + {// FOR N points on the fiber N-1 edges are possible + N = E[i]->n_pts(); + totalEdgeVertices = totalEdgeVertices + N- 1; + } + return totalEdgeVertices; } - - ///adds an edge from two std::lists with positions and radii. - /// NOT NECESSARY. - void - addEdge(std::list< stim::vec > pos, std::list radii) + // sum of all the fiber lengths + float + lengthOfNetwork() { - edge a = edge(pos,radii); - E.push_back(a); + float networkLength=0;float N=0; + for (int i=0; i < sizeE(); i ++) + { + N = E[i]->length(); + networkLength = networkLength + N; + } + return networkLength; } - - ///returns the forest as a std::string. For testing only. + /// get index of a node on a fiber + // by matching the node on fiber to already set vertices (both strings) + // used in obj file conversion + int + getIndexes(std::string* input, std::string searched, int sizeV) + { + int result = 0; + for (int i = 0; i < sizeV; i++) + { + if (input[i] == searched) + { + result = i + 1; + } + } + return result; + } + // strObj returns a string of fiber indices corresponding to vectors of positions in the fiber including intermediate nodes std::string - edges_to_str() + strObj(std::string* strArray, int sizeV) { std::stringstream ss; - for(unsigned int i = 0; i < E.size(); i++) + std::stringstream oss; + for(unsigned int i = 0; i < N; i++) { - ss << i << ": " << E[i]->str() << std::endl; + ss.str(std::string()); + for(unsigned int d = 0; d < 3; d++) + { + ss< > @@ -541,8 +579,158 @@ class network{ ss << "[from Node " << E[i] -> Nodes[1] << " to " << E[i] -> Nodes[2] << "]"; return ss.str(); } - + //load an obj file into a network class + stim::network + objToNetwork(stim::obj objInput) + { + stim::network nwc; + // function to convert an obj file loaded using stim/visualization/obj.h + // to a 3D network class using network.h methods. + std::vector< stim::vec > fiberPositions; //initialize positions on the fiber to be added to network + objInput.getLine(1, fiberPositions); int numDim = fiberPositions[0].size();//find dimensions of the vertices in obj file + // to verify if the nodes are already pushed into node list + std::vector validList; + validList.assign(objInput.numV(), false); + // go through each of the lines "l followed by indices in obj and add all start and end vertices as the nodes + // using addNode function for adding nodes + // and the line as an edge(belongs to fiber class) using addEdge function + std::vector vertexIndices(objInput.numV()); + std::vector< stim::vec > allVerticesList = objInput.getAllV(vertexIndices); + for (unsigned i =1; i< objInput.numL() + 1; i++) + { + // edges/fibers could be added to the network class + std::vector< stim::vec > fiberPositions; + + objInput.getLine(i, fiberPositions); + // finding size to allocate radii + int numPointsOnFiber = fiberPositions.size(); + // to extend it to a 3D postion if it is a 1D/2D vertex in space + std::vector< stim::vec > fiberPositions1(numPointsOnFiber); + // 2D case append and assign the last dimension to zero + if (numDim == 2) + {// 2D case append and assign the last dimension to zero repeating it for all the points on fiber + for (int j = 0; j < numPointsOnFiber; j ++) + { + fiberPositions1[j][numDim - 2] = fiberPositions[j][0]; + fiberPositions1[j][numDim -1 ] = fiberPositions[j][1]; + fiberPositions1[j][numDim] = 0; + } + } + // 1D case append and assign last two dimensions to zero + else if (numDim == 1) + { + for (int j = 0; j < numPointsOnFiber; j ++) + { + fiberPositions1[j][numDim - 2] = fiberPositions[j][0]; + fiberPositions1[j][numDim -1 ] = 0; + fiberPositions1[j][numDim] = 0; + } + } + else + { + fiberPositions1 = fiberPositions; + } + std::vector radii(numPointsOnFiber); // allocating radii to zero + // then add edge + // edge* a = new edge(fiberPositions1,radii); + //E.push_back(a); + nwc.addEdge(fiberPositions1, radii); + // add nodes + stim::vec v0(3);stim::vec vN(3); + int endId = numPointsOnFiber -1; + v0[0] = fiberPositions1[0][0];v0[1] = fiberPositions1[0][1];v0[2] = fiberPositions1[0][2]; + vN[0] = fiberPositions1[endId][0];vN[1] = fiberPositions1[endId][1];vN[2] = fiberPositions1[endId][2]; + // VISITED INDEXES OF the nodes are set to true + if(!validList[objInput.getIndex(allVerticesList, v0)]) + { + //V.push_back(node(v0)); + nwc.addNode(v0); + validList[objInput.getIndex(allVerticesList, v0)] = true; + } + if(!validList[objInput.getIndex(allVerticesList, vN)]) + { + //V.push_back(node(vN)); + nwc.addNode(vN); + validList[objInput.getIndex(allVerticesList, vN)] = true; + } + } + return nwc; + } + // convert ground and test case to network ,kd trees + boost::tuple< ANNkd_tree*, ANNkd_tree*, stim::network, stim::network > + LoadNetworks(std::string gold_filename, std::string test_filename) + { + ANNkd_tree* kdTree1;ANNkd_tree* kdTree2; + using namespace stim; + network network1;network network2; + stim::obj objFile1(gold_filename); + network1 = objToNetwork(objFile1); + kdTree1 = generateKdTreeFromNetwork(network1); + std::cout<<"Gold Standard:network and kdtree generated"< objFile2(test_filename); + network2 = objToNetwork(objFile2); + kdTree2 = generateKdTreeFromNetwork(network2); + std::cout<<"Test Network:network and kdtree generated"< p0, std::vector p1) + { + double sum = 0; + for(unsigned int d = 0; d < 3; d++) + sum += p0[d] * p1[d]; + return sqrt(sum); + } + // sum of elements in a vector + double sum(stim::vec metricList) + { + float sumMetricList = 0; + for (int count=0; count network) //kd-tree stores all points in the fiber for fast searching + { + ANNkd_tree* kdt; + double **c; //centerline (array of double pointers) + float* r; // array of fiber radii + unsigned int n_data = network.sizeV(); //set the number of points + c = (double**) malloc(sizeof(double*) * n_data); //allocate the array pointer + for(unsigned int i = 0; i < n_data; i++) //allocate space for each point + {c[i] = (double*) malloc(sizeof(double) * 3);} + r = (float*) malloc(sizeof(float) * n_data); //allocate space for the radii + stim::vec node; + for(int i = 0; i < n_data; i++) + { + node = network.V[i].getPosition(); + //convert_to_double + for(unsigned int d = 0; d < 3; d++){ //for each dimension + c[i][d] = double(node[d]); //copy the coordinate + } + r[i] = r[i]; //copy the radius + } + ANNpointArray pts = (ANNpointArray)c; //create an array of data points + kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree + return kdt; + } + + ///@param pos: std::vector of stim vecs with the positions of the point in the fiber. + ///@param radii: std::vector of floats with the radii of the fiber at positions in pos. + ///adds an edge from two std::vectors with positions and radii. + void + addEdge(std::vector< stim::vec > pos, std::vector radii) + { + edge *a = new edge(pos,radii); + E.push_back(a); + } + void + addNode(stim::vec nodes) + { + V.push_back(nodes); + } ///exports the graph. void to_csv() @@ -575,6 +763,65 @@ class network{ } ofs.close(); } + // gaussian function + float gaussianFunction(float x, float std=25) + { + float normalization = multfactor * std; + float evaluate = 1.0 - (exp(-x/(2*std*std))); + return evaluate; + } + boost::tuple< float, float > + compareSkeletons(boost::tuple< ANNkd_tree*, ANNkd_tree*, stim::network, stim::network > networkKdtree) + { + float gFPR, gFNR; + gFPR = CompareNetKD(networkKdtree.get<0>(), networkKdtree.get<3>()); + gFNR = CompareNetKD(networkKdtree.get<1>(), networkKdtree.get<2>()); + return boost::make_tuple(gFPR, gFNR); + } + // gaussian distance of points on network to Kdtree + float + CompareNetKD(ANNkd_tree* kdTreeGroundtruth, stim::network networkTruth) + { + std::vector< stim::vec > fiberPoints; + float netmetsMetric=0; + double eps = 0; // error bound + ANNdistArray dists1;dists1 = new ANNdist[1]; // near neighbor distances + ANNdistArray dists2; dists2 = new ANNdist[1];// near neighbor distances + ANNidxArray nnIdx1; nnIdx1 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices + ANNidxArray nnIdx2; nnIdx2 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices + int N; int numQueryPoints = networkTruth.totalEdges(); + float totalNetworkLength = networkTruth.lengthOfNetwork(); + stim::vec fiberMetric(networkTruth.sizeE()); + //for each fiber + for (unsigned i=0; i < networkTruth.sizeE(); i ++) + { + std::vector p1(3); std::vector p2(3);//temporary variables to store point positions + fiberPoints = networkTruth.E[i]->centerline(); + N = networkTruth.E[i]->n_pts(); + stim::vec segmentMetric(N-1); + // for each segment in the fiber + for (unsigned j = 0; j < N - 1; j++) + { + ANNpoint queryPt1; queryPt1 = annAllocPt(3); + ANNpoint queryPt2; queryPt2 = annAllocPt(3); + //for each dimension of the point + for(unsigned int d = 0; d < 3; d++) + { + queryPt1[d] = double(fiberPoints[j][d]); + p1[d] = double(fiberPoints[j][d]); + queryPt2[d] = double(fiberPoints[j + 1][d]); + p2[d] = double(fiberPoints[j + 1][d]);// for each dimension + } + kdTreeGroundtruth->annkSearch( queryPt1, 1, nnIdx1, dists1, eps); // error bound + kdTreeGroundtruth->annkSearch( queryPt2, 1, nnIdx2, dists2, eps); // error bound + float dist1 = gaussianFunction(dists1[0]);float dist2 = gaussianFunction(dists2[0]); + segmentMetric[j] = (((dist1 + dist2) / 2) * dist(p1, p2)) ; + } + fiberMetric[i] = sum(segmentMetric); + } + netmetsMetric = sum(fiberMetric)/totalNetworkLength ; + return netmetsMetric; + } void removeCharsFromString(std::string &str, char* charsToRemove ) { for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) { str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() ); @@ -588,10 +835,10 @@ class network{ ofs.open("Graph.obj", std::ofstream::out | std::ofstream::app); int num = V.size(); string *strArray = new string[num]; - for(int i = 0; i < V.size(); i++) + for(int i = 0; i < allVerticesList.size(); i++) { std::string str; - str = V[i].str(); + str = allVerticesList[i].str(); removeCharsFromString(str, "[],"); ofs << "v " << str << "\n"; removeCharsFromString(str," "); -- libgit2 0.21.4