matlab编写L1中值点云骨架提取程序

Question

想用matlab编写L1中值点云骨架提取程序，参考文献是L1-Medial Skeleton of Point Cloud

Answer 1

L1中值算法提取比较全的可以看github的程序，一个是Python的程序：https://github.com/MarcSchotman/skeletons-from-poincloud/blob/master/L1-skeleton.py%EF%BC%8C

import numpy as np
from scipy.spatial import distance
import math, random, sys
from utils import *

def get_thetas(r,h):
    """
    matrix of values r
    """
    thetas =  np.exp((-r**2)/((h/2)**2)) 
    #Clip to JUST not zero
    thetas =  np.clip(thetas, 10**-323, None)
    return thetas

def get_alphas(x,points, h):
    """
    INPUTS:
        x: 1x3 center we of interest, np.ndarray 
        points: Nx3 array of all the points, np.ndarray 
        h: size of local neighboorhood, float 
    """
    r = np.linalg.norm(x - points, axis = 1) + 10**-10
    theta = get_thetas(r, h)

    alphas = theta/r
    return alphas

def get_betas(x,points, h):
    """
    INPUTS:
        x: 1x3 center we of interest, np.ndarray 
        points: Nx3 array of all the points, np.ndarray 
        h: size of local neighboorhood, float 
    """
    r = np.linalg.norm( x - points, axis = 1) + 10**-10
    theta = get_thetas(r,h)

    betas = theta/r**2

    return np.array(betas)

def get_density_weights(points, h0, for_center=False, center = [0,0,0]):
    """
    INPUTS:
        x: 1x3 center we of interest, np.ndarray 
        points: Nx3 array of all the points, np.ndarray 
        h: size of local neighboorhood, float 
    RETURNS:
        - np.array Nx1 of density weights assoiscated to each point
    """
    density_weights = []

    if for_center:
        r = points - center
        r2 = np.einsum('ij,ij->i',r, r)
        density_weights = np.einsum('i->', np.exp((-r2)/((h0/4)**2)))
    else:

        for point in points:
            r = point - points
            r2 = np.einsum('ij,ij->i',r, r)
            #This calculation includes the point itself thus one entry will be zero resultig in the needed + 1 in formula dj = 1+ sum(theta(p_i - p_j))
            density_weight = np.einsum('i->', np.exp((-r2)/((h0/4)**2)))
            density_weights.append(density_weight)

    return np.array(density_weights)

def get_term1(center, points, h, density_weights):
    """
    INPUTS:
        center: 1x3 center we of interest, np.ndarray 
        points: Nx3 array of all the points, np.ndarray 
        h: size of local neighboorhood, float 
        h0: size of first local neighboorhood, float 
    RETURNS:
        - term1 of the equation as float
    """
    t1_t = time.perf_counter()

    r = points - center
    r2 = np.einsum('ij,ij->i',r, r)

    thetas =  np.exp( -r2 / ((h/2)**2)) 
    #Clip to JUST not zero
    # thetas =  np.clip(thetas, 10**-323, None)

    #DIFFERS FROM PAPER
    # r_norm = np.sqrt(r_norm, axis = 1)
    # alphas = thetas/r_norm

    alphas = thetas/density_weights

    denom = np.einsum('i->',alphas)
    if denom  > 10**-20:
        # term1 = np.sum((points.T*alphas).T, axis = 0)/denom
        term1 = np.einsum('j,jk->k',alphas, points) / denom
    else: 
        term1 = np.array(False)
    
    t2_t = time.perf_counter()
    tt = round(t2_t - t1_t, 5)

    return term1, tt

def get_term2(center, centers, h):
    """
    INPUTS:
        center: 1x3 center we of interest, np.ndarray 
        centers: Nx3 array of all the centers (excluding the current center), np.ndarray 
        h: size of local neighboorhood, float 
    RETURNS:
        - term2 of the equation as float
    """
    t1 = time.perf_counter()
    
    x = center - centers
    r2 = np.einsum('ij,ij->i',x, x)
    r = 1/np.sqrt(r2)
    # r3 = np.sum(r**1.2, axis = 1)
    thetas =  np.exp((-r2)/((h/2)**2)) 

    # r_norm = np.linalg.norm(r,axis = 1)
    #DIFFERS FROM PAPER
    #betas =np.einsum('i,i->i', thetas, density_weights)# / r2
    betas = np.einsum('i,i->i',thetas,r)
    
    denom = np.einsum('i->',betas)
    
    if denom > 10**-20:
        num = np.einsum('j,jk->k',betas, x)

        term2 = num/denom
    else:
        term2 = np.array(False)
    
    t2 = time.perf_counter()
    tt = round(t2-t1, 4)
    return term2, tt

