From 81529de3657e8853576cc3bc096241aff3ee8d0a Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 11 Oct 2024 16:19:23 +0100 Subject: [PATCH] Add a wrapper function to dela with the replicants and other FV. (#1554) * Add a center option to solids and export them * Add a gmt_illuminate() wrapper. * Type annotation. * Type annotation. * Add a "replicants" function to replicate several copies of a 3D FV and paint them. * Change a solids test. * Add a wrapper function to dela with the replicants and other FV. FINALY MAY HAVE FOUND AND FIXED THE HEISENBUG. --- src/GMT.jl | 2 + src/choropleth_utils.jl | 10 ++-- src/common_options.jl | 4 +- src/custom_symb_funs.jl | 2 +- src/libgmt.jl | 14 ++++- src/psxy.jl | 96 +++++++++++++++++++++++++++++---- src/solids.jl | 116 ++++++++++++++++++++++++++-------------- src/utils_types.jl | 8 +-- test/test_solidos.jl | 2 +- 9 files changed, 193 insertions(+), 61 deletions(-) diff --git a/src/GMT.jl b/src/GMT.jl index 75eb9b2cf..f315b6ba2 100644 --- a/src/GMT.jl +++ b/src/GMT.jl @@ -184,6 +184,8 @@ export lazinfo, lazread, lazwrite, lasread, laswrite, + cube, dodecahedron, icosahedron, sphere, octahedron, tetrahedron, + df2ds, ODE2ds, sprintf, @?, @dir diff --git a/src/choropleth_utils.jl b/src/choropleth_utils.jl index a6a399ccf..5025a46ea 100644 --- a/src/choropleth_utils.jl +++ b/src/choropleth_utils.jl @@ -35,13 +35,15 @@ function cpt4dcw(codes::Vector{<:AbstractString}, vals::Vector{<:Real}; kwargs.. P::Ptr{GMT_PALETTE} = palette_init(G_API[1], C) # A pointer to a GMT CPT Ccat::GMTcpt = gmt("makecpt -T"*join(_codes, ",")) - rgb = [0.0, 0.0, 0.0]; + rgb = [0.0, 0.0, 0.0, 0.0]; + _rgb = [0.0, 0.0, 0.0, 0.0]; inc = (C.minmax[2] - C.minmax[1]) / size(Ccat.cpt, 1) for k = 1:size(Ccat.cpt, 1) - @GC.preserve C gmt_get_rgb_from_z(G_API[1], P, _vals[k], rgb) - [Ccat.colormap[k, n] = copy(rgb[n]) for n = 1:3] - [Ccat.cpt[k, n+1] = copy(rgb[(n % 3) + 1]) for n = 0:5] # cpt = [rgb rgb] + gmt_get_rgb_from_z(G_API[1], P, _vals[k], rgb) + Ccat.colormap[k, 1], Ccat.colormap[k, 2], Ccat.colormap[k, 3] = rgb[1], rgb[2], rgb[3] + _rgb[1], _rgb[2], _rgb[3] = rgb[1], rgb[2], rgb[3] + [Ccat.cpt[k, n+1] = copy(_rgb[(n % 3) + 1]) for n = 0:5] # cpt = [rgb rgb] Ccat.range[k,1] = C.minmax[1] + (k-1) * inc Ccat.range[k,2] = C.minmax[1] + k * inc end diff --git a/src/common_options.jl b/src/common_options.jl index 7c708f063..5b5d782a0 100644 --- a/src/common_options.jl +++ b/src/common_options.jl @@ -4664,13 +4664,13 @@ function put_in_legend_bag(d::Dict, cmd, arg, O::Bool=false, opt_l::String="") extra_opt *= ((t = scan_opt(cmd_[1], "-G", true)) != "") ? t : "" for k = 1:numel(pens) pens[k] *= extra_opt end if ((ind = findfirst(arg.colnames .== "Zcolor")) !== nothing) - rgb = [0.0, 0.0, 0.0] + rgb = [0.0, 0.0, 0.0, 0.0] P::Ptr{GMT_PALETTE} = palette_init(G_API[1], CURRENT_CPT[1]) # A pointer to a GMT CPT gmt_get_rgb_from_z(G_API[1], P, arg[gindex[1],ind], rgb) cmd_[1] *= " -G" * arg2str(rgb.*255) for k = 1:numel(pens) gmt_get_rgb_from_z(G_API[1], P, arg[gindex[k+1],ind]+10eps(), rgb) - pens[k] *= " -G" * arg2str(rgb.*255) + pens[k] *= @sprintf("-G%.0f/%.0f/%.0f", rgb[1]*255, rgb[2]*255, rgb[3]*255) end end end diff --git a/src/custom_symb_funs.jl b/src/custom_symb_funs.jl index 5764dbc3b..2c5b030fd 100644 --- a/src/custom_symb_funs.jl +++ b/src/custom_symb_funs.jl @@ -51,7 +51,7 @@ function helper_cusymb(cache, name, fmt, format, subdir="") fmt = ".eps" end (fmt == "" && format != "") && (fmt = format) - _fmt = (fmt == "" && format == "") ? ".eps" : string(fmt) + _fmt::String = (fmt == "" && format == "") ? ".eps" : string(fmt) _fmt[1] != '.' && (_fmt = "." * _fmt) return _fmt, name end diff --git a/src/libgmt.jl b/src/libgmt.jl index 6b8ba03fa..7d33338e3 100644 --- a/src/libgmt.jl +++ b/src/libgmt.jl @@ -414,10 +414,22 @@ function GMT_Put_Strings(API::Ptr{Cvoid}, family::Integer, object::Ptr{Cvoid}, t ccall((:GMT_Put_Strings, libgmt), Cint, (Cstring, UInt32, Cstring, Ptr{Ptr{UInt8}}), API, family, object, txt) end -function gmt_get_rgb_from_z(API::Ptr{Cvoid}, P::Ptr{GMT_PALETTE}, value::Cdouble, rgb::Vector{Float64}) +function gmt_get_rgb_from_z(API::Ptr{Cvoid}, P::Ptr{GMT_PALETTE}, value, rgb) ccall((:gmt_get_rgb_from_z, libgmt), Cint, (Cstring, Ptr{Cvoid}, Cdouble, Ptr{Cdouble}), GMT_Get_Ctrl(API), P, value, rgb) end +function gmt_get_index(API::Ptr{Cvoid}, P::Ptr{GMT_PALETTE}, value) + ccall((:gmt_get_index, libgmt), Cint, (Cstring, Ptr{Cvoid}, Ptr{Cdouble}), GMT_Get_Ctrl(API), P, value) +end + +function gmt_get_rgb_lookup(API::Ptr{Cvoid}, P::Ptr{GMT_PALETTE}, ind, value, rgb) + ccall((:gmt_get_rgb_lookup, libgmt), Cint, (Cstring, Ptr{Cvoid}, Cint, Cdouble, Ptr{Cdouble}), GMT_Get_Ctrl(API), P, ind, value, rgb) +end + +function gmt_illuminate(API::Ptr{Cvoid}, intensity, rgb::Vector{Float64}) + ccall((:gmt_illuminate, libgmt), Cvoid, (Cstring, Cdouble, Ptr{Cdouble}), GMT_Get_Ctrl(API), intensity, rgb) +end + function gmt_getrgb(API::Ptr{Cvoid}, colorname::String, rgb::Vector{Float64}) ccall((:gmt_getrgb, libgmt), Cuchar, (Cstring, Ptr{UInt8}, Ptr{Cdouble}), GMT_Get_Ctrl(API), colorname, rgb) end diff --git a/src/psxy.jl b/src/psxy.jl index 8817976d7..54157bd68 100644 --- a/src/psxy.jl +++ b/src/psxy.jl @@ -56,7 +56,7 @@ function common_plot_xyz(cmd0::String, arg1, caller::String, first::Bool, is3D:: (opt_p == "") ? (opt_p = CURRENT_VIEW[1]; cmd *= opt_p) : (CURRENT_VIEW[1] = opt_p) # Save for eventual use in other modules. if (is3D && isFV(arg1)) # case of 3D faces - arg1 = deal_faceverts(arg1, d, opt_p) + arg1 = deal_faceverts(arg1, d) !haskey(d, :aspect3) && (d[:aspect3] = "equal") # Needs thinking elseif is_gridtri arg1 = sort_visible_triangles(arg1) @@ -126,8 +126,6 @@ function common_plot_xyz(cmd0::String, arg1, caller::String, first::Bool, is3D:: end if (isa(arg1, GDtype) && !contains(opt_f, "T") && !contains(opt_f, "t") && !contains(opt_R, "T") && !contains(opt_R, "t")) - #isa(arg1, GMTdataset) && (cmd = set_fT(arg1, cmd, opt_f)) - #isa(arg1, Vector{<:GMTdataset}) && (cmd = set_fT(arg1[1], cmd, opt_f)) cmd = isa(arg1, GMTdataset) ? set_fT(arg1, cmd, opt_f) : set_fT(arg1[1], cmd, opt_f) end @@ -367,12 +365,34 @@ function plt_txt_attrib!(D::GDtype, d::Dict, _cmd::Vector{String}) end # --------------------------------------------------------------------------------------------------- -function deal_faceverts(arg1, d, opt_p) +# replicate=(centers=[], zcolor=[], cmap=..., scales=..., ) +function deal_replicants(arg1, d) + (val = find_in_dict(d, [:replicate])[1]) === nothing && return arg1 + ((xyz = get(val, :centers, nothing)) === nothing) && error("Can't replicate without centers") + _xyz = (promote_type(eltype(arg1[1].data), eltype(xyz))).(xyz) + ((zcolor = get(val, :zcolor, nothing)) === nothing) && (zcolor = collect(1.0:size(xyz, 1))) + cpt::GMTcpt = get(val, :cmap, GMTcpt()) + (isempty(cpt)) && (cpt = gmt("makecpt -T1/$(length(zcolor))")) + scales = get(val, :scales, nothing) + spl = split(CURRENT_VIEW[1][4:end], '/') + azim = parse(Float64, spl[1]) + elev = length(spl) > 1 ? parse(Float64, spl[2]) : 90.0 + arg1, = sort_visible_faces(arg1, azim, elev) # Sort & kill invisible + replicant(arg1, _xyz, azim, elev, zcolor, cpt; scales=scales) +end + +# --------------------------------------------------------------------------------------------------- +""" + deal_faceverts(arg1, d) +""" +function deal_faceverts(arg1, d) # Deal with the situation where we are plotting 3D FV's - spl = split(opt_p[4:end], '/') + (is_in_dict(d, [:replicate]) !== nothing) && return deal_replicants(arg1, d) + + spl = split(CURRENT_VIEW[1][4:end], '/') azim = parse(Float64, spl[1]) elev = length(spl) > 1 ? parse(Float64, spl[2]) : 90.0 - arg1, dotprod = sort_visible_faces(arg1, azim, elev) + arg1, dotprod = sort_visible_faces(arg1, azim, elev) # Sort & kill invisible if (is_in_dict(d, [:G :fill]) === nothing) # If fill not set we use the dotprod and a gray CPT to set the fill is_in_dict(d, [:Z :level :levels]) === nothing && (d[:Z] = dotprod) (is_in_dict(d, CPTaliases) === nothing) && (d[:C] = gmt("makecpt -T0/1 -C150,210")) # Users may still set a CPT @@ -1424,7 +1444,7 @@ end # --------------------------------------------------------------------------------------------------- """ - FV = sort_visible_faces(FV::Vector{<:GMTdataset}, azim, elev) -> ::Vector{GMTdataset} + FV, projs = sort_visible_faces(FV::Vector{<:GMTdataset}, azim, elev) -> Tuple{Vector{GMTdataset}, Vector{Float64}} Take a Faces-Vertices dataset and delete the invisible faces from view vector. Next sort them by distance so that the furthest faces are drawn on first and hence do not hide the others. @@ -1462,6 +1482,17 @@ function sort_visible_faces(FV::Vector{<:GMTdataset}, azim, elev)::Tuple{Vector{ end # --------------------------------------------------------------------------------------------------- +""" + Dv = sort_visible_triangles(Dv::Vector{<:GMTdataset}; del_hidden=false, zfact=1.0) -> Vector{GMTdataset} + +Take a Faces-Vertices dataset produced by grid2tri() and sort them by distance so that the furthest faces are +drawn on first and hence do not hide others. Optionally remove the invisible triangles. This not true by default +because we may want to see the inside of a top surface. + +`zfact` is a factor use to bring the z units to the same of the x,y units as for example when `z` is in km +and `x` and `y` are in degrees, case in which an automatic projection takes place, or in meters. We need this +if we want that the normal compuations makes sense. +""" function sort_visible_triangles(Dv::Vector{<:GMTdataset}; del_hidden=false, zfact=1.0)::Vector{GMTdataset} azim, elev = parse.(Float64, split(CURRENT_VIEW[1][4:end], '/')) sin_az, cos_az, sin_el = sind(azim), cosd(azim), sind(elev) @@ -1477,7 +1508,7 @@ function sort_visible_triangles(Dv::Vector{<:GMTdataset}; del_hidden=false, zfac end # ---------------------- Now sort by distance to the viewer ---------------------- - Dc = gmtspatial(Dv, Q=true, o="0,1") + Dc = gmtspatial(Dv, Q=true, o="0,1") # Should this directly in Julia and avoid calling GMT (=> a copy of Dv) dists = [(Dc.data[1,1] * sin_az + Dc.data[1,2] * cos_az, (Dv[1].bbox[5] + Dv[1].bbox[6]) / 2 * sin_el)] for k = 2:size(Dc, 1) push!(dists, (Dc.data[k,1] * sin_az + Dc.data[k,2] * cos_az, (Dv[k].bbox[5] + Dv[k].bbox[6]) / 2 * sin_el)) @@ -1491,7 +1522,7 @@ function sort_visible_triangles(Dv::Vector{<:GMTdataset}; del_hidden=false, zfac return Dv end -# --------------------------------------------------------------------------------------------------- +#= --------------------------------------------------------------------------------------------------- function tri_normals(Dv::Vector{<:GMTdataset}; zfact=1.0) t = isgeog(Dv) ? mapproject(Dv, J="t$((Dv[1].ds_bbox[1] + Dv[1].ds_bbox[2])/2)/1:1", C=true, F=true) : Dv Dc = gmtspatial(Dv, Q=true) @@ -1515,3 +1546,50 @@ function quiver3(Dv::Vector{<:GMTdataset}; first=true, zfact=1.0, kwargs...) D = mat2ds(mat, geom=wkbLineStringZ) common_plot_xyz("", D, "plot3d", first, true, kwargs...) end +=# + +# --------------------------------------------------------------------------------------------------- +function replicant(FV::Vector{<:GMTdataset}, xyz::Matrix{<:Real}, azim, elev, cval=1:size(xyz, 1), + cpt::GMTcpt=GMTcpt(); scales=nothing) + + @assert length(cval) == size(xyz, 1) + if (scales !== nothing) + (length(scales) == 1) && (scales = fill(scales, length(cval))) + @assert length(scales) == length(cval) + end + FV, normals = sort_visible_faces(FV::Vector{<:GMTdataset}, azim, elev) + V::GMTdataset{Float64,2}, F::GMTdataset{Int,2} = FV[1], FV[2] + + n_rows = size(F.data, 2) # Number of rows (vertexes of the polygon) + n_cols = size(V.data, 2) # Number of columns (2 for x,y; 3 for x,y,z) + n_faces = size(F.data, 1) # Number of faces of base body + + D1 = Vector{GMTdataset}(undef, n_faces) + t = zeros(eltype(V.data), n_rows, n_cols) + for face = 1:n_faces # Loop over faces + for c = 1:n_cols, r = 1:n_rows + t[r,c] = V.data[F.data[face, r], c] + end + D1[face] = GMTdataset(data=copy(t)) + end + + cpt::GMTcpt = (isempty(cpt)) ? gmt(T="1/$(length(cval))") : cpt + P::Ptr{GMT_PALETTE} = palette_init(G_API[1], cpt) # A pointer to a GMT CPT + rgb = [0.0, 0.0, 0.0, 0.0] + cor = [0.0, 0.0, 0.0] + D2 = Vector{GMTdataset}(undef, size(xyz, 1) * n_faces) + + _scales = (scales === nothing) ? fill(one(eltype(V.data)), length(cval)) : # Try to make _scales the same type as the V points + (eltype(scales) == Float64 && eltype(V.data) == Float32) ? Float32.(scales) : scales + + for k = 1:size(xyz, 1) # Loop over number of positions. For each of these we have a new body + gmt_get_rgb_from_z(G_API[1], P, cval[k], rgb) + for face = 1:n_faces # Loop over number of faces of the base body + cor[1], cor[2], cor[3] = rgb[1], rgb[2], rgb[3] + gmt_illuminate(G_API[1], normals[face], cor) + txt = @sprintf("-G%.0f/%.0f/%.0f", cor[1]*255, cor[2]*255, cor[3]*255) + D2[(k-1)*n_faces+face] = GMTdataset(data=(D1[face].data .+ xyz[k:k,1:3]) * _scales[k], header=txt) + end + end + return set_dsBB(D2) +end diff --git a/src/solids.jl b/src/solids.jl index 4c3fc74f9..44324b6cf 100644 --- a/src/solids.jl +++ b/src/solids.jl @@ -1,28 +1,34 @@ # This file contains modified versions of some functions from the Comodo.jl package (Apache 2.0 license). # https://github.com/COMODO-research/Comodo.jl """ - FV = icosahedron(r=1.0) + FV = icosahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) -Creates an icosahedron mesh with radius `r. +Creates an icosahedron mesh with radius `r`. + +- `radius`: the keyword `radius` is an alternative to the positional argument `r`. +- `origin`: A tuple of three numbers defining the origin of the body. Default is `(0.0, 0.0, 0.0)`. """ -function icosahedron(r=1.0) +function icosahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) + + _r::Float64 = (radius != 1.0) ? radius : r # If spelled, the `radius` kwarg take precedence + o::Matrix{Float64} = [origin[1] origin[2] origin[3]] ϕ = Base.MathConstants.golden # (1.0+sqrt(5.0))/2.0, Golden ratio - s = r / sqrt(ϕ + 2.0) + s = _r / sqrt(ϕ + 2.0) t = ϕ * s - V = [0.0 -s -t # The Vertices - 0.0 -s t - 0.0 s t - 0.0 s -t - -s -t 0.0 - -s t 0.0 - s t 0.0 - s -t 0.0 - -t 0.0 -s - t 0.0 -s - t 0.0 s - -t 0.0 s] + V =[0 -s -t # The Vertices + 0 -s t + 0 s t + 0 s -t + -s -t 0 + -s t 0 + s t 0 + s -t 0 + -t 0 -s + t 0 -s + t 0 s + -t 0 s] F = [9 4 1 # The Faces 1 5 9 @@ -45,26 +51,32 @@ function icosahedron(r=1.0) 6 7 4 6 4 9] + (o[1] != 0 || o[2] != 0 || o[3] != 0) && (V .+= o) DV = GMTdataset(data=V, geom=wkbPointZ); set_dsBB!(DV) return [DV, GMTdataset(data=F)] end # -------------------------------------------------------- """ - FV = octahedron(r=1.0) + FV = octahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) Creates an octahedron mesh with radius `r`. + +- `radius`: the keyword `radius` is an alternative to the positional argument `r`. +- `origin`: A tuple of three numbers defining the origin of the body. Default is `(0.0, 0.0, 0.0)`. """ -function octahedron(r=1.0) +function octahedron(r=1.0; radius=0.0, origin=(0.0, 0.0, 0.0)) - s = r/sqrt(2.0) + _r::Float64 = (radius != 1.0) ? radius : r # If spelled, the `radius` kwarg take precedence + o::Matrix{Float64} = [origin[1] origin[2] origin[3]] + s = _r / sqrt(2.0) V = [-s -s 0.0 # The Vertices s -s 0.0 s s 0.0 -s s 0.0 - 0.0 0.0 -r - 0.0 0.0 r] + 0.0 0.0 -_r + 0.0 0.0 _r] F = [5 2 1 # The Faces 5 3 2 @@ -75,20 +87,27 @@ function octahedron(r=1.0) 6 3 4 6 4 1] + (o[1] != 0 || o[2] != 0 || o[3] != 0) && (V .+= o) DV = GMTdataset(data=V, geom=wkbPointZ); set_dsBB!(DV) return [DV, GMTdataset(data=F)] end # ---------------------------------------------------------------------------- """ - FV = dodecahedron(r=1.0) + FV = dodecahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) Creates an dodecahedron mesh with radius `r`. + +- `radius`: the keyword `radius` is an alternative to the positional argument `r`. +- `origin`: A tuple of three numbers defining the origin of the body. Default is `(0.0, 0.0, 0.0)`. """ -function dodecahedron(r=1.0) +function dodecahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) + + _r::Float64 = (radius != 1.0) ? radius : r # If spelled, the `radius` kwarg take precedence + o::Matrix{Float64} = [origin[1] origin[2] origin[3]] ϕ = Base.MathConstants.golden # (1.0+sqrt(5.0))/2.0, Golden ratio - s = r/sqrt(3.0) + s = _r / sqrt(3.0) t = ϕ*s w = (ϕ-1.0)*s @@ -125,26 +144,33 @@ function dodecahedron(r=1.0) 1 2 8 17 4 1 4 19 15 20 12 18 5 15 19] - + + (o[1] != 0 || o[2] != 0 || o[3] != 0) && (V .+= o) DV = GMTdataset(data=V, geom=wkbPointZ); set_dsBB!(DV) return [DV, GMTdataset(data=F)] end # ---------------------------------------------------------------------------- """ - FV = tetrahedron(r=1.0) + FV = tetrahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) Creates a tetrahedron mesh with radius `r`. + +- `radius`: the keyword `radius` is an alternative to the positional argument `r`. +- `origin`: A tuple of three numbers defining the origin of the body. Default is `(0.0, 0.0, 0.0)`. """ -function tetrahedron(r=1.0) +function tetrahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) + + _r::Float64 = (radius != 1.0) ? radius : r # If spelled, the `radius` kwarg take precedence + o::Matrix{Float64} = [origin[1] origin[2] origin[3]] a = r*sqrt(2.0)/sqrt(3.0) - b = -r*sqrt(2.0)/3.0 - c = -r/3.0 + b = -_r * sqrt(2.0)/3.0 + c = -_r / 3.0 V = [-a b c # The Vertices a b c - 0.0 0.0 r + 0.0 0.0 _r 0.0 -2.0*b c] F = [1 2 3 # The Faces @@ -152,19 +178,26 @@ function tetrahedron(r=1.0) 4 3 2 4 1 3] + (o[1] != 0 || o[2] != 0 || o[3] != 0) && (V .+= o) DV = GMTdataset(data=V, geom=wkbPointZ); set_dsBB!(DV) return [DV, GMTdataset(data=F)] end # ---------------------------------------------------------------------------- """ - FV = cube(r=1.0) + FV = cube(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) Creates a cube mesh with radius `r`. + +- `radius`: the keyword `radius` is an alternative to the positional argument `r`. +- `origin`: A tuple of three numbers defining the origin of the body. Default is `(0.0, 0.0, 0.0)`. """ -function cube(r=1.0) +function cube(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0)) - s = r/sqrt(3.0) + _r::Float64 = (radius != 1.0) ? radius : r # If spelled, the `radius` kwarg take precedence + o::Matrix{Float64} = [origin[1] origin[2] origin[3]] + + s = _r / sqrt(3.0) V = [-s -s -s # The Vertices -s s -s @@ -182,13 +215,14 @@ function cube(r=1.0) 7 8 4 3 8 5 1 4] + (o[1] != 0 || o[2] != 0 || o[3] != 0) && (V .+= o) DV = GMTdataset(data=V, geom=wkbPointZ); set_dsBB!(DV) return [DV, GMTdataset(data=F)] end # ---------------------------------------------------------------------------- """ - FV = geosphere(n, r=1; radius=1.0) + FV = sphere(r=1; n=1, radius=1.0, center=(0.0, 0.0, 0.0)) Generate a geodesic sphere triangulation based on the number of refinement iterations `n` and the radius `r`. Geodesic spheres (aka Buckminster-Fuller spheres) are triangulations @@ -197,25 +231,29 @@ icosahedron. Next this icosahedron is refined `n` times, while nodes are pushed surface with radius `r` at each iteration. - `radius`: the keyword `radius` is an alternative to the positional argument `r`. -- `n`: is the number of times of iterations used to obtain the sphere from the icosahedron. +- `n`: is the number of iterations used to obtain the sphere from the icosahedron. +- `center`: A tuple of three numbers defining the center of the sphere. ### Returns A two elements vector of GMTdataset where first contains the vertices and the second the indices that define the faces. """ -function geosphere(n, r=1; radius=1.0) - (radius != 1 && r == 1) && (r = radius) - FV = icosahedron(r) +function sphere(r=1; n=1, radius=1.0, center=(0.0, 0.0, 0.0)) + _r::Float64 = (radius != 1.0) ? radius : r # If spelled, the `radius` kwarg take precedence + o::Matrix{Float64} = [center[1] center[2] center[3]] + + FV = icosahedron(_r) # Can't use the center here because subTriSplit() screws it up # Sub-triangulate the icosahedron for geodesic sphere triangulation if (n > 0) # If refinement is requested Fn::Matrix{Int} = FV[2].data; Vn::Matrix{Float64} = FV[1].data; # Initialise Fn, Vn as input F and V for q = 1:n # iteratively refine triangulation and push back radii to be r Vn, Fn = subTriSplit(Vn, Fn, 1) # Sub-triangulate T, P, R = cart2sph(view(Vn,:,1),view(Vn,:,2),view(Vn,:,3)) # Convert to spherical coordinates - x, y, z = sph2cart(T, P, r.*ones(size(R))) # Push back radii + x, y, z = sph2cart(T, P, _r .* ones(size(R))) # Push back radii Vn[:,1] .= x; Vn[:,2] .= y; Vn[:,3] .= z end end + (o[1] != 0 || o[2] != 0 || o[3] != 0) && (Vn .+= o) DV = GMTdataset(data=Vn, geom=wkbPointZ); set_dsBB!(DV) return [DV, GMTdataset(data=Fn)] end diff --git a/src/utils_types.jl b/src/utils_types.jl index c4eb3b0c4..50284c8de 100644 --- a/src/utils_types.jl +++ b/src/utils_types.jl @@ -930,11 +930,11 @@ function line2multiseg(M::Matrix{<:Real}; is3D::Bool=false, color::GMTcpt=GMTcpt if (!isempty(color)) z_col = color_col - rgb = [0.0, 0.0, 0.0]; + rgb = [0.0, 0.0, 0.0, 0.0] P::Ptr{GMT_PALETTE} = palette_init(G_API[1], color); # A pointer to a GMT CPT for k = 1:n_ds z = (use_row_number) ? z4color[k] : M[k, z_col] - @GC.preserve color gmt_get_rgb_from_z(G_API[1], P, z, rgb) + gmt_get_rgb_from_z(G_API[1], P, z, rgb) t = @sprintf(",%.0f/%.0f/%.0f", rgb[1]*255, rgb[2]*255, rgb[3]*255) _hdr[k] = (first) ? " -W"*t : _hdr[k] * t end @@ -2219,9 +2219,9 @@ end mksymbol(f::Function, cmd0::String="", arg1=nothing; kwargs...) """ function mksymbol(f::Function, cmd0::String="", arg1=nothing; kwargs...) - # Make a fig and convert it to EPS so it can be used as a custom symbol is plot(3) + # Make a fig and convert it to EPS so it can be used as a custom symbol in plot(3) d = KW(kwargs) - t = ((val = find_in_dict(d, [:symbname :symb_name :symbol])[1]) !== nothing) ? string(val) : "GMTsymbol" + t::String = ((val = find_in_dict(d, [:symbname :symb_name :symbol])[1]) !== nothing) ? string(val) : "GMTsymbol" (t == "GMTsymbol" && (f == flower_minho || f == matchbox)) && (t = string(f)) !haskey(d, :name) && (d[:name] = t * ".eps") if (f == flower_minho || f == matchbox) # Special case for the Flower Minho symbol diff --git a/test/test_solidos.jl b/test/test_solidos.jl index a20a31d12..f2e60a70b 100644 --- a/test/test_solidos.jl +++ b/test/test_solidos.jl @@ -5,5 +5,5 @@ GMT.dodecahedron() GMT.tetrahedron() GMT.cube() - GMT.geosphere(2) + GMT.sphere() end