diff --git a/FCWTAPI/GaussianSmoothing.cs b/FCWTAPI/GaussianSmoothing.cs index 19fac1b..8f7c1b1 100644 --- a/FCWTAPI/GaussianSmoothing.cs +++ b/FCWTAPI/GaussianSmoothing.cs @@ -105,6 +105,107 @@ public class GaussianSmoothing } return dim2Smoothing; } + /// + /// Performs 2D Gaussian smoothing using an elliptic Gaussian kernel + /// + /// + /// Deviation of the convolving Gaussian in the frequency axis + /// Deviation of the convolving Gaussian in the time axis + /// + /// + public static double[,] EllipticGaussianConvolution(double[,] inputData, double frequencyDeviation, double timeDeviation) + { + if (inputData.GetLength(0) < (frequencyDeviation * 6 + 1) || inputData.GetLength(1) < (timeDeviation * 6 + 1)) + { + throw new ArgumentException("Matrix may not be smaller than the convoluting Gaussian kernel"); + } + double[,] frequencyKernel = CalculateNormalized1DSampleKernel(frequencyDeviation); + double[,] timeKernel = CalculateNormalized1DSampleKernel(timeDeviation); + // Calculates the 1D kernel to be used for smoothing + double[,] frequencySmoothing = new double[inputData.GetLength(0), inputData.GetLength(1)]; + double[,] timeSmoothing = new double[inputData.GetLength(0), inputData.GetLength(1)]; + // x-direction smoothing + for (int i = 0; i < inputData.GetLength(0); i++) + { + // Calculates each point in the x direction + for (int j = 0; j < inputData.GetLength(1); j++) + frequencySmoothing[i, j] = ProcessPoint(inputData, i, j, frequencyKernel, 0); + } + //y-direction smoothing + for (int i = 0; i < inputData.GetLength(0); i++) + { + //Calculates each point in the y direction from the x smoothed array res 1 + for (int j = 0; j < inputData.GetLength(1); j++) + timeSmoothing[i, j] = ProcessPoint(frequencySmoothing, i, j, timeKernel, 1); + } + return timeSmoothing; + } + /// + /// Performs Elliptic Gaussian smoothing on individual rows or columns of an array + /// + /// + /// + /// Deviation of the convolving Gaussian in the frequency axis + /// Deviation of the convolving Gaussian in the time axis + /// + /// + public static double[] SliceEllipticGaussianConvolution(double[,] inputData, double frequencyDeviation, double timeDeviation, int sliceIndex, int dimension) + { + if (inputData.GetLength(0) < (frequencyDeviation * 6 + 1) || inputData.GetLength(1) < (timeDeviation * 6 + 1)) + { + throw new ArgumentException("Matrix may not be smaller than the convoluting Gaussian kernel"); + } + if (dimension != 0 && dimension != 1) + { + throw new ArgumentException("Dimension can only be 0 and 1", nameof(dimension)); + } + double[,] frequencyKernel = CalculateNormalized1DSampleKernel(frequencyDeviation); + double[,] timeKernel = CalculateNormalized1DSampleKernel(timeDeviation); + double[,] sliceDirectionSmoothing; + double[] smoothedSlice; + // Calculates the 1D kernel to be used for smoothing + if (dimension == 0) + { + // Calculates each point in the x direction + int freqSize = (int)Math.Ceiling(frequencyDeviation * 3) * 2 + 1; + sliceDirectionSmoothing = new double[freqSize, inputData.GetLength(1)]; + smoothedSlice = new double[inputData.GetLength(1)]; + for (int i = 0; i < freqSize; i++) + { + int dataIndex = sliceIndex + i - freqSize / 2; + for (int j = 0; j < inputData.GetLength(1); j++) + { + sliceDirectionSmoothing[i, j] = ProcessPoint(inputData, dataIndex, j, frequencyKernel, 0); + } + } + for (int j = 0; j < inputData.GetLength(1); j++) + { + smoothedSlice[j] = ProcessPoint(sliceDirectionSmoothing, freqSize / 2, j, timeKernel, 1); + } + return smoothedSlice; + } + else + { + // Calculates each point in the x direction + int timeSize = (int)Math.Ceiling(timeDeviation * 3) * 2 + 1; + sliceDirectionSmoothing = new double[inputData.GetLength(0), timeSize]; + smoothedSlice = new double[inputData.GetLength(0)]; + for (int i = 0; i < inputData.GetLength(0); i++) + { + for (int j = 0; j < timeSize; j++) + { + int dataIndex = sliceIndex + j - timeSize / 2; + sliceDirectionSmoothing[i, j] = ProcessPoint(inputData, i, dataIndex, frequencyKernel, 0); + } + } + for (int i = 0; i < inputData.GetLength(0); i++) + { + smoothedSlice[i] = ProcessPoint(sliceDirectionSmoothing, i, timeSize / 2, timeKernel, 1); + } + return smoothedSlice; + } + } + /// /// Method to process each individual point of a given double[,] array /// diff --git a/FCWTAPI/MinMaxIdentification.cs b/FCWTAPI/MinMaxIdentification.cs index d4d4970..8859e2c 100644 --- a/FCWTAPI/MinMaxIdentification.cs +++ b/FCWTAPI/MinMaxIdentification.cs @@ -3,11 +3,258 @@ using System.Linq; using System.Text; using System.Threading.Tasks; +using FCWTNET; -namespace fCWT.NET +namespace FCWTNET { public class MinMaxIdentification { + // This constructor is a WIP I'll fix it when I start finding relative maxima + public double[,] InputData; + public double FrequencySmoothingWidth; + public double TimeSmoothingWidth; + //public double MinFrequencyWidth; + //public int? RowSpacings; + public int TimeDerivativeDistance; + //public double IntensityThreshold; + public string WorkingPath; + public MinMaxIdentification(double[,] inputData, int timeDerivativeDistance, string workingPath) + { + InputData = inputData; + TimeDerivativeDistance = timeDerivativeDistance; + WorkingPath = workingPath; + } + /// + /// Method to take derivatives of a 1D array with Δx based on the index + /// This method exists entirely for testing purposes + /// + /// Input 1D array + /// Distance on each side of the target point to select initial and final points + /// + public static double[] StandardArrayDerivative(double[] inputSlice, int derivativeDistance) + { + if (inputSlice.Length / 2 < derivativeDistance) + { + throw new ArgumentException("derivativeDistance must be at least less than half the length of inputSlice", nameof(derivativeDistance)); + } + double[] outputArray = new double[inputSlice.Length]; + for (int i = 0; i < inputSlice.Length; i++) + { + double intensityDifference; + double axisDifference; + // Edge case where the target point is closer than 1 derivativeDistance to the left edge + if ( i - derivativeDistance < 0) + { + intensityDifference = inputSlice[i + derivativeDistance] - inputSlice[0]; + axisDifference = i + derivativeDistance; + } + // Edge case where the target point is closer than 1 derivativeDistance to the right edge + else if (i + derivativeDistance >= inputSlice.Length) + { + intensityDifference = inputSlice[^1] - inputSlice[i - derivativeDistance]; + axisDifference = inputSlice.Length - i + derivativeDistance - 1; + } + // Standard case + else + { + intensityDifference = inputSlice[i + derivativeDistance] - inputSlice[i - derivativeDistance]; + axisDifference = 2 * derivativeDistance; + } + outputArray[i] = intensityDifference / axisDifference; + } + return outputArray; + } + /// + /// Method for taking a downsampled derivative of data in the SortedDictionary format + /// This takes a derivative at 1 point for every 2 * derivativeDistance Points + /// Do not use this method multiple times in a row to take higher derivatives, + /// use it once if downsampling is required, then use StandardDerivative + /// This method is likely not parallelizable + /// + /// + /// Distance on each side of the target point to select initial and final points + /// + /// + public static SortedDictionary DownsampledSliceDerivative(SortedDictionary filteredData, int derivativeDistance) + { + int prevIndex = -1; + int counter = -1; + SortedDictionary derivativePoints = new SortedDictionary(); + foreach (var point in filteredData) + { + if (prevIndex == -1 || prevIndex == point.Key - 1) + { + counter++; + prevIndex = point.Key; + if (counter % (2 * derivativeDistance + 1) == 2 * derivativeDistance) + { + double startValue = filteredData[point.Key - (2 * derivativeDistance)]; + double endValue = point.Value; + double derivativeValue = (endValue - startValue) / (2 * derivativeDistance); + int derivativeIndex = point.Key - derivativeDistance; + derivativePoints.Add(derivativeIndex, derivativeValue); + } + } + else + { + prevIndex = -1; + counter = 0; + } + } + if (derivativePoints.Count == 0) + { + throw new Exception("No derivative points could be calculated, derivativeDistance may be too large relative to the data"); + } + return derivativePoints; + } + /// + /// Method to take a derivative of data in the SortedDictionary format without downsampling. + /// This loses a single point at the edge of each coherent packet which may need to be fixed + /// however, it won't cause problems for the application of peak identification. + /// + /// + /// If the data was previously downsampled during differentiation, the derivativeDistance used must be given here + /// + public static SortedDictionary StandardDerivative(SortedDictionary inputData, int? originalDerivativeDistance = null) + { + int prevIndex = -1; + SortedDictionary derivativePoints = new SortedDictionary(); + int indexSpacing; + KeyValuePair oneEarlierPoint = new KeyValuePair(-1, 0); + KeyValuePair twoEarlierPoint = new KeyValuePair(-1, 0); + if (!originalDerivativeDistance.HasValue) + { + indexSpacing = 1; + } + else + { + indexSpacing = (originalDerivativeDistance.Value * 2) + 1; + } + foreach(var point in inputData) + { + // Checks that we're in a packet of evenly spaced data with the anticipated spacing between indices + if (prevIndex == -1 || prevIndex == point.Key - indexSpacing) + { + prevIndex = point.Key; + // Checks that there is a point to the left and right of oneEarlierPoint + // the point which we are taking the derivative at. + if (oneEarlierPoint.Key != -1 && twoEarlierPoint.Key != -1) + { + double derivativeValue = (point.Value - twoEarlierPoint.Value) / (point.Key - twoEarlierPoint.Key); + int derivativeIndex = oneEarlierPoint.Key; + derivativePoints.Add(derivativeIndex, derivativeValue); + } + // Updates oneEarlierPoint and twoEarlierPoint + twoEarlierPoint = oneEarlierPoint; + oneEarlierPoint = point; + + } + // If we're not in a packet of evenly spaced data, resets to begin identifying a new packet + else + { + prevIndex = -1; + oneEarlierPoint = point; + twoEarlierPoint = new KeyValuePair(-1, 0); + } + } + if (derivativePoints.Count == 0) + { + throw new ArgumentException("No regions of properly spaced data were identified, originalDerivativeDistance is likely incorrect", nameof(originalDerivativeDistance)); + } + return derivativePoints; + + } + /// + /// Smooths a downsampledDerivative using a smoothing deviation scaled with that downsampling + /// This method is designed to work with downsampled and intensity filtered data + /// If a sequential region does not have enough points for the Gaussian smoothing, it is not smoothed and is left alone. + /// + /// Input downsampled derivative + /// Deviation of the Gaussian kernel used to smooth in terms of the original size of the data + /// Derivative distance used when calculating the downsampled derivative + /// + public static SortedDictionary DownsampledDerivativeSmoothing(SortedDictionary downsampledDerivative, double derivativeSmoothingDeviation, int derivativeDistance) + { + int prevIndex = -1; + SortedDictionary smoothedDerivativePoints = new SortedDictionary(); + List clusterIndexList = new List(); + List clusterValueList = new List(); + bool smoothingAppliedOnce = false; + foreach(var point in downsampledDerivative) + { + if ((prevIndex == -1 || prevIndex == point.Key - (2 * derivativeDistance + 1)) && point.Key != downsampledDerivative.Keys.Last()) + { + clusterIndexList.Add(point.Key); + clusterValueList.Add(point.Value); + prevIndex = point.Key; + } + else + { + if (point.Key == downsampledDerivative.Keys.Last()) + { + clusterIndexList.Add(point.Key); + clusterValueList.Add(point.Value); + } + if(clusterIndexList.Count > (6 * (derivativeSmoothingDeviation/ (1 + 2 * (double)derivativeDistance))) + 1) + { + double[] derivativeArray = clusterValueList.ToArray(); + double[] smoothedDerivative = GaussianSmoothing.GaussianSmoothing1D(derivativeArray, derivativeSmoothingDeviation/(1 + 2 * (double)derivativeDistance)); + for (int i = 0; i < smoothedDerivative.Length; i++) + { + smoothedDerivativePoints.Add(clusterIndexList[i], smoothedDerivative[i]); + } + if (!smoothingAppliedOnce) + { + smoothingAppliedOnce = true; + } + } + else + { + for(int i = 0; i < clusterIndexList.Count; i++) + { + smoothedDerivativePoints.Add(clusterIndexList[i], clusterValueList[i]); + } + } + clusterIndexList = new List(); + clusterValueList = new List(); + clusterIndexList.Add(point.Key); + clusterValueList.Add(point.Value); + prevIndex = -1; + } + } + if (smoothedDerivativePoints.Count == 0) + { + throw new ArgumentException("No regions of properly spaced data were identified, derivativeDistance is likely incorrect", nameof(derivativeDistance)); + } + if (!smoothingAppliedOnce) + { + throw new ArgumentException("No sequential regions exist large enough for smoothing to be applied. Smoothing Kernel is too big", nameof(derivativeSmoothingDeviation)); + } + return smoothedDerivativePoints; + } + public int[] FrequencySliceMaxIdentification() + { + throw new NotImplementedException(); + } + public static SortedDictionary IndexLinkedIntensityFiltering(double[] inputSlice, double intensityThreshold) + { + if(intensityThreshold > 1 || intensityThreshold < 0) + { + throw new ArgumentException("intensityThreshold must be greater than 0 and less than or equal to 1", nameof(intensityThreshold)); + } + double maxValue = inputSlice.Max(); + double minValue = inputSlice.Min(); + double minIntensity = Math.Abs(maxValue) > Math.Abs(minValue) ? Math.Abs(maxValue) * intensityThreshold: Math.Abs(minValue) * intensityThreshold; + SortedDictionary intensityFilteredData = new SortedDictionary(); + for (int i = 0; i < inputSlice.Length; i++) + { + if(Math.Abs(inputSlice[i]) > minIntensity) + { + intensityFilteredData.Add(i, inputSlice[i]); + } + } + return intensityFilteredData; + } } } diff --git a/FCWTAPI/PlottingUtils.cs b/FCWTAPI/PlottingUtils.cs index b2589d7..962a5ca 100644 --- a/FCWTAPI/PlottingUtils.cs +++ b/FCWTAPI/PlottingUtils.cs @@ -58,8 +58,16 @@ public static PlotModel GenerateCWTHeatMap(double[,] data, string plotTitle, dou } public static PlotModel GenerateCWTHeatMap(double[,] data, string plotTitle) { - double[] defaultTimeAxis = new double[data.GetLength(1)]; - double[] defaultFrequencyAxis = new double[data.GetLength(0)]; + double[] defaultTimeAxis = new double[data.GetLength(0)]; + double[] defaultFrequencyAxis = new double[data.GetLength(1)]; + for (int i = 0; i < data.GetLength(1); i++) + { + defaultFrequencyAxis[i] = Math.Pow(2, 1 + (0.001* i)); + } + for (int i = 0; i < data.GetLength(0); i++) + { + defaultTimeAxis[i] = i; + } return GenerateCWTHeatMap(data, plotTitle, defaultTimeAxis, defaultFrequencyAxis); } /// @@ -101,11 +109,18 @@ public static PlotModel GenerateCWTContourPlot(double[,] data, string plotTitle, } public static PlotModel GenerateCWTContourPlot(double[,] data, string plotTitle) { - double[] defaultTimeAxis = new double[data.GetLength(1)]; - double[] defaultFrequencyAxis = new double[data.GetLength(0)]; + double[] defaultTimeAxis = new double[data.GetLength(0)]; + double[] defaultFrequencyAxis = new double[data.GetLength(1)]; + for (int i = 0; i < data.GetLength(1); i++) + { + defaultFrequencyAxis[i] = Math.Pow(2, 1 + (0.001 * i)); + } + for (int i = 0; i < data.GetLength(0); i++) + { + defaultTimeAxis[i] = i; + } return GenerateCWTContourPlot(data, plotTitle, defaultTimeAxis, defaultFrequencyAxis); } - public static PlotModel GenerateXYPlotCWT(double[,] data, int[] rowIndices, double[] timeAxis, double[] freqAxis, PlotTitles plotTitle, XYPlotOptions mode, string? customTitle = null) { string actualTitle; @@ -211,10 +226,86 @@ public static PlotModel GenerateXYPlotCWT(double[,] data, int[] rowIndices, doub } public static PlotModel GenerateXYPlotCWT(double[,] data, int[] rowIndices, PlotTitles plotTitle, XYPlotOptions mode, string? customTitle = null) { - double[] defaultTimeAxis = new double[data.GetLength(1)]; - double[] defaultFrequencyAxis = new double[data.GetLength(0)]; - return GenerateXYPlotCWT(data, rowIndices, defaultTimeAxis, defaultFrequencyAxis, plotTitle, mode); + double[] defaultTimeAxis = new double[data.GetLength(0)]; + double[] defaultFrequencyAxis = new double[data.GetLength(1)]; + for (int i = 0; i < data.GetLength(1); i++) + { + defaultFrequencyAxis[i] = Math.Pow(2, 1 + (0.001 * i)); + } + for (int i = 0; i < data.GetLength(0); i++) + { + defaultTimeAxis[i] = i; + } + return GenerateXYPlotCWT(data, rowIndices, defaultTimeAxis, defaultFrequencyAxis, plotTitle, mode, customTitle); + } + public static PlotModel Plot1DArray(double[] data) + { + var plotModel = new PlotModel(); + plotModel.Axes.Add(new LinearAxis + { + Position = AxisPosition.Bottom, + MinimumPadding = 0.05, + MaximumPadding = 0.05, + Title = "Index", + FontSize = 14, + TitleFontSize = 16, + AxisTitleDistance = 10, + Minimum = 0, + Maximum = data.Length - 1, + }); + plotModel.Axes.Add(new LinearAxis + { + Position = AxisPosition.Left, + MinimumPadding = 0.05, + MaximumPadding = 0.05, + Title = "Intensity", + FontSize = 14, + TitleFontSize = 16, + AxisTitleDistance = 10 + }); + LineSeries dataSeries = new LineSeries(); + for (int i = 0; i < data.Length; i++) + { + dataSeries.Points.Add(new DataPoint(i, data[i])); + } + plotModel.Series.Add(dataSeries); + return plotModel; + + } + public static PlotModel PlotSortedPointDictionary(SortedDictionary data) + { + var plotModel = new PlotModel(); + plotModel.Axes.Add(new LinearAxis + { + Position = AxisPosition.Bottom, + MinimumPadding = 0.05, + MaximumPadding = 0.05, + Title = "Index", + FontSize = 14, + TitleFontSize = 16, + AxisTitleDistance = 10, + Minimum = data.Keys.First(), + Maximum = data.Keys.Last() - 1, + }); + plotModel.Axes.Add(new LinearAxis + { + Position = AxisPosition.Left, + MinimumPadding = 0.05, + MaximumPadding = 0.05, + Title = "Intensity", + FontSize = 14, + TitleFontSize = 16, + AxisTitleDistance = 10 + }); + LineSeries dataSeries = new LineSeries(); + foreach(var point in data) + { + dataSeries.Points.Add(new DataPoint(point.Key, point.Value)); + } + plotModel.Series.Add(dataSeries); + return plotModel; } + /// /// Exports generated plots to PDF files /// Expansion to this method should include options for other file formats @@ -224,6 +315,7 @@ public static PlotModel GenerateXYPlotCWT(double[,] data, int[] rowIndices, Plot /// Width of the plot /// Height of the plot /// + public static void ExportPlotPDF(PlotModel plotModel, string filePath, int plotWidth = 700, int plotHeight = 600) { if (Path.GetExtension(filePath) == null) diff --git a/TestFCWTAPI/GaussianSmoothingTests.cs b/TestFCWTAPI/GaussianSmoothingTests.cs index 7523609..fb36aa7 100644 --- a/TestFCWTAPI/GaussianSmoothingTests.cs +++ b/TestFCWTAPI/GaussianSmoothingTests.cs @@ -119,7 +119,7 @@ public void TestProcessPoint() } } [Test] - public void TestGaussianConvolution() + public static void TestGaussianConvolution() { double[,] invalid2DArray = new double[,] { @@ -169,11 +169,11 @@ public void TestGaussianConvolution() // Deviation of our test 2D Gaussian double testDeviation = 1; // Calculates expected value of a point at the center of the 2D Gaussian - double centerPoint = 1 / (2 * Math.PI * testDeviation); + double centerPoint = 1 / (2 * Math.PI * Math.Pow(testDeviation, 2)); // Calculates expected value of a point 1 unit off from the center of the 2D Gaussian in any direction - double adjacentPoint = 1 / (2 * Math.PI * testDeviation) * Math.Exp(-(1 * 1) / (2 * testDeviation * testDeviation)); + double adjacentPoint = 1 / (2 * Math.PI * Math.Pow(testDeviation, 2)) * Math.Exp(-(1 * 1) / (2 * testDeviation * testDeviation)); // Calculates expected value of a point 1 unit off from the center in in x and y direction - double diagonalPoint = 1 / (2 * Math.PI * testDeviation) * Math.Exp(-2 / (2 * testDeviation * testDeviation)); + double diagonalPoint = 1 / (2 * Math.PI * Math.Pow(testDeviation, 2)) * Math.Exp(-2 / (2 * testDeviation * testDeviation)); // Calculats the Gaussian convolution of our test 2D Gaussian with tesd2DArray double[,] testBlurredArray = GaussianSmoothing.GaussianConvolution(test2DArray, testDeviation); for (int p = 0; p < 4; p++) @@ -191,6 +191,154 @@ public void TestGaussianConvolution() * fixed later. The unit test for this new function should include a module to ensure * that the sum of the original array equals the sum of the blurred array */ + } + [Test] + public static void TestEllipticGaussianConvolution() + { + double[,] invalid2DArray = new double[,] + { + {1, 2, 3, 4, 5 }, + {4, 5, 6, 7, 8 }, + {9, 10, 11, 12, 13 }, + }; + Assert.Throws(() => GaussianSmoothing.EllipticGaussianConvolution(invalid2DArray, 1, 2)); + double testDim1Deviation = 1; + double testDim2Deviation = 2; + double[,] test2DArray = new double[51, 51]; + // List of points to add to a test array + var pointArray = new (int, int, double)[] + { + (13, 13, 1), // Generic point + (27, 13, -1), // Sample point with a negative value + (2, 43, 1), // Sample point next to an edge + (50, 50, 1) // Sample point in a corner + }; + double centerPoint = 1 / (2 * Math.PI * testDim1Deviation * testDim2Deviation); + // Generates a 2D array containing all points from pointArray + for (int i = 0; i < 51; i++) + { + for (int j = 0; j < 51; j++) + { + for (int t = 0; t < 4; t++) + if (i == pointArray[t].Item1 && j == pointArray[t].Item2) + { + test2DArray[i, j] = pointArray[t].Item3; + } + } + } + + // Generates a list of coordinates to points adjacent to the points in pointArray off by 1 unit in the higher dimension + var adjacentDim1Array = new (int, int, double)[] + { + (12, 13, 1), + (28, 13, -1), + (3, 43, 1), + (49, 50, 1) + }; + double dim1AdjacentPoint = 1 / (2 * Math.PI * testDim1Deviation * testDim2Deviation) * Math.Exp(-1 / (2 * Math.Pow(testDim1Deviation, 2))); + + // Generates a list of coordinates to points adjacent to the points in pointArray off by 1 unit in the lower dimension + var adjacentDim2Array = new (int, int, double)[] + { + (13, 14, 1), + (27, 12, -1), + (2, 42, 1), + (50, 49, 1) + + }; + double dim2AdjacentPoint = 1 / (2 * Math.PI * testDim1Deviation * testDim2Deviation) * Math.Exp(-1 / (2 * Math.Pow(testDim2Deviation, 2))); + + // Generates a list of coordinates to points diagonal to the points in pointArray + var diagonalArray = new (int, int, double)[] + { + (14, 14, 1), + (26, 14, -1), + (1, 42, 1), + (49, 49, 1) + }; + double diagonalPoint = 1 / (2 * Math.PI * testDim1Deviation * testDim2Deviation) * Math.Exp(-1 * ((1 / (2 * Math.Pow(testDim1Deviation, 2))) + (1 / (2* Math.Pow(testDim2Deviation, 2))))); + double[,] testAsymmetricalBlurArray = GaussianSmoothing.EllipticGaussianConvolution(test2DArray, testDim1Deviation, testDim2Deviation); + for(int p = 0; p < 4; p++) + { + // Checks all the center point values for correctness + Assert.AreEqual(centerPoint * pointArray[p].Item3, testAsymmetricalBlurArray[pointArray[p].Item1, pointArray[p].Item2], 0.001); + // Checks all the adjacent dim1 point values for correctness + Assert.AreEqual(dim1AdjacentPoint * pointArray[p].Item3, testAsymmetricalBlurArray[adjacentDim1Array[p].Item1, adjacentDim1Array[p].Item2], 0.001); + // Checks all the adjacent dim2 point values for correctness + Assert.AreEqual(dim2AdjacentPoint * pointArray[p].Item3, testAsymmetricalBlurArray[adjacentDim2Array[p].Item1, adjacentDim2Array[p].Item2], 0.001); + // Checks all the diagonal point values for correctness + Assert.AreEqual(diagonalPoint * pointArray[p].Item3, testAsymmetricalBlurArray[diagonalArray[p].Item1, diagonalArray[p].Item2], 0.001); + } + } + + [Test] + public static void TestSliceEllipticGaussianConvolution() + { + double[,] invalid2DArray = new double[,] + { + {1, 2, 3, 4, 5 }, + {4, 5, 6, 7, 8 }, + {9, 10, 11, 12, 13 }, + }; + Assert.Throws(() => GaussianSmoothing.SliceEllipticGaussianConvolution(invalid2DArray, 1, 2, 1, 0)); + double testDim1Deviation = 1; + double testDim2Deviation = 2; + double[,] test2DArray = new double[51, 51]; + // List of points to add to a test array + var pointArray = new (int, int, double)[] + { + (20, 20, 1), // Generic middle point + (20, 2, 1), // Sample point next to a dim 2 edge + (21, 46, 1), // Off slice point for dim 1 smoothing + (3, 46, 1), // Off slice point next to a dim 1 edge for dim 1 smoothing + (2, 20, 1), // Sample point next to a dim 1 edge + (46, 19, 1), // Off slice point for dim 2 smoothing + (46, 3, 1) // off slice point next to a dim 2 edge for dim 2 smoothing + }; + double centerPoint = 1 / (2 * Math.PI * testDim1Deviation * testDim2Deviation); + double dim1AdjacentPoint = 1 / (2 * Math.PI * testDim1Deviation * testDim2Deviation) * Math.Exp(-1 / (2 * Math.Pow(testDim1Deviation, 2))); + double dim2AdjacentPoint = 1 / (2 * Math.PI * testDim1Deviation * testDim2Deviation) * Math.Exp(-1 / (2 * Math.Pow(testDim2Deviation, 2))); + double diagonalPoint = 1 / (2 * Math.PI * testDim1Deviation * testDim2Deviation) * Math.Exp(-1 * ((1 / (2 * Math.Pow(testDim1Deviation, 2))) + (1 / (2 * Math.Pow(testDim2Deviation, 2))))); + // Generates a 2D array containing all points from pointArray + for (int i = 0; i < 51; i++) + { + for (int j = 0; j < 51; j++) + { + for (int t = 0; t < pointArray.Length; t++) + if (i == pointArray[t].Item1 && j == pointArray[t].Item2) + { + test2DArray[i, j] = pointArray[t].Item3; + } + } + } + double[] centerDim1Smoothing = GaussianSmoothing.SliceEllipticGaussianConvolution(test2DArray, testDim1Deviation, testDim2Deviation, 20, 0); + Assert.AreEqual(centerPoint, centerDim1Smoothing[2], 0.001); + Assert.AreEqual(centerPoint, centerDim1Smoothing[20], 0.001); + Assert.AreEqual(dim2AdjacentPoint, centerDim1Smoothing[3], 0.001); + Assert.AreEqual(dim2AdjacentPoint, centerDim1Smoothing[19], 0.001); + Assert.AreEqual(dim1AdjacentPoint, centerDim1Smoothing[46], 0.001); + Assert.AreEqual(diagonalPoint, centerDim1Smoothing[45], 0.001); + + double[] edgeDim1Smoothing = GaussianSmoothing.SliceEllipticGaussianConvolution(test2DArray, testDim1Deviation, testDim2Deviation, 2, 0); + Assert.AreEqual(centerPoint, edgeDim1Smoothing[20], 0.001); + Assert.AreEqual(dim2AdjacentPoint, edgeDim1Smoothing[21], 0.001); + Assert.AreEqual(dim1AdjacentPoint, edgeDim1Smoothing[46], 0.001); + Assert.AreEqual(diagonalPoint, edgeDim1Smoothing[45], 0.001); + + + double[] centerDim2Smoothing = GaussianSmoothing.SliceEllipticGaussianConvolution(test2DArray, testDim1Deviation, testDim2Deviation, 20, 1); + Assert.AreEqual(centerPoint, centerDim2Smoothing[2], 0.001); + Assert.AreEqual(centerPoint, centerDim2Smoothing[20], 0.001); + Assert.AreEqual(dim1AdjacentPoint, centerDim2Smoothing[3], 0.001); + Assert.AreEqual(dim1AdjacentPoint, centerDim2Smoothing[21], 0.001); + Assert.AreEqual(dim2AdjacentPoint, centerDim2Smoothing[46], 0.001); + Assert.AreEqual(diagonalPoint, centerDim2Smoothing[45], 0.001); + + double[] edgeDim2Smoothing = GaussianSmoothing.SliceEllipticGaussianConvolution(test2DArray, testDim1Deviation, testDim2Deviation, 2, 1); + Assert.AreEqual(centerPoint, edgeDim2Smoothing[20], 0.001); + Assert.AreEqual(dim1AdjacentPoint, edgeDim2Smoothing[21], 0.001); + Assert.AreEqual(dim2AdjacentPoint, edgeDim2Smoothing[46], 0.001); + Assert.AreEqual(diagonalPoint, edgeDim2Smoothing[45], 0.001); } [Test] @@ -206,13 +354,13 @@ public static void TestGaussianSmoothing1D() double[] smoothedArray = GaussianSmoothing.GaussianSmoothing1D(testArray, 1); double centerPoint = 1 / Math.Sqrt(2 * Math.PI * testDeviation); double adjacentPoint = 1 / Math.Sqrt(2 * Math.PI * testDeviation) * Math.Exp(-(1 * 1) / (2 * testDeviation * testDeviation)); - double sum = 0; Assert.AreEqual(centerPoint, smoothedArray[2], 0.001); Assert.AreEqual(centerPoint, smoothedArray[20], 0.001); Assert.AreEqual(centerPoint, smoothedArray[49], 0.001); Assert.AreEqual(adjacentPoint, smoothedArray[1], 0.001); Assert.AreEqual(adjacentPoint, smoothedArray[21], 0.001); Assert.AreEqual(adjacentPoint, smoothedArray[48], 0.001); + } } } \ No newline at end of file diff --git a/TestFCWTAPI/MinMaxIdentificationFullDataTests.cs b/TestFCWTAPI/MinMaxIdentificationFullDataTests.cs new file mode 100644 index 0000000..274ccd1 --- /dev/null +++ b/TestFCWTAPI/MinMaxIdentificationFullDataTests.cs @@ -0,0 +1,176 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using FCWTNET; +using System.IO; +using NUnit.Framework; + +namespace TestFCWTAPI +{ + public class MinMaxIdentificationFullDataTests + { + public CWTObject MyoglobinTestCWT { get; set; } + public CWTObject ComplexTestCWT { get; set; } + [OneTimeSetUp] + public void TestDataSetup() + { + // Create test data from real transients + string myoZ13TransientPath = Path.Combine(TestContext.CurrentContext.TestDirectory, "TestDataFiles", "SampleTransient_Myo_z13.csv"); + TransientData myoZ13Transient = new(myoZ13TransientPath); + myoZ13Transient.ImportTransientData(800000); + string testPath = Path.Combine(TestContext.CurrentContext.TestDirectory, "TestDataFiles"); + var myoZ13CWT = new CWTObject(myoZ13Transient.Data, 2, 3, 125, (float)(2 * Math.Tau), 4, false, myoZ13Transient.SamplingRate, myoZ13Transient.CalibrationFactor, testPath); + myoZ13CWT.PerformCWT(); + myoZ13CWT.CalculateTimeAxis(); + myoZ13CWT.CalculateFrequencyAxis(); + MyoglobinTestCWT = myoZ13CWT; + // Plot modulus and real heatmap for basic visualization + string complexTransientPath = Path.Combine(TestContext.CurrentContext.TestDirectory, "TestDataFiles", "ComplexTransient_protein_mixture_512ms.csv"); + TransientData complexTransient = new(complexTransientPath); + complexTransient.ImportTransientData(700000); + var complexTransientCWT = new CWTObject(complexTransient.Data, 2, 3, 125, (float)(15 * Math.Tau), 4, false, complexTransient.SamplingRate, complexTransient.CalibrationFactor, testPath); + complexTransientCWT.PerformCWT(); + complexTransientCWT.CalculateTimeAxis(); + complexTransientCWT.CalculateFrequencyAxis(); + ComplexTestCWT = complexTransientCWT; + } + [Ignore("Long test time")] + [Test] + // Test method for generating plots to carefully examine the CWT of the sample MS transient data + public void PlotGeneration() + { + MyoglobinTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "myoZ13_modulus.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "myo z13"); + MyoglobinTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "myoZ13_modulus_Zoomed.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "myo z13", 0.075, 0.095, 2.1, 3); + MyoglobinTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "myoZ13_modulus_HiZoom.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "myo z13", 0.082, 0.083, 2.4, 2.5); + MyoglobinTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Real, "myoZ13_Real.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "myo z13"); + ComplexTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "complexTransient_modulus.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "complexTransient"); + ComplexTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "zoomed_Band2complexTransient_modulus.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "complexTransient", 0.06, 0.08, 19, 30); + ComplexTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "zoomed_Band1complexTransient_modulus.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "complexTransient", 0.02, 0.04, 19, 30); + } + [Ignore("Long test time")] + [Test] + // Test method to play around with smoothing + public void SmoothingTests() + { + double[,] data = ComplexTestCWT.OutputCWT.ModulusCalculation(); + double[,] timeWindowedData; + double[] dummyTimeAx; + CWTExtensions.TimeWindowing(0.05, 0.09, ComplexTestCWT.TimeAxis, data, out dummyTimeAx, out timeWindowedData); + data = timeWindowedData; + double[,] freqWindowedData; + double[] windowedFreqAxis; + ComplexTestCWT.FrequencyAxis.FrequencyWindowing(21, 23.5, ComplexTestCWT.FrequencyAxis.WaveletCenterFrequencies, data, CWTFrequencies.FrequencyUnits.WaveletFrequency, out windowedFreqAxis, out freqWindowedData); + data = freqWindowedData; + double[,] smoothedData = GaussianSmoothing.EllipticGaussianConvolution(data, 1, 700); + //double[,] smoothedData = data; + double[,] xyReflectedData = new double[smoothedData.GetLength(1), smoothedData.GetLength(0)]; + for (int i = 0; i < smoothedData.GetLength(0); i++) + { + for (int j = 0; j < smoothedData.GetLength(1); j++) + { + xyReflectedData[j, i] = smoothedData[i, j]; + } + } + var smoothedHeatMap = PlottingUtils.GenerateCWTHeatMap(xyReflectedData, "SmoothedBand2"); + string heatMapFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "smoothedHeatMap1freq700time.pdf"); + PlottingUtils.ExportPlotPDF(smoothedHeatMap, heatMapFilePath); + string beatInBeatFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "beatinbeatSlice1f700t.pdf"); + int[] targetBeatInBeatIndex = new int[] { -1 * Array.BinarySearch(windowedFreqAxis, 22.2) - 1 }; + var xyPlotBeatInBeat = PlottingUtils.GenerateXYPlotCWT(xyReflectedData, targetBeatInBeatIndex, PlottingUtils.PlotTitles.Custom, PlottingUtils.XYPlotOptions.Single, "Slice at 22.2"); + PlottingUtils.ExportPlotPDF(xyPlotBeatInBeat, beatInBeatFilePath); + string beatStdBeatFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "standardBeatSlice1f700t.pdf"); + int[] targetStandardBeatIndex = new int[] { -1 * Array.BinarySearch(windowedFreqAxis, 22.6) - 1 }; + var xyPlotStandardBeat = PlottingUtils.GenerateXYPlotCWT(xyReflectedData, targetStandardBeatIndex, PlottingUtils.PlotTitles.Custom, PlottingUtils.XYPlotOptions.Single, "Slice at 22.6"); + PlottingUtils.ExportPlotPDF(xyPlotStandardBeat, beatStdBeatFilePath); + + } + [Ignore("Long test time")] + [Test] + public void SlicePlottingTest() + { + double[,] data = ComplexTestCWT.OutputCWT.ModulusCalculation(); + double[,] timeWindowedData; + double[] dummyTimeAx; + CWTExtensions.TimeWindowing(0.04, 0.11, ComplexTestCWT.TimeAxis, data, out dummyTimeAx, out timeWindowedData); + data = timeWindowedData; + double[,] freqWindowedData; + double[] windowedFreqAxis; + ComplexTestCWT.FrequencyAxis.FrequencyWindowing(21, 23.5, ComplexTestCWT.FrequencyAxis.WaveletCenterFrequencies, data, CWTFrequencies.FrequencyUnits.WaveletFrequency, out windowedFreqAxis, out freqWindowedData); + data = freqWindowedData; + int targetBeatInBeatIndex = -1 * Array.BinarySearch(windowedFreqAxis, 22.2) - 1; + int targetStandardBeatIndex = -1 * Array.BinarySearch(windowedFreqAxis, 22.6) - 1; + double[] rawBeatingSlice = new double[data.GetLength(1)]; + double[] rawStandardSlice = new double[data.GetLength(1)]; + for (int i = 0; i < data.GetLength(1); i++) + { + rawStandardSlice[i] = data[targetStandardBeatIndex, i]; + rawBeatingSlice[i] = data[targetBeatInBeatIndex, i]; + } + double[] smoothedBeatingSlice = GaussianSmoothing.SliceEllipticGaussianConvolution(data, 1, 700, targetBeatInBeatIndex, 0); + double[] smoothedStandardSlice = GaussianSmoothing.SliceEllipticGaussianConvolution(data, 1, 700, targetStandardBeatIndex, 0); + var smoothedBeatingSlicePlot = PlottingUtils.Plot1DArray(smoothedBeatingSlice); + var smoothedStandardSlicePlot = PlottingUtils.Plot1DArray(smoothedStandardSlice); + string beatingSmoothSliceFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "beatingSliceSmoothed.pdf"); + string standardSmoothSliceFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "standardSliceSmoothed.pdf"); + PlottingUtils.ExportPlotPDF(smoothedStandardSlicePlot, standardSmoothSliceFilePath); + PlottingUtils.ExportPlotPDF(smoothedBeatingSlicePlot, beatingSmoothSliceFilePath); + + double[] beatingSliceDeriv = MinMaxIdentification.StandardArrayDerivative(smoothedBeatingSlice, 5); + var beatingSliceDerivPlot = PlottingUtils.Plot1DArray(beatingSliceDeriv); + string beatingSliceDerivFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "beatingSliceDeriv5dist.pdf"); + PlottingUtils.ExportPlotPDF(beatingSliceDerivPlot, beatingSliceDerivFilePath); + + double[] standardSliceDeriv = MinMaxIdentification.StandardArrayDerivative(smoothedStandardSlice, 5); + var standardSliceDerivPlot = PlottingUtils.Plot1DArray(standardSliceDeriv); + string standardSliceDerivFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "standardSliceDeriv5dist.pdf"); + PlottingUtils.ExportPlotPDF(standardSliceDerivPlot, standardSliceDerivFilePath); + double[] smoothStdSliceDeriv = GaussianSmoothing.GaussianSmoothing1D(standardSliceDeriv, 1000); + var stdSmoothSliceDerivPlot = PlottingUtils.Plot1DArray(smoothStdSliceDeriv); + string smoothStdSliceDerivPath = Path.Combine(ComplexTestCWT.WorkingPath, "smooth1000StDSliceDeriv5dist.pdf"); + PlottingUtils.ExportPlotPDF(stdSmoothSliceDerivPlot, smoothStdSliceDerivPath); + + } + [Ignore("Long test time")] + [Test] + public void TestDerivitiveAnalysis() + { + double[,] data = ComplexTestCWT.OutputCWT.ModulusCalculation(); + double[,] timeWindowedData; + double[] dummyTimeAx; + CWTExtensions.TimeWindowing(0.04, 0.11, ComplexTestCWT.TimeAxis, data, out dummyTimeAx, out timeWindowedData); + data = timeWindowedData; + double[,] freqWindowedData; + double[] windowedFreqAxis; + ComplexTestCWT.FrequencyAxis.FrequencyWindowing(21, 23.5, ComplexTestCWT.FrequencyAxis.WaveletCenterFrequencies, data, CWTFrequencies.FrequencyUnits.WaveletFrequency, out windowedFreqAxis, out freqWindowedData); + data = freqWindowedData; + int targetBeatInBeatIndex = -1 * Array.BinarySearch(windowedFreqAxis, 22.2) - 1; + int targetStandardBeatIndex = -1 * Array.BinarySearch(windowedFreqAxis, 22.6) - 1; + double[] rawBeatingSlice = new double[data.GetLength(1)]; + double[] rawStandardSlice = new double[data.GetLength(1)]; + for (int i = 0; i < data.GetLength(1); i++) + { + rawStandardSlice[i] = data[targetStandardBeatIndex, i]; + rawBeatingSlice[i] = data[targetBeatInBeatIndex, i]; + } + double[] smoothedStandardSlice = GaussianSmoothing.SliceEllipticGaussianConvolution(data, 1, 700, targetStandardBeatIndex, 0); + var intFilteredStandardSlice = MinMaxIdentification.IndexLinkedIntensityFiltering(smoothedStandardSlice, 0.2); + var intFiltStdPlot = PlottingUtils.PlotSortedPointDictionary(intFilteredStandardSlice); + string standardSliceFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "intFiltStdSlicePlot.pdf"); + PlottingUtils.ExportPlotPDF(intFiltStdPlot, standardSliceFilePath); + var downsampledStdDeriv = MinMaxIdentification.DownsampledSliceDerivative(intFilteredStandardSlice, 200); + var dsampledDerivPlot = PlottingUtils.PlotSortedPointDictionary(downsampledStdDeriv); + string standardDerivFilePath = Path.Combine(ComplexTestCWT.WorkingPath, "StdSliceDownsampledDeriv200Deriv.pdf"); + PlottingUtils.ExportPlotPDF(dsampledDerivPlot, standardDerivFilePath); + var smoothedStdDeriv = MinMaxIdentification.DownsampledDerivativeSmoothing(downsampledStdDeriv, 1000, 200); + var smoothedDerivPlot = PlottingUtils.PlotSortedPointDictionary(smoothedStdDeriv); + var smoothedStdDerivFile = Path.Combine(ComplexTestCWT.WorkingPath, "smoothed1000dsampled200derivStd.pdf"); + PlottingUtils.ExportPlotPDF(smoothedDerivPlot, smoothedStdDerivFile); + var secondStdDeriv = MinMaxIdentification.StandardDerivative(smoothedStdDeriv, 200); + var secondDerivPlot = PlottingUtils.PlotSortedPointDictionary(secondStdDeriv); + var secondDerivFile = Path.Combine(ComplexTestCWT.WorkingPath, "smoothedDerivandsecondderiv.pdf"); + PlottingUtils.ExportPlotPDF(secondDerivPlot, secondDerivFile); + } + } +} diff --git a/TestFCWTAPI/MinMaxIdentificationSmallDataTests.cs b/TestFCWTAPI/MinMaxIdentificationSmallDataTests.cs new file mode 100644 index 0000000..a068ad7 --- /dev/null +++ b/TestFCWTAPI/MinMaxIdentificationSmallDataTests.cs @@ -0,0 +1,161 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using FCWTNET; +using System.IO; +using NUnit.Framework; + +namespace TestFCWTAPI +{ + public class MinMaxIdentificationSmallDataTests + { + [Test] + public static void StandardArrayDerivative() + { + double[] testArray = new double[15]; + for (int i = 1; i <= testArray.Length; i++) + { + testArray[i - 1] = i * i; + } + double[] testDerivative = MinMaxIdentification.StandardArrayDerivative(testArray, 2); + Assert.AreEqual(12, testDerivative[5]); + Assert.AreEqual(5, testDerivative[1]); + Assert.AreEqual(28, testDerivative[14]); + Assert.Throws(() => MinMaxIdentification.StandardArrayDerivative(testArray, 9)); + } + [Test] + public static void TestIndexLinkedIntensityFiltering() + { + double[] testArray = new double[16]; + for (int i = 0; i < testArray.Length; i++) + { + testArray[i] = Math.Sin(Math.PI * i / 8); + } + SortedDictionary testIntensityFiltering = MinMaxIdentification.IndexLinkedIntensityFiltering(testArray, 0.5); + foreach(var pair in testIntensityFiltering) + { + Assert.IsTrue(Math.Abs(pair.Value) >= 0.5); + Assert.AreEqual(testArray[pair.Key], pair.Value); + } + Assert.Throws(() => MinMaxIdentification.IndexLinkedIntensityFiltering(testArray, -1)); + Assert.Throws(() => MinMaxIdentification.IndexLinkedIntensityFiltering(testArray, 2)); + } + [Test] + public static void TestDownsampledSliceDerivative() + { + double[] testArray = new double[32]; + for (int i = 0; i < testArray.Length; i++) + { + testArray[i] = Math.Sin(Math.PI * i / 16); + } + SortedDictionary filteredPointSet = MinMaxIdentification.IndexLinkedIntensityFiltering(testArray, 0.2); + int testDerivativeDistance = 1; + SortedDictionary downsampledDeriv = MinMaxIdentification.DownsampledSliceDerivative(filteredPointSet, testDerivativeDistance); + Assert.AreEqual(filteredPointSet.Count / (2 * testDerivativeDistance + 1), downsampledDeriv.Count); + foreach(var point in downsampledDeriv) + { + double expectedDeriv = (double)(filteredPointSet[point.Key + 1] - filteredPointSet[point.Key - 1]) / 2; + Assert.AreEqual(expectedDeriv, point.Value, 0.001); + } + Assert.Throws(() => MinMaxIdentification.DownsampledSliceDerivative(filteredPointSet, 33)); + } + [Test] + public static void TestDownsampledDerivativeSmoothing() + { + double[] testArray = new double[128]; + for (int i = 0; i < testArray.Length; i++) + { + if (i < 32 || (i > 64 && i < 96)) + { + testArray[i] = Math.Sin(Math.PI * i / 32); + } + else + { + testArray[i] = 0.3 * Math.Sin(Math.PI * i / 32); + } + + } + SortedDictionary filteredPointSet = MinMaxIdentification.IndexLinkedIntensityFiltering(testArray, 0.2); + int testDerivativeDistance = 1; + SortedDictionary downsampledDeriv = MinMaxIdentification.DownsampledSliceDerivative(filteredPointSet, testDerivativeDistance); + List SmoothableCluster = new List(); + List UnsmoothableCluster = new List(); + + foreach(var point in downsampledDeriv) + { + if (point.Key < 32) + { + SmoothableCluster.Add(point.Value); + } + else if(point.Key < 64) + { + UnsmoothableCluster.Add(point.Value); + } + } + double[] smoothedCluster = GaussianSmoothing.GaussianSmoothing1D(SmoothableCluster.ToArray(), 1); + SortedDictionary testSmoothedDeriv = MinMaxIdentification.DownsampledDerivativeSmoothing(downsampledDeriv, 3, testDerivativeDistance); + int firstSmoothableClusterCounter = 0; + int secondSmoothableClusterCounter = 0; + int firstUnsmoothedClusterCounter = 0; + int secondUnsmoothedClusterCounter = 0; + foreach (var point in testSmoothedDeriv) + { + if (point.Key < 32) + { + Assert.AreEqual(smoothedCluster[firstSmoothableClusterCounter], point.Value, 0.00001); + firstSmoothableClusterCounter++; + } + else if (point.Key < 64) + { + Assert.AreEqual(UnsmoothableCluster[firstUnsmoothedClusterCounter], point.Value, 0.00001); + firstUnsmoothedClusterCounter++; + } + else if (point.Key < 96) + { + Assert.AreEqual(smoothedCluster[secondSmoothableClusterCounter], point.Value, 0.00001); + secondSmoothableClusterCounter++; + } + else + { + Assert.AreEqual(UnsmoothableCluster[secondUnsmoothedClusterCounter], point.Value, 0.00001); + secondUnsmoothedClusterCounter++; + } + } + + } + [Test] + public static void TestStandardDerivative() + { + SortedDictionary sampleNoDerivativeDistance = new SortedDictionary(); + SortedDictionary sampleWithDerivativeDistance = new SortedDictionary(); + int derivativeDistance = 1; + for(int i = 0; i < 10; i++) + { + sampleNoDerivativeDistance.Add(i, i * i); + sampleWithDerivativeDistance.Add((i * (2 * derivativeDistance + 1)) + derivativeDistance, i * i); + } + for (int i = 13; i < 20; i++) + { + sampleNoDerivativeDistance.Add(i, i * i); + sampleWithDerivativeDistance.Add((i * (2 * derivativeDistance + 1)) + derivativeDistance, i * i); + } + SortedDictionary derivativeNoDerivativeDistance = MinMaxIdentification.StandardDerivative(sampleNoDerivativeDistance); + SortedDictionary derivativeWithDerivativeDistance = MinMaxIdentification.StandardDerivative(sampleWithDerivativeDistance, 1); + Assert.AreEqual(sampleWithDerivativeDistance.Count - 4, derivativeWithDerivativeDistance.Count); + Assert.AreEqual(sampleNoDerivativeDistance.Count - 4, derivativeNoDerivativeDistance.Count); + foreach(var point in derivativeNoDerivativeDistance) + { + double expectedDerivative = (sampleNoDerivativeDistance[point.Key + 1] - sampleNoDerivativeDistance[point.Key - 1]) / 2; + Assert.AreEqual(expectedDerivative, point.Value); + } + foreach(var point in derivativeWithDerivativeDistance) + { + double expectedDerivative = (sampleWithDerivativeDistance[point.Key + 2 * derivativeDistance + 1] - sampleWithDerivativeDistance[point.Key - (2 * derivativeDistance + 1)]) / (2 * (2 * derivativeDistance + 1)); + Assert.AreEqual(expectedDerivative, point.Value); + } + } + + } +} diff --git a/TestFCWTAPI/MinMaxIdentificationTest.cs b/TestFCWTAPI/MinMaxIdentificationTest.cs deleted file mode 100644 index 1682817..0000000 --- a/TestFCWTAPI/MinMaxIdentificationTest.cs +++ /dev/null @@ -1,54 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using System.Threading.Tasks; -using FCWTNET; -using System.IO; -using NUnit.Framework; - -namespace TestFCWTAPI -{ - public class MinMaxIdentificationTest - { - public CWTObject MyoglobinTestCWT { get; set; } - public CWTObject ComplexTestCWT { get; set; } - [OneTimeSetUp] - public void TestDataSetup() - { - // Create test data from real transients - string myoZ13TransientPath = Path.Combine(TestContext.CurrentContext.TestDirectory, "TestDataFiles", "SampleTransient_Myo_z13.csv"); - TransientData myoZ13Transient = new(myoZ13TransientPath); - myoZ13Transient.ImportTransientData(800000); - string testPath = Path.Combine(TestContext.CurrentContext.TestDirectory, "TestDataFiles"); - var myoZ13CWT = new CWTObject(myoZ13Transient.Data, 2, 3, 125, (float)(2 * Math.Tau), 4, false, myoZ13Transient.SamplingRate, myoZ13Transient.CalibrationFactor, testPath); - myoZ13CWT.PerformCWT(); - myoZ13CWT.CalculateTimeAxis(); - myoZ13CWT.CalculateFrequencyAxis(); - MyoglobinTestCWT = myoZ13CWT; - // Plot modulus and real heatmap for basic visualization - string complexTransientPath = Path.Combine(TestContext.CurrentContext.TestDirectory, "TestDataFiles", "ComplexTransient_protein_mixture_512ms.csv"); - TransientData complexTransient = new(complexTransientPath); - complexTransient.ImportTransientData(700000); - var complexTransientCWT = new CWTObject(complexTransient.Data, 2, 3, 125, (float)(15 * Math.Tau), 4, false, complexTransient.SamplingRate, complexTransient.CalibrationFactor, testPath); - complexTransientCWT.PerformCWT(); - complexTransientCWT.CalculateTimeAxis(); - complexTransientCWT.CalculateFrequencyAxis(); - ComplexTestCWT = complexTransientCWT; - } - [Test] - // Test method for generating plots to carefully examine the CWT of the sample MS transient data - public void PlotGeneration() - { - MyoglobinTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "myoZ13_modulus.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "myo z13"); - MyoglobinTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "myoZ13_modulus_Zoomed.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "myo z13", 0.075, 0.095, 2.1, 3); - MyoglobinTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "myoZ13_modulus_HiZoom.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "myo z13", 0.082, 0.083, 2.4, 2.5); - MyoglobinTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Real, "myoZ13_Real.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "myo z13"); - ComplexTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "complexTransient_modulus.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "complexTransient"); - ComplexTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "zoomed_Band2complexTransient_modulus.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "complexTransient", 0.06, 0.08, 19, 30); - ComplexTestCWT.GenerateHeatMap(CWTObject.CWTFeatures.Modulus, "zoomed_Band1complexTransient_modulus.pdf", CWTFrequencies.FrequencyUnits.WaveletFrequency, "complexTransient", 0.02, 0.04, 19, 30); - - - } - } -} diff --git a/TestFCWTAPI/PlotingTests.cs b/TestFCWTAPI/PlotingTests.cs index 06e2bbc..4b83e2c 100644 --- a/TestFCWTAPI/PlotingTests.cs +++ b/TestFCWTAPI/PlotingTests.cs @@ -65,6 +65,8 @@ public void TestGenerateCWTHeatMap() } var testHeatMap = PlottingUtils.GenerateCWTHeatMap(testData, "Test Heatmap", timeAxis, freqAxis); PlottingUtils.ExportPlotPDF(testHeatMap, "testCWTheatmap.pdf"); + var noAxesHeatMap = PlottingUtils.GenerateCWTHeatMap(testData, "Test Heatmap"); + PlottingUtils.ExportPlotPDF(noAxesHeatMap, "testNoAxesCWTheatmap.pdf"); } [Test] public void TestGenerateCWTContourPlot() @@ -79,6 +81,8 @@ public void TestGenerateCWTContourPlot() } var testContour = PlottingUtils.GenerateCWTContourPlot(testData, "Test Contour", timeAxis, freqAxis); PlottingUtils.ExportPlotPDF(testContour, "testCWTcontour.pdf"); + var noAxesContourMap = PlottingUtils.GenerateCWTContourPlot(testData, "Test Contour"); + PlottingUtils.ExportPlotPDF(testContour, "testNoAxesCWTcontour.pdf"); } [Test] public void TestGenerateXYPlotCWT() @@ -97,6 +101,7 @@ public void TestGenerateXYPlotCWT() var testEvolutionPlot = PlottingUtils.GenerateXYPlotCWT(testData, testMultiple, timeAxis, freqAxis, PlottingUtils.PlotTitles.Evolution, PlottingUtils.XYPlotOptions.Evolution); var testSinglePlot = PlottingUtils.GenerateXYPlotCWT(testData, testSingle, timeAxis, freqAxis, PlottingUtils.PlotTitles.Single, PlottingUtils.XYPlotOptions.Evolution); var testCustomTitle = PlottingUtils.GenerateXYPlotCWT(testData, testSingle, timeAxis, freqAxis, PlottingUtils.PlotTitles.Custom, PlottingUtils.XYPlotOptions.Single, "Custom Title"); + var testNoAxesEvolution = PlottingUtils.GenerateXYPlotCWT(testData, testMultiple, PlottingUtils.PlotTitles.Custom, PlottingUtils.XYPlotOptions.Evolution, "NoAxesEvolution"); Assert.Throws(() => PlottingUtils.GenerateXYPlotCWT(testData, testSingle, timeAxis, freqAxis, PlottingUtils.PlotTitles.Custom, PlottingUtils.XYPlotOptions.Single)); double[] badTimeAxis = new double[980]; double[] badFreqAxis = new double[980]; @@ -106,7 +111,7 @@ public void TestGenerateXYPlotCWT() PlottingUtils.ExportPlotPDF(testEvolutionPlot, "testevolutionCWT.pdf"); PlottingUtils.ExportPlotPDF(testSinglePlot, "testsingleCWT.pdf"); PlottingUtils.ExportPlotPDF(testCustomTitle, "customtitleCWT.pdf"); - + PlottingUtils.ExportPlotPDF(testNoAxesEvolution, "No_Axes_XY_Plot.pdf"); } }