Commit 685a889c859634e1f0cf5a4b75ccdb3b82e0d13c

Authored by Jiaming Guo
2 parents 5068402b 56737905

Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

Showing 34 changed files with 1034 additions and 339 deletions   Show diff stats
matlab/cls_ConfusionMatrix.m 0 → 100644
  1 +function M = cls_ConfusionMatrix(GT, T)
  2 +
  3 +%calculate the classes (unique elements in the GT array)
  4 +C = unique(GT);
  5 +nc = length(C); %calculate the number of classes
  6 +M = zeros(nc); %allocate space for the confusion matrix
  7 +
  8 +%for each class
  9 +for ci = 1:nc
  10 + for cj = 1:nc
  11 + M(ci, cj) = nnz((GT == C(ci)) .* (T == C(cj)));
  12 + end
  13 +end
0 \ No newline at end of file 14 \ No newline at end of file
matlab/cls_MeanClassFeatures.m 0 → 100644
  1 +function S = cls_MeanClassFeatures(F, T)
  2 +%Calculates the mean set of features for each class given the feature matrix F and targets T
  3 +
  4 +C = unique(T); %get the class IDs
  5 +nc = length(C);
  6 +
  7 +S = zeros(nc, size(F, 2)); %allocate space for the mean feature vectors
  8 +for c = 1:nc %for each class
  9 + S(c, :) = mean(F(T == C(c), :)); %calculate the mean feature vector for class c
  10 +end
  11 +
  12 +S = S';
0 \ No newline at end of file 13 \ No newline at end of file
matlab/cls_PlotConfusionMatrix.m 0 → 100644
  1 +function cls_PlotConfusionMatrix(M)
  2 +
  3 +
  4 +%normalize each row by its column
  5 +sum_cols = repmat(sum(M, 1), size(M, 1), 1);
  6 +Mc = M ./ sum_cols;
  7 +subplot(2, 1, 1),
  8 +bar(Mc');
  9 +
  10 +sum_rows = repmat(sum(M, 2), 1, size(M, 2));
  11 +Mr = M ./ sum_rows;
  12 +subplot(2, 1, 2),
  13 +bar(Mr);
0 \ No newline at end of file 14 \ No newline at end of file
matlab/stim_images2matrix.m 0 → 100644
  1 +function S = stim_images2matrix(filemask)
  2 +%This function loads a set of images as a 3D matrix. Color images are
  3 +%converted to grayscale when loaded, so the resulting matrix is always 3D
  4 +%with size X x Y x Z, where:
  5 +% X is the size of the images along the X axis
  6 +% Y is the size of the images along the Y axis
  7 +% Z is the number of images
  8 +%
  9 +% all images are assumed to be the same size (though they do not have to
  10 +% be the same file format or number of bits per pixel
  11 +
  12 + files = dir(filemask);
  13 +
  14 + %figure out the file size
  15 + I = imread([files(1).folder '/' files(1).name]);
  16 + X = size(I, 1);
  17 + Y = size(I, 2);
  18 + Z = length(files);
  19 +
  20 + S = zeros(X, Y, Z, 'uint8');
  21 +
  22 + h = waitbar(0, ['Loading ' num2str(Z) ' images...']);
  23 + for i = 1:Z
  24 + I = rgb2gray(imread([files(1).folder '/' files(1).name]));
  25 + S(:, :, i) = I;
  26 + waitbar(i/Z, h);
  27 + end
  28 + close(h);
  29 +end
  30 +
  31 +
python/structen.py 0 → 100644
  1 +# -*- coding: utf-8 -*-
  2 +"""
  3 +Created on Sun Mar 12 21:54:40 2017
  4 +
  5 +@author: david
  6 +"""
  7 +
  8 +import numpy
  9 +import scipy.ndimage
  10 +import progressbar
  11 +import glob
  12 +
  13 +def st2(I, s=1, dtype=numpy.float32):
  14 +
  15 + #calculates the 2D structure tensor for an image using a gaussian window with standard deviation s
  16 +
  17 + #calculate the gradient
  18 + dI = numpy.gradient(I)
  19 +
  20 + #calculate the dimensions of the tensor field
  21 + field_dims = [dI[0].shape[0], dI[0].shape[1], 3]
  22 +
  23 + #allocate space for the tensor field
  24 + Tg = numpy.zeros(field_dims, dtype=dtype)
  25 +
  26 + #calculate the gradient components of the tensor
  27 + ti = 0
  28 + for i in range(2):
  29 + for j in range(i + 1):
  30 + Tg[:, :, ti] = dI[j] * dI[i]
  31 + ti = ti + 1
  32 +
  33 + #blur the tensor field
  34 + T = numpy.zeros(field_dims, dtype=dtype)
  35 +
  36 + for i in range(3):
  37 + T[:, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, i], [s, s])
  38 +
  39 +
  40 + return T
  41 +
  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
  44 + #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 +
  52 + print("\nCalculating gradient...")
  53 + dI = numpy.gradient(I)
  54 + #calculate the dimensions of the tensor field
  55 + field_dims = [dI[0].shape[0], dI[0].shape[1], dI[0].shape[2], 6]
  56 +
  57 + #allocate space for the tensor field
  58 + Tg = numpy.zeros(field_dims, dtype=numpy.float32)
  59 +
  60 + #calculate the gradient components of the tensor
  61 + ti = 0
  62 + print("Calculating tensor components...")
  63 + bar = progressbar.ProgressBar()
  64 + bar.max_value = 6
  65 + for i in range(3):
  66 + for j in range(i + 1):
  67 + Tg[:, :, :, ti] = dI[j] * dI[i]
  68 + ti = ti + 1
  69 + bar.update(ti)
  70 +
  71 + #blur the tensor field
  72 + T = numpy.zeros(field_dims, dtype=numpy.float32)
  73 +
  74 + print("\nConvolving tensor field...")
  75 + bar = progressbar.ProgressBar()
  76 + bar.max_value = 6
  77 + for i in range(6):
  78 + T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], s)
  79 + bar.update(i+1)
  80 +
  81 + return T
  82 +
  83 +def st(I, s=1):
  84 + if I.ndim == 3:
  85 + return st3(I, s)
  86 + elif I.ndim == 2:
  87 + return st2(I, s)
  88 + else:
  89 + print("Image must be 2D or 3D")
  90 + return
  91 +
  92 +
  93 +
  94 +def sym2mat(T):
  95 + #Calculate the full symmetric matrix from a list of lower triangular elements.
  96 + #The lower triangular components in the input field T are assumed to be the
  97 + # final dimension of the input matrix.
  98 +
  99 + # | 1 2 4 7 |
  100 + # | 0 3 5 8 |
  101 + # | 0 0 6 9 |
  102 + # | 0 0 0 10 |
  103 +
  104 + in_s = T.shape
  105 +
  106 + #get the number of tensor elements
  107 + n = in_s[T.ndim - 1]
  108 +
  109 + #calculate the dimension of the symmetric matrix
  110 + d = int(0.5 * (numpy.sqrt(8. * n + 1.) - 1.))
  111 +
  112 + #calculate the dimensions of the output field
  113 + out_s = list(in_s)[:-1] + [d] + [d]
  114 +
  115 + #allocate space for the output field
  116 + R = numpy.zeros(out_s)
  117 +
  118 + ni = 0
  119 + for i in range(d):
  120 + for j in range(i + 1):
  121 + R[..., i, j] = T[..., ni]
  122 + if i != j:
  123 + R[..., j, i] = T[..., ni]
  124 + ni = ni + 1
  125 +
  126 + return R
  127 +
  128 +def st2vec(S, vector='largest'):
  129 + #Create a color image from a 2D or 3D structure tensor slice
  130 +
  131 + #convert the field to a full rank-2 tensor
  132 + T = sym2mat(S);
  133 + del(S)
  134 +
  135 + #calculate the eigenvectors and eigenvalues
  136 + l, v = numpy.linalg.eig(T)
  137 +
  138 + #get the dimension of the tensor field
  139 + d = T.shape[2]
  140 +
  141 + #allocate space for the vector field
  142 + V = numpy.zeros([T.shape[0], T.shape[1], 3])
  143 +
  144 + idx = l.argsort()
  145 +
  146 + for di in range(d):
  147 + if vector == 'smallest':
  148 + b = idx[:, :, 0] == di
  149 + elif vector == 'largest':
  150 + b = idx[:, :, d-1] == di
  151 + else:
  152 + b = idx[:, :, 1] == di
  153 + V[b, 0:d] = v[b, :, di]
  154 +
  155 + return V
  156 +
  157 +def loadstack(filemask):
  158 + #Load an image stack as a 3D grayscale array
  159 +
  160 + #get a list of all files matching the given mask
  161 + files = [file for file in glob.glob(filemask)]
  162 +
  163 + #calculate the size of the output stack
  164 + I = scipy.misc.imread(files[0])
  165 + X = I.shape[0]
  166 + Y = I.shape[1]
  167 + Z = len(files)
  168 +
  169 + #allocate space for the image stack
  170 + M = numpy.zeros([X, Y, Z]).astype('float32')
  171 +
  172 + #create a progress bar
  173 + bar = progressbar.ProgressBar()
  174 + bar.max_value = Z
  175 +
  176 + #for each file
  177 + for z in range(Z):
  178 + #load the file and save it to the image stack
  179 + M[:, :, z] = scipy.misc.imread(files[z], flatten="True").astype('float32')
  180 + bar.update(z+1)
  181 + return M
  182 +
  183 +def anisotropy(S):
  184 +
  185 + Sf = sym2mat(S)
  186 +
  187 + #calculate the eigenvectors and eigenvalues
  188 + l, v = numpy.linalg.eig(Sf)
  189 +
  190 + #store the sorted eigenvalues
  191 + ls = numpy.sort(l)
  192 + l0 = ls[:, :, 0]
  193 + l1 = ls[:, :, 1]
  194 + l2 = ls[:, :, 2]
  195 +
  196 + #calculate the linear anisotropy
  197 + Cl = (l2 - l1)/(l2 + l1 + l0)
  198 +
  199 + #calculate the planar anisotropy
  200 + Cp = 2 * (l1 - l0) / (l2 + l1 + l0)
  201 +
  202 + #calculate the spherical anisotropy
  203 + Cs = 3 * l0 / (l2 + l1 + l0)
  204 +
  205 + #calculate the fractional anisotropy
  206 + l_hat = (l0 + l1 + l2)/3
  207 + fa_num = (l2 - l_hat) ** 2 + (l1 - l_hat) ** 2 + (l0 - l_hat) ** 2;
  208 + fa_den = l0 ** 2 + l1 ** 2 + l2 ** 2
  209 + FA = numpy.sqrt(3./2.) * numpy.sqrt(fa_num) / numpy.sqrt(fa_den)
  210 +
  211 + return FA, Cl, Cp, Cs
  212 +
  213 +def st2amira(filename, T):
  214 + #generates a tensor field that can be imported into Amira
  215 +
  216 + # 0 dx dx ----> 0
  217 + # 1 dx dy ----> 1
  218 + # 2 dy dy ----> 3
  219 + # 3 dx dz ----> 2
  220 + # 4 dy dz ----> 4
  221 + # 5 dz dz ----> 5
  222 +
  223 + #swap the 2nd and 3rd tensor components
  224 + A = numpy.copy(T)
  225 + A[..., 3] = T[..., 2]
  226 + A[..., 2] = T[..., 3]
  227 +
  228 + #roll the tensor axis so that it is the leading component
  229 + #A = numpy.rollaxis(A, A.ndim - 1)
  230 + A.tofile(filename)
  231 + print("\n", A.shape)
  232 +
  233 +def resample3(T, s=2):
  234 + #resample a tensor field by an integer factor s
  235 + #This function first convolves the field with a box filter and then
  236 + # re-samples to create a smaller field
  237 +
  238 + #check the format for the window size
  239 + if type(s) is not list:
  240 + s = [s] * 3
  241 + elif len(s) == 1:
  242 + s = s * 3
  243 + elif len(s) == 2:
  244 + s.insert(1, s[0])
  245 + s = numpy.array(s)
  246 +
  247 + bar = progressbar.ProgressBar()
  248 + bar.max_value = T.shape[3]
  249 +
  250 + #blur with a uniform box filter of size r
  251 + for t in range(T.shape[3]):
  252 + T[..., t] = scipy.ndimage.filters.uniform_filter(T[..., t], 2 * s)
  253 + bar.update(t+1)
  254 +
  255 + #resample at a rate of r
  256 + R = T[::s[0], ::s[1], ::s[2], :]
  257 + return R
  258 +
  259 +def color3(prefix, T, vector='largest', aniso=True):
  260 + #Saves a stack of color images corresponding to the eigenvector and optionally scaled by anisotropy
  261 +
  262 + bar = progressbar.ProgressBar()
  263 + bar.max_value = T.shape[2]
  264 +
  265 + #for each z-axis slice
  266 + for z in range(T.shape[2]):
  267 + S = T[:, :, z, :] #get the slice
  268 + V = st2vec(S, vector='smallest') #calculate the vector
  269 + C = numpy.absolute(V) #calculate the absolute value
  270 +
  271 + if aniso == True: #if the image is scaled by anisotropy
  272 + FA, Cl, Cp, Cs = anisotropy(S) #calculate the anisotropy of the slice
  273 + if vector == 'largest':
  274 + A = Cl
  275 + elif vector == 'smallest':
  276 + A = Cp
  277 + else: #otherwise just scale by 1
  278 + A = numpy.ones(T.shape[0], T.shape[1])
  279 + image = C * numpy.expand_dims(A, 3)
  280 +
  281 + filename = prefix + str(z).zfill(4) + ".bmp"
  282 + scipy.misc.imsave(filename, image)
  283 + bar.update(z + 1)
