Skip to content
This repository was archived by the owner on Feb 26, 2025. It is now read-only.

Commit 4c82008

Browse files
seiriosseirios
committed
Many improvements to metrics
+ Fix bugs + Update codes + Handle multiple resolutions + Add missing plots
1 parent ac77646 commit 4c82008

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

48 files changed

+556
-131
lines changed

‎flatmap/code/metrics/coverage.py‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,19 @@
1-
import voxcell as vc
1+
from voxcell import VoxelData
22
import numpy as np
33
import sys
44

5-
flatmap_file = sys.argv[1]
6-
rdepth_file = sys.argv[2]
5+
import flatmap_util as fmutil
76

8-
fmap_vd = vc.VoxelData.load_nrrd(flatmap_file)
9-
fmap = fmap_vd.raw
7+
flatmap_nrrd = sys.argv[1]
8+
rdepth_nrrd = sys.argv[2]
109

11-
rdepth_vd = vc.VoxelData.load_nrrd(rdepth_file)
12-
mask = ~np.isnan(rdepth_vd.raw)
10+
fmap_vd, fmask = fmutil.load_flatmap(flatmap_nrrd)
1311

14-
is_unmapped = lambda x: x[0] == -1 or x[1] == -1
12+
rdepth_vd = VoxelData.load_nrrd(rdepth_nrrd)
13+
mask = ~np.isnan(rdepth_vd.raw)
1514

1615
n_total = np.sum(mask)
17-
n_not_mapped = np.sum([is_unmapped(x) for x in fmap[mask]])
16+
n_mapped = np.sum(fmask)
17+
n_not_mapped = n_total - n_mapped
1818

19-
print('{} {}'.format(n_not_mapped,1.0 - n_not_mapped / n_total))
19+
print('{} {}'.format(n_not_mapped,n_mapped / n_total))

‎flatmap/code/metrics/laplacian.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@
1111
rdepth = VoxelData.load_nrrd(rdepth_nrrd)
1212
mask = VoxelData.load_nrrd(mask_nrrd)
1313

14-
lapl = laplace(rdepth.raw,mode='nearest')
14+
lapl = laplace(rdepth.raw, mode='nearest')
1515
lapl[mask.raw != 1] = np.nan # keep only interior values
1616

1717
rdepth.with_data(lapl.astype(np.float32)).save_nrrd(output_nrrd)

‎flatmap/code/metrics/orthogonality.py‎

Lines changed: 12 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -8,21 +8,21 @@
88
norm_eps = 1e-12
99

1010
flatmap_nrrd = sys.argv[1]
11-
orient_x_file = sys.argv[2]
12-
orient_y_file = sys.argv[3]
13-
orient_z_file = sys.argv[4]
11+
orient_x_nrrd = sys.argv[2]
12+
orient_y_nrrd = sys.argv[3]
13+
orient_z_nrrd = sys.argv[4]
1414
mask_nrrd = sys.argv[5]
1515
output_nrrd = sys.argv[6]
1616

1717
# load flat map
18-
fmap_vd = VoxelData.load_nrrd(flatmap_nrrd)
18+
fmap_vd, fmsk = fmutil.load_flatmap(flatmap_nrrd)
1919
fmap = fmap_vd.raw
2020
vshape = fmap.shape[0:3]
2121

2222
# load orientation vector
23-
nx_vd = VoxelData.load_nrrd(orient_x_file)
24-
ny_vd = VoxelData.load_nrrd(orient_y_file)
25-
nz_vd = VoxelData.load_nrrd(orient_z_file)
23+
nx_vd = VoxelData.load_nrrd(orient_x_nrrd)
24+
ny_vd = VoxelData.load_nrrd(orient_y_nrrd)
25+
nz_vd = VoxelData.load_nrrd(orient_z_nrrd)
2626

2727
# load mask
2828
mask_vd = VoxelData.load_nrrd(mask_nrrd)
@@ -31,7 +31,6 @@
3131
# get flat coordinates
3232
fx = fmap[:,:,:,0]
3333
fy = fmap[:,:,:,1]
34-
fmsk = fx > -1
3534

3635
## we use 0 as background value for gradient computation
3736
fx[~fmsk] = 0.0
@@ -50,13 +49,13 @@
5049
dy_z = grad_yflat[2][fmsk]
5150

