kineBox/example_starauto.m at master · adsoto/kineBox · 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
151
function example_starauto
% Sample code for running automatic kinematics of sea star tube feet
% Working with Andres' blue led images


%% Preliminaries

% Extract directories
paths = givePaths;

frInterval = 500;

% Number of frames used for mean image
numFrames = 50;

numThresh = 4;

% % Path to sample zebrafish video
% vid_path = [paths.vid_root filesep 'Seastars' filesep 'Star Prints' ...
%             filesep 'C0423'];

% Path to sample zebrafish video
vid_path = [paths.vid_root filesep 'Seastars' filesep 'Star Prints' ...
            filesep 'C0423'];


% Load video info (v)
v = defineVidObject(vid_path);


%% Prompt user for input

% Read next image
im = getFrame(vid_path,v,1);

% Convert to grayscale, enhance contrast
%imGray = imadjust(adapthisteq(rgb2gray(im)));

% Fine approximate threshold value
%tVal = graythresh(imGray);
%thresh = multithresh(imGray,4);
%imSeg = imquantize(imGray,thresh);
%RGB = label2rgb(imSeg);

%imshow(RGB);
% Binary image
%imBW = im2bw(imGray,tVal);


frSkip = floor(frInterval/numFrames) - 1;

% Loop thru frames
for i = ceil(frInterval/2):(v.UserData.NumFrames-floor(frInterval/2))

    frames = (ceil(i-frInterval/2)+1):frSkip:floor(i+frInterval/2);

    im = rgb2gray(getFrame(vid_path,v,i));

    imMean = meanImage(vid_path,v,'enhance contrast',frames);

    %im1 = imsubtract(im,imMean);
    %imMean = adapthisteq(stasisImage(vid_path,v,'enhance contrast',frames));

    thresh = multithresh(imMean,numThresh);

    imSeg = imquantize(imMean,thresh);
    %RGB = label2rgb(imSeg);


    for j = 1:numThresh
        imFeet = showFeet(imSeg==j,[30 120]);

    end

%    [centers, radii] = imfindcircles(adapthisteq(imMean),[50 120],'ObjectPolarity','bright', ...
%    'Sensitivity',.9);


    if 1
        imshow(imMean,'InitialMagnification','fit');
        plot(centers(:,1),centers(:,2),'+r')
        %h = viscircles(centers,radii);
        sss=2
    end

    % Update status
    disp(['Done ' num2str(i) ' of ' ...
          num2str((v.UserData.NumFrames-frInterval)) ' frames']);
end

% Display
figure
subplot(1,2,1)
imshow(im,'InitialMagnification','fit');
subplot(1,2,2)
imshow(imGray,'InitialMagnification','fit');

delete(v)


function imOut = showFeet(im,dia_range)

%propDiff

minArea = 0.2*min(dia_range)^2;
minArea2 = min(dia_range)^2;
maxArea = max(dia_range)^2;

se = strel('disk',round(min(dia_range)/4));


im = bwareaopen(im,minArea);

im1 = imdilate(im,se);

im2 = imfill(im1,'holes');

%se2 = strel('disk',round(min(dia_range)/8));

im3 = imerode(im2,se);

im4 = bwareaopen(im3,minArea2);


% Get area and centroid
    props = regionprops(logical(im4),'Centroid','Area',...
                        'MajorAxisLength','MinorAxisLength');

    cCent = [];

    for i = 1:length(props)
        if props(i).Area < maxArea && ...
                props(i).MajorAxisLength
            % Circle centroid
            cCent = [cCent;props(i).Centroid(1) props(i).Centroid(2)];
        end
    end

    % Choose largest blob
    if ~isempty(cCent)
        im5 = bwselect(im4,cCent(:,1),cCent(:,2));
    else
        im5 = im4;
    end

imOut = im5;