0 \ No newline at end of file 284 \ No newline at end of file
stim/cuda/ivote.cuh deleted
1 -#ifndef STIM_CUDA_IVOTE_H  
2 -#define STIM_CUDA_IVOTE_H  
3 -  
4 -#include <stim/cuda/ivote/down_sample.cuh>  
5 -#include <stim/cuda/ivote/local_max.cuh>  
6 -#include <stim/cuda/ivote/update_dir.cuh>  
7 -#include <stim/cuda/ivote/vote.cuh>  
8 -  
9 -namespace stim{  
10 - namespace cuda{  
11 -  
12 - }  
13 -}  
14 -  
15 -  
16 -  
17 -#endif  
18 \ No newline at end of file 0 \ No newline at end of file
stim/cuda/ivote/down_sample.cuh deleted
1 -#ifndef STIM_CUDA_DOWN_SAMPLE_H  
2 -#define STIM_CUDA_DOWN_SAMPLE_H  
3 -  
4 -#include <iostream>  
5 -#include <cuda.h>  
6 -#include <stim/cuda/cudatools.h>  
7 -#include <stim/cuda/templates/gaussian_blur.cuh>  
8 -  
9 -namespace stim{  
10 - namespace cuda{  
11 -  
12 - template<typename T>  
13 - __global__ void down_sample(T* gpuI, T* gpuI0, T resize, unsigned int x, unsigned int y){  
14 -  
15 - unsigned int sigma_ds = 1/resize;  
16 - unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));  
17 - unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));  
18 -  
19 -  
20 - // calculate the 2D coordinates for this current thread.  
21 - int xi = blockIdx.x * blockDim.x + threadIdx.x;  
22 - int yi = blockIdx.y;  
23 - // convert 2D coordinates to 1D  
24 - int i = yi * x_ds + xi;  
25 -  
26 - if(xi< x_ds && yi< y_ds){  
27 -  
28 - int x_org = xi * sigma_ds ;  
29 - int y_org = yi * sigma_ds ;  
30 - int i_org = y_org * x + x_org;  
31 - gpuI[i] = gpuI0[i_org];  
32 - }  
33 -  
34 - }  
35 -  
36 -  
37 - /// Applies a Gaussian blur to a 2D image stored on the GPU  
38 - template<typename T>  
39 - void gpu_down_sample(T* gpuI, T* gpuI0, T resize, size_t x, size_t y){  
40 -  
41 -  
42 - unsigned int sigma_ds = (unsigned int)(1.0f/resize);  
43 - size_t x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));  
44 - size_t y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));  
45 -  
46 - //get the number of pixels in the image  
47 -// unsigned int pixels_ds = x_ds * y_ds;  
48 -  
49 - unsigned int max_threads = stim::maxThreadsPerBlock();  
50 - dim3 threads(max_threads, 1);  
51 - dim3 blocks(x_ds/threads.x + (x_ds %threads.x == 0 ? 0:1) , y_ds);  
52 -  
53 - stim::cuda::gpu_gaussian_blur2<float>(gpuI0, sigma_ds,x ,y);  
54 -  
55 - //resample the image  
56 - down_sample<float> <<< blocks, threads >>>(gpuI, gpuI0, resize, x, y);  
57 -  
58 - }  
59 -  
60 - /// Applies a Gaussian blur to a 2D image stored on the CPU  
61 - template<typename T>  
62 - void cpu_down_sample(T* re_img, T* image, T resize, unsigned int x, unsigned int y){  
63 -  
64 - //get the number of pixels in the image  
65 - unsigned int pixels = x * y;  
66 - unsigned int bytes = sizeof(T) * pixels;  
67 -  
68 - unsigned int sigma_ds = 1/resize;  
69 - unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));  
70 - unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));  
71 - unsigned int bytes_ds = sizeof(T) * x_ds * y_ds;  
72 -  
73 -  
74 -  
75 - //allocate space on the GPU for the original image  
76 - T* gpuI0;  
77 - cudaMalloc(&gpuI0, bytes);  
78 -  
79 -  
80 - //copy the image data to the GPU  
81 - cudaMemcpy(gpuI0, image, bytes, cudaMemcpyHostToDevice);  
82 -  
83 - //allocate space on the GPU for the down sampled image  
84 - T* gpuI;  
85 - cudaMalloc(&gpuI, bytes_ds);  
86 -  
87 - //run the GPU-based version of the algorithm  
88 - gpu_down_sample<T>(gpuI, gpuI0, resize, x, y);  
89 -  
90 - //copy the image data to the GPU  
91 - cudaMemcpy(re_img, gpuI, bytes_ds, cudaMemcpyHostToDevice);  
92 -  
93 - cudaFree(gpuI0);  
94 - cudeFree(gpuI);  
95 - }  
96 -  
97 - }  
98 -}  
99 -  
100 -#endif  
stim/cuda/ivote/re_sample.cuh deleted
1 -#ifndef STIM_CUDA_RE_SAMPLE_H  
2 -#define STIM_CUDA_RE_SAMPLE_H  
3 -  
4 -#include <iostream>  
5 -#include <cuda.h>  
6 -#include <stim/cuda/cudatools.h>  
7 -#include <stim/cuda/templates/gaussian_blur.cuh>  
8 -  
9 -namespace stim{  
10 - namespace cuda{  
11 -  
12 - template<typename T>  
13 - __global__ void cuda_re_sample(T* gpuI, T* gpuI0, T resize, unsigned int x, unsigned int y){  
14 -  
15 - unsigned int sigma_ds = 1/resize;  
16 - unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));  
17 - unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));  
18 -  
19 -  
20 - // calculate the 2D coordinates for this current thread.  
21 - int xi = blockIdx.x * blockDim.x + threadIdx.x;  
22 - int yi = blockIdx.y;  
23 - // convert 2D coordinates to 1D  
24 - int i = yi * x + xi;  
25 -  
26 - if(xi< x && yi< y){  
27 - if(xi%sigma_ds==0){  
28 - if(yi%sigma_ds==0){  
29 - gpuI[i] = gpuI0[(yi/sigma_ds)*x_ds + xi/sigma_ds];  
30 - }  
31 - }  
32 - else gpuI[i] = 0;  
33 -  
34 - //int x_org = xi * sigma_ds ;  
35 - //int y_org = yi * sigma_ds ;  
36 - //int i_org = y_org * x + x_org;  
37 - //gpuI[i] = gpuI0[i_org];  
38 - }  
39 -  
40 - }  
41 -  
42 -  
43 - /// Applies a Gaussian blur to a 2D image stored on the GPU  
44 - template<typename T>  
45 - void gpu_re_sample(T* gpuI, T* gpuI0, T resize, unsigned int x, unsigned int y){  
46 -  
47 -  
48 - //unsigned int sigma_ds = 1/resize;  
49 - //unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));  
50 - //unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));  
51 -  
52 - //get the number of pixels in the image  
53 - //unsigned int pixels_ds = x_ds * y_ds;  
54 -  
55 - unsigned int max_threads = stim::maxThreadsPerBlock();  
56 - dim3 threads(max_threads, 1);  
57 - dim3 blocks(x/threads.x + (x %threads.x == 0 ? 0:1) , y);  
58 -  
59 - //stim::cuda::gpu_gaussian_blur2<float>(gpuI0, sigma_ds,x ,y);  
60 -  
61 - //resample the image  
62 - cuda_re_sample<float> <<< blocks, threads >>>(gpuI, gpuI0, resize, x, y);  
63 -  
64 - }  
65 -  
66 - /// Applies a Gaussian blur to a 2D image stored on the CPU  
67 - template<typename T>  
68 - void cpu_re_sample(T* out, T* in, T resize, unsigned int x, unsigned int y){  
69 -  
70 - //get the number of pixels in the image  
71 - unsigned int pixels = x*y;  
72 - unsigned int bytes = sizeof(T) * pixels;  
73 -  
74 - unsigned int sigma_ds = 1/resize;  
75 - unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));  
76 - unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));  
77 - unsigned int bytes_ds = sizeof(T) * x_ds * y_ds;  
78 -  
79 -  
80 -  
81 - //allocate space on the GPU for the original image  
82 - T* gpuI0;  
83 - cudaMalloc(&gpuI0, bytes_ds);  
84 -  
85 -  
86 - //copy the image data to the GPU  
87 - cudaMemcpy(gpuI0, in, bytes_ds, cudaMemcpyHostToDevice);  
88 -  
89 - //allocate space on the GPU for the down sampled image  
90 - T* gpuI;  
91 - cudaMalloc(&gpuI, bytes);  
92 -  
93 - //run the GPU-based version of the algorithm  
94 - gpu_re_sample<T>(gpuI, gpuI0, resize, x, y);  
95 -  
96 - //copy the image data to the GPU  
97 - cudaMemcpy(re_img, gpuI, bytes_ds, cudaMemcpyHostToDevice);  
98 -  
99 - cudaFree(gpuI0);  
100 - cudeFree(gpuI);  
101 - }  
102 -  
103 - }  
104 -}  
105 -  
106 -#endif  
107 \ No newline at end of file 0 \ No newline at end of file
stim/cuda/ivote_atomic_bb.cuh deleted
1 -#ifndef STIM_CUDA_IVOTE_ATOMIC_BB_H  
2 -#define STIM_CUDA_IVOTE_ATOMIC_BB_H  
3 -  
4 -extern bool DEBUG;  
5 -#include <stim/cuda/ivote/down_sample.cuh>  
6 -#include <stim/cuda/ivote/local_max.cuh>  
7 -#include <stim/cuda/ivote/update_dir_bb.cuh>  
8 -#include <stim/cuda/ivote/vote_atomic_bb.cuh>  
9 -  
10 -namespace stim{  
11 - namespace cuda{  
12 -  
13 - }  
14 -}  
15 -  
16 -  
17 -  
18 -#endif  
19 \ No newline at end of file 0 \ No newline at end of file
stim/envi/agilent_binary.h
@@ -243,13 +243,13 @@ public: @@ -243,13 +243,13 @@ public:
243 if (device >= 0) { //if a CUDA device is specified 243 if (device >= 0) { //if a CUDA device is specified
244 int dev_count; 244 int dev_count;
245 HANDLE_ERROR(cudaGetDeviceCount(&dev_count)); //get the number of CUDA devices 245 HANDLE_ERROR(cudaGetDeviceCount(&dev_count)); //get the number of CUDA devices
246 - std::cout << "Number of CUDA devices: " << dev_count << std::endl; //output the number of CUDA devices 246 + //std::cout << "Number of CUDA devices: " << dev_count << std::endl; //output the number of CUDA devices
247 cudaDeviceProp prop; 247 cudaDeviceProp prop;
248 - std::cout << "CUDA devices----" << std::endl; 248 + //std::cout << "CUDA devices----" << std::endl;
249 for (int d = 0; d < dev_count; d++) { //for each CUDA device 249 for (int d = 0; d < dev_count; d++) { //for each CUDA device
250 cudaGetDeviceProperties(&prop, d); //get the property of the first device 250 cudaGetDeviceProperties(&prop, d); //get the property of the first device
251 //float cc = prop.major + prop.minor / 10.0f; //calculate the compute capability 251 //float cc = prop.major + prop.minor / 10.0f; //calculate the compute capability
252 - std::cout << d << ": [" << prop.major << "." << prop.minor << "] " << prop.name << std::endl; //display the device information 252 + //std::cout << d << ": [" << prop.major << "." << prop.minor << "] " << prop.name << std::endl; //display the device information
253 //if(cc > best_device_cc){ 253 //if(cc > best_device_cc){
254 // best_device_cc = cc; //if this is better than the previous device, use it 254 // best_device_cc = cc; //if this is better than the previous device, use it
255 // best_device_id = d; 255 // best_device_id = d;
@@ -258,7 +258,7 @@ public: @@ -258,7 +258,7 @@ public:
258 if (dev_count > 0 && dev_count > device) { //if the first device is not an emulator 258 if (dev_count > 0 && dev_count > device) { //if the first device is not an emulator
259 cudaGetDeviceProperties(&prop, device); //get the property of the requested CUDA device 259 cudaGetDeviceProperties(&prop, device); //get the property of the requested CUDA device
260 if (prop.major != 9999) { 260 if (prop.major != 9999) {
261 - std::cout << "Using device " << device << std::endl; 261 + //std::cout << "Using device " << device << std::endl;
262 HANDLE_ERROR(cudaSetDevice(device)); 262 HANDLE_ERROR(cudaSetDevice(device));
263 } 263 }
264 } 264 }
@@ -1368,6 +1368,39 @@ public: @@ -1368,6 +1368,39 @@ public:
1368 return false; 1368 return false;
1369 } 1369 }
1370 1370
  1371 + void band_bounds(double wavelength, size_t& low, size_t& high) {
  1372 + if (header.interleave == envi_header::BSQ) { //if the infile is bsq file
  1373 + if (header.data_type == envi_header::float32)
  1374 + ((bsq<float>*)file)->band_bounds(wavelength, low, high);
  1375 + else if (header.data_type == envi_header::float64)
  1376 + ((bsq<double>*)file)->band_bounds(wavelength, low, high);
  1377 + else {
  1378 + std::cout << "ERROR: unidentified data type" << std::endl;
  1379 + exit(1);
  1380 + }
  1381 + }
  1382 + else if (header.interleave == envi_header::BIL) {
  1383 + if (header.data_type == envi_header::float32)
  1384 + ((bil<float>*)file)->band_bounds(wavelength, low, high);
  1385 + else if (header.data_type == envi_header::float64)
  1386 + ((bil<double>*)file)->band_bounds(wavelength, low, high);
  1387 + else {
  1388 + std::cout << "ERROR: unidentified data type" << std::endl;
  1389 + exit(1);
  1390 + }
  1391 + }
  1392 + else if (header.interleave == envi_header::BIP) {
  1393 + if (header.data_type == envi_header::float32)
  1394 + ((bip<float>*)file)->band_bounds(wavelength, low, high);
  1395 + else if (header.data_type == envi_header::float64)
  1396 + ((bip<double>*)file)->band_bounds(wavelength, low, high);
  1397 + else {
  1398 + std::cout << "ERROR: unidentified data type" << std::endl;
  1399 + exit(1);
  1400 + }
  1401 + }
  1402 + }
  1403 +
