Skip to content

Commit

Permalink
swc clean up minor bug fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@MKRangers
MKRangers committed Jan 20, 2021
1 parent c15b7e2 commit 8d2d516
Show file tree
Hide file tree
Showing 5 changed files with 74 additions and 48 deletions.
26 changes: 21 additions & 5 deletions hackathon/MK/NeuronStructNavigator/NeuronStructUtilities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -358,7 +358,8 @@ void NeuronStructUtil::somaCleanUp(NeuronTree& inputTree)

void NeuronStructUtil::removeRedunNodes(profiledTree& inputProfiledTree)
{
vector<ptrdiff_t> delLocs;
boost::container::flat_set<int> delLocs;
boost::container::flat_set<int> nodeIDMarked;
for (auto& nodeID : inputProfiledTree.node2childLocMap)
{
const NeuronSWC& currNode = inputProfiledTree.tree.listNeuron.at(inputProfiledTree.node2LocMap.at(nodeID.first));
Expand All @@ -367,15 +368,30 @@ void NeuronStructUtil::removeRedunNodes(profiledTree& inputProfiledTree)
const NeuronSWC& childNode = inputProfiledTree.tree.listNeuron.at(childLoc);
if (childNode.x == currNode.x && childNode.y == currNode.y && childNode.z == currNode.z)
{
delLocs.push_back(childLoc);
for (auto& node : inputProfiledTree.tree.listNeuron)
if (node.parent == childNode.n) node.parent = currNode.n;

delLocs.insert(childLoc);
nodeIDMarked.insert(childNode.n);
if (inputProfiledTree.node2childLocMap.find(childNode.n) != inputProfiledTree.node2childLocMap.end())
for (auto& grandChildLoc : inputProfiledTree.node2childLocMap.at(childNode.n)) inputProfiledTree.tree.listNeuron[grandChildLoc].parent = currNode.n;
{
if (nodeIDMarked.find(currNode.n) == nodeIDMarked.end())
{
for (auto& grandChildLoc : inputProfiledTree.node2childLocMap.at(childNode.n)) inputProfiledTree.tree.listNeuron[grandChildLoc].parent = currNode.n;
}
else
{
int paNodeID = currNode.parent;
while (nodeIDMarked.find(paNodeID) != nodeIDMarked.end()) paNodeID = inputProfiledTree.tree.listNeuron.at(inputProfiledTree.node2LocMap.at(paNodeID)).parent;
for (auto& grandChildLoc : inputProfiledTree.node2childLocMap.at(childNode.n)) inputProfiledTree.tree.listNeuron[grandChildLoc].parent = paNodeID;
}
}
}
}
}

sort(delLocs.rbegin(), delLocs.rend());
for (auto& delLoc : delLocs) inputProfiledTree.tree.listNeuron.erase(inputProfiledTree.tree.listNeuron.begin() + delLoc);
for (boost::container::flat_set<int>::reverse_iterator rit = delLocs.rbegin(); rit != delLocs.rend(); ++rit)
inputProfiledTree.tree.listNeuron.erase(inputProfiledTree.tree.listNeuron.begin() + *rit);
profiledTreeReInit(inputProfiledTree);
}

Expand Down
10 changes: 6 additions & 4 deletions hackathon/MK/NeuronStructNavigator/integratedDataTypes.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -626,7 +626,7 @@ int integratedDataTypes::profiledTree::findNearestSegEndNodeID(const CellAPO inp
vector<int> nodeIDs;
if (this->nodeTileMap.find(targetNodeTileKey) != this->nodeTileMap.end())
{
//cout << targetNodeTileKey << endl;
//cout << "soma node tile: " << targetNodeTileKey << endl;
nodeIDs = this->nodeTileMap.at(targetNodeTileKey);
for (auto& nodeID : nodeIDs)
{
Expand All @@ -647,9 +647,11 @@ int integratedDataTypes::profiledTree::findNearestSegEndNodeID(const CellAPO inp
}
}
}
return outputNodeID;

if (outputNodeID != 0) return outputNodeID;
}
else

