From 8be1ab930be4cf927f46b79c654a102b1c9bd8fc Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Tue, 7 Apr 2015 13:24:08 -0500 Subject: [PATCH] initial commit of the new Matlab repository --- rtsApplyBrewer.m | 41 +++++++++++++++++++++++++++++++++++++++++ rtsBaselineCorrect.m | 61 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtsClass2Color.m | 15 +++++++++++++++ rtsColorMapField.m | 30 ++++++++++++++++++++++++++++++ rtsColorMapRainbow.m | 37 +++++++++++++++++++++++++++++++++++++ rtsCropMask.m | 21 +++++++++++++++++++++ rtsEnviClassify.m | 39 +++++++++++++++++++++++++++++++++++++++ rtsEnviClose.m | 3 +++ rtsEnviID.m | 8 ++++++++ rtsEnviLoadClasses.m~ | 11 +++++++++++ rtsEnviLoadHeader.m | 78 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtsEnviLoadImage.m | 44 ++++++++++++++++++++++++++++++++++++++++++++ rtsEnviLoadTraining.m | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtsEnviOpen.m | 35 +++++++++++++++++++++++++++++++++++ rtsEnviRandomForest2C_apply.m | 41 +++++++++++++++++++++++++++++++++++++++++ rtsEnviRandomForest2C_train.m | 76 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtsEnviRandomForest2C_validate.m | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtsEnviRead.m | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtsEnviSaveImage.m | 41 +++++++++++++++++++++++++++++++++++++++++ rtsGaussianBlurND.m | 30 ++++++++++++++++++++++++++++++ rtsHSVtoRGB.m | 29 +++++++++++++++++++++++++++++ rtsInt2Str.m | 20 ++++++++++++++++++++ rtsLoadImageStack.m | 33 +++++++++++++++++++++++++++++++++ rtsMaskReconstruct.m | 35 +++++++++++++++++++++++++++++++++++ rtsMonteCarloSphere.m | 29 +++++++++++++++++++++++++++++ rtsOffsetPlots.m | 11 +++++++++++ rtsRAW2Images.m | 8 ++++++++ rtsROC.m | 35 +++++++++++++++++++++++++++++++++++ rtsRandomizeMask.m | 22 ++++++++++++++++++++++ rtsSaveRAW.m | 11 +++++++++++ stimBrewerColormap.m | 28 ++++++++++++++++++++++++++++ 31 files changed, 1062 insertions(+), 0 deletions(-) create mode 100644 rtsApplyBrewer.m create mode 100644 rtsBaselineCorrect.m create mode 100644 rtsClass2Color.m create mode 100644 rtsColorMapField.m create mode 100644 rtsColorMapRainbow.m create mode 100644 rtsCropMask.m create mode 100644 rtsEnviClassify.m create mode 100644 rtsEnviClose.m create mode 100644 rtsEnviID.m create mode 100644 rtsEnviLoadClasses.m~ create mode 100644 rtsEnviLoadHeader.m create mode 100644 rtsEnviLoadImage.m create mode 100644 rtsEnviLoadTraining.m create mode 100644 rtsEnviOpen.m create mode 100644 rtsEnviRandomForest2C_apply.m create mode 100644 rtsEnviRandomForest2C_train.m create mode 100644 rtsEnviRandomForest2C_validate.m create mode 100644 rtsEnviRead.m create mode 100644 rtsEnviSaveImage.m create mode 100644 rtsGaussianBlurND.m create mode 100644 rtsHSVtoRGB.m create mode 100644 rtsInt2Str.m create mode 100644 rtsLoadImageStack.m create mode 100644 rtsMaskReconstruct.m create mode 100644 rtsMonteCarloSphere.m create mode 100644 rtsOffsetPlots.m create mode 100644 rtsRAW2Images.m create mode 100644 rtsROC.m create mode 100644 rtsRandomizeMask.m create mode 100644 rtsSaveRAW.m create mode 100644 stimBrewerColormap.m diff --git a/rtsApplyBrewer.m b/rtsApplyBrewer.m new file mode 100644 index 0000000..ee140cc --- /dev/null +++ b/rtsApplyBrewer.m @@ -0,0 +1,41 @@ +function result = rtsApplyBrewer(field, minVal, maxVal) + +if nargin < 3 + maxVal = max(field(:)); +end +if nargin < 2 + minVal = min(field(:)); +end + +fieldMaskMin = field >= minVal; +fieldMaskMax = field <= maxVal; + +fieldMask = min(fieldMaskMin, fieldMaskMax); + +ctrlPts = zeros(11, 3); + +ctrlPts(1, :) = [0.192157, 0.211765, 0.584314]; +ctrlPts(2, :) = [0.270588, 0.458824, 0.705882]; +ctrlPts(3, :) = [0.454902, 0.678431, 0.819608]; +ctrlPts(4, :) = [0.670588, 0.85098, 0.913725]; +ctrlPts(5, :) = [0.878431, 0.952941, 0.972549]; +ctrlPts(6, :) = [1, 1, 0.74902]; +ctrlPts(7, :) = [0.996078, 0.878431, 0.564706]; +ctrlPts(8, :) = [0.992157, 0.682353, 0.380392]; +ctrlPts(9, :) = [0.956863, 0.427451, 0.262745]; +ctrlPts(10, :) = [0.843137, 0.188235, 0.152941]; +ctrlPts(11, :) = [0.647059, 0, 0.14902]; + +%X = 0:1/10:1; +X = minVal:(maxVal - minVal)/10:maxVal; + +R = interp1(X, ctrlPts(:, 1), field); +G = interp1(X, ctrlPts(:, 2), field); +B = interp1(X, ctrlPts(:, 3), field); + +R = R.* fieldMask; +G = G.* fieldMask; +B = B.* fieldMask; + +result = cat(ndims(field)+1, R, G, B); + diff --git a/rtsBaselineCorrect.m b/rtsBaselineCorrect.m new file mode 100644 index 0000000..43ba602 --- /dev/null +++ b/rtsBaselineCorrect.m @@ -0,0 +1,61 @@ +function R = rtsBaselineCorrect(A, b, dim) + +%This function baselines an n-dimensional array given a series of points b +%B = baseline-corrected data +%A = raw data +%b = baseline points in the format [b0 b1 ... bn] +%dim = dimension along which the baseline correction will occur + + +%allocate space for the sub-component indices +sz = size(A); +inds = repmat({1},1,ndims(A)); + +%allocate space for the result +R = zeros(size(A)); + +%for each baseline point +for bp = 1:length(b) - 1 + + %extract the data at the first baseline point + n = b(bp); + inds = getIndices(A, n, dim); + BP1 = A(inds{:}); + + %extract the data at the second baseline point + n = b(bp+1); + inds = getIndices(A, n, dim); + BP2 = A(inds{:}); + + %for each point in between the baseline points + for p = b(bp):b(bp+1) + + %compute the weighting of the linear function between BP1 and BP2 + a2 = (p - b(bp))/(b(bp+1) - b(bp)); + a1 = 1.0 - a2; + + n = p; + inds = getIndices(A, n, dim); + + baseFunc = a1 * BP1 + a2 * BP2; + R(inds{:}) = A(inds{:}) - baseFunc; + end +end + + +function I = getIndices(A, n, dim) + +%this function returns indices for all elements in S at +%point n along dimension dim + +sz = size(A); +I = repmat({1},1,ndims(A)); + +for d = 1:ndims(A) + if d ~= dim + I{d} = 1:sz(d); + else + I{d} = n; + end +end + diff --git a/rtsClass2Color.m b/rtsClass2Color.m new file mode 100644 index 0000000..7377873 --- /dev/null +++ b/rtsClass2Color.m @@ -0,0 +1,15 @@ +function C = rtsClass2Color(classes, classList, colorList) + +C = zeros(length(classes), 3); + +for i = 1:length(classes) + + for c = 1:length(classList) + + if classes(i) == classList(c) + C(i, :) = colorList(c, :); + end + + end + +end \ No newline at end of file diff --git a/rtsColorMapField.m b/rtsColorMapField.m new file mode 100644 index 0000000..d947377 --- /dev/null +++ b/rtsColorMapField.m @@ -0,0 +1,30 @@ +function Result = rtsColorMapField(F, range) +%creates a color map for the field F and returns an array with +%an additional size-3 color dimension + +R = zeros(size(F)); +G = zeros(size(F)); +B = zeros(size(F)); + +%if no range is specified, estimate it using min and max values +if nargin == 1 + range(1) = min(F(:)); + range(2) = max(F(:)); +end +range_size = range(2) - range(1); + +num_elements = size(F(:)); +for i =1:num_elements + if(F(i) > range(1) && F(i) < range(2)) + hue_degrees = (1 - (F(i) - range(1))/range_size)*240; + HSVval = [hue_degrees 1.0 1.0]; + RGBval = rtsHSVtoRGB(HSVval); + else + RGBval = [0 0 F(i)/range(1)]; + end + R(i) = RGBval(1); + G(i) = RGBval(2); + B(i) = RGBval(3); +end +Result = cat(ndims(F)+1, R, G, B); +imshow(Result); diff --git a/rtsColorMapRainbow.m b/rtsColorMapRainbow.m new file mode 100644 index 0000000..ff38316 --- /dev/null +++ b/rtsColorMapRainbow.m @@ -0,0 +1,37 @@ +function result = rtsColorMapRainbow(field, minVal, maxVal) + +if nargin < 3 + maxVal = max(field(:)); +end +if nargin < 2 + minVal = min(field(:)); +end + +fieldMaskMin = field >= minVal; +fieldMaskMax = field <= maxVal; + +fieldMask = fieldMaskMin;%min(fieldMaskMin, fieldMaskMax); + +ctrlPts = zeros(3, 3); + +ctrlPts(1, :) = [0.0, 0.0, 0.5]; +ctrlPts(2, :) = [0.0, 0.0, 1.0]; +ctrlPts(3, :) = [0.0, 1.0, 1.0]; +ctrlPts(4, :) = [0.0, 1.0, 0.0]; +ctrlPts(5, :) = [1.0, 1.0, 0.0]; +ctrlPts(6, :) = [1.0, 0.0, 0.0]; +ctrlPts(7, :) = [0.5, 0.0, 0.0]; + +%X = 0:1/10:1; +X = minVal:(maxVal - minVal)/6:maxVal; + +R = interp1(X, ctrlPts(:, 1), field); +G = interp1(X, ctrlPts(:, 2), field); +B = interp1(X, ctrlPts(:, 3), field); + +%R = R.* fieldMask; +%G = G.* fieldMask; +%B = B.* fieldMask; + +result = cat(ndims(field)+1, R, G, B); + diff --git a/rtsCropMask.m b/rtsCropMask.m new file mode 100644 index 0000000..0c36960 --- /dev/null +++ b/rtsCropMask.m @@ -0,0 +1,21 @@ +function new_mask = rtsCropMask(mask, max_pixels) + +%convert to logical +%mask = mask > 0; + +if nnz(mask) < max_pixels + new_mask = mask; +else + j = 0; + for i = 1:length(mask(:)) + if mask(i) > 0 + j = j + 1; + if j == max_pixels + last_pixel = i; + end + end + end + new_mask = zeros(size(mask, 1), size(mask, 2)); + new_mask(1:last_pixel) = mask(1:last_pixel); +end + \ No newline at end of file diff --git a/rtsEnviClassify.m b/rtsEnviClassify.m new file mode 100644 index 0000000..3c3dfb3 --- /dev/null +++ b/rtsEnviClassify.m @@ -0,0 +1,39 @@ +function C = rtsEnviClassify(B, envi_filename, mask, nFeatures) + +% B = TreeBagger class +% nFeatures = number of features cropped out (0 means use all features) + +if nargin < 4 + nFeatures = 0; +end + + +batchSize = 100000; + + + +Ns = nnz(mask) %get the number of samples to be classified +Nb = ceil(Ns / batchSize); %get the number of batches + +C(Ns, 1) = char(0); + +%open the ENVI file +envi = rtsEnviOpen(envi_filename, [envi_filename '.hdr'], mask); + +for b = 1:Nb + + %load a batch + batch = rtsEnviRead(envi, batchSize); + + %crop out the number of features specified (if it is specified) + if nFeatures ~= 0 + batch = batch(1:nFeatures, :); + end + + %classify the batch + C( (b-1)*batchSize + 1 : (b-1)*batchSize + size(batch, 2) ) = cell2mat(predict(B, batch')); + + disp([num2str( round( b / Nb * 100 ) ) '%']); +end + +rtsEnviClose(envi); \ No newline at end of file diff --git a/rtsEnviClose.m b/rtsEnviClose.m new file mode 100644 index 0000000..3ee3283 --- /dev/null +++ b/rtsEnviClose.m @@ -0,0 +1,3 @@ +function rtsEnviClose(enviFID) + +fclose(enviFID.fid); \ No newline at end of file diff --git a/rtsEnviID.m b/rtsEnviID.m new file mode 100644 index 0000000..10f559e --- /dev/null +++ b/rtsEnviID.m @@ -0,0 +1,8 @@ +classdef rtsEnviID < handle + properties + fid + header + mask + fpos + end +end \ No newline at end of file diff --git a/rtsEnviLoadClasses.m~ b/rtsEnviLoadClasses.m~ new file mode 100644 index 0000000..ee30933 --- /dev/null +++ b/rtsEnviLoadClasses.m~ @@ -0,0 +1,11 @@ +function [F, C] = rtsEnviLoadClasses(filenames, classnames, ppc) + +Nc = size(filenames, 1); + +for c = 1:Nc + %load the masks + image = imread(filenames{c}); + mask = image(:, :, 1)' > 0; + + +end \ No newline at end of file diff --git a/rtsEnviLoadHeader.m b/rtsEnviLoadHeader.m new file mode 100644 index 0000000..839a9e4 --- /dev/null +++ b/rtsEnviLoadHeader.m @@ -0,0 +1,78 @@ +function s = rtsEnviLoadHeader(filename) +%create a cell array of fields +s = struct; + +fid = fopen(filename, 'r', 'n', 'US-ASCII'); + +%read the first field and make sure it is "ENVI" +fname = GetFieldName(fid); +if strcmp(fname, 'ENVI') == 0 + disp('Not an ENVI header file'); + return; +end + +while feof(fid) == 0 + fname = GetFieldName(fid); + if feof(fid) == 1 + return; + end + [value, valid] = ReadField(fid, fname); + if valid == 1 + s = setfield(s, fname, value); + end +end +fclose(fid); + +function t = GetFieldName(fid) +string_struct = textscan(fid, '%s', 1, 'Delimiter', '='); +if feof(fid) == 1 + t = []; + return; +end +t = string_struct{1}{1}; +t = strtrim(t); +t(t==' ') = '_'; + +function [v, valid] = ReadField(fid, field_name) +valid = 1; +stringFields = {'file_type', 'interleave', 'sensor_type', 'wavelength_units'}; +intFields = {'samples', 'lines', 'bands', 'header_offset', 'data_type', 'byte_order'}; + +%if the field is "description", read between the brackets +if strcmp(field_name, 'description') == 1 + textscan(fid, '%[{]', 1); + string_struct = textscan(fid, '%[^}]', 1, 'Whitespace', ''); + textscan(fid, '%[}]', 1); + v = string_struct{1}{1}; + v = strtrim(v); + return; +end +if max(strcmp(field_name, intFields)) ~= 0 + v = fscanf(fid, '%d'); + return; +end +if max(strcmp(field_name, stringFields)) ~= 0 + string_struct = textscan(fid, '%s', 1, 'Whitespace', '\n'); + v = string_struct{1}{1}; + v = strtrim(v); + return; +end + +%read and return the wavelength values +if strcmp(field_name, 'wavelength') == 1 + v = []; + textscan(fid, '%[{]', 1); + c = ' '; + while c ~= '}' + new = fscanf(fid, '%f'); + v = [v new]; + c = fscanf(fid, '%c', 1); + end + return; +end + +%if it doesn't match anything, just read until the end of the line +%string_struct = textscan(fid, '%s', 1, 'Whitespace', '\n'); +string_struct = textscan(fid, '%s', 1, 'Delimiter', '\n'); +v = ''; +valid = 0; diff --git a/rtsEnviLoadImage.m b/rtsEnviLoadImage.m new file mode 100644 index 0000000..dd6510d --- /dev/null +++ b/rtsEnviLoadImage.m @@ -0,0 +1,44 @@ +function [image, header] = rtsEnviLoadImage(filename, headername) +%enter the file name (no .hdr extension) + +header_file = headername; +header = rtsEnviLoadHeader(header_file); + +if strcmp(header.interleave, 'bsq') == 1 + image = zeros(header.samples, header.lines, header.bands); +end +if strcmp(header.interleave, 'bip') == 1 + image = zeros(header.bands, header.samples, header.lines); +end +if strcmp(header.interleave, 'bil') == 1 + image = zeros(header.samples, header.bands, header.lines); +end + + +file_bytes = header.samples*header.lines*header.bands; +%check to make sure there is enough memory available +if ispc + [~, sys] = memory; + if sys.PhysicalMemory.Available < file_bytes*header.data_type && strcmp(header.interleave, 'bsq') == 0 + strResponse = input('The data set will require virtual memory for permutation. Continue? (y/n) ', 's'); + if strcmp(strResponse, 'n') == 1 + return; + end + end +else + disp('Sorry, I can''t check available memory for this OS. Keep your fingers crossed!'); +end + +fid = fopen(filename); +%skip the header - start reading header_offset bytes from the beginning +fseek(fid, header.header_offset, 'bof'); + +%read the data +image(:) = fread(fid, file_bytes, 'float32'); + +if strcmp(header.interleave, 'bip') == 1 + image = permute(image, [2 3 1]); +end +if strcmp(header.interleave, 'bil') == 1 + image = permute(image, [1 3 2]); +end \ No newline at end of file diff --git a/rtsEnviLoadTraining.m b/rtsEnviLoadTraining.m new file mode 100644 index 0000000..1788a08 --- /dev/null +++ b/rtsEnviLoadTraining.m @@ -0,0 +1,57 @@ +function [F, C] = rtsEnviLoadTraining(filenames, classnames, binfile, ppc) + +Nc = length(filenames); %calculate the number of classes + +if nargin < 3 + ppc = 0; +end + +%---------Load the Masks----------- +%first, get the total number of pixels +Ns = 0; +for c = 1:Nc + %load the mask + image = imread(filenames{c}); + + if c == 1 + mask = boolean(zeros(size(image, 2), size(image, 1), Nc)); + end + + if ppc ~= 0 + mask(:, :, c) = rtsRandomizeMask(image(:, :, 1)' > 0, ppc); + else + mask(:, :, c) = image(:, :, 1)' > 0; + end + +end + +%calculate the number of samples +Ns = Ns + nnz(mask); + +%-----------Load the Training Data------------- +disp(['loading ' num2str(Ns) ' samples']); +ci = 1; +for c = 1:Nc + + disp(['loading class ' classnames(c) '...']); + batch = nnz(mask(:, :, c)); + %create an ENVI file object + envi = rtsEnviOpen(binfile, [binfile '.hdr'], mask(:, :, c)); + + if c == 1 + F = zeros(Ns, envi.header.bands); + C(Ns, 1) = char(0); + end + + %read a bunch of spectra + F(ci:ci+batch - 1, :) = rtsEnviRead(envi, batch)'; + C(ci:ci+batch - 1, 1) = repmat([classnames(c)], [batch 1]); + + %update the current index + ci = ci+batch; + + %close the ENVI file + rtsEnviClose(envi); + + disp('done'); +end \ No newline at end of file diff --git a/rtsEnviOpen.m b/rtsEnviOpen.m new file mode 100644 index 0000000..badc640 --- /dev/null +++ b/rtsEnviOpen.m @@ -0,0 +1,35 @@ +function EID = rtsEnviOpen(filename, headername, mask) + +%opens an ENVI file and returns an ID structure +% filename = name of ENVI file +% headername = name of ENVI header +% FID = structure containing file information +% mask = binary image specifying valid spectra + +%assign the mask variable +if(nargin < 3) + mask = 1; +end + +%open the file and save the file ID +fid = fopen(filename, 'r'); + +%open the ENVI header file +header = rtsEnviLoadHeader(headername); + +%this is currently only valid for BIP files +if(~strcmp(header.interleave, 'bip')) + disp('Error: This function only works for BIP interleave files. Load in ENVI and convert.'); +end +if(header.data_type ~= 4) + disp('Error: This function only works for floating-point files.'); +end + +EID = rtsEnviID; +EID.fid = fid; +EID.header = header; +EID.mask = mask; +EID.fpos = 0; + +%EID = struct('fid', fid, 'header', header, 'mask', mask, 'fpos', 0); + diff --git a/rtsEnviRandomForest2C_apply.m b/rtsEnviRandomForest2C_apply.m new file mode 100644 index 0000000..e1d2a7c --- /dev/null +++ b/rtsEnviRandomForest2C_apply.m @@ -0,0 +1,41 @@ +function P = rtsEnviRandomForest2C_apply(RF, EnviFileName, EnviHeaderName, TissueMask, reference) + +%Applies a 2-class random forest classifier to the given image. Only +%pixels specified in TissueMask are classified. The estimated rule data is +%returned as an image sequence. + +%load the tissue mask +maskImage = imread(TissueMask); +maskBinary = (maskImage(:, :, 1) > 0)'; + +nPixels = nnz(maskBinary); +disp(['Number of pixels: ' num2str(nPixels)]); + + +%load the data +disp(['Loading Features: ' EnviFileName]); +tLoadTime = tic; +fid = rtsEnviOpen(EnviFileName, EnviHeaderName, maskBinary); +F = rtsEnviRead(fid, nPixels); +rtsEnviClose(fid); +disp(['Time: ' num2str(toc(tLoadTime)) 's']); + +%transpose +F = F'; + +%apply the reference +if nargin == 5 + Fnorm = repmat(F(:, reference), 1, size(F, 2)); + F = F./Fnorm; + F(:, reference) = []; +end + +%classify the data +disp('Classifying...'); +tClassTime = tic; +T = predict(RF, F); +disp(['Time: ' num2str(toc(tClassTime)) 's']); +%Cl = Cpost > threshold; + +%reconstruct the image from the mask +P = rtsMaskReconstruct(T', maskBinary, -1); \ No newline at end of file diff --git a/rtsEnviRandomForest2C_train.m b/rtsEnviRandomForest2C_train.m new file mode 100644 index 0000000..8bf13d7 --- /dev/null +++ b/rtsEnviRandomForest2C_train.m @@ -0,0 +1,76 @@ +function RF = rtsEnviRandomForest2C_train(EnviFileName, EnviHeaderName, RoiImages, reference) + +%Creates a 2-class random forest classifier for an Envi BIP file using the +%provided ROI images. The resulting classifier uses regression to map +%class A to 1 and class B to 0. This allows the creation of ROC curves. +% +% EnviFilename - Name of an ENVI BIP file +% RoiImages - Cell array containing names of ROI images (masks) + +%default parameters +maxPixels = 200000; +threshold = 0.5; +nTrees = 100; +nThreads = 8; +%trainPixels = 200000; + +%determine the number of classes +nClasses = length(RoiImages); +%make sure that there are only two classes +if nClasses > 2 + disp('This classifier only supports 2 classes.'); + return; +end + +%for each class, load the training data +T = []; +F = []; +for c = 1:nClasses + + %load the class mask + maskImage = imread(RoiImages{c}); + maskBinary = (maskImage(:, :, 1) > 0)'; + + disp('------------------------------------------------------'); + %load epithelium spectra + disp(['Loading Training Class ' num2str(c) ' pixels: ' EnviFileName]); + tLoadTime = tic; + fid = rtsEnviOpen(EnviFileName, EnviHeaderName, maskBinary); + Fc = rtsEnviRead(fid, maxPixels); + rtsEnviClose(fid); + + if c == 1 + Tc = ones(size(Fc, 2), 1); + else + Tc = zeros(size(Fc, 2), 1); + end + + disp(['Time: ' num2str(toc(tLoadTime)) 's']); + + %add features and targets to the final vectors + T = [T; Tc]; + F = [F; Fc']; +end + +%apply the reference +if nargin == 4 + Fnorm = repmat(F(:, reference), 1, size(F, 2)); + F = F./Fnorm; + F(:, reference) = []; +end + +%parallelize if specified +if nThreads > 1 + matlabpool('open', nThreads); + paraoptions = statset('UseParallel', 'always'); +end + +%train the classifier +disp('Creating a Random Forest classifier...'); +tTrainTime = tic; +RF = TreeBagger(nTrees, F, T, 'method', 'regression', 'Options',paraoptions); +disp(['Time: ' num2str(toc(tTrainTime)) 's']); + +if nThreads > 1 + matlabpool close +end \ No newline at end of file diff --git a/rtsEnviRandomForest2C_validate.m b/rtsEnviRandomForest2C_validate.m new file mode 100644 index 0000000..626a03c --- /dev/null +++ b/rtsEnviRandomForest2C_validate.m @@ -0,0 +1,65 @@ +function [Croc, Cauc] = rtsEnviRandomForest2C_validate(RF, EnviFileName, EnviHeaderName, RoiImages, reference) + +%Validates a 2-class random forest classifier, creating ROC curves. + +%parameters +maxPixels = 200000; + +%determine the number of classes +nClasses = length(RoiImages); +%make sure that there are only two classes +if nClasses > 2 + disp('This classifier only supports 2 classes.'); + return; +end + +%for each class, load the validation data +T = []; +F = []; +for c = 1:nClasses + + %load the class mask + maskImage = imread(RoiImages{c}); + maskBinary = (maskImage(:, :, 1) > 0)'; + + disp('------------------------------------------------------'); + %load epithelium spectra + disp(['Loading Validation Class ' num2str(c) ' pixels: ' EnviFileName]); + tLoadTime = tic; + fid = rtsEnviOpen(EnviFileName, EnviHeaderName, maskBinary); + Fc = rtsEnviRead(fid, maxPixels); + rtsEnviClose(fid); + + if c == 1 + Tc = ones(size(Fc, 2), 1); + else + Tc = zeros(size(Fc, 2), 1); + end + + disp(['Time: ' num2str(toc(tLoadTime)) 's']); + + %add features and targets to the final vectors + T = [T; Tc]; + F = [F; Fc']; +end + +%apply the reference +if nargin == 5 + Fnorm = repmat(F(:, reference), 1, size(F, 2)); + F = F./Fnorm; + F(:, reference) = []; +end + +%classify the data (get the estimated posterior probability) +disp('Validating...'); +tClassTime = tic; +Tpost = predict(RF, F); +disp(['Time: ' num2str(toc(tClassTime)) 's']); + +%calculate and display the ROC curve +disp('Calculating ROC...'); +tRocTime = tic; +[Croc, Cauc, threshList] = rtsROC(Tpost, T); +disp(['Time: ' num2str(toc(tRocTime)) 's']); +disp(['AUC = ' num2str(Cauc)]); +plot(Croc(:, 1), Croc(:, 2)); \ No newline at end of file diff --git a/rtsEnviRead.m b/rtsEnviRead.m new file mode 100644 index 0000000..d46db6c --- /dev/null +++ b/rtsEnviRead.m @@ -0,0 +1,68 @@ +function S = rtsEnviRead(fid, batchSize) + +%This function reads a batch of spectra from the given ENVI file ID +% fid = ENVI file ID created using rtsEnviOpen +% batchSize = number of spectra to read + +%if there is no mask, just load each spectrum in order +if(fid.mask == 1) + + %compute the new batch size in case we are near the eof + nSpectra = fid.header.samples * fid.header.lines; + remainingSpectra = nSpectra - fid.fpos; + + if(batchSize > remainingSpectra) + batchSize = remainingSpectra; + end + + S = fread(fid.fid, [fid.header.bands batchSize], '*float32'); + + %increment the fid counter + fid.fpos = fid.fpos + batchSize; + +%otherwise only load valid spectra +else + + %compute the new batch size in case we are near the eof + if(fid.fpos == 0) + remainingSpectra = nnz(fid.mask); + else + nSpectra = nnz(fid.mask); + maskRead = fid.mask(1:fid.fpos+1); + remainingSpectra = nSpectra - nnz(maskRead); + end + + if(batchSize > remainingSpectra) + batchSize = remainingSpectra; + end + + %allocate space for the spectra + S = zeros(fid.header.bands, batchSize); + + %for each spectrum in the batch + for s = 1:batchSize + + %while the current spectrum is invalid + skip = 0; + while (~fid.mask(fid.fpos + 1)) + %read the invalid spectrum + %invalid = fread(fid.fid, [fid.header.bands 1], '*float32'); + skip = skip + 1; + + %increment the file position + fid.fpos = fid.fpos + 1; + end + fseek(fid.fid, skip * fid.header.bands * 4, 'cof'); + + test = fread(fid.fid, [fid.header.bands 1], '*float32'); + if size(test) ~= size(S(:, s)) + size(test) + size(S(:, s)) + end + S(:, s) = test; + fid.fpos = fid.fpos + 1; + + end + +end + diff --git a/rtsEnviSaveImage.m b/rtsEnviSaveImage.m new file mode 100644 index 0000000..a22402b --- /dev/null +++ b/rtsEnviSaveImage.m @@ -0,0 +1,41 @@ +function rtsEnviSaveImage(S, fileName, headerName, W) + +% S = 3D hyperspectral data +% W = list of wavelength values + +%save the image +rtsSaveRAW(S, fileName, 'float32', 'l'); + +%create the header file +headerFile = [fileName '.hdr']; +if nargin == 3 + headerFile = headerName; +end + +fid = fopen(headerFile, 'w'); + +%write the header and description +fprintf(fid, 'ENVI\n'); +fprintf(fid, 'description = {\n'); +fprintf(fid, ' File Saved using Matlab.}\n'); +fprintf(fid, 'samples = %d\n', size(S, 1)); +fprintf(fid, 'lines = %d\n', size(S, 2)); +fprintf(fid, 'bands = %d\n', size(S, 3)); +fprintf(fid, 'header offset = 0\n'); +fprintf(fid, 'file type = ENVI Standard\n'); +fprintf(fid, 'data type = 4\n'); +fprintf(fid, 'interleave = bsq\n'); +fprintf(fid, 'sensor type = Unknown\n'); +fprintf(fid, 'byte order = 0\n'); + +if nargin == 4 + %output wavelength units + fprintf(fid, 'wavelength = {\n'); + for i = 1:length(W)-1 + fprintf(fid, ' %f,\n', W(i)); + end + fprintf(fid, ' %f\n', W(length(W))); + fprintf(fid, '}\n'); +end + +fclose(fid); diff --git a/rtsGaussianBlurND.m b/rtsGaussianBlurND.m new file mode 100644 index 0000000..549d15c --- /dev/null +++ b/rtsGaussianBlurND.m @@ -0,0 +1,30 @@ +function R = rtsGaussianBlurND(input, sigma) +%This function convolves an ND array with a gaussian kernel specified by +%sigma. This takes advantage of the separability of the kernel and is +%therefore very fast. + +%disp('begin***********************************************'); + +%get the number of dimensions for the convolution +dim = ndims(input); + +%initialize the result to the source +R = input; + + +for d=1:dim + %if the dimension is greater than 1 and the sigma not a dirac + if size(input, d) > 1 && sigma(d) > 0 + %create a 1D kernel + kernelshape = ones(1, dim); + kernelshape(d) = sigma(d)*6; + kernel = zeros(kernelshape); + kernel(:) = ndgauss(sigma(d)*6, sigma(d)); + + %perform the convolution + R = convn(R, kernel, 'same'); + end +end + + + \ No newline at end of file diff --git a/rtsHSVtoRGB.m b/rtsHSVtoRGB.m new file mode 100644 index 0000000..6524ea9 --- /dev/null +++ b/rtsHSVtoRGB.m @@ -0,0 +1,29 @@ +function RGBval = rtsHSVtoRGB(HSVval) + +H = HSVval(1); +S = HSVval(2); +V = HSVval(3); + +C = V*S; +Hprime = H/60; +X = C*(1 - abs(mod(Hprime, 2) - 1)); + +RGBval = [0 0 0]; +if Hprime >= 0 && Hprime < 1 + RGBval = [C X 0]; +end +if Hprime >= 1 && Hprime < 2 + RGBval = [X C 0]; +end +if Hprime >= 2 && Hprime < 3 + RGBval = [0 C X]; +end +if Hprime >= 3 && Hprime < 4 + RGBval = [0 X C]; +end +if Hprime >= 4 && Hprime < 5 + RGBval = [X 0 C]; +end +if Hprime >= 5 && Hprime < 6 + RGBval = [C 0 X]; +end \ No newline at end of file diff --git a/rtsInt2Str.m b/rtsInt2Str.m new file mode 100644 index 0000000..ac1d025 --- /dev/null +++ b/rtsInt2Str.m @@ -0,0 +1,20 @@ +function result = rtsInt2Str(num, digits) + + +%get the number of digits in the current number +temp = num; +current_digits = 0; +while temp >= 1 + temp = temp/10; + current_digits = current_digits+1; +end + +if current_digits > digits + result = int2str(num); +else + result = []; + for i = 1:digits - current_digits + result = [result '0']; + end + result = [result int2str(num)]; +end \ No newline at end of file diff --git a/rtsLoadImageStack.m b/rtsLoadImageStack.m new file mode 100644 index 0000000..9a8c3b9 --- /dev/null +++ b/rtsLoadImageStack.m @@ -0,0 +1,33 @@ +function R = rtsLoadImageStack(directoryName) + +fileList = dir(directoryName); +fileList = fileList(3:length(fileList)); + +nFiles = length(fileList); + +%enable the progress bar +gui_active(1); +h = progressbar([], 0, ['Loading ' num2str(nFiles) ' Images...'], 'Load Stack'); + + +%load each file into a volume +for i=1:nFiles + fileName = [directoryName '\' fileList(i).name]; + image = imread(fileName); + + %allocate space + if i == 1 + R = zeros(size(image, 1), size(image, 2), nFiles, 'uint8'); + end + R(:, :, i) = image(:, :, 1); + + h = progressbar(h, 1/(nFiles)); + if ~gui_active + progressbar(h, -1); + break; + end +end + +progressbar(h, -1); + + diff --git a/rtsMaskReconstruct.m b/rtsMaskReconstruct.m new file mode 100644 index 0000000..eed102e --- /dev/null +++ b/rtsMaskReconstruct.m @@ -0,0 +1,35 @@ +function I = rtsMaskReconstruct(data, mask, fillVal) + +%This function reconstructs an array of data created from a mask +% data = a VxD matrix, where D = positions and V = values +% mask = a binary image with D nonzero values + +%compute the number of values at each pixel +nV = size(data, 1); + +%allocate space +nX = size(mask, 1); +nY = size(mask, 2); +I = ones(nV, nX*nY) * fillVal; + +d=1; +%for each position +for p = 1:nX*nY + + %if the mask value at the position is true, copy values + if(mask(p) && d < size(data, 2)) + I(:, p) = data(:, d); + + %increment the data pointer + d = d + 1; + + %if the mask value at the position is zero, set all values to zero + else + I(:, p) = ones(nV, 1) * fillVal; + end + +end + +%reshape the array into an image +I = reshape(I', [nX nY nV]); + diff --git a/rtsMonteCarloSphere.m b/rtsMonteCarloSphere.m new file mode 100644 index 0000000..c749360 --- /dev/null +++ b/rtsMonteCarloSphere.m @@ -0,0 +1,29 @@ +function S = rtsMonteCarloSphere(sqrtN, startAngle, endAngle) + +%allocate space for the samples +S = zeros(sqrtN^2, 2); + +if sqrtN == 1 + S = [0 0]; + return; +end +cosStart = cos(startAngle); +cosEnd = cos(endAngle); + +i=1; % array index +oneoverN = 1.0/sqrtN; +for a = 0:sqrtN-1 + for b = 0:sqrtN-1 + %generate unbiased distribution of spherical coords + U1 = rand(1, 1); + U2 = rand(1, 1); + %x = 1.0 - (a + U1) * oneoverN * (1.0 - cosEnd) + x = cosStart - (a + U1) * oneoverN * (1.0 - cosEnd); + y = (b + U2) * oneoverN; + theta = acos(x); + phi = 2.0 * pi * y; + S(i, 1) = theta; + S(i, 2) = phi; + i = i+1; + end +end diff --git a/rtsOffsetPlots.m b/rtsOffsetPlots.m new file mode 100644 index 0000000..eebc22e --- /dev/null +++ b/rtsOffsetPlots.m @@ -0,0 +1,11 @@ +function R = rtsOffsetPlots(D, s) +%This function offsets a series of plots in D by intervals of s +%R = the resulting offset data +%D = the data as a 2D matrix +%s = the interval used as an offset + + nPlots = size(D, 2); + + offsetMat = repmat(1:nPlots, size(D, 1), 1) * s; + + R = D + offsetMat; \ No newline at end of file diff --git a/rtsRAW2Images.m b/rtsRAW2Images.m new file mode 100644 index 0000000..334f360 --- /dev/null +++ b/rtsRAW2Images.m @@ -0,0 +1,8 @@ +function rtsRAW2Images(source, destination_directory) + +for i=1:size(source, 3) + filename = [destination_directory '\' rtsInt2Str(i, 4) '.bmp']; + imwrite(source(:, :, i)', filename); + disp(filename); +end + \ No newline at end of file diff --git a/rtsROC.m b/rtsROC.m new file mode 100644 index 0000000..82721b1 --- /dev/null +++ b/rtsROC.m @@ -0,0 +1,35 @@ +function [ROC, AUC, T] = rtsROC(thresholds, class) + +%Computes the ROC curves given a set of thresholds and class assignments +% class = boolean array describing actual class membership +% thresholds = threshold values for class assignment +% ROC = ROC curve as column vectors (1 = 1-specificity, 2 = sensitivity) +% AUC = area under the ROC curve computed using the trapezoid method +% T = threshold values associated with each 1-specificity value + +nValues = length(thresholds); +ROC = zeros(nValues, 2); + +%sort both arrays +[T, ix] = sort(thresholds, 'descend'); +C = class(ix); + +%compute the necessary global values +P = nnz(class); +N = nValues - P; + +for i = 1:nValues + TP = nnz(C(1:i)); + sensitivity = TP/P; + + FP = i - TP; + FPR = FP/N; + specificity = 1 - FPR; + + ROC(i, 1) = 1 - specificity; + ROC(i, 2) = sensitivity; +end + +AUC = trapz(ROC(:, 1), ROC(:, 2)); + + \ No newline at end of file diff --git a/rtsRandomizeMask.m b/rtsRandomizeMask.m new file mode 100644 index 0000000..62866d3 --- /dev/null +++ b/rtsRandomizeMask.m @@ -0,0 +1,22 @@ +function R = rtsRandomizeMask(mask, N) + +%error checking +if N > nnz(mask) + N = nnz(mask); +end + +if N < 0 + N = 0; +end + +%get the indices of all nonzero values in the mask +ind = find(mask); + +%randomize the indices +rind = ind(randperm(size(ind, 1))); + +%create the new randomized mask (random subset of the old mask) +R = zeros(size(mask)); +R(rind(1:N)) = 1; + +R = R > 0; \ No newline at end of file diff --git a/rtsSaveRAW.m b/rtsSaveRAW.m new file mode 100644 index 0000000..a6745c1 --- /dev/null +++ b/rtsSaveRAW.m @@ -0,0 +1,11 @@ +function rtsSaveRAW(data, filename, type, endian) +%Loads a RAW image file. +%LoadRAW(filename, x, y, z, c, type) returns a [CxXxYxZ] array representing +%RAW data loaded from disk. 'type' defines the data type to be loaded (ex. +%'uint8'. 'endian' defines the byte order: 'b' = big endian, 'l' = little +%endian. + +file_id = fopen(filename, 'w', endian); +vol = zeros(size(data, 1), size(data, 2), size(data, 3), size(data, 4)); +vol(:) = fwrite(file_id, data, type); +fclose(file_id); \ No newline at end of file diff --git a/stimBrewerColormap.m b/stimBrewerColormap.m new file mode 100644 index 0000000..0a2c521 --- /dev/null +++ b/stimBrewerColormap.m @@ -0,0 +1,28 @@ +function result = stimBrewerColormap(R) + +%returns a Brewer colormap with the specified resolution R + +ctrlPts = zeros(11, 3); + +ctrlPts(1, :) = [0.192157, 0.211765, 0.584314]; +ctrlPts(2, :) = [0.270588, 0.458824, 0.705882]; +ctrlPts(3, :) = [0.454902, 0.678431, 0.819608]; +ctrlPts(4, :) = [0.670588, 0.85098, 0.913725]; +ctrlPts(5, :) = [0.878431, 0.952941, 0.972549]; +ctrlPts(6, :) = [1, 1, 0.74902]; +ctrlPts(7, :) = [0.996078, 0.878431, 0.564706]; +ctrlPts(8, :) = [0.992157, 0.682353, 0.380392]; +ctrlPts(9, :) = [0.956863, 0.427451, 0.262745]; +ctrlPts(10, :) = [0.843137, 0.188235, 0.152941]; +ctrlPts(11, :) = [0.647059, 0, 0.14902]; + +X = 1:11; + +r = 1:11/R:11; + +R = interp1(X, ctrlPts(:, 1), r); +G = interp1(X, ctrlPts(:, 2), r); +B = interp1(X, ctrlPts(:, 3), r); + +result = [R' G' B']; + -- libgit2 0.21.4