Skeleton3D.jl

# from https://www.mathworks.com/matlabcentral/fileexchange/43400-skeleton3d



function FillEulerLUT()

    LUT=(zeros(Int8,255,1));
    
    LUT[1]=  1;
    LUT[3]= -1;
    LUT[5]= -1;
    LUT[7]=  1;
    LUT[9]= -3;
    LUT[11] = -1;
    LUT[13] = -1;
    LUT[15] =  1;
    LUT[17] = -1;
    LUT[19] =  1;
    LUT[21] =  1;
    LUT[23] = -1;
    LUT[25] =  3;
    LUT[27] =  1;
    LUT[29] =  1;
    LUT[31] = -1;
    LUT[33] = -3;
    LUT[35] = -1;
    LUT[37] =  3;
    LUT[39] =  1;
    LUT[41] =  1;
    LUT[43] = -1;
    LUT[45] =  3;
    LUT[47] =  1;
    LUT[49] = -1;
    LUT[51] =  1;
    
    LUT[53] =  1;
    LUT[55] = -1;
    LUT[57] =  3;
    LUT[59] =  1;
    LUT[61] =  1;
    LUT[63] = -1;
    LUT[65] = -3;
    LUT[67] =  3;
    LUT[69] = -1;
    LUT[71] =  1;
    LUT[73] =  1;
    LUT[75] =  3;
    LUT[77] = -1;
    LUT[79] =  1;
    LUT[81] = -1;
    LUT[83] =  1;
    LUT[85] =  1;
    LUT[87] = -1;
    LUT[89] =  3;
    LUT[91] =  1;
    LUT[93] =  1;
    LUT[95] = -1;
    LUT[97] =  1;
    LUT[99] =  3;
    LUT[101] =  3;
    LUT[103] =  1;
    
    LUT[105] =  5;
    LUT[107] =  3;
    LUT[109] =  3;
    LUT[111] =  1;
    LUT[113] = -1;
    LUT[115] =  1;
    LUT[117] =  1;
    LUT[119] = -1;
    LUT[121] =  3;
    LUT[123] =  1;
    LUT[125] =  1;
    LUT[127] = -1;
    LUT[129] = -7;
    LUT[131] = -1;
    LUT[133] = -1;
    LUT[135] =  1;
    LUT[137] = -3;
    LUT[139] = -1;
    LUT[141] = -1;
    LUT[143] =  1;
    LUT[145] = -1;
    LUT[147] =  1;
    LUT[149] =  1;
    LUT[151] = -1;
    LUT[153] =  3;
    LUT[155] =  1;
    
    LUT[157] =  1;
    LUT[159] = -1;
    LUT[161] = -3;
    LUT[163] = -1;
    LUT[165] =  3;
    LUT[167] =  1;
    LUT[169] =  1;
    LUT[171] = -1;
    LUT[173] =  3;
    LUT[175] =  1;
    LUT[177] = -1;
    LUT[179] =  1;
    LUT[181] =  1;
    LUT[183] = -1;
    LUT[185] =  3;
    LUT[187] =  1;
    LUT[189] =  1;
    LUT[191] = -1;
    LUT[193] = -3;
    LUT[195] =  3;
    LUT[197] = -1;
    LUT[199] =  1;
    LUT[201] =  1;
    LUT[203] =  3;
    LUT[205] = -1;
    LUT[207] =  1;
    
    LUT[209] = -1;
    LUT[211] =  1;
    LUT[213] =  1;
    LUT[215] = -1;
    LUT[217] =  3;
    LUT[219] =  1;
    LUT[221] =  1;
    LUT[223] = -1;
    LUT[225] =  1;
    LUT[227] =  3;
    LUT[229] =  3;
    LUT[231] =  1;
    LUT[233] =  5;
    LUT[235] =  3;
    LUT[237] =  3;
    LUT[239] =  1;
    LUT[241] = -1;
    LUT[243] =  1;
    LUT[245] =  1;
    LUT[247] = -1;
    LUT[249] =  3;
    LUT[251] =  1;
    LUT[253] =  1;
    LUT[255] = -1;
    return LUT
