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

@@ -184,6 +184,8 @@ export
 
 	lazinfo, lazread, lazwrite, lasread, laswrite,
 
+	cube, dodecahedron, icosahedron, sphere, octahedron, tetrahedron,
+
 	df2ds, ODE2ds,
 	sprintf,
 	@?, @dir

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

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

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

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

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