1371 // Retrieve a spectrum at the specified 1D location 1404 // Retrieve a spectrum at the specified 1D location
1372 1405
1373 /// @param ptr is a pointer to pre-allocated memory of size B*sizeof(T) 1406 /// @param ptr is a pointer to pre-allocated memory of size B*sizeof(T)
@@ -62,31 +62,6 @@ protected: @@ -62,31 +62,6 @@ protected:
62 return (T)((1.0 - alpha) * low_v + alpha * high_v); //interpolate 62 return (T)((1.0 - alpha) * low_v + alpha * high_v); //interpolate
63 } 63 }
64 64
65 - /// Gets the two band indices surrounding a given wavelength  
66 - void band_bounds(double wavelength, size_t& low, size_t& high){  
67 - size_t B = Z();  
68 - for(high = 0; high < B; high++){  
69 - if(w[high] > wavelength) break;  
70 - }  
71 - low = 0;  
72 - if(high > 0)  
73 - low = high-1;  
74 - }  
75 -  
76 - /// Get the list of band numbers that bound a list of wavelengths  
77 - void band_bounds(std::vector<double> wavelengths,  
78 - std::vector<unsigned long long>& low_bands,  
79 - std::vector<unsigned long long>& high_bands){  
80 -  
81 - unsigned long long W = w.size(); //get the number of wavelengths in the list  
82 - low_bands.resize(W); //pre-allocate space for the band lists  
83 - high_bands.resize(W);  
84 -  
85 - for(unsigned long long wl = 0; wl < W; wl++){ //for each wavelength  
86 - band_bounds(wavelengths[wl], low_bands[wl], high_bands[wl]); //find the low and high bands  
87 - }  
88 - }  
89 -  
90 /// Returns the interpolated in the given spectrum based on the given wavelength 65 /// Returns the interpolated in the given spectrum based on the given wavelength
91 66
92 /// @param s is the spectrum in main memory of length Z() 67 /// @param s is the spectrum in main memory of length Z()
@@ -139,6 +114,31 @@ protected: @@ -139,6 +114,31 @@ protected:
139 } 114 }
140 115
141 public: 116 public:
  117 +
  118 + /// Gets the two band indices surrounding a given wavelength
  119 + void band_bounds(double wavelength, size_t& low, size_t& high) {
  120 + size_t B = Z();
  121 + for (high = 0; high < B; high++) {
  122 + if (w[high] > wavelength) break;
  123 + }
  124 + low = 0;
  125 + if (high > 0)
  126 + low = high - 1;
  127 + }
  128 +
  129 + /// Get the list of band numbers that bound a list of wavelengths
  130 + void band_bounds(std::vector<double> wavelengths,
  131 + std::vector<unsigned long long>& low_bands,
  132 + std::vector<unsigned long long>& high_bands) {
  133 +
  134 + unsigned long long W = w.size(); //get the number of wavelengths in the list
  135 + low_bands.resize(W); //pre-allocate space for the band lists
  136 + high_bands.resize(W);
  137 +
  138 + for (unsigned long long wl = 0; wl < W; wl++) { //for each wavelength
  139 + band_bounds(wavelengths[wl], low_bands[wl], high_bands[wl]); //find the low and high bands
  140 + }
  141 + }