end

function p_EulerInv(img,LUT)
    #img=logical, LUT=int8, EulerInv=logical
    if numel(LUT) > 255
        display("skeleton3D:p_EulerInv:LUTwithTooManyElems","LUT with 255 elements expected");
    end
    # Calculate Euler characteristic for each octant and sum up
    eulerChar = int8(zeros(size(img,1),1));
    
    n=ones(Int8,size(img,1),1);#pre-make
    # Octant SWU
    bitorTable = uint8[128; 64; 32; 16; 8; 4; 2];
    n[:]=1;#restore to ones
    n[img[:,25]==1] = bitor[n[img[:,25]==1], bitorTable[1]];
    n[img[:,26]==1] = bitor[n[img[:,26]==1], bitorTable[2]];
    n[img[:,16]==1] = bitor[n[img[:,16]==1], bitorTable[3]];
    n[img[:,17]==1] = bitor[n[img[:,17]==1], bitorTable[4]];
    n[img[:,22]==1] = bitor[n[img[:,22]==1], bitorTable[5]];
    n[img[:,23]==1] = bitor[n[img[:,23]==1], bitorTable[6]];
    n[img[:,13]==1] = bitor[n[img[:,13]==1], bitorTable[7]];
    eulerChar = eulerChar + LUT[n];
    # Octant SEU
    n[:]=1;#restore to ones
    n[img[:,27]==1] = bitor[n[img[:,27]==1], bitorTable[1]];
    n[img[:,24]==1] = bitor[n[img[:,24]==1], bitorTable[2]];
    n[img[:,18]==1] = bitor[n[img[:,18]==1], bitorTable[3]];
    n[img[:,15]==1] = bitor[n[img[:,15]==1], bitorTable[4]];
    n[img[:,26]==1] = bitor[n[img[:,26]==1], bitorTable[5]];
    n[img[:,23]==1] = bitor[n[img[:,23]==1], bitorTable[6]];
    n[img[:,17]==1] = bitor[n[img[:,17]==1], bitorTable[7]];
    eulerChar = eulerChar + LUT[n];
    # Octant NWU
    n[:]=1;#restore to ones
    n[img[:,19]==1] = bitor[n[img[:,19]==1], bitorTable[1]];
    n[img[:,22]==1] = bitor[n[img[:,22]==1], bitorTable[2]];
    n[img[:,10]==1] = bitor[n[img[:,10]==1], bitorTable[3]];
    n[img[:,13]==1] = bitor[n[img[:,13]==1], bitorTable[4]];
    n[img[:,20]==1] = bitor[n[img[:,20]==1], bitorTable[5]];
    n[img[:,23]==1] = bitor[n[img[:,23]==1], bitorTable[6]];
    n[img[:,11]==1] = bitor[n[img[:,11]==1], bitorTable[7]];
    eulerChar = eulerChar + LUT[n];
    # Octant NEU
    n[:]=1;#restore to ones
    n[img[:,21]==1] = bitor[n[img[:,21]==1], bitorTable[1]];
    n[img[:,24]==1] = bitor[n[img[:,24]==1], bitorTable[2]];
    n[img[:,20]==1] = bitor[n[img[:,20]==1], bitorTable[3]];
    n[img[:,23]==1] = bitor[n[img[:,23]==1], bitorTable[4]];
    n[img[:,12]==1] = bitor[n[img[:,12]==1], bitorTable[5]];
    n[img[:,15]==1] = bitor[n[img[:,15]==1], bitorTable[6]];
    n[img[:,11]==1] = bitor[n[img[:,11]==1], bitorTable[7]];
    eulerChar = eulerChar + LUT[n];
    # Octant SWB
    n[:]=1;#restore to ones
    n[img[:, 7]==1] = bitor[n[img[:, 7]==1], bitorTable[1]];
    n[img[:,16]==1] = bitor[n[img[:,16]==1], bitorTable[2]];
    n[img[:, 8]==1] = bitor[n[img[:, 8]==1], bitorTable[3]];
    n[img[:,17]==1] = bitor[n[img[:,17]==1], bitorTable[4]];
    n[img[:, 4]==1] = bitor[n[img[:, 4]==1], bitorTable[5]];
    n[img[:,13]==1] = bitor[n[img[:,13]==1], bitorTable[6]];
    n[img[:, 5]==1] = bitor[n[img[:, 5]==1], bitorTable[7]];
    eulerChar = eulerChar + LUT[n];
    # Octant SEB
    n[:]=1;#restore to ones
    n[img[:, 9]==1] = bitor[n[img[:, 9]==1], bitorTable[1]];
    n[img[:, 8]==1] = bitor[n[img[:, 8]==1], bitorTable[2]];
    n[img[:,18]==1] = bitor[n[img[:,18]==1], bitorTable[3]];
    n[img[:,17]==1] = bitor[n[img[:,17]==1], bitorTable[4]];
    n[img[:, 6]==1] = bitor[n[img[:, 6]==1], bitorTable[5]];
    n[img[:, 5]==1] = bitor[n[img[:, 5]==1], bitorTable[6]];
    n[img[:,15]==1] = bitor[n[img[:,15]==1], bitorTable[7]];
    eulerChar = eulerChar + LUT[[n]];
    # Octant NWB
    n[:]=1;#restore to ones
    n[img[:, 1]==1] = bitor[n[img[:, 1]==1], bitorTable[1]];
    n[img[:,10]==1] = bitor[n[img[:,10]==1], bitorTable[2]];
    n[img[:, 4]==1] = bitor[n[img[:, 4]==1], bitorTable[3]];
    n[img[:,13]==1] = bitor[n[img[:,13]==1], bitorTable[4]];
    n[img[:, 2]==1] = bitor[n[img[:, 2]==1], bitorTable[5]];
    n[img[:,11]==1] = bitor[n[img[:,11]==1], bitorTable[6]];
    n[img[:, 5]==1] = bitor[n[img[:, 5]==1], bitorTable[7]];
    eulerChar = eulerChar + LUT[n];
    # Octant NEB
    n[:]=1;#restore to ones
    n[img[:, 3]==1] = bitor[n[img[:, 3]==1], bitorTable[1]];
    n[img[:, 2]==1] = bitor[n[img[:, 2]==1], bitorTable[2]];
    n[img[:,12]==1] = bitor[n[img[:,12]==1], bitorTable[3]];
    n[img[:,11]==1] = bitor[n[img[:,11]==1], bitorTable[4]];
    n[img[:, 6]==1] = bitor[n[img[:, 6]==1], bitorTable[5]];
    n[img[:, 5]==1] = bitor[n[img[:, 5]==1], bitorTable[6]];
    n[img[:,15]==1] = bitor[n[img[:,15]==1], bitorTable[7]];
    eulerChar = eulerChar + LUT[n];
    
    EulerInv= eulerChar==0 ;
    return EulerInv
