/*
* Procedures for generating myelinated axon from original, unmyelinated morphology
* Preserves original geometry, replaces axon with new sections named Myelin, Node, and Unmyelin
* AUTHOR: Aman Aberra, Duke University
* CONTACT: aman.aberra@duke.edu
*/
// input scale factor and section list, scales pt3d diam info
proc scale_diam2() { local f,i,ii localobj diams, diams2, diams_sec,scale_seclist
f = $1 // scale factor
scale_seclist = $o2
diams = new Vector()
forsec scale_seclist {
diams_sec = getpt3d(5)
diams.append(diams_sec) // append this sections diam vector
}
diams2 = diams.mul(f)
i = 0
forsec scale_seclist {
for ii = 0, n3d() - 1{
pt3dchange(ii,diams2.x[i])
i += 1
}
}
}
// input myelin section list, scales based on diameter dependent g_ratio taken from (Micheva 2016)
proc scale_diam3() { local f, i, ii, nn, p1,p2,p3,p4 localobj diams, diams2, diams_sec, scale_seclist, g_ratios,g_ratios_sec, ones
scale_seclist = $o1
diams = new Vector()
g_ratios_sec = new Vector()
p1 = 0.6425 // polynomial curve fit of d vs. g_ratio from Micheva 2016 data
p2 = -1.778
p3 = 1.749
p4 = 0.1518
g_ratios = new Vector()
forsec scale_seclist {
diams_sec = getpt3d(5)
diams.append(diams_sec) // append this sections diam vector
g_ratios_sec = diams_sec.c.pow(3).mul(p1)
g_ratios_sec = g_ratios_sec.c.add(diams_sec.c.pow(2).mul(p2))
g_ratios_sec = g_ratios_sec.c.add(diams_sec.c.mul(p3))
g_ratios_sec = g_ratios_sec.c.add(p4)
for nn = 0, g_ratios_sec.size()-1 {
if (g_ratios_sec.x[nn] > 0.8) g_ratios_sec.x[nn] = 0.8 // set max g-ratio
if (g_ratios_sec.x[nn] < 0.40) g_ratios_sec.x[nn] = 0.40 // set min g-ratio
}
g_ratios.append(g_ratios_sec)
}
ones = new Vector(g_ratios.size())
ones = ones.c.fill(1)
diams2 = diams.c.mul(ones.c.div(g_ratios)) // scale each compartment diameter by its specific g_ratio
i = 0
forsec scale_seclist {
for ii = 0, n3d() - 1{
pt3dchange(ii,diams2.x[i])
i += 1
}
}
printf("Scaled diameter of myelin sections using variable g-ratio\n")
}
// Skip axon[0] convert every section to myelin + node
INL_ratio = 100
INL_ratio_term = 70 // INL ratio for terminations
nodeL = 1
min_myelinL = 20 // 20 minimum myelinated section length (Waxman 1970)
min_myelinD = 0.2 // 0.2 µm (Hildebrand 1993)
min_PMAS = 50 // minimum premyelin axon segment (PMAS), i.e. axon hillock + initial segment length
myelinL_error = 0.1 // allowed error in myelin length
nodeL_error = 0.1 // allowed error in node length
max_myelin_order = 0 // max_order - max_myelin_order is max branch order to myelinate. 0 myelinates all branches (same convention as prune_ax)
create Myelin[2], Node[2],Unmyelin[2]
objref iseg_secList, Node_secList,Myelin_secList,Unmyelin_secList, axonal
objref myelinCnts
load_file("myelinBiophysics.hoc") // load after section list objects instantiated
proc myelinate_axon() { localobj axon_secList
axon_secList = $o1
add_myelin(axon_secList)
//scale_diam2(1/g_ratio,Myelin_secList) // use constant g_ratio
scale_diam3(Myelin_secList) // use diam dependent g ratio using 3° polynomial fit to Micheva 2016 data
//scale_diam4(Node_secList) // use diam dependent g ratio to scale down node diameters
geom_nseg(40,axonal) // assign nseg to each section (2 per 40 um)
myelin_biophys()
forsec axonal cell.all.append() // add to cell sectionlist
forsec iseg_secList cell.axonal.append() // add initial segment to cell.axonal()
forsec axonal cell.axonal.append() // add rest of new myelinated axon to cell.axonal()
//forsec iseg_secList axonal.append() // add here to avoid re-adding to cell.all
}
// input axon section list, e.g. cell.axonal
proc add_myelin() { local cnt, myelinL, nn, include_PMAS_myelin, myelinL_sec, max_order localobj secx, secy, secz, length, diamvec, \
axon_secList, parent_SecList, children_SecList, child_SecRef
//pt3dconst(1)
axon_secList = $o1
iseg_secList = new SectionList()
Myelin_secList = new SectionList()
Node_secList = new SectionList()
Unmyelin_secList = new SectionList()
axonal = new SectionList()
forsec "axon[0]" iseg_secList.append()
axon_secList.remove(iseg_secList) // remove axon[0] from list
forsec iseg_secList {
if (L >= (min_PMAS + min_myelinL) ) { // add myelin/node before 1st bifurcation
numMyelin=1 // account for first myelin and nodal sections
numNode=1
include_PMAS_myelin = 1
} else { // don't add myelin/node before 1st bifurcation
numMyelin=0 // start count at 0 for rest of arbor
numNode=0
include_PMAS_myelin = 0
}
//axonal.append() // don't add iseg to axonal seclist until adding axonal to all
}
numAxonal = 0
numUnmyelin = 0
objref myelinCnts
myelinCnts = new Vector() // number of internodal sections to replace each existing axonal section
// calculate number of new myelin, nodal, and unmyelinated sections to create
max_order = get_max_order(axon_secList) // get maximum branch order of original axon
forsec axon_secList {
numAxonal += 1
if (type_xtra(1)==2||type_xtra(1)==5){
myelinL=diam(0)*INL_ratio_term
} else {
myelinL = diam(0)*INL_ratio
}
// 3 conditions for myelination: 1) above min diam 2) above min length 3) below max branch order
// must be larger than minimum diameter for myelinated fibers
if (diam(0) >= min_myelinD) {
// must be below maximum branch order, e.g. if max_order = 10, max_myelin_order = 1 => myelinate up to order = 9
// unless section is in main axon, which is always myelinated (if diam is below min_myelinD)
if (check_in_secList(main_ax_list) || order_xtra < max_order + 1 - max_myelin_order) {
numMyelin_sec = int(L/(myelinL + nodeL))
if (numMyelin_sec == 0){ // int rounds down to 0 if < 1, no myelin if L < min_myelinL + nodeL
numMyelin_sec = int(L/(min_myelinL + nodeL)) // if section is shorter than internodal length calculated by INL_ratio*diam, use minimum myelin length
//numMyelin_sec = 1 // if section is shorter than INL from INL_ratio*diam, make full section myelinated
}
} else { // section is above max order and not part of main axon
numMyelin_sec = 0
}
} else {
numMyelin_sec = 0
}
numMyelin += numMyelin_sec
numNode += numMyelin_sec
if (numMyelin_sec==0) numUnmyelin+=1 // if section still too thin/short for minimum myelin L, leave unmyelinated
myelinCnts.append(numMyelin_sec) // keeps number of myelin sections per existing section
}
// create new axonal sections and add to SectionLists
create Myelin[numMyelin], Node[numNode]
if (numUnmyelin >= 1) {
create Unmyelin[numUnmyelin]
//printf("Number of unmyelinated sections = %g\n",numUnmyelin)
}
forsec "Myelin" {
Myelin_secList.append()
axonal.append()
}
forsec "Node" {
Node_secList.append()
axonal.append()
}
forsec "Unmyelin" {
Unmyelin_secList.append()
axonal.append()
}
printf("Myelinating axon: Replacing %g Axonal sections w/ %g Myelin, %g Node, %g Unmyelin sections\n",numAxonal,numMyelin,numNode,numUnmyelin)
// connect Myelin[0] to iseg if exists before 1st bifurcation
if (include_PMAS_myelin==1) {
print "Adding myelin before the 1st bifurcation"
forsec iseg_secList {
children_SecList = getchildren() // get existing children sections
forsec children_SecList {disconnect()} // disconnect children from cell.axon[0]
connect Myelin[0](0), 1 // connect 1st myelin section to iseg
//printf("Iseg: connected Myelin[0] to %s\n",secname())
// get coordinates of original iseg
secx = getpt3d(1)
secy = getpt3d(2)
secz = getpt3d(3)
length = getpt3d(4)
diamvec = getpt3d(5)
last_pt3d_i = length.indwhere(">",min_PMAS)-1 // until 1st point that's farther than the PMAS length
// remove points after last_pt3d_i, assign them to 1st myelin/node
while (arc3d(n3d()-1) > min_PMAS) {
pt3dremove(n3d()-1)
}
// get myelinL
myelinL = length.x[length.size()-1] - (min_PMAS + nodeL) // get new myelin length to fit all myelin/node sections in old section
}
Myelin[0] {
first_pt3d_i = last_pt3d_i
last_pt3d_i = length.indwhere(">",min_PMAS+myelinL) // until 1st point that's farther than myelinL
if (last_pt3d_i > first_pt3d_i) last_pt3d_i -= 1 // point before point that's farther than myelinL
//printf("Iseg: Myelin[%g], 1st pt = %g. Last point = %g\n",mye_cnt,first_pt3d_i,last_pt3d_i)
last_pt3d_i = add_new_points(first_pt3d_i,last_pt3d_i,secx,secy,secz,length,diamvec,myelinL,myelinL_error,1)
}
connect Node[0](0), Myelin[0](1) // connect 1st node to 1st myelin
Node[0] {
first_pt3d_i = last_pt3d_i
//length = get_arc3d(secx,secy,secz)
//last_pt3d_i = length.indwhere(">",myelinL + nodeL)
last_pt3d_i = length.size()-1 // set to last coordinate
if (last_pt3d_i > first_pt3d_i) {
last_pt3d_i -= 1
} else { // if less than or equal to 1st pt, use same point
last_pt3d_i = first_pt3d_i
}
//printf("Iseg: Node[0], 1st pt = %g. Last point = %g\n",first_pt3d_i,last_pt3d_i)
/*
if (first_pt3d_i == last_pt3d_i){
//printf(" Adding 1st interpolated point\n")
add_interp_pt(last_pt3d_i,secx,secy,secz,diamvec,nodeL,-1) // add additional point along same direction as pt3d with distance nodeL
}
*/
//assign_pts(first_pt3d_i,last_pt3d_i,secx,secy,secz,diamvec) // add 2nd point using original coordinates
last_pt3d_i = add_new_points(first_pt3d_i,last_pt3d_i,secx,secy,secz,length,diamvec,nodeL,nodeL_error,1)
//printf("Iseg: Node[0].L = %.2f\n",L)
}
forsec children_SecList { // reattach children of iseg
//printf("Disconnecting %s from parent\n",secname())
child_SecRef = new SectionRef()
disconnect() // disconnect from original parent (current section)
connect child_SecRef.sec(0), Node[0](1)
//printf("Reconnecting to Node[%g]\n",0)
}
mye_cnt = 1 // myelin counter starts at 1, because 0 already added
//print "Added myelin before 1st bifurcation"
} else {
// otherwise, first myelinated sections will be attached to Node[0], start at 0
mye_cnt = 0
print "No myelin before 1st bifurcation"
}
// Assign parameters for Myelin[0] if necessary
// Deal with rest of axon
sec_cnt = 0 // section counter, myelin section counter above (mye_cnt)
unmye_cnt = 0 // unmyelinated section counter
forsec axon_secList { // replace each section with (myelin+node) units
// get coordinates of current section
secx = getpt3d(1)
secy = getpt3d(2)
secz = getpt3d(3)
length = getpt3d(4) // outputs lengths of 3d points (not normalized)
diamvec = getpt3d(5) // outputs diameter of 3d points
// get children and parent sections in SectionLists
children_SecList = getchildren()
parent_SecList = getparent()
numMyelin_sec = myelinCnts.x[sec_cnt] // get number of myelin sections to create
//node_diam = diam(0) // same node diameter for all nodes in this section
//myelin_diam = node_diam/g_ratio // same myelin diameter for all myelin in this section
// Now that we've gotten all info we need, delete current section
delete_section()
myelinL_sec = 0 // start at 0 for every axonal section being myelinated
// Start connecting new sections replacing original axonal section
if (numMyelin_sec >=1 ) {
//printf("First pt sec (%f,%f,%f) \n",secx.x[0],secy.x[0],secz.x[0])
forsec parent_SecList {
connect Myelin[mye_cnt](0), 1 // connect first myelin to end of original parent
// start coordinates at end of parent (in case original parent has already been replaced)
secx.x[0] = x3d(n3d()-1)
secy.x[0] = y3d(n3d()-1)
secz.x[0] = z3d(n3d()-1)
diamvec.x[0] = diam3d(n3d()-1)
length = get_arc3d(secx,secy,secz) // get new length vector with updated x,y,z vectors
//printf("New first pt sec (%f,%f,%f) \n",secx.x[0],secy.x[0],secz.x[0])
myelinL = (length.x[length.size()-1] - numMyelin_sec*nodeL)/numMyelin_sec // get new myelin length to fit all myelin/node sections in old section
// connect first myelin and nodal section to original parent
//printf("Creating %g myelin (L=%.3f um) and nodes. n3d() = %g\n",numMyelin_sec,myelinL,length.size())
//printf("Connecting 1st Myelin[%g] to its parent: %s\n",mye_cnt,secname())
}
Myelin[mye_cnt] connect Node[mye_cnt](0), 1 // connect first node to end of first myelin
// assign coordinates and diameters for 1st myelinated section
Myelin[mye_cnt] {
first_pt3d_i = 0 // start from 1st pt3d
last_pt3d_i = length.indwhere(">",myelinL) // until 1st point that's farther than myelinL
if (last_pt3d_i > first_pt3d_i) {
last_pt3d_i -= 1 // point before point that's farther than myelinL
}
//printf("Myelin[%g], 1st pt = %g. Last point = %g\n",mye_cnt,first_pt3d_i,last_pt3d_i)
//***** Add new point (before or after last_pt3d_i) or use existing ones to trace out myelin path
last_pt3d_i = add_new_points(first_pt3d_i,last_pt3d_i,secx,secy,secz,length,diamvec,myelinL,myelinL_error,1)
//*****
myelinL_sec = L // add length of 1st myelin section to section's running total
//printf("L = %.3f. myelinL_sec = %.2f\n",L,myelinL_sec)
}
// loop through remaining nodal/myelin sections, connect and give pt3d coordinates
if (numMyelin_sec >= 2) {
for nn = 0, numMyelin_sec-1 {
//printf("Adding myelin+node section #%g\n",nn)
Node[mye_cnt+nn] { // start with 1st node - already connected to 1st myelin
first_pt3d_i = last_pt3d_i // start at last point of previous section
//first_pt3d_i = length.indwhere(">",myelinL_sec)
last_pt3d_i = length.indwhere(">",myelinL_sec + nodeL)
if (last_pt3d_i > first_pt3d_i) {
last_pt3d_i -= 1
} else { // if less than or equal to 1st pt, use same point
last_pt3d_i = first_pt3d_i
}
//last_pt3d_i = first_pt3d_i + 1 // use next point for interpolation
//last_pt3d_i = length.indwhere(">",(myelinL + nodeL)*(nn+1))
//printf("Node[%g],1st pt = %g. Last point = %g\n",mye_cnt+nn,first_pt3d_i,last_pt3d_i)
//assign_pts(first_pt3d_i,first_pt3d_i,secx,secy,secz,diamvec) // add first point from original coordinates
// Add new point 1 µm in the direction of original axon, modify coordinate vectors
last_pt3d_i = add_new_points(first_pt3d_i,last_pt3d_i,secx,secy,secz,length,diamvec,nodeL,nodeL_error,-1)
//assign_pts(first_pt3d_i,last_pt3d_i,secx,secy,secz,diamvec)
// If adding another myelinated section, connect it to end of node
if (nn < numMyelin_sec-1) connect Myelin[mye_cnt+nn+1](0), 1 // connect 1 end of current node (parent) to 0 end of next myelin (child)
myelinL_sec += L
//printf("Node L=%.2f. myelinL_sec = %f\n",L,myelinL_sec)
}
if (nn < numMyelin_sec-1) { // end section on node
Myelin[mye_cnt+nn+1] {
first_pt3d_i = last_pt3d_i
//first_pt3d_i = length.indwhere(">",myelinL_sec + nodeL*(nn+1)) // start at length after current total myelin length
//first_pt3d_i = length.indwhere(">",myelinL_sec) // start at length after current total myelin length (incl node)
//last_pt3d_i = length.indwhere(">",myelinL+myelinL_sec + nodeL*(nn+1))
last_pt3d_i = length.indwhere(">",myelinL_sec+myelinL) // running total incl node
//printf("last_i = %g. (%f,%f,%f). length(last_i) = %f. myelinL_sec = %.3f + myelinL = %.3f\n",last_pt3d_i,secx.x[last_pt3d_i],secy.x[last_pt3d_i],secz.x[last_pt3d_i],length.x[last_pt3d_i],myelinL_sec,myelinL)
//last_pt3d_i = length.indwhere(">",myelinL*(nn+2) + nodeL*(nn+1))
if (last_pt3d_i > first_pt3d_i) {
last_pt3d_i -= 1 // point before point that's farther than myelinL
} else if (last_pt3d_i < 0) {
myelinL = length.x[length.size()-1] - length.x[first_pt3d_i] - nodeL // readjust myelinL for remaining part of section
last_pt3d_i = length.size() - 2 // use second to last point
}
//printf("Myelin[%g],1st pt = %g. Last point = %g\n",mye_cnt+nn+1,first_pt3d_i,last_pt3d_i)
last_pt3d_i = add_new_points(first_pt3d_i,last_pt3d_i,secx,secy,secz,length,diamvec,myelinL,myelinL_error,1)
/*
if (first_pt3d_i == last_pt3d_i){
printf(" Adding interpolated point\n")
// add additional point along same direction as pt3d with distance myelinL
add_interp_pt(last_pt3d_i,secx,secy,secz,diamvec,myelinL,1) // use negative to interpolate forward
}
assign_pts(first_pt3d_i,last_pt3d_i,secx,secy,secz,diamvec)
*/
myelinL_sec += L // add to running total
//printf("L = %.2f. myelinL_sec = %.2f\n",L, myelinL_sec)
connect Node[mye_cnt+nn+1](0), 1 // connect 0 end of next node (child) to 1 end of current myelin
}
}
}
} else {
Node[mye_cnt] {
//first_pt3d_i = last_pt3d_i+1
//last_pt3d_i = length.indwhere(">",(myelinL + nodeL)*(nn+1))
first_pt3d_i = last_pt3d_i
last_pt3d_i = length.size()-1 // use last coordinate
/*
if (last_pt3d_i > first_pt3d_i) {
last_pt3d_i -= 1
} else { // if less than or equal to 1st pt, use same point
last_pt3d_i = first_pt3d_i
}
if (first_pt3d_i == last_pt3d_i){
printf(" Adding 1st interpolated point\n")
add_interp_pt(last_pt3d_i,secx,secy,secz,diamvec,nodeL,-1) // add additional point along same direction as pt3d with distance nodeL
}
assign_pts(first_pt3d_i,last_pt3d_i,secx,secy,secz,diamvec)
*/
last_pt3d_i = add_new_points(first_pt3d_i,last_pt3d_i,secx,secy,secz,length,diamvec,nodeL,nodeL_error,1)
}
}
// connect existing children to last node
forsec children_SecList {
//printf("Disconnecting %s from parent\n",secname())
child_SecRef = new SectionRef()
disconnect() // disconnect from original parent (current section)
connect child_SecRef.sec(0), Node[mye_cnt+numMyelin_sec-1](1)
//printf("Reconnecting to Node[%g]\n",mye_cnt+numMyelin_sec-1)
}
mye_cnt += numMyelin_sec
} else { // leave section unmyelinated
//printf("Leaving section unmyelinated\n")
// connect to parent
forsec parent_SecList connect Unmyelin[unmye_cnt](0), 1
// assign coordinates
Unmyelin[unmye_cnt] {
assign_pts(0,secx.size()-1,secx,secy,secz,diamvec)
}
// connect to children
forsec children_SecList {
//printf("Disconnecting %s from parent\n",secname())
child_SecRef = new SectionRef()
disconnect() // disconnect from original parent (current section)
connect child_SecRef.sec(0), Unmyelin[unmye_cnt](1) // connect end of unmyelianted section to original children
//printf("Reconnecting to Unmyelin[%g]\n",unmye_cnt)
}
unmye_cnt += 1
}
sec_cnt += 1 // increment section counter 1
}
}
obfunc getpt3d() { local dim, nn, ii localobj vec
// argument is 1, 2, 3, or 4, for x, y, z, or arc3d
dim = $1
nn = n3d()
vec = new Vector(nn)
if (dim == 1){
for ii = 0, nn-1 vec.x[ii] = x3d(ii)
} else if (dim == 2){
for ii = 0, nn-1 vec.x[ii] = y3d(ii)
} else if (dim == 3){
for ii = 0, nn-1 vec.x[ii] = z3d(ii)
} else if (dim == 4){
for ii = 0, nn-1 vec.x[ii] = arc3d(ii)
//vec.div(vec.x[nn-1]) // normalize length
} else if (dim == 5){
for ii = 0, nn-1 vec.x[ii] = diam3d(ii)
}
return vec
}
// Returns section list with single entry for currently accessed section's parent
obfunc getparent() { localobj current_sec, parent
current_sec = new SectionRef()
parent = new SectionList()
current_sec.parent() {
parent.append()
// print "Parent: ", secname()
}
return parent
}
// Returns section list with entries for currently accessed section's parent
obfunc getchildren() { local i localobj current_sec, children
children = new SectionList()
children.children() // append children
//forsec children print "Child: ", secname()
return children
}
// assign pt3d points from original section to current new section
// input first, last indices and coordinate vectors
// assign_pts(i1,i2,x,y,z)
proc assign_pts() { local i, i1, i2 localobj x, y, z, diamvec
i1 = $1
i2 = $2
x = $o3
y = $o4
z = $o5
diamvec = $o6
for i = i1,i2 {
pt3dadd(x.x[i], y.x[i],z.x[i], diamvec.x[i]) // add selected points from vectors to current section
}
//printf("Added points from %g to %g. (%.3f,%.3f,%.3f) to (%.3f,%.3f,%.3f)\n",i1, i2, x.x[i1],y.x[i1],z.x[i1],x.x[i2],y.x[i2],z.x[i2])
}
// adds additional pt3d point using direction of current point and either next or previous point
// at interpolated distance using input length
// dir is -1 or 1
// add_interp_pt(last_pt3d_i,secx,secy,secz,diamvec,len,dir)
obfunc add_interp_pt() { local i1, len, xu, yu, zu,xn,yn,zn,dir localobj x,y,z,diams, inds, interp_pt
i1 = $1 // current coordinate
x = $o2
y = $o3
z = $o4
diams = $o5
len = $6
dir = $7
inds = new Vector()
if (dir == 1){
inds.append(i1,i1+1) // extract i1 and next index to interpolate forward
x = x.ind(inds) // extract current coordinate and previous one
y = y.ind(inds)
z = z.ind(inds)
xu = x.x[1] - x.x[0] // x displacement (forward)
yu = y.x[1] - y.x[0]
zu = z.x[1] - z.x[0]
r = sqrt(xu^2 + yu^2 + zu^2) // total displacement
// add point at distance len and direction -<xu,yu,zu> from i2th coordinate (should be within section)
xn = x.x[0] + len*xu/r // new x point
yn = y.x[0] + len*yu/r
zn = z.x[0] + len*zu/r
} else if (dir == -1) {
inds.append(i1-1,i1) // extract previous index and i2 to interpolate from last point
x = x.ind(inds) // extract current coordinate and previous one
y = y.ind(inds)
z = z.ind(inds)
xu = x.x[1] - x.x[0] // x displacement (using previous point)
yu = y.x[1] - y.x[0]
zu = z.x[1] - z.x[0]
r = sqrt(xu^2 + yu^2 + zu^2) // total displacement
// add point at distance len and direction -<xu,yu,zu> from i2th coordinate (should be within section)
xn = x.x[1] + len*xu/r // new x point
yn = y.x[1] + len*yu/r
zn = z.x[1] + len*zu/r
}
pt3dadd(xn, yn, zn, diams.x[i1])
if (dir==1){
dist = sqrt( (xn - x.x[0] )^2 + (yn - y.x[0] )^2 + (zn - z.x[0] )^2 )
//printf("+1 Added point to %s (total %g): Start (%.5f,%.5f,%.5f). New (%.5f,%.5f,%.5f). Dist = %.3f\n",secname(),n3d(), x.x[0],y.x[0],z.x[0],xn,yn,zn,dist)
} else{
dist = sqrt( (xn - x.x[1])^2 + (yn - y.x[1])^2 + (zn - z.x[1])^2 )
//printf("-1 Added point to %s (total %g): Start (%.5f,%.5f,%.5f). New (%.5f,%.5f,%.5f). Dist = %.3f\n",secname(),n3d(), x.x[1],y.x[1],z.x[1],xn,yn,zn,dist)
}
interp_pt = new Vector()
interp_pt.append(xn,yn,zn) // output coordinates of new point
return interp_pt
}
// add_new_points(secx,secy,secz,length,diamvec,interp_pt)
func add_new_points() { local first_pt3d_i, last_pt3d_i, secL, secerr, secL_add, dir localobj secx, secy, secz, length, diamvec, interp_pt
first_pt3d_i = $1
last_pt3d_i = $2
secx = $o3
secy = $o4
secz = $o5
length = $o6
diamvec = $o7
secL = $8
secerr = $9
dir = $10
if (length.x[last_pt3d_i] - length.x[first_pt3d_i] >= secL + secerr ) {
// pt3d coordinates would make section too long
//printf("dist: %f - %f = %f. secL = %f. secerr = %f\n",length.x[last_pt3d_i],length.x[first_pt3d_i],length.x[last_pt3d_i] - length.x[first_pt3d_i],secL,secerr)
last_pt3d_i -= 1 // use 2nd to last coordinate
assign_pts(first_pt3d_i,last_pt3d_i,secx,secy,secz, diamvec) // add existing points (shortened)
secL_add = secL - (length.x[last_pt3d_i] - length.x[first_pt3d_i]) // distance of interpolated point from last existing point
//printf("Adding point to make section shorter at %g with secL_add = %.3f\n",last_pt3d_i,secL_add)
interp_pt = add_interp_pt(last_pt3d_i,secx,secy,secz,diamvec, secL_add, dir)
// replace old point to coordinate, length, and diameter vectors
/*
secx.x[last_pt3d_i] = interp_pt.x[0]
secy.x[last_pt3d_i] = interp_pt.x[1]
secz.x[last_pt3d_i] = interp_pt.x[2]
length.x[last_pt3d_i] = length.x[last_pt3d_i-1] + sqrt( (secx.x[last_pt3d_i] - interp_pt.x[0])^2 + (secy.x[last_pt3d_i]-interp_pt.x[1])^2 + (secz.x[last_pt3d_i]-interp_pt.x[2])^2 )
*/
// don't change diamvec entry, since this wasn't changed
secx.insrt(last_pt3d_i+1,interp_pt.x[0])
secy.insrt(last_pt3d_i+1,interp_pt.x[1])
secz.insrt(last_pt3d_i+1,interp_pt.x[2])
length.insrt(last_pt3d_i+1,length.x[last_pt3d_i]+sqrt((secx.x[last_pt3d_i] - interp_pt.x[0])^2 + (secy.x[last_pt3d_i]-interp_pt.x[1])^2 + (secz.x[last_pt3d_i]-interp_pt.x[2])^2) ) // insert distnace of new point from first point
//printf("Inserting length: %f + %f = %f\n",length.x[last_pt3d_i], sqrt( (secx.x[last_pt3d_i] - interp_pt.x[0])^2 + (secy.x[last_pt3d_i]-interp_pt.x[1])^2 + (secz.x[last_pt3d_i]-interp_pt.x[2])^2 ),length.x[last_pt3d_i] + sqrt( (secx.x[last_pt3d_i] - interp_pt.x[0])^2 + (secy.x[last_pt3d_i]-interp_pt.x[1])^2 + (secz.x[last_pt3d_i]-interp_pt.x[2])^2 ) )
diamvec.insrt(last_pt3d_i+1,diamvec.x[last_pt3d_i]) // insert same diameter of previous point
return last_pt3d_i+1 // position of new point in modified coordinate vectors
} else if (length.x[last_pt3d_i] - length.x[first_pt3d_i] <= (secL - secerr) ) {
// pt3d coordinates would make section too short
// add additional point using existing coordinates direction
//printf("dist2: %f - %f = %f. secL = %f. secerr = %f\n",length.x[last_pt3d_i],length.x[first_pt3d_i],length.x[last_pt3d_i] - length.x[first_pt3d_i],secL,secerr)
assign_pts(first_pt3d_i,last_pt3d_i,secx,secy,secz, diamvec) // add existing points (shortened)
secL_add = secL - (length.x[last_pt3d_i] - length.x[first_pt3d_i]) // distance of interpolated point from last existing point to make L=myelinL
//printf("Adding point to make section longer at %g with secL_add = %.3f\n",last_pt3d_i,secL_add)
interp_pt = add_interp_pt(last_pt3d_i,secx,secy,secz,diamvec,secL_add,dir)
// insert new point to coordinate, length, and diameter vectors
secx.insrt(last_pt3d_i+1,interp_pt.x[0])
secy.insrt(last_pt3d_i+1,interp_pt.x[1])
secz.insrt(last_pt3d_i+1,interp_pt.x[2])
//printf("Inserting length: %f + %f = %f\n",length.x[last_pt3d_i], sqrt( (secx.x[last_pt3d_i] - interp_pt.x[0])^2 + (secy.x[last_pt3d_i]-interp_pt.x[1])^2 + (secz.x[last_pt3d_i]-interp_pt.x[2])^2 ),length.x[last_pt3d_i] + sqrt( (secx.x[last_pt3d_i] - interp_pt.x[0])^2 + (secy.x[last_pt3d_i]-interp_pt.x[1])^2 + (secz.x[last_pt3d_i]-interp_pt.x[2])^2 ) )
length.insrt(last_pt3d_i+1, length.x[last_pt3d_i] + sqrt( (secx.x[last_pt3d_i] - interp_pt.x[0])^2 + (secy.x[last_pt3d_i]-interp_pt.x[1])^2 + (secz.x[last_pt3d_i]-interp_pt.x[2])^2 )) // insert distnace of new point from first point
diamvec.insrt(last_pt3d_i+1,diamvec.x[last_pt3d_i]) // insert same diameter of previous point
return last_pt3d_i+1 // position of new point in modified coordinate vectors
} else {
// just use existing points
//printf("dist: %f - %f = %f. secL = %f. secerr = %f\n",length.x[last_pt3d_i],length.x[first_pt3d_i],length.x[last_pt3d_i] - length.x[first_pt3d_i],secL,secerr)
assign_pts(first_pt3d_i,last_pt3d_i,secx,secy,secz,diamvec)
//printf("Used existing points, giving L= %.2f\n",L)
return last_pt3d_i // position of new point in modified coordinate vectors
}
}
// sets 2 comp per chunkSize in specified sectionList
// geom_nseg(chunkSize,secList)
proc geom_nseg() { local secIndex, chunkSize localobj seclist
chunkSize = 40
if( numarg() > 0 ) {
chunkSize = $1
}
seclist = $o2
secIndex=0
forsec seclist {
nseg = 1 + 2*int(L/chunkSize)
secIndex += 1
}
}
obfunc get_arc3d() {local i, dist_i localobj x, y, z, length
x = $o1
y = $o2
z = $o3
length = new Vector(x.size())
length.x[0] = 0 // first point starts at 0
for i = 1, x.size() -1 {
dist_i = sqrt( (x.x[i] - x.x[i-1])^2 + (y.x[i] - y.x[i-1])^2 + (z.x[i] - z.x[i-1])^2 )
length.x[i] = length.x[i-1] + dist_i
}
return length
}
// input sectionlist and check if currently accessed section is a member
// in_secList = check_in_secList(SectionList)
// returns 1 if in sectionlist
func check_in_secList() { localobj temp_secList
// copy into new sectionlist
temp_secList = new SectionList()
forsec $o1 temp_secList.append()
temp_secList.append() // append currently accessed section
return (temp_secList.unique > 0)
}
// max_order = get_max_order(SectionList)
// gets maximum branch order of input SectionList
// SectionList should have have xtra inserted and should have order_xtra defined (setpointers())
func get_max_order() { local max_order localobj seclist
seclist = $o1
max_order = 0
forsec seclist {
if (ismembrane("xtra")){
if (order_xtra > max_order) max_order = order_xtra // get max order
} else {
print "xtra not inserted in ", secname()
}
}
return max_order
}