def get_sigma(center, centers, h):

    t1 = time.perf_counter()
    #These are the weights
    r = centers - center
    r2 = np.einsum('ij,ij->i',r, r)
    thetas = np.exp((-r2)/((h/2)**2))     

    # thetas = get_thetas(r,h)
    #Thetas are further clipped to a minimum value to prevent infinite covariance
    # weights = np.clip(thetas, 10**-10, None)
    #substract mean then calculate variance\
    cov = np.einsum('j,jk,jl->kl',thetas,r,r)
    # cov = np.zeros((3,3))
    # for index in range(len(r)):
    #     cov += weights[index]*np.outer(r[index],r[index])
    # centers -= np.mean(centers, axis = 0)
    # # print(centers)
    # cov = np.cov(centers.T, aweights=weights)

    #Get eigenvalues and eigenvectors
    values, vectors = np.linalg.eig(cov)

    if np.iscomplex(values).any():
        values = np.real(values)
        vectors = np.real(vectors)
        vectors_norm = np.sqrt(np.einsum('ij,ij->i',vectors, vectors))
        vectors = vectors/vectors_norm

    #Argsort always works from low --> to high so taking the negative values will give us high --> low indices
    sorted_indices = np.argsort(-values)

    values_sorted = values[sorted_indices]
    vectors_sorted = vectors[:,sorted_indices]

    sigma = values_sorted[0]/np.sum(values_sorted)

    t2 = time.perf_counter()

    return sigma,  vectors_sorted, t2-t1

def get_h0(points):

    x_max = points[:,0].max(); x_min = points[:,0].min()

    y_max = points[:,1].max(); y_min = points[:,1].min()

    z_max = points[:,2].max(); z_min = points[:,2].min()
    print("BB values: \n\tx:",x_max - x_min,"\n\ty:",y_max -y_min,"\n\tz:",z_max - z_min)

    diagonal = ((x_max-x_min)**2 + (y_max-y_min)**2+ (z_max-z_min)**2)**.5

    Npoints = len(points)

    return 2*diagonal/(Npoints**(1./3))

class myCenter:
    def __init__(self, center, h, index):
        self.center = center
        self.h = h
        self.label = "non_branch_point"
        self.index = index
        self.connections = []
        self.bridge_connections = None
        self.closest_neighbours = np.array([])
        self.head_tail = False
        self.branch_number = None
        self.eigen_vectors = np.zeros((3,3))
        self.sigma = 0.5

    def set_non_branch(self):
        if self.label != 'branch_point' and self.label !='removed':
            self.set_label('non_branch_point')
            self.connections = []
            self.bridge_connections = None
            self.head_tail = False
            self.branch_number = None

    def set_as_bridge_point(self, key, connection):
        if self.label != 'removed':
            self.set_non_branch()
            self.set_label('bridge_point')
            self.bridge_connections = connection
            self.branch_number = key

    def set_as_branch_point(self, key):
    
        self.connections = []
        self.bridge_connections = None
        self.head_tail = False
        self.branch_number = None
        self.label = 'branch_point'
        self.branch_number = key

    def set_eigen_vectors(self,eigen_vectors):
        if self.label == "non_branch_point":
            self.eigen_vectors = eigen_vectors

    def set_sigma(self,sigma):
        if self.label != "branch_point":
            self.sigma = sigma 

    def set_closest_neighbours(self, closest_neighbours):
        """
        """
        self.closest_neighbours = closest_neighbours

    def set_label(self, label):
        if self.label !='removed':
            self.label = label

    def set_center(self,center):

        if self.label != "branch_point":
            self.center = center

    def set_h(self,h):
        if self.label != "branch_point":
            self.h = h

