diff --git a/curve/CurveUtils.cs b/curve/CurveUtils.cs index bff793cc..bc5a9db5 100644 --- a/curve/CurveUtils.cs +++ b/curve/CurveUtils.cs @@ -146,6 +146,56 @@ public static void InPlaceSmooth(IList vertices, int iStart, int iEnd, + /// + /// smooth set of vertices using extra buffer + /// + public static void IterativeSmooth(IList vertices, double alpha, int nIterations, bool bClosed) + { + IterativeSmooth(vertices, 0, vertices.Count, alpha, nIterations, bClosed); + } + /// + /// smooth set of vertices using extra buffer + /// + public static void IterativeSmooth(IList vertices, int iStart, int iEnd, double alpha, int nIterations, bool bClosed, Vector3d[] buffer = null) + { + int N = vertices.Count; + if (buffer == null || buffer.Length < N) + buffer = new Vector3d[N]; + if (bClosed) { + for (int iter = 0; iter < nIterations; ++iter) { + for (int ii = iStart; ii < iEnd; ++ii) { + int i = (ii % N); + int iPrev = (ii == 0) ? N - 1 : ii - 1; + int iNext = (ii + 1) % N; + Vector3d prev = vertices[iPrev], next = vertices[iNext]; + Vector3d c = (prev + next) * 0.5f; + buffer[i] = (1 - alpha) * vertices[i] + (alpha) * c; + } + for (int ii = iStart; ii < iEnd; ++ii) { + int i = (ii % N); + vertices[i] = buffer[i]; + } + } + } else { + for (int iter = 0; iter < nIterations; ++iter) { + for (int i = iStart; i <= iEnd; ++i) { + if (i == 0 || i >= N - 1) + continue; + Vector3d prev = vertices[i - 1], next = vertices[i + 1]; + Vector3d c = (prev + next) * 0.5f; + buffer[i] = (1 - alpha) * vertices[i] + (alpha) * c; + } + for (int ii = iStart; ii < iEnd; ++ii) { + int i = (ii % N); + vertices[i] = buffer[i]; + } + } + } + } + + + + } diff --git a/curve/CurveUtils2.cs b/curve/CurveUtils2.cs index de6729c0..ed1bf710 100644 --- a/curve/CurveUtils2.cs +++ b/curve/CurveUtils2.cs @@ -173,14 +173,25 @@ public static void LaplacianSmoothConstrained(Polygon2d poly, double alpha, int for (int i = 0; i < N; ++i ) { Vector2d curpos = poly[i]; Vector2d smoothpos = (poly[(i + N - 1) % N] + poly[(i + 1) % N]) * 0.5; - bool contained = true; - if (bAllowShrink == false || bAllowGrow == false) - contained = origPoly.Contains(smoothpos); + bool do_smooth = true; - if (bAllowShrink && contained == false) - do_smooth = false; - if (bAllowGrow && contained == true) - do_smooth = false; + if (bAllowShrink == false || bAllowGrow == false) { + bool is_inside = origPoly.Contains(smoothpos); + if (is_inside == true) + do_smooth = bAllowShrink; + else + do_smooth = bAllowGrow; + } + + // [RMS] this is old code...I think not correct? + //bool contained = true; + //if (bAllowShrink == false || bAllowGrow == false) + // contained = origPoly.Contains(smoothpos); + //bool do_smooth = true; + //if (bAllowShrink && contained == false) + // do_smooth = false; + //if (bAllowGrow && contained == true) + // do_smooth = false; if ( do_smooth ) { Vector2d newpos = beta * curpos + alpha * smoothpos; @@ -246,6 +257,37 @@ public static void Split(List objects, out List set1, out List set2, + public static Polygon2d SplitToTargetLength(Polygon2d poly, double length) + { + Polygon2d result = new Polygon2d(); + result.AppendVertex(poly[0]); + for (int j = 0; j < poly.VertexCount; ++j) { + int next = (j + 1) % poly.VertexCount; + double len = poly[j].Distance(poly[next]); + if (len < length) { + result.AppendVertex(poly[next]); + continue; + } + + int steps = (int)Math.Ceiling(len / length); + for (int k = 1; k < steps; ++k) { + double t = (double)(k) / (double)steps; + Vector2d v = (1.0 - t) * poly[j] + (t) * poly[next]; + result.AppendVertex(v); + } + + if (j < poly.VertexCount - 1) { + Util.gDevAssert(poly[j].Distance(result.Vertices[result.VertexCount - 1]) > 0.0001); + result.AppendVertex(poly[next]); + } + } + + return result; + } + + + + /// /// Remove polygons and polygon-holes smaller than minArea /// diff --git a/curve/PolyLine2d.cs b/curve/PolyLine2d.cs index 5ff1aac5..a0d201e5 100644 --- a/curve/PolyLine2d.cs +++ b/curve/PolyLine2d.cs @@ -340,6 +340,43 @@ public PolyLine2d Transform(ITransform2 xform) } + + static public PolyLine2d MakeBoxSpiral(Vector2d center, double len, double spacing) + { + double d = spacing / 2; + PolyLine2d pline = new PolyLine2d(); + pline.AppendVertex(center); + + Vector2d c = center; + c.x += spacing / 2; + pline.AppendVertex(c); + c.y += spacing; + pline.AppendVertex(c); + double accum = spacing / 2 + spacing; + + double w = spacing / 2; + double h = spacing; + + double sign = -1.0; + while (accum < len) { + w += spacing; + c.x += sign * w; + pline.AppendVertex(c); + accum += w; + + h += spacing; + c.y += sign * h; + pline.AppendVertex(c); + accum += h; + + sign *= -1.0; + } + + return pline; + } + + + } diff --git a/geometry3Sharp.csproj b/geometry3Sharp.csproj index 17d9ff96..f2318307 100644 --- a/geometry3Sharp.csproj +++ b/geometry3Sharp.csproj @@ -104,6 +104,7 @@ + @@ -294,6 +295,7 @@ + diff --git a/implicit/ImplicitFieldSampler3d.cs b/implicit/ImplicitFieldSampler3d.cs new file mode 100644 index 00000000..8f01b3fa --- /dev/null +++ b/implicit/ImplicitFieldSampler3d.cs @@ -0,0 +1,91 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace g3 +{ + /// + /// Sample implicit fields into a dense grid + /// + public class ImplicitFieldSampler3d + { + public DenseGrid3f Grid; + public double CellSize; + public Vector3d GridOrigin; + public ShiftGridIndexer3 Indexer; + public AxisAlignedBox3i GridBounds; + + public float BackgroundValue; + + + public enum CombineModes + { + DistanceMinUnion = 0 + } + public CombineModes CombineMode = CombineModes.DistanceMinUnion; + + + public ImplicitFieldSampler3d(AxisAlignedBox3d fieldBounds, double cellSize) + { + CellSize = cellSize; + GridOrigin = fieldBounds.Min; + Indexer = new ShiftGridIndexer3(GridOrigin, CellSize); + + Vector3d max = fieldBounds.Max; max += cellSize; + int ni = (int)((max.x - GridOrigin.x) / CellSize) + 1; + int nj = (int)((max.y - GridOrigin.y) / CellSize) + 1; + int nk = (int)((max.z - GridOrigin.z) / CellSize) + 1; + + GridBounds = new AxisAlignedBox3i(0, 0, 0, ni, nj, nk); + + BackgroundValue = (float)((ni + nj + nk) * CellSize); + Grid = new DenseGrid3f(ni, nj, nk, BackgroundValue); + } + + + + public DenseGridTrilinearImplicit ToImplicit() { + return new DenseGridTrilinearImplicit(Grid, GridOrigin, CellSize); + } + + + + public void Clear(float f) + { + BackgroundValue = f; + Grid.assign(BackgroundValue); + } + + + + public void Sample(BoundedImplicitFunction3d f, double expandRadius = 0) + { + AxisAlignedBox3d bounds = f.Bounds(); + + Vector3d expand = expandRadius * Vector3d.One; + Vector3i gridMin = Indexer.ToGrid(bounds.Min-expand), + gridMax = Indexer.ToGrid(bounds.Max+expand) + Vector3i.One; + gridMin = GridBounds.ClampExclusive(gridMin); + gridMax = GridBounds.ClampExclusive(gridMax); + + AxisAlignedBox3i gridbox = new AxisAlignedBox3i(gridMin, gridMax); + switch (CombineMode) { + case CombineModes.DistanceMinUnion: + sample_min(f, gridbox.IndicesInclusive()); + break; + } + } + + + void sample_min(BoundedImplicitFunction3d f, IEnumerable indices) + { + gParallel.ForEach(indices, (idx) => { + Vector3d v = Indexer.FromGrid(idx); + double d = f.Value(ref v); + Grid.set_min(ref idx, (float)d); + }); + } + + } +} diff --git a/math/AxisAlignedBox3d.cs b/math/AxisAlignedBox3d.cs index 96f1e02f..ed41045d 100644 --- a/math/AxisAlignedBox3d.cs +++ b/math/AxisAlignedBox3d.cs @@ -100,7 +100,6 @@ public Vector3d Center { get { return new Vector3d(0.5 * (Min.x + Max.x), 0.5 * (Min.y + Max.y), 0.5 * (Min.z + Max.z)); } } - public static bool operator ==(AxisAlignedBox3d a, AxisAlignedBox3d b) { return a.Min == b.Min && a.Max == b.Max; } @@ -138,6 +137,17 @@ public Vector3d Corner(int i) return new Vector3d(x, y, z); } + /// + /// Returns point on face/edge/corner. For each coord value neg==min, 0==center, pos==max + /// + public Vector3d Point(int xi, int yi, int zi) + { + double x = (xi < 0) ? Min.x : ((xi == 0) ? (0.5*(Min.x + Max.x)) : Max.x); + double y = (yi < 0) ? Min.y : ((yi == 0) ? (0.5*(Min.y + Max.y)) : Max.y); + double z = (zi < 0) ? Min.z : ((zi == 0) ? (0.5*(Min.z + Max.z)) : Max.z); + return new Vector3d(x, y, z); + } + // TODO ////! 0 == bottom-left, 1 = bottom-right, 2 == top-right, 3 == top-left diff --git a/math/AxisAlignedBox3f.cs b/math/AxisAlignedBox3f.cs index c3f8b7ee..356483ac 100644 --- a/math/AxisAlignedBox3f.cs +++ b/math/AxisAlignedBox3f.cs @@ -150,6 +150,18 @@ public Vector3f Corner(int i) } + /// + /// Returns point on face/edge/corner. For each coord value neg==min, 0==center, pos==max + /// + public Vector3f Point(int xi, int yi, int zi) + { + float x = (xi < 0) ? Min.x : ((xi == 0) ? (0.5f * (Min.x + Max.x)) : Max.x); + float y = (yi < 0) ? Min.y : ((yi == 0) ? (0.5f * (Min.y + Max.y)) : Max.y); + float z = (zi < 0) ? Min.z : ((zi == 0) ? (0.5f * (Min.z + Max.z)) : Max.z); + return new Vector3f(x, y, z); + } + + //! value is subtracted from min and added to max public void Expand(float fRadius) { diff --git a/math/IndexUtil.cs b/math/IndexUtil.cs index 62269e54..8891eb2d 100644 --- a/math/IndexUtil.cs +++ b/math/IndexUtil.cs @@ -268,6 +268,25 @@ public static void sort_indices(ref Index3i tri) + public static Vector3i ToGrid3Index(int idx, int nx, int ny) + { + int x = idx % nx; + int y = (idx / nx) % ny; + int z = idx / (nx * ny); + return new Vector3i(x, y, z); + } + + public static int ToGrid3Linear(int i, int j, int k, int nx, int ny) { + return i + nx * (j + ny * k); + } + public static int ToGrid3Linear(Vector3i ijk, int nx, int ny) { + return ijk.x + nx * (ijk.y + ny * ijk.z); + } + public static int ToGrid3Linear(ref Vector3i ijk, int nx, int ny) { + return ijk.x + nx * (ijk.y + ny * ijk.z); + } + + /// /// Filter out invalid entries in indices[] list. Will return indices itself if diff --git a/mesh/MeshEditor.cs b/mesh/MeshEditor.cs index 117d744b..8f201332 100644 --- a/mesh/MeshEditor.cs +++ b/mesh/MeshEditor.cs @@ -246,10 +246,10 @@ public virtual int[] StitchUnorderedEdges(List EdgePairs, int group_id /// vertex ordering must reslut in appropriate orientation (which is...??) /// [TODO] check and fail on bad orientation /// - public virtual int[] StitchSpan(int[] vspan1, int[] vspan2, int group_id = -1) + public virtual int[] StitchSpan(IList vspan1, IList vspan2, int group_id = -1) { - int N = vspan1.Length; - if (N != vspan2.Length) + int N = vspan1.Count; + if (N != vspan2.Count) throw new Exception("MeshEditor.StitchSpan: spans are not the same length!!"); N--; diff --git a/queries/MeshQueries.cs b/queries/MeshQueries.cs index 0789ea6f..0f56c1c2 100644 --- a/queries/MeshQueries.cs +++ b/queries/MeshQueries.cs @@ -26,11 +26,11 @@ public static DistPoint3Triangle3 TriangleDistance(DMesh3 mesh, int ti, Vector3d /// Find point-normal(Z) frame at closest point to queryPoint on mesh. /// Returns interpolated vertex-normal frame if available, otherwise tri-normal frame. /// - public static Frame3f NearestPointFrame(DMesh3 mesh, ISpatial spatial, Vector3d queryPoint) + public static Frame3f NearestPointFrame(DMesh3 mesh, ISpatial spatial, Vector3d queryPoint, bool bForceFaceNormal = false) { int tid = spatial.FindNearestTriangle(queryPoint); Vector3d surfPt = TriangleDistance(mesh, tid, queryPoint).TriangleClosest; - if (mesh.HasVertexNormals) + if (mesh.HasVertexNormals && bForceFaceNormal == false) return SurfaceFrame(mesh, tid, surfPt); else return new Frame3f(surfPt, mesh.GetTriNormal(tid)); @@ -113,7 +113,7 @@ public static DistTriangle3Triangle3 TrianglesDistance(DMesh3 mesh1, int ti, DMe /// Find point-normal frame at ray-intersection point on mesh, or return false if no hit. /// Returns interpolated vertex-normal frame if available, otherwise tri-normal frame. /// - public static bool RayHitPointFrame(DMesh3 mesh, ISpatial spatial, Ray3d ray, out Frame3f hitPosFrame) + public static bool RayHitPointFrame(DMesh3 mesh, ISpatial spatial, Ray3d ray, out Frame3f hitPosFrame, bool bForceFaceNormal = false) { hitPosFrame = new Frame3f(); int tid = spatial.FindNearestHitTriangle(ray); @@ -123,7 +123,7 @@ public static bool RayHitPointFrame(DMesh3 mesh, ISpatial spatial, Ray3d ray, ou if (isect.Result != IntersectionResult.Intersects) return false; Vector3d surfPt = ray.PointAt(isect.RayParameter); - if (mesh.HasVertexNormals) + if (mesh.HasVertexNormals && bForceFaceNormal == false) hitPosFrame = SurfaceFrame(mesh, tid, surfPt); // TODO isect has bary-coords already!! else hitPosFrame = new Frame3f(surfPt, mesh.GetTriNormal(tid)); @@ -135,7 +135,7 @@ public static bool RayHitPointFrame(DMesh3 mesh, ISpatial spatial, Ray3d ray, ou /// Get point-normal frame on surface of mesh. Assumption is that point lies in tID. /// returns interpolated vertex-normal frame if available, otherwise tri-normal frame. /// - public static Frame3f SurfaceFrame(DMesh3 mesh, int tID, Vector3d point) + public static Frame3f SurfaceFrame(DMesh3 mesh, int tID, Vector3d point, bool bForceFaceNormal = false) { if (!mesh.IsTriangle(tID)) throw new Exception("MeshQueries.SurfaceFrame: triangle " + tID + " does not exist!"); @@ -143,7 +143,7 @@ public static Frame3f SurfaceFrame(DMesh3 mesh, int tID, Vector3d point) mesh.GetTriVertices(tID, ref tri.V0, ref tri.V1, ref tri.V2); Vector3d bary = tri.BarycentricCoords(point); point = tri.PointAt(bary); - if (mesh.HasVertexNormals) { + if (mesh.HasVertexNormals && bForceFaceNormal == false) { Vector3d normal = mesh.GetTriBaryNormal(tID, bary.x, bary.y, bary.z); return new Frame3f(point, normal); } else diff --git a/spatial/DenseGrid2.cs b/spatial/DenseGrid2.cs new file mode 100644 index 00000000..84eefdc1 --- /dev/null +++ b/spatial/DenseGrid2.cs @@ -0,0 +1,280 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + + +namespace g3 +{ + + /// + /// 2D dense grid of floating-point scalar values. + /// + public class DenseGrid2f + { + public float[] Buffer; + public int ni, nj; + + public DenseGrid2f() + { + ni = nj = 0; + } + + public DenseGrid2f(int ni, int nj, float initialValue) + { + resize(ni, nj); + assign(initialValue); + } + + public DenseGrid2f(DenseGrid2f copy) + { + Buffer = new float[copy.Buffer.Length]; + Array.Copy(copy.Buffer, Buffer, Buffer.Length); + ni = copy.ni; nj = copy.nj; + } + + public void swap(DenseGrid2f g2) + { + Util.gDevAssert(ni == g2.ni && nj == g2.nj); + var tmp = g2.Buffer; + g2.Buffer = this.Buffer; + this.Buffer = tmp; + } + + public int size { get { return ni * nj; } } + + public void resize(int ni, int nj) + { + Buffer = new float[ni * nj]; + this.ni = ni; this.nj = nj; + } + + public void assign(float value) + { + for (int i = 0; i < Buffer.Length; ++i) + Buffer[i] = value; + } + + public void assign_border(float value, int rings) + { + for ( int j = 0; j < rings; ++j ) { + int jb = nj - 1 - j; + for ( int i = 0; i < ni; ++i ) { + Buffer[i + ni * j] = value; + Buffer[i + ni * jb] = value; + } + } + int stop = nj - 1 - rings; + for ( int j = rings; j < stop; ++j ) { + for ( int i = 0; i < rings; ++i ) { + Buffer[i + ni * j] = value; + Buffer[(ni - 1 - i) + ni * j] = value; + } + } + } + + + public void clear() { + Array.Clear(Buffer, 0, Buffer.Length); + } + + public void copy(DenseGrid2f copy) + { + Array.Copy(copy.Buffer, this.Buffer, this.Buffer.Length); + } + + public float this[int i] { + get { return Buffer[i]; } + set { Buffer[i] = value; } + } + + public float this[int i, int j] { + get { return Buffer[i + ni * j]; } + set { Buffer[i + ni * j] = value; } + } + + public float this[Vector2i ijk] { + get { return Buffer[ijk.x + ni * ijk.y]; } + set { Buffer[ijk.x + ni * ijk.y] = value; } + } + + public void get_x_pair(int i0, int j, out double a, out double b) + { + int offset = ni * j; + a = Buffer[offset + i0]; + b = Buffer[offset + i0 + 1]; + } + + public void apply(Func f) + { + for (int j = 0; j < nj; j++ ) { + for ( int i = 0; i < ni; i++ ) { + int idx = i + ni * j; + Buffer[idx] = f(Buffer[idx]); + } + } + } + + public void set_min(DenseGrid2f grid2) + { + for (int k = 0; k < Buffer.Length; ++k) + Buffer[k] = Math.Min(Buffer[k], grid2.Buffer[k]); + } + public void set_max(DenseGrid2f grid2) + { + for (int k = 0; k < Buffer.Length; ++k) + Buffer[k] = Math.Max(Buffer[k], grid2.Buffer[k]); + } + + public AxisAlignedBox2i Bounds { + get { return new AxisAlignedBox2i(0, 0, ni, nj); } + } + + + public IEnumerable Indices() + { + for (int y = 0; y < nj; ++y) { + for (int x = 0; x < ni; ++x) + yield return new Vector2i(x, y); + } + } + + + public IEnumerable InsetIndices(int border_width) + { + int stopy = nj - border_width, stopx = ni - border_width; + for (int y = border_width; y < stopy; ++y) { + for (int x = border_width; x < stopx; ++x) + yield return new Vector2i(x, y); + } + } + + + } + + + + + + + /// + /// 2D dense grid of integers. + /// + public class DenseGrid2i + { + public int[] Buffer; + public int ni, nj; + + public DenseGrid2i() + { + ni = nj = 0; + } + + public DenseGrid2i(int ni, int nj, int initialValue) + { + resize(ni, nj); + assign(initialValue); + } + + public DenseGrid2i(DenseGrid2i copy) + { + resize(copy.ni, copy.nj); + Array.Copy(copy.Buffer, this.Buffer, this.Buffer.Length); + } + + public int size { get { return ni * nj; } } + + public void resize(int ni, int nj) + { + Buffer = new int[ni * nj]; + this.ni = ni; this.nj = nj; + } + + public void clear() { + Array.Clear(Buffer, 0, Buffer.Length); + } + + + public void copy(DenseGrid2i copy) + { + Array.Copy(copy.Buffer, this.Buffer, this.Buffer.Length); + } + + public void assign(int value) + { + for (int i = 0; i < Buffer.Length; ++i) + Buffer[i] = value; + } + + public int this[int i] { + get { return Buffer[i]; } + set { Buffer[i] = value; } + } + + public int this[int i, int j] { + get { return Buffer[i + ni * j]; } + set { Buffer[i + ni * j] = value; } + } + + public int this[Vector2i ijk] { + get { return Buffer[ijk.x + ni * ijk.y]; } + set { Buffer[ijk.x + ni * ijk.y] = value; } + } + + public void increment(int i, int j) + { + Buffer[i + ni * j]++; + } + public void decrement(int i, int j) + { + Buffer[i + ni * j]--; + } + + public void atomic_increment(int i, int j) + { + System.Threading.Interlocked.Increment(ref Buffer[i + ni * j]); + } + + public void atomic_decrement(int i, int j) + { + System.Threading.Interlocked.Decrement(ref Buffer[i + ni * j]); + } + + public void atomic_incdec(int i, int j, bool decrement = false) { + if ( decrement ) + System.Threading.Interlocked.Decrement(ref Buffer[i + ni * j]); + else + System.Threading.Interlocked.Increment(ref Buffer[i + ni * j]); + } + + public int sum() { + int sum = 0; + for (int i = 0; i < Buffer.Length; ++i) + sum += Buffer[i]; + return sum; + } + + + public IEnumerable Indices() + { + for (int y = 0; y < nj; ++y) { + for (int x = 0; x < ni; ++x) + yield return new Vector2i(x, y); + } + } + + + public IEnumerable InsetIndices(int border_width) + { + int stopy = nj - border_width, stopx = ni - border_width; + for (int y = border_width; y < stopy; ++y) { + for (int x = border_width; x < stopx; ++x) + yield return new Vector2i(x, y); + } + } + + } + + + +} diff --git a/spatial/DenseGrid3.cs b/spatial/DenseGrid3.cs index 6ce6041d..cbe6050c 100644 --- a/spatial/DenseGrid3.cs +++ b/spatial/DenseGrid3.cs @@ -6,12 +6,42 @@ namespace g3 { - + /// + /// 3D dense grid of floating-point scalar values. + /// public class DenseGrid3f { public float[] Buffer; public int ni, nj, nk; + public DenseGrid3f() + { + ni = nj = nk = 0; + } + + public DenseGrid3f(int ni, int nj, int nk, float initialValue) + { + resize(ni, nj, nk); + assign(initialValue); + } + + public DenseGrid3f(DenseGrid3f copy) + { + Buffer = new float[copy.Buffer.Length]; + Array.Copy(copy.Buffer, Buffer, Buffer.Length); + ni = copy.ni; nj = copy.nj; nk = copy.nk; + } + + public void swap(DenseGrid3f g2) + { + Util.gDevAssert(ni == g2.ni && nj == g2.nj && nk == g2.nk); + var tmp = g2.Buffer; + g2.Buffer = this.Buffer; + this.Buffer = tmp; + } + + public int size { get { return ni * nj * nk; } } + public void resize(int ni, int nj, int nk) { Buffer = new float[ni * nj * nk]; @@ -24,6 +54,24 @@ public void assign(float value) Buffer[i] = value; } + public void set_min(ref Vector3i ijk, float f) + { + int idx = ijk.x + ni * (ijk.y + nj * ijk.z); + if (f < Buffer[idx]) + Buffer[idx] = f; + } + public void set_max(ref Vector3i ijk, float f) + { + int idx = ijk.x + ni * (ijk.y + nj * ijk.z); + if (f > Buffer[idx]) + Buffer[idx] = f; + } + + public float this[int i] { + get { return Buffer[i]; } + set { Buffer[i] = value; } + } + public float this[int i, int j, int k] { get { return Buffer[i + ni * (j + nj * k)]; } set { Buffer[i + ni * (j + nj * k)] = value; } @@ -53,6 +101,49 @@ public void apply(Func f) } } + + public DenseGrid2f get_slice(int slice_i, int dimension) + { + DenseGrid2f slice; + if (dimension == 0) { + slice = new DenseGrid2f(nj, nk, 0); + for (int k = 0; k < nk; ++k) + for (int j = 0; j < nj; ++j) + slice[j, k] = Buffer[slice_i + ni * (j + nj * k)]; + } else if (dimension == 1) { + slice = new DenseGrid2f(ni, nk, 0); + for (int k = 0; k < nk; ++k) + for (int i = 0; i < ni; ++i) + slice[i, k] = Buffer[i + ni * (slice_i + nj * k)]; + } else { + slice = new DenseGrid2f(ni, nj, 0); + for (int j = 0; j < nj; ++j) + for (int i = 0; i < ni; ++i) + slice[i, j] = Buffer[i + ni * (j + nj * slice_i)]; + } + return slice; + } + + + public void set_slice(DenseGrid2f slice, int slice_i, int dimension) + { + if (dimension == 0) { + for (int k = 0; k < nk; ++k) + for (int j = 0; j < nj; ++j) + Buffer[slice_i + ni * (j + nj * k)] = slice[j, k]; + } else if (dimension == 1) { + for (int k = 0; k < nk; ++k) + for (int i = 0; i < ni; ++i) + Buffer[i + ni * (slice_i + nj * k)] = slice[i, k]; + } else { + for (int j = 0; j < nj; ++j) + for (int i = 0; i < ni; ++i) + Buffer[i + ni * (j + nj * slice_i)] = slice[i, j]; + } + } + + + public AxisAlignedBox3i Bounds { get { return new AxisAlignedBox3i(0, 0, 0, ni, nj, nk); } } @@ -81,26 +172,53 @@ public IEnumerable InsetIndices(int border_width) } + public Vector3i to_index(int idx) { + int x = idx % ni; + int y = (idx / ni) % nj; + int z = idx / (ni * nj); + return new Vector3i(x, y, z); + } + public int to_linear(int i, int j, int k) + { + return i + ni * (j + nj * k); + } + public int to_linear(ref Vector3i ijk) + { + return ijk.x + ni * (ijk.y + nj * ijk.z); + } + public int to_linear(Vector3i ijk) + { + return ijk.x + ni * (ijk.y + nj * ijk.z); + } } - + /// + /// 3D dense grid of integers. + /// public class DenseGrid3i { public int[] Buffer; public int ni, nj, nk; + public DenseGrid3i() + { + ni = nj = nk = 0; + } + public DenseGrid3i(int ni, int nj, int nk, int initialValue) { resize(ni, nj, nk); assign(initialValue); } + public int size { get { return ni * nj * nk; } } + public void resize(int ni, int nj, int nk) { Buffer = new int[ni * nj * nk]; @@ -113,11 +231,21 @@ public void assign(int value) Buffer[i] = value; } + public int this[int i] { + get { return Buffer[i]; } + set { Buffer[i] = value; } + } + public int this[int i, int j, int k] { get { return Buffer[i + ni * (j + nj * k)]; } set { Buffer[i + ni * (j + nj * k)] = value; } } + public int this[Vector3i ijk] { + get { return Buffer[ijk.x + ni * (ijk.y + nj * ijk.z)]; } + set { Buffer[ijk.x + ni * (ijk.y + nj * ijk.z)] = value; } + } + public void increment(int i, int j, int k) { Buffer[i + ni * (j + nj * k)]++; @@ -143,6 +271,69 @@ public void atomic_incdec(int i, int j, int k, bool decrement = false) { else System.Threading.Interlocked.Increment(ref Buffer[i + ni * (j + nj * k)]); } + + + + public DenseGrid2i get_slice(int slice_i, int dimension) + { + DenseGrid2i slice; + if ( dimension == 0 ) { + slice = new DenseGrid2i(nj, nk, 0); + for (int k = 0; k < nk; ++k) + for (int j = 0; j < nj; ++j) + slice[j, k] = Buffer[slice_i + ni * (j + nj * k)]; + } else if (dimension == 1) { + slice = new DenseGrid2i(ni, nk, 0); + for (int k = 0; k < nk; ++k) + for (int i = 0; i < ni; ++i) + slice[i, k] = Buffer[i + ni * (slice_i + nj * k)]; + } else { + slice = new DenseGrid2i(ni, nj, 0); + for (int j = 0; j < nj; ++j) + for (int i = 0; i < ni; ++i) + slice[i, j] = Buffer[i + ni * (j + nj * slice_i)]; + } + return slice; + } + + + /// + /// convert to binary bitmap + /// + public Bitmap3 get_bitmap(int thresh = 0) + { + Bitmap3 bmp = new Bitmap3(new Vector3i(ni, nj, nk)); + for (int i = 0; i < Buffer.Length; ++i) + bmp[i] = (Buffer[i] > thresh) ? true : false; + return bmp; + } + + + public IEnumerable Indices() + { + for (int z = 0; z < nk; ++z) { + for (int y = 0; y < nj; ++y) { + for (int x = 0; x < ni; ++x) + yield return new Vector3i(x, y, z); + } + } + } + + + public IEnumerable InsetIndices(int border_width) + { + int stopy = nj - border_width, stopx = ni - border_width; + for (int z = border_width; z < nk - border_width; ++z) { + for (int y = border_width; y < stopy; ++y) { + for (int x = border_width; x < stopx; ++x) + yield return new Vector3i(x, y, z); + } + } + } + + + + }