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