codebase / stimlib

Commit 08ada8c2fbcefbfe7f9ff8b94d6db7deb93be8d1

Authored by David Mayerich
1 parent 42ddd92c

updated image to allow ID index access

Showing 6 changed files with 275 additions and 59 deletions   Show diff stats
Inline Side-by-side
python/enviProcess.py deleted
View file @42ddd92
1   -#!/usr/bin/python3
2   -
3   -#import system processes
4   -import subprocess, sys
5   -
6   -if len(sys.argv) > 1:
7   - infile = int(sys.argv[1])
8   -
9   -basefile = infile + "-base"
10   -normfile = infile + "-norm"
11   -
12   -runcommand = "hsiproc " + infile + basefile + " --baseline baseline.txt"
13   -subprocess.call(runcommand, shell=True)
14 0 \ No newline at end of file
python/structen.py
  Show/Hide comments   View file @08ada8c
... ... @@ -5,23 +5,26 @@ Created on Sun Mar 12 21:54:40 2017
5 5 @author: david
6 6 """
7 7  
8   -import numpy
  8 +import numpy as np
  9 +import scipy as sp
9 10 import scipy.ndimage
10 11 import progressbar
11 12 import glob
12 13  
13   -def st2(I, s=1, dtype=numpy.float32):
  14 +def st2(I, s=1, dtype=np.float32):
14 15  
15 16 #calculates the 2D structure tensor for an image using a gaussian window with standard deviation s
16 17  
17 18 #calculate the gradient
18   - dI = numpy.gradient(I)
  19 + dI = np.gradient(I.astype(dtype))
  20 +
  21 + #print(dI[1][20, 164])
19 22  
20 23 #calculate the dimensions of the tensor field
21 24 field_dims = [dI[0].shape[0], dI[0].shape[1], 3]
22 25  
23 26 #allocate space for the tensor field
24   - Tg = numpy.zeros(field_dims, dtype=dtype)
  27 + Tg = np.zeros(field_dims, dtype=dtype)
25 28  
26 29 #calculate the gradient components of the tensor
27 30 ti = 0
... ... @@ -31,7 +34,7 @@ def st2(I, s=1, dtype=numpy.float32):
31 34 ti = ti + 1
32 35  
33 36 #blur the tensor field
34   - T = numpy.zeros(field_dims, dtype=dtype)
  37 + T = np.zeros(field_dims, dtype=dtype)
35 38  
36 39 for i in range(3):
37 40 T[:, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, i], [s, s])
... ... @@ -39,23 +42,18 @@ def st2(I, s=1, dtype=numpy.float32):
39 42  
40 43 return T
41 44  
42   -def st3(I, s=1):
43   - #calculate the structure tensor field for the 3D input image I given the window size s in 3D
  45 +def st3(I, s=1, v=[1, 1, 1], dtype=np.float32):
  46 + #calculate the structure tensor field for the 3D input image I given the window size s and voxel size v
44 47 #check the format for the window size
45   - if type(s) is not list:
46   - s = [s] * 3
47   - elif len(s) == 1:
48   - s = s * 3
49   - elif len(s) == 2:
50   - s.insert(1, s[0])
51 48  
  49 + v = np.array(v)
52 50 print("\nCalculating gradient...")
53   - dI = numpy.gradient(I)
  51 + dI = np.gradient(I.astype(dtype), v[0], v[1], v[2])
54 52 #calculate the dimensions of the tensor field
55 53 field_dims = [dI[0].shape[0], dI[0].shape[1], dI[0].shape[2], 6]
56 54  
57 55 #allocate space for the tensor field
58   - Tg = numpy.zeros(field_dims, dtype=numpy.float32)
  56 + Tg = np.zeros(field_dims, dtype=np.float32)
59 57  
60 58 #calculate the gradient components of the tensor
61 59 ti = 0
... ... @@ -69,13 +67,15 @@ def st3(I, s=1):
69 67 bar.update(ti)
70 68  
71 69 #blur the tensor field
72   - T = numpy.zeros(field_dims, dtype=numpy.float32)
  70 + T = np.zeros(field_dims, dtype=np.float32)
73 71  
74 72 print("\nConvolving tensor field...")
75 73 bar = progressbar.ProgressBar()
76 74 bar.max_value = 6
  75 + sigma = s / v
  76 + print(sigma)
77 77 for i in range(6):
78   - T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], s)
  78 + T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], sigma)
79 79 bar.update(i+1)
80 80  
81 81 return T
... ... @@ -107,13 +107,13 @@ def sym2mat(T):
107 107 n = in_s[T.ndim - 1]
108 108  
109 109 #calculate the dimension of the symmetric matrix
110   - d = int(0.5 * (numpy.sqrt(8. * n + 1.) - 1.))
  110 + d = int(0.5 * (np.sqrt(8. * n + 1.) - 1.))
111 111  
112 112 #calculate the dimensions of the output field
113 113 out_s = list(in_s)[:-1] + [d] + [d]
114 114  
115 115 #allocate space for the output field
116   - R = numpy.zeros(out_s)
  116 + R = np.zeros(out_s)
117 117  
118 118 ni = 0
119 119 for i in range(d):
... ... @@ -125,9 +125,8 @@ def sym2mat(T):
125 125  
126 126 return R
127 127  
128   -def st2vec(S, vector='largest'):
129   - #Create a color image from a 2D or 3D structure tensor slice
130   -
  128 +def st2vec(S, vector=0):
  129 +
131 130 if(S.ndim != 3):
132 131 print("ERROR: a 2D slice is expected")
133 132 return
... ... @@ -137,27 +136,23 @@ def st2vec(S, vector='largest'):
137 136 del(S)
138 137  
139 138 #calculate the eigenvectors and eigenvalues
140   - l, v = numpy.linalg.eig(T)
  139 + l, v = np.linalg.eig(T)
141 140  
142 141 #get the dimension of the tensor field
143 142 d = T.shape[2]
144 143  
145 144 #allocate space for the vector field
146   - V = numpy.zeros([T.shape[0], T.shape[1], 3])
147   -
  145 + V = np.zeros([T.shape[0], T.shape[1], 3])
  146 +
  147 + #arrange the indices for each pixel from smallest to largest eigenvalue
148 148 idx = l.argsort()
149 149  
150 150 for di in range(d):
151   - if vector == 'smallest':
152   - b = idx[:, :, 0] == di
153   - elif vector == 'largest':
154   - b = idx[:, :, d-1] == di
155   - else:
156   - b = idx[:, :, 1] == di
  151 + b = idx[:, :, -1-vector] == di
157 152 V[b, 0:d] = v[b, :, di]
158 153  
159 154 return V
160   -
  155 +
161 156 def loadstack(filemask):
162 157 #Load an image stack as a 3D grayscale array
163 158  
... ... @@ -171,7 +166,7 @@ def loadstack(filemask):
171 166 Z = len(files)
172 167  
173 168 #allocate space for the image stack
174   - M = numpy.zeros([X, Y, Z]).astype('float32')
  169 + M = np.zeros([X, Y, Z]).astype('float32')
175 170  
176 171 #create a progress bar
177 172 bar = progressbar.ProgressBar()
... ... @@ -184,15 +179,16 @@ def loadstack(filemask):
184 179 bar.update(z+1)
185 180 return M
186 181  
187   -def anisotropy(S):
  182 +#calculate the anisotropy of a structure tensor given the tensor field S
  183 +def anisotropy3(S):
188 184  
189 185 Sf = sym2mat(S)
190 186  
191 187 #calculate the eigenvectors and eigenvalues
192   - l, v = numpy.linalg.eig(Sf)
  188 + l, v = np.linalg.eig(Sf)
193 189  
194 190 #store the sorted eigenvalues
195   - ls = numpy.sort(l)
  191 + ls = np.sort(l)
196 192 l0 = ls[:, :, 0]
197 193 l1 = ls[:, :, 1]
198 194 l2 = ls[:, :, 2]
... ... @@ -210,10 +206,38 @@ def anisotropy(S):
210 206 l_hat = (l0 + l1 + l2)/3
211 207 fa_num = (l2 - l_hat) ** 2 + (l1 - l_hat) ** 2 + (l0 - l_hat) ** 2;
212 208 fa_den = l0 ** 2 + l1 ** 2 + l2 ** 2
213   - FA = numpy.sqrt(3./2.) * numpy.sqrt(fa_num) / numpy.sqrt(fa_den)
  209 + FA = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den)
214 210  
215 211 return FA, Cl, Cp, Cs
216 212  
  213 +#calculate the fractional anisotropy
  214 +def fa(S):
  215 + Sf = sym2mat(S)
  216 +
  217 + #calculate the eigenvectors and eigenvalues
  218 + l, v = np.linalg.eig(Sf)
  219 +
  220 + #store the sorted eigenvalues
  221 + ls = np.sort(l)
  222 + l0 = ls[:, :, 0]
  223 + l1 = ls[:, :, 1]
  224 +
  225 + #if this is a 2D tensor, calculate and return the coherence
  226 + if(S.shape[2] == 3):
  227 + C = ((l0 - l1) / (l0 + l1)) ** 2
  228 + return C
  229 +
  230 + #if this is a 3D tensor
  231 + elif(S.shape[2] == 6):
  232 + l2 = ls[:, :, 2]
  233 +
  234 + #calculate the fractional anisotropy
  235 + l_hat = (l0 + l1 + l2)/3
  236 + fa_num = (l2 - l_hat) ** 2 + (l1 - l_hat) ** 2 + (l0 - l_hat) ** 2;
  237 + fa_den = l0 ** 2 + l1 ** 2 + l2 ** 2
  238 + FA = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den)
  239 + return FA
  240 +
217 241 def st2amira(filename, T):
218 242 #generates a tensor field that can be imported into Amira
219 243  
... ... @@ -225,7 +249,7 @@ def st2amira(filename, T):
225 249 # 5 dz dz ----> 5
226 250  
227 251 #swap the 2nd and 3rd tensor components
228   - A = numpy.copy(T)
  252 + A = np.copy(T)
229 253 A[..., 3] = T[..., 2]
230 254 A[..., 2] = T[..., 3]
231 255  
... ... @@ -246,7 +270,7 @@ def resample3(T, s=2):
246 270 s = s * 3
247 271 elif len(s) == 2:
248 272 s.insert(1, s[0])
249   - s = numpy.array(s)
  273 + s = np.array(s)
250 274  
251 275 bar = progressbar.ProgressBar()
252 276 bar.max_value = T.shape[3]
... ... @@ -270,7 +294,7 @@ def color3(prefix, T, vector='largest', aniso=True):
270 294 for z in range(T.shape[2]):
271 295 S = T[:, :, z, :] #get the slice
272 296 V = st2vec(S, vector='smallest') #calculate the vector
273   - C = numpy.absolute(V) #calculate the absolute value
  297 + C = np.absolute(V) #calculate the absolute value
274 298  
275 299 if aniso == True: #if the image is scaled by anisotropy
276 300 FA, Cl, Cp, Cs = anisotropy(S) #calculate the anisotropy of the slice
... ... @@ -279,9 +303,121 @@ def color3(prefix, T, vector='largest', aniso=True):
279 303 elif vector == 'smallest':
280 304 A = Cp
281 305 else: #otherwise just scale by 1
282   - A = numpy.ones(T.shape[0], T.shape[1])
283   - image = C * numpy.expand_dims(A, 3)
  306 + A = np.ones(T.shape[0], T.shape[1])
  307 + image = C * np.expand_dims(A, 3)
284 308  
285 309 filename = prefix + str(z).zfill(4) + ".bmp"
286 310 scipy.misc.imsave(filename, image)
287   - bar.update(z + 1)
288 311 \ No newline at end of file
  312 + bar.update(z + 1)
  313 +
  314 +def st2stack(T, outfile, **kwargs):
  315 + eigenvector = False #initialize the colormap flags to false
  316 + aniso_color = False
  317 + aniso_alpha = False
  318 + #image = False
  319 + aniso_pwr = 1
  320 + cimage_pwr = 1
  321 + aimage_pwr = 1
  322 + anisostretch = 1 #set the contrast stretch parameter
  323 + alpha_channel = False
  324 + alpha_image = False
  325 + color_image = False
  326 +
  327 + for k,v in kwargs.items(): #for each argument
  328 + if(k == "ev"): #if the user wants a colormap based on an eigenvector
  329 + eigenvector = True #set the eigenvector flag to true
  330 + ev = v #save the desired eigenvector
  331 + if(k == "aniso"): #if the user wants to factor in anisotropy
  332 + aniso = True #set the anisotropy flag to true
  333 + aniso_channel = v #save the anisotropy channel
  334 + if(k == "aniso_color"):
  335 + aniso_color = v
  336 + if(k == "aniso_alpha"):
  337 + aniso_alpha = v
  338 + if(k == "apwr"): #if the user wants to amplify the anisotropy
  339 + aniso_pwr = v #set the anisotropy exponent
  340 + if(k == "cipwr"): #if the user specifies the image power
  341 + cimage_pwr = v
  342 + if(k == "aipwr"):
  343 + aimage_pwr = v
  344 + if(k == "alphaimage"):
  345 + Ia = v
  346 + alpha_image = True
  347 + if(k == "colorimage"):
  348 + Ic = v
  349 + color_image = True
  350 + if(k == "anisostretch"):
  351 + anisostretch = v
  352 + if(k == "alpha"):
  353 + alpha_channel = v
  354 +
  355 + bar = progressbar.ProgressBar()
  356 + bar.max_value = T.shape[2]
  357 + for i in range(0, T.shape[2]):
  358 + #for i in range(0, 50):
  359 +
  360 + if(alpha_image or alpha_channel):
  361 + img = np.ones([T.shape[0], T.shape[1], 4])
  362 + else:
  363 + img = np.ones([T.shape[0], T.shape[1], 3])
  364 + if(eigenvector):
  365 + V = st2vec(T[:, :, i], ev) #get the vector field for slice i corresponding to eigenvector ev
  366 + img[:, :, 0:3] = V #update the image with the vector field information
  367 + if(aniso): #if the user is requesting anisotropy be incorporated into the image
  368 + FA, Cl, Cp, Cs = anisotropy(T[:, :, i]) #calculate the anisotropy of the tensors in slice i
  369 + if(aniso_channel == "fa"):
  370 + A = FA
  371 + elif(aniso_channel == "l"):
  372 + A = Cl
  373 + elif(aniso_channel == "p"):
  374 + A = Cp
  375 + else:
  376 + A = Cs
  377 + if(aniso_alpha):
  378 + print("rendering anisotropy to the alpha channel")
  379 + img[:, :, 3] = A ** aniso_pwr * anisostretch
  380 + if(aniso_color):
  381 + print("rendering anisotropy to the color channel")
  382 + img[:, :, 0:3] = img[:, :, 0:3] * np.expand_dims(A ** aniso_pwr, 3) * anisostretch
  383 + if(alpha_image):
  384 + img[:, :, 3] = Ia[:, :, i]/255 ** aimage_pwr
  385 + if(color_image):
  386 + img[:, :, 0:3] = img[:, :, 0:3] * (np.expand_dims(Ic[:, :, i], 3)/255) ** cimage_pwr #multiply the vector field by the image intensity
  387 + #outname = outfile + str(i).zfill(3) + ".bmp" #get the file name to be saved
  388 + outname = outfile.replace("*", str(i).zfill(3))
  389 +
  390 + sp.misc.imsave(outname, np.ndarray.astype(np.abs(img)*255, "uint8")) #save the output image
  391 + bar.update(i+1)
  392 +
  393 +
  394 +#this function takes a 3D image and turns it into a stack of color images based on the structure tensor
  395 +def img2stack(I, outfile, **kwargs):
  396 +
  397 + vs = [1, 1, 1] #set the default voxel size to 1
  398 + w = 5
  399 +
  400 + for k,v in kwargs.items(): #for each argument
  401 + if(k == "voxelsize"): #if the voxel size is specified
  402 + if(len(v) == 1): #if the user just specifies one value
  403 + vs = [v] * 3 #assume that the voxels are isotropic, create a list of 3 v's
  404 + elif(len(v) == 2): #if the user specifies two values
  405 + vs[0] = v[0] #assume that the voxels are isotropic along (x, y) and anisotropic along z
  406 + vs[1] = v[0]
  407 + vs[2] = v[1]
  408 + elif(len(v) == 3):
  409 + vs = v
  410 + if(k == "window"): #if the user specifies a window size
  411 + w = v
  412 +
  413 + T = st3(I, w, vs) #calculate the structure tensor
  414 +
  415 + st2stack(T, outfile, **kwargs)
  416 +
  417 +def stack2stack(infile_mask, outfile, **kwargs):
  418 +
  419 + I = loadstack(infile_mask) #load the file mask
  420 + for k,v in kwargs.items(): #for each argument
  421 + if(k == "ipwr"):
  422 + img = I
  423 +
  424 + img2stack(I, outfile, image=img, **kwargs) #call the function to convert the image to an output ST stack
289 425 \ No newline at end of file
... ...
stim/image/image.h
  Show/Hide comments   View file @08ada8c
... ... @@ -157,6 +157,24 @@ public:
157 157 }
158 158 #endif
159 159  
  160 + //determines if a filename represents a valid file format that can be loaded/saved
  161 + static bool test_filename(std::string f) {
  162 + stim::filename fname = f;
  163 + std::string ext = fname.extension();
  164 +#ifdef USING_OPENCV
  165 + if (ext == "bmp" ||
  166 + ext == "jpg" ||
  167 + ext == "png" ||
  168 + ext == "pbm" ||
  169 + ext == "tif" )
  170 + return true;
  171 +#else
  172 + if (ext == "pbm" || ext == "bmp")
  173 + return true;
  174 +#endif
  175 + return false;
  176 + }
  177 +
160 178 //save a Netpbm file
161 179 void load_netpbm(std::string filename) {
162 180 std::ifstream infile(filename.c_str(), std::ios::in | std::ios::binary); //open an output file
... ... @@ -423,6 +441,11 @@ public:
423 441 return img[idx(x, y, c)];
424 442 }
425 443  
  444 + /// This function returns a pixel reference based on a 1D index into the image
  445 + T& operator()(size_t i) {
  446 + return img[i];
  447 + }
  448 +
426 449 /// Set all elements in the image to a given scalar value
427 450  
428 451 /// @param v is the value used to set all values in the image
... ...
stim/optics/scalarfield.h
  Show/Hide comments   View file @08ada8c
... ... @@ -428,7 +428,8 @@ public:
428 428 }
429 429 }
430 430  
431   - /// Propagate the field along its orthogonal direction by a distance d
  431 + /// Apply a low pass filter preserving all frequencies lower than or equal to "highest"
  432 + // @param highest is the highest frequency passed
432 433 void lowpass(T highest){
433 434 cpu_scalar_lowpass(E, E, X.len(), Y.len(), highest, R[0], R[1]);
434 435 }
... ...
stim/visualization/camera.h
  Show/Hide comments   View file @08ada8c
... ... @@ -3,6 +3,7 @@
3 3 #include <stim/math/matrix_sq.h>
4 4  
5 5 #include <ostream>
  6 +#include <iostream>
6 7  
7 8 #ifndef STIM_CAMERA_H
8 9 #define STIM_CAMERA_H
... ... @@ -71,6 +72,7 @@ public:
71 72 delta = focus;
72 73  
73 74 focus -= delta;
  75 + std::cout << focus << std::endl;
74 76  
75 77 Dolly(d*delta);
76 78 }
... ...
stim/visualization/glut_template.cpp
  Show/Hide comments   View file @08ada8c
... ... @@ -11,6 +11,72 @@
11 11 #include <stim/visualization/camera.h>
12 12 #include <GL/glut.h>
13 13  
  14 +void render_scene();
  15 +
  16 +void draw_colorcube() {
  17 + glBegin(GL_LINES); //start a line series
  18 + glColor3f(0, 0, 0);
  19 + glVertex3f(-0.5, -0.5, -0.5);
  20 + glColor3f(1, 0, 0);
  21 + glVertex3f(0.5, -0.5, -0.5);
  22 +
  23 + glColor3f(0, 0, 0);
  24 + glVertex3f(-0.5, -0.5, -0.5);
  25 + glColor3f(0, 1, 0);
  26 + glVertex3f(-0.5, 0.5, -0.5);
  27 +
  28 + glColor3f(0, 0, 0);
  29 + glVertex3f(-0.5, -0.5, -0.5);
  30 + glColor3f(0, 0, 1);
  31 + glVertex3f(-0.5, -0.5, 0.5);
  32 +
  33 + glColor3f(1, 1, 1);
  34 + glVertex3f(0.5, 0.5, 0.5);
  35 + glColor3f(1, 1, 0);
  36 + glVertex3f(0.5, 0.5, -0.5);
  37 +
  38 + glColor3f(1, 1, 1);
  39 + glVertex3f(0.5, 0.5, 0.5);
  40 + glColor3f(1, 0, 1);
  41 + glVertex3f(0.5, -0.5, 0.5);
  42 +
  43 + glColor3f(1, 1, 1);
  44 + glVertex3f(0.5, 0.5, 0.5);
  45 + glColor3f(0, 1, 1);
  46 + glVertex3f(-0.5, 0.5, 0.5);
  47 +
  48 + glColor3f(1, 0, 0);
  49 + glVertex3f(0.5, -0.5, -0.5);
  50 + glColor3f(1, 1, 0);
  51 + glVertex3f(0.5, 0.5, -0.5);
  52 +
  53 + glColor3f(1, 0, 0);
  54 + glVertex3f(0.5, -0.5, -0.5);
  55 + glColor3f(1, 0, 1);
  56 + glVertex3f(0.5, -0.5, 0.5);
  57 +
  58 + glColor3f(0, 1, 0);
  59 + glVertex3f(-0.5, 0.5, -0.5);
  60 + glColor3f(1, 1, 0);
  61 + glVertex3f(0.5, 0.5, -0.5);
  62 +
  63 + glColor3f(0, 1, 0);
  64 + glVertex3f(-0.5, 0.5, -0.5);
  65 + glColor3f(0, 1, 1);
  66 + glVertex3f(-0.5, 0.5, 0.5);
  67 +
  68 + glColor3f(0, 0, 1);
  69 + glVertex3f(-0.5, -0.5, 0.5);
  70 + glColor3f(1, 0, 1);
  71 + glVertex3f(0.5, -0.5, 0.5);
  72 +
  73 + glColor3f(0, 0, 1);
  74 + glVertex3f(-0.5, -0.5, 0.5);
  75 + glColor3f(0, 1, 1);
  76 + glVertex3f(-0.5, 0.5, 0.5);
  77 + glEnd();
  78 +}}
  79 +
14 80 struct glut_template_struct {
15 81 float theta_scale = 0.01f;
16 82 float phi_scale = 0.01f;
... ... @@ -53,10 +119,11 @@ void glut_display() {
53 119 glLoadIdentity(); //set it to the identity matrix
54 120 stim::vec3<float> p = gt.cam.getPosition(); //get the camera parameters
55 121 stim::vec3<float> u = gt.cam.getUp();
56   - stim::vec3<float> d = gt.cam.getDirection();
57   - gluLookAt(p[0], p[1], p[2], d[0], d[1], d[2], u[0], u[1], u[2]); //specify the camera parameters to OpenGL
  122 + stim::vec3<float> c = gt.cam.getLookAt();
  123 + gluLookAt(p[0], p[1], p[2], c[0], c[1], c[2], u[0], u[1], u[2]); //specify the camera parameters to OpenGL
58 124  
59 125 render_axes(); //render the axes if the user requests them
  126 + draw_colorcube();
60 127 glutSwapBuffers(); //swap in the back buffer (double-buffering is used to prevent tearing)
61 128 }
62 129  
... ...