end

function p_is_simple(N)
    #is_simple=logical, N=logical
    
    # copy neighbors for labeling
    n_p = size(N,1);
    # is_simple = logical(ones(n_p, 1)); ##ok<LOGL>
    is_simple=fill(true,n_p, 1)
    
    cube = zeros(Int8,n_p, 26);
    cube[:, 1:13]=N[:, 1:13];
    cube[:, 14:26]=N[:,15:27];
    
    label = 2*(ones(Int8,n_p, 1));
    
    # for all points in the neighborhood
    for i=1:26
        
        idx = cube[:,i]==1 & is_simple;
        
        if sum(idx .!=0)>0
            
            # start recursion with any octant that contains the point i
            if i==1 ||i==2||i==4||i==5||i==10||i==11||i==13
                cube[idx,:] = p_oct_label(1, label, cube[idx,:] );
            elseif i==3 ||i==6||i==12||i==14
                cube[idx,:] = p_oct_label(2, label, cube[idx,:] );
            elseif i==7 ||i==8||i==15||i==16
                cube[idx,:] = p_oct_label(3, label, cube[idx,:] );
            elseif i==9 ||i==17
                cube[idx,:] = p_oct_label(4, label, cube[idx,:] );
            elseif i==18 ||i==19||i==21||i==22
                cube[idx,:] = p_oct_label(5, label, cube[idx,:] );
            elseif i==20 ||i==23
                cube[idx,:] = p_oct_label(6, label, cube[idx,:] );
            elseif i==24 ||i==25
                cube[idx,:] = p_oct_label(7, label, cube[idx,:] );
            elseif i==26
                cube[idx,:] = p_oct_label(8, label, cube[idx,:] );
            end
    
            label[idx] = label[idx]+1;
            del_idx = label>=4;
            
            if sum(del_idx .!=0)>0
                is_simple[del_idx] = false;
            end
        end
    end

    return is_simple
