NAS/plotStaticMNMorphology3.m at main · vitorchaud/NAS · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150


function [upperAxes, surface_handle, visualizationBoxLength] = plotStaticMNMorphology3 ( morphologicProp_struct, d_min, d_max, azimute, elevation, baseColor, sectionColor,...
                                                                    figureFileName, geometryType, distanceMetric, isShowSectionOnly, isColoredTrees, showOnlyTree)


    numOfSegments = morphologicProp_struct.('numOfSegments');
    type = morphologicProp_struct.('data').('local').('type');
    parent = morphologicProp_struct.('data').('local').('parent');
    radius = morphologicProp_struct.('data').('local').('radius');
    x = morphologicProp_struct.('data').('local').('x');
    y = morphologicProp_struct.('data').('local').('y');
    z = morphologicProp_struct.('data').('local').('z');

    proximal_position_x = morphologicProp_struct.('data').('local').('proximal_position_x');
    proximal_position_y = morphologicProp_struct.('data').('local').('proximal_position_y');
    proximal_position_z = morphologicProp_struct.('data').('local').('proximal_position_z');

    sectionIndexes =  getCompartmentIndexesBySection5(morphologicProp_struct, d_min, d_max, 0, distanceMetric);

    if isColoredTrees
        zeroOrderSegments = morphologicProp_struct.('data').('zeroOrderSegments');
        stemParent = morphologicProp_struct.('data').('local').('stemParent');
        cmap = hsv(length(zeroOrderSegments));
    end

    visualizationBoxLength = 0;


    for i = 1:numOfSegments

        if type(i)==1 && parent(i)==1 %% Skip false compartments
            continue;
        end

        if isShowSectionOnly && ~sectionIndexes(i)
            continue;
        end


        compartment_pos_x = x(i);
        compartment_pos_y = y(i);
        compartment_pos_z = z(i);

        if(type(i) == 1)
            [Xs,Ys,Zs] = sphere(20); % Compartment is soma
            surface_handle(i) = surf(radius(i)*Xs + compartment_pos_x, radius(i)*Ys + compartment_pos_y, radius(i)*Zs + compartment_pos_z);
            set(surface_handle(i),'FaceColor',[1 0 0],'FaceAlpha',0.5);
            hold on;
        else

            if(abs(compartment_pos_x) > visualizationBoxLength)
                visualizationBoxLength = abs(compartment_pos_x);
            end
            if(abs(compartment_pos_y) > visualizationBoxLength)
                visualizationBoxLength = abs(compartment_pos_y);
            end
            if(abs(compartment_pos_z) > visualizationBoxLength)
                visualizationBoxLength = abs(compartment_pos_z);
            end

            if strcmp(geometryType, 'lines')
                hl = plot3([proximal_position_x(i) compartment_pos_x],...
                    [proximal_position_y(i) compartment_pos_y],...
                    [proximal_position_z(i) compartment_pos_z]);
                set(hl, 'Color', [0.039, 0.039, 0.498]);
            else
                dx = compartment_pos_x - proximal_position_x(i);
                dy = compartment_pos_y - proximal_position_y(i);
                dz = compartment_pos_z - proximal_position_z(i);

                d = sqrt(dx^2 + dy^2 + dz^2);

                % orthogonal vector
                vx = d*dy;
                vy = -d*dx;

                angle = -(180/pi) * acos(dz/d) + 180;

                if strcmp(geometryType, 'cylinder')
                    [X,Y,Z] = cylinder(radius(i), 10);
                end
                if strcmp(geometryType, 'frustum')
                    [X,Y,Z] = cylinder([radius(i) radius(parent(i))], 10);
                end

                surface_handle(i) = surf((X + compartment_pos_x),(Y + compartment_pos_y),(d*Z + compartment_pos_z));
                rotate(surface_handle(i), [vx vy 0], angle, [compartment_pos_x, compartment_pos_y, compartment_pos_z]);
            end
        end

        if ~strcmp(geometryType, 'lines')
            set(surface_handle(i),'FaceColor',baseColor,'EdgeColor', baseColor,'FaceAlpha',0.7);
        end

        if isColoredTrees
            if type(i)==1
                color=[0,0,0];
            else
                if nnz(zeroOrderSegments == stemParent(i))
                    color = cmap(zeroOrderSegments == stemParent(i),:);
                else
                    color = [1 0 0]; % in case the tree was excluded
                end
            end
            set(surface_handle(i),'FaceColor',color,'EdgeColor', color,'FaceAlpha',0.8,'EdgeAlpha',0.8);
        end
        if showOnlyTree
            if zeroOrderSegments(showOnlyTree) == stemParent(i)
                set(surface_handle(i),'FaceColor',color,'EdgeColor', color,'FaceAlpha',0.8,'EdgeAlpha',0.8);
            else
                set(surface_handle(i),'FaceColor',color,'EdgeColor', color,'FaceAlpha',0.005,'EdgeAlpha',0.005);
            end
        end

        if isShowSectionOnly
            if(sectionIndexes(i))
                set(surface_handle(i),'FaceColor',sectionColor,'EdgeColor', sectionColor,'FaceAlpha',0.2);
            end
        end


        hold on
    end


    upperAxes = gca;
    axis square
%     visualizationBoxLength = d_max;
    xlim([-visualizationBoxLength visualizationBoxLength]);
    ylim([-visualizationBoxLength visualizationBoxLength]);
    zlim([-visualizationBoxLength visualizationBoxLength]);
    xlabel('x');
    ylabel('y');
    zlabel('z');

    view(azimute, elevation);

%     frame = getframe(gcf);
%     im = frame2im(frame);
%     [imind,cm] = rgb2ind(im,256);
%     imwrite(imind,cm,figureFileName,'png');


end