142 /// Get a mask that has all pixels with inf or NaN values masked out (false) 142 /// Get a mask that has all pixels with inf or NaN values masked out (false)
143 void mask_finite(unsigned char* out_mask, unsigned char* mask, bool PROGRESS = false){ 143 void mask_finite(unsigned char* out_mask, unsigned char* mask, bool PROGRESS = false){
144 size_t XY = X() * Y(); 144 size_t XY = X() * Y();
stim/iVote/ivote2.cuh 0 → 100644
  1 +#ifndef STIM_IVOTE2_CUH
  2 +#define STIM_IVOTE2_CUH
  3 +
  4 +#include <iostream>
  5 +#include <fstream>
  6 +#include <stim/cuda/cudatools/error.h>
  7 +#include <stim/cuda/templates/gradient.cuh>
  8 +#include <stim/cuda/arraymath.cuh>
  9 +#include <stim/iVote/ivote2/iter_vote2.cuh>
  10 +#include <stim/iVote/ivote2/local_max.cuh>
  11 +#include <stim/math/constants.h>
  12 +#include <stim/math/vector.h>
  13 +#include <stim/visualization/colormap.h>
  14 +
  15 +
  16 +namespace stim {
  17 + // this function precomputes the atan2 values
  18 + template<typename T>
  19 + void atan_2(T* cpuTable, unsigned int rmax) {
  20 + int xsize = 2 * rmax + 1; //initialize the width and height of the window which atan2 are computed in.
  21 + int ysize = 2 * rmax + 1;
  22 + int yi = rmax; // assign the center coordinates of the atan2 window to yi and xi
  23 + int xi = rmax;
  24 + for (int xt = 0; xt < xsize; xt++) { //for each element in the atan2 table
  25 + for (int yt = 0; yt < ysize; yt++) {
  26 + int id = yt * xsize + xt; //convert the current 2D coordinates to 1D
  27 + int xd = xi - xt; // calculate the distance between the pixel and the center of the atan2 window
  28 + int yd = yi - yt;
  29 + T atan_2d = atan2((T)yd, (T)xd); // calculate the angle between the pixel and the center of the atan2 window and store the result.
  30 + cpuTable[id] = atan_2d;
  31 + }
  32 + }
  33 + }
  34 +
  35 + //this kernel invert the 2D image
  36 + template<typename T>
  37 + __global__ void cuda_invert(T* gpuI, size_t x, size_t y) {
  38 + // calculate the 2D coordinates for this current thread.
  39 + size_t xi = blockIdx.x * blockDim.x + threadIdx.x;
  40 + size_t yi = blockIdx.y * blockDim.y + threadIdx.y;
  41 +
  42 + if (xi >= x || yi >= y) return;
  43 + size_t i = yi * x + xi; // convert 2D coordinates to 1D
  44 + gpuI[i] = 255 - gpuI[i]; //invert the pixel intensity
  45 + }
  46 +
  47 +
  48 +
  49 + //this function calculate the threshold using OTSU method
  50 + template<typename T>
  51 + T th_otsu(T* pts, size_t pixels, unsigned int th_num = 20) {
  52 + T Imax = pts[0]; //initialize the maximum value to the first one
  53 + T Imin = pts[0]; //initialize the maximum value to the first on
  54 +
  55 + for (size_t n = 0; n < pixels; n++) { //for every value
  56 + if (pts[n] > Imax) { //if the value is higher than the current max
  57 + Imax = pts[n];
  58 + }
  59 + }
  60 + for (size_t n = 0; n< pixels; n++) { //for every value
  61 + if (pts[n] < Imin) { //if the value is higher than the current max
  62 + Imin = pts[n];
  63 + }
  64 + }
  65 +
  66 + T th_step = ((Imax - Imin) / th_num);
  67 + vector<T> var_b;
  68 + for (unsigned int t0 = 0; t0 < th_num; t0++) {
  69 + T th = t0 * th_step + Imin;
  70 + unsigned int n_b(0), n_o(0); //these variables save the number of elements that are below and over the threshold
  71 + T m_b(0), m_o(0); //these variables save the mean value for each cluster
  72 + for (unsigned int idx = 0; idx < pixels; idx++) {
  73 + if (pts[idx] <= th) {
  74 + m_b += pts[idx];
  75 + n_b += 1;
  76 + }
  77 + else {
  78 + m_o += pts[idx];
  79 + n_o += 1;
  80 + }
  81 + }
  82 +
  83 + m_b = m_b / n_b; //calculate the mean value for the below threshold cluster
  84 + m_o = m_o / n_o; //calculate the mean value for the over threshold cluster
  85 +
  86 + var_b.push_back(n_b * n_o * pow((m_b - m_o), 2));
  87 + }
  88 +
  89 + vector<float>::iterator max_var = std::max_element(var_b.begin(), var_b.end()); //finding maximum elements in the vector
  90 + size_t th_idx = std::distance(var_b.begin(), max_var);
  91 + T threshold = Imin + (T)(th_idx * th_step);
  92 + return threshold;
  93 + }
  94 +
  95 + //this function performs the 2D iterative voting algorithm on the image stored in the gpu
  96 + template<typename T>
  97 + void gpu_ivote2(T* gpuI, unsigned int rmax, size_t x, size_t y, bool invert = false, T t = 0, std::string outname_img = "out.bmp", std::string outname_txt = "out.txt",
  98 + int iter = 8, T phi = 15.0f * (float)stim::PI / 180, int conn = 8, bool debug = false) {
  99 +
  100 + size_t pixels = x * y; //compute the size of input image
  101 + //
  102 + if (invert) { //if inversion is required call the kernel to invert the image
  103 + unsigned int max_threads = stim::maxThreadsPerBlock();
  104 + dim3 threads((unsigned int)sqrt(max_threads), (unsigned int)sqrt(max_threads));
  105 + dim3 blocks((unsigned int)x / threads.x + 1, (unsigned int)y / threads.y + 1);
  106 + cuda_invert << <blocks, threads >> > (gpuI, x, y);
  107 + }
  108 + //
  109 + size_t table_bytes = (size_t)(pow(2 * rmax + 1, 2) * sizeof(T)); // create the atan2 table
  110 + T* cpuTable = (T*)malloc(table_bytes); //assign memory on the cpu for atan2 table
  111 + atan_2<T>(cpuTable, rmax); //call the function to precompute the atan2 table
  112 + T* gpuTable; HANDLE_ERROR(cudaMalloc(&gpuTable, table_bytes));
  113 + HANDLE_ERROR(cudaMemcpy(gpuTable, cpuTable, table_bytes, cudaMemcpyHostToDevice)); //copy atan2 table to the gpu
  114 +
  115 + size_t bytes = pixels* sizeof(T); //calculate the bytes of the input
  116 + float dphi = phi / iter; //change in phi for each iteration
  117 +
  118 + float* gpuGrad; HANDLE_ERROR(cudaMalloc(&gpuGrad, bytes * 2)); //allocate space to store the 2D gradient
  119 + float* gpuVote; HANDLE_ERROR(cudaMalloc(&gpuVote, bytes)); //allocate space to store the vote image
  120 +
  121 + stim::cuda::gpu_gradient_2d<float>(gpuGrad, gpuI, x, y); //calculate the 2D gradient
  122 + stim::cuda::gpu_cart2polar<float>(gpuGrad, x, y); //convert cartesian coordinate of gradient to the polar
  123 +
  124 + for (int i = 0; i < iter; i++) { //for each iteration
  125 + cudaMemset(gpuVote, 0, bytes); //reset the vote image to 0
  126 + stim::cuda::gpu_vote<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y, debug); //perform voting
  127 + stim::cuda::gpu_update_dir<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y, debug); //update the voter directions
  128 + phi = phi - dphi; //decrement phi
  129 + }
  130 + stim::cuda::gpu_local_max<float>(gpuI, gpuVote, conn, x, y); //calculate the local maxima
  131 +
  132 + T* pts = (T*)malloc(bytes); //allocate memory on the cpu to store the output of iterative voting
  133 + HANDLE_ERROR(cudaMemcpy(pts, gpuI, bytes, cudaMemcpyDeviceToHost)); //copy the output from gpu to the cpu memory
  134 +
  135 + T threshold;
  136 + if (t == 0) threshold = stim::th_otsu<T>(pts, pixels); //if threshold value is not set call the function to compute the threshold
  137 + else threshold = t;
  138 +
  139 + std::ofstream output; //save the thresholded detected seeds in a text file
  140 + output.open(outname_txt);
  141 + output << "X" << " " << "Y" << " " << "threshold" << "\n";
  142 + size_t ind;
  143 + for (size_t ix = 0; ix < x; ix++) {
  144 + for (size_t iy = 0; iy < y; iy++) {
  145 + ind = iy * x + ix;
  146 + if (pts[ind] > threshold) {
  147 + output << ix << " " << iy << " " << pts[ind] << "\n";
  148 + pts[ind] = 1;
  149 + }
  150 + else pts[ind] = 0;
  151 + }
  152 + }
  153 + output.close();
  154 +
  155 + HANDLE_ERROR(cudaMemcpy(gpuI, pts, bytes, cudaMemcpyHostToDevice)); //copy the points to the gpu
  156 + stim::cpu2image(pts, outname_img, x, y); //output the image
  157 +
  158 + }
  159 +
  160 +
  161 + template<typename T>
  162 + void cpu_ivote2(T* cpuI, unsigned int rmax, size_t x, size_t y, bool invert = false, T t = 0, std::string outname_img = "out.bmp", std::string outname_txt = "out.txt",
  163 + int iter = 8, T phi = 15.0f * (float)stim::PI / 180, int conn = 8, bool debug = false) {
  164 + size_t bytes = x*y * sizeof(T);
  165 + T* gpuI; //allocate space on the gpu to save the input image
  166 + HANDLE_ERROR(cudaMalloc(&gpuI, bytes));
  167 + HANDLE_ERROR(cudaMemcpy(gpuI, cpuI, bytes, cudaMemcpyHostToDevice)); //copy the image to the gpu
  168 + stim::gpu_ivote2<T>(gpuI, rmax, x, y, invert, t, outname_img, outname_txt, iter, phi, conn, debug); //call the gpu version of the ivote
  169 + HANDLE_ERROR(cudaMemcpy(cpuI, gpuI, bytes, cudaMemcpyDeviceToHost)); //copy the output to the cpu
  170 + }
  171 +}
  172 +#endif
