codebase / stim-matlab

Browse Code »

Commit 8be1ab930be4cf927f46b79c654a102b1c9bd8fc

Authored by David Mayerich 2015-04-07 13:24:08 -0500

0 parents

initial commit of the new Matlab repository

Showing 31 changed files with 1062 additions and 0 deletions Show diff stats

Inline Side-by-side

rtsApplyBrewer.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsApplyBrewer.m
	1	+function result = rtsApplyBrewer(field, minVal, maxVal)
	2	+
	3	+if nargin < 3
	4	+ maxVal = max(field(:));
	5	+end
	6	+if nargin < 2
	7	+ minVal = min(field(:));
	8	+end
	9	+
	10	+fieldMaskMin = field >= minVal;
	11	+fieldMaskMax = field <= maxVal;
	12	+
	13	+fieldMask = min(fieldMaskMin, fieldMaskMax);
	14	+
	15	+ctrlPts = zeros(11, 3);
	16	+
	17	+ctrlPts(1, :) = [0.192157, 0.211765, 0.584314];
	18	+ctrlPts(2, :) = [0.270588, 0.458824, 0.705882];
	19	+ctrlPts(3, :) = [0.454902, 0.678431, 0.819608];
	20	+ctrlPts(4, :) = [0.670588, 0.85098, 0.913725];
	21	+ctrlPts(5, :) = [0.878431, 0.952941, 0.972549];
	22	+ctrlPts(6, :) = [1, 1, 0.74902];
	23	+ctrlPts(7, :) = [0.996078, 0.878431, 0.564706];
	24	+ctrlPts(8, :) = [0.992157, 0.682353, 0.380392];
	25	+ctrlPts(9, :) = [0.956863, 0.427451, 0.262745];
	26	+ctrlPts(10, :) = [0.843137, 0.188235, 0.152941];
	27	+ctrlPts(11, :) = [0.647059, 0, 0.14902];
	28	+
	29	+%X = 0:1/10:1;
	30	+X = minVal:(maxVal - minVal)/10:maxVal;
	31	+
	32	+R = interp1(X, ctrlPts(:, 1), field);
	33	+G = interp1(X, ctrlPts(:, 2), field);
	34	+B = interp1(X, ctrlPts(:, 3), field);
	35	+
	36	+R = R.* fieldMask;
	37	+G = G.* fieldMask;
	38	+B = B.* fieldMask;
	39	+
	40	+result = cat(ndims(field)+1, R, G, B);
	41	+
...	...

rtsBaselineCorrect.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsBaselineCorrect.m
	1	+function R = rtsBaselineCorrect(A, b, dim)
	2	+
	3	+%This function baselines an n-dimensional array given a series of points b
	4	+%B = baseline-corrected data
	5	+%A = raw data
	6	+%b = baseline points in the format [b0 b1 ... bn]
	7	+%dim = dimension along which the baseline correction will occur
	8	+
	9	+
	10	+%allocate space for the sub-component indices
	11	+sz = size(A);
	12	+inds = repmat({1},1,ndims(A));
	13	+
	14	+%allocate space for the result
	15	+R = zeros(size(A));
	16	+
	17	+%for each baseline point
	18	+for bp = 1:length(b) - 1
	19	+
	20	+ %extract the data at the first baseline point
	21	+ n = b(bp);
	22	+ inds = getIndices(A, n, dim);
	23	+ BP1 = A(inds{:});
	24	+
	25	+ %extract the data at the second baseline point
	26	+ n = b(bp+1);
	27	+ inds = getIndices(A, n, dim);
	28	+ BP2 = A(inds{:});
	29	+
	30	+ %for each point in between the baseline points
	31	+ for p = b(bp):b(bp+1)
	32	+
	33	+ %compute the weighting of the linear function between BP1 and BP2
	34	+ a2 = (p - b(bp))/(b(bp+1) - b(bp));
	35	+ a1 = 1.0 - a2;
	36	+
	37	+ n = p;
	38	+ inds = getIndices(A, n, dim);
	39	+
	40	+ baseFunc = a1 * BP1 + a2 * BP2;
	41	+ R(inds{:}) = A(inds{:}) - baseFunc;
	42	+ end
	43	+end
	44	+
	45	+
	46	+function I = getIndices(A, n, dim)
	47	+
	48	+%this function returns indices for all elements in S at
	49	+%point n along dimension dim
	50	+
	51	+sz = size(A);
	52	+I = repmat({1},1,ndims(A));
	53	+
	54	+for d = 1:ndims(A)
	55	+ if d ~= dim
	56	+ I{d} = 1:sz(d);
	57	+ else
	58	+ I{d} = n;
	59	+ end
	60	+end
	61	+
...	...

rtsClass2Color.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsClass2Color.m
	1	+function C = rtsClass2Color(classes, classList, colorList)
	2	+
	3	+C = zeros(length(classes), 3);
	4	+
	5	+for i = 1:length(classes)
	6	+
	7	+ for c = 1:length(classList)
	8	+
	9	+ if classes(i) == classList(c)
	10	+ C(i, :) = colorList(c, :);
	11	+ end
	12	+
	13	+ end
	14	+
	15	+end
0	16	\ No newline at end of file
...	...

