Created
February 2, 2024 15:21
-
-
Save danlooo/758b40019bda15a84de97b15672d049b to your computer and use it in GitHub Desktop.
View sattelite data in Makie.jl: Globe at lower zoom levels, only load data required for the current camera view
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using GeometryBasics, LinearAlgebra, GLMakie, FileIO, CoordinateTransformations, ColorSchemes | |
function intersect_line_sphere(p::Point3f, d::Point3f, c::Point3f, sphere_radius::Float32) | |
# Use indexing to access point/vector components | |
# Calculate coefficients of the quadratic equation (a*t^2 + b*t + c = 0) | |
a = dot(d, d) | |
b = 2 * dot(d, p .- c) | |
c = dot(p .- c, p .- c) .- sphere_radius^2 | |
# Calculate discriminant | |
discriminant = b^2 - 4 * a * c | |
if discriminant < 0 | |
return nothing | |
elseif discriminant == 0 | |
# One intersection point (tangent) | |
t = -b / (2 * a) | |
intersection_point = p + t * d | |
return Point3f(intersection_point...) | |
else | |
# Two intersection points | |
t1 = (-b + sqrt(discriminant)) / (2 * a) | |
t2 = (-b - sqrt(discriminant)) / (2 * a) | |
intersection_point1 = p + t1 * d | |
intersection_point2 = p + t2 * d | |
return Point3f(intersection_point1...), Point3f(intersection_point2...) | |
end | |
end | |
function get_frustrum(cam) | |
rect_ps = [ | |
Point3f(1, 1, 1), Point3f(-1, 1, 1), | |
Point3f(1, -1, 1), Point3f(-1, -1, 1)] | |
inv_pv = inv(cam.projectionview[]) | |
return map(rect_ps) do p | |
p = inv_pv * to_ndim(Point4f, p, 1) | |
return p[Vec(1, 2, 3)] / p[4] | |
end | |
end | |
fig = Figure() | |
globe_scene = LScene(fig[1, 1]; show_axis=false) | |
globe_scene.scene.clear = true | |
globe_cam = Makie.Camera3D( | |
globe_scene.scene, | |
show_axis=false, | |
projectiontype=Makie.Perspective, | |
clipping_mode=:static, | |
# disable translation to prevent skewing | |
mouse_translationspeed=0, | |
keyboard_translationspeed=0 | |
) | |
globe_cam.near[] = 0.0001f0 | |
earth_texture = load(Makie.assetpath("earth.png")) | |
earth = mesh!(globe_scene, Sphere(Point3f(0), 1.0f0), color=earth_texture) | |
center!(globe_scene.scene) | |
cam = globe_scene.scene.camera | |
eyeposition = globe_cam.eyeposition | |
lookat = globe_cam.lookat | |
frustrum = map(pv -> get_frustrum(cam), cam.projectionview) | |
cam_scene = LScene(fig[1, 2]) | |
center!(cam_scene.scene) | |
meshscatter!(cam_scene, frustrum, color=:blue) | |
intersections = map(frustrum, eyeposition) do frustrum, eypos | |
result = Point3f[] | |
for point in frustrum | |
res = intersect_line_sphere(Point(eypos), Point(eypos .- point), Point3f(0), 1.0f0) | |
if !isnothing(res) | |
if res isa Tuple | |
res = sort!([res...], by=(x -> norm(x - eypos)))[1] | |
end | |
push!(result, res) | |
end | |
end | |
return result | |
end | |
visible = map(intersections) do intersections | |
if length(intersections) == 4 | |
earth.visible = false | |
else | |
earth.visible = true | |
end | |
end | |
plane = map(intersections) do points | |
if length(points) == 4 | |
return GeometryBasics.Mesh(meta(points; | |
normals=fill(Vec3f(0), 4), | |
uv=[Vec2f(0, 1), Vec2f(1, 1), Vec2f(0, 0), Vec2f(1, 0)] | |
), | |
[GLTriangleFace(1, 2, 3), GLTriangleFace(2, 3, 4)]) | |
else | |
return GeometryBasics.Mesh(meta([Point3f(0)]; normals=[Vec3f(0)], uv=[Vec2f(NaN)]), GLTriangleFace[]) | |
end | |
end | |
sphere_to_cart_trans = SphericalFromCartesian() | |
plane_bbox = map(intersections) do points | |
if length(points) == 4 | |
# move date line from meridian to common date line | |
p = sphere_to_cart_trans.(points) .|> x -> (x.θ > 0 ? x.θ - π : x.θ + π, x.ϕ) .|> rad2deg | |
return p | |
else | |
return [(-180.0, 180.0), (180.0, -180.0), (-180.0, -180.0), (180.0, 180.0)] | |
end | |
end | |
# TODO: Replace mock texture with actual data inside the calculated bounding box | |
s = 100 | |
values = reshape(range(1, 255, s * s), (s, s)) .|> round .|> Integer | |
plane_texture = map(x -> ColorSchemes.viridis[x], values) | |
scatter!(cam_scene, eyeposition, color=:black) | |
meshscatter!(cam_scene, intersections, color=:white) | |
wireframe!(cam_scene, Sphere(Point3f(0), 1.0f0)) | |
mesh!(cam_scene, plane, color=plane_texture, shading=NoShading) | |
mesh!(globe_scene, plane, color=plane_texture, shading=NoShading) | |
scatter(fig[1, 3], plane_bbox; axis=(ylabel="latitude [°]", xlabel="longitude [°]", title="Bounding Box")) | |
fig |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
The frustum forms a plane that intersects with the sphere resulting into a circle map. However, the viewport is rectangular. We can fill corner points with those being even further away e.g. using an orthographic projection.