0 \ No newline at end of file 173 \ No newline at end of file
stim/iVote/ivote2/iter_vote2.cuh 0 → 100644
  1 +#ifndef STIM_CUDA_ITER_VOTE2_H
  2 +#define STIM_CUDA_ITER_VOTE2_H
  3 +
  4 +//extern bool DEBUG;
  5 +
  6 +#include "update_dir_bb.cuh"
  7 +#include "vote_atomic_bb.cuh"
  8 +
  9 +namespace stim{
  10 + namespace cuda{
  11 +
  12 + }
  13 +}
  14 +
  15 +
  16 +
  17 +#endif
0 \ No newline at end of file 18 \ No newline at end of file
stim/cuda/ivote/local_max.cuh renamed to stim/iVote/ivote2/local_max.cuh
stim/cuda/ivote/update_dir.cuh renamed to stim/iVote/ivote2/update_dir.cuh
stim/cuda/ivote/update_dir_bb.cuh renamed to stim/iVote/ivote2/update_dir_bb.cuh
@@ -97,7 +97,7 @@ namespace stim{ @@ -97,7 +97,7 @@ namespace stim{
97 } 97 }
98 98
99 template<typename T> 99 template<typename T>
100 - void gpu_update_dir(T* gpuVote, T* gpuGrad, T* gpuTable, T phi, unsigned int rmax, size_t x, size_t y){ 100 + void gpu_update_dir(T* gpuVote, T* gpuGrad, T* gpuTable, T phi, unsigned int rmax, size_t x, size_t y, bool DEBUG = false){
101 101
102 //calculate the number of bytes in the array 102 //calculate the number of bytes in the array
103 size_t bytes = x * y * sizeof(T); 103 size_t bytes = x * y * sizeof(T);
stim/cuda/ivote/update_dir_shared.cuh renamed to stim/iVote/ivote2/update_dir_shared.cuh
stim/cuda/ivote/update_dir_threshold_global.cuh renamed to stim/iVote/ivote2/update_dir_threshold_global.cuh
stim/cuda/ivote/vote.cuh renamed to stim/iVote/ivote2/vote.cuh
stim/cuda/ivote/vote_atomic.cuh renamed to stim/iVote/ivote2/vote_atomic.cuh
stim/cuda/ivote/vote_atomic_bb.cuh renamed to stim/iVote/ivote2/vote_atomic_bb.cuh
@@ -87,7 +87,7 @@ namespace stim{ @@ -87,7 +87,7 @@ namespace stim{
87 /// @param x and y are the spatial dimensions of the gradient image 87 /// @param x and y are the spatial dimensions of the gradient image
88 /// @param gradmag defines whether or not the gradient magnitude is taken into account during the vote 88 /// @param gradmag defines whether or not the gradient magnitude is taken into account during the vote
89 template<typename T> 89 template<typename T>
90 - void gpu_vote(T* gpuVote, T* gpuGrad, T* gpuTable, T phi, unsigned int rmax, size_t x, size_t y, bool gradmag = true){ 90 + void gpu_vote(T* gpuVote, T* gpuGrad, T* gpuTable, T phi, unsigned int rmax, size_t x, size_t y, bool DEBUG = false, bool gradmag = true){
91 unsigned int max_threads = stim::maxThreadsPerBlock(); 91 unsigned int max_threads = stim::maxThreadsPerBlock();
92 dim3 threads( (unsigned int)sqrt(max_threads), (unsigned int)sqrt(max_threads) ); 92 dim3 threads( (unsigned int)sqrt(max_threads), (unsigned int)sqrt(max_threads) );
93 dim3 blocks((unsigned int)x/threads.x + 1, (unsigned int)y/threads.y + 1); 93 dim3 blocks((unsigned int)x/threads.x + 1, (unsigned int)y/threads.y + 1);
@@ -96,7 +96,7 @@ namespace stim{ @@ -96,7 +96,7 @@ namespace stim{
96 if (DEBUG) std::cout<<"Shared Memory required: "<<shared_mem_req<<std::endl; 96 if (DEBUG) std::cout<<"Shared Memory required: "<<shared_mem_req<<std::endl;
97 size_t shared_mem = stim::sharedMemPerBlock(); 97 size_t shared_mem = stim::sharedMemPerBlock();
98 if(shared_mem_req > shared_mem){ 98 if(shared_mem_req > shared_mem){
99 - std::cout<<"Error: insufficient shared memory for this implementation of cuda_update_dir()."<<std::endl; 99 + std::cout<<"Error: insufficient shared memory for this implementation of cuda_vote()."<<std::endl;
100 exit(1); 100 exit(1);
101 } 101 }
102 102
stim/cuda/ivote/vote_atomic_shared.cuh renamed to stim/iVote/ivote2/vote_atomic_shared.cuh
stim/cuda/ivote/vote_shared.cuh renamed to stim/iVote/ivote2/vote_shared.cuh
stim/cuda/ivote/vote_shared_32-32.cuh renamed to stim/iVote/ivote2/vote_shared_32-32.cuh
stim/cuda/ivote/vote_threshold_global.cuh renamed to stim/iVote/ivote2/vote_threshold_global.cuh
stim/image/image.h
@@ -53,6 +53,10 @@ class image{ @@ -53,6 +53,10 @@ class image{
53 void allocate(){ 53 void allocate(){
54 unalloc(); 54 unalloc();
55 img = (T*) malloc( sizeof(T) * R[0] * R[1] * R[2] ); //allocate memory 55 img = (T*) malloc( sizeof(T) * R[0] * R[1] * R[2] ); //allocate memory
  56 + if (img == NULL) {
  57 + std::cout << "stim::image ERROR - failed to allocate memory for image" << std::endl;
  58 + exit(1);
  59 + }
56 } 60 }
57 61
58 void allocate(size_t x, size_t y, size_t c){ //allocate memory based on the resolution 62 void allocate(size_t x, size_t y, size_t c){ //allocate memory based on the resolution
@@ -228,6 +232,14 @@ public: @@ -228,6 +232,14 @@ public:
228 } 232 }
229 } 233 }
230 #endif 234 #endif
  235 + //Copy N data points from source to dest, casting while doing so
  236 + template<typename S, typename D>
  237 + void type_copy(S* source, D* dest, size_t N) {
  238 + if (typeid(S) == typeid(D)) //if both types are the same
  239 + memcpy(dest, source, N * sizeof(S)); //just use a memcpy
  240 + for (size_t n = 0; n < N; n++) //otherwise, iterate through each element
  241 + dest[n] = (D)source[n]; //copy and cast
  242 + }
231 /// Load an image from a file 243 /// Load an image from a file
232 void load(std::string filename){ 244 void load(std::string filename){
233 #ifdef USING_OPENCV 245 #ifdef USING_OPENCV
@@ -236,13 +248,15 @@ public: @@ -236,13 +248,15 @@ public:
236 std::cout<<"ERROR stim::image::load() - unable to find image "<<filename<<std::endl; 248 std::cout<<"ERROR stim::image::load() - unable to find image "<<filename<<std::endl;
237 exit(1); 249 exit(1);
238 } 250 }
  251 + int cv_type = cvImage.type();
239 int cols = cvImage.cols; 252 int cols = cvImage.cols;
240 int rows = cvImage.rows; 253 int rows = cvImage.rows;
241 int channels = cvImage.channels(); 254 int channels = cvImage.channels();
242 allocate(cols, rows, channels); //allocate space for the image 255 allocate(cols, rows, channels); //allocate space for the image
  256 + size_t img_bytes = bytes();
243 unsigned char* cv_ptr = (unsigned char*)cvImage.data; 257 unsigned char* cv_ptr = (unsigned char*)cvImage.data;
244 - if(C() == 1) //if this is a single-color image, just copy the data  
245 - memcpy(img, cv_ptr, bytes()); 258 + if (C() == 1) //if this is a single-color image, just copy the data
  259 + type_copy<unsigned char, T>(cv_ptr, img, size());
246 if(C() == 3) //if this is a 3-color image, OpenCV uses BGR interleaving 260 if(C() == 3) //if this is a 3-color image, OpenCV uses BGR interleaving
247 from_opencv(cv_ptr, X(), Y()); 261 from_opencv(cv_ptr, X(), Y());
248 #else 262 #else
stim/math/matrix.h
@@ -33,32 +33,58 @@ namespace stim{ @@ -33,32 +33,58 @@ namespace stim{
33 } 33 }
34 } 34 }
35 35
  36 + //class encapsulates a mat4 file, and can be used to write multiple matrices to a single mat4 file
  37 + class mat4file {
  38 + std::ofstream matfile;
  39 +
  40 + public:
  41 + /// Constructor opens a mat4 file for writing
  42 + mat4file(std::string filename) {
  43 + matfile.open(filename, std::ios::binary);
  44 + }
  45 +
  46 + bool is_open() {
  47 + return matfile.is_open();
  48 + }
  49 +
  50 + void close() {
  51 + matfile.close();
  52 + }
  53 +
  54 + bool writemat(char* data, std::string varname, size_t sx, size_t sy, mat4Format format) {
  55 + //save the matrix file here (use the mat4 function above)
  56 + //data format: https://maxwell.ict.griffith.edu.au/spl/matlab-page/matfile_format.pdf (page 32)
  57 +
  58 + int MOPT = 0; //initialize the MOPT type value to zero
  59 + int m = 0; //little endian
  60 + int o = 0; //reserved, always 0
  61 + int p = format;
  62 + int t = 0;
  63 + MOPT = m * 1000 + o * 100 + p * 10 + t; //calculate the type value
  64 + int mrows = (int)sx;
  65 + int ncols = (int)sy;
  66 + int imagf = 0; //assume real (for now)
  67 + varname.push_back('\0'); //add a null to the string
  68 + int namlen = (int)varname.size(); //calculate the name size
  69 +
  70 + size_t bytes = sx * sy * mat4Format_size(format);
  71 + matfile.write((char*)&MOPT, 4);
  72 + matfile.write((char*)&mrows, 4);
  73 + matfile.write((char*)&ncols, 4);
  74 + matfile.write((char*)&imagf, 4);
  75 + matfile.write((char*)&namlen, 4);
  76 + matfile.write((char*)&varname[0], namlen);
  77 + matfile.write((char*)data, bytes); //write the matrix data
  78 + return is_open();
  79 + }
  80 + };
  81 +
36 static void save_mat4(char* data, std::string filename, std::string varname, size_t sx, size_t sy, mat4Format format){ 82 static void save_mat4(char* data, std::string filename, std::string varname, size_t sx, size_t sy, mat4Format format){
37 - //save the matrix file here (use the mat4 function above)  
38 - //data format: https://maxwell.ict.griffith.edu.au/spl/matlab-page/matfile_format.pdf (page 32)  
39 -  
40 - int MOPT = 0; //initialize the MOPT type value to zero  
41 - int m = 0; //little endian  
42 - int o = 0; //reserved, always 0  
43 - int p = format;  
44 - int t = 0;  
45 - MOPT = m * 1000 + o * 100 + p * 10 + t; //calculate the type value  
46 - int mrows = (int)sx;  
47 - int ncols = (int)sy;  
48 - int imagf = 0; //assume real (for now)  
49 - varname.push_back('\0'); //add a null to the string  
50 - int namlen = (int)varname.size(); //calculate the name size  
51 -  
52 - size_t bytes = sx * sy * mat4Format_size(format);  
53 - std::ofstream outfile(filename, std::ios::binary);  
54 - outfile.write((char*)&MOPT, 4);  
55 - outfile.write((char*)&mrows, 4);  
56 - outfile.write((char*)&ncols, 4);  
57 - outfile.write((char*)&imagf, 4);  
58 - outfile.write((char*)&namlen, 4);  
59 - outfile.write((char*)&varname[0], namlen);  
60 - outfile.write((char*)data, bytes); //write the matrix data  
61 - outfile.close(); 83 + mat4file outfile(filename); //create a mat4 file object
  84 + if (outfile.is_open()) { //if the file is open
  85 + outfile.writemat(data, varname, sx, sy, format); //write the matrix
  86 + outfile.close(); //close the file
  87 + }
62 } 88 }
63 89
64 template <class T> 90 template <class T>
@@ -409,8 +435,29 @@ public: @@ -409,8 +435,29 @@ public:
409 } 435 }
410 } 436 }
411 437
412 - // saves the matrix as a Level-4 MATLAB file  
413 - void mat4(std::string filename, std::string name = std::string("unknown"), mat4Format format = mat4_float) { 438 + void raw(std::string filename) {
  439 + ofstream out(filename, std::ios::binary);
  440 + if (out) {
  441 + out.write((char*)data(), rows() * cols() * sizeof(T));
  442 + out.close();
  443 + }
  444 + }
  445 +
  446 + void mat4(stim::mat4file& file, std::string name = std::string("unknown"), mat4Format format = mat4_float) {
  447 + //make sure the matrix name is valid (only numbers and letters, with a letter at the beginning
  448 + for (size_t c = 0; c < name.size(); c++) {
  449 + if (name[c] < 48 || //if the character isn't a number or letter, replace it with '_'
  450 + (name[c] > 57 && name[c] < 65) ||
  451 + (name[c] > 90 && name[c] < 97) ||
  452 + (name[c] > 122)) {
  453 + name[c] = '_';
  454 + }
  455 + }
  456 + if (name[0] < 65 ||
  457 + (name[0] > 91 && name[0] < 97) ||
  458 + name[0] > 122) {
  459 + name = std::string("m") + name;
  460 + }
414 if (format == mat4_float) { 461 if (format == mat4_float) {
415 if (sizeof(T) == 4) format = mat4_float32; 462 if (sizeof(T) == 4) format = mat4_float32;
416 else if (sizeof(T) == 8) format = mat4_float64; 463 else if (sizeof(T) == 8) format = mat4_float64;
@@ -419,7 +466,40 @@ public: @@ -419,7 +466,40 @@ public:
419 exit(1); 466 exit(1);
420 } 467 }
421 } 468 }
422 - stim::save_mat4((char*)M, filename, name, rows(), cols(), format); 469 + //the name is now valid
  470 +
  471 + //if the size of the array is more than 100,000,000 elements, the matrix isn't supported
  472 + if (rows() * cols() > 100000000) { //break the matrix up into multiple parts
  473 + //mat4file out(filename); //create a mat4 object to write the matrix
  474 + if (file.is_open()) {
  475 + if (rows() < 100000000) { //if the size of the row is less than 100,000,000, split the matrix up by columns
  476 + size_t ncols = 100000000 / rows(); //calculate the number of columns that can fit in one matrix
  477 + size_t nmat = (size_t)std::ceil((double)cols() / (double)ncols); //calculate the number of matrices required
  478 + for (size_t m = 0; m < nmat; m++) { //for each matrix
  479 + std::stringstream ss;
  480 + ss << name << "_part_" << m + 1;
  481 + if (m == nmat - 1)
  482 + file.writemat((char*)(data() + m * ncols * rows()), ss.str(), rows(), cols() - m * ncols, format);
  483 + else
  484 + file.writemat((char*)(data() + m * ncols * rows()), ss.str(), rows(), ncols, format);
  485 + }
  486 + }
  487 + }
  488 + }
  489 + //call the mat4 subroutine
  490 + else
  491 + //stim::save_mat4((char*)M, filename, name, rows(), cols(), format);
  492 + file.writemat((char*)data(), name, rows(), cols(), format);
  493 + }
  494 +
  495 + // saves the matrix as a Level-4 MATLAB file
  496 + void mat4(std::string filename, std::string name = std::string("unknown"), mat4Format format = mat4_float) {
  497 + stim::mat4file matfile(filename);
  498 +
  499 + if (matfile.is_open()) {
  500 + mat4(matfile, name, format);
  501 + matfile.close();
  502 + }
423 } 503 }
424 }; 504 };
425 505
@@ -243,7 +243,7 @@ public: @@ -243,7 +243,7 @@ public:
243 return false; 243 return false;
244 } 244 }
245 245
246 -//#ifndef __NVCC__ 246 +#ifndef __NVCC__
247 /// Outputs the vector as a string 247 /// Outputs the vector as a string
248 std::string str() const{ 248 std::string str() const{
249 std::stringstream ss; 249 std::stringstream ss;
@@ -261,7 +261,7 @@ public: @@ -261,7 +261,7 @@ public:
261 261
262 return ss.str(); 262 return ss.str();
263 } 263 }
264 -//#endif 264 +#endif
265 265
266 size_t size(){ return 3; } 266 size_t size(){ return 3; }
267 267
stim/parser/filename.h
@@ -89,6 +89,11 @@ protected: @@ -89,6 +89,11 @@ protected:
89 absolute.push_back(relative[i]); 89 absolute.push_back(relative[i]);
90 } 90 }
91 } 91 }
  92 + else {
  93 + if (relative[0] == ".")
  94 + relative = std::vector<std::string>(relative.begin() + 1, relative.end());
  95 + absolute = relative;
  96 + }
