diff --git a/HemoglobinSpectra.mat b/HemoglobinSpectra.mat
new file mode 100755
index 0000000000000000000000000000000000000000..5837407318d6e6e196b1333602aaca7226e30822
Binary files /dev/null and b/HemoglobinSpectra.mat differ
diff --git a/LICENSE b/LICENSE
new file mode 100755
index 0000000000000000000000000000000000000000..1c670dba37c07ba8265a9151e99d3ffe2a0a4a5c
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,15 @@
+MIT License
+
+Copyright (c) 2020, University of Texas Southwestern Medical Center.
+
+All rights reserved.
+
+Contributors: Devin O'Kelly
+
+Department: Department of Radiology
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
diff --git a/MakeSampleData.m b/MakeSampleData.m
new file mode 100644
index 0000000000000000000000000000000000000000..bc04f342f14aa88ae0d22c4c00c1bd725877dd8e
--- /dev/null
+++ b/MakeSampleData.m
@@ -0,0 +1,529 @@
+
+%% MakeSampleData generates datasets of a Shepp-Logan phantom according to concentration time course
+% of provided spectral endmembers.
+%
+%
+% ABOUT:
+% author - Devin O'Kelly
+% date - 28 Nov 2019
+% last update - 19 May 2020
+%
+%
+%
+
+
+
+ %% Global dataset parameters.
+clearvars; close all; clc;
+ rng(172969105)
+ datasetName = 'testData_';
+
+ saveFolder = './MultispectralSheppLoganData';
+
+ timecourseLength = 'short'; % 'long'
+ NWL = 5; % Number of wavelengths to sample between 680 and 980.
+ oversample = 1; % number of sequential times to repeat acquisition of a particular channel.
+ dynamics = 'local'; % 'local','global'
+ isbreathing = 1; % 1
+
+ dt = 0.1;
+ imageSize = 128;
+ WLs = [700 730 756 800 850];
+
+
+ % Creating the name.
+ LENG = [timecourseLength,'_'];
+ NCHAN = 'SelectedWavelengths_';
+ OSAMP = ['OS',num2str(oversample),'_'];
+ DYN = [dynamics,'dynam_'];
+ if isbreathing
+ BRTH = 'WBrth';
+ else
+ BRTH = 'NoBrth';
+ end
+
+ datasetName = [datasetName LENG NCHAN OSAMP DYN BRTH];
+ switch timecourseLength
+ case 'short'
+ Nt = 256;
+ case 'long'
+ Nt = 2048;
+ end
+ time = (1:Nt).*dt;
+ phantomType = 'Positive Shepp-Logan';
+
+
+
+ jitter = 0.0; % Variation in breathing rate.
+ BreathPeriod = 1; %Sec
+ SampleRate = 1./dt; %Hz
+ ImagingDuration = max(time); %Sec
+
+ %% Phantom distortion
+ if isbreathing
+ SpikeMeanHeight = 0.15; % Arbitrary unit.
+ SpikeVarHeight = 0.002;
+ else
+ SpikeMeanHeight = 0.0001; % Very small displacement values, when interpolated, will not cause an appreciable amount of variation in the edge profile.
+ SpikeVarHeight = 0.0000000001;
+ end
+ % TODO: Add random wavelength selection.
+ [Ny,Nx,N] = deal(imageSize);
+ WLorder = 1:numel(WLs);
+
+ if boolean(oversample)
+ WLorder = repelem(WLorder,oversample); %#ok<UNRCH>
+ WLorder = WLorder(:);
+ end
+
+
+
+ %% Ellipses/phantom shape and composition over time.
+ % For N ellipses and T timepoints and C components, need an N x (C x T)
+ % array to describe it.
+ % So for 10 ellipses, 4 components, and 250 time points, that's 10000
+ % entries.
+ switch phantomType
+ case 'Shepp-Logan'
+ case {'Modified Shepp-Logan','Positive Shepp-Logan'}
+
+ % SO2 values at each 'plateau' of the SO2 timecourse.
+ switch dynamics
+ case 'local' % Rather arbitrary values of change between each segment.
+ SO2Values = [0.4 0.6 0.4;
+ 1.0 0.5 0.0
+ 0.2 0.5 0.1
+ 0.5 0.5 0.6
+ 0.2 0.9 0.2
+ 0.1 0.1 0.9
+ 0.9 0.8 0.7
+ 1.0 1.0 1.0
+ 0.7 0.6 0.2
+ 0.9 0.2 0.9];
+
+ % Case for no dynamics.
+ case 'no'
+ SO2Values = [0.4 0.4 0.4; % Case of no variation or dynamics of the underlying signal.
+ 0.5 0.5 0.5
+ 0.1 0.1 0.1
+ 0.5 0.5 0.5
+ 0.9 0.9 0.9
+ 0.2 0.2 0.2
+ 0.7 0.7 0.7
+ 1.0 1.0 1.0
+ 0.6 0.6 0.6
+ 0.3 0.3 0.3];
+ % Case for global dynamics.
+ case 'global'
+ SO2Values = [0.3 0.7 0.4;
+ 0.3 0.7 0.4;
+ 0.3 0.7 0.4;
+ 0.3 0.7 0.4;
+ 0.3 0.7 0.4;
+ 0.3 0.7 0.4;
+ 0.3 0.7 0.4;
+ 0.3 0.7 0.4;
+ 0.3 0.7 0.4;
+ 0.3 0.7 0.4];
+ otherwise
+ error();
+ end
+
+ % SO2 switching times, relative to the overall challenge as [0,1]
+ % TODO: Randomize.
+ % TODO: What happens if they're 0 or 1 or out of order?
+
+ SO2Times = [1/3 2/3;
+ 1/3 2/3;
+ 1/3 2/3;
+ 1/3 2/3;
+ 1/3 2/3;
+ 1/3 2/3;
+ 1/3 2/3;
+ 1/3 2/3;
+ 1/3 2/3;
+ 1/3 2/3];
+
+ % Time constants of response. Constant for each component.
+ TimeConstants = [0.2
+ 0.3
+ 0.02
+ 0.04
+ 0.2
+ 0.1
+ 0.05
+ 0.4
+ 0.2
+ 0.1];
+ % TODO: Base time constants off of the number of components.
+% TimeConstants=2*rand(10,1);
+
+ end
+
+
+
+%% Assign into parameters.
+sampleData.parameters.sampling.Wavelengths = WLs;
+sampleData.parameters.sampling.TimeStep = dt;
+sampleData.parameters.sampling.Fs = SampleRate;
+
+sampleData.parameters.image.Ny = Ny;
+sampleData.parameters.image.Nx = Nx;
+
+sampleData.parameters.breathing.SpikeHeightMean = SpikeMeanHeight;
+sampleData.parameters.breathing.SpikeHeightVar = SpikeVarHeight;
+sampleData.parameters.breathing.TimeJitter = jitter;
+sampleData.parameters.breathing.BreathPeriod = BreathPeriod;
+
+
+sampleData.parameters.SO2.SO2Values = SO2Values;
+sampleData.parameters.SO2.SO2Times = SO2Times;
+sampleData.parameters.SO2.TimeConstants = TimeConstants;
+
+sampleData.timeCourses.Time = time;
+
+%% Generate ellipses and the overall image.
+ [phantomImage,ellipses,phantomLayers]=phantomEllipses(Ny,...
+ phantomType);
+
+sampleData.parameters.image.PhantomImage = phantomImage;
+sampleData.parameters.image.Ellipses = ellipses;
+sampleData.parameters.image.PhantomLayers = phantomLayers;
+
+
+%% Create the timecourse for each component.
+nEllipses = size(phantomLayers,3);
+SO2Times = [SO2Times ones(nEllipses,1)]; % Closes SO2 times on the right.
+timeCourses = zeros(nEllipses,numel(time)); % Storage for SO2 timecourses.
+
+
+% Allocation.
+ tcIm = zeros(N,N,numel(time));
+ spectralImageTimeCourse_distorted = zeros(N,N,numel(WLs),numel(time));
+ tcImStore = cell(nEllipses,1);
+ spectralCourse = zeros(numel(WLs),numel(time),nEllipses);
+
+% Create smoothed timecourses.
+for k=1:nEllipses
+ % Allocation.
+ oxygenTC = zeros(size(time));
+
+ % Assign the values of the oxygenation between the time switches.
+ % TODO: Clarify logic here.
+ for switchInd=size(SO2Times,2):-1:1
+ oxygenTC(time<=(SO2Times(k,switchInd).*max(time)))=SO2Values(k,switchInd);
+ end
+
+ % Generate this timecourse's convolution kernel.
+ kern=exp(-TimeConstants(k)*time);
+ kern=kern./sum(kern); % normalize
+ kern=[zeros(size(kern)) kern]; % Pad up front. Necessary?
+ kern = smoothdata(kern);
+ oxygenTC=convWFrontPad(oxygenTC,kern);
+
+ % Store the TC and its spectrum in the library.
+ timeCourses(k,:)=oxygenTC;
+ spectralCourse(:,:,k)=makeSpectralSet(oxygenTC,WLs);
+end
+
+sampleData.timeCourses.SO2TimeCourses = timeCourses;
+sampleData.timeCourses.SpectralTimeCourses = spectralCourse;
+
+%% Prepare for simulating movement.
+[distMap,zz] = makeDistortion(N);
+ interpolator=makeInterpolator(zz,phantomLayers); % TODO: check that the interpolator works for both stacks and not-stacks.
+ nEllipses=size(ellipses,1);
+
+sampleData.parameters.distortion.DistortionMap = distMap;
+sampleData.parameters.distortion.DistortionCoordinates = zz;
+
+% Make grid vectors and interpolating grids.
+xVec = real(zz(1,:));
+yVec = imag(zz(:,1));
+zVec = permute(1:nEllipses,[1 3 2]); % Puts the non-singleton dimension in dim 3.
+
+[X,Y,Z] = ndgrid(xVec,yVec,zVec);
+
+Zmix = X+1i*Y;
+
+
+%% Make breathing timecourse
+
+% TODO: Make this simpler and more debuggable.
+q2 = makeBarTimeCourse(BreathPeriod,jitter,time,SampleRate,SpikeVarHeight,SpikeMeanHeight);
+q = imgaussfilt(q2);
+
+sampleData.timeCourses.BreathingTimeCourse = q;
+
+trueSpectralImage_distorted = zeros(N,N,nEllipses,numel(WLs),numel(time));
+trueSpectralImage = zeros(N,N,nEllipses,numel(WLs),numel(time));
+
+
+%% Distort images and combine into spectral image. Also record the undistorted version.
+
+ for t=1:numel(time)
+ % grab the breathing timecourse and apply the distortion.
+ dx=q(t);
+ Zdistort=Zmix+dx*1i*distMap; % TODO: Check the math here.
+ distortedIm=interpolator(real(Zdistort),imag(Zdistort),Z);
+
+ % grab the spectral timecourse for each component.
+ expander = squeeze(spectralCourse(:,t,:));
+ expander = permute(expander,[4 3 2 1]);
+
+ % Take tensor product of vectors, create representative image.
+ totalIm_distorted = distortedIm.*expander;
+ totalIm_true = phantomLayers.*expander;
+
+ trueSpectralImage_distorted(:,:,:,:,t) = totalIm_distorted;
+ trueSpectralImage(:,:,:,:,t) = totalIm_true;
+
+ spectralImage_distorted =squeeze(sum(totalIm_distorted,3));
+ spectralImageTimeCourse_distorted(:,:,:,t)=spectralImage_distorted;
+
+
+ end
+
+%% Save component time-course image.
+sampleData.images.CompleteData_distorted = trueSpectralImage_distorted;
+sampleData.images.CompleteData = trueSpectralImage;
+
+
+%% Simulate the process of acquiring data from the object distribution.
+
+% Allocation.
+sampledImageTimeCourse=zeros(N,N,numel(time));
+
+% Set the wavelength sampling index.
+WLindex=1;
+WLlog=zeros(size(time)); % Record which wavelength was used when.
+
+
+%% Loop over the whole time series and assign the sampled image.
+for k=1:numel(time)
+
+
+ if WLindex>numel(WLorder) % If we've already cycled through the indices...
+ WLindex=1; % Bring back to 1.
+ end
+
+ % Assign the image to the data store.
+ % TODO: Simulate noise (acoustic and electrical).
+ % TODO: With noise, get SNR information.
+ sampledImageTimeCourse(:,:,k)=spectralImageTimeCourse_distorted(:,:,WLorder(WLindex),k);
+
+ % Store the wavelength we acquired with.
+ WLlog(k)=WLs(WLorder(WLindex));
+ WLindex=WLindex+1;
+end
+
+sampleData.parameters.sampling.SamplingOrder = WLorder;
+sampleData.parameters.sampling.SamplingWavelengths = WLs(WLorder);
+
+sampleData.images.SampledData = sampledImageTimeCourse;
+sampleData.timeCourses.Wavelengths = WLlog;
+
+%% Save the object distribution, as well as the sampled distribution.
+
+CD_distorted = squeeze(sum(sum(sampleData.images.CompleteData_distorted,3),4));
+CD = squeeze(sum(sum(sampleData.images.CompleteData,3),4));
+
+
+% make the local directory (just use a temp folder)
+% tdir_local=tempname(fullfile(pwd,'genData'));
+tdir_local = fullfile(saveFolder,datasetName);
+mkdir(tdir_local);
+cd(tdir_local);
+
+% make MAT files for storage
+% TODO: Consolidate into a structure.
+save('testData.mat','sampleData','-v7.3');
+
+%% Create a video of the sampled dataset.
+
+% Quick fun thing: Can alias a deep handle reference by making an anonymous
+% function.
+samp = @() sampleData.images.SampledData; % Neat!
+
+v = VideoWriter([datasetName '.avi']);
+v.FrameRate = SampleRate;
+open(v);
+
+imToWrite = normalizeData(permute(samp(),[1 2 4 3]));
+
+writeVideo(v,imToWrite);
+
+close(v);
+
+
+%% Go home.
+% cd('..\..');
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% HELPER FUNCTIONS %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% Normalize data to [0,1]
+function dat = normalizeData(input)
+
+ minX = min(input(:));
+ maxX = max(input(:));
+
+ dat = (input-minX)./(maxX-minX);
+
+end
+
+
+%% Overlay perlin noise on the image provided.
+function im = perlin_noise(im)
+
+ [n, m] = size(im);
+ i = 0;
+ w = sqrt(n*m);
+
+ while w > 3
+ i = i + 1;
+ d = interp2(randn(ceil((n-1)/(2^(i-1))+1),ceil((m-1)/(2^(i-1))+1))*4, i-1, 'spline');
+ im = im + i * d(1:n, 1:m);
+ w = w - ceil(w/2 - 1);
+ end
+end
+
+%% Makes the spectral timeseries corresponding to the given SO2 timeseries.
+function oxyCourse = makeSpectralSet(SO2Course,WLs)
+
+ load HemoglobinSpectra.mat
+
+ Hb_lambda = Hb.spectra(:,1); % Get wavelengths used for Hb spectra (every 2nm)
+ ox_abs = Hb.spectra(:,2); % Get extinction vals of HbO spectra (every 2nm)
+ dox_abs = Hb.spectra(:,3); % Get extinction vals of HbD spectra (every 2nm)
+
+ % Resample all of the spectra to give us just the wavelengths we're
+ % interested in.
+ ox_spec = interp1(Hb_lambda,ox_abs,WLs);
+ dox_spec = interp1(Hb_lambda,dox_abs,WLs);
+
+ specMat = [ox_spec(:),dox_spec(:)]; % Concatenate into a mixing matrix
+
+ HbOConc = SO2Course;
+ HbConc = 1 - SO2Course;
+
+ concMat = [HbOConc(:).';HbConc(:).'];
+ oxyCourse = specMat*concMat;
+
+end
+
+
+
+%% Returns the time-stabilized timecourse of oxygenation.
+function output = convWFrontPad(target,kern)
+
+
+ N = numel(target);
+
+ % Expand out the timecourse by padding the front.
+ target2 = [ones(size(target))*target(1), target,ones(size(target))*target(end)];
+
+ tempconv = conv(target2,kern,'same');
+
+ output=tempconv((N+1):(end-N));
+end
+
+
+%% Creates a distortion map of the given image.
+% Based on an elliptical Mobius transform.
+% TODO: Stabilize a little better, figure out the conjugacy/transpose issue.
+function [distortionMap,zz] = makeDistortion(Npix)
+
+ [xx,yy]=meshgrid(linspace(-2,2,Npix));
+ zz=xx+1i*yy;
+
+ Mmat=[-exp(1i*pi/4) 0.1i;
+ 0 exp(1i*pi/4)]
+
+ distortionMap = mobiusTransform(zz,Mmat)+2i-real(zz);
+ distortionMap = -distortionMap./abs(distortionMap);
+
+end
+
+%% Makes a GriddedInterpolant object to allow resampling of image coordinates.
+function interpolator = makeInterpolator(zz,im)
+ isStack = (ndims(im)==3);
+ if isStack
+ stackSize = size(im,3);
+
+ zVec = permute([1:stackSize],[1 3 2]);
+
+ gridVecs = {real(zz(1,:)),imag(zz(:,1)),zVec};
+
+ interpolator = griddedInterpolant(gridVecs,im);
+ else
+ interpolator = griddedInterpolant(real(zz)',imag(zz)',im');
+ end
+end
+
+%% Applies the Mobius transform to the given z array.
+function output = mobiusTransform(z,xfrm)
+
+ a=xfrm(1,1);
+ b=xfrm(1,2);
+ c=xfrm(2,1);
+ d=xfrm(2,2);
+
+ output = (a.*z+b)./(c.*z+d);
+end
+
+%% Creates a timecourse of breathing for simulating motion.
+function [tc] = makeBarTimeCourse(BreathPeriod,jitter,time,SampleRate,SpikeVarHeight,SpikeMeanHeight)
+
+ % First localize the variables.
+ T = BreathPeriod;
+ fs = SampleRate;
+ tc = zeros(size(time));
+
+ % Get parameters for inter-breath time. This is index, so convert.
+ meanInterTime = BreathPeriod*fs;
+ if jitter == 0
+ varInterTime = meanInterTime./1e6;
+ else
+ varInterTime = jitter*(fs.^2);
+ end
+
+ b_time = varInterTime./meanInterTime;
+ a_time = meanInterTime./b_time;
+
+ % Get parameters for breathing motion.
+ meanMotion = SpikeMeanHeight;
+ if SpikeVarHeight == 0
+ varMotion = SpikeMeanHeight./1e6;
+ else
+ varMotion = SpikeVarHeight;
+ end
+
+ b_motion = varMotion./meanMotion;
+ a_motion = meanMotion./b_motion;
+
+ % Initial time point. Want to not deal weird edge artifacts or
+ % negative time.
+ tind = round(gamrnd(a_time,b_time));
+
+ while tind < max(time)*fs
+
+ % Put a spike of random height at the timepoint given.
+ tc(tind) = gamrnd(a_motion,b_motion);
+
+ % advance to the next time point.
+ tind = tind + round(gamrnd(a_time,b_time));
+ end
+end
+
+
+
+
+
+
+
+
+
+
diff --git a/README.md b/README.md
old mode 100644
new mode 100755
index 5e4d721f07ae7cdf7a2d83dacdef20f5673457f9..dd8a7ed9830533db887865b921bce0495cd584bc
--- a/README.md
+++ b/README.md
@@ -1,3 +1,13 @@
# Kalata Filters - MATLAB
-Implementations of the Kalata filters described in
\ No newline at end of file
+Implementations of the Kalata filters described in:
+
+D.S. O'Kelly, Y. Guo, R.P. Mason. Evaluating online filtering algorithms to enhance dynamic multispectral optoacoustic tomography, Photoacoustics (2020)
+
+INSTRUCTIONS FOR USE:
+
+1. Run MakeSampleData.m to create test data for running the Kalata filter.
+
+2. Run filterEx.m with varying parameters to test the Kalata filtering algorithm.
+
+If any substantial portion(s) of this code, not including the use of functions derived from phantomEllipses.m, is/are used in another publication or further development, please cite the above publication.
\ No newline at end of file
diff --git a/filterEx.m b/filterEx.m
new file mode 100644
index 0000000000000000000000000000000000000000..66204f6c50e5f37a3305361c23a8dff1bd34aa23
--- /dev/null
+++ b/filterEx.m
@@ -0,0 +1,77 @@
+
+% Given an Ny x Nx x Nt matrix, an Nt x 1 vector of channel numbers, and
+% measures of dt,ow, and on, produce a filtered time-series.
+
+
+NOISE_STD = 500;
+
+burn_in = 3;
+filt_order = 0;
+proc_noise_var = 300^2;
+meas_noise_var = 700^2;
+
+
+dataFolder = './MultispectralSheppLoganData/testData_short_SelectedWavelengths_OS1_localdynam_WBrth';
+
+if ~exist('sampleData','var')||~(strcmp(pwd,dataFolder))
+% cd(dataFolder);
+ load(fullfile(dataFolder,'testData.mat'));
+else
+ clearvars -except sampleData dataFolder
+end
+
+
+% The scanData is the one-wavelength-at-a-time train of data.
+scanData = sampleData.images.SampledData;
+scanData = scanData + NOISE_STD * randn(size(scanData));
+
+% Get the whole dataset as well as whatever meta is necessary.
+completeDataset = sampleData.images.CompleteData;
+scanMeta = sampleData.parameters;
+channelSet = scanMeta.sampling.Wavelengths;
+Wavelengths = channelSet;
+frameSampledWLs = sampleData.timeCourses.Wavelengths;
+
+
+
+timeList = sampleData.timeCourses.Time;
+[~,channelIndexList] = ismember(frameSampledWLs,Wavelengths);
+
+
+SW_filteredData = kalataFilter(scanData,channelIndexList,timeList,proc_noise_var,meas_noise_var,0,burn_in);
+alpha_filteredData = kalataFilter(scanData,channelIndexList,timeList,proc_noise_var,meas_noise_var,1,burn_in);
+alphabeta_filteredData = kalataFilter(scanData,channelIndexList,timeList,proc_noise_var,meas_noise_var,2,burn_in);
+trueData = squeeze(sum(completeDataset,3));
+%%
+for k = 1:numel(timeList)
+ ssim_SW(k) = calcWholeImageSSIM(SW_filteredData(:,:,:,k),trueData(:,:,:,k));
+ ssim_a(k) = calcWholeImageSSIM(alpha_filteredData(:,:,:,k),trueData(:,:,:,k));
+ ssim_ab(k) = calcWholeImageSSIM(alphabeta_filteredData(:,:,:,k),trueData(:,:,:,k));
+
+end
+figure;
+plot(ssim_SW); hold on; plot(ssim_a); plot(ssim_ab);
+
+
+
+
+
+% Functions as supplement to MATLAB's built-in SSIM calculations.
+function wholessim = calcWholeImageSSIM(im,ref)
+ C1 = 0.02;
+ C2 = 0.01;
+
+ mean_ref = mean(ref(:));
+ mean_im = mean(im(:));
+
+ covMat = cov(ref(:),im(:));
+
+ Nt1 = (2*mean_ref*mean_im + C1);
+ Nt2 = (2*covMat(1,2) + C2);
+ Dt1 = (mean_ref.^2 + mean_im.^2 + C1);
+ Dt2 = (covMat(1,1) + covMat(2,2) + C2);
+
+ wholessim = (Nt1.*Nt2)./(Dt1.*Dt2);
+
+
+end
diff --git a/kalataFilter.m b/kalataFilter.m
new file mode 100644
index 0000000000000000000000000000000000000000..e2c52f19532d91bf03b170198a64e589b215a7c3
--- /dev/null
+++ b/kalataFilter.m
@@ -0,0 +1,263 @@
+
+function [filteredOutput] = kalataFilter(inputImages,channelIndexList,timeArray,sigma_w,sigma_n,filt_order,burn_in)
+%kalataFilter.m Apply Kalata filters to an input time-series of images
+%acquired at single channels.
+%
+%
+%
+% INPUTS:
+% inputImages : An Ny*Nx*Nt array of images, with each Ny*Nx image acquired at a specific channel
+% recorded in channelIndexList, such that inputImages(:,:,k) is the image recorded at
+% channelIndexList(k)
+% channelIndexList: An Nt*1 array of doubles, one entry for each image, with a unique linear index for each channel.
+% Values are expected to lie between 1 and NChannels.
+% timeArray: An Nt*1 array of doubles, corresponding to the time points at which each image was acquired.
+% sigma_w: Process noise variance
+% sigma_n: Measurement noise variance
+% filt_order: Order of the Kalata filter, with values:
+% 0 - Sliding Window (SW) filter
+% 1 - Alpha filter
+% 2 - Alpha-Beta filter
+% burn_in: Number of frames of each channel to record before beginning to apply the kinematic model.
+% Integer-valued
+%
+% OUTPUTS:
+% filteredOutput: An Ny*Nx*Nc*Nt array of filtered multispectral images, such that each Ny*Nx*Nc multi-channel image is the estimate of the true multi-channel image at timepoint k.
+%
+%
+% author - Devin O'Kelly
+% date - 28 Nov 2019
+% last update - 19 May 2020
+%
+
+ % Preparing for buffer creation.
+ Nw = numel(unique(channelIndexList));
+ [Ny,Nx,NFrames] = size(inputImages);
+
+ filteredOutput = ones(Ny,Nx,Nw,NFrames);
+ lastTime = 0;
+ channelTimeList = zeros(Nw,1);
+ nMeas = zeros(Nw,1);
+
+ switch filt_order
+ case 0 % Sliding window
+ est.state = zeros(Ny,Nx,Nw);
+ case 1 % Alpha
+ est.state = zeros(Ny,Nx,Nw);
+ case 2 % Alpha-beta
+ est.state = zeros(Ny,Nx,Nw);
+ est.d_state = zeros(Ny,Nx,Nw);
+ otherwise
+ error('Invalid filter order!');
+ end
+
+
+ %%
+ for k = 1:NFrames
+ %%% Simulate data acquisition, along with channel index and timestamp
+ frameData = inputImages(:,:,k);
+ channelIndex = channelIndexList(k);
+ curTime = timeArray(k);
+
+
+ dt = curTime - lastTime; % Time update for kinematic component.
+ dT = curTime - channelTimeList(channelIndex); % Sampling time since the last time this channel was observed.
+
+
+ % Kinematic update / predict step
+ if nMeas(channelIndex) < burn_in
+ switch filt_order
+ case 0
+ est.state = est.state;
+ case 1
+ est.state = est.state;
+ case 2
+ est.state = est.state;
+ est.d_state = est.d_state.*0;
+ end
+ else
+ switch filt_order
+ case 0
+ est.state = est.state;
+ case 1
+ est.state = est.state;
+ case 2
+ est.state = est.state + dt.*est.d_state;
+ est.d_state = est.d_state;
+ end
+
+ end
+
+ % Calculate tracking index based on time since last observation.
+ trackingIndex = calcLambda(dT,sigma_w,sigma_n);
+
+ switch filt_order
+ case 0
+ alpha = 1;
+ case 1
+ italpha = calcABG_from_lambda(trackingIndex);
+ [alpha] = calc_scheduling(nMeas(channelIndex),italpha);
+ case 2
+ [italpha,itbeta] = calcABG_from_lambda(trackingIndex);
+ [alpha,beta] = calc_scheduling(nMeas(channelIndex),italpha,itbeta);
+ end
+
+
+ % Take measurement residual.
+ channelDiff = frameData - est.state(:,:,channelIndex);
+
+
+ if nMeas(channelIndex) == 0 % Let the first estimate take on the first value.
+ est.state(:,:,channelIndex) = frameData;
+ switch filt_order
+ case {2}
+ est.d_state(:,:,channelIndex) = est.d_state(:,:,channelIndex).*0;
+ otherwise
+ end
+ else
+ switch filt_order
+ case 0
+ est.state(:,:,channelIndex) = frameData;
+ case 1
+ est.state(:,:,channelIndex) = est.state(:,:,channelIndex) + alpha.*channelDiff;
+ case 2
+ est.state(:,:,channelIndex) = est.state(:,:,channelIndex) + alpha.*channelDiff;
+ est.d_state(:,:,channelIndex) = est.d_state(:,:,channelIndex) + (beta./dT).*channelDiff;
+ end
+ end
+
+ channelTimeList(channelIndex) = curTime;
+ lastTime = curTime;
+ nMeas(channelIndex) = nMeas(channelIndex) + 1;
+
+ filteredOutput(:,:,:,k) = est.state;
+
+
+ if rem(k,50) == 0
+ disp([num2str(k) ' completed']);
+ end
+
+ end
+end
+
+
+
+
+
+
+
+% calculating tracking index
+% Expression from:
+% Paul R. Kalata: The tracking index: A generalized parameter for alpha-beta and alpha-beta-gamma target trackers.
+% IEEE Transactions on Aerospace and Electronic Systems, AES-20(2):174-181, March 1984.
+function lambda = calcLambda(dt,process_var,noise_var)
+
+ sw = sqrt(process_var);
+ sn = sqrt(noise_var);
+
+ lambda = sw.*(dt.^2)./sn;
+
+end
+
+% Calculation of alpha, alpha-beta, and alpha-beta-gamma coefficients.
+%
+% Expressions for A and A-B coefficients from:
+% Paul R. Kalata: The tracking index: A generalized parameter for alpha-beta and alpha-beta-gamma target trackers.
+% IEEE Transactions on Aerospace and Electronic Systems, AES-20(2):174-181, March 1984.
+%
+% Expressions for A-B-G coefficients from:
+% J. E. Gray and W. J. Murray: A derivation of an analytic expression for the tracking index for the alpha-beta-gamma filter.
+% IEEE Trans. on Aerospace and Electronic Systems, 29:1064-1065, 1993.
+%
+% Input: A single scalar value (the tracking index) output from
+% calcLambda()
+%
+% Output: Coefficients for 0th, 1st, or 2nd-order Kalata filters depending
+% on the number of outputs denoted.
+
+function [alpha, beta, gamma] = calcABG_from_lambda(lambda)
+
+ switch nargout
+ case 1
+ alpha = (-(lambda.^2) + sqrt(lambda.^4 + 16.*(lambda.^2)))./8;
+ case 2
+ r = (4 + lambda - sqrt(8.*lambda + lambda.^2))./(4);
+ alpha = 1 - r.^2;
+ beta = 2.*(2-alpha) - 4.*sqrt(1-alpha);
+
+ case 3
+ b = lambda/2 - 3;
+ c = lambda/2 +3;
+ d = -1;
+ p = c - (b.^2)/3;
+ q = (2.*b.^3)./27 - (b.*c)/3 + d;
+ v = sqrt(q.^2 + (4.*p.^3)/27);
+ z = -(q + v/2).^(1/3);
+ s = z - p./(3.*z)-b./3;
+ alpha = 1 - s.^2;
+ beta = 2.*(1-s).^2;
+ gamma = (beta.^2)./(2.*alpha);
+ otherwise
+ disp('Number of output arguments specifies filter order.');
+ end
+
+end
+
+
+% Calculating the scheduling based on the sample index.
+% Functions based on expressions from:
+% Fixed-Gain, Two-Stage Estimators For Tracking Maneuvering Targets, W.D.
+% Blair, NSWC Dahlgren Division: https://apps.dtic.mil/dtic/tr/fulltext/u2/a255832.pdf
+% Link current as of 2020 May 19
+function [varargout] = calc_scheduling(k,varargin)
+
+ nvarg = numel(varargin);
+ switch nvarg
+ case 1
+ [alpha] = varargin{:};
+ if isnan(alpha)
+ alpha = -Inf;
+ end
+
+ itValue_alpha = 1./(k+1);
+ alpha_k = max(itValue_alpha,alpha);
+ varargout = {alpha_k};
+ case 2
+ [alpha, beta] = varargin{:};
+ if isnan(alpha)
+ alpha = -Inf;
+ end
+ if isnan(beta)
+ beta = -Inf;
+ end
+
+ itValue_alpha = (2.*(2.*k + 1)./((k+1).*(k+2)));
+ itValue_beta = (6./((k+1).*(k+2)));
+
+ alpha_k = max(itValue_alpha,alpha);
+ beta_k = max(itValue_beta,beta);
+ varargout = {alpha_k beta_k};
+ case 3
+ [alpha, beta, gamma] = varargin{:};
+ if isnan(alpha)
+ alpha = -Inf;
+ end
+ if isnan(beta)
+ beta = -Inf;
+ end
+ if isnan(gamma)
+ gamma = -Inf;
+ end
+
+ itValue_alpha = (3.*(3.*k.^2 + 3.*k + 2)./((k+1).*(k+2).*(k+3)));
+ itValue_beta = (18.*(2.*k + 1)./((k+1).*(k+2).*(k+3)));
+ itValue_gamma = (60)./((k+1).*(k+2).*(k+3));
+
+ alpha_k = max(itValue_alpha,alpha);
+ beta_k = max(itValue_beta,beta);
+ gamma_k = max(itValue_gamma,gamma);
+ varargout = {alpha_k beta_k gamma_k};
+ otherwise
+ end
+end
+
diff --git a/phantomEllipses.m b/phantomEllipses.m
new file mode 100755
index 0000000000000000000000000000000000000000..1499dca9d83a18f6fa37b4f0afd73972ee80667c
--- /dev/null
+++ b/phantomEllipses.m
@@ -0,0 +1,201 @@
+function [p,ellipse,holderp]=phantomEllipses(varargin)
+%PHANTOM Generate a head phantom image.
+% P = PHANTOM(DEF,N) generates an image of a head phantom that can
+% be used to test the numerical accuracy of RADON and IRADON or other
+% 2-D reconstruction algorithms. P is a grayscale intensity image that
+% consists of one large ellipse (representing the brain) containing
+% several smaller ellipses (representing features in the brain).
+%
+% DEF is a string that specifies the type of head phantom to generate.
+% Valid values are:
+%
+% 'Shepp-Logan' A test image used widely by researchers in
+% tomography
+% 'Modified Shepp-Logan' (default) A variant of the Shepp-Logan phantom
+% in which the contrast is improved for better
+% visual perception.
+%
+% N is a scalar that specifies the number of rows and columns in P.
+% If you omit the argument, N defaults to 256.
+%
+% P = PHANTOM(E,N) generates a user-defined phantom, where each row
+% of the matrix E specifies an ellipse in the image. E has six columns,
+% with each column containing a different parameter for the ellipses:
+%
+% Column 1: A the additive intensity value of the ellipse
+% Column 2: a the length of the horizontal semi-axis of the ellipse
+% Column 3: b the length of the vertical semi-axis of the ellipse
+% Column 4: x0 the x-coordinate of the center of the ellipse
+% Column 5: y0 the y-coordinate of the center of the ellipse
+% Column 6: phi the angle (in degrees) between the horizontal semi-axis
+% of the ellipse and the x-axis of the image
+%
+% For purposes of generating the phantom, the domains for the x- and
+% y-axes span [-1,1]. Columns 2 through 5 must be specified in terms
+% of this range.
+%
+% [P,E] = PHANTOM(...) returns the matrix E used to generate the phantom.
+%
+% Class Support
+% -------------
+% All inputs must be of class double. All outputs are of class double.
+%
+% Remarks
+% -------
+% For any given pixel in the output image, the pixel's value is equal to the
+% sum of the additive intensity values of all ellipses that the pixel is a
+% part of. If a pixel is not part of any ellipse, its value is 0.
+%
+% The additive intensity value A for an ellipse can be positive or negative;
+% if it is negative, the ellipse will be darker than the surrounding pixels.
+% Note that, depending on the values of A, some pixels may have values outside
+% the range [0,1].
+%
+% See also RADON, IRADON.
+
+% Copyright 1993-2003 The MathWorks, Inc.
+% $Revision: 1.13.4.2 $ $Date: 2003/08/01 18:09:35 $
+
+
+% References:
+% A. K. Jain, "Fundamentals of Digital Image Processing", p. 439.
+% P. A. Toft, "The Radon Transform, Theory and Implementation" (unpublished
+% dissertation), p. 199.
+
+[ellipse,n] = parse_inputs(varargin{:});
+
+p = zeros(n);
+pe= zeros(n);
+holderp=zeros(n,n,size(ellipse,1));
+
+xax = ( (0:n-1)-(n-1)/2 ) / ((n-1)/2);
+xg = repmat(xax, n, 1); % x coordinates, the y coordinates are rot90(xg)
+
+for k = 1:size(ellipse,1)
+ asq = ellipse(k,2)^2; % a^2
+ bsq = ellipse(k,3)^2; % b^2
+ phi = ellipse(k,6)*pi/180; % rotation angle in radians
+ x0 = ellipse(k,4); % x offset
+ y0 = ellipse(k,5); % y offset
+ A = ellipse(k,1); % Amplitude change for this ellipse
+ x=xg-x0; % Center the ellipse
+ y=rot90(xg)-y0;
+ cosp = cos(phi);
+ sinp = sin(phi);
+ idx=find(((x.*cosp + y.*sinp).^2)./asq + ((y.*cosp - x.*sinp).^2)./bsq <= 1);
+ p(idx) = p(idx) + A;
+ pe(idx)=p(idx)+A;
+ holderp(:,:,k)=pe;
+ pe(:)=0;
+end
+
+
+function [e,n] = parse_inputs(varargin)
+% e is the m-by-6 array which defines ellipses
+% n is the size of the phantom brain image
+
+n=256; % The default size
+e = [];
+defaults = {'shepp-logan', 'modified shepp-logan','positive shepp-logan'};
+
+for i=1:nargin
+ if ischar(varargin{i}) % Look for a default phantom
+ def = lower(varargin{i});
+ idx = strmatch(def, defaults);
+ if isempty(idx)
+ eid = sprintf('Images:%s:unknownPhantom',mfilename);
+ msg = 'Unknown default phantom selected.';
+ error(eid,'%s',msg);
+ end
+ switch defaults{idx}
+ case 'shepp-logan'
+ e = shepp_logan;
+ case 'modified shepp-logan'
+ e = modified_shepp_logan;
+ case 'positive shepp-logan'
+ e = positive_shepp_logan;
+ end
+ elseif numel(varargin{i})==1
+ n = varargin{i}; % a scalar is the image size
+ elseif ndims(varargin{i})==2 && size(varargin{i},2)==6
+ e = varargin{i}; % user specified phantom
+ else
+ eid = sprintf('Images:%s:invalidInputArgs',mfilename);
+ msg = 'Invalid input arguments.';
+ error(eid,'%s',msg);
+ end
+end
+
+if isempty(e) % ellipse is not yet defined
+ e = modified_shepp_logan;
+end
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Default head phantoms:
+%
+
+function shep=shepp_logan
+%
+% This is the default head phantom, taken from AK Jain, 439.
+%
+% A a b x0 y0 phi
+% ---------------------------------
+shep = [ 1 .69 .92 0 0 0
+ -.98 .6624 .8740 0 -.0184 0
+ -.02 .1100 .3100 .22 0 -18
+ -.02 .1600 .4100 -.22 0 18
+ .01 .2100 .2500 0 .35 0
+ .01 .0460 .0460 0 .1 0
+ .01 .0460 .0460 0 -.1 0
+ .01 .0460 .0230 -.08 -.605 0
+ .01 .0230 .0230 0 -.606 0
+ .01 .0230 .0460 .06 -.605 0 ];
+
+
+function toft=modified_shepp_logan
+%
+% This head phantom is the same as the Shepp-Logan except
+% the intensities are changed to yield higher contrast in
+% the image. Taken from Toft, 199-200.
+%
+% A a b x0 y0 phi
+% ---------------------------------
+toft = [ 1 .69 .92 0 0 0
+ -.8 .6624 .8740 0 -.0184 0
+ -.2 .1100 .3100 .22 0 -18
+ -.2 .1600 .4100 -.22 0 18
+ .1 .2100 .2500 0 .35 0
+ .1 .0460 .0460 0 .1 0
+ .1 .0460 .0460 0 -.1 0
+ .1 .0460 .0230 -.08 -.605 0
+ .1 .0230 .0230 0 -.606 0
+ .1 .0230 .0460 .06 -.605 0 ];
+
+
+ function toft=positive_shepp_logan
+%
+% This head phantom is the same as the Shepp-Logan except
+% the intensities are changed to yield higher contrast in
+% the image. Taken from Toft, 199-200.
+%
+% A a b x0 y0 phi
+% ---------------------------------
+toft = [ 1 .69 .92 0 0 0
+ 0.8 .6624 .8740 0 -.0184 0
+ 0.2 .1100 .3100 .22 0 -18
+ 0.2 .1600 .4100 -.22 0 18
+ .1 .2100 .2500 0 .35 0
+ .1 .0460 .0460 0 .1 0
+ .1 .0460 .0460 0 -.1 0
+ .1 .0460 .0230 -.08 -.605 0
+ .1 .0230 .0230 0 -.606 0
+ .1 .0230 .0460 .06 -.605 0 ];
+
+
+
+
+
+