classdef TotalFocusingMethod < handle
	% TotalFocusingMethod	Performs the total focusing method (TFM) on full matrix capture (FMC) data (version 0.1020140721).
	% Simplifies coverting full matrix capture (FMC) data to an image of the sample using the total focusing method (TFM).
	% This documentation is not complete. If adding to the documentation, please follow the format described in:
	% http://www.mathworks.co.uk/help/matlab/matlab_prog/create-help-for-classes.html
	%
	% TotalFocusingMethod	Copyright:
	% 
	% Copyright (c) 2014, Phillip Anthony Petcher
	% All rights reserved.
	% 
	% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
	% 
	% - Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
	% - Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
	% 
	% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
	% 
	% TotalFocusingMethod	Properties:
	%	Image - (get only) - the last total focusing method (TFM) image created.
	%	DataAnalytic - (get only) - the input signal after it has been band pass filtered, interpolated to a higher sampling frequency, and converted to the analytic signal form.
	%	Speed - (get only) - the wave propagation speed in m/s.
	%	DataLength - (get only) - the number of time points for a given capture for any of the receivers.
	%	StartTime - (get only) - the time value in s of the first time point in any given capture, which must be greater than zero.
	%	HighPass - (get only) - the high pass filter frequency in Hz (not MHz), which must be greater than zero.
	%	LowPass - (get only) - the low pass filter frequency in Hz (not MHz), which must be greater than HighPass, and less than a quarter of the sampling frequency.
	%	SamplingFreq - (get only) - the sampling frequency in Hz (not MHz) that the data was captured at.
	%	SamplingFreqInterp - (get only) - the sampling frequency in Hz (not MHz) that the data is interpolated to.
	%	PadLengthFreq - (get only) - the size to pad the frequency domain signal to when calculating the analytic signal.
	%	PadLengthTime - (get only) - the size to pad the time domain signal to when calculating the analytic signal.
	%	FilterStart - (get only) - the last point before the band pass region of the signal when filtering it.
	%	FilterEnd - (get only) - the first point after the band pass region of the signal when filtering it.
	%	LastPositive - (get only) - the last positive point in the FFT when calculating the analytic signal.
	%	PositionsX - (get only) - the X positions of the transducer elements in m.  Note that for a linear array, this is usually all that is necessary.
	%	PositionsY - (get only) - the Y positions of the transducer elements in m.  Note that this is the position in the vertical (depth) direction, and usually will just be zeros, or left blank.
	%	PositionsZ - (get only) - the Z positions of the transducer elements in m.
	%	Elements - (get only) - the total number of transducer elements.
	%	RegionX - (get only) - the X points that are to be imaged with the total focusing method (TFM).
	%	RegionY - (get only) - the Y points that are to be imaged with the total focusing method (TFM).
	%	RegionZ - (get only) - the Z points that are to be imaged with the total focusing method (TFM).  This will not be used for a two dimensional image, and usually will just be zeros, or left blank.
	%	Paths - (get only) - the distances in m from each transducer to each point that is to be imaged.
	%	RefreshPaths = true - (get only) - if a parameter change means that the path lengths need to be calculated again.
	%	RefreshAnalytic = true - (get only) - if a parameter change means that the analytic signals need to be calculated again.
	%
	% TotalFocusingMethod	Methods:
	%	SetDataParameters(DataLength,StartTime,SamplingFreq,HighPass,LowPass,Speed) - set the basic parameters which describe the data to be processed.
	%	SetTransducerPosition(PositionsX,PositionsY,PositionsZ) - set the transducer positions, usually relative to the centre of the transducer, but any reference point can be used as long as it is consistent.
	%	SetImageRegion(XStart,XEnd,XIncrement,YStart,YEnd,YIncrement,ZStart,ZEnd,ZIncrement) - set the region that is to be imaged with the total focusing method (TFM).
	%	CalculatePaths() - calculates the paths from each transducer to each point that is to be imaged.
	%	CalculatePathsWedge(SpeedWedge,Interface) - calculates the paths from each transducer to each point that is to be imaged, with a wedge attached to the transducer. Currently, the calculation is for two dimensions only.
	%	AnalyticSignal(Data) - to the input data, this method applies a band pass frequency filter, interpolates to a higher sampling frequency, and coverts the input signal to the analytic signal (takes the real signal and makes it complex).
	%	Tfm(Data) - Calculates the total focusing method (TFM) image.
	%	FigPlot = PlotResults(ExistingFigure) - displays the total focusing method (TFM) image.
	%	FigPlot = PlotPaths(ExistingFigure,Transducer) - displays the time for a wave to reach a point in the TFM image for a single transducer element.
	%	ClearCurrentData(ClearImage,ClearDataAnalytic,ClearPaths) - clears data within the class.

	%% MicroPulse5 properties - only get methods available.
	properties (SetAccess = private)
		Image;		% (NOT USER SETTABLE - USE Tfm METHOD) The last total focusing method (TFM) image created.
		DataAnalytic;		% (NOT USER SETTABLE - USE AnalyticSignal METHOD) The input signal after it has been band pass filtered, interpolated to a higher sampling frequency, and converted to the analytic signal form.
		Speed;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The wave propagation speed in m/s.
		DataLength;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The number of time points for a given capture for any of the receivers.
		StartTime;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The time value in s of the first time point in any given capture.
		HighPass;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The high pass filter frequency in Hz (not MHz), which must be greater than zero.
		LowPass;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The low pass filter frequency in Hz (not MHz), which must be greater than HighPass, and less than a quarter of the sampling frequency.
		SamplingFreq;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The sampling frequency in Hz (not MHz) that the data was captured at.
		SamplingFreqInterp;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The sampling frequency in Hz (not MHz) that the data is interpolated to.
		PadLengthFreq;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The size to pad the frequency domain signal to when calculating the analytic signal.
		PadLengthTime;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The size to pad the time domain signal to when calculating the analytic signal.
		FilterStart;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The last point before the band pass region of the signal when filtering it.
		FilterEnd;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The first point after the band pass region of the signal when filtering it.
		LastPositive;		% (NOT USER SETTABLE - USE SetDataParameters METHOD) The last positive point in the FFT when calculating the analytic signal.
		PositionsX;		% (NOT USER SETTABLE - USE SetTransducerPosition METHOD) The X positions of the transducer elements in m.  Note that for a linear array, this is usually all that is necessary.
		PositionsY;		% (NOT USER SETTABLE - USE SetTransducerPosition METHOD) The Y positions of the transducer elements in m.  Note that this is the position in the vertical (depth) direction, and usually will just be zeros, or left blank.
		PositionsZ;		% (NOT USER SETTABLE - USE SetTransducerPosition METHOD) The Z positions of the transducer elements in m.
		Elements;		% (NOT USER SETTABLE - USE SetTransducerPosition METHOD) The total number of transducer elements.
		RegionX;		% (NOT USER SETTABLE - USE SetImageRegion METHOD) The X points that are to be imaged with the total focusing method (TFM).
		RegionY;		% (NOT USER SETTABLE - USE SetImageRegion METHOD) The Y points that are to be imaged with the total focusing method (TFM).
		RegionZ;		% (NOT USER SETTABLE - USE SetImageRegion METHOD) The Z points that are to be imaged with the total focusing method (TFM).  This will not be used for a two dimensional image, and usually will just be zeros, or left blank.
		Paths;		% (NOT USER SETTABLE - USE SetImageRegion METHOD) The distances in m from each transducer to each point that is to be imaged.
		SpeedWedge;		% (NOT USER SETTABLE - USE CalculatePathsWedge METHOD) The wave propagation speed through the wedge in m/s.
		Interface;		% (NOT USER SETTABLE - USE CalculatePathsWedge METHOD) The points describing the wedge/sample interface, currently just in two dimensions, as [XStart YStart; XEnd YEnd].
		EdgeOfWedge;		% (NOT USER SETTABLE - USE CalculatePathsWedge METHOD) If the optimum path for the propagating wave may lie outside the wedge/sample interface for the inspection region chosen.
		RefreshPaths = true;		% (NOT USER SETTABLE) If a parameter change means that the path lengths need to be calculated again.
		RefreshAnalytic = true;		% (NOT USER SETTABLE) If a parameter change means that the analytic signals need to be calculated again.
	end

	%% Public methods (accessible by users of the class)
	methods

		%% Public methods

		function SetDataParameters(obj,DataLength,StartTime,SamplingFreq,HighPass,LowPass,Speed)
			% Set the basic parameters which describe the data to be processed.
			% DataLength is the number of time points for a given capture for any of the receivers.
			% StartTime is the time value in s of the first time point in any given capture.
			% SamplingFreq is the sampling frequency in Hz (not MHz) that the data was captured at.
			% HighPass is the high pass filter frequency in Hz (not MHz), and must be greater than zero.
			% LowPass is the low pass filter frequency in Hz (not MHz), and must be greater than HighPass, and less than a quarter of the sampling frequency.
			% Speed is the wave propagation speed in m/s.
			if (DataLength >= 2) && (SamplingFreq >= 0.0) && (HighPass >= 0.0) && (LowPass > HighPass) && (LowPass < SamplingFreq / 4) && (Speed > 0.0)
				obj.DataLength = DataLength;
				obj.StartTime = StartTime;
				obj.SamplingFreq = SamplingFreq;
				obj.HighPass = HighPass;
				obj.LowPass = LowPass;
				obj.Speed = Speed;
			else
				error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
			end
			TimeScale = 2;
			FreqScale = max(1,pow2(nextpow2(10 * obj.LowPass / obj.SamplingFreq)));