class myCenters:

    def set_my_non_branch_centers(self):

        my_non_branch_centers = []

        for center in self.myCenters:
            if center.label =='non_branch_point' or center.label == 'bridge_point':
                my_non_branch_centers.append(center)
        self.my_non_branch_centers = my_non_branch_centers

    def get_nearest_neighbours(self):

        distances =distance.squareform(distance.pdist(self.centers ))
        self.closest  = np.argsort(distances, axis =1 )

        for center in self.myCenters:
            # center.set_closest_neighbours(self.closest[center.index,1:])
            closest = self.closest[center.index, :].copy()         
            sorted_local_distances = distances[center.index, closest]**2

            #Returns zero if ALL values are within the range
            in_neighboorhood = np.argmax(sorted_local_distances >= (center.h)**2)
            if in_neighboorhood == 0:
                in_neighboorhood = -1

            center.set_closest_neighbours( closest[1:in_neighboorhood])

    def __init__(self,centers, h0, maxPoints):
        self.centers = centers + 10**-20 #Making sure centers are never the same as the actual points which can lead to bad things
        self.myCenters=[]
        self.my_non_branch_centers=[]
        index = 0
        for center in centers:
            self.myCenters.append(myCenter(center, h, index))
            index+=1
        self.skeleton = {}
        self.closest = []
        self.sigmas = np.array([None] * len(centers))
        self.h0 = h0
        self.h = h0       
        self.eigen_vectors = [None] * len(centers)
        self.branch_points = [None] * len(centers)
        self.non_branch_points = [None] * len(centers)
        self.maxPoints = maxPoints
        self.get_nearest_neighbours()
        self.set_my_non_branch_centers()
        self.Nremoved = 0

        #From the official code
        self.search_distance = .4
        self.too_close_threshold = 0.01
        self.allowed_branch_length = 5


    def remove_centers(self,indices):
        """
        Removes a center completely
        """
        if not isinstance(indices,list):
            indices = list([indices])

        for index in sorted(indices, reverse=True):
            center = self.myCenters[index]
            center.set_label("removed")
            self.centers[center.index] = [9999,9999,9999]
        self.set_my_non_branch_centers()
        self.Nremoved += len(indices)

    def get_non_branch_points(self):

        non_branch_points = []
        for center in self.myCenters:
            if center.label != "branch_point" and center.label != "removed":
                non_branch_points.append(center.index)

        return non_branch_points

    def get_bridge_points(self):

        bridge_points = []
        for key in self.skeleton:
            head = self.skeleton[key]['head_bridge_connection']
            tail = self.skeleton[key]['tail_bridge_connection']

            if head[0] and head[1] != None:
                if not head[1] in bridge_points:
                    bridge_points.append(head[1])
            if tail[0] and tail[1] != None:
                if not tail[1] in bridge_points:
                    bridge_points.append(tail[1])

        return bridge_points

    def update_sigmas(self):

        k = 5
        
        new_sigmas = []
    
        for center in self.my_non_branch_centers:

            index = center.index

            indices = np.array(self.closest[index,:k]).astype(int)
            
            sigma_nearest_k_neighbours = self.sigmas[indices]

            mean_sigma = np.mean(sigma_nearest_k_neighbours)
            
            new_sigmas.append(mean_sigma)
            

        index = 0
        for center in self.my_non_branch_centers:
            
            center.set_sigma(new_sigmas[index])

            self.sigmas[center.index] = new_sigmas[index]
            index +=1

    def update_properties(self):
        
        self.set_my_non_branch_centers()

        for center in self.myCenters:
            index = center.index
            self.centers[index] = center.center
            self.eigen_vectors[index] = center.eigen_vectors
            self.sigmas[index] = center.sigma

        
        self.get_nearest_neighbours()
        self.update_sigmas()

    def update_labels_connections(self):
        """
        Update all the labels of all the centers
            1) goes through all the branches and checks if the head has a bridge connection or a branch connection
                - If bridge connection this is still the head/tail of the branch
                - If it has a branch connection it is simply connected to another branch --> It is no head/tail anymore
            2) Checks if bridges are still bridges
            3) Sets all other points to simple non_branch_points
        """

        updated_centers = []
        for key in self.skeleton:

            branch = self.skeleton[key]

            head = self.myCenters[branch['branch'][0]];  tail = self.myCenters[branch['branch'][-1]]

            #This is either a None value (for not having found a bridge point / connected branch) or this is an integer index
            head_connection = branch['head_bridge_connection'][1]
            tail_connection = branch['tail_bridge_connection'][1]

            if head_connection != None:

                head_connection = self.myCenters[head_connection]

                if branch['head_bridge_connection'][0]:
                    head_connection.set_as_bridge_point(key, head.index)
                    head.head_tail = True
                else:
                    head_connection.set_as_branch_point(key)
                    head.head_tail = False

                head.set_as_branch_point(key)
                head.connections = [head_connection.index, branch['branch'][1]]

                updated_centers.append(head_connection.index)  
                updated_centers.append(head.index)
            else:
                head.set_as_branch_point(key)
                head.head_tail = True
                
            if tail_connection != None:
                
                tail_connection = self.myCenters[tail_connection]

                if branch['tail_bridge_connection'][0]:
                    tail.head_tail = True
                    tail_connection.set_as_bridge_point(key, tail.index)
                else:
                    tail.head_tail = False
                    tail_connection.set_as_branch_point(key)

                tail.set_as_branch_point(key)
                tail.connections = [tail_connection.index, branch['branch'][-2]]
                updated_centers.append(tail_connection.index)
                updated_centers.append(tail.index)
            else:
                tail.set_as_branch_point(key)
                tail.head_tail = True

            # 1) Go through the branch list and set each center t branch_point and set the head_tail value appropriately
            # 2) Set the connections
            index = 1
            for center in branch['branch'][1:-1]:

                center = self.myCenters[center]

                center.set_as_branch_point(key)

                center.connections.append(branch['branch'][index-1])
                center.connections.append(branch['branch'][index+1])
                center.head_tail = False

                updated_centers.append(center.index)
                index+=1

        for center in self.myCenters:

            if center.index in updated_centers:
                continue
            center.set_non_branch()


        for key in self.skeleton:
            branch = self.skeleton[key]

            for index in branch['branch']:
                if branch['branch'].count(index) > 1:
                    print("ERROR: This branch has multiple counts of 1 index...", branch['branch'])
                    break

    def contract(self, points, local_indices, h, density_weights, mu = 0.35):
        """
        Updates the centers by the algorithm suggested in "L1-medial skeleton of Point Cloud 2010"
        INPUT:
            - Centers
            - points belonging to centers
            - local neighbourhood h0
            - mu factor for force between centers (preventing them from clustering)
        OUTPUT: 
            - New centers
            - Sigmas (indicator for the strength of dominant direction)
            - The eigenvectors of the points belonging to the centers
        """
        self.h = h

        t1_total = time.perf_counter(); term1_t = 0; term2_t = 0; sigma_t = 0
        t_pre =0; t_post = 0

        error_center = 0; N = 0;
        for myCenter in self.myCenters:
            
            t1 = time.perf_counter()
            #Get the closest 50 centers to do calculations with
            centers_indices = myCenter.closest_neighbours
            #Get the density weight of these centers
            centers_in = np.array(self.centers[centers_indices])

            my_local_indices = local_indices[myCenter.index]
            local_points = points[my_local_indices]

            t2 = time.perf_counter()
            t_pre += t2-t1
            #Check if we have enough points and centers
            shape = local_points.shape
            if len(shape) ==1:
                continue
            elif shape[0] > 2 and len(centers_in) > 1:                
                
                density_weights_points =density_weights[my_local_indices]

                term1, delta_t1 = get_term1(myCenter.center, local_points,  h, density_weights_points)
            
                term2, delta_t2 = get_term2(myCenter.center, centers_in,  h)
                
                term1_t += delta_t1; term2_t += delta_t2

                if term1.any() and term2.any():

                    sigma, vecs, delta_ts = get_sigma(myCenter.center, centers_in, h)   
                    # sigma = np.clip(sigma, 0 ,1.)          
                    sigma_t += delta_ts

                    #DIFFERS FROM PAPER
                    # mu = mu_length/sigma_length * (sigma - min_sigma)
                    # if mu < 0:
                    #     continue

                    # mu_average +=mu 

                    t1 = time.perf_counter()
                    new_center = term1 + mu*sigma*term2

                    error_center+= np.linalg.norm(myCenter.center - new_center);  N+=1
                    
                    #Update this center object
                    myCenter.set_center(new_center)
                    myCenter.set_eigen_vectors(vecs)
                    myCenter.set_sigma(sigma)
                    myCenter.set_h(h)               
                    t2 = time.perf_counter()

                    t_post += t2- t1
        t2_total = time.perf_counter(); total_time = round(t2_total - t1_total,4);

这边是这里一部分的，安装一下环境改改。
另外一个是原始这个文献上的，要使用linux系统，还有qt比较麻烦。
要考虑原始文章的话还是得用原版的程序，毕竟简短的写个matlab的公式肯定不一样，就不是原本的意思了。