rtsColorMapField.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsColorMapField.m
	1	+function Result = rtsColorMapField(F, range)
	2	+%creates a color map for the field F and returns an array with
	3	+%an additional size-3 color dimension
	4	+
	5	+R = zeros(size(F));
	6	+G = zeros(size(F));
	7	+B = zeros(size(F));
	8	+
	9	+%if no range is specified, estimate it using min and max values
	10	+if nargin == 1
	11	+ range(1) = min(F(:));
	12	+ range(2) = max(F(:));
	13	+end
	14	+range_size = range(2) - range(1);
	15	+
	16	+num_elements = size(F(:));
	17	+for i =1:num_elements
	18	+ if(F(i) > range(1) && F(i) < range(2))
	19	+ hue_degrees = (1 - (F(i) - range(1))/range_size)*240;
	20	+ HSVval = [hue_degrees 1.0 1.0];
	21	+ RGBval = rtsHSVtoRGB(HSVval);
	22	+ else
	23	+ RGBval = [0 0 F(i)/range(1)];
	24	+ end
	25	+ R(i) = RGBval(1);
	26	+ G(i) = RGBval(2);
	27	+ B(i) = RGBval(3);
	28	+end
	29	+Result = cat(ndims(F)+1, R, G, B);
	30	+imshow(Result);
...	...

rtsColorMapRainbow.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsColorMapRainbow.m
	1	+function result = rtsColorMapRainbow(field, minVal, maxVal)
	2	+
	3	+if nargin < 3
	4	+ maxVal = max(field(:));
	5	+end
	6	+if nargin < 2
	7	+ minVal = min(field(:));
	8	+end
	9	+
	10	+fieldMaskMin = field >= minVal;
	11	+fieldMaskMax = field <= maxVal;
	12	+
	13	+fieldMask = fieldMaskMin;%min(fieldMaskMin, fieldMaskMax);
	14	+
	15	+ctrlPts = zeros(3, 3);
	16	+
	17	+ctrlPts(1, :) = [0.0, 0.0, 0.5];
	18	+ctrlPts(2, :) = [0.0, 0.0, 1.0];
	19	+ctrlPts(3, :) = [0.0, 1.0, 1.0];
	20	+ctrlPts(4, :) = [0.0, 1.0, 0.0];
	21	+ctrlPts(5, :) = [1.0, 1.0, 0.0];
	22	+ctrlPts(6, :) = [1.0, 0.0, 0.0];
	23	+ctrlPts(7, :) = [0.5, 0.0, 0.0];
	24	+
	25	+%X = 0:1/10:1;
	26	+X = minVal:(maxVal - minVal)/6:maxVal;
	27	+
	28	+R = interp1(X, ctrlPts(:, 1), field);
	29	+G = interp1(X, ctrlPts(:, 2), field);
	30	+B = interp1(X, ctrlPts(:, 3), field);
	31	+
	32	+%R = R.* fieldMask;
	33	+%G = G.* fieldMask;
	34	+%B = B.* fieldMask;
	35	+
	36	+result = cat(ndims(field)+1, R, G, B);
	37	+
...	...

rtsCropMask.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsCropMask.m
	1	+function new_mask = rtsCropMask(mask, max_pixels)
	2	+
	3	+%convert to logical
	4	+%mask = mask > 0;
	5	+
	6	+if nnz(mask) < max_pixels
	7	+ new_mask = mask;
	8	+else
	9	+ j = 0;
	10	+ for i = 1:length(mask(:))
	11	+ if mask(i) > 0
	12	+ j = j + 1;
	13	+ if j == max_pixels
	14	+ last_pixel = i;
	15	+ end
	16	+ end
	17	+ end
	18	+ new_mask = zeros(size(mask, 1), size(mask, 2));
	19	+ new_mask(1:last_pixel) = mask(1:last_pixel);
	20	+end
	21	+
0	22	\ No newline at end of file
...	...