%			FractionError = 0.1;
%			FreqScale = max(1,pow2(nextpow2((2.0 * pi * obj.LowPass / asin(FractionError)) / obj.SamplingFreq)));
			obj.PadLengthTime = pow2(nextpow2(DataLength * TimeScale));
			obj.PadLengthFreq = obj.PadLengthTime * FreqScale;
			obj.SamplingFreqInterp = obj.SamplingFreq * FreqScale;
			obj.FilterStart = floor(obj.HighPass * obj.PadLengthTime / obj.SamplingFreq);
			obj.FilterEnd = ceil(obj.LowPass * obj.PadLengthTime / obj.SamplingFreq + 2);
			obj.LastPositive = obj.PadLengthTime / 2 + 1;
			obj.RefreshPaths = true;
			obj.RefreshAnalytic = true;
		end

		function SetTransducerPosition(obj,PositionsX,PositionsY,PositionsZ)
			% Set the transducer positions, usually relative to the centre of the transducer, but any reference point can be used as long as it is consistent.
			% PositionsX is the X positions of the transducer elements in m.  Note that for a linear array, this is usually all that is necessary.
			% PositionsY is the Y positions of the transducer elements in m.  Note that this is the position in the vertical (depth) direction, and usually will just be zeros, or left blank.
			% PositionsZ is the Z positions of the transducer elements in m.
			if isempty(PositionsX)
				error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
			end
			obj.PositionsX = PositionsX(:);
			obj.PositionsY = [];
			obj.PositionsZ = [];
			if (nargin > 2) && ~isempty(PositionsY)
				if all(size(PositionsX) == size(PositionsY))
					obj.PositionsY = PositionsY(:);
				else
					error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
				end
				if nargin > 3
					if all(size(PositionsX) == size(PositionsZ))
						obj.PositionsZ = PositionsZ(:);
					else
						error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
					end
				end
			end
			obj.Elements = length(obj.PositionsX);
			obj.RefreshPaths = true;
		end

		function SetImageRegion(obj,XStart,XEnd,XIncrement,YStart,YEnd,YIncrement,ZStart,ZEnd,ZIncrement)
			% Set the region that is to be imaged with the total focusing method (TFM).
			% XStart, XEnd, and XIncrement, are the variables that describe the start, end, and increment for the image axes in the X direction in m.
			% YStart, YEnd, and YIncrement, are the variables that describe the start, end, and increment for the image axes in the Y direction in m.
			% ZStart, ZEnd, and ZIncrement, are the variables that describe the start, end, and increment for the image axes in the Z direction in m.  For a linear array, this will usually be left blank.
			obj.RegionX = [];
			obj.RegionY = [];
			obj.RegionZ = [];
			if (XEnd > XStart) && (XIncrement < (XEnd - XStart))
				obj.RegionX(1,:) = XStart:XIncrement:XEnd;
			else
				error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
			end
			if (YEnd > YStart) && (YIncrement < (YEnd - YStart))
				obj.RegionY(:,1) = YStart:YIncrement:YEnd;
			else
				error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
			end
			if nargin > 9
				if (ZEnd > ZStart) && (ZIncrement < (ZEnd - ZStart))
					obj.RegionZ(1,1,:) = ZStart:ZIncrement:ZEnd;
				else
					error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
				end
			end
			obj.RefreshPaths = true;
		end

		function CalculatePaths(obj)
			% Calculates the paths from each transducer to each point that is to be imaged.
			if isempty(obj.Speed) || isempty(obj.StartTime) || isempty(obj.RegionX) || isempty(obj.RegionY) || isempty(obj.PositionsX)
				error('TotalFocusingMethod:OutOfOrder','The methods that must be called before this one have not yet been called.');
			end
			Factor1 = obj.SamplingFreqInterp / obj.Speed;
			Factor2 = obj.SamplingFreqInterp * obj.StartTime - 0.75;		% The factor of 0.75 means that 1.5 (0.75 from each part) does not need to be added when calculating the total path during the TFM.
			if isempty(obj.RegionZ)
				obj.Paths = zeros(length(obj.RegionY),length(obj.RegionX),1,length(obj.PositionsX));
			else
				obj.Paths = zeros(length(obj.RegionY),length(obj.RegionX),length(obj.RegionZ),length(obj.PositionsX));
			end
