using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Threading.Tasks; using CaeMesh; using CaeGlobals; using DynamicTypeDescriptor; using System.Collections.Concurrent; using Priority_Queue; namespace CaeResults { public class CloudInterpolator { // Variables private int _nx; private int _ny; private int _nz; private int _nxy; private int _numOfValues; private double _deltaX; private double _deltaY; private double _deltaZ; private CloudPoint[] _cloudPoints; private BoundingBox _sourceBox; private BoundingBox[] _regionBoxes; // Tag of each bounding box contains dictionary // Variables public int NumOfValues { get { return _numOfValues; } } // number of values/dimensions per point public CloudPoint[] CloudPoints { get { return _cloudPoints; } } // Constructor public CloudInterpolator(CloudPoint[] cloudPoints) { _cloudPoints = cloudPoints; // _sourceBox = ComputeAllPointsBoundingBox(_cloudPoints); _sourceBox.InflateIfThinn(1E-6); // double l = _sourceBox.GetDiagonal() / 100; _nx = (int)Math.Ceiling(_sourceBox.GetXSize() / l); _ny = (int)Math.Ceiling(_sourceBox.GetYSize() / l); _nz = (int)Math.Ceiling(_sourceBox.GetZSize() / l); // int currNumBoxes = _nx * _ny * _nz; int maxNumBoxes = 10_000_000; if (currNumBoxes > maxNumBoxes) { double factor = Math.Pow((double)maxNumBoxes / currNumBoxes, 0.333333); l /= factor; } // _nx = (int)Math.Ceiling(_sourceBox.GetXSize() / l); _ny = (int)Math.Ceiling(_sourceBox.GetYSize() / l); _nz = (int)Math.Ceiling(_sourceBox.GetZSize() / l); // _nxy = _nx * _ny; _deltaX = _sourceBox.GetXSize() / _nx; _deltaY = _sourceBox.GetYSize() / _ny; _deltaZ = _sourceBox.GetZSize() / _nz; // _regionBoxes = AssignPointsToRegions(_cloudPoints, _sourceBox, _nx, _ny, _nz); } public void InterpolateAt(double[] point, CloudInterpolatorEnum interpolator, double radius, out double[] distance, out double[] values) { if (interpolator == CloudInterpolatorEnum.ClosestPoint) { InterpolateByClosestPoint(point, out distance, out values); } else if (interpolator == CloudInterpolatorEnum.Gauss) { InterpolateByGauss(point, radius, out distance, out values); } else if (interpolator == CloudInterpolatorEnum.Shepard) { InterpolateByShepard(point, radius, out distance, out values); } else throw new NotSupportedException(); // values = values.ToArray(); // copy } public void InterpolateByClosestPoint(double[] point, out double[] distance, out double[] values) { Dictionary regions = GetClosestRegions(point, -1); // double absX; double absY; double absZ; double d; double minD = double.MaxValue; CloudPoint bestPoint = new CloudPoint(); foreach (var regionEntry in regions) { if (regionEntry.Value.IsMaxOutsideDistance2SmallerThan(point, minD * minD)) { foreach (var cloudPoint in (HashSet)regionEntry.Value.Tag) { absX = Math.Abs(cloudPoint.Coor[0] - point[0]); if (absX < minD) { absY = Math.Abs(cloudPoint.Coor[1] - point[1]); if (absY < minD) { absZ = Math.Abs(cloudPoint.Coor[2] - point[2]); if (absZ < minD) { d = Math.Sqrt(absX * absX + absY * absY + absZ * absZ); // if (d < minD) { minD = d; bestPoint = cloudPoint; } } } } } } } // distance = new double[] { bestPoint.Coor[0] - point[0], bestPoint.Coor[1] - point[1], bestPoint.Coor[2] - point[2]}; // values = bestPoint.Values; // copy values since multiply with magnitude } public void InterpolateByGauss(double[] point, double radius, out double[] distance, out double[] values) { Dictionary regions = GetClosestRegions(point, radius); // double absX; double absY; double absZ; double d; double s = 1; double w; double radius2 = radius * radius; Dictionary cloudPointWeight = new Dictionary(); foreach (var regionEntry in regions) { if (regionEntry.Value.IsMaxOutsideDistance2SmallerThan(point, radius2)) { foreach (var cloudPoint in (HashSet)regionEntry.Value.Tag) { absX = Math.Abs(cloudPoint.Coor[0] - point[0]); if (absX < radius) { absY = Math.Abs(cloudPoint.Coor[1] - point[1]); if (absY < radius) { absZ = Math.Abs(cloudPoint.Coor[2] - point[2]); if (absZ < radius) { d = Math.Sqrt(absX * absX + absY * absY + absZ * absZ); // w = Math.Exp(Math.Pow(-s * d / radius, 2)); cloudPointWeight.Add(cloudPoint, w); } } } } } } // if (cloudPointWeight.Count > 0) { double wSum = 0; foreach (var entry in cloudPointWeight) wSum += entry.Value; // values = new double[_numOfValues]; // CloudPoint cp; foreach (var entry in cloudPointWeight) { cp = entry.Key; for (int i = 0; i < _numOfValues; i++) { values[i] += (entry.Value / wSum) * cp.Values[i]; } } // distance = new double[] { 0, 0, 0 }; } else { values = new double[_numOfValues]; distance = new double[] { radius, radius, radius }; } // } public void InterpolateByShepard(double[] point, double radius, out double[] distance, out double[] values) { Dictionary regions = GetClosestRegions(point, radius); // bool foundExactPosition = false; double absX; double absY; double absZ; double d; double w; double radius2 = radius * radius; Dictionary cloudPointWeight = new Dictionary(); foreach (var regionEntry in regions) { if (regionEntry.Value.IsMaxOutsideDistance2SmallerThan(point, radius2)) { foreach (var cloudPoint in (HashSet)regionEntry.Value.Tag) { absX = Math.Abs(cloudPoint.Coor[0] - point[0]); if (absX < radius) { absY = Math.Abs(cloudPoint.Coor[1] - point[1]); if (absY < radius) { absZ = Math.Abs(cloudPoint.Coor[2] - point[2]); if (absZ < radius) { d = Math.Sqrt(absX * absX + absY * absY + absZ * absZ); // if (d == 0) { cloudPointWeight.Clear(); cloudPointWeight.Add(cloudPoint, 1); foundExactPosition = true; break; } else { w = 1 / Math.Pow(d, 2); cloudPointWeight.Add(cloudPoint, w); } } } } } } if (foundExactPosition) break; } // if (cloudPointWeight.Count > 0) { double wSum = 0; foreach (var entry in cloudPointWeight) wSum += entry.Value; // values = new double[_numOfValues]; // CloudPoint cp; foreach (var entry in cloudPointWeight) { cp = entry.Key; for (int i = 0; i < _numOfValues; i++) { values[i] += (entry.Value / wSum) * cp.Values[i]; } } // distance = new double[] { 0, 0, 0 }; } else { values = new double[_numOfValues]; distance = new double[] { radius, radius, radius }; } // } private Dictionary GetClosestRegions(double[] point, double radius) { int i; int j; int k; int mini; int maxi; int minj; int maxj; int mink; int maxk; int index; BoundingBox bb; Dictionary regions = new Dictionary(); int num; int layer; // i = (int)Math.Floor((point[0] - _sourceBox.MinX) / _deltaX); j = (int)Math.Floor((point[1] - _sourceBox.MinY) / _deltaY); k = (int)Math.Floor((point[2] - _sourceBox.MinZ) / _deltaZ); if (i < 0) i = 0; else if (i >= _nx) i = _nx - 1; if (j < 0) j = 0; else if (j >= _ny) j = _ny - 1; if (k < 0) k = 0; else if (k >= _nz) k = _nz - 1; index = k * _nxy + j * _nx + i; bb = _regionBoxes[index]; if (bb != null) regions.Add(index, bb); // layer = 0; num = bb == null ? 0 : ((HashSet)bb.Tag).Count; // Add next layer of regions - at least one while (num == 0 || layer < 1) { layer++; if (radius <= 0) { mini = i - layer; maxi = i + layer; minj = j - layer; maxj = j + layer; mink = k - layer; maxk = k + layer; } else { mini = (int)Math.Floor((point[0] - radius - _sourceBox.MinX) / _deltaX); maxi = (int)Math.Ceiling((point[0] + radius - _sourceBox.MinX) / _deltaX); minj = (int)Math.Floor((point[1] - radius - _sourceBox.MinY) / _deltaY); maxj = (int)Math.Ceiling((point[1] + radius - _sourceBox.MinY) / _deltaY); mink = (int)Math.Floor((point[2] - radius - _sourceBox.MinZ) / _deltaZ); maxk = (int)Math.Ceiling((point[2] + radius - _sourceBox.MinZ) / _deltaZ); num = 1; } // if (mini < 0) mini = 0; if (mini >= _nx) mini = _nx - 1; if (maxi < 0) maxi = 0; if (maxi >= _nx) maxi = _nx - 1; // if (minj < 0) minj = 0; if (minj >= _ny) minj = _ny - 1; if (maxj < 0) maxj = 0; if (maxj >= _ny) maxj = _ny - 1; // if (mink < 0) mink = 0; if (mink >= _nz) mink = _nz - 1; if (maxk < 0) maxk = 0; if (maxk >= _nz) maxk = _nz - 1; // for (int kk = mink; kk <= maxk; kk++) { for (int jj = minj; jj <= maxj; jj++) { for (int ii = mini; ii <= maxi; ii++) { index = kk * _nxy + jj * _nx + ii; if (!regions.ContainsKey(index)) { bb = _regionBoxes[index]; // if (bb != null && ((HashSet)bb.Tag).Count > 0) { regions.Add(index, bb); num += ((HashSet)bb.Tag).Count; } } } } } } // return regions; } // private BoundingBox ComputeAllPointsBoundingBox(CloudPoint[] cloudPoints) { _numOfValues = -1; BoundingBox bb = new BoundingBox(); bb.IncludeFirstCoor(cloudPoints[0].Coor); for (int i = 1; i < cloudPoints.Length; i++) { bb.IncludeCoorFast(cloudPoints[i].Coor); if (_numOfValues == -1) _numOfValues = cloudPoints[i].Values.Length; if (_numOfValues != cloudPoints[i].Values.Length) throw new CaeException("The interpolator does not contain the same number of values at each cloud point."); } return bb; } private BoundingBox[] AssignPointsToRegions(CloudPoint[] cloudPoints, BoundingBox sourceBox, int nx, int ny, int nz) { int nxy = nx * ny; double deltaX = sourceBox.GetXSize() / nx; double deltaY = sourceBox.GetYSize() / ny; double deltaZ = sourceBox.GetZSize() / nz; // BoundingBox bb; BoundingBox[] regions = new BoundingBox[nxy * nz]; // int pointI; int pointJ; int pointK; int regionIndex; // If cell box max value is on the border of the region division, the cell will be a member of both space regions //foreach (var cloudPoint in cloudPoints) for (int i = 0; i < cloudPoints.Length; i++) { pointI = (int)Math.Floor((cloudPoints[i].Coor[0] - sourceBox.MinX) / deltaX); pointJ = (int)Math.Floor((cloudPoints[i].Coor[1] - sourceBox.MinY) / deltaY); pointK = (int)Math.Floor((cloudPoints[i].Coor[2] - sourceBox.MinZ) / deltaZ); // if (pointI < 0) pointI = 0; else if (pointI >= _nx) pointI = _nx - 1; if (pointJ < 0) pointJ = 0; else if (pointJ >= _ny) pointJ = _ny - 1; if (pointK < 0) pointK = 0; else if (pointK >= _nz) pointK = _nz - 1; // regionIndex = pointK * nxy + pointJ * nx + pointI; bb = regions[regionIndex]; if (bb == null) { bb = new BoundingBox(); bb.MinX = sourceBox.MinX + pointI * deltaX; bb.MaxX = bb.MinX + deltaX; bb.MinY = sourceBox.MinY + pointJ * deltaY; bb.MaxY = bb.MinY + deltaY; bb.MinZ = sourceBox.MinZ + pointK * deltaZ; bb.MaxZ = bb.MinZ + deltaZ; bb.Tag = new HashSet(); regions[regionIndex] = bb; } ((HashSet)bb.Tag).Add(cloudPoints[i]); } // return regions; } } }