Skip to content

Commit

Permalink
NeuronRecon_Manager update
Browse files Browse the repository at this point in the history
  • Loading branch information
@MKRangers
MKRangers committed Jan 17, 2021
1 parent 79e4ed2 commit c15b7e2
Show file tree
Hide file tree
Showing 8 changed files with 232 additions and 58 deletions.
2 changes: 1 addition & 1 deletion hackathon/MK/NeuronStructNavigator/NeuronStructExplorer.cpp
View file
Expand Up @@ -1160,7 +1160,7 @@ boost::container::flat_map<int, profiledTree> NeuronStructExplorer::groupGeoConn
{
set<int> groupedSegIDs;
groupedSegIDs.insert(*unGroupedSegIDs.begin());
cout << *unGroupedSegIDs.begin() << endl;
//cout << *unGroupedSegIDs.begin() << endl;
NeuronStructExplorer::rc_findConnectedSegs(inputProfiledTree, groupedSegIDs, *unGroupedSegIDs.begin());

NeuronTree outputTree;
Expand Down
10 changes: 10 additions & 0 deletions hackathon/MK/NeuronStructNavigator/NeuronStructExplorer.h
View file
Expand Up @@ -52,6 +52,7 @@ class NeuronStructExplorer
static inline string getNodeTileKey_noZratio(const NeuronSWC& inputNode, float nodeTileLength = NODE_TILE_LENGTH);
static inline string getNodeTileKey(const float nodeCoords[], float nodeTileLength = NODE_TILE_LENGTH);
static inline string getNodeTileKey(const ImageMarker& inputMarker, float nodeTileLength = NODE_TILE_LENGTH);
static inline string getNodeTileKey(const CellAPO& inputAPO, float nodeTileLength = NODE_TILE_LENGTH);

// Returns the corresponding string segment tile key with the given node or marker.
static inline string getSegTileKey(const NeuronSWC& inputNode, float segTileLength = SEGtileXY_LENGTH);
Expand Down Expand Up @@ -213,6 +214,15 @@ inline string NeuronStructExplorer::getNodeTileKey(const ImageMarker& inputMarke
return keyLabel;
}

inline string NeuronStructExplorer::getNodeTileKey(const CellAPO& inputAPO, float nodeTileLength)
{
string xLabel = to_string(int(inputAPO.x / nodeTileLength));
string yLabel = to_string(int(inputAPO.y / nodeTileLength));
string zLabel = to_string(int(inputAPO.z / (nodeTileLength / zRATIO)));
string keyLabel = xLabel + "_" + yLabel + "_" + zLabel;
return keyLabel;
}

inline string NeuronStructExplorer::getSegTileKey(const NeuronSWC& inputNode, float segTileLength)
{
string xLabel = to_string(int(inputNode.x / segTileLength));
Expand Down
143 changes: 129 additions & 14 deletions hackathon/MK/NeuronStructNavigator/integratedDataTypes.cpp
View file
Expand Up @@ -613,20 +613,118 @@ void integratedDataTypes::profiledTree::nodeTileResize(float nodeTileLength)
}
}

