property_spatial_image.py 文件源码

def compute_image_property(self, property_name, min_contact_area=None, sub_factor=8.):
        """
        """
        computable_properties = ['barycenter','volume','neighborhood_size','layer','mean_curvature','gaussian_curvature']
        try:
            assert property_name in computable_properties
        except:
            print "Property \""+property_name+"\" can not be computed on image"
            print "Try with one of the following :"
            for p in computable_properties:
                print "  * "+p
        else:
            if self._image is not None:
                if property_name in ['barycenter','volume']:
                    graph = graph_from_image(self._image,background=self.background,spatio_temporal_properties=[property_name],ignore_cells_at_stack_margins=self.ignore_cells_at_stack_margins)
                    property_dict = graph.vertex_property(property_name)
                elif property_name == 'neighborhood_size':
                    neighbors = [self.image_graph.neighbors(l) for l in self.labels]
                    property_dict = dict(zip(self.labels,map(len,neighbors)))
                elif property_name == 'layer':
                    if min_contact_area is None:
                        min_contact_area = self.min_contact_area
                    graph = graph_from_image(self._image,background=self.background,spatio_temporal_properties=['L1'],ignore_cells_at_stack_margins=self.ignore_cells_at_stack_margins, min_contact_area=min_contact_area)
                    first_layer = graph.vertex_property('L1')
                    second_layer_cells = [v for v in graph.vertices() if np.any([first_layer[n] for n in graph.neighbors(v)]) and not first_layer[v]]
                    second_layer = dict(zip(list(graph.vertices()),[v in second_layer_cells for v in graph.vertices()]))
                    property_dict = dict(zip(self.labels,[1 if first_layer[l] else 2 if second_layer[l] else 3 for l in self.labels]))
                elif property_name in ['mean_curvature','gaussian_curvature']:
                    if not self.has_image_property('barycenter'):
                        self.compute_image_property('barycenter')
                    if not self.has_image_property('layer'):
                        print "--> Computing layer property"
                        self.compute_image_property('layer')

                    cell_centers = self.image_property('barycenter')
                    L1_cells = self.labels[self.image_property('layer').values()==1]

                    from openalea.cellcomplex.property_topomesh.utils.implicit_surfaces import implicit_surface_topomesh
                    from openalea.cellcomplex.property_topomesh.property_topomesh_analysis import compute_topomesh_property, compute_topomesh_vertex_property_from_faces
                    from openalea.cellcomplex.property_topomesh.property_topomesh_optimization import property_topomesh_vertices_deformation, topomesh_triangle_split

                    sub_binary_image = (self._image!=self.background).astype(float)[::sub_factor,::sub_factor,::sub_factor]
                    surface_topomesh = implicit_surface_topomesh(sub_binary_image,np.array(sub_binary_image.shape),sub_factor*np.array(self._image.voxelsize),center=False)
                    property_topomesh_vertices_deformation(surface_topomesh,iterations=10)

                    compute_topomesh_property(surface_topomesh,'barycenter',2)
                    compute_topomesh_property(surface_topomesh,'normal',2,normal_method='orientation')

                    compute_topomesh_vertex_property_from_faces(surface_topomesh,'normal',adjacency_sigma=2,neighborhood=5)
                    compute_topomesh_property(surface_topomesh,'mean_curvature',2)
                    compute_topomesh_vertex_property_from_faces(surface_topomesh,property_name,adjacency_sigma=2,neighborhood=5)

                    surface_cells = L1_cells[vq(surface_topomesh.wisp_property('barycenter',0).values(),cell_centers.values(L1_cells))[0]]
                    surface_topomesh.update_wisp_property('label',0,array_dict(surface_cells,list(surface_topomesh.wisps(0))))

                    L1_cell_property = nd.sum(surface_topomesh.wisp_property(property_name,0).values(),surface_cells,index=L1_cells)/nd.sum(np.ones_like(surface_cells),surface_cells,index=L1_cells)
                    L1_cell_property = array_dict(L1_cell_property,L1_cells)    
                    property_dict = array_dict([L1_cell_property[l] if (l in L1_cells) and (not np.isnan(L1_cell_property[l])) else 0 for l in self.labels],self.labels)

                property_data = [property_dict[l] for l in self.labels]
                self.update_image_property(property_name,property_data)