diff --git a/common/Domain.cpp b/common/Domain.cpp index 945a6c6a..14dc17dd 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -217,17 +217,23 @@ void Domain::read_swc(const std::string &Filename) { double z = k*voxel_length + start_z; double distance; - double s = ((x-xp)*alpha+(y-yp)*beta+(z-zp)*gamma) / (alpha*alpha + beta*beta + gamma*gamma); + double s = ((x-xp)*alpha+(y-yp)*beta+(z-zp)*gamma) / (alpha*alpha + beta*beta + gamma*gamma); + + double di = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi)); + double dp = rp - sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp)); + if (s > length){ - distance = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi)); + distance = di; } else if (s < 0.0){ - distance = rp - sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp)); + distance = dp; } else { // linear variation for radius double radius = rp + (ri - rp)*s/length; distance = radius - sqrt((x-xp-alpha*s)*(x-xp-alpha*s) + (y-yp-beta*s)*(y-yp-beta*s) + (z-zp-gamma*s)*(z-zp-gamma*s)); + if (distance < di) distance = di; + if (distance < dp) distance = dp; } if ( distance > 0.0 ){ /* label the voxel */