if (this->nodeTileMap.find(targetNodeTileKey) == this->nodeTileMap.end() || outputNodeID == 0)
{
for (int k = -1; k <= 1; ++k)
{
Expand All @@ -670,7 +672,7 @@ int integratedDataTypes::profiledTree::findNearestSegEndNodeID(const CellAPO inp
int newYkey = stoi(tileKeyStrings[1]) - j;
int newZkey = stoi(tileKeyStrings[2]) - k;
string newTileKey = to_string(newXkey) + "_" + to_string(newYkey) + "_" + to_string(newZkey);
//cout << newTileKey << endl;
//cout << "surrounding tile: " << newTileKey << endl;
if (this->nodeTileMap.find(newTileKey) != this->nodeTileMap.end())
{
nodeIDs.clear();
Expand Down
51 changes: 30 additions & 21 deletions hackathon/MK/NeuronStructNavigator/integratedDataTypes.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,18 +162,14 @@ namespace integratedDataTypes
NeuronTree tree;
map<int, size_t> node2LocMap;
map<int, vector<size_t>> node2childLocMap;
map<string, vector<int>> nodeTileMap; // tile label -> node ID

NeuronTree sortedTree;
map<string, float> morphFeatureMap;
void getMorphFeatures();
map<string, vector<int>> nodeTileMap; // tile label -> node IDs

/* ================= Segment-related Data Members and Functions ================= */
map<int, segUnit> segs; // key = seg ID
map<string, vector<int>> segHeadMap; // tile label -> seg ID
map<string, vector<int>> segTailMap; // tile label -> seg ID

// -- NOT INITIALIZED WHEN PROFILEDTREE IS CONSTRUCTED -- //
// --------- NOT INITIALIZED WHEN PROFILEDTREE IS CONSTRUCTED --------- //
boost::container::flat_map<int, boost::container::flat_set<int>> segHeadClusters; // key is ordered cluster number label; cluster number -> all seg IDs with heads in the cluster
boost::container::flat_map<int, boost::container::flat_set<int>> segTailClusters; // key is ordered cluster number label; cluster number -> all seg IDs with tails in the cluster
boost::container::flat_map<int, int> headSeg2ClusterMap; // segment ID -> the cluster in which the segment head is located
Expand All @@ -186,31 +182,44 @@ namespace integratedDataTypes
boost::container::flat_map<string, integratedDataTypes::overlappedCoord> overlappedCoordMap; // node coord key -> profiled overlapped coordinate

boost::container::flat_map<int, boost::container::flat_set<vector<float>>> segEndClusterNodeMap; // segEnd cluster ID -> all seg end nodes' coordinates in the cluster
boost::container::flat_map<int, vector<float>> segEndClusterCentroidMap; // segEnd cluster ID -> the coordiate of the centroid of all nodes in the cluster
boost::container::flat_map<int, vector<float>> segEndClusterCentroidMap; // segEnd cluster ID -> the coordiate of the centroid of all nodes in the cluster
boost::container::flat_map<int, vector<segPairProfile>> cluster2segPairMap; // segEnd cluster -> all possible seg pair combinations in the cluster

map<int, segUnit> seg2MiddleBranchingMap; // original segment ID -> rearranged segment with head node in the middle (self branching)
// -------------------------------------------------------------------- //

boost::container::flat_map<int, vector<segPairProfile>> cluster2segPairMap; // segEnd cluster -> all possible seg pair combinations in the cluster
// ------------------------------------------------------ //
void getSegEndClusterNodeMap();
void getSegEndClusterCentoirds();

void nodeSegMapGen();
void nodeCoordKeySegMapGen();
void segEndCoordKeySegMapGen();
void nodeCoordKeyNodeIDmapGen();
void overlappedCoordMapGen();

map<int, segUnit> seg2MiddleBranchingMap;
int findNearestSegEndNodeID(const CellAPO inputAPO);
void assembleSegs2singleTree(int rootNodeID);
void getSegEndClusterNodeMap(); // -> this->segEndClusterNodeMap
void getSegEndClusterCentoirds(); // -> this->segEndClusterCentroidMap

void nodeSegMapGen(); // -> this->node2segMap
void nodeCoordKeySegMapGen(); // -> this->nodeCoordKey2segMap
void segEndCoordKeySegMapGen(); // -> this->segEndCoordKey2segMap
void nodeCoordKeyNodeIDmapGen(); // -> this->nodeCoordKey2nodeIDMap
void overlappedCoordMapGen(); // -> overlappedCoordMap

// ----------------- Segment Assembling Methods ----------------- //
// Pick up the nearst node to the given soma coordinate. The node has to be a singular node that doesn't overlap any other nodes.
// If no node is found within the soma node tile and its 26 surrounding tiles, then a random singular will be picked.
int findNearestSegEndNodeID(const CellAPO inputAPO);

// Clean up and organize the 'mess' of manually traced segments, so that the segments will be "combed through" with a given starting node.
// This method is required before running [this->assembleSegs2singleTree] to transform segments into a single tree.
void combSegs(int rootNodeID);

// This method takes "combed" segments and assembles them to a single tree. (Only 1 segment)
void assembleSegs2singleTree(int rootNodeID);
// -------------------------------------------------------------- //

private:
void nodeSegMapGen(const map<int, segUnit>& segMap, boost::container::flat_map<int, int>& node2segMap);
void nodeCoordKeySegMapGen(const map<int, segUnit>& segMap, boost::container::flat_multimap<string, int>& nodeCoordKey2segMap);
void segEndCoordKeySegMapGen(const map<int, segUnit>& segMap, boost::container::flat_multimap<string, int>& segEndCoordKey2segMap);
void nodeCoordKeyNodeIDmapGen(const QList<NeuronSWC>& nodeList, boost::container::flat_multimap<string, int>& nodeCoordKey2nodeIDmap);

// This method is called by [this->combSegs] to recursively examine every segment and correct the heading direction.
void rc_reverseSegs(const int leadingSegID, const int startingEndNodeID, set<int>& checkedSegIDs);

// This method is called when it needs to be converted to a segment that self-branches from the node to which other segment's end is attached.
pair<int, segUnit> splitSegWithMiddleHead(const segUnit& inputSeg, int newHeadID);
/* ============ END of [Segment-related Data Members and Functions] ============ */

Expand Down
14 changes: 12 additions & 2 deletions hackathon/MK/swcRename/ReconOperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,11 +116,13 @@ void ReconOperator::removeDupedNodes()
QStringList fileList = inputFolder.entryList();

QString outputFolderQ = this->rootPath + "\\allCleanedUp\\";
QString remainingFolderQ = this->rootPath + "\\stillMoreThan1Root\\";
QDir outputDir(outputFolderQ);
QDir remainingDir(remainingFolderQ);
if (!outputDir.exists()) outputDir.mkpath(".");

/*QString subTreeFolderQ = this->rootPath + "\\subTrees\\";
QDir subTreeDir(subTreeFolderQ);
if (!subTreeDir.exists()) subTreeDir.mkpath(".");*/

for (auto& file : fileList)
{
QString baseName;
Expand Down Expand Up @@ -150,13 +152,20 @@ void ReconOperator::removeDupedNodes()
NeuronStructUtil::removeRedunNodes(inputProfiledTree);
if (NeuronStructUtil::multipleSegsCheck(inputTree))
{
/*QString treesFolderQ = subTreeFolderQ + baseName + "\\";
QDir treesFolderDir(treesFolderQ);
if (!treesFolderDir.exists()) treesFolderDir.mkpath(".");*/

clock_t start = clock();
boost::container::flat_map<int, profiledTree> connectedTrees = NeuronStructExplorer::groupGeoConnectedTrees(inputProfiledTree.tree);
cout << endl << "-- " << connectedTrees.size() << " separate trees identified." << endl << endl;
map<int, int> tree2HeadNodeMap;
int minNodeID, maxNodeID;
for (auto& connectedTree : connectedTrees)
{
//writeSWC_file(treesFolderQ + QString::number(int(connectedTrees.find(connectedTree.first) - connectedTrees.begin()) + 1) + ".swc", connectedTree.second.tree);

NeuronStructUtil::removeRedunNodes(inputProfiledTree);
cout << "Processing tree " << int(connectedTrees.find(connectedTree.first) - connectedTrees.begin()) + 1 << "..." << endl;
minNodeID = 10000000;
maxNodeID = 0;
Expand All @@ -172,6 +181,7 @@ void ReconOperator::removeDupedNodes()
connectedTree.second.assembleSegs2singleTree(nearestNodeID);
tree2HeadNodeMap.insert({ connectedTree.first, nearestNodeID });
connectedTree.second.tree.listNeuron[connectedTree.second.node2LocMap.at(nearestNodeID)].type = 1;
cout << "-----> Finish with tree " << int(connectedTrees.find(connectedTree.first) - connectedTrees.begin()) + 1 << endl << endl;
}
else
{
Expand Down
21 changes: 5 additions & 16 deletions hackathon/MK/tester/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ int main(int argc, char* argv[])
paras.push_back(paraString);
}

//string funcName = "groupTree";
//string funcName = "removeRedunNode";
/************************************/

ImgTester myImgTester;
Expand Down Expand Up @@ -71,9 +71,12 @@ int main(int argc, char* argv[])
else if (!funcName.compare("removeRedunNode"))
{
NeuronTree inputTree = readSWC_file(QString::fromStdString(paras.at(0)));
//QString inputFileNameQ = "C:\\Users\\hkuo9\\Desktop\\the_cell_for_test_preprocess_\\dupSegs_composite4\\191812_3801-X5383-Y20989_finalized.ano.eswc";
//NeuronTree inputTree = readSWC_file(inputFileNameQ);
profiledTree inputProfiledTree(inputTree);
NeuronStructUtil::removeRedunNodes(inputProfiledTree);
writeSWC_file(QString::fromStdString(paras.at(0)) + "_redunRemoved.swc", inputProfiledTree.tree);
//writeSWC_file("C:\\Users\\hkuo9\\Desktop\\the_cell_for_test_preprocess_\\dupSegs_composite4\\191812_3801-X5383-Y20989_finalized.ano.eswc_dupRemoved.swc", inputProfiledTree.tree);
}
else if (!funcName.compare("dupNodeRemove"))
{
Expand All @@ -98,16 +101,8 @@ int main(int argc, char* argv[])
clock_t start = clock();
//profiledTree inputProfiledTree(readSWC_file("C:\\Users\\hkuo9\\Desktop\\test1\\" + file));
profiledTree inputProfiledTree(readSWC_file(QString::fromStdString(paras.at(0)) + "\\" + file));
/*for (auto& seg : inputProfiledTree.segs)
{
cout << seg.first << ": " << endl;
for (auto& node : seg.second.nodes) cout << "(" << node.x << ", " << node.y << ", " << node.z << ")" << " ";
cout << endl;
}
cout << endl;*/
//inputProfiledTree.segEndCoordKeySegMapGen();
//for (auto& coordKey : inputProfiledTree.segEndCoordKey2segMap) cout << coordKey.first << " " << coordKey.second << endl;
boost::container::flat_map<int, profiledTree> groupedProfiledTrees = NeuronStructExplorer::groupGeoConnectedTrees(readSWC_file(QString::fromStdString(paras.at(0)) + "\\" + file));
//boost::container::flat_map<int, profiledTree> groupedProfiledTrees = NeuronStructExplorer::groupGeoConnectedTrees(inputProfiledTree1.tree);
//vector<profiledTree> groupedProfiledTrees = NeuronStructExplorer::groupGeoConnectedTrees(readSWC_file("C:\\Users\\hkuo9\\Desktop\\test1\\" + file));
clock_t end = clock();

Expand All @@ -119,12 +114,6 @@ int main(int argc, char* argv[])
QDir outputDir(outputFolderQ);
if (!outputDir.exists()) outputDir.mkpath(".");

/*vector<>
for (boost::container::flat_map<int, profiledTree>::reverse_iterator rit = groupedProfiledTrees.rbegin(); rit != groupedProfiledTrees.rend(); ++rit)
{
}*/

for (boost::container::flat_map<int, profiledTree>::iterator it = groupedProfiledTrees.begin(); it != groupedProfiledTrees.end(); ++it)
{
QString treeNameQ = outputFolderQ + "\\" + QString::number(int(it - groupedProfiledTrees.begin() + 1)) + ".swc";
Expand Down

0 comments on commit 8d2d516

Please sign in to comment.