diff --git a/src/DIPOL-UF/DipolUfApp.cs b/src/DIPOL-UF/DipolUfApp.cs index 1b5ade6b..a3bd110e 100644 --- a/src/DIPOL-UF/DipolUfApp.cs +++ b/src/DIPOL-UF/DipolUfApp.cs @@ -1,7 +1,7 @@  #define HOST_SERVER -// #define IN_PROCESS +#define IN_PROCESS using System; using System.Diagnostics; @@ -31,6 +31,7 @@ private static int Main() var host = connStr is null ? null : new DipolHost(new Uri(connStr)); host?.Open(); #else + var connStr = UiSettingsProvider.Settings.GetArray("RemoteLocations")?.FirstOrDefault(); var pInfo = new ProcessStartInfo() diff --git a/src/DIPOL-UF/Helper.cs b/src/DIPOL-UF/Helper.cs index e652a951..7c7c3efe 100644 --- a/src/DIPOL-UF/Helper.cs +++ b/src/DIPOL-UF/Helper.cs @@ -26,7 +26,7 @@ namespace DIPOL_UF { public static class Helper { - + private static readonly double DEpsilon = Math.Pow(2, -53); public static Dispatcher UiDispatcher => Application.Current?.Dispatcher ?? throw new InvalidOperationException("Dispatcher is unavailable."); @@ -410,7 +410,61 @@ DescriptionAttribute attr : key); - + + internal static double Interpolate(double x1, double x2, double y1, double y2, double x0) + { + if (Math.Abs(x2 - x1) < 1e-15) + { + return 0.0; + } + return y1 + (y2 - y1) / (x2 - x1) * (x0 - x1); + } + + internal static double Percentile(ReadOnlySpan data, double p) + { + if (data.IsEmpty || p is < 0 or > 1) + { + return double.NaN; + } + + if (Math.Abs(p) < DEpsilon) + { + return data[0]; + } + + if (Math.Abs(1 - p) < DEpsilon) + { + return data[data.Length - 1]; + } + + + var buffer = ArrayPool.Shared.Rent(data.Length); + try + { + data.CopyTo(buffer); + Array.Sort(buffer, 0, data.Length); + + // Strictly non-negative value between `0` and `n - 1` + var x = (data.Length - 1) * p; + + // Integral part of `x` + var i = (int)Math.Floor(x); + // Non-integral part of `x` + var g = x - i; + + if (i == data.Length - 1) + { + return buffer[i]; + } + + return (1 - g) * buffer[i] + g * buffer[i + 1]; + } + finally + { + ArrayPool.Shared.Return(buffer); + } + } + } internal class CameraTupleOrderComparer : IComparer<(string Id, IDevice Camera)> diff --git a/src/DIPOL-UF/Models/DipolImagePresenter.cs b/src/DIPOL-UF/Models/DipolImagePresenter.cs index 9709eb39..aa1c5bb1 100644 --- a/src/DIPOL-UF/Models/DipolImagePresenter.cs +++ b/src/DIPOL-UF/Models/DipolImagePresenter.cs @@ -40,31 +40,22 @@ private enum GeometryLayer public struct GaussianFitResults : IEquatable { - private static double Const { get; } = 2 * Math.Sqrt(2 * Math.Log(2)); - - public double Baseline { get; } - public double Scale { get; } public double Center { get; } - public double Sigma { get; } - + public double FWHM { get; } public int Origin { get; set; } - public double FWHM => Const * Sigma; - - public bool IsValid => Baseline >= 0 && Scale > 0 && Sigma > 0; - public GaussianFitResults(double baseLine, double scale, double center, double sigma) => - (Baseline, Scale, Center, Sigma, Origin) = (baseLine, scale, center, sigma, 0); + public GaussianFitResults(double center, double fwhm) => + (Center, FWHM, Origin) = (center, fwhm, 0); public override bool Equals(object? obj) => obj is GaussianFitResults other && Equals(other); public bool Equals(GaussianFitResults other) => - (Baseline, Scale, Center, Sigma, Origin) - == (other.Baseline, other.Scale, other.Center, other.Sigma, other.Origin); + (Center, FWHM, Origin) == (other.Center, other.FWHM, other.Origin); public override int GetHashCode() => - HashCode.Combine(Baseline, Scale, Center, Sigma); + HashCode.Combine(Center, FWHM, Origin); public static bool operator ==(GaussianFitResults left, GaussianFitResults right) => left.Equals(right); @@ -670,7 +661,7 @@ private Point GetImageScale(Point p) GetImageScale(p.Y)) : p; - private (GaussianFitResults Row, GaussianFitResults Column) GetImageStatistics(Point clickPosition) + private (GaussianFitResults Row, GaussianFitResults Column) GetImageStatistics() { if (LastKnownImageControlSize.IsEmpty || DisplayedImage == null) { @@ -681,7 +672,7 @@ private Point GetImageScale(Point p) var sizePix = GetPixelScale(SamplerGeometry!.Size); var halfSizePix = GetPixelScale(SamplerGeometry!.HalfSize); var image = _sourceImage!; - + var rowStart = Math.Max(0, (int) (centerPix.Y - halfSizePix.Height)); var colStart = Math.Max(0, (int) (centerPix.X - halfSizePix.Width)); var width = Math.Min((int) sizePix.Width, image.Width - colStart); @@ -692,9 +683,11 @@ private Point GetImageScale(Point p) { buffer = ArrayPool.Shared.Rent(width * height); var tempView = new Span2D(buffer, height, width); - CastViewToDouble(image, colStart, rowStart, width, height, tempView); - var stats = ComputeFullWidthHalfMax(tempView); + + var backGround = Helper.Percentile(buffer.AsSpan(0, width * height), 0.025); + + var stats = ComputeFullWidthHalfMax(tempView, backGround); stats.Column.Origin = colStart; stats.Row.Origin = rowStart; return stats; @@ -705,12 +698,17 @@ private Point GetImageScale(Point p) { ArrayPool.Shared.Return(buffer); } + + } } - private List<(int X, int Y)> GetPixelsInArea(GeometryLayer layer) + private IReadOnlyList<(int X, int Y)> GetPixelsInArea(GeometryLayer layer) { - if (!LastKnownImageControlSize.IsEmpty && DisplayedImage != null) + if ( + !LastKnownImageControlSize.IsEmpty && + DisplayedImage is not null + ) { var centerPix = GetPixelScale(SamplerCenterPos); @@ -743,6 +741,10 @@ private Point GetImageScale(Point p) switch (layer) { case GeometryLayer.Annulus: + if (SamplerGeometry is null || GapGeometry is null) + { + return new List<(int X, int Y)> { (0, 0) }; + } var annulusPixels = new List<(int X, int Y)>( (pixelXLims.Max - pixelXLims.Min) * @@ -751,16 +753,26 @@ private Point GetImageScale(Point p) for (var xPix = pixelXLims.Min; xPix <= pixelXLims.Max; xPix++) for (var yPix = pixelYLims.Min; yPix <= pixelYLims.Max; yPix++) { - if ((SamplerGeometry?.IsInsideChecker(xPix, yPix, centerPix, halfSizePix, - thcknssPix) ?? false) && - !(GapGeometry?.IsInsideChecker(xPix, yPix, centerPix, gapHalfSizePix, - thcknssPix) ?? true)) + if (SamplerGeometry.IsInsideChecker( + xPix, yPix, centerPix, halfSizePix, + thcknssPix + ) && + !GapGeometry.IsInsideChecker( + xPix, yPix, centerPix, gapHalfSizePix, + thcknssPix + )) + { annulusPixels.Add((X: xPix, Y: yPix)); + } } return annulusPixels; case GeometryLayer.Gap: + if (ApertureGeometry is null || GapGeometry is null) + { + return new List<(int X, int Y)> { (0, 0) }; + } var gapPixels = new List<(int X, int Y)>( (gapPixelXLims.Max - gapPixelXLims.Min) * @@ -769,16 +781,26 @@ private Point GetImageScale(Point p) for (var xPix = gapPixelXLims.Min; xPix <= gapPixelXLims.Max; xPix++) for (var yPix = gapPixelYLims.Min; yPix <= gapPixelYLims.Max; yPix++) { - if ((GapGeometry?.IsInsideChecker(xPix, yPix, centerPix, halfSizePix, - thcknssPix) ?? false) && - !(ApertureGeometry?.IsInsideChecker(xPix, yPix, centerPix, gapHalfSizePix, - thcknssPix) ?? true)) + if (GapGeometry.IsInsideChecker( + xPix, yPix, centerPix, halfSizePix, + thcknssPix + ) && + !ApertureGeometry.IsInsideChecker( + xPix, yPix, centerPix, gapHalfSizePix, + thcknssPix + )) + { gapPixels.Add((X: xPix, Y: yPix)); + } } return gapPixels; case GeometryLayer.Aperture: + if (ApertureGeometry is null) + { + return new List<(int X, int Y)> { (0, 0) }; + } var aperturePixels = new List<(int X, int Y)>( (aperturePixelXLims.Max - aperturePixelXLims.Min) * @@ -787,9 +809,13 @@ private Point GetImageScale(Point p) for (var xPix = aperturePixelXLims.Min; xPix <= aperturePixelXLims.Max; xPix++) for (var yPix = aperturePixelYLims.Min; yPix <= aperturePixelYLims.Max; yPix++) { - if (ApertureGeometry?.IsInsideChecker(xPix, yPix, centerPix, - apertureHalfSizePix, thcknssPix) ?? false) + if (ApertureGeometry.IsInsideChecker( + xPix, yPix, centerPix, + apertureHalfSizePix, thcknssPix + )) + { aperturePixels.Add((X: xPix, Y: yPix)); + } } return aperturePixels; @@ -825,7 +851,7 @@ private void ImageDoubleClickCommandExecute(MouseEventArgs args) private (GaussianFitResults Row, GaussianFitResults Column) ImageRightClickCommandExecute(MouseEventArgs args) { - if (!(args.Source is FrameworkElement elem)) + if (args.Source is not FrameworkElement elem) { return default; } @@ -839,7 +865,7 @@ private void ImageDoubleClickCommandExecute(MouseEventArgs args) #endif try { - return GetImageStatistics(pos); + return GetImageStatistics(); } catch (Exception e) { @@ -870,7 +896,7 @@ GeometryDescriptor CommonRectangle(double size, double thickness) var path = new List>>(4) { Tuple.Create>( - new Point(0, 0), null), + new Point(0, 0), (_, _) => { }), Tuple.Create>( new Point(size, 0), (cont, pt) => cont.LineTo(pt, true, false)), Tuple.Create>( @@ -1040,68 +1066,69 @@ private static void CastViewToDouble(Image image, int col, int row, int width, i } - private static (GaussianFitResults Row, GaussianFitResults Column) ComputeFullWidthHalfMax(ReadOnlySpan2D data) + private static (GaussianFitResults Row, GaussianFitResults Column) ComputeFullWidthHalfMax(ReadOnlySpan2D data, double bkg) { var center = FindBrightestPixel(data); ReadOnlySpan row = IK.ILSpanCasts.SpanExtensions.GetRow(data, center.Row); - var rowBuff = new double[data.Width]; - row.CopyTo(rowBuff); - var colBuff = new double[data.Height]; - ReadOnlySpan2D column = data.Slice(0, center.Column, 1, data.Height); - for (var i = 0; i < colBuff.Length; i++) - { - colBuff[i] = column[i, 0]; - } - var args = new double[data.Width]; - for (var i = 0; i < args.Length; i++) + Span rowView = data.Width > 128 + ? new double[data.Width] + : stackalloc double[data.Width]; + + Span colView = data.Height > 128 + ? new double[data.Height] + : stackalloc double[data.Height]; + + ReadOnlySpan2D column = data.Slice(0, center.Column, width: 1, height: data.Height); + + row.CopyTo(rowView); + + for (var i = 0; i < colView.Length; i++) { - args[i] = i; + colView[i] = column[i, 0]; } - var rowStats = ComputeFullWidthHalfMax(args, rowBuff); - var colStats = ComputeFullWidthHalfMax(args, colBuff); + var rowStats = + new GaussianFitResults(center.Row, ComputeFullWidthHalfMax(rowView, center.Column, bkg)); + var colStats = + new GaussianFitResults(center.Column, ComputeFullWidthHalfMax(colView, center.Row, bkg)); return (rowStats, colStats); } - private static GaussianFitResults ComputeFullWidthHalfMax(double[] arg, double[] data) + private static double ComputeFullWidthHalfMax(ReadOnlySpan data, int pos, double bkg) { - var (baseLine, max) = Helper.MinMax(data); - - // Gaussian with non-zero baseline - static double ExpFunc2(double scale, double sigma, double center, double background, double x) => - background + scale * Math.Exp(-(x - center) * (x - center) / 2 / sigma / sigma); - - // Distance between prediction and actual data - double DistFunc(double scale, double sigma, double center, double background) => - Distance.Euclidean( - // Prediction - Generate.Map(arg, t => ExpFunc2(scale, sigma, center, background, t)), - // Data - data - ); - - // Minimizes Euclidean distance - MinimizationResult minimizationResult = NelderMeadSimplex.Minimum( - // Vector -> double distance function - ObjectiveFunction.Value(v => DistFunc(v[0], v[1], v[2], v[3])), - // Initial guess - CreateVector.Dense( - new[] - { - max, - data.Length / 4.0, - data.Length / 2.0, - baseLine - } - ) - ); + if (pos > data.Length || pos < 0) + { + return default; + } + + var max = data[pos] - bkg; + + (double left, double right) = (0, data.Length - 1); + + for(var i = pos; i >= 0; i--) + { + if (data[i] - bkg < 0.5 * max) + { + left = Helper.Interpolate(data[i] - bkg, data[i + 1] - bkg, i, i + 1, 0.5 * max); + break; + } + } + + for (var i = pos; i < data.Length; i++) + { + if (data[i] - bkg < 0.5 * max) + { + right = Helper.Interpolate(data[i - 1] - bkg, data[i] - bkg, i - 1, i, 0.5 * max); + break; + } + } - Vector prediction = minimizationResult.MinimizingPoint; - return new GaussianFitResults(prediction[3], prediction[0], prediction[2], prediction[1]); + return right - left; } + private static (int Row, int Column) FindBrightestPixel(ReadOnlySpan2D data) { var result = (Row: data.Height / 2, Column: data.Width / 2); diff --git a/src/DIPOL-UF/Properties/AssemblyInfo.cs b/src/DIPOL-UF/Properties/AssemblyInfo.cs index 2bdd064e..958372a3 100644 --- a/src/DIPOL-UF/Properties/AssemblyInfo.cs +++ b/src/DIPOL-UF/Properties/AssemblyInfo.cs @@ -51,5 +51,5 @@ // You can specify all the values or you can default the Build and Revision Numbers // by using the '*' as shown below: // [assembly: AssemblyVersion("1.0.*")] -[assembly: AssemblyVersion("2.6.0")] +[assembly: AssemblyVersion("2.6.1")] //[assembly: AssemblyFileVersion("1.0.*")] diff --git a/src/DIPOL-UF/ViewModels/DipolImagePresenterViewModel.cs b/src/DIPOL-UF/ViewModels/DipolImagePresenterViewModel.cs index 4ca5e6ad..c20708ca 100644 --- a/src/DIPOL-UF/ViewModels/DipolImagePresenterViewModel.cs +++ b/src/DIPOL-UF/ViewModels/DipolImagePresenterViewModel.cs @@ -300,11 +300,6 @@ private void DisplayFWHMEstimates( DipolImagePresenter.GaussianFitResults col ) { - // Handles incorrect fits - if (!row.IsValid || !col.IsValid) - { - (row, col) = (default, default); - } var builder = new StringBuilder(); //builder.AppendLine($"{Properties.Localization.FWHM_Center}\t({col.Origin + col.Center:F2}, {row.Origin + row.Center:F2})"); builder.AppendFormat(Properties.Localization.FWHM_Center, col.Origin + col.Center, row.Origin + row.Center);