-
Notifications
You must be signed in to change notification settings - Fork 0
/
ampel.m
93 lines (83 loc) · 3.32 KB
/
ampel.m
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
#!/bin/octave
function[ centroid_red, centroid_blue_transformed, centroid_blue, red, blue, rotation, translation ] = ampel( pos, channel_order, n, method );
## This function "ampel" splices the given table including 3d positions
## etc from images of twocolor fluorescence microscopy into red & blue
## color chanels based on the frame number. It uses "dbscan" to cluster
## the images, and sverages the positions within each
## cluster. Depending on "method" either finds best rigid transform, or
## best affine transform correlating the red and blue sets via least
## squares (singular-value-decomposition).
## @ moritz siegel
global nwd
global rf
global minpts
global dist
## sort red/blue/empty of many frames.
red = blue = empty = zeros( size( pos ));
nb = nr = ne = 1;
for k = 1 : size( pos, 1 )
## exclude blank lines (text).
if ( all( pos( k, : ) == 0 ) )
disp( sprintf( 'warning: line %d is empty; skipping.', k ));
fprintf( rf, 'warning: empty line; skipping.\n' );
continue
endif
## red pill / blue pill?
channel = mod( pos( k, 2 ), 3 );
switch ( channel )
case channel_order(1)
blue( nb, : ) = pos( k, : );
nb = nb + 1;
case channel_order(2)
empty( ne, : ) = pos( k, : );
ne = ne + 1;
case channel_order(3)
red( nr, : ) = pos( k, : );
nr = nr + 1;
endswitch
endfor
disp( sprintf( 'sorted localisations:\n %d red\n %d blue\n %d empty', nr, nb, ne ) );
fprintf( rf, 'sorted localisations:\n %d red\n %d blue\n %d empty\n', nr, nb, ne );
## analyse clusters & average centroids.
[ assignments_blue, c_blue ] = dbscan( blue( 1:300, 3:4 ), minpts, dist );
for k = 1 : c_blue
idx = find( assignments_blue == k );
centroid_blue( k, 1:2 ) = mean( blue( idx, 3:4 ), 1 );
endfor
[ assignments_red, c_red ] = dbscan( red( 1:300, 3:4 ), minpts, dist );
for k = 1 : c_red
idx = find( assignments_red == k );
centroid_red( k, 1:2 ) = mean( red( idx, 3:4 ), 1 );
endfor
## need to ditch clusters if the sets are not equally sized.
assert( c_red == c_blue, "clusters are not equally sized" );
## recover transform.
switch( method )
case "rigid"
## we only have 2d data currently, introduce linear dependent z component.
centroid_red = [ centroid_red, zeros( c_red, 1 ) ]
centroid_blue = [ centroid_blue, zeros( c_blue, 1 ) ]
[ rotation, translation ] = rig( centroid_blue, centroid_red )
case "affine"
## we only have 2d data currently, introduce noise for z component,
## matrix cant be singular, else cholesky decomposition fails.
z = 1 + 1e-3 * rand( c_red, 1 );
centroid_red = [ centroid_red, z ];
centroid_blue = [ centroid_blue, z ];
[ rotation, translation ] = affine( centroid_blue, centroid_red );
endswitch
## transform blue set to red set.
centroid_blue_transformed = ( rotation * centroid_blue' )' + translation;
## save matrices.
swd = sprintf( "%s/position%d", nwd, n ); # concat working dir
mkdir( swd )
chdir( swd )
save red.mat red
save blue.mat blue
save empty.mat empty
save centroid_red.mat centroid_red
save centroid_blue.mat centroid_blue
save rotation.mat rotation
save translation.mat translation
save centroid_blue_transformed.mat centroid_blue_transformed
endfunction