from_bellaiche.m 3.58 KB
Newer Older
1
2
3
4
5
function obj = from_bellaiche( ...
    cells, ...
    frame, ...
    sides, ...
    vertices, ...
6
    mesh_file_path, ...
7
8
9
10
11
12
13
14
15
    units )
%FROM_BELLAICHE Returns a deproj object built from the results of the tool
%from Yohannes Bellaiche lab.

    %% Get information about source image.

    width = frame.imageSize( 1 );
    pixel_size = frame.scale1D;
    
16
17
18
19
20
    %% Read the mesh file.
    
    fprintf('Loading and processing the tissue mesh file.\n' )
    tic
    
21
    [ V, F ] = deproj.ply_read( mesh_file_path );
22
23
24
25
    si = scatteredInterpolant( V(:,1), V(:,2), V(:,3), 'linear', 'nearest' );
    
    fprintf('Loaded %d vertices in %.1f seconds.\n', size(V, 1), toc )
    
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
    %% Compute local curvature from the height-map.
    
    tic
    fprintf('Computing tissue local curvature.\n' )
    [ H, min_y, min_x ] = deproj.mesh_to_heightmap( V, pixel_size );
    % We need to compute a smoothing scale for H, lest the curvature values
    % will have strong discontinuities along the edges of the mesh.    
    med_area = median( face_areas( V, F ) ); % um^2    
    faces_scale = 2 * sqrt( med_area )/ pi * pixel_size; % pixel units
    [ curvMean, curvGauss, curvK1, curvK2 ] = deproj.compute_curvatures( ...
        H, ...
        faces_scale, ...
        pixel_size );
    
    fprintf('Done in %.1f seconds.\n', toc )
41
    
42
43
44
45
46
47
    %% Build the junction graph.
    
    tic
    fprintf('Building the junction graph.\n' )
    
    edge_table = table();
48
49
50
    
    % Prune duplicate edges.
    edge_table.EndNodes = unique( sides.vertices, 'rows' );
51
52
53
    
    node_table = table();
    jxy = vertices.XYs;
54
55

    jz = si( jxy );
56
57
58
59
60
61
    node_table.Centroid = [ jxy jz ];
    node_table.ID = vertices.numbers;
    
    junction_graph = graph( edge_table, node_table, 'omitselfloops' );
    
    fprintf('Built the junction graph with %d nodes and %d edges in %.1f seconds.\n', ...
62
        junction_graph.numnodes, junction_graph.numedges, toc )
63
64
65
66
    
    %% Convert YB structure to epicells.
    
    n_cells = numel( cells.contour_indices );
67
68
69

    fprintf('Computing morphological descriptors of %d objects.\n', n_cells-1 )
    tic
70
71
72
73
74
75
76
77
78
79
80
81
    
    % We don't go for the first one, which contains the contour of the tissue.
    % Get boundary.
    for i = n_cells : -1 : 2
        
        % Get the boundary.
        % It is about converting linear indices to XY coordinates. We do
        % this by integer division since we know the image width.
        ids = int32( cells.contour_indices{ i } );
        x = idivide( ids, width);
        y = rem( ids, width);
        P = double( [ x y ] );
82
        P = deproj.find_countour( P );
83
        P = P * pixel_size;
84
        center = mean( P );
85
        
86
87
        % Get Z from the mesh.
        z = si( P );
88
89
90
91
92
        boundary = [ P z ];
        
        % Get the junction ids.
        junction_ids = cells.vertices{ i };
        
93
94
95
96
97
        % Get center position with respect to crvature image.
        xc = round( center( 1 ) / pixel_size ) - min_x;
        yc = round( center( 2 ) / pixel_size ) - min_y;
        curvatures = [ curvMean( yc, xc ), curvGauss( yc, xc ), curvK1( yc, xc), curvK2( yc, xc ) ];
        
98
        id = cells.numbers( i );
99
        epicells( i-1 ) = epicell( boundary, junction_ids, id, curvatures );
100
101
102
103
        
        
    end
    
104
105
    fprintf('Done in %.1f seconds.\n', toc )
    
106
107
    obj = deproj( epicells, junction_graph, units );
    
108
    
109
110
111
112
113
114
115
116
117
    
    function  areas = face_areas( V, F )
        
        vec_a = V( F(:, 2), :) - V( F(:, 1), :);
        vec_b = V( F(:, 3), :) - V( F(:, 1), :);
        vec_c = cross( vec_a, vec_b, 2);
        areas = 1/2 * sum( sqrt( sum( vec_c .* vec_c, 2 ) ) );
        
    end
118
119
end