5251
# compute cross product of X and Y gradients
53-
dx = np.vstack((dx_x,dx_y,dx_z)).T
54-
dy = np.vstack((dy_x,dy_y,dy_z)).T
52+
dx = np.vstack((dx_x, dx_y, dx_z)).T
53+
dy = np.vstack((dy_x, dy_y, dy_z)).T
5554

5655
cross = np.cross(dx,dy)
57-
norms = np.linalg.norm(cross,axis=1)
56+
norms = np.linalg.norm(cross, axis=1)
5857
fmsk_norms = norms > norm_eps
59-
cross[fmsk_norms] = cross[fmsk_norms] / norms[fmsk_norms,None] # unit norm
58+
cross[fmsk_norms] = cross[fmsk_norms] / norms[fmsk_norms, None] # unit norm
6059

6160
# get orientation vectors
6261
nx = nx_vd.raw[fmsk] * np.sign(nx_vd.voxel_dimensions[0]) # X axis direction
@@ -69,7 +68,7 @@
6968
orient[fmsk_norms] = orient[fmsk_norms] / norms[fmsk_norms,None] # unit norm
7069

7170
# compute alignment between cross product and orientation vector
72-
ortho = np.einsum('ij,ij->i',orient,cross) # dot product
71+
ortho = np.einsum('ij,ij->i', orient, cross) # dot product
7372

7473
# setup result image
7574
ortho3d = np.full(vshape, np.nan, dtype=np.float32)

‎flatmap/code/metrics/pairwise_preimage_continuity.py‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,13 @@
1717
if len(sys.argv) > 4:
1818
ncpu = int(sys.argv[4])
1919

20-
# load flat map
21-
fmap_vd = VoxelData.load_nrrd(flatmap_nrrd)
20+
# load and discretize flat map
21+
fmap_vd, fmask = fmutil.load_flatmap(flatmap_nrrd)
22+
fmap_d = fmutil.discretize_flatmap(fmap_vd, fmask, pixel_res)
23+
fmap_vd = fmap_vd.with_data(fmap_d)
24+
del fmask
2225

26+
# get pre-images
2327
voxels = fmutil.get_preimages_vox(fmap_vd, pixel_res)
2428
del fmap_vd # FREE
2529

‎flatmap/code/metrics/per_voxel.R‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
library(dplyr)
2+
library(ggplot2)
3+
4+
if(!exists("argv")) argv <- commandArgs(trailingOnly = T)
5+
argc <- length(argv)
6+
7+
div_txt = argv[1]
8+
lap_txt = argv[2]
9+
ort_txt = argv[3]
10+
out_png = argv[4]
11+
12+
div <- read.table(div_txt, na.strings = c("nan"))
13+
lap <- read.table(lap_txt, na.strings = c("nan"))
14+
ort <- read.table(ort_txt, na.strings = c("nan"))
15+
16+
div$metric <- as.factor("divergence")
17+
lap$metric <- as.factor("laplacian")
18+
ort$metric <- as.factor("orthogonality")
19+
20+
dat <- bind_rows(ort,lap,div)
21+
22+
p <- ggplot(dat, aes(V1, metric)) +
23+
geom_boxplot(aes(colour = metric),
24+
outlier.size = 0.5,
25+
outlier.alpha = 0.01,
26+
outlier.stroke = 0.01) +
27+
xlab("") +
28+
ylab("") +
29+
theme_classic() +
30+
theme(axis.line.y = element_blank(),
31+
axis.ticks.y = element_blank(),
32+
legend.position = "none")
33+
34+
ggsave(out_png, p,
35+
units = "in",
36+
width = 5,
37+
height = 2,
38+
dpi = 300)

‎flatmap/code/metrics/pixel_metric.gnuplot‎

