diff --git a/python/enviProcess.py b/python/enviProcess.py deleted file mode 100644 index 1e655b5..0000000 --- a/python/enviProcess.py +++ /dev/null @@ -1,13 +0,0 @@ -#!/usr/bin/python3 - -#import system processes -import subprocess, sys - -if len(sys.argv) > 1: - infile = int(sys.argv[1]) - -basefile = infile + "-base" -normfile = infile + "-norm" - -runcommand = "hsiproc " + infile + basefile + " --baseline baseline.txt" -subprocess.call(runcommand, shell=True) \ No newline at end of file diff --git a/python/structen.py b/python/structen.py index f284ff6..672ee5e 100644 --- a/python/structen.py +++ b/python/structen.py @@ -5,23 +5,26 @@ Created on Sun Mar 12 21:54:40 2017 @author: david """ -import numpy +import numpy as np +import scipy as sp import scipy.ndimage import progressbar import glob -def st2(I, s=1, dtype=numpy.float32): +def st2(I, s=1, dtype=np.float32): #calculates the 2D structure tensor for an image using a gaussian window with standard deviation s #calculate the gradient - dI = numpy.gradient(I) + dI = np.gradient(I.astype(dtype)) + + #print(dI[1][20, 164]) #calculate the dimensions of the tensor field field_dims = [dI[0].shape[0], dI[0].shape[1], 3] #allocate space for the tensor field - Tg = numpy.zeros(field_dims, dtype=dtype) + Tg = np.zeros(field_dims, dtype=dtype) #calculate the gradient components of the tensor ti = 0 @@ -31,7 +34,7 @@ def st2(I, s=1, dtype=numpy.float32): ti = ti + 1 #blur the tensor field - T = numpy.zeros(field_dims, dtype=dtype) + T = np.zeros(field_dims, dtype=dtype) for i in range(3): T[:, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, i], [s, s]) @@ -39,23 +42,18 @@ def st2(I, s=1, dtype=numpy.float32): return T -def st3(I, s=1): - #calculate the structure tensor field for the 3D input image I given the window size s in 3D +def st3(I, s=1, v=[1, 1, 1], dtype=np.float32): + #calculate the structure tensor field for the 3D input image I given the window size s and voxel size v #check the format for the window size - if type(s) is not list: - s = [s] * 3 - elif len(s) == 1: - s = s * 3 - elif len(s) == 2: - s.insert(1, s[0]) + v = np.array(v) print("\nCalculating gradient...") - dI = numpy.gradient(I) + dI = np.gradient(I.astype(dtype), v[0], v[1], v[2]) #calculate the dimensions of the tensor field field_dims = [dI[0].shape[0], dI[0].shape[1], dI[0].shape[2], 6] #allocate space for the tensor field - Tg = numpy.zeros(field_dims, dtype=numpy.float32) + Tg = np.zeros(field_dims, dtype=np.float32) #calculate the gradient components of the tensor ti = 0 @@ -69,13 +67,15 @@ def st3(I, s=1): bar.update(ti) #blur the tensor field - T = numpy.zeros(field_dims, dtype=numpy.float32) + T = np.zeros(field_dims, dtype=np.float32) print("\nConvolving tensor field...") bar = progressbar.ProgressBar() bar.max_value = 6 + sigma = s / v + print(sigma) for i in range(6): - T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], s) + T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], sigma) bar.update(i+1) return T @@ -107,13 +107,13 @@ def sym2mat(T): n = in_s[T.ndim - 1] #calculate the dimension of the symmetric matrix - d = int(0.5 * (numpy.sqrt(8. * n + 1.) - 1.)) + d = int(0.5 * (np.sqrt(8. * n + 1.) - 1.)) #calculate the dimensions of the output field out_s = list(in_s)[:-1] + [d] + [d] #allocate space for the output field - R = numpy.zeros(out_s) + R = np.zeros(out_s) ni = 0 for i in range(d): @@ -125,9 +125,8 @@ def sym2mat(T): return R -def st2vec(S, vector='largest'): - #Create a color image from a 2D or 3D structure tensor slice - +def st2vec(S, vector=0): + if(S.ndim != 3): print("ERROR: a 2D slice is expected") return @@ -137,27 +136,23 @@ def st2vec(S, vector='largest'): del(S) #calculate the eigenvectors and eigenvalues - l, v = numpy.linalg.eig(T) + l, v = np.linalg.eig(T) #get the dimension of the tensor field d = T.shape[2] #allocate space for the vector field - V = numpy.zeros([T.shape[0], T.shape[1], 3]) - + V = np.zeros([T.shape[0], T.shape[1], 3]) + + #arrange the indices for each pixel from smallest to largest eigenvalue idx = l.argsort() for di in range(d): - if vector == 'smallest': - b = idx[:, :, 0] == di - elif vector == 'largest': - b = idx[:, :, d-1] == di - else: - b = idx[:, :, 1] == di + b = idx[:, :, -1-vector] == di V[b, 0:d] = v[b, :, di] return V - + def loadstack(filemask): #Load an image stack as a 3D grayscale array @@ -171,7 +166,7 @@ def loadstack(filemask): Z = len(files) #allocate space for the image stack - M = numpy.zeros([X, Y, Z]).astype('float32') + M = np.zeros([X, Y, Z]).astype('float32') #create a progress bar bar = progressbar.ProgressBar() @@ -184,15 +179,16 @@ def loadstack(filemask): bar.update(z+1) return M -def anisotropy(S): +#calculate the anisotropy of a structure tensor given the tensor field S +def anisotropy3(S): Sf = sym2mat(S) #calculate the eigenvectors and eigenvalues - l, v = numpy.linalg.eig(Sf) + l, v = np.linalg.eig(Sf) #store the sorted eigenvalues - ls = numpy.sort(l) + ls = np.sort(l) l0 = ls[:, :, 0] l1 = ls[:, :, 1] l2 = ls[:, :, 2] @@ -210,10 +206,38 @@ def anisotropy(S): l_hat = (l0 + l1 + l2)/3 fa_num = (l2 - l_hat) ** 2 + (l1 - l_hat) ** 2 + (l0 - l_hat) ** 2; fa_den = l0 ** 2 + l1 ** 2 + l2 ** 2 - FA = numpy.sqrt(3./2.) * numpy.sqrt(fa_num) / numpy.sqrt(fa_den) + FA = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den) return FA, Cl, Cp, Cs +#calculate the fractional anisotropy +def fa(S): + Sf = sym2mat(S) + + #calculate the eigenvectors and eigenvalues + l, v = np.linalg.eig(Sf) + + #store the sorted eigenvalues + ls = np.sort(l) + l0 = ls[:, :, 0] + l1 = ls[:, :, 1] + + #if this is a 2D tensor, calculate and return the coherence + if(S.shape[2] == 3): + C = ((l0 - l1) / (l0 + l1)) ** 2 + return C + + #if this is a 3D tensor + elif(S.shape[2] == 6): + l2 = ls[:, :, 2] + + #calculate the fractional anisotropy + l_hat = (l0 + l1 + l2)/3 + fa_num = (l2 - l_hat) ** 2 + (l1 - l_hat) ** 2 + (l0 - l_hat) ** 2; + fa_den = l0 ** 2 + l1 ** 2 + l2 ** 2 + FA = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den) + return FA + def st2amira(filename, T): #generates a tensor field that can be imported into Amira @@ -225,7 +249,7 @@ def st2amira(filename, T): # 5 dz dz ----> 5 #swap the 2nd and 3rd tensor components - A = numpy.copy(T) + A = np.copy(T) A[..., 3] = T[..., 2] A[..., 2] = T[..., 3] @@ -246,7 +270,7 @@ def resample3(T, s=2): s = s * 3 elif len(s) == 2: s.insert(1, s[0]) - s = numpy.array(s) + s = np.array(s) bar = progressbar.ProgressBar() bar.max_value = T.shape[3] @@ -270,7 +294,7 @@ def color3(prefix, T, vector='largest', aniso=True): for z in range(T.shape[2]): S = T[:, :, z, :] #get the slice V = st2vec(S, vector='smallest') #calculate the vector - C = numpy.absolute(V) #calculate the absolute value + C = np.absolute(V) #calculate the absolute value if aniso == True: #if the image is scaled by anisotropy FA, Cl, Cp, Cs = anisotropy(S) #calculate the anisotropy of the slice @@ -279,9 +303,121 @@ def color3(prefix, T, vector='largest', aniso=True): elif vector == 'smallest': A = Cp else: #otherwise just scale by 1 - A = numpy.ones(T.shape[0], T.shape[1]) - image = C * numpy.expand_dims(A, 3) + A = np.ones(T.shape[0], T.shape[1]) + image = C * np.expand_dims(A, 3) filename = prefix + str(z).zfill(4) + ".bmp" scipy.misc.imsave(filename, image) - bar.update(z + 1) \ No newline at end of file + bar.update(z + 1) + +def st2stack(T, outfile, **kwargs): + eigenvector = False #initialize the colormap flags to false + aniso_color = False + aniso_alpha = False + #image = False + aniso_pwr = 1 + cimage_pwr = 1 + aimage_pwr = 1 + anisostretch = 1 #set the contrast stretch parameter + alpha_channel = False + alpha_image = False + color_image = False + + for k,v in kwargs.items(): #for each argument + if(k == "ev"): #if the user wants a colormap based on an eigenvector + eigenvector = True #set the eigenvector flag to true + ev = v #save the desired eigenvector + if(k == "aniso"): #if the user wants to factor in anisotropy + aniso = True #set the anisotropy flag to true + aniso_channel = v #save the anisotropy channel + if(k == "aniso_color"): + aniso_color = v + if(k == "aniso_alpha"): + aniso_alpha = v + if(k == "apwr"): #if the user wants to amplify the anisotropy + aniso_pwr = v #set the anisotropy exponent + if(k == "cipwr"): #if the user specifies the image power + cimage_pwr = v + if(k == "aipwr"): + aimage_pwr = v + if(k == "alphaimage"): + Ia = v + alpha_image = True + if(k == "colorimage"): + Ic = v + color_image = True + if(k == "anisostretch"): + anisostretch = v + if(k == "alpha"): + alpha_channel = v + + bar = progressbar.ProgressBar() + bar.max_value = T.shape[2] + for i in range(0, T.shape[2]): + #for i in range(0, 50): + + if(alpha_image or alpha_channel): + img = np.ones([T.shape[0], T.shape[1], 4]) + else: + img = np.ones([T.shape[0], T.shape[1], 3]) + if(eigenvector): + V = st2vec(T[:, :, i], ev) #get the vector field for slice i corresponding to eigenvector ev + img[:, :, 0:3] = V #update the image with the vector field information + if(aniso): #if the user is requesting anisotropy be incorporated into the image + FA, Cl, Cp, Cs = anisotropy(T[:, :, i]) #calculate the anisotropy of the tensors in slice i + if(aniso_channel == "fa"): + A = FA + elif(aniso_channel == "l"): + A = Cl + elif(aniso_channel == "p"): + A = Cp + else: + A = Cs + if(aniso_alpha): + print("rendering anisotropy to the alpha channel") + img[:, :, 3] = A ** aniso_pwr * anisostretch + if(aniso_color): + print("rendering anisotropy to the color channel") + img[:, :, 0:3] = img[:, :, 0:3] * np.expand_dims(A ** aniso_pwr, 3) * anisostretch + if(alpha_image): + img[:, :, 3] = Ia[:, :, i]/255 ** aimage_pwr + if(color_image): + img[:, :, 0:3] = img[:, :, 0:3] * (np.expand_dims(Ic[:, :, i], 3)/255) ** cimage_pwr #multiply the vector field by the image intensity + #outname = outfile + str(i).zfill(3) + ".bmp" #get the file name to be saved + outname = outfile.replace("*", str(i).zfill(3)) + + sp.misc.imsave(outname, np.ndarray.astype(np.abs(img)*255, "uint8")) #save the output image + bar.update(i+1) + + +#this function takes a 3D image and turns it into a stack of color images based on the structure tensor +def img2stack(I, outfile, **kwargs): + + vs = [1, 1, 1] #set the default voxel size to 1 + w = 5 + + for k,v in kwargs.items(): #for each argument + if(k == "voxelsize"): #if the voxel size is specified + if(len(v) == 1): #if the user just specifies one value + vs = [v] * 3 #assume that the voxels are isotropic, create a list of 3 v's + elif(len(v) == 2): #if the user specifies two values + vs[0] = v[0] #assume that the voxels are isotropic along (x, y) and anisotropic along z + vs[1] = v[0] + vs[2] = v[1] + elif(len(v) == 3): + vs = v + if(k == "window"): #if the user specifies a window size + w = v + + T = st3(I, w, vs) #calculate the structure tensor + + st2stack(T, outfile, **kwargs) + +def stack2stack(infile_mask, outfile, **kwargs): + + I = loadstack(infile_mask) #load the file mask + for k,v in kwargs.items(): #for each argument + if(k == "ipwr"): + img = I + + img2stack(I, outfile, image=img, **kwargs) #call the function to convert the image to an output ST stack \ No newline at end of file diff --git a/stim/image/image.h b/stim/image/image.h index 1307b48..cd88de2 100644 --- a/stim/image/image.h +++ b/stim/image/image.h @@ -157,6 +157,24 @@ public: } #endif + //determines if a filename represents a valid file format that can be loaded/saved + static bool test_filename(std::string f) { + stim::filename fname = f; + std::string ext = fname.extension(); +#ifdef USING_OPENCV + if (ext == "bmp" || + ext == "jpg" || + ext == "png" || + ext == "pbm" || + ext == "tif" ) + return true; +#else + if (ext == "pbm" || ext == "bmp") + return true; +#endif + return false; + } + //save a Netpbm file void load_netpbm(std::string filename) { std::ifstream infile(filename.c_str(), std::ios::in | std::ios::binary); //open an output file @@ -423,6 +441,11 @@ public: return img[idx(x, y, c)]; } + /// This function returns a pixel reference based on a 1D index into the image + T& operator()(size_t i) { + return img[i]; + } + /// Set all elements in the image to a given scalar value /// @param v is the value used to set all values in the image diff --git a/stim/optics/scalarfield.h b/stim/optics/scalarfield.h index df049d0..6aebc8c 100644 --- a/stim/optics/scalarfield.h +++ b/stim/optics/scalarfield.h @@ -428,7 +428,8 @@ public: } } - /// Propagate the field along its orthogonal direction by a distance d + /// Apply a low pass filter preserving all frequencies lower than or equal to "highest" + // @param highest is the highest frequency passed void lowpass(T highest){ cpu_scalar_lowpass(E, E, X.len(), Y.len(), highest, R[0], R[1]); } diff --git a/stim/visualization/camera.h b/stim/visualization/camera.h index c55351c..6eb3c5e 100644 --- a/stim/visualization/camera.h +++ b/stim/visualization/camera.h @@ -3,6 +3,7 @@ #include #include +#include #ifndef STIM_CAMERA_H #define STIM_CAMERA_H @@ -71,6 +72,7 @@ public: delta = focus; focus -= delta; + std::cout << focus << std::endl; Dolly(d*delta); } diff --git a/stim/visualization/glut_template.cpp b/stim/visualization/glut_template.cpp index 40be42e..784bc18 100644 --- a/stim/visualization/glut_template.cpp +++ b/stim/visualization/glut_template.cpp @@ -11,6 +11,72 @@ #include #include +void render_scene(); + +void draw_colorcube() { + glBegin(GL_LINES); //start a line series + glColor3f(0, 0, 0); + glVertex3f(-0.5, -0.5, -0.5); + glColor3f(1, 0, 0); + glVertex3f(0.5, -0.5, -0.5); + + glColor3f(0, 0, 0); + glVertex3f(-0.5, -0.5, -0.5); + glColor3f(0, 1, 0); + glVertex3f(-0.5, 0.5, -0.5); + + glColor3f(0, 0, 0); + glVertex3f(-0.5, -0.5, -0.5); + glColor3f(0, 0, 1); + glVertex3f(-0.5, -0.5, 0.5); + + glColor3f(1, 1, 1); + glVertex3f(0.5, 0.5, 0.5); + glColor3f(1, 1, 0); + glVertex3f(0.5, 0.5, -0.5); + + glColor3f(1, 1, 1); + glVertex3f(0.5, 0.5, 0.5); + glColor3f(1, 0, 1); + glVertex3f(0.5, -0.5, 0.5); + + glColor3f(1, 1, 1); + glVertex3f(0.5, 0.5, 0.5); + glColor3f(0, 1, 1); + glVertex3f(-0.5, 0.5, 0.5); + + glColor3f(1, 0, 0); + glVertex3f(0.5, -0.5, -0.5); + glColor3f(1, 1, 0); + glVertex3f(0.5, 0.5, -0.5); + + glColor3f(1, 0, 0); + glVertex3f(0.5, -0.5, -0.5); + glColor3f(1, 0, 1); + glVertex3f(0.5, -0.5, 0.5); + + glColor3f(0, 1, 0); + glVertex3f(-0.5, 0.5, -0.5); + glColor3f(1, 1, 0); + glVertex3f(0.5, 0.5, -0.5); + + glColor3f(0, 1, 0); + glVertex3f(-0.5, 0.5, -0.5); + glColor3f(0, 1, 1); + glVertex3f(-0.5, 0.5, 0.5); + + glColor3f(0, 0, 1); + glVertex3f(-0.5, -0.5, 0.5); + glColor3f(1, 0, 1); + glVertex3f(0.5, -0.5, 0.5); + + glColor3f(0, 0, 1); + glVertex3f(-0.5, -0.5, 0.5); + glColor3f(0, 1, 1); + glVertex3f(-0.5, 0.5, 0.5); + glEnd(); +}} + struct glut_template_struct { float theta_scale = 0.01f; float phi_scale = 0.01f; @@ -53,10 +119,11 @@ void glut_display() { glLoadIdentity(); //set it to the identity matrix stim::vec3 p = gt.cam.getPosition(); //get the camera parameters stim::vec3 u = gt.cam.getUp(); - stim::vec3 d = gt.cam.getDirection(); - 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 + stim::vec3 c = gt.cam.getLookAt(); + 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 render_axes(); //render the axes if the user requests them + draw_colorcube(); glutSwapBuffers(); //swap in the back buffer (double-buffering is used to prevent tearing) } -- libgit2 0.21.4