end

function p_oct_label(octant, label, cube)
    #cube=uint8, octant=double integer, label=uint8, cube=uint8
    
    # check if there are points in the octant with value 1
    if ( octant==1 )
        
        # set points in this octant to current label
        # and recurseive labeling of adjacent octants
        idx = cube[:,1] == 1;
        if sum(idx .!=0)>0
            cube[idx,1] = label[idx];
        end
        
        idx = cube[:,2] == 1;
        if sum(idx .!=0)>0
            cube[idx,2] = label[idx];
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,4] == 1;
        if sum(idx .!=0)>0
            cube[idx,4] = label[idx];
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,5] == 1;
        if sum(idx .!=0)>0
            cube[idx,5] = label[idx];
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,10] == 1;
        if sum(idx .!=0)>0
            cube[idx,10] = label[idx];
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,11] == 1;
        if sum(idx .!=0)>0
            cube[idx,11] = label[idx];
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,13] == 1;
        if sum(idx .!=0)>0
            cube[idx,13] = label[idx];
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
        end
        
    end
    
    if ( octant==2 )
        
        idx = cube[:,2] == 1;
        if sum(idx .!=0)>0
            cube[idx,2] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,5] == 1;
        if sum(idx .!=0)>0
            cube[idx,5] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,11] == 1;
        if sum(idx .!=0)>0
            cube[idx,11] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,3] == 1;
        if sum(idx .!=0)>0
            cube[idx,3] = label[idx];
        end
    
        idx = cube[:,6] == 1;
        if sum(idx .!=0)>0
            cube[idx,6] = label[idx];
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,12] == 1;
        if sum(idx .!=0)>0
            cube[idx,12] = label[idx];
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,14] == 1;
        if sum(idx .!=0)>0
            cube[idx,14] = label[idx];
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
    
    end
    
    if( octant==3 )
        
        idx = cube[:,4] == 1;
        if sum(idx .!=0)>0
            cube[idx,4] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,5] == 1;
        if sum(idx .!=0)>0
            cube[idx,5] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,13] == 1;
        if sum(idx .!=0)>0
            cube[idx,13] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,7] == 1;
        if sum(idx .!=0)>0
            cube[idx,7] = label[idx];
        end
    
        idx = cube[:,8] == 1;
        if sum(idx .!=0)>0
            cube[idx,8] = label[idx];
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,15] == 1;
        if sum(idx .!=0)>0
            cube[idx,15] = label[idx];
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,16] == 1;
        if sum(idx .!=0)>0
            cube[idx,16] = label[idx];
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
        
    end
    
    if( octant==4 )
        
        idx = cube[:,5] == 1;
        if sum(idx .!=0)>0
            cube[idx,5] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,6] == 1;
        if sum(idx .!=0)>0
            cube[idx,6] = label[idx];
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,14] == 1;
        if sum(idx .!=0)>0
            cube[idx,14] = label[idx];
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,8] == 1;
        if sum(idx .!=0)>0
            cube[idx,8] = label[idx];
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,16] == 1;
        if sum(idx .!=0)>0
            cube[idx,16] = label[idx];
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,9] == 1;
        if sum(idx .!=0)>0
            cube[idx,9] = label[idx];
        end
    
        idx = cube[:,17] == 1;
        if sum(idx .!=0)>0
            cube[idx,17] = label[idx];
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
    
    end
    
    if ( octant==5 )
        
        idx = cube[:,10] == 1;
        if sum(idx .!=0)>0
            cube[idx,10] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,11] == 1;
        if sum(idx .!=0)>0
            cube[idx,11] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,13] == 1;
        if sum(idx .!=0)>0
            cube[idx,13] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,18] == 1;
        if sum(idx .!=0)>0
            cube[idx,18] = label[idx];
        end
    
        idx = cube[:,19] == 1;
        if sum(idx .!=0)>0
            cube[idx,19] = label[idx];
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,21] == 1;
        if sum(idx .!=0)>0
            cube[idx,21] = label[idx];
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,22] == 1;
        if sum(idx .!=0)>0
            cube[idx,22] = label[idx];
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
    
    end
    
    if( octant==6 )
        
        idx = cube[:,11] == 1;
        if sum(idx .!=0)>0
            cube[idx,11] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,12] == 1;
        if sum(idx .!=0)>0
            cube[idx,12] = label[idx];
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,14] == 1;
        if sum(idx .!=0)>0
            cube[idx,14] = label[idx];
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,19] == 1;
        if sum(idx .!=0)>0
            cube[idx,19] = label[idx];
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
        end
    
    
        idx = cube[:,22] == 1;
        if sum(idx .!=0)>0
            cube[idx,22] = label[idx];
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,20] == 1;
        if sum(idx .!=0)>0
            cube[idx,20] = label[idx];
        end
    
        idx = cube[:,23] == 1;
        if sum(idx .!=0)>0
            cube[idx,23] = label[idx];
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
    
    end
    
    if( octant==7 )
        
        idx = cube[:,13] == 1;
        if sum(idx .!=0)>0
            cube[idx,13] = label[idx];
            cube[idx,:] = p_oct_label(1,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,15] == 1;
        if sum(idx .!=0)>0
            cube[idx,15] = label[idx];
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,16] == 1;
        if sum(idx .!=0)>0
            cube[idx,16] = label[idx];
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(8,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,21] == 1;
        if sum(idx .!=0)>0
            cube[idx,21] = label[idx];
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,22] == 1;
        if sum(idx .!=0)>0
            cube[idx,22] = label[idx];
            cube[idx,:] = p_oct_label([5,label[idx],cube[idx,:]]);
            cube[idx,:] = p_oct_label([6,label[idx],cube[idx,:]]);
            cube[idx,:] = p_oct_label([8,label[idx],cube[idx,:]]);
        end
    
        idx = cube[:,24] == 1;
        if sum(idx .!=0)>0
            cube[idx,24] = label[idx];
        end
        
        idx = cube[:,25] == 1;
        if sum(idx .!=0)>0
            cube[idx,25] = label[idx];
            cube[idx,:] = p_oct_label([8,label[idx],cube[idx,:]]);
        end
    end
    
    if ( octant==8 )
        
        idx = cube[:,14] == 1;
        if sum(idx .!=0)>0
            cube[idx,14] = label[idx];
            cube[idx,:] = p_oct_label(2,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
        end
    
        idx = cube[:,16] == 1;
        if sum(idx .!=0)>0
            cube[idx,16] = label[idx];
            cube[idx,:] = p_oct_label(3,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,17] == 1;
        if sum(idx .!=0)>0
            cube[idx,17] = label[idx];
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,22] == 1;
        if sum(idx .!=0)>0
            cube[idx,22] = label[idx];
            cube[idx,:] = p_oct_label(5,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,17] == 1;
        if sum(idx .!=0)>0
            cube[idx,17] = label[idx];
            cube[idx,:] = p_oct_label(4,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,23] == 1;
        if sum(idx .!=0)>0
            cube[idx,23] = label[idx];
            cube[idx,:] = p_oct_label(6,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,25] == 1;
        if sum(idx .!=0)>0
            cube[idx,25] = label[idx];
            cube[idx,:] = p_oct_label(7,label[idx],cube[idx,:]);
        end
        
        idx = cube[:,26] == 1;
        if sum(idx .!=0)>0
            cube[idx,26] = label[idx];
        end
    end
    return cube