Lines changed: 188 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
1+
#!/gpfs/bbp.cscs.ch/home/bolanos/bin/gnuplot -persist
2+
#
3+
#
4+
# G N U P L O T
5+
# Version 5.4 patchlevel 1 last modified 2021-01-02
6+
#
7+
# Copyright (C) 1986-1993, 1998, 2004, 2007-2020
8+
# Thomas Williams, Colin Kelley and many others
9+
#
10+
# gnuplot home: http://www.gnuplot.info
11+
# mailing list: [email protected]
12+
# faq, bugs, etc: type "help FAQ"
13+
# immediate help: type "help" (plot window: hit 'h')
14+
# set terminal x11 noreplotonresize nopersist
15+
# set output
16+
unset clip points
17+
set clip one
18+
unset clip two
19+
unset clip radial
20+
set errorbars front 1.000000
21+
set border 31 front lt black linewidth 1.000 dashtype solid
22+
set zdata
23+
set ydata
24+
set xdata
25+
set y2data
26+
set x2data
27+
set boxwidth
28+
set boxdepth 0
29+
set style fill empty border
30+
set style rectangle back fc bgnd fillstyle solid 1.00 border lt -1
31+
set style circle radius graph 0.02
32+
set style ellipse size graph 0.05, 0.03 angle 0 units xy
33+
set dummy x, y
34+
set format x "% h"
35+
set format y "% h"
36+
set format x2 "% h"
37+
set format y2 "% h"
38+
set format z "% h"
39+
set format cb "% h"
40+
set format r "% h"
41+
set ttics format "% h"
42+
set timefmt "%d/%m/%y,%H:%M"
43+
set angles radians
44+
set tics back
45+
unset grid
46+
unset raxis
47+
set theta counterclockwise right
48+
set style parallel front lt black linewidth 2.000 dashtype solid
49+
set key notitle
50+
set key fixed left top vertical Right noreverse enhanced autotitle nobox
51+
set key noinvert samplen 4 spacing 1 width 0 height 0
52+
set key maxcolumns 0 maxrows 0
53+
set key noopaque
54+
unset label
55+
unset arrow
56+
unset style line
57+
unset style arrow
58+
set style histogram clustered gap 2 title textcolor lt -1
59+
unset object
60+
unset walls
61+
set style textbox transparent margins 1.0, 1.0 border lt -1 linewidth 1.0
62+
set offsets 0, 0, 0, 0
63+
set pointsize 1
64+
set pointintervalbox 1
65+
set encoding default
66+
unset polar
67+
unset parametric
68+
unset spiderplot
69+
unset decimalsign
70+
unset micro
71+
unset minussign
72+
set view 60, 30, 1, 1
73+
set view azimuth 0
74+
set view equal xyz
75+
set rgbmax 255
76+
set samples 100, 100
77+
set isosamples 10, 10
78+
set surface
79+
unset contour
80+
set cntrlabel format '%8.3g' font '' start 5 interval 20
81+
set mapping cartesian
82+
set datafile separator whitespace
83+
set datafile nocolumnheaders
84+
unset hidden3d
85+
set cntrparam order 4
86+
set cntrparam linear
87+
set cntrparam levels 5
88+
set cntrparam levels auto
89+
set cntrparam firstlinetype 0 unsorted
90+
set cntrparam points 5
91+
set size ratio 0 1,1
92+
set origin 0,0
93+
set style data points
94+
set style function lines
95+
unset xzeroaxis
96+
unset yzeroaxis
97+
unset zzeroaxis
98+
unset x2zeroaxis
99+
unset y2zeroaxis
100+
set xyplane relative 0
101+
set tics scale 1, 0.5, 1, 1, 1
102+
set mxtics default
103+
set mytics default
104+
set mztics default
105+
set mx2tics default
106+
set my2tics default
107+
set mcbtics default
108+
set mrtics default
109+
set nomttics
110+
set xtics border in scale 1,0.5 mirror norotate autojustify
111+
set xtics norangelimit autofreq
112+
set ytics border in scale 1,0.5 mirror norotate autojustify
113+
set ytics norangelimit autofreq
114+
set ztics border in scale 1,0.5 nomirror norotate autojustify
115+
set ztics norangelimit autofreq
116+
unset x2tics
117+
unset y2tics
118+
set cbtics border in scale 1,0.5 mirror norotate autojustify
119+
set cbtics norangelimit autofreq
120+
set rtics axis in scale 1,0.5 nomirror norotate autojustify
121+
set rtics norangelimit autofreq
122+
unset ttics
123+
set title ""
124+
set title font "" textcolor lt -1 norotate
125+
set timestamp bottom
126+
set timestamp ""
127+
set timestamp font "" textcolor lt -1 norotate
128+
set trange [ * : * ] noreverse nowriteback
129+
set urange [ * : * ] noreverse nowriteback
130+
set vrange [ * : * ] noreverse nowriteback
131+
set xlabel "Flat pixel size [arb. unit]"
132+
set xlabel font "" textcolor lt -1 norotate
133+
set x2label ""
134+
set x2label font "" textcolor lt -1 norotate
135+
set xrange [ * : * ] noreverse writeback
136+
set x2range [ * : * ] noreverse writeback
137+
set ylabel "Pre-image radius [um]"
138+
set ylabel font "" textcolor lt -1 rotate
139+
set y2label ""
140+
set y2label font "" textcolor lt -1 rotate
141+
set yrange [ * : * ] noreverse writeback
142+
set y2range [ * : * ] noreverse writeback
143+
set zlabel ""
144+
set zlabel font "" textcolor lt -1 norotate
145+
set zrange [ * : * ] noreverse writeback
146+
set cblabel ""
147+
set cblabel font "" textcolor lt -1 rotate
148+
set cbrange [ * : * ] noreverse writeback
149+
set rlabel ""
150+
set rlabel font "" textcolor lt -1 norotate
151+
set rrange [ * : * ] noreverse writeback
152+
unset logscale
153+
unset jitter
154+
set zero 1e-08
155+
set lmargin -1
156+
set bmargin -1
157+
set rmargin -1
158+
set tmargin -1
159+
set locale "en_US.UTF-8"
160+
set pm3d explicit at s
161+
set pm3d scansautomatic
162+
set pm3d interpolate 1,1 flush begin noftriangles noborder corners2color mean
163+
set pm3d clip z
164+
set pm3d nolighting
165+
set palette positive nops_allcF maxcolors 0 gamma 1.5 color model RGB
166+
set palette rgbformulae 7, 5, 15
167+
set colorbox default
168+
set colorbox vertical origin screen 0.9, 0.2 size screen 0.05, 0.6 front noinvert bdefault
169+
set style boxplot candles range 1.50 outliers pt 7 separation 1 labels auto unsorted
170+
set loadpath
171+
set fontpath
172+
set psdir
173+
set fit nolog brief errorvariables nocovariancevariables errorscaling prescale nowrap v5
174+
GNUTERM = "x11"
175+
I = {0.0, 1.0}
176+
VoxelDistance = 0.0
177+
set border 3
178+
set xtics nomirror out
179+
set ytics nomirror out
180+
set xrange [0:*] noext
181+
set yrange [0:*] noext
182+
## Last datafile plotted: "pixel_metric.txt"
183+
fit a*x infile u (1.0/$1):2:3 zerr via a
184+
set term pdfcairo
185+
set out outfile
186+
plot infile u (1.0/$1):2:3 title '' w yerr pt 4 lc rgbcolor 'red',\
187+
a*x title sprintf("%.0f +/- %.0f x",a,a_err) lc 0
188+
# EOF