int end2bodyCount = 0;
void integratedDataTypes::profiledTree::combSegs(int rootNodeID)
int integratedDataTypes::profiledTree::findNearestSegEndNodeID(const CellAPO inputAPO)
{
if (this->node2segMap.empty()) this->nodeSegMapGen();
if (this->nodeCoordKey2segMap.empty()) this->nodeCoordKeySegMapGen();
if (this->segEndCoordKey2segMap.empty()) this->segEndCoordKeySegMapGen();
if (this->nodeCoordKey2nodeIDMap.empty()) this->nodeCoordKeyNodeIDmapGen();
int leadingSegID = this->node2segMap.at(rootNodeID);
set<int> checkedSegIDs = { this->node2segMap.at(rootNodeID) };
this->rc_reverseSegs(leadingSegID, rootNodeID, checkedSegIDs);
this->segEndClusterNodeMap.clear();
this->segEndCoordKeySegMapGen();
this->nodeCoordKey2segMap.clear();
this->nodeCoordKeySegMapGen();

cout << "original segment number: " << this->segs.size() << endl;
cout << "checked segment number: " << checkedSegIDs.size() << endl;
cout << "total end to body count: " << end2bodyCount << endl;
string targetNodeTileKey = NeuronStructExplorer::getNodeTileKey(inputAPO);
float dist = 100000;
int outputNodeID = 0;
vector<int> nodeIDs;
if (this->nodeTileMap.find(targetNodeTileKey) != this->nodeTileMap.end())
{
//cout << targetNodeTileKey << endl;
nodeIDs = this->nodeTileMap.at(targetNodeTileKey);
for (auto& nodeID : nodeIDs)
{
const NeuronSWC& node = this->tree.listNeuron.at(this->node2LocMap.at(nodeID));
string nodeCoordKey = to_string(node.x) + "_" + to_string(node.y) + "_" + to_string(node.z);
if (this->segEndCoordKey2segMap.find(nodeCoordKey) != this->segEndCoordKey2segMap.end())
{
pair<boost::container::flat_multimap<string, int>::iterator, boost::container::flat_multimap<string, int>::iterator> range = this->segEndCoordKey2segMap.equal_range(nodeCoordKey);
pair<boost::container::flat_multimap<string, int>::iterator, boost::container::flat_multimap<string, int>::iterator> range2 = this->nodeCoordKey2segMap.equal_range(nodeCoordKey);
if (range.second - range.first == 1 && range2.second - range2.first == 1)
{
float currNodeDist = sqrtf((node.x - inputAPO.x) * (node.x - inputAPO.x) + (node.y - inputAPO.y) * (node.y - inputAPO.y) + (node.z - inputAPO.z) * (node.z - inputAPO.z));
if (currNodeDist < dist)
{
outputNodeID = nodeID;
dist = currNodeDist;
}
}
}
}
return outputNodeID;
}
else
{
for (int k = -1; k <= 1; ++k)
{
for (int j = -1; j <= 1; ++j)
{
for (int i = -1; i <= 1; ++i)
{
vector<string> tileKeyStrings;
stringstream ss(targetNodeTileKey);
while (ss.good())
{
string subStr;
getline(ss, subStr, '_');
tileKeyStrings.push_back(subStr);
}

int newXkey = stoi(tileKeyStrings[0]) - i;
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;
if (this->nodeTileMap.find(newTileKey) != this->nodeTileMap.end())
{
nodeIDs.clear();
nodeIDs = this->nodeTileMap.at(newTileKey);
for (auto& nodeID : nodeIDs)
{
const NeuronSWC& node = this->tree.listNeuron.at(this->node2LocMap.at(nodeID));
string nodeCoordKey = to_string(node.x) + "_" + to_string(node.y) + "_" + to_string(node.z);
if (this->segEndCoordKey2segMap.find(nodeCoordKey) != this->segEndCoordKey2segMap.end())
{
pair<boost::container::flat_multimap<string, int>::iterator, boost::container::flat_multimap<string, int>::iterator> range = this->segEndCoordKey2segMap.equal_range(nodeCoordKey);
pair<boost::container::flat_multimap<string, int>::iterator, boost::container::flat_multimap<string, int>::iterator> range2 = this->nodeCoordKey2segMap.equal_range(nodeCoordKey);
if (range.second - range.first == 1 && range2.second - range2.first == 1)
{
float currNodeDist = sqrtf((node.x - inputAPO.x) * (node.x - inputAPO.x) + (node.y - inputAPO.y) * (node.y - inputAPO.y) + (node.z - inputAPO.z) * (node.z - inputAPO.z));
if (currNodeDist < dist)
{
outputNodeID = nodeID;
dist = currNodeDist;
}
}
}
}
}
}
}
}

if (outputNodeID != 0) return outputNodeID;
else
{
for (auto& segEndCoord : this->segEndCoordKey2segMap)
{
pair<boost::container::flat_multimap<string, int>::iterator, boost::container::flat_multimap<string, int>::iterator> range = this->segEndCoordKey2segMap.equal_range(segEndCoord.first);
pair<boost::container::flat_multimap<string, int>::iterator, boost::container::flat_multimap<string, int>::iterator> range2 = this->nodeCoordKey2segMap.equal_range(segEndCoord.first);
if (range.second - range.first == 1 && range2.second - range2.first == 1)
{
const segUnit& currSeg = this->segs.at(segEndCoord.second);
const NeuronSWC& currHeadNode = currSeg.nodes.at(currSeg.seg_nodeLocMap.at(currSeg.head));
string headCoordKey = to_string(currHeadNode.x) + "_" + to_string(currHeadNode.y) + "_" + to_string(currHeadNode.z);
if (!headCoordKey.compare(segEndCoord.first)) return currSeg.head;
else
{
for (auto& tailID : currSeg.tails)
{
const NeuronSWC& tailNode = currSeg.nodes.at(currSeg.seg_nodeLocMap.at(tailID));
string tailCoordKey = to_string(tailNode.x) + "_" + to_string(tailNode.y) + "_" + to_string(tailNode.z);
if (!tailCoordKey.compare(segEndCoord.first)) return tailID;
}
}
}
}
}
}
}