92 } 97 }
93 98
94 /// Parses a directory string into a drive (NULL if not Windows) and list of hierarchical directories 99 /// Parses a directory string into a drive (NULL if not Windows) and list of hierarchical directories
stim/visualization/aabb3.h
@@ -6,14 +6,14 @@ @@ -6,14 +6,14 @@
6 6
7 namespace stim{ 7 namespace stim{
8 8
9 - template<typename T>  
10 - using aabb3 = aabbn<T, 3>;  
11 -/*/// Structure for a 3D axis aligned bounding box 9 + //template<typename T>
  10 + //using aabb3 = aabbn<T, 3>;
  11 +/// Structure for a 3D axis aligned bounding box
12 template<typename T> 12 template<typename T>
13 struct aabb3 : public aabbn<T, 3>{ 13 struct aabb3 : public aabbn<T, 3>{
14 14
15 - aabb3() : aabbn() {}  
16 - aabb3(T x0, T y0, T z0, T x1, T y1, T z1){ 15 + CUDA_CALLABLE aabb3() : aabbn() {}
  16 + CUDA_CALLABLE aabb3(T x0, T y0, T z0, T x1, T y1, T z1){
17 low[0] = x0; 17 low[0] = x0;
18 low[1] = y0; 18 low[1] = y0;
19 low[2] = z0; 19 low[2] = z0;
@@ -22,11 +22,39 @@ struct aabb3 : public aabbn&lt;T, 3&gt;{ @@ -22,11 +22,39 @@ struct aabb3 : public aabbn&lt;T, 3&gt;{
22 high[2] = x2; 22 high[2] = x2;
23 } 23 }
24 24
25 - aabb3 aabbn<T, 3>() { 25 + CUDA_CALLABLE aabb3(T x, T y, T z) {
  26 + low[0] = high[0] = x;
  27 + low[1] = high[1] = y;
  28 + low[2] = high[2] = z;
  29 + }
  30 +
  31 + CUDA_CALLABLE void insert(T x, T y, T z) {
  32 + T p[3];
  33 + p[0] = x;
  34 + p[1] = y;
  35 + p[2] = z;
  36 + aabbn<T, 3>::insert(p);
  37 + }
26 38
  39 + CUDA_CALLABLE void trim_low(T x, T y, T z) {
  40 + T p[3];
  41 + p[0] = x;
  42 + p[1] = y;
  43 + p[2] = z;
  44 + aabbn<T, 3>::trim_low(p);
27 } 45 }
28 46
29 -};*/ 47 + CUDA_CALLABLE void trim_high(T x, T y, T z) {
  48 + T p[3];
  49 + p[0] = x;
  50 + p[1] = y;
  51 + p[2] = z;
  52 + aabbn<T, 3>::trim_high(p);
  53 + }
  54 +
  55 +
  56 +
  57 +};
30 58
31 } 59 }
32 60
stim/visualization/aabbn.h
@@ -25,26 +25,58 @@ struct aabbn{ @@ -25,26 +25,58 @@ struct aabbn{
25 init(i); 25 init(i);
26 } 26 }
27 27
  28 + /// For even inputs to the constructor, the input could be one point or a set of pairs of points
28 CUDA_CALLABLE aabbn(T x0, T x1) { 29 CUDA_CALLABLE aabbn(T x0, T x1) {
29 - low[0] = x0;  
30 - high[0] = x1; 30 + if (D == 1) {
  31 + low[0] = x0;
  32 + high[0] = x1;
  33 + }
  34 + else if (D == 2) {
  35 + low[0] = high[0] = x0;
  36 + low[1] = high[1] = x1;
  37 + }
31 } 38 }
32 39
  40 + /// In the case of 3 inputs, this must be a 3D bounding box, so initialize to a box of size 0 at (x, y, z)
  41 + /*CUDA_CALLABLE aabbn(T x, T y, T z) {
  42 + low[0] = high[0] = x;
  43 + low[1] = high[1] = y;
  44 + low[2] = high[2] = z;
  45 + }*/
  46