‎flatmap/code/metrics/preimage_connectedness.py‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,13 @@
1717
if len(sys.argv) > 4:
1818
ncpu = int(sys.argv[4])
1919

20-
# load flat map
21-
fmap_vd = VoxelData.load_nrrd(flatmap_nrrd)
20+
# load and discretize flat map
21+
fmap_vd, fmask = fmutil.load_flatmap(flatmap_nrrd)
22+
fmap_d = fmutil.discretize_flatmap(fmap_vd, fmask, pixel_res)
23+
fmap_vd = VoxelData(fmap_d, (1, 1))
24+
del fmask
2225

26+
# get pre-images
2327
voxels = fmutil.get_preimages_vox(fmap_vd, pixel_res)
2428
del fmap_vd # FREE
2529

‎flatmap/code/metrics/preimage_geometry.py‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,10 +26,14 @@
2626
if len(sys.argv) > 7:
2727
ncpu = int(argv[7])
2828

29-
# load flat map
30-
fmap_vd = VoxelData.load_nrrd(flatmap_nrrd)
29+
# load and discretize flat map
30+
fmap_vd, fmask = fmutil.load_flatmap(flatmap_nrrd)
3131
vshape = fmap_vd.shape[0:3]
32+
fmap_d = fmutil.discretize_flatmap(fmap_vd, fmask, pixel_res)
33+
fmap_vd = fmap_vd.with_data(fmap_d)
34+
del fmask
3235

36+
# get pre-images
3337
points = fmutil.get_preimages_pts(fmap_vd, pixel_res)
3438
del fmap_vd # FREE
3539

0 commit comments

Comments
 (0)