void integratedDataTypes::profiledTree::assembleSegs2singleTree(int rootNodeID)
Expand All @@ -642,9 +740,26 @@ void integratedDataTypes::profiledTree::assembleSegs2singleTree(int rootNodeID)
}
NeuronStructUtil::removeDupHeads(this->tree);
profiledTreeReInit(*this);
cout << " --> segment number: " << this->segs.size() << endl;
}
}

//int end2bodyCount = 0;
void integratedDataTypes::profiledTree::combSegs(int rootNodeID)
{
if (this->node2segMap.empty()) this->nodeSegMapGen();
if (this->nodeCoordKey2segMap.empty()) this->nodeCoordKeySegMapGen();
if (this->segEndCoordKey2segMap.empty()) this->segEndCoordKeySegMapGen();
if (this->nodeCoordKey2nodeIDMap.empty()) this->nodeCoordKeyNodeIDmapGen();
int leadingSegID = this->node2segMap.at(rootNodeID);
set<int> checkedSegIDs = { this->node2segMap.at(rootNodeID) };
this->rc_reverseSegs(leadingSegID, rootNodeID, checkedSegIDs);

cout << "original segment number: " << this->segs.size() << endl;
cout << "checked segment number: " << checkedSegIDs.size() << endl;
//cout << "total end to body count: " << end2bodyCount << endl;
}

