1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
/*
* // Copyright (c) 2021 Feng Yang
* //
* // I am making my contributions/submissions to this project solely in my
* // personal capacity and am not conveying any rights to any intellectual
* // property of any third parties.
*/
/// Returns true if \p phi is inside the implicit surface (< 0).
/// - Parameter phi: The level set value.
/// - Returns: True if inside the implicit surface, false otherwise.
pub fn is_inside_sdf(phi: f64) -> bool {
return phi < 0.0;
}
/// Returns smeared Heaviside function.
///
/// This function returns smeared (or smooth) Heaviside (or step) function
/// between 0 and 1. If \p phi is less than -1.5, it will return 0. If \p phi
/// is greater than 1.5, it will return 1. Between -1.5 and 1.5, the function
/// will return smooth profile between 0 and 1. Derivative of this function is
/// smeared_delta_sdf.
/// - Parameter phi: The level set value.
/// - Returns: Smeared Heaviside function.
pub fn smeared_heaviside_sdf(phi: f64) -> f64 {
return if phi > 1.5 {
1.0
} else {
if phi < -1.5 {
0.0
} else {
0.5 + phi / 3.0 + 0.5 * crate::constants::K_INV_PI_D * f64::sin(crate::constants::K_PI_D * phi / 1.5)
}
};
}
/// Returns smeared delta function.
///
/// This function returns smeared (or smooth) delta function between 0 and 1.
/// If \p phi is less than -1.5, it will return 0. If \p phi is greater than
/// 1.5, it will also return 0. Between -1.5 and 1.5, the function will return
/// smooth delta function. Integral of this function is smeared_heaviside_sdf.
/// - Parameter phi: The level set value.
/// - Returns: Smeared delta function.
pub fn smeared_delta_sdf(phi: f64) -> f64 {
return if f64::abs(phi) > 1.5 {
0.0
} else {
1.0 / 3.0 + 1.0 / 3.0 * f64::cos(crate::constants::K_PI_D * phi / 1.5)
};
}
/// Returns the fraction occupied by the implicit surface.
///
/// The input parameters, \p phi0 and \p phi1, are the level set values,
/// measured from two nearby points. This function computes how much the
/// implicit surface occupies the line between two points. For example, if both
/// \p phi0 and \p phi1 are negative, it means the points are both inside the
/// surface, thus the function will return 1. If both are positive, it will
/// return 0 because both are outside the surface. If the signs are different,
/// then only one of the points is inside the surface and the function will
/// return a value between 0 and 1.
/// - Parameters:
/// - phi0: The level set value from the first point.
/// - phi1: The level set value from the second point.
/// - Returns: The fraction occupied by the implicit surface.
pub fn fraction_inside_sdf(phi0: f64, phi1: f64) -> f64 {
return if is_inside_sdf(phi0) && is_inside_sdf(phi1) {
1.0
} else if is_inside_sdf(phi0) && !is_inside_sdf(phi1) {
phi0 / (phi0 - phi1)
} else if !is_inside_sdf(phi0) && is_inside_sdf(phi1) {
phi1 / (phi1 - phi0)
} else {
0.0
};
}
fn cycle_array(arr: &mut [f64; 4]) {
let t = arr[0];
for i in 0..arr.len() {
arr[i] = arr[i + 1]
}
if let Some(last) = arr.last_mut() {
*last = t;
}
}
/// Returns the fraction occupied by the implicit surface.
///
/// Given four signed distance values (square corners), determine what fraction
/// of the square is "inside". The original implementation can be found from
/// Christopher Batty's variational fluid code at
/// https://github.com/christopherbatty/Fluid3D.
/// - Parameters:
/// - phi_bottom_left: The level set value on the bottom-left corner.
/// - phi_bottom_right: The level set value on the bottom-right corner.
/// - phi_top_left: The level set value on the top-left corner.
/// - phi_top_right: The level set value on the top-right corner.
/// - Returns: The fraction occupied by the implicit surface.
pub fn fraction_inside(phi_bottom_left: f64, phi_bottom_right: f64,
phi_top_left: f64, phi_top_right: f64) -> f64 {
let inside_count = match phi_bottom_left < 0.0 {
true => 1,
false => 0
} + match phi_top_left < 0.0 {
true => 1,
false => 0
} + match phi_bottom_right < 0.0 {
true => 1,
false => 0
} + match phi_top_right < 0.0 {
true => 1,
false => 0
};
let mut list = [phi_bottom_left, phi_bottom_right, phi_top_right, phi_top_left];
return if inside_count == 4 {
1.0
} else if inside_count == 3 {
// rotate until the positive value is in the first position
while list[0] < 0.0 {
cycle_array(&mut list);
}
// Work out the area of the exterior triangle
let side0 = 1.0 - fraction_inside_sdf(list[0], list[3]);
let side1 = 1.0 - fraction_inside_sdf(list[0], list[1]);
1.0 - 0.5 * side0 * side1
} else if inside_count == 2 {
// rotate until a negative value is in the first position, and the next
// negative is in either slot 1 or 2.
while list[0] >= 0.0 || !(list[1] < 0.0 || list[2] < 0.0) {
cycle_array(&mut list);
}
if list[1] < 0.0 { // the matching signs are adjacent
let side_left = fraction_inside_sdf(list[0], list[3]);
let side_right = fraction_inside_sdf(list[1], list[2]);
0.5 * (side_left + side_right)
} else { // matching signs are diagonally opposite
// determine the centre point's sign to disambiguate this case
let middle_point = 0.25 * (list[0] + list[1] + list[2] + list[3]);
if middle_point < 0.0 {
let mut area: f64 = 0.0;
// first triangle (top left)
let side1 = 1.0 - fraction_inside_sdf(list[0], list[3]);
let side3 = 1.0 - fraction_inside_sdf(list[2], list[3]);
area += 0.5 * side1 * side3;
// second triangle (top right)
let side2 = 1.0 - fraction_inside_sdf(list[2], list[1]);
let side0 = 1.0 - fraction_inside_sdf(list[0], list[1]);
area += 0.5 * side0 * side2;
1.0 - area
} else {
let mut area: f64 = 0.0;
// first triangle (bottom left)
let side0 = fraction_inside_sdf(list[0], list[1]);
let side1 = fraction_inside_sdf(list[0], list[3]);
area += 0.5 * side0 * side1;
// second triangle (top right)
let side2 = fraction_inside_sdf(list[2], list[1]);
let side3 = fraction_inside_sdf(list[2], list[3]);
area += 0.5 * side2 * side3;
area
}
}
} else if inside_count == 1 {
// rotate until the negative value is in the first position
while list[0] >= 0.0 {
cycle_array(&mut list);
}
// Work out the area of the interior triangle, and subtract from 1.
let side0 = fraction_inside_sdf(list[0], list[3]);
let side1 = fraction_inside_sdf(list[0], list[1]);
0.5 * side0 * side1
} else {
0.0
};
}
pub fn distance_to_zero_level_set(phi0: f64, phi1: f64) -> f64 {
return if f64::abs(phi0) + f64::abs(phi1) > f64::EPSILON {
f64::abs(phi0) / (f64::abs(phi0) + f64::abs(phi1))
} else {
0.5
};
}