rtsEnviClassify.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviClassify.m
	1	+function C = rtsEnviClassify(B, envi_filename, mask, nFeatures)
	2	+
	3	+% B = TreeBagger class
	4	+% nFeatures = number of features cropped out (0 means use all features)
	5	+
	6	+if nargin < 4
	7	+ nFeatures = 0;
	8	+end
	9	+
	10	+
	11	+batchSize = 100000;
	12	+
	13	+
	14	+
	15	+Ns = nnz(mask) %get the number of samples to be classified
	16	+Nb = ceil(Ns / batchSize); %get the number of batches
	17	+
	18	+C(Ns, 1) = char(0);
	19	+
	20	+%open the ENVI file
	21	+envi = rtsEnviOpen(envi_filename, [envi_filename '.hdr'], mask);
	22	+
	23	+for b = 1:Nb
	24	+
	25	+ %load a batch
	26	+ batch = rtsEnviRead(envi, batchSize);
	27	+
	28	+ %crop out the number of features specified (if it is specified)
	29	+ if nFeatures ~= 0
	30	+ batch = batch(1:nFeatures, :);
	31	+ end
	32	+
	33	+ %classify the batch
	34	+ C( (b-1)batchSize + 1 : (b-1)batchSize + size(batch, 2) ) = cell2mat(predict(B, batch'));
	35	+
	36	+ disp([num2str( round( b / Nb * 100 ) ) '%']);
	37	+end
	38	+
	39	+rtsEnviClose(envi);
0	40	\ No newline at end of file
...	...

rtsEnviClose.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviClose.m
	1	+function rtsEnviClose(enviFID)
	2	+
	3	+fclose(enviFID.fid);
0	4	\ No newline at end of file
...	...

rtsEnviID.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviID.m
	1	+classdef rtsEnviID < handle
	2	+ properties
	3	+ fid
	4	+ header
	5	+ mask
	6	+ fpos
	7	+ end
	8	+end
0	9	\ No newline at end of file
...	...

rtsEnviLoadClasses.m~ 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviLoadClasses.m~
	1	+function [F, C] = rtsEnviLoadClasses(filenames, classnames, ppc)
	2	+
	3	+Nc = size(filenames, 1);
	4	+
	5	+for c = 1:Nc
	6	+ %load the masks
	7	+ image = imread(filenames{c});
	8	+ mask = image(:, :, 1)' > 0;
	9	+
	10	+
	11	+end
0	12	\ No newline at end of file
...	...

rtsEnviLoadHeader.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviLoadHeader.m
	1	+function s = rtsEnviLoadHeader(filename)
	2	+%create a cell array of fields
	3	+s = struct;
	4	+
	5	+fid = fopen(filename, 'r', 'n', 'US-ASCII');
	6	+
	7	+%read the first field and make sure it is "ENVI"
	8	+fname = GetFieldName(fid);
	9	+if strcmp(fname, 'ENVI') == 0
	10	+ disp('Not an ENVI header file');
	11	+ return;
	12	+end
	13	+
	14	+while feof(fid) == 0
	15	+ fname = GetFieldName(fid);
	16	+ if feof(fid) == 1
	17	+ return;
	18	+ end
	19	+ [value, valid] = ReadField(fid, fname);
	20	+ if valid == 1
	21	+ s = setfield(s, fname, value);
	22	+ end
	23	+end
	24	+fclose(fid);
	25	+
	26	+function t = GetFieldName(fid)
	27	+string_struct = textscan(fid, '%s', 1, 'Delimiter', '=');
	28	+if feof(fid) == 1
	29	+ t = [];
	30	+ return;
	31	+end
	32	+t = string_struct{1}{1};
	33	+t = strtrim(t);
	34	+t(t==' ') = '_';
	35	+
	36	+function [v, valid] = ReadField(fid, field_name)
	37	+valid = 1;
	38	+stringFields = {'file_type', 'interleave', 'sensor_type', 'wavelength_units'};
	39	+intFields = {'samples', 'lines', 'bands', 'header_offset', 'data_type', 'byte_order'};
	40	+
	41	+%if the field is "description", read between the brackets
	42	+if strcmp(field_name, 'description') == 1
	43	+ textscan(fid, '%[{]', 1);
	44	+ string_struct = textscan(fid, '%[^}]', 1, 'Whitespace', '');
	45	+ textscan(fid, '%[}]', 1);
	46	+ v = string_struct{1}{1};
	47	+ v = strtrim(v);
	48	+ return;
	49	+end
	50	+if max(strcmp(field_name, intFields)) ~= 0
	51	+ v = fscanf(fid, '%d');
	52	+ return;
	53	+end
	54	+if max(strcmp(field_name, stringFields)) ~= 0
	55	+ string_struct = textscan(fid, '%s', 1, 'Whitespace', '\n');
	56	+ v = string_struct{1}{1};
	57	+ v = strtrim(v);
	58	+ return;
	59	+end
	60	+
	61	+%read and return the wavelength values
	62	+if strcmp(field_name, 'wavelength') == 1
	63	+ v = [];
	64	+ textscan(fid, '%[{]', 1);
	65	+ c = ' ';
	66	+ while c ~= '}'
	67	+ new = fscanf(fid, '%f');
	68	+ v = [v new];
	69	+ c = fscanf(fid, '%c', 1);
	70	+ end
	71	+ return;
	72	+end
	73	+
	74	+%if it doesn't match anything, just read until the end of the line
	75	+%string_struct = textscan(fid, '%s', 1, 'Whitespace', '\n');
	76	+string_struct = textscan(fid, '%s', 1, 'Delimiter', '\n');
	77	+v = '';
	78	+valid = 0;
...	...

rtsEnviLoadImage.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviLoadImage.m
	1	+function [image, header] = rtsEnviLoadImage(filename, headername)
	2	+%enter the file name (no .hdr extension)
	3	+
	4	+header_file = headername;
	5	+header = rtsEnviLoadHeader(header_file);
	6	+
	7	+if strcmp(header.interleave, 'bsq') == 1
	8	+ image = zeros(header.samples, header.lines, header.bands);
	9	+end
	10	+if strcmp(header.interleave, 'bip') == 1
	11	+ image = zeros(header.bands, header.samples, header.lines);
	12	+end
	13	+if strcmp(header.interleave, 'bil') == 1
	14	+ image = zeros(header.samples, header.bands, header.lines);
	15	+end
	16	+
	17	+
	18	+file_bytes = header.samplesheader.linesheader.bands;
	19	+%check to make sure there is enough memory available
	20	+if ispc
	21	+ [~, sys] = memory;
	22	+ if sys.PhysicalMemory.Available < file_bytes*header.data_type && strcmp(header.interleave, 'bsq') == 0
	23	+ strResponse = input('The data set will require virtual memory for permutation. Continue? (y/n) ', 's');
	24	+ if strcmp(strResponse, 'n') == 1
	25	+ return;
	26	+ end
	27	+ end
	28	+else
	29	+ disp('Sorry, I can''t check available memory for this OS. Keep your fingers crossed!');
	30	+end
	31	+
	32	+fid = fopen(filename);
	33	+%skip the header - start reading header_offset bytes from the beginning
	34	+fseek(fid, header.header_offset, 'bof');
	35	+
	36	+%read the data
	37	+image(:) = fread(fid, file_bytes, 'float32');
	38	+
	39	+if strcmp(header.interleave, 'bip') == 1
	40	+ image = permute(image, [2 3 1]);
	41	+end
	42	+if strcmp(header.interleave, 'bil') == 1
	43	+ image = permute(image, [1 3 2]);
	44	+end
0	45	\ No newline at end of file
...	...

rtsEnviLoadTraining.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviLoadTraining.m
	1	+function [F, C] = rtsEnviLoadTraining(filenames, classnames, binfile, ppc)
	2	+
	3	+Nc = length(filenames); %calculate the number of classes
	4	+
	5	+if nargin < 3
	6	+ ppc = 0;
	7	+end
	8	+
	9	+%---------Load the Masks-----------
	10	+%first, get the total number of pixels
	11	+Ns = 0;
	12	+for c = 1:Nc
	13	+ %load the mask
	14	+ image = imread(filenames{c});
	15	+
	16	+ if c == 1
	17	+ mask = boolean(zeros(size(image, 2), size(image, 1), Nc));
	18	+ end
	19	+
	20	+ if ppc ~= 0
	21	+ mask(:, :, c) = rtsRandomizeMask(image(:, :, 1)' > 0, ppc);
	22	+ else
	23	+ mask(:, :, c) = image(:, :, 1)' > 0;
	24	+ end
	25	+
	26	+end
	27	+
	28	+%calculate the number of samples
	29	+Ns = Ns + nnz(mask);
	30	+
	31	+%-----------Load the Training Data-------------
	32	+disp(['loading ' num2str(Ns) ' samples']);
	33	+ci = 1;
	34	+for c = 1:Nc
	35	+
	36	+ disp(['loading class ' classnames(c) '...']);
	37	+ batch = nnz(mask(:, :, c));
	38	+ %create an ENVI file object
	39	+ envi = rtsEnviOpen(binfile, [binfile '.hdr'], mask(:, :, c));
	40	+
	41	+ if c == 1
	42	+ F = zeros(Ns, envi.header.bands);
	43	+ C(Ns, 1) = char(0);
	44	+ end
	45	+
	46	+ %read a bunch of spectra
	47	+ F(ci:ci+batch - 1, :) = rtsEnviRead(envi, batch)';
	48	+ C(ci:ci+batch - 1, 1) = repmat([classnames(c)], [batch 1]);
	49	+
	50	+ %update the current index
	51	+ ci = ci+batch;
	52	+
	53	+ %close the ENVI file
	54	+ rtsEnviClose(envi);
	55	+
	56	+ disp('done');
	57	+end
0	58	\ No newline at end of file
...	...

rtsEnviOpen.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviOpen.m
	1	+function EID = rtsEnviOpen(filename, headername, mask)
	2	+
	3	+%opens an ENVI file and returns an ID structure
	4	+% filename = name of ENVI file
	5	+% headername = name of ENVI header
	6	+% FID = structure containing file information
	7	+% mask = binary image specifying valid spectra
	8	+
	9	+%assign the mask variable
	10	+if(nargin < 3)
	11	+ mask = 1;
	12	+end
	13	+
	14	+%open the file and save the file ID
	15	+fid = fopen(filename, 'r');
	16	+
	17	+%open the ENVI header file
	18	+header = rtsEnviLoadHeader(headername);
	19	+
	20	+%this is currently only valid for BIP files
	21	+if(~strcmp(header.interleave, 'bip'))
	22	+ disp('Error: This function only works for BIP interleave files. Load in ENVI and convert.');
	23	+end
	24	+if(header.data_type ~= 4)
	25	+ disp('Error: This function only works for floating-point files.');
	26	+end
	27	+
	28	+EID = rtsEnviID;
	29	+EID.fid = fid;
	30	+EID.header = header;
	31	+EID.mask = mask;
	32	+EID.fpos = 0;
	33	+
	34	+%EID = struct('fid', fid, 'header', header, 'mask', mask, 'fpos', 0);
	35	+
...	...

rtsEnviRandomForest2C_apply.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviRandomForest2C_apply.m
	1	+function P = rtsEnviRandomForest2C_apply(RF, EnviFileName, EnviHeaderName, TissueMask, reference)
	2	+
	3	+%Applies a 2-class random forest classifier to the given image. Only
	4	+%pixels specified in TissueMask are classified. The estimated rule data is
	5	+%returned as an image sequence.
	6	+
	7	+%load the tissue mask
	8	+maskImage = imread(TissueMask);
	9	+maskBinary = (maskImage(:, :, 1) > 0)';
	10	+
	11	+nPixels = nnz(maskBinary);
	12	+disp(['Number of pixels: ' num2str(nPixels)]);
	13	+
	14	+
	15	+%load the data
	16	+disp(['Loading Features: ' EnviFileName]);
	17	+tLoadTime = tic;
	18	+fid = rtsEnviOpen(EnviFileName, EnviHeaderName, maskBinary);
	19	+F = rtsEnviRead(fid, nPixels);
	20	+rtsEnviClose(fid);
	21	+disp(['Time: ' num2str(toc(tLoadTime)) 's']);
	22	+
	23	+%transpose
	24	+F = F';
	25	+
	26	+%apply the reference
	27	+if nargin == 5
	28	+ Fnorm = repmat(F(:, reference), 1, size(F, 2));
	29	+ F = F./Fnorm;
	30	+ F(:, reference) = [];
	31	+end
	32	+
	33	+%classify the data
	34	+disp('Classifying...');
	35	+tClassTime = tic;
	36	+T = predict(RF, F);
	37	+disp(['Time: ' num2str(toc(tClassTime)) 's']);
	38	+%Cl = Cpost > threshold;
	39	+
	40	+%reconstruct the image from the mask
	41	+P = rtsMaskReconstruct(T', maskBinary, -1);
0	42	\ No newline at end of file
...	...

rtsEnviRandomForest2C_train.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviRandomForest2C_train.m
	1	+function RF = rtsEnviRandomForest2C_train(EnviFileName, EnviHeaderName, RoiImages, reference)
	2	+
	3	+%Creates a 2-class random forest classifier for an Envi BIP file using the
	4	+%provided ROI images. The resulting classifier uses regression to map
	5	+%class A to 1 and class B to 0. This allows the creation of ROC curves.
	6	+%
	7	+% EnviFilename - Name of an ENVI BIP file
	8	+% RoiImages - Cell array containing names of ROI images (masks)
	9	+
	10	+%default parameters
	11	+maxPixels = 200000;
	12	+threshold = 0.5;
	13	+nTrees = 100;
	14	+nThreads = 8;
	15	+%trainPixels = 200000;
	16	+
	17	+%determine the number of classes
	18	+nClasses = length(RoiImages);
	19	+%make sure that there are only two classes
	20	+if nClasses > 2
	21	+ disp('This classifier only supports 2 classes.');
	22	+ return;
	23	+end
	24	+
	25	+%for each class, load the training data
	26	+T = [];
	27	+F = [];
	28	+for c = 1:nClasses
	29	+
	30	+ %load the class mask
	31	+ maskImage = imread(RoiImages{c});
	32	+ maskBinary = (maskImage(:, :, 1) > 0)';
	33	+
	34	+ disp('------------------------------------------------------');
	35	+ %load epithelium spectra
	36	+ disp(['Loading Training Class ' num2str(c) ' pixels: ' EnviFileName]);
	37	+ tLoadTime = tic;
	38	+ fid = rtsEnviOpen(EnviFileName, EnviHeaderName, maskBinary);
	39	+ Fc = rtsEnviRead(fid, maxPixels);
	40	+ rtsEnviClose(fid);
	41	+
	42	+ if c == 1
	43	+ Tc = ones(size(Fc, 2), 1);
	44	+ else
	45	+ Tc = zeros(size(Fc, 2), 1);
	46	+ end
	47	+
	48	+ disp(['Time: ' num2str(toc(tLoadTime)) 's']);
	49	+
	50	+ %add features and targets to the final vectors
	51	+ T = [T; Tc];
	52	+ F = [F; Fc'];
	53	+end
	54	+
	55	+%apply the reference
	56	+if nargin == 4
	57	+ Fnorm = repmat(F(:, reference), 1, size(F, 2));
	58	+ F = F./Fnorm;
	59	+ F(:, reference) = [];
	60	+end
	61	+
	62	+%parallelize if specified
	63	+if nThreads > 1
	64	+ matlabpool('open', nThreads);
	65	+ paraoptions = statset('UseParallel', 'always');
	66	+end
	67	+
	68	+%train the classifier
	69	+disp('Creating a Random Forest classifier...');
	70	+tTrainTime = tic;
	71	+RF = TreeBagger(nTrees, F, T, 'method', 'regression', 'Options',paraoptions);
	72	+disp(['Time: ' num2str(toc(tTrainTime)) 's']);
	73	+
	74	+if nThreads > 1
	75	+ matlabpool close
	76	+end
0	77	\ No newline at end of file
...	...

rtsEnviRandomForest2C_validate.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviRandomForest2C_validate.m
	1	+function [Croc, Cauc] = rtsEnviRandomForest2C_validate(RF, EnviFileName, EnviHeaderName, RoiImages, reference)
	2	+
	3	+%Validates a 2-class random forest classifier, creating ROC curves.
	4	+
	5	+%parameters
	6	+maxPixels = 200000;
	7	+
	8	+%determine the number of classes
	9	+nClasses = length(RoiImages);
	10	+%make sure that there are only two classes
	11	+if nClasses > 2
	12	+ disp('This classifier only supports 2 classes.');
	13	+ return;
	14	+end
	15	+
	16	+%for each class, load the validation data
	17	+T = [];
	18	+F = [];
	19	+for c = 1:nClasses
	20	+
	21	+ %load the class mask
	22	+ maskImage = imread(RoiImages{c});
	23	+ maskBinary = (maskImage(:, :, 1) > 0)';
	24	+
	25	+ disp('------------------------------------------------------');
	26	+ %load epithelium spectra
	27	+ disp(['Loading Validation Class ' num2str(c) ' pixels: ' EnviFileName]);
	28	+ tLoadTime = tic;
	29	+ fid = rtsEnviOpen(EnviFileName, EnviHeaderName, maskBinary);
	30	+ Fc = rtsEnviRead(fid, maxPixels);
	31	+ rtsEnviClose(fid);
	32	+
	33	+ if c == 1
	34	+ Tc = ones(size(Fc, 2), 1);
	35	+ else
	36	+ Tc = zeros(size(Fc, 2), 1);
	37	+ end
	38	+
	39	+ disp(['Time: ' num2str(toc(tLoadTime)) 's']);
	40	+
	41	+ %add features and targets to the final vectors
	42	+ T = [T; Tc];
	43	+ F = [F; Fc'];
	44	+end
	45	+
	46	+%apply the reference
	47	+if nargin == 5
	48	+ Fnorm = repmat(F(:, reference), 1, size(F, 2));
	49	+ F = F./Fnorm;
	50	+ F(:, reference) = [];
	51	+end
	52	+
	53	+%classify the data (get the estimated posterior probability)
	54	+disp('Validating...');
	55	+tClassTime = tic;
	56	+Tpost = predict(RF, F);
	57	+disp(['Time: ' num2str(toc(tClassTime)) 's']);
	58	+
	59	+%calculate and display the ROC curve
	60	+disp('Calculating ROC...');
	61	+tRocTime = tic;
	62	+[Croc, Cauc, threshList] = rtsROC(Tpost, T);
	63	+disp(['Time: ' num2str(toc(tRocTime)) 's']);
	64	+disp(['AUC = ' num2str(Cauc)]);
	65	+plot(Croc(:, 1), Croc(:, 2));
0	66	\ No newline at end of file
...	...

rtsEnviRead.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviRead.m
	1	+function S = rtsEnviRead(fid, batchSize)
	2	+
	3	+%This function reads a batch of spectra from the given ENVI file ID
	4	+% fid = ENVI file ID created using rtsEnviOpen
	5	+% batchSize = number of spectra to read
	6	+
	7	+%if there is no mask, just load each spectrum in order
	8	+if(fid.mask == 1)
	9	+
	10	+ %compute the new batch size in case we are near the eof
	11	+ nSpectra = fid.header.samples * fid.header.lines;
	12	+ remainingSpectra = nSpectra - fid.fpos;
	13	+
	14	+ if(batchSize > remainingSpectra)
	15	+ batchSize = remainingSpectra;
	16	+ end
	17	+
	18	+ S = fread(fid.fid, [fid.header.bands batchSize], '*float32');
	19	+
	20	+ %increment the fid counter
	21	+ fid.fpos = fid.fpos + batchSize;
	22	+
	23	+%otherwise only load valid spectra
	24	+else
	25	+
	26	+ %compute the new batch size in case we are near the eof
	27	+ if(fid.fpos == 0)
	28	+ remainingSpectra = nnz(fid.mask);
	29	+ else
	30	+ nSpectra = nnz(fid.mask);
	31	+ maskRead = fid.mask(1:fid.fpos+1);
	32	+ remainingSpectra = nSpectra - nnz(maskRead);
	33	+ end
	34	+
	35	+ if(batchSize > remainingSpectra)
	36	+ batchSize = remainingSpectra;
	37	+ end
	38	+
	39	+ %allocate space for the spectra
	40	+ S = zeros(fid.header.bands, batchSize);
	41	+
	42	+ %for each spectrum in the batch
	43	+ for s = 1:batchSize
	44	+
	45	+ %while the current spectrum is invalid
	46	+ skip = 0;
	47	+ while (~fid.mask(fid.fpos + 1))
	48	+ %read the invalid spectrum
	49	+ %invalid = fread(fid.fid, [fid.header.bands 1], '*float32');
	50	+ skip = skip + 1;
	51	+
	52	+ %increment the file position
	53	+ fid.fpos = fid.fpos + 1;
	54	+ end
	55	+ fseek(fid.fid, skip * fid.header.bands * 4, 'cof');
	56	+
	57	+ test = fread(fid.fid, [fid.header.bands 1], '*float32');
	58	+ if size(test) ~= size(S(:, s))
	59	+ size(test)
	60	+ size(S(:, s))
	61	+ end
	62	+ S(:, s) = test;
	63	+ fid.fpos = fid.fpos + 1;
	64	+
	65	+ end
	66	+
	67	+end
	68	+
...	...

rtsEnviSaveImage.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsEnviSaveImage.m
	1	+function rtsEnviSaveImage(S, fileName, headerName, W)
	2	+
	3	+% S = 3D hyperspectral data
	4	+% W = list of wavelength values
	5	+
	6	+%save the image
	7	+rtsSaveRAW(S, fileName, 'float32', 'l');
	8	+
	9	+%create the header file
	10	+headerFile = [fileName '.hdr'];
	11	+if nargin == 3
	12	+ headerFile = headerName;
	13	+end
	14	+
	15	+fid = fopen(headerFile, 'w');
	16	+
	17	+%write the header and description
	18	+fprintf(fid, 'ENVI\n');
	19	+fprintf(fid, 'description = {\n');
	20	+fprintf(fid, ' File Saved using Matlab.}\n');
	21	+fprintf(fid, 'samples = %d\n', size(S, 1));
	22	+fprintf(fid, 'lines = %d\n', size(S, 2));
	23	+fprintf(fid, 'bands = %d\n', size(S, 3));
	24	+fprintf(fid, 'header offset = 0\n');
	25	+fprintf(fid, 'file type = ENVI Standard\n');
	26	+fprintf(fid, 'data type = 4\n');
	27	+fprintf(fid, 'interleave = bsq\n');
	28	+fprintf(fid, 'sensor type = Unknown\n');
	29	+fprintf(fid, 'byte order = 0\n');
	30	+
	31	+if nargin == 4
	32	+ %output wavelength units
	33	+ fprintf(fid, 'wavelength = {\n');
	34	+ for i = 1:length(W)-1
	35	+ fprintf(fid, ' %f,\n', W(i));
	36	+ end
	37	+ fprintf(fid, ' %f\n', W(length(W)));
	38	+ fprintf(fid, '}\n');
	39	+end
	40	+
	41	+fclose(fid);
...	...

rtsGaussianBlurND.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsGaussianBlurND.m
	1	+function R = rtsGaussianBlurND(input, sigma)
	2	+%This function convolves an ND array with a gaussian kernel specified by
	3	+%sigma. This takes advantage of the separability of the kernel and is
	4	+%therefore very fast.
	5	+
	6	+%disp('begin***********************************************');
	7	+
	8	+%get the number of dimensions for the convolution
	9	+dim = ndims(input);
	10	+
	11	+%initialize the result to the source
	12	+R = input;
	13	+
	14	+
	15	+for d=1:dim
	16	+ %if the dimension is greater than 1 and the sigma not a dirac
	17	+ if size(input, d) > 1 && sigma(d) > 0
	18	+ %create a 1D kernel
	19	+ kernelshape = ones(1, dim);
	20	+ kernelshape(d) = sigma(d)*6;
	21	+ kernel = zeros(kernelshape);
	22	+ kernel(:) = ndgauss(sigma(d)*6, sigma(d));
	23	+
	24	+ %perform the convolution
	25	+ R = convn(R, kernel, 'same');
	26	+ end
	27	+end
	28	+
	29	+
	30	+
0	31	\ No newline at end of file
...	...

rtsHSVtoRGB.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsHSVtoRGB.m
	1	+function RGBval = rtsHSVtoRGB(HSVval)
	2	+
	3	+H = HSVval(1);
	4	+S = HSVval(2);
	5	+V = HSVval(3);
	6	+
	7	+C = V*S;
	8	+Hprime = H/60;
	9	+X = C*(1 - abs(mod(Hprime, 2) - 1));
	10	+
	11	+RGBval = [0 0 0];
	12	+if Hprime >= 0 && Hprime < 1
	13	+ RGBval = [C X 0];
	14	+end
	15	+if Hprime >= 1 && Hprime < 2
	16	+ RGBval = [X C 0];
	17	+end
	18	+if Hprime >= 2 && Hprime < 3
	19	+ RGBval = [0 C X];
	20	+end
	21	+if Hprime >= 3 && Hprime < 4
	22	+ RGBval = [0 X C];
	23	+end
	24	+if Hprime >= 4 && Hprime < 5
	25	+ RGBval = [X 0 C];
	26	+end
	27	+if Hprime >= 5 && Hprime < 6
	28	+ RGBval = [C 0 X];
	29	+end
0	30	\ No newline at end of file
...	...

rtsInt2Str.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsInt2Str.m
	1	+function result = rtsInt2Str(num, digits)
	2	+
	3	+
	4	+%get the number of digits in the current number
	5	+temp = num;
	6	+current_digits = 0;
	7	+while temp >= 1
	8	+ temp = temp/10;
	9	+ current_digits = current_digits+1;
	10	+end
	11	+
	12	+if current_digits > digits
	13	+ result = int2str(num);
	14	+else
	15	+ result = [];
	16	+ for i = 1:digits - current_digits
	17	+ result = [result '0'];
	18	+ end
	19	+ result = [result int2str(num)];
	20	+end
0	21	\ No newline at end of file
...	...

rtsLoadImageStack.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsLoadImageStack.m
	1	+function R = rtsLoadImageStack(directoryName)
	2	+
	3	+fileList = dir(directoryName);
	4	+fileList = fileList(3:length(fileList));
	5	+
	6	+nFiles = length(fileList);
	7	+
	8	+%enable the progress bar
	9	+gui_active(1);
	10	+h = progressbar([], 0, ['Loading ' num2str(nFiles) ' Images...'], 'Load Stack');
	11	+
	12	+
	13	+%load each file into a volume
	14	+for i=1:nFiles
	15	+ fileName = [directoryName '\' fileList(i).name];
	16	+ image = imread(fileName);
	17	+
	18	+ %allocate space
	19	+ if i == 1
	20	+ R = zeros(size(image, 1), size(image, 2), nFiles, 'uint8');
	21	+ end
	22	+ R(:, :, i) = image(:, :, 1);
	23	+
	24	+ h = progressbar(h, 1/(nFiles));
	25	+ if ~gui_active
	26	+ progressbar(h, -1);
	27	+ break;
	28	+ end
	29	+end
	30	+
	31	+progressbar(h, -1);
	32	+
	33	+
...	...

rtsMaskReconstruct.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsMaskReconstruct.m
	1	+function I = rtsMaskReconstruct(data, mask, fillVal)
	2	+
	3	+%This function reconstructs an array of data created from a mask
	4	+% data = a VxD matrix, where D = positions and V = values
	5	+% mask = a binary image with D nonzero values
	6	+
	7	+%compute the number of values at each pixel
	8	+nV = size(data, 1);
	9	+
	10	+%allocate space
	11	+nX = size(mask, 1);
	12	+nY = size(mask, 2);
	13	+I = ones(nV, nXnY) fillVal;
	14	+
	15	+d=1;
	16	+%for each position
	17	+for p = 1:nX*nY
	18	+
	19	+ %if the mask value at the position is true, copy values
	20	+ if(mask(p) && d < size(data, 2))
	21	+ I(:, p) = data(:, d);
	22	+
	23	+ %increment the data pointer
	24	+ d = d + 1;
	25	+
	26	+ %if the mask value at the position is zero, set all values to zero
	27	+ else
	28	+ I(:, p) = ones(nV, 1) * fillVal;
	29	+ end
	30	+
	31	+end
	32	+
	33	+%reshape the array into an image
	34	+I = reshape(I', [nX nY nV]);
	35	+
...	...

rtsMonteCarloSphere.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsMonteCarloSphere.m
	1	+function S = rtsMonteCarloSphere(sqrtN, startAngle, endAngle)
	2	+
	3	+%allocate space for the samples
	4	+S = zeros(sqrtN^2, 2);
	5	+
	6	+if sqrtN == 1
	7	+ S = [0 0];
	8	+ return;
	9	+end
	10	+cosStart = cos(startAngle);
	11	+cosEnd = cos(endAngle);
	12	+
	13	+i=1; % array index
	14	+oneoverN = 1.0/sqrtN;
	15	+for a = 0:sqrtN-1
	16	+ for b = 0:sqrtN-1
	17	+ %generate unbiased distribution of spherical coords
	18	+ U1 = rand(1, 1);
	19	+ U2 = rand(1, 1);
	20	+ %x = 1.0 - (a + U1) * oneoverN * (1.0 - cosEnd)
	21	+ x = cosStart - (a + U1) * oneoverN * (1.0 - cosEnd);
	22	+ y = (b + U2) * oneoverN;
	23	+ theta = acos(x);
	24	+ phi = 2.0 * pi * y;
	25	+ S(i, 1) = theta;
	26	+ S(i, 2) = phi;
	27	+ i = i+1;
	28	+ end
	29	+end
...	...

rtsOffsetPlots.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsOffsetPlots.m
	1	+function R = rtsOffsetPlots(D, s)
	2	+%This function offsets a series of plots in D by intervals of s
	3	+%R = the resulting offset data
	4	+%D = the data as a 2D matrix
	5	+%s = the interval used as an offset
	6	+
	7	+ nPlots = size(D, 2);
	8	+
	9	+ offsetMat = repmat(1:nPlots, size(D, 1), 1) * s;
	10	+
	11	+ R = D + offsetMat;
0	12	\ No newline at end of file
...	...

rtsRAW2Images.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsRAW2Images.m
	1	+function rtsRAW2Images(source, destination_directory)
	2	+
	3	+for i=1:size(source, 3)
	4	+ filename = [destination_directory '\' rtsInt2Str(i, 4) '.bmp'];
	5	+ imwrite(source(:, :, i)', filename);
	6	+ disp(filename);
	7	+end
	8	+
0	9	\ No newline at end of file
...	...

rtsROC.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsROC.m
	1	+function [ROC, AUC, T] = rtsROC(thresholds, class)
	2	+
	3	+%Computes the ROC curves given a set of thresholds and class assignments
	4	+% class = boolean array describing actual class membership
	5	+% thresholds = threshold values for class assignment
	6	+% ROC = ROC curve as column vectors (1 = 1-specificity, 2 = sensitivity)
	7	+% AUC = area under the ROC curve computed using the trapezoid method
	8	+% T = threshold values associated with each 1-specificity value
	9	+
	10	+nValues = length(thresholds);
	11	+ROC = zeros(nValues, 2);
	12	+
	13	+%sort both arrays
	14	+[T, ix] = sort(thresholds, 'descend');
	15	+C = class(ix);
	16	+
	17	+%compute the necessary global values
	18	+P = nnz(class);
	19	+N = nValues - P;
	20	+
	21	+for i = 1:nValues
	22	+ TP = nnz(C(1:i));
	23	+ sensitivity = TP/P;
	24	+
	25	+ FP = i - TP;
	26	+ FPR = FP/N;
	27	+ specificity = 1 - FPR;
	28	+
	29	+ ROC(i, 1) = 1 - specificity;
	30	+ ROC(i, 2) = sensitivity;
	31	+end
	32	+
	33	+AUC = trapz(ROC(:, 1), ROC(:, 2));
	34	+
	35	+
0	36	\ No newline at end of file
...	...

rtsRandomizeMask.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsRandomizeMask.m
	1	+function R = rtsRandomizeMask(mask, N)
	2	+
	3	+%error checking
	4	+if N > nnz(mask)
	5	+ N = nnz(mask);
	6	+end
	7	+
	8	+if N < 0
	9	+ N = 0;
	10	+end
	11	+
	12	+%get the indices of all nonzero values in the mask
	13	+ind = find(mask);
	14	+
	15	+%randomize the indices
	16	+rind = ind(randperm(size(ind, 1)));
	17	+
	18	+%create the new randomized mask (random subset of the old mask)
	19	+R = zeros(size(mask));
	20	+R(rind(1:N)) = 1;
	21	+
	22	+R = R > 0;
0	23	\ No newline at end of file
...	...

rtsSaveRAW.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/rtsSaveRAW.m
	1	+function rtsSaveRAW(data, filename, type, endian)
	2	+%Loads a RAW image file.
	3	+%LoadRAW(filename, x, y, z, c, type) returns a [CxXxYxZ] array representing
	4	+%RAW data loaded from disk. 'type' defines the data type to be loaded (ex.
	5	+%'uint8'. 'endian' defines the byte order: 'b' = big endian, 'l' = little
	6	+%endian.
	7	+
	8	+file_id = fopen(filename, 'w', endian);
	9	+vol = zeros(size(data, 1), size(data, 2), size(data, 3), size(data, 4));
	10	+vol(:) = fwrite(file_id, data, type);
	11	+fclose(file_id);
0	12	\ No newline at end of file
...	...

stimBrewerColormap.m 0 → 100644

Wrap text Show/Hide comments View file @8be1ab9

	1	+++ a/stimBrewerColormap.m
	1	+function result = stimBrewerColormap(R)
	2	+
	3	+%returns a Brewer colormap with the specified resolution R
	4	+
	5	+ctrlPts = zeros(11, 3);
	6	+
	7	+ctrlPts(1, :) = [0.192157, 0.211765, 0.584314];
	8	+ctrlPts(2, :) = [0.270588, 0.458824, 0.705882];
	9	+ctrlPts(3, :) = [0.454902, 0.678431, 0.819608];
	10	+ctrlPts(4, :) = [0.670588, 0.85098, 0.913725];
	11	+ctrlPts(5, :) = [0.878431, 0.952941, 0.972549];
	12	+ctrlPts(6, :) = [1, 1, 0.74902];
	13	+ctrlPts(7, :) = [0.996078, 0.878431, 0.564706];
	14	+ctrlPts(8, :) = [0.992157, 0.682353, 0.380392];
	15	+ctrlPts(9, :) = [0.956863, 0.427451, 0.262745];
	16	+ctrlPts(10, :) = [0.843137, 0.188235, 0.152941];
	17	+ctrlPts(11, :) = [0.647059, 0, 0.14902];
	18	+
	19	+X = 1:11;
	20	+
	21	+r = 1:11/R:11;
	22	+
	23	+R = interp1(X, ctrlPts(:, 1), r);
	24	+G = interp1(X, ctrlPts(:, 2), r);
	25	+B = interp1(X, ctrlPts(:, 3), r);
	26	+
	27	+result = [R' G' B'];
	28	+
...	...