void integratedDataTypes::profiledTree::rc_reverseSegs(const int leadingSegID, const int startingEndNodeID, set<int>& checkedSegIDs)
{
checkedSegIDs.insert(leadingSegID);
Expand Down Expand Up @@ -691,7 +806,7 @@ void integratedDataTypes::profiledTree::rc_reverseSegs(const int leadingSegID, c
cout << "---------" << endl;
cout << "head to body: " << headNode.n << " -> " << it->second << endl;
cout << "---------" << endl << endl;
++end2bodyCount;
//++end2bodyCount;
this->rc_reverseSegs(this->node2segMap.at(it->second), it->second, checkedSegIDs);
}
}
Expand All @@ -708,7 +823,7 @@ void integratedDataTypes::profiledTree::rc_reverseSegs(const int leadingSegID, c
{
cout << "---------" << endl;
cout << "tail to body: " << tailNode.n << " -> " << it->second << endl;
++end2bodyCount;
//++end2bodyCount;
this->seg2MiddleBranchingMap.insert(this->splitSegWithMiddleHead(this->segs.at(this->node2segMap.at(it->second)), it->second));
cout << "---------" << endl << endl;
this->segs[this->node2segMap.at(it->second)].to_be_deleted = true;
Expand Down
6 changes: 3 additions & 3 deletions hackathon/MK/NeuronStructNavigator/integratedDataTypes.h
View file
Expand Up @@ -199,11 +199,11 @@ namespace integratedDataTypes
void nodeCoordKeyNodeIDmapGen();
void overlappedCoordMapGen();

set<set<int>> groupedTreeSegs; // seg IDs of each geometrically grouped trees
map<int, segUnit> seg2MiddleBranchingMap;
void combSegs(int rootNodeID);
int findNearestSegEndNodeID(const CellAPO inputAPO);
void assembleSegs2singleTree(int rootNodeID);

void combSegs(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);
Expand Down
108 changes: 85 additions & 23 deletions hackathon/MK/swcRename/ReconOperator.cpp
View file
Expand Up @@ -101,7 +101,7 @@ void ReconOperator::denAxonCombine(bool dupRemove)
trees.push_back(readSWC_file(axonFullName));
NeuronTree outputTree = NeuronStructUtil::swcCombine(trees);

if (dupRemove) outputTree = NeuronStructUtil::removeDupNodes(outputTree);
//if (dupRemove) outputTree = NeuronStructUtil::removeDupNodes(outputTree);
QString outputSWCfullName = outputFolderName + axonFile;
writeSWC_file(outputSWCfullName, outputTree);
}
Expand All @@ -113,41 +113,103 @@ void ReconOperator::removeDupedNodes()
{
QDir inputFolder(this->rootPath);
inputFolder.setFilter(QDir::Files | QDir::NoDotAndDotDot);
QStringList swcFileList = inputFolder.entryList();
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(".");

for (auto& file : swcFileList)
for (auto& file : fileList)
{
QString inputSWCFullName = this->rootPath + "\\" + file;
NeuronTree inputTree = readSWC_file(inputSWCFullName);
if (NeuronStructUtil::multipleSegsCheck(inputTree))
QString baseName;
if (file.endsWith(".swc")) baseName = file.left(file.length() - 4);
else if (file.endsWith(".eswc")) baseName = file.left(file.length() - 5);

bool apoFound = false;
CellAPO somaAPO;
QString targetAPOfileName = baseName + ".apo";
for (auto& file2 : fileList)
{
QString outputSWCfullName;
NeuronTree dupRemovedTree = NeuronStructUtil::removeDupNodes(inputTree);
if (!NeuronStructUtil::multipleSegsCheck(dupRemovedTree))
if (file2 == targetAPOfileName)
{
if (file.endsWith("eswc")) outputSWCfullName = outputFolderQ + file.left(file.length() - 4) + "swc";
else outputSWCfullName = outputFolderQ + file;
writeSWC_file(outputSWCfullName, dupRemovedTree);
}
else
{
if (file.endsWith("eswc")) outputSWCfullName = remainingFolderQ + file.left(file.length() - 4) + "swc";
else outputSWCfullName = remainingFolderQ + file;
if (!remainingDir.exists()) remainingDir.mkpath(".");
writeSWC_file(outputSWCfullName, dupRemovedTree);
somaAPO = *readAPO_file(this->rootPath + "\\" + file2).begin();
apoFound = true;
break;
}
}
else

if (apoFound)
{
QString oldName = this->rootPath + "\\" + file;
QString newName = outputFolderQ + file;
QFile::copy(oldName, newName);
cout << "Input soma coordinate: " << somaAPO.x << " " << somaAPO.y << " " << somaAPO.z << endl;
QString inputSWCFullName = this->rootPath + "\\" + file;
NeuronTree inputTree = readSWC_file(inputSWCFullName);
NeuronTree noDupSegTree = NeuronStructUtil::removeDupSegs(inputTree);
profiledTree inputProfiledTree(noDupSegTree);
NeuronStructUtil::removeRedunNodes(inputProfiledTree);
if (NeuronStructUtil::multipleSegsCheck(inputTree))
{
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)
{
cout << "Processing tree " << int(connectedTrees.find(connectedTree.first) - connectedTrees.begin()) + 1 << "..." << endl;
minNodeID = 10000000;
maxNodeID = 0;
if (connectedTree.second.node2LocMap.begin()->first < minNodeID) minNodeID = connectedTree.second.node2LocMap.begin()->first;
if (connectedTree.second.node2LocMap.rbegin()->first > maxNodeID) maxNodeID = connectedTree.second.node2LocMap.rbegin()->first;

int nearestNodeID = connectedTree.second.findNearestSegEndNodeID(somaAPO);
cout << "ID of the nearest node to soma coordinates " << nearestNodeID << " -> ";
const NeuronSWC& nearestNode = connectedTree.second.tree.listNeuron.at(connectedTree.second.node2LocMap.at(nearestNodeID));
float dist = sqrtf((nearestNode.x - somaAPO.x) * (nearestNode.x - somaAPO.x) + (nearestNode.y - somaAPO.y) * (nearestNode.y - somaAPO.y) + (nearestNode.z - somaAPO.z) * (nearestNode.z - somaAPO.z));
if (dist <= 35)
{
connectedTree.second.assembleSegs2singleTree(nearestNodeID);
tree2HeadNodeMap.insert({ connectedTree.first, nearestNodeID });
connectedTree.second.tree.listNeuron[connectedTree.second.node2LocMap.at(nearestNodeID)].type = 1;
}
else
{
cout << "Tree " << int(connectedTrees.find(connectedTree.first) - connectedTrees.begin()) + 1 << " Not in the soma range." << endl;
if (connectedTree.second.segs.size() > 1)
{
cout << " More than 1 segments are identified. Still needs to assemble segments." << endl;
connectedTree.second.assembleSegs2singleTree(nearestNodeID);
}
cout << " Change node type to 7." << endl;
for (auto& node : connectedTree.second.tree.listNeuron) node.type = 7;
cout << "-----> Finish with tree " << int(connectedTrees.find(connectedTree.first) - connectedTrees.begin()) + 1 << endl << endl;
}
}
clock_t end = clock();
float duration = float(end - start) / CLOCKS_PER_SEC;
cout << "--> All trees processed. " << duration << " seconds elapsed." << endl;

NeuronSWC somaNode;
if (minNodeID > 1) somaNode.n = minNodeID - 1;
else somaNode.n = maxNodeID + 1;

somaNode.x = somaAPO.x;
somaNode.y = somaAPO.y;
somaNode.z = somaAPO.z;
somaNode.parent = -1;
somaNode.type = 1;
for (auto& treeID : tree2HeadNodeMap)
{
profiledTree& currTree = connectedTrees[treeID.first];
currTree.tree.listNeuron[currTree.node2LocMap.at(treeID.second)].parent = somaNode.n;
}

NeuronTree outputTree;
for (auto& connectedTree : connectedTrees) outputTree.listNeuron.append(connectedTree.second.tree.listNeuron);
outputTree.listNeuron.push_front(somaNode);
writeSWC_file(outputFolderQ + baseName + ".swc", outputTree);
}
}
}
}
1 change: 1 addition & 0 deletions hackathon/MK/swcRename/ReconOperator.h
View file
Expand Up @@ -6,6 +6,7 @@
#include <qfile.h>

#include "NeuronStructUtilities.h"
#include "NeuronStructExplorer.h"

class ReconOperator
{
Expand Down
2 changes: 1 addition & 1 deletion hackathon/MK/swcRename/renameSWC.ui
View file
Expand Up @@ -823,7 +823,7 @@
</size>
</property>
<property name="title">
<string>Remove duplicated nodes:</string>
<string>Assemble segments to single tree:</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout_11">
<item>
Expand Down

0 comments on commit c15b7e2

Please sign in to comment.