end

function pk_get_nh(img,i,x,y,z)
    #nhood=logical, img=logical, i=double integer
    sz=size(img);
    # x,y,z=ind2sub(sz,i);

    # nhood = logical(zeros(length(i),27)); ##ok<LOGL>
    nhood=fill(false,length(i),27)

    for xx=1:3
        for yy=1:3
            for zz=1:3
                w=sub2ind( [3 3 3],xx,yy,zz);
                idx = sub2ind(sz,x+xx-2,y+yy-2,z+zz-2);
                nhood[:,w]=img[idx];
            end
        end
    end
    return nhood
end



function Skeleton3D(skel,spare)
    # # SKELETON3D Calculate the 3D skeleton of an arbitrary binary volume using parallel medial axis thinning.
    # #
    # # skel = SKELETON3D(img) returns the skeleton of the binary volume 'img'
    # # skel = SKELETON3D(img,mask) preserves foreground voxels in 'mask'
    # #
    # # MATLAB vectorized implementation of the algorithm by Lee, Kashyap and Chu
    # # "Building skeleton models via 3-D medial surface/axis thinning algorithms."
    # # Computer Vision, Graphics, and Image Processing, 56(6):462???478, 1994.
    # #
    # # Inspired by the ITK implementation by Hanno Homann
    # # http://hdl.handle.net/1926/1292
    # # and the Fiji/ImageJ plugin by Ignacio Arganda-Carreras
    # # http://fiji.sc/wiki/index.php/Skeletonize3D
    # #
    # # Philip Kollmannsberger (philipk@gmx.net)
    # #
    # # For more information, see <a
    # # href="matlab:web('http://www.mathworks.com/matlabcentral/fileexchange/43400-skeleton3d')">Skeleton3D</a> at the MATLAB File Exchange.
    
    # #Edited to not require the IPT (by replacing the padarray function) and
    # #to make compatible with ML6.5 (by changing all pre-allocation that involve
    # #the __,'like',p syntax).
    # #Further edits were made to be more strict with data types.
    
    # if ~isa(skel,'logical') || sum(size(skel)>1)~=3
    #     error('first input must be a 3D logical array')
    # end
    # if nargin==2
    #     if ~isa(spare,'logical') || ~isequal(size(spare),size(skel))
    #         error('both inputs must be 3D logical arrays of the same size')
    #     end
    # end
    
    # pad volume with zeros to avoid edge effects
    # tmp=skel;
    # skel=logical(zeros(size(tmp)+[2 2 2] )); ##ok<LOGL>
    # skel[2:end-1,2:end-1,2:end-1]=tmp;
    # if (nargin==2)
    #     tmp=spare;
    #     spare=logical(zeros(size(tmp)+[2 2 2] )); ##ok<LOGL>
    #     spare[2:end-1,2:end-1,2:end-1]=tmp;
    # end

    skel=padarray(skel,Pad(1, 1,1));

    # if (nargin==2)
    #     spare=padarray(spare,[1 1 1]);
    # end

    
    # if (nargin==2)
    #     spare=padarray(spare,Pad(1, 1,1));
    # end


    # fill lookup table
    eulerLUT = FillEulerLUT();

    width = size(skel,1);
    height = size(skel,2);
    depth = size(skel,3);

    unchangedBorders = 0;
        
    while (unchangedBorders < 6 )  # loop until no change for all six border types
        unchangedBorders = 0;
        for currentBorder=1:6 # loop over all 6 directions
            cands=fill(false,width,height,depth); ##ok<LOGL>
            # cands= false(width,height,depth, "like", skel);
            if currentBorder==4
                x=2:size(skel,1); # identify border voxels as candidates
                cands[x,:,:]=((skel[x,:,:] - skel[x.-1,:,:]) .!= 0);
            elseif currentBorder==3
                x=1:size(skel,1)-1;
                cands[x,:,:]=((skel[x,:,:] - skel[x.+1,:,:]) .!= 0);
            elseif currentBorder==1
                y=2:size(skel,2);
                cands[:,y,:]=((skel[:,y,:] - skel[:,y.-1,:]) .!= 0);
            elseif currentBorder==2
                y=1:(size(skel,2)-1);
                cands[:,y,:]=((skel[:,y,:] - skel[:,y.+1,:]) .!= 0);
            elseif currentBorder==6
                z=2:size(skel,3);
                cands[:,:,z]=((skel[:,:,z] - skel[:,:,z.-1]) .!= 0);
            elseif currentBorder==5
                z=1:size(skel,3)-1;
                cands[:,:,z]=((skel[:,:,z] - skel[:,:,z.+1]) .!= 0);
            end
            
            # if excluded voxels were passed, remove them from candidates
            # if (nargin==2)
            #     cands = cands & ~spare;
            # end
            
            # # make sure all candidates are indeed foreground voxels
            # cands = cands[:] & skel[:];
            # make sure all candidates are indeed foreground voxels
            cands[:]=cands[:] .==1 .& skel[:] .==1
            cands1=cands[:] .==1 .& skel[:] .==1
            
            noChange = true;
                        
            if sum(cands1)>0
                cands = findall(n -> n != 0, cands)
                x=map(x->x[1],cands)
                y=map(x->x[2],cands)
                z=map(x->x[3],cands)
                
                # get 26-neighbourhood of candidates in volume
                nhood = pk_get_nh(skel,cands,x,y,z);
                
                # remove all endpoints (exactly one nb) from list
                di1 = sum(nhood,dims=2)==2;
                nhood[di1,:]=[];
                cands[di1]=[];
                x[di1]=[];
                y[di1]=[];
                z[di1]=[];

                
                # remove all non-Euler-invariant points from list
                di2 = p_EulerInv(nhood, eulerLUT);
                nhood[di2,:]=[];
                cands[di2]=[];
                x[di2]=[];
                y[di2]=[];
                z[di2]=[];
                
                # remove all non-simple points from list
                di3 = p_is_simple(nhood);
                x[di3]=[];
                y[di3]=[];
                z[di3]=[];
                
                
                # if any candidates left: divide into 8 independent subvolumes
                if (isempty(x))
                    x1 = logical(mod(x,2));
                    x2 = ~x1;
                    y1 = logical(mod(y,2));
                    y2 = ~y1;
                    z1 = logical(mod(z,2));
                    z2 = ~z1;
                    ilst[1].l = x1 & y1 & z1;
                    ilst[2].l = x2 & y1 & z1;
                    ilst[3].l = x1 & y2 & z1;
                    ilst[4].l = x2 & y2 & z1;
                    ilst[5].l = x1 & y1 & z2;
                    ilst[6].l = x2 & y1 & z2;
                    ilst[7].l = x1 & y2 & z2;
                    ilst[8].l = x2 & y2 & z2;
                    
                    #do parallel re-checking for all points in each subvolume
                    # With euler criteria and simple point --> keep topology
                    for i = 1:8                    
                        if any(ilst[i].l)
                            idx = ilst[i].l;
                            li = sub2ind( [width height depth],x[idx],y[idx],z[idx] );
                                                  
                            skel[li]=0; # remove points
                            nh = pk_get_nh(skel,li,x,y,z);
                            di_rc_euler = p_EulerInv(nh, eulerLUT);
                            di_rc = p_is_simple(nh);                  
    
                             # Check for the simple plint criteria
                            no_single_change = true ;
                            if sum(di_rc)>0 # if topology changed: revert
                                skel[li[di_rc]] = true;
                                no_single_change = false;
                            end
                            
                            # Check for the euler criteria
                            if sum(di_rc_euler)>0
                                skel[li[di_rc_euler]] = true;
                                no_single_change = false;
                            end
                   
                            if no_single_change
                                noChange = false; # at least one voxel removed
                            end
                        end
                    end
                end
            end
            
            if (noChange)
                unchangedBorders = unchangedBorders + 1;
            end
            
        end
        

    end
    
    # get rid of padded zeros
    skel = skel[2:end-1,2:end-1,2:end-1];
    return skel
end