% 			for Transducer = 1:1:obj.Elements
% 				obj.Paths(:,:,:,Transducer) = sqrt(bsxfun(@plus,bsxfun(@plus,(obj.RegionX - obj.PositionsX(Transducer)) .^ 2,(obj.RegionY - obj.PositionsY(Transducer)) .^2),(obj.RegionZ - obj.PositionsZ(Transducer)) .^2)) .* Factor1 - Factor2;
% 			end
			if isempty(obj.PositionsY)
				for Transducer = 1:1:obj.Elements
					obj.Paths(:,:,1,Transducer) = bsxfun(@plus,(obj.RegionX - obj.PositionsX(Transducer)) .^ 2,obj.RegionY .^ 2);
				end
			else
				for Transducer = 1:1:obj.Elements
					obj.Paths(:,:,1,Transducer) = bsxfun(@plus,(obj.RegionX - obj.PositionsX(Transducer)) .^ 2,(obj.RegionY - obj.PositionsY(Transducer)) .^ 2);
				end
			end
			if ~isempty(obj.PositionsZ) || ~isempty(obj.RegionZ)
				if isempty(obj.PositionsZ)
					for Transducer = 1:1:obj.Elements
						obj.Paths(:,:,:,Transducer) = bsxfun(@plus,obj.Paths(:,:,:,Transducer),obj.RegionZ .^ 2);
					end
				elseif isempty(obj.RegionZ)
					for Transducer = 1:1:obj.Elements
						obj.Paths(:,:,:,Transducer) = bsxfun(@plus,obj.Paths(:,:,:,Transducer),obj.PositionsZ(Transducer) * obj.PositionsZ(Transducer));
					end
				else
					for Transducer = 1:1:obj.Elements
						obj.Paths(:,:,:,Transducer) = bsxfun(@plus,obj.Paths(:,:,:,Transducer),(obj.RegionZ - obj.PositionsZ(Transducer)) .^ 2);
					end
				end
			end
			obj.Paths = sqrt(obj.Paths) .* Factor1 - Factor2;
			if any(any(any(any(obj.Paths > (obj.DataLength * obj.PadLengthFreq / obj.PadLengthTime)))))
				warning('TotalFocusingMethod:SignalTooShort','The received signal is too short for the propagation distance - TFM cannot be performed with the current settings.');
			else
				obj.RefreshPaths = false;
			end
		end

		function CalculatePathsWedge(obj,SpeedWedge,Interface)
			% Calculates the paths from each transducer to each point that is to be imaged, with a wedge attached to the transducer. Currently, the calculation is for two dimensions only.
			% SpeedWedge is the speed of the ultrasound wave in the wedge material.
			% Interface is the coordinates describing the interface, which so far is just [FirstPointX FirstPointY; SecondPointX SecondPointY] as thus assumes a linear interace.
			if isempty(obj.Speed) || isempty(obj.StartTime) || isempty(obj.RegionX) || isempty(obj.RegionY) || isempty(obj.PositionsX)
				error('TotalFocusingMethod:OutOfOrder','The methods that must be called before this one have not yet been called.');
			end
			if isempty(obj.PositionsY) || all(size(Interface,1) ~= [2 2])
				error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
			end
			if ~isempty(obj.PositionsZ) || ~isempty(obj.RegionZ)
				error('TotalFocusingMethod:IncompleteImplementation','A failure has occurred due to the incomplete implementation of the TotalFocusingMethod class.');
			end
			obj.SpeedWedge = SpeedWedge;
			obj.Interface = Interface;
			obj.EdgeOfWedge = false;
			%ErrorMaximum = 10 ^ floor(log10(0.5 * max(obj.Speed,obj.SpeedWedge) / obj.SamplingFreqInterp));
			ErrorMaximum = 0.5 * max(obj.Speed,obj.SpeedWedge) / obj.SamplingFreqInterp;
			InterfacePoints = ceil(sqrt((Interface(2,1) - Interface(1,1)) ^ 2 + (Interface(2,2) - Interface(1,2)) ^ 2) / ErrorMaximum);
			InterfaceX(1,1,:) = (0:1:InterfacePoints) .* ((Interface(2,1) - Interface(1,1)) / InterfacePoints) + Interface(1,1);
			if Interface(2,2) == Interface(1,2)
				InterfaceY = Interface(1,2);
				if (obj.SpeedWedge < obj.Speed)% && (max(obj.PositionsX) < min(obj.RegionX))
					HorizontalInterface = true;
				else
					HorizontalInterface = false;
					warning('TotalFocusingMethod:HorizontalInterface','The imaging conditions are not suitable for the horizontal interface condition. A better implementation might fix this.');
				end
			else
				InterfaceY(1,1,:) = (0:1:InterfacePoints) .* ((Interface(2,2) - Interface(1,2)) / InterfacePoints) + Interface(1,2);
				HorizontalInterface = false;
			end
			obj.Paths = zeros(length(obj.RegionY),length(obj.RegionX),1,length(obj.PositionsX));
			InverseSpeed = 1.0 / obj.Speed;
			InverseSpeedWedge = 1.0 / obj.SpeedWedge;
			for Transducer = 1:1:length(obj.PositionsX)
				if HorizontalInterface
					%LastX = obj.PositionsX(Transducer) + abs(obj.PositionsY(Transducer) * tan(asin(obj.SpeedWedge / obj.Speed)));
					LastXIndex = InterfacePoints + 1;%min(ceil((LastX - Interface(1,1)) / ((Interface(2,1) - Interface(1,1)) / InterfacePoints) + 1),InterfacePoints + 1);
					FirstXIndex = 1;%max(floor((obj.PositionsX(Transducer) - Interface(1,1)) .* (InterfacePoints / (Interface(2,1) - Interface(1,1))) + 1),1);
					LastYIndex = 1;
					FirstYIndex = 1;
				else
					LastXIndex = InterfacePoints + 1;
					FirstXIndex = 1;
				end
				TransducerToInterface = sqrt((InterfaceX(FirstXIndex:1:LastXIndex) - obj.PositionsX(Transducer)) .^ 2 + (InterfaceY(FirstYIndex:1:LastYIndex) - obj.PositionsY(Transducer)) .^ 2) .* InverseSpeedWedge;
				for CurrentX = 1:1:length(obj.RegionX)
					InterfaceToPoint = sqrt(bsxfun(@plus,(InterfaceX(FirstXIndex:1:LastXIndex) - obj.RegionX(CurrentX)) .^ 2,bsxfun(@minus,InterfaceY(FirstYIndex:1:LastYIndex),obj.RegionY) .^ 2)) .* InverseSpeed;
					[obj.Paths(:,CurrentX,1,Transducer), Indices] = min(bsxfun(@plus,TransducerToInterface,InterfaceToPoint),[],3);
					%figure; plot([obj.PositionsX(Transducer) obj.PositionsX(Transducer); InterfaceX(Indices(1) + FirstXIndex - 1) InterfaceX(Indices(end) + FirstXIndex - 1); obj.RegionX(CurrentX) obj.RegionX(CurrentX)],[obj.PositionsY(Transducer) obj.PositionsY(Transducer); InterfaceY InterfaceY; obj.RegionY(1) obj.RegionY(end)]); axis('ij'); axis('equal');
					%Difference = obj.SpeedWedge / obj.Speed - squeeze(abs(obj.PositionsX(Transducer) - InterfaceX(Indices)) ./ sqrt((obj.PositionsX(Transducer) - InterfaceX(Indices)) .^ 2 + obj.PositionsY(Transducer) ^ 2)) ./ (squeeze(abs(obj.RegionX(CurrentX) - InterfaceX(Indices))) ./ sqrt(squeeze(obj.RegionX(CurrentX) - InterfaceX(Indices)) .^ 2 + obj.RegionY .^ 2));		% Check against Snell's law (this is only for an interface that is along the x-axis).
					%plot(Difference .* (100.0 * obj.Speed / obj.SpeedWedge));
					%if (Transducer == 1) && (CurrentX == 1)
					%	figure; plot([obj.PositionsX(Transducer) obj.PositionsX(Transducer); InterfaceX(Indices(1)) InterfaceX(Indices(end)); obj.RegionX(CurrentX) obj.RegionX(CurrentX)],[obj.PositionsY(Transducer) obj.PositionsY(Transducer); InterfaceY(Indices(1)) InterfaceY(Indices(end)); obj.RegionY(1) obj.RegionY(end)]);
					%elseif (Transducer == length(obj.PositionsX)) && (CurrentX == length(obj.RegionX))
					%	hold('on'); line([obj.PositionsX(Transducer) obj.PositionsX(Transducer); InterfaceX(Indices(1)) InterfaceX(Indices(end)); obj.RegionX(CurrentX) obj.RegionX(CurrentX)],[obj.PositionsY(Transducer) obj.PositionsY(Transducer); InterfaceY(Indices(1)) InterfaceY(Indices(end)); obj.RegionY(1) obj.RegionY(end)]); hold('off');
					%end
					if any(Indices == 1) || any(Indices == (InterfacePoints - FirstXIndex + 2))
						obj.EdgeOfWedge = true;
					end
				end
			end
			obj.Paths = obj.Paths .* obj.SamplingFreqInterp - (obj.SamplingFreqInterp * obj.StartTime - 0.75);		% The factor of 0.75 means that 1.5 (0.75 from each part) does not need to be added when calculating the total path during the TFM.
			if obj.EdgeOfWedge
				warning('TotalFocusingMethod:InspectionRegionTooLarge','The optimum path for the propagating wave may lie outside the wedge/sample interface for the inspection region chosen.');
			end
			if any(any(any(any(obj.Paths > (obj.DataLength * obj.PadLengthFreq / obj.PadLengthTime)))))
				warning('TotalFocusingMethod:SignalTooShort','The received signal is too short for the propagation distance - TFM cannot be performed with the current settings.');
			else
				obj.RefreshPaths = false;
			end
		end

		function LoadPaths(obj,Paths)
			% Loads previously calculated paths from each transducer to each point that is to be imaged.
			% Paths are the paths to be loaded - the user must check that the data is suitable.
			obj.Paths = Paths;
			if any(any(any(any(obj.Paths > (obj.DataLength * obj.PadLengthFreq / obj.PadLengthTime)))))
				warning('TotalFocusingMethod:SignalTooShort','The received signal is too short for the propagation distance - TFM cannot be performed with the current settings.');
			else
				obj.RefreshPaths = false;
			end
		end

		function AnalyticSignal(obj,Data)
			% To the input data, this method applies a band pass frequency filter, interpolates to a higher sampling frequency, and coverts the input signal to the analytic signal (takes the real signal and makes it complex).
			% Data is the input data in the form Data(Time,Receiver,Transmitter).  The Receiver and Transmitter dimensions should be the same.
			if size(Data,2) ~= size(Data,3)
				error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
			end
