spatial_image_analysis.py 文件源码

def inertia_axis(self, labels = None, real = True, verbose=False):
        """
        Return the inertia axis of cells, also called the shape main axis.
        Return 3 (3D-oriented) vectors by rows and 3 (length) values.
        """
        # Check the provided `labels`:
        labels = self.label_request(labels)

        # results
        inertia_eig_vec = []
        inertia_eig_val = []
        N = len(labels); percent=0
        for i,label in enumerate(labels):
            if verbose and i*100/float(N) >= percent: print "{}%...".format(percent),; percent += 10
            if verbose and i+1==N: print "100%"
            slices = self.boundingbox(label, real=False)
            center = copy.copy(self.center_of_mass(label, real=False))
            # project center into the slices sub_image coordinate
            if slices is not None:
                for i,slice in enumerate(slices):
                    center[i] = center[i] - slice.start
                label_image = (self.image[slices] == label)
            else:
                print 'No boundingbox found for label {}'.format(label)
                label_image = (self.image == label)

            # compute the indices of voxel with adequate label
            xyz = label_image.nonzero()
            if len(xyz)==0:
                continue # obviously no reasons to go further !
            coord = coordinates_centering3D(xyz, center)
            # compute the variance-covariance matrix (1/N*P.P^T):
            cov = compute_covariance_matrix(coord)
            # Find the eigen values and vectors.
            eig_val, eig_vec = eigen_values_vectors(cov)
            # convert to real-world units if asked:
            if real:
                for i in xrange(3):
                    eig_val[i] *= np.linalg.norm( np.multiply(eig_vec[i],self._voxelsize) )

            inertia_eig_vec.append(eig_vec)
            inertia_eig_val.append(eig_val)

        if len(labels)==1 :
            return return_list_of_vectors(inertia_eig_vec[0]), inertia_eig_val[0]
        else:
            return self.convert_return(return_list_of_vectors(inertia_eig_vec),labels), self.convert_return(inertia_eig_val,labels)