diff --git a/edt/edt.go b/edt/edt.go new file mode 100644 index 0000000..db2460b --- /dev/null +++ b/edt/edt.go @@ -0,0 +1,293 @@ +package edt + +import "math" + +type Bounds struct { + X int + Y int +} + +type Coordinate struct { + X int + Y int +} + +/* +Represents a Euclidean Distance Transform of dimensions Dimensions with the transform +stored in D and a lookup in L. The lookup represents the nearest occupied coordinate (x, y) +relative to the current coordinate. +*/ +type EDT struct { + Dimensions Bounds //spatial dimensions + D []float64 //distance transform + L []Coordinate //transform lookup + gx []float64 //buffer for our vector function, per row + gy []float64 //buffer for our vector function, per column +} + +func NewEDT(dimensions Bounds) *EDT { + edt := &EDT{ + Dimensions: dimensions, + } + + edt.InitializeD() + return edt +} + +func (edt *EDT) InitializeD() { + + n := edt.Dimensions.X * edt.Dimensions.Y + if n <= 0 { + return + } + + edt.D = make([]float64, n) + edt.L = make([]Coordinate, n) + edt.gx = make([]float64, edt.Dimensions.X) + edt.gy = make([]float64, edt.Dimensions.Y) + + for i := range n { + edt.D[i] = math.Inf(+1) + } +} + +/* +Reinitializes EDT using occupancy values. Occupancy dimensions must match EDT. +*/ +func (edt *EDT) AssignOccupancy(occupancy []bool) { + + n := edt.Dimensions.X * edt.Dimensions.Y + if len(occupancy) != n { + return + } + + //assign occupancy values + for i := range n { + if occupancy[i] { + edt.D[i] = 0 + } else { + edt.D[i] = math.Inf(+1) + } + } +} + +/* +EDT must already have been initialized with occupancy assigned +*/ +func (edt *EDT) ComputeDistanceTransform() { + if len(edt.D) == 0 { + return + } + //compute edt for each row + for j := range edt.Dimensions.Y { + //prepare our row vector + edt.ConstructGx(j) + + //compute distance transform for this row + //rowD, rowL := edt.DistanceTransform1D(edt.gx) + rowD, rowL := edt1D_with_labels(edt.gx) + + //update corresponding row + edt.UpdateRow(j, rowD, rowL) + } + + //compute edt for each column + for i := range edt.Dimensions.X { + //prepare our column vector + edt.ConstructGy(i) + + //compute distance transform for this column + //colD, colL := edt.DistanceTransform1D(edt.gy) + colD, colL := edt1D_with_labels(edt.gy) + + //update corresponding column + edt.UpdateCol(i, colD, colL) + } +} + +/*Computes EDT and assigns lookup coordinates for given vector g*/ +func (edt *EDT) DistanceTransform1D(g []float64) (d []float64, l []int) { + //perform 1d distance transform + n := len(g) + d = make([]float64, n) + l = make([]int, n) + + v := make([]int, n) + // envelope sites (seed indices) + z := make([]float64, n+1) // changeover abscissae + k := 0 + + // top + v[0] = 0 + z[0] = math.Inf(-1) + z[1] = math.Inf(+1) + + // initialize first valid seed + k = -1 + for q := 0; q < n; q++ { + if g[q] == math.Inf(+1) { + continue + } + + // push q: pop while it starts before current segment + var s float64 + for { + if k < 0 { + k = 0 + v[0] = q + z[0] = math.Inf(-1) + z[1] = math.Inf(+1) + break + } + i := v[k] + // intersection s(i,q) + s = ((float64(q*q) + g[q]) - (float64(i*i) + g[i])) / float64(2*(q-i)) + if s <= z[k] { + k-- + if k < 0 { + continue + } // force a push with -inf + } else { + k++ + v[k] = q + z[k] = s + z[k+1] = math.Inf(+1) + break + } + } + } + + // evaluate + k = 0 + for x := 0; x < n; x++ { + for z[k+1] <= float64(x) { + k++ + } + i := v[k] + l[x] = i + dx := float64(x - i) + d[x] = dx*dx + g[i] + } + + return d, l +} + +func edt1D_with_labels(g []float64) ([]float64, []int) { + n := len(g) + d := make([]float64, n) + l := make([]int, n) + + // --- 1. Gather valid seeds --- + type site struct { + i int + gi float64 + } + seeds := make([]site, 0, n) + for i := 0; i < n; i++ { + if !math.IsInf(g[i], +1) { // it's a seed + seeds = append(seeds, site{i, g[i]}) + } + } + + // --- 2. Handle empty case directly --- + if len(seeds) == 0 { + for i := 0; i < n; i++ { + d[i] = math.Inf(+1) + l[i] = -1 + } + return d, l + } + + // --- 3. Lower-envelope construction (Felzenszwalb–Huttenlocher) --- + v := make([]int, len(seeds)) + z := make([]float64, len(seeds)+1) + k := 0 + v[0] = seeds[0].i + z[0] = math.Inf(-1) + z[1] = math.Inf(+1) + + for q := 1; q < len(seeds); q++ { + i := seeds[q].i + gi := seeds[q].gi + for { + j := v[k] + gj := g[j] + s := ((float64(i*i) + gi) - (float64(j*j) + gj)) / float64(2*(i-j)) + if s <= z[k] { + k-- + if k < 0 { + k = 0 + break + } + continue + } + k++ + v[k] = i + z[k] = s + z[k+1] = math.Inf(+1) + break + } + } + + // --- 4. Evaluate --- + k = 0 + for x := 0; x < n; x++ { + for z[k+1] <= float64(x) { + k++ + } + i := v[k] + l[x] = i + dx := float64(x - i) + d[x] = dx*dx + g[i] + } + + return d, l +} + +/* +Prepares row vector (gx) from specified column in D +*/ +func (edt *EDT) ConstructGx(col int) { + if col < edt.Dimensions.Y { + //we can be very efficient due to slice referencing + startIdx := col * edt.Dimensions.Y + edt.gx = edt.D[startIdx : startIdx+edt.Dimensions.X] + } +} + +/* +Prepares column vector (gy) from specified row in D +*/ +func (edt *EDT) ConstructGy(row int) { + if row < edt.Dimensions.X { + for j := range edt.Dimensions.Y { + edt.gy[j] = edt.D[j*edt.Dimensions.X+row] + } + } +} + +/* +Writes data and lookup values into EDT row specified by rowidx +*/ +func (edt *EDT) UpdateRow(rowidx int, data []float64, lookup []int) { + if rowidx < edt.Dimensions.Y { + startidx := rowidx * edt.Dimensions.X + for i := range edt.Dimensions.X { + edt.D[startidx+i] = data[i] + edt.L[startidx+i].X = lookup[i] + } + } +} + +/* +Writes data and lookup values into EDT column specified by colidx +*/ +func (edt *EDT) UpdateCol(colidx int, data []float64, lookup []int) { + if colidx < edt.Dimensions.X { + for j := range edt.Dimensions.Y { + edt.D[j*edt.Dimensions.X+colidx] = data[j] + edt.L[j*edt.Dimensions.X+colidx].X = edt.L[lookup[j]*edt.Dimensions.X+colidx].X + edt.L[j*edt.Dimensions.X+colidx].Y = lookup[j] + } + } +} diff --git a/elements/flaskRoundedBottom.go b/elements/flaskRoundedBottom.go index 5a7f413..d335456 100644 --- a/elements/flaskRoundedBottom.go +++ b/elements/flaskRoundedBottom.go @@ -12,6 +12,8 @@ import ( const ( RoundedBottomFlaskWidth = 32 //pixels RoundedBottomFlaskHeight = 32 //pixels + RBFlaskMaskWidth = 46 //pixels + RBFlaskMaskHeight = 46 //pixels RoundedBottomFlaskFluidRadius = 9 //pixels: represents spherical portion of the flask where fluid will be contained RoundedBottomFlaskFluidOriginX = 16 //pixels: x origin of fluid area RoundedBottomFlaskFluidOriginY = 20 //pixels: y origin of fluid area @@ -22,27 +24,38 @@ const ( type RoundedBottomFlask struct { MappedEntityBase - fluid *fluid.Fluid //our physical representation of the fluid - fluidbuff *ebiten.Image //predraw for the fluid - fluidcellbuff *ebiten.Image //persistent fluid sprite, we redraw this everywhere we want to represent fluid - flaskbase *ebiten.Image //flask background (container) - flaskhighlight *ebiten.Image //flask foreground (glassware highlights) - fieldscale gamedata.Vector //used for transforming from fluid-space to sprite-space - angle float64 + fluid *fluid.Fluid //our physical representation of the fluid + fluidbuffS *ebiten.Image //predraw for the fluid, static + fluidbuffD *ebiten.Image //predraw for the fluid, dynamic + fluidcellbuff *ebiten.Image //persistent fluid sprite, we redraw this everywhere we want to represent fluid + flaskbase *ebiten.Image //flask background (container) + flaskhighlight *ebiten.Image //flask foreground (glassware highlights) + flaskboundarymask *ebiten.Image //flask boundary mask + fieldscaleStatic gamedata.Vector //used for transforming from fluid-space to sprite-space: static + fieldscaleDyanmic gamedata.Vector //used for transforming from fluid-space to sprite-space: dynamic + angle float64 //fluid color business fluidcolor color.RGBA //premultiplied fluid color, as set by external sources fluidcolorF []float32 //for caching of individual color values, we compute rarely + + //boundary + boundaryinitialized bool + container *fluid.Boundary } func NewRoundedBottomFlask() *RoundedBottomFlask { flask := &RoundedBottomFlask{ - fluidbuff: ebiten.NewImage(RoundedBottomFlaskFluidRadius*2, RoundedBottomFlaskFluidRadius*2), - fluidcellbuff: ebiten.NewImage(1, 1), - flaskbase: ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskBase]), - flaskhighlight: ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskHighlights]), - angle: 0, - fluidcolorF: make([]float32, 4), //one field each for R,G,B,A + fluidbuffS: ebiten.NewImage(RoundedBottomFlaskFluidRadius*2, RoundedBottomFlaskFluidRadius*2), + fluidbuffD: ebiten.NewImage(RBFlaskMaskWidth, RBFlaskMaskHeight), + fluidcellbuff: ebiten.NewImage(1, 1), + flaskbase: ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskBase]), + flaskhighlight: ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskHighlights]), + flaskboundarymask: ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskBoundaryMap46]), + //flaskboundarymask: ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskBoundaryMap]), + angle: 0, + fluidcolorF: make([]float32, 4), //one field each for R,G,B,A + boundaryinitialized: false, } flask.Initialize() return flask @@ -66,17 +79,23 @@ func (flask *RoundedBottomFlask) Initialize() { flask.fluid.Initialize() flask.fluid.Block = false //rounded flask, not a rect volume - //compute fieldscale using newly created fluid - flask.fieldscale = gamedata.Vector{ + //compute fieldscales using newly created fluid + flask.fieldscaleStatic = gamedata.Vector{ X: RoundedBottomFlaskFluidRadius * 2 / float64(flask.fluid.Field.Nx-1), Y: RoundedBottomFlaskFluidRadius * 2 / float64(flask.fluid.Field.Ny-1), } + flask.fieldscaleDyanmic = gamedata.Vector{ + X: RBFlaskMaskWidth / float64(flask.fluid.Field.Nx-2), + Y: RBFlaskMaskHeight / float64(flask.fluid.Field.Ny-2), + } + //setup default fluid color flask.SetFluidColor(color.RGBA{R: 0x0, G: 0x0, B: 0xff, A: 0xff}) } func (flask *RoundedBottomFlask) Update() { + if flask.paused { return } @@ -91,7 +110,11 @@ func (flask *RoundedBottomFlask) Draw() { flask.Sprite.DrawImage(flask.flaskbase, nil) //render fluid - flask.RenderFluid() + if flask.boundaryinitialized { + flask.RenderFluidDynamic() + } else { + flask.RenderFluidStatic() + } //render flask foreground flask.Sprite.DrawImage(flask.flaskhighlight, nil) @@ -107,8 +130,10 @@ func (flask *RoundedBottomFlask) GetAngle() float64 { return flask.angle } -func (flask *RoundedBottomFlask) RenderFluid() { - flask.fluidbuff.Clear() +func (flask *RoundedBottomFlask) RenderFluidDynamic() { + flask.fluidbuffD.Clear() + //flask.fluidbuffD.Fill(color.White) + //vector.StrokeRect(flask.fluidbuffD, 0, 0, 46, 46, 1, color.White, true) //construct fluid buffer from fluid simulation for i := range flask.fluid.Field.Nx { @@ -125,17 +150,56 @@ func (flask *RoundedBottomFlask) RenderFluid() { celldensity = 1 }*/ - ox := float64(i) * flask.fieldscale.X - oy := float64(j) * flask.fieldscale.Y + ox := float64(i) * flask.fieldscaleDyanmic.X + oy := float64(j) * flask.fieldscaleDyanmic.Y op := &ebiten.DrawImageOptions{} op.GeoM.Translate(-.5, -.5) - op.GeoM.Scale(flask.fieldscale.X, flask.fieldscale.Y) + op.GeoM.Scale(flask.fieldscaleDyanmic.X, flask.fieldscaleDyanmic.Y) op.GeoM.Translate(ox, oy) op.ColorScale.ScaleAlpha(celldensity) op.ColorScale.Scale(flask.fluidcolorF[0], flask.fluidcolorF[1], flask.fluidcolorF[2], flask.fluidcolorF[3]) - // op.ColorM.Scale(0, 0, 1, 1) - //flask.Sprite.DrawImage(flask.fluidcellbuff, op) - flask.fluidbuff.DrawImage(flask.fluidcellbuff, op) + flask.fluidbuffD.DrawImage(flask.fluidcellbuff, op) + } + } + + //transform buffer for our flask space + + s := float64(RoundedBottomFlaskWidth) / RBFlaskMaskWidth * 67. / 104 + op := &ebiten.DrawImageOptions{} + op.GeoM.Translate(-RBFlaskMaskWidth/2, -RBFlaskMaskHeight/2) + op.GeoM.Scale(s, -s*1.15) + op.GeoM.Translate(RoundedBottomFlaskWidth/2, RoundedBottomFlaskHeight/2+2) + //op.GeoM.Translate(RoundedBottomFlaskFluidOriginX, RoundedBottomFlaskFluidOriginY) + flask.Sprite.DrawImage(flask.fluidbuffD, op) +} + +func (flask *RoundedBottomFlask) RenderFluidStatic() { + flask.fluidbuffS.Clear() + + //construct fluid buffer from fluid simulation + for i := range flask.fluid.Field.Nx { + for j := range flask.fluid.Field.Ny { + + idx := i*flask.fluid.Field.Ny + j + + if flask.fluid.Field.CellType[idx] != fluid.CellTypeFluid { + continue + } + + celldensity := flask.fluid.ParticleDensity[idx] / flask.fluid.ParticleRestDensity + /*if celldensity > 0.8 { + celldensity = 1 + }*/ + + ox := float64(i) * flask.fieldscaleStatic.X + oy := float64(j) * flask.fieldscaleStatic.Y + op := &ebiten.DrawImageOptions{} + op.GeoM.Translate(-.5, -.5) + op.GeoM.Scale(flask.fieldscaleStatic.X, flask.fieldscaleStatic.Y) + op.GeoM.Translate(ox, oy) + op.ColorScale.ScaleAlpha(celldensity) + op.ColorScale.Scale(flask.fluidcolorF[0], flask.fluidcolorF[1], flask.fluidcolorF[2], flask.fluidcolorF[3]) + flask.fluidbuffS.DrawImage(flask.fluidcellbuff, op) } } @@ -145,7 +209,7 @@ func (flask *RoundedBottomFlask) RenderFluid() { op.GeoM.Translate(-RoundedBottomFlaskFluidRadius, -RoundedBottomFlaskFluidRadius) op.GeoM.Scale(1, -1) op.GeoM.Translate(RoundedBottomFlaskFluidOriginX, RoundedBottomFlaskFluidOriginY) - flask.Sprite.DrawImage(flask.fluidbuff, op) + flask.Sprite.DrawImage(flask.fluidbuffS, op) } func (flask *RoundedBottomFlask) SetFluidColor(c color.RGBA) { @@ -155,3 +219,47 @@ func (flask *RoundedBottomFlask) SetFluidColor(c color.RGBA) { flask.fluidcolorF[2] = float32(flask.fluidcolor.B) / float32(flask.fluidcolor.A) flask.fluidcolorF[3] = float32(flask.fluidcolor.A) / 0xff } + +func (flask *RoundedBottomFlask) InitializeBoundary() { + + //prepare the dimensions of our boundary map + dimensions := fluid.BoundaryDimensions{ + X: RBFlaskMaskWidth, + Y: RBFlaskMaskHeight, + } + + //instantiate the boundary structure + flask.container = fluid.NewBoundary(dimensions) + + //load pixel data from boundary + var pixels []byte = make([]byte, dimensions.X*dimensions.Y*4) + flask.flaskboundarymask.ReadPixels(pixels) + + //populate our boundary map based on the pixel data + for i := 0; i < len(pixels); i += 4 { + cellidx := i / 4 + boundary := pixels[i] == 0 + flask.container.Cells[cellidx] = !boundary + } + + //apply to fluid simulation + flask.fluid.SetBoundary(flask.container) + flask.boundaryinitialized = true +} + +// set new boundary mask and reinitialize associated data, including computing and +// passing that into the associated fluid simulation +func (flask *RoundedBottomFlask) SetBoundaryMap(img *ebiten.Image) { + flask.flaskboundarymask = img + flask.InitializeBoundary() +} + +func (flask *RoundedBottomFlask) ToggleBoundaryMask() { + if flask.boundaryinitialized { + flask.fluid.SetBoundary(nil) + flask.boundaryinitialized = false + } else { + flask.InitializeBoundary() + flask.boundaryinitialized = true + } +} diff --git a/fluid/boundary.go b/fluid/boundary.go new file mode 100644 index 0000000..0d4269f --- /dev/null +++ b/fluid/boundary.go @@ -0,0 +1,19 @@ +package fluid + +type BoundaryDimensions struct { + X int + Y int +} + +type Boundary struct { + Dimensions BoundaryDimensions + Cells []bool +} + +func NewBoundary(dimensions BoundaryDimensions) *Boundary { + b := &Boundary{ + Dimensions: dimensions, + Cells: make([]bool, dimensions.X*dimensions.Y), + } + return b +} diff --git a/fluid/fluid.go b/fluid/fluid.go index 9fd1179..ebe8d22 100644 --- a/fluid/fluid.go +++ b/fluid/fluid.go @@ -4,6 +4,7 @@ package fluid import ( + "fluids/edt" "fluids/utils" "math" ) @@ -41,6 +42,9 @@ type Fluid struct { flipPicRatio float32 numSubSteps int //number of simulation substeps Block bool //rectangular or circular field container + + //for arbitrary boundaries + edt *edt.EDT } func NewFluid(dimensions FieldVector, spacing float32) *Fluid { @@ -52,6 +56,7 @@ func NewFluid(dimensions FieldVector, spacing float32) *Fluid { flipPicRatio: FluidDefaultFlipPicRatio, numSubSteps: FluidDefaultSubSteps, Block: true, + edt: nil, } f.Initialize() @@ -268,6 +273,11 @@ func (f *Fluid) HandleParticleCollisions() { minDist2 := minDist * minDist */ + if f.edt != nil { + f.HandleBoundaryCollisions() + return + } + if f.Block { minX := f.Field.H + f.particleRadius maxX := float32(f.Field.Nx-1)*f.Field.H - f.particleRadius @@ -651,3 +661,53 @@ func (f *Fluid) SolveIncompressibility() { func (f *Fluid) ToggleShape() { f.Block = !f.Block } + +func (f *Fluid) SetBoundary(b *Boundary) { + + if b != nil { + dim := edt.Bounds{ + X: b.Dimensions.X, + Y: b.Dimensions.Y, + } + + f.edt = edt.NewEDT(dim) + f.edt.AssignOccupancy(b.Cells) + f.edt.ComputeDistanceTransform() + } else { + f.edt = nil + } +} + +// we perform our boundary resolution on particles according to our defined fluid boundary +func (f *Fluid) HandleBoundaryCollisions() { + + //grid spacing within our distance field + dx := f.Field.Dimensions.X / float32(f.edt.Dimensions.X) + dy := f.Field.Dimensions.Y / float32(f.edt.Dimensions.Y) + + //find cell of distance field for which particles belong, check if it's in bounds + for i := range f.Particles { + p := &f.Particles[i] + + xi := utils.Clamp(int(math.Floor(float64(p.Position.X/dx))), 0, f.edt.Dimensions.X-1) + yi := utils.Clamp(int(math.Floor(float64(p.Position.Y/dy))), 0, f.edt.Dimensions.Y-1) + + cellidx := xi + yi*f.edt.Dimensions.X + + //if our current cell isn't a boundary, then we skip + if f.edt.D[cellidx] == 0 { + continue + } + + //find where the particle actually belongs + pos := f.edt.L[cellidx] + newx := (float32(pos.X-xi) + f.particleRadius) * dx + newy := (float32(pos.Y-yi) + f.particleRadius) * dy + p.Position.X += newx + p.Position.Y += newy + p.Velocity.X = 0 + p.Velocity.Y = 0 + + } + +} diff --git a/game/game.go b/game/game.go index 99f767f..9c30096 100644 --- a/game/game.go +++ b/game/game.go @@ -335,6 +335,10 @@ func (g *Game) ManageFlaskInputs() { g.mfangle = g.mfangle + (angle - g.mdangle) } + + if inpututil.IsKeyJustPressed(ebiten.KeyM) { + g.flask.ToggleBoundaryMask() + } } func (g *Game) UpdateFlask() { diff --git a/resources/images.go b/resources/images.go index 1f712f9..cb8259a 100644 --- a/resources/images.go +++ b/resources/images.go @@ -13,8 +13,10 @@ import ( type ImageName string const ( - RoundedBottomFlaskBase ImageName = "RoundedBottomFlaskBase" - RoundedBottomFlaskHighlights ImageName = "RoundedBottomFlaskHighlights" + RoundedBottomFlaskBase ImageName = "RoundedBottomFlaskBase" + RoundedBottomFlaskHighlights ImageName = "RoundedBottomFlaskHighlights" + RoundedBottomFlaskBoundaryMap ImageName = "RoundedBottomFlaskBoundaryMap" + RoundedBottomFlaskBoundaryMap46 ImageName = "RoundedBottomFlaskBoundaryMap46" ) var ( @@ -23,13 +25,19 @@ var ( //go:embed rounded_bottom_flask_base.png rounded_bottom_flask_base []byte //go:embed rounded_bottom_flask_highlights.png - rounded_bottom_flask_highlitsh []byte + rounded_bottom_flask_highlight []byte + //go:embed rounded_bottom_flask_boundary_map.png + rounded_bottom_flask_boundary_map []byte + //go:embed rounded_bottom_flask_boundary_map_46.png + rounded_bottom_flask_boundary_map_46 []byte ) func LoadImages() { ImageBank = make(map[ImageName]*ebiten.Image) ImageBank[RoundedBottomFlaskBase] = LoadImagesFatal(rounded_bottom_flask_base) - ImageBank[RoundedBottomFlaskHighlights] = LoadImagesFatal(rounded_bottom_flask_highlitsh) + ImageBank[RoundedBottomFlaskHighlights] = LoadImagesFatal(rounded_bottom_flask_highlight) + ImageBank[RoundedBottomFlaskBoundaryMap] = LoadImagesFatal(rounded_bottom_flask_boundary_map) + ImageBank[RoundedBottomFlaskBoundaryMap46] = LoadImagesFatal(rounded_bottom_flask_boundary_map_46) } func LoadImagesFatal(b []byte) *ebiten.Image { diff --git a/resources/rounded_bottom_flask_boundary_map.png b/resources/rounded_bottom_flask_boundary_map.png new file mode 100644 index 0000000..862c6eb Binary files /dev/null and b/resources/rounded_bottom_flask_boundary_map.png differ diff --git a/resources/rounded_bottom_flask_boundary_map_46.png b/resources/rounded_bottom_flask_boundary_map_46.png new file mode 100644 index 0000000..4af5d1b Binary files /dev/null and b/resources/rounded_bottom_flask_boundary_map_46.png differ