%			MemoryUsageLimit = 1 * 1024 ^ 3;
%			MaxPadLengthFreq = max(obj.PadLengthTime,2 ^ floor(log2(MemoryUsageLimit / (8 * size(Data,2) * size(Data,3)))));
%			obj.PadLengthFreq = min(obj.PadLengthFreq,MaxPadLengthFreq);
			DataFft = fft(Data,obj.PadLengthTime);
			if obj.FilterStart > 0
				DataFft(1:1:obj.FilterStart,:,:) = 0.0;
			else
				DataFft(1,:,:) = DataFft(1,:,:) .* 0.5;
			end
			if obj.FilterEnd > obj.LastPositive
				LastPoint = obj.LastPositive;
				DataFft(obj.LastPositive,:,:) = DataFft(obj.LastPositive,:,:) .* 0.5;
			else
				LastPoint = obj.FilterEnd;
			end
			obj.DataAnalytic = ifft(DataFft(1:1:LastPoint,:,:),obj.PadLengthFreq);
			obj.RefreshAnalytic = false;
		end

		function Tfm(obj,Data)
			% Calculates the total focusing method (TFM) image.
			% Data is the input data in the form Data(Time,Receiver,Transmitter).  The Receiver and Transmitter dimensions should be the same.
			if obj.RefreshPaths
				obj.CalculatePaths();
			end
			if nargin > 1
				obj.AnalyticSignal(Data);
			elseif obj.RefreshAnalytic
				error('TotalFocusingMethod:ParameterMissing','A parameter was missing.');
			end
			obj.DataAnalytic(end,:,:) = 0.0;
			SizePaths = size(obj.Paths);
			obj.Image = zeros(SizePaths(1)*SizePaths(2)*SizePaths(3),1);
			for Transmitter = 1:1:obj.Elements
				for Receiver = Transmitter:1:obj.Elements
					TotalPath = floor(obj.Paths(:,:,:,Transmitter) + obj.Paths(:,:,:,Receiver));
					TotalPath(TotalPath < 1) = obj.PadLengthFreq;
					if Transmitter == Receiver
						obj.Image = obj.Image + obj.DataAnalytic(TotalPath,Receiver,Transmitter);
					else
						obj.Image = obj.Image + obj.DataAnalytic(TotalPath,Receiver,Transmitter) + obj.DataAnalytic(TotalPath,Transmitter,Receiver);
					end
				end
			end
			obj.Image = reshape(obj.Image,SizePaths(1),SizePaths(2),SizePaths(3));
		end

		function FigPlot = PlotResults(obj,ExistingFigure,Normalise)
			% Displays the total focusing method (TFM) image.
			% ExistingFigure is if there is an existing figure that the data should be plotted to, and can left blank if a new figure is to be created by the method.
			if isempty(obj.Image)
				error('TotalFocusingMethod:OutOfOrder','The methods that must be called before this one have not yet been called.');
			end
			ColourMap = 1.0 - gray(256);
			if (nargin > 1) && ~isempty(ExistingFigure)
				FigPlot = ExistingFigure;
			else
				FigPlot = figure('Color','w','Colormap',ColourMap,'InvertHardcopy','off');
			end
			ImageMagnitude = abs(obj.Image);
			if (nargin < 3) || Normalise
				imagesc(obj.RegionX([1 end]) .* 1000.0,obj.RegionY([1 end]) .* 1000.0,ImageMagnitude(:,:,1) .* (1.0 / max(max(max(ImageMagnitude)))));
			else
				imagesc(obj.RegionX([1 end]) .* 1000.0,obj.RegionY([1 end]) .* 1000.0,ImageMagnitude(:,:,1));
			end
			TheXLabel = xlabel('X Position (mm)'); TheYLabel = ylabel('Y Position (mm)'); ColourBar = colorbar; ColourBarLabel = ylabel(ColourBar,'Response (arbitrary)');
			set(get(FigPlot,'CurrentAxes'),'FontName','Arial','FontSize',14);
			set([TheXLabel TheYLabel ColourBarLabel],'FontName','Arial','FontSize',16);
			set(get(FigPlot,'CurrentAxes'),'TickDir','out','TickLength',[0.008 0.008],'XMinorTick','on','YMinorTick','on','XColor',[0.0 0.0 0.0],'YColor',[0.0 0.0 0.0],'LineWidth',0.5,'Box','on');
			set(ColourBar,'TickDir','in','TickLength',[0.008 0.008],'YMinorTick','on','XColor',[0.0 0.0 0.0],'YColor',[0.0 0.0 0.0],'LineWidth',0.5,'Box','on');
		end

		function FigPlot = PlotPaths(obj,ExistingFigure,Transducer)
			% Displays the time for a wave to reach a point in the TFM image for a single transducer element.
			% ExistingFigure is if there is an existing figure that the data should be plotted to, and can left blank if a new figure is to be created by the method.
			% Transducer is the element in the phased array for which the times to the imaging region are to be plotted.
			if isempty(obj.Paths)
				error('TotalFocusingMethod:OutOfOrder','The methods that must be called before this one have not yet been called.');
			end
			ColourMap = 1.0 - gray(256);
			if (nargin > 1) && ~isempty(ExistingFigure)
				FigPlot = ExistingFigure;
			else
				FigPlot = figure('Color','w','Colormap',ColourMap,'InvertHardcopy','off');
			end
			if (nargin > 2) && ~isempty(Transducer)
				if Transducer > obj.Elements
					error('TotalFocusingMethod:ParameterIncorrect','A parameter was incorrect.');
				end
			else
				Transducer = floor(obj.Elements * 0.5 + 0.5);
			end
			PathDistances = (obj.Paths(:,:,floor(size(obj.Paths,3) * 0.5 + 0.5),Transducer) + (obj.SamplingFreqInterp * obj.StartTime - 0.75)) .* (1.0e6 / obj.SamplingFreqInterp);
			imagesc(obj.RegionX([1 end]) .* 1000.0,obj.RegionY([1 end]) .* 1000.0,PathDistances);
			TheXLabel = xlabel('X Position (mm)'); TheYLabel = ylabel('Y Position (mm)'); ColourBar = colorbar; ColourBarLabel = ylabel(ColourBar,'Time of flight ({\mu}s)');
			set(get(FigPlot,'CurrentAxes'),'FontName','Arial','FontSize',14);
			set([TheXLabel TheYLabel ColourBarLabel],'FontName','Arial','FontSize',16);
			set(get(FigPlot,'CurrentAxes'),'TickDir','out','TickLength',[0.008 0.008],'XMinorTick','on','YMinorTick','on','XColor',[0.0 0.0 0.0],'YColor',[0.0 0.0 0.0],'LineWidth',0.5,'Box','on');
			set(ColourBar,'TickDir','in','TickLength',[0.008 0.008],'YMinorTick','on','XColor',[0.0 0.0 0.0],'YColor',[0.0 0.0 0.0],'LineWidth',0.5,'Box','on');
		end

		function ClearCurrentData(obj,ClearImage,ClearDataAnalytic,ClearPaths)
			% Clears data within the class.
			% ClearImage, ClearDataAnalytic, and ClearPaths should be set to true to clear the data for each of those variables.
			if (nargin > 1) && ~isempty(ClearImage) && ClearImage
				obj.Image = [];
			end
			if (nargin > 2) && ~isempty(ClearDataAnalytic) && ClearDataAnalytic
				obj.DataAnalytic = [];
			end
			if (nargin > 3) && ~isempty(ClearPaths) && ClearPaths
				obj.Paths = [];
			end
		end

	end

end