C++ ZSwcTree类代码示例

ZSwcTree* ZSwcGenerator::createSwc(
    const ZObject3dScan &blockObj, int z, const ZDvidInfo &dvidInfo)
{
#ifdef _FLYEM_
  ZObject3dScan slice = blockObj.getSlice(z);
  size_t stripeNumber = slice.getStripeNumber();

  ZSwcTree *tree = new ZSwcTree;
  for (size_t s = 0; s < stripeNumber; ++s) {
    const ZObject3dStripe &stripe = slice.getStripe(s);
    int nseg = stripe.getSegmentNumber();
    int y = stripe.getY();
    int z = stripe.getZ();
    for (int i = 0; i < nseg; ++i) {
      int x0 = stripe.getSegmentStart(i);
      int x1 = stripe.getSegmentEnd(i);
      for (int x = x0; x <= x1; ++x) {
        ZIntCuboid cuboid = dvidInfo.getBlockBox(x, y, z);
        tree->merge(createBoxSwc(cuboid));
      }
    }
  }

  return tree;
#else
  UNUSED_PARAMETER(&blockObj);
  UNUSED_PARAMETER(z);
  UNUSED_PARAMETER(&dvidInfo);
  return NULL;
#endif
}

ZSwcTree* ZSwcGenerator::createSwc(const ZObject3dScan &obj)
{
  if (obj.isEmpty()) {
    return NULL;
  }

  ZSwcTree *tree = new ZSwcTree();
  tree->forceVirtualRoot();
  Swc_Tree_Node *root = tree->root();

  size_t stripeNumber = obj.getStripeNumber();
  for (size_t i = 0; i < stripeNumber; ++i) {
    const ZObject3dStripe &stripe = obj.getStripe(i);
    int segNumber = stripe.getSegmentNumber();
    int y = stripe.getY();
    int z = stripe.getZ();
    for (int j = 0; j < segNumber; ++j) {
      Swc_Tree_Node *tn =
          SwcTreeNode::makePointer(stripe.getSegmentStart(j), y, z, 2.0);
      SwcTreeNode::setFirstChild(root, tn);
      Swc_Tree_Node *tn2 =
          SwcTreeNode::makePointer(stripe.getSegmentEnd(j), y, z, 2.0);
      SwcTreeNode::setFirstChild(tn, tn2);
    }
  }

  tree->resortId();

  return tree;
}

ZSwcTree* ZSwcGenerator::createSurfaceSwc(
    const ZObject3dScan &obj, int sparseLevel)
{
  size_t volume = obj.getBoundBox().getVolume();

  int intv = 0;
  if (volume > MAX_INT32) {
    intv = iround(Cube_Root((double) volume / MAX_INT32));
  }

  ZStack *stack = NULL;
  std::cout << "Creating object mask ..." << "ds: " << intv <<  std::endl;
  if (intv > 0) {
    ZObject3dScan obj2 = obj;
    obj2.downsampleMax(intv, intv, intv);
    stack = obj2.toStackObject();
  } else {
    stack = obj.toStackObject();
  }

  ZSwcTree *tree = NULL;
  if (stack != NULL) {
    tree = createSurfaceSwc(*stack, sparseLevel);
    tree->setColor(obj.getColor());
    tree->rescale(intv + 1, intv + 1, intv + 1);
    delete stack;
  }

  return tree;
}

ZSwcTree* ZSwcGenerator::createSwc(const ZFlyEmNeuronRange &range)
{
  if (range.isEmpty()) {
    return NULL;
  }

  double minZ = range.getMinZ();
  double maxZ = range.getMaxZ();

  double dz = (maxZ - minZ) / 50.0;

  if (dz == 0.0) { //Avoid dead loop
    dz = 1.0;
  }

  ZSwcTree *tree = createVirtualRootSwc();

  for (double z = minZ; z <= maxZ; z += dz) {
    double r = range.getRadius(z);
    ZSwcTree *subtree = createCircleSwc(0, 0, z, r);
    tree->merge(subtree, true);
  }

  tree->resortId();

  return tree;
}

void ZSwcNodeObjsModel::setupModelData(ZObjsItem *parent)
{
  QList<QVariant> data;

  m_typeToRow.clear();
  m_swcTreeNodeToRow.clear();
  m_swcTreeNodeToType.clear();

  m_typeToRow[SwcTreeNode::TERMINAL] = 0;
  m_typeToRow[SwcTreeNode::BRANCH_POINT] = 1;

  data.clear();
  data << "Termini" << "id" << "type" << "radius" << "x" << "y" << "z"
       << "label";
  ZObjsItem *terminalItem = new ZObjsItem(data, NULL, parent);
  terminalItem->setCheckState(Qt::Checked);
  parent->appendChild(terminalItem);

  data.clear();
  data << "Branch Points" << "id" << "type" << "radius" << "x" << "y" << "z"
       << "label";
  ZObjsItem *branchPointItem = new ZObjsItem(data, NULL, parent);
  branchPointItem->setCheckState(Qt::Checked);
  parent->appendChild(branchPointItem);

  int terminalRow = 0;
  int branchPointRow = 0;

  QList<ZSwcTree*> swcList = m_doc->getSwcList();
  for (int i=0; i<swcList.size(); i++) {
    data.clear();
    ZSwcTree *swcTree = swcList.at(i);

    //ZObjsItem *nodeParent = new ZObjsItem(data, swcTree, parent);
    //nodeParent->setCheckState(swcTree->isVisible() ? Qt::Checked : Qt::Unchecked);
    //nodeParent->setToolTip(QString("source: %1").arg(QString::fromStdString(swcTree->source())));
    //parent->appendChild(nodeParent);

    swcTree->updateIterator(SWC_TREE_ITERATOR_DEPTH_FIRST);   //depth first
    for (Swc_Tree_Node *tn = swcTree->begin(); tn != swcTree->end(); tn = swcTree->next()) {
      if (!SwcTreeNode::isVirtual(tn)) {
        data.clear();
        data << "" << tn->node.id << tn->node.type << tn->node.d << tn->node.x
             << tn->node.y << tn->node.z  << tn->node.label << "";
        if (SwcTreeNode::isBranchPoint(tn)) {
          m_swcTreeNodeToType[tn] = SwcTreeNode::BRANCH_POINT;
          m_swcTreeNodeToRow[tn] = branchPointRow++;
          ZObjsItem *node = new ZObjsItem(data, tn, branchPointItem);
          branchPointItem->appendChild(node);
        } else if (SwcTreeNode::isRoot(tn) || SwcTreeNode::isLeaf(tn)) {
          m_swcTreeNodeToType[tn] = SwcTreeNode::TERMINAL;
          m_swcTreeNodeToRow[tn] = terminalRow++;
          ZObjsItem *node = new ZObjsItem(data, tn, terminalItem);
          terminalItem->appendChild(node);
        }
      }
    }
  }
}

Swc_Tree_Node* ZStackOperator::getHitObject<Swc_Tree_Node>() const
{
  Swc_Tree_Node *tn = NULL;
  ZSwcTree *tree = getHitObject<ZSwcTree>();
  if (tree != NULL) {
    tn = tree->getHitNode();
  }

  return tn;
}

ZSwcTree* ZSwcGenerator::createSwcByRegionSampling(
    const ZVoxelArray &voxelArray, double radiusAdjustment)
{
#ifdef _DEBUG_2
  voxelArray.print();
#endif

  ZDoubleVector voxelSizeArray(voxelArray.size());

  const std::vector<ZVoxel> &voxelData = voxelArray.getInternalData();
  //Retrieve voxel size
  for (size_t i = 0; i < voxelSizeArray.size(); ++i) {
    voxelSizeArray[i] = -voxelData[i].value();
  }

  std::vector<int> indexArray;
  voxelSizeArray.sort(indexArray);

  std::vector<bool> sampled(voxelArray.size(), true);

  for (size_t i = 1; i < voxelArray.size(); ++i) {
    size_t currentVoxelIndex = indexArray[i];
    const ZVoxel &currentVoxel = voxelData[currentVoxelIndex];
    for (size_t j = 0; j < i; ++j) {
      size_t prevVoxelIndex = indexArray[j];
      if (sampled[prevVoxelIndex]) {
        const ZVoxel &prevVoxel = voxelData[prevVoxelIndex];
        double dist = currentVoxel.distanceTo(prevVoxel);
        if (dist < prevVoxel.value()) {
          sampled[currentVoxelIndex] = false;
          break;
        }
      }
    }
  }

  Swc_Tree_Node *prevTn = NULL;

  for (size_t i = 0; i < voxelArray.size(); ++i) {
    if (sampled[i]) {
      Swc_Tree_Node *tn = SwcTreeNode::makePointer();
      SwcTreeNode::setPos(
            tn, voxelData[i].x(), voxelData[i].y(), voxelData[i].z());
      SwcTreeNode::setRadius(
            tn, voxelData[i].value() + radiusAdjustment);
      Swc_Tree_Node_Set_Parent(tn, prevTn);
      prevTn = tn;
    }
  }

  ZSwcTree *tree = new ZSwcTree;
  tree->setDataFromNodeRoot(prevTn);

  return tree;
}

Geo3d_Scalar_Field* ZNeuronTracer::extractSeedSkel(const Stack *mask)
{
  Stack *skel = Stack_Bwthin(mask, NULL);

  /* alloc <dist> */
  Stack *dist = Stack_Bwdist_L_U16(mask, NULL, 0);


  ZStackProcessor::RemoveBranchPoint(skel, 26);

  Stack *skel_proc = C_Stack::clone(skel);
 Geo3d_Scalar_Field *field1 = extractLineSeed(skel_proc, dist);
 C_Stack::kill(skel_proc);

 for (int i = 0; i <field1->size; ++i) {
   int x = field1->points[i][0];
   int y = field1->points[i][1];
   int z = field1->points[i][2];
   Set_Stack_Pixel(skel, x, y, z, 0, 0);
 }

#ifdef _DEBUG_2
  C_Stack::write(GET_TEST_DATA_DIR + "/test.tif", skel);
#endif


 Geo3d_Scalar_Field *field2 = extractLineSeed(skel, dist, 0);

 Geo3d_Scalar_Field *field = Geo3d_Scalar_Field_Merge(field1, field2, NULL);

 Kill_Geo3d_Scalar_Field(field1);
 Kill_Geo3d_Scalar_Field(field2);

#ifdef _DEBUG_2
  ZSwcTree tree;
  tree.forceVirtualRoot();
  for (int i = 0; i <field->size; ++i) {
    int x = field->points[i][0];
    int y = field->points[i][1];
    int z = field->points[i][2];
    double radius = field->values[i];
    SwcTreeNode::setFirstChild(
          tree.root(), SwcTreeNode::makePointer(x, y, z, radius));
  }
  tree.save(GET_TEST_DATA_DIR + "/test.swc");
#endif

  /* free <dist> */
  C_Stack::kill(dist);

  C_Stack::kill(skel);

  return field;
}

ZSwcTree *ZNeuronConstructor::reconstruct(
    std::vector<Locseg_Chain*> &chainArray)
{
  ZSwcTree *tree = NULL;

  if (!chainArray.empty()) {
    int chain_number = chainArray.size();
    /* <neuronComponent> allocated */
    Neuron_Component *neuronComponent =
        Make_Neuron_Component_Array(chain_number);

    for (int i = 0; i < chain_number; i++) {
      Set_Neuron_Component(neuronComponent + i,
                           NEUROCOMP_TYPE_LOCSEG_CHAIN,
                           chainArray[i]);
    }

    /* reconstruct neuron */
    /* alloc <ns> */
    double zscale = 1.0;
    Neuron_Structure *ns = Locseg_Chain_Comp_Neurostruct(
          neuronComponent, chain_number, m_signal, zscale, m_connWorkspace);

    Process_Neuron_Structure(ns);

    if (m_connWorkspace->crossover_test == TRUE) {
      Neuron_Structure_Crossover_Test(ns, zscale);
    }

    /* alloc <ns2> */
    Neuron_Structure* ns2=
        Neuron_Structure_Locseg_Chain_To_Circle_S(ns, 1.0, 1.0);

    Neuron_Structure_To_Tree(ns2);

    tree = new ZSwcTree;
    tree->setData(Neuron_Structure_To_Swc_Tree_Circle_Z(ns2, 1.0, NULL));
    tree->resortId();

    /* free <ns2> */
    Kill_Neuron_Structure(ns2);
    /* free <ns> */
    ns->comp = NULL;
    Kill_Neuron_Structure(ns);

    /* free <neuronComponent> */
    Clean_Neuron_Component_Array(neuronComponent, chain_number);
    free(neuronComponent);
  }

  return tree;
}

ZSwcTree* ZNeuronTracer::trace(ZStack *stack, bool doResampleAfterTracing)
{
  ZSwcTree *tree = NULL;

  if (stack != NULL) {
    tree = trace(stack->c_stack(), doResampleAfterTracing);
    if (tree != NULL) {
      tree->translate(stack->getOffset());
    }
  }

  return tree;
}

std::vector<double> ZSwcTrunkSizeFeatureAnalyzer::computeFeature(Swc_Tree_Node *tn)
{
  ZSwcTree tree;
  tree.setDataFromNodeRoot(tn);

  tree.labelTrunkLevel(m_trunkAnalyzer);

  vector<double> featureArray(1);
  featureArray[0] =
      SwcTreeNode::downstreamSize(tn, SwcTreeNode::labelDifference);

  return featureArray;
}

std::vector<Swc_Tree_Node*> ZSwcRangeAnalyzer::getOutsideNode(
    const ZSwcTree &host, const ZSwcTree &targetTree)
{
  std::vector<Swc_Tree_Node*> nodeArray;

  std::vector<ZEllipsoid> range = computeLayerRange(host);
  ZSwcNodeCompositeSelector selector;
  for (std::vector<ZEllipsoid>::iterator iter = range.begin();
       iter != range.end(); ++iter) {
    iter->setCenterX(0);
    iter->setCenterY(0);
    ZSwcNodeEllipsoidRangeSelector *childSelector =
        new ZSwcNodeEllipsoidRangeSelector;
    childSelector->setExpandFactor(3.0);
    //iter->print();
    /*
    iter->setSize(iter->getXRadius() + sqrt(iter->getXRadius()),
                  iter->getYRadius() + sqrt(iter->getYRadius()),
                  iter->getZRadius() + sqrt(iter->getZRadius()) * 5);
                  */

    childSelector->setRange(*iter);
    selector.addSelector(childSelector);
  }

  std::vector<ZEllipsoid> range2 = computeLayerRange(targetTree);
  std::vector<ZPoint> pointArray = computerLayerCentroid(targetTree);

  Swc_Tree_Node bufferNode;
  targetTree.updateIterator(SWC_TREE_ITERATOR_DEPTH_FIRST);
  for (Swc_Tree_Node *tn = targetTree.begin(); tn != NULL;
       tn = targetTree.next()) {
    if (SwcTreeNode::isRegular(tn)) {
      bufferNode = *tn;
      for (size_t i = 0; i < range2.size(); ++i) {
        if (range2[i].containsPoint(
              SwcTreeNode::x(tn), SwcTreeNode::y(tn), SwcTreeNode::z(tn))) {
          SwcTreeNode::translate(&bufferNode, -pointArray[i].x(),
                                 -pointArray[i].y(), 0.0);
          break;
        }
      }
      if (!selector.isSelected(&bufferNode)) {
        nodeArray.push_back(tn);
      }
    }
  }

  return nodeArray;
}

ZSwcTree* ZObject3dPlayer::getSwcDecoration() const
{
  const ZObject3d *obj = getCompleteData();

  ZSwcTree *tree = NULL;
  if (obj != NULL) {
    if (!obj->isEmpty()) {
      tree = ZSwcGenerator::createSwc(*obj, 1.0, 3);
      tree->setColor(obj->getColor());
    }
  }

  return tree;
}

pair<Swc_Tree_Node*, Swc_Tree_Node*> ZSwcConnector::identifyConnection(
      const ZSwcTree &hook, const ZSwcTree &loop)
{
  pair<Swc_Tree_Node*, Swc_Tree_Node*> conn;

  loop.updateIterator(SWC_TREE_ITERATOR_DEPTH_FIRST);
  m_dist = Infinity;

  ZSwcTree::TerminalIterator hookIter(&hook);

  ZSwcTree::DepthFirstIterator loopIter(&loop);
  while (loopIter.hasNext()) {
    Swc_Tree_Node *tn = loopIter.next();
    if (SwcTreeNode::isRegular(tn)) {
      hookIter.begin();
      while (hookIter.hasNext()) {
        Swc_Tree_Node *head = hookIter.next();
        double d = SwcTreeNode::scaledDistance(
              head, tn, m_resolution.voxelSizeX(), m_resolution.voxelSizeY(),
              m_resolution.voxelSizeZ());
        if (d < m_dist) {
          m_dist = d;
          conn.first = head;
          conn.second = const_cast<Swc_Tree_Node*>(tn);
        }
      }
    }
  }

  return conn;
}

std::vector<ZEllipsoid> ZSwcRangeAnalyzer::computeLayerRange(const ZSwcTree &tree)
{
  std::vector<ZEllipsoid> range;

  ZCuboid treeBound = tree.boundBox();

  double zStart = treeBound.firstCorner().z();
  double zEnd = treeBound.lastCorner().z();

  ZSwcNodeZRangeSelector selector;
  for (double z = zStart; z <= zEnd; z += m_zStep) {
    selector.setZRange(z - m_zMargin, z + m_zMargin);
    std::vector<Swc_Tree_Node*> nodeSet = selector.selectFrom(tree);
    if (!nodeSet.empty()) {
      ZCuboid bound = SwcTreeNode::boundBox(nodeSet.begin(), nodeSet.end());
      ZEllipsoid ellipsoid;
      ellipsoid.setCenter(bound.center());
      double rxy = sqrt(bound.width() * bound.width() +
                        bound.height() * bound.height()) / 2.0;
      ellipsoid.setSize(rxy, rxy, bound.depth() / 2.0);
      range.push_back(ellipsoid);
    }
  }

  return range;
}