s(i, :));
    p1 = center - 30*cart;
    p2 = center + 30*cart;
    cyl = [p1 p2 10];

    <span class="comment">% generation of 3D image</span>
    img = discreteCylinder(lx, ly, lz, cyl);

    <span class="comment">% display image isosurface</span>
    clf; set(gca, <span class="string">'fontsize'</span>, 14);
    patch(isosurface(img, .5), <span class="string">'FaceColor'</span>, <span class="string">'g'</span>, <span class="string">'LineStyle'</span>, <span class="string">'none'</span>);

    <span clas