Merge pull request #1507 from GenericMappingTools/FV-in-plot

Deal with the FaceVertices data (hide invisible faces) in plot3.

Author: joa-quim

src/common_options.jl

@@ -386,7 +386,18 @@ function parse_JZ(d::Dict, cmd::String, del::Bool=true; O::Bool=false, is3D::Boo
 		(o == "") && (o = " ")		# When scaning a string with only -J (i.e., no args), do not error
 		(o[1] != 'X' || o[end] == 'd') &&  @warn("aspect3 works only in linear projections (and no geog), ignoring it.") 
 		if (o[1] == 'X' && o[end] != 'd' && length(o) > 1)
-			opt_J = " -JZ" * split(o[2:end],'/')[1];		seek_JZ = false
+			s_val::String = string(val)
+			if (contains(s_val, ':'))						# An explicit aspect ratio, e.g. 3:2
+				dims = parse.(Float64, split(s_val, ':'))
+				t = split(o[2:end],'/')[1]
+				isletter(t[end]) && (t = t[1:end-1])		# Remove the letter 'x' which hopefully a 'c'
+				w = parse(Float64, t)
+				h = w * dims[2] / dims[1]					# Apply the aspect ratio
+				opt_J = string(" -JZ", h)
+			else
+				opt_J = " -JZ" * split(o[2:end],'/')[1];
+			end
+			seek_JZ = false
 			cmd *= opt_J
 		end
 	end
@@ -602,9 +613,12 @@ function append_figsize(d::Dict, opt_J::String, width::String="", scale::Bool=fa
 	return opt_J
 end
 
-set_aspect_ratio(aspect::Nothing, width::String, def_fig::Bool=false, is_aspect3::Bool=false)::String = set_aspect_ratio("", width, def_fig, is_aspect3)
-set_aspect_ratio(aspect::Symbol, width::String, def_fig::Bool=false, is_aspect3::Bool=false)::String = set_aspect_ratio(string(aspect), width, def_fig, is_aspect3)
-set_aspect_ratio(aspect::Real, width::String, def_fig::Bool=false, is_aspect3::Bool=false)::String = set_aspect_ratio(string(aspect, ":1"), width, def_fig, is_aspect3)
+set_aspect_ratio(aspect::Nothing, width::String, def_fig::Bool=false, is_aspect3::Bool=false)::String =
+	set_aspect_ratio("", width, def_fig, is_aspect3)
+set_aspect_ratio(aspect::Symbol, width::String, def_fig::Bool=false, is_aspect3::Bool=false)::String =
+	set_aspect_ratio(string(aspect), width, def_fig, is_aspect3)
+set_aspect_ratio(aspect::Real, width::String, def_fig::Bool=false, is_aspect3::Bool=false)::String =
+	set_aspect_ratio(string(aspect, ":1"), width, def_fig, is_aspect3)
 function set_aspect_ratio(aspect::String, width::String, def_fig::Bool=false, is_aspect3::Bool=false)::String
 	# Set the aspect ratio. ASPECT can be "equal", "eq"; "square", "sq" or a ratio in the form "4:3", "16:12", etc.
 	def_fig && (width = split(DEF_FIG_SIZE, '/')[1])
@@ -3711,6 +3725,7 @@ function _read_data(d::Dict, cmd::String, arg, opt_R::String="", is3D::Bool=fals
 			if (!occursin("?", opt_R) && !is_onecol)		# is_onecol is true for DateTime data
 				dx::Float64 = (wesn_f64[2] - wesn_f64[1]) * 0.005;	dy::Float64 = (wesn_f64[4] - wesn_f64[3]) * 0.005;
 				wesn_f64[1] -= dx;	wesn_f64[2] += dx;	wesn_f64[3] -= dy;	wesn_f64[4] += dy;
+				length(wesn_f64) == 6 && (dz::Float64 = (wesn_f64[6] - wesn_f64[5]) * 0.005; wesn_f64[5] -= dz;	wesn_f64[6] += dz;)
 				wesn_f64 = round_wesn(wesn_f64, is_geo)		# Add a pad if not-tight
 				if (isGMTdataset(arg))						# Needed for the guess_proj case
 					if ((wesn_f64[3] < -90 || wesn_f64[4] > 90) || ((wesn_f64[2] - wesn_f64[1]) > 360))
@@ -3779,8 +3794,10 @@ function round_wesn(_wesn::Vector{Float64}, geo::Bool=false, pad=zeros(2))::Vect
 		dy = (wesn[4] - wesn[3]) * pad[2]
 		wesn[1:4] = [wesn[1]-dx, wesn[2]+dx, wesn[3]-dy, wesn[4]+dy]
 	end
-	set::Vector{Bool} = zeros(Bool, 2)
-	range::Vector{Float64} = [0.0, 0.0]
+	set = zeros(Bool, 3)
+	range = [0.0, 0.0, 0.0]
+	n_axe = (length(wesn) == 6) ? 3 : 2
+
 	if (wesn[1] == wesn[2])
 		wesn[1] -= abs(wesn[1]) * 0.05;	wesn[2] += abs(wesn[2]) * 0.05
 		if (wesn[1] == wesn[2])  wesn[1] = -0.1;	wesn[2] = 0.1;	end		# x was = 0
@@ -3791,6 +3808,7 @@ function round_wesn(_wesn::Vector{Float64}, geo::Bool=false, pad=zeros(2))::Vect
 	end
 	range[1] = wesn[2] - wesn[1]
 	range[2] = wesn[4] - wesn[3]
+	(n_axe == 3) && (range[3] = wesn[6] - wesn[5])
 	if (geo) 					# Special checks due to periodicity
 		if (range[1] > 306.0) 	# If within 15% of a full 360 we promote to 360
 			if ((wesn[1] + wesn[2]) / 2 < 100)  wesn[1] = -180.;	wesn[2] = 180.
@@ -3806,14 +3824,14 @@ function round_wesn(_wesn::Vector{Float64}, geo::Bool=false, pad=zeros(2))::Vect
 
 	_log10(x) = log(x) / 2.30258509299	# Compute log10 with ln because JET & SnoopCompile acuse it of "failed to optimize"
 	item = 1
-	for side = 1:2
+	for side = 1:n_axe
 		set[side] && continue			# Done above */
 		mag::Float64 = round(_log10(range[side])) - 1.0
 		inc = exp10(mag)
 		if ((range[side] / inc) > 10.0) inc *= 2.0	end	# Factor of 2 in the rounding
 		if ((range[side] / inc) > 10.0) inc *= 2.5	end	# Factor of 5 in the rounding
 		s = 1.0
-		if (geo) 	# Use arc integer minutes or seconds if possible
+		if (geo && side < 3) 	# Use arc integer minutes or seconds if possible
 			if (inc < 1.0 && inc > 0.05) 				# Nearest arc minute
 				s, inc = 60.0, 1.0
 				if ((s * range[side] / inc) > 10.0) inc *= 2.0	end		# 2 arcmin
@@ -3897,12 +3915,13 @@ end
 
 # ---------------------------------------------------------------------------------------------------
 """
-    isFV(D::Vector{<:GMTdataset})::Bool
+    isFV(D)::Bool
 
-Check if D, a vector of GMTdataset, is a Face-Vertices ensemble.
+Check if D is a Face-Vertices ensemble (a 2 elements vector of GMTdataset).
 """
-function isFV(D::Vector{<:GMTdataset})::Bool
-	length(D) == 2 && (D[1].geom == wkbPoint || D[1].geom == wkbPointZ) && eltype(D[2].data) <: Integer
+function isFV(D)::Bool
+	(isa(D, Vector{<:GMTdataset}) && length(D) > 1) &&
+		(D[1].geom == wkbPoint || D[1].geom == wkbPointZ) && eltype(D[2].data) <: Integer
 end
 
 # ---------------------------------------------------------------------------------------------------

src/gmt_main.jl

@@ -1094,25 +1094,18 @@ function dataset_init_FV(API::Ptr{Nothing}, FV)::Ptr{GMT_MATRIX}
 	n_segs = size(F.data, 1)				# Number of segments or faces (polygons)
 	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)
-	view_vec = [sind(200) * cosd(30), cosd(200) * cosd(30), sind(30)]
-
-	view_proj = triage_faces(V, F, view_vec)
-	n_visible_faces = sum(view_proj .> 0)
-	dim = [1, n_visible_faces, n_rows, n_cols]		# [1, GMT_SEG+1, GMT_ROW+1, GMT_COL+1]
+	dim = [1, n_segs, n_rows, n_cols]		# [1, GMT_SEG+1, GMT_ROW+1, GMT_COL+1]
 
 	pdim = pointer(dim)
 	D = convert(Ptr{GMT_DATASET}, GMT_Create_Data(API, GMT_IS_DATASET, GMT_IS_PLP, GMT_NO_STRINGS, pdim, NULL, NULL, 0, 0, NULL))
 	DS::GMT_DATASET = unsafe_load(D)
 	DT = unsafe_load(unsafe_load(DS.table))			# GMT_DATATABLE
 
-	n_records, count_vis = 0, 0
+	n_records = 0
 	tmp = zeros(n_rows, n_cols)
 
 	for seg = 1:n_segs 								# Each row in F (a face) is a new data segment (a polygon)
-		view_proj[seg] <= 0 && continue
-		count_vis += 1
-
-		DSv = convert(Ptr{Nothing}, unsafe_load(DT.segment, count_vis))		# DT.segment = Ptr{Ptr{GMT_DATASEGMENT}}
+		DSv = convert(Ptr{Nothing}, unsafe_load(DT.segment, seg))		# DT.segment = Ptr{Ptr{GMT_DATASEGMENT}}
 		S = GMT_Alloc_Segment(API, GMT_NO_STRINGS, n_rows, n_cols, "", DSv) # Ptr{GMT_DATASEGMENT}
 		Sb = unsafe_load(S)							# GMT_DATASEGMENT;		Sb.data -> Ptr{Ptr{Float64}}
 
@@ -1126,7 +1119,7 @@ function dataset_init_FV(API::Ptr{Nothing}, FV)::Ptr{GMT_MATRIX}
 		n_records += n_rows							# Must manually keep track of totals
 		DS.type_ = GMT_READ_DATA
 		unsafe_store!(S, Sb)
-		unsafe_store!(DT.segment, S, count_vis)
+		unsafe_store!(DT.segment, S, seg)
 	end
 	DT.n_records, DS.n_records = n_records, n_records	# They are equal because our GMT_DATASET have only one table
 	Dt = unsafe_load(DS.table)
@@ -1137,21 +1130,6 @@ function dataset_init_FV(API::Ptr{Nothing}, FV)::Ptr{GMT_MATRIX}
 	return D
 end
 
-function triage_faces(V, F, view_vec)
-	# Compute the dot product between the view vector and the normal of each face
-	n_faces = size(F.data, 1)		# Number of segments or faces (polygons)
-	n_rows = size(F.data, 2)		# Number of rows (vertexes of the polygon)
-	tmp = zeros(n_rows, 3)
-	proj = zeros(n_faces)
-	for face = 1:n_faces 			# Each row in F (a face) is a new data segment (a polygon)
-		for c = 1:3, r = 1:n_rows
-			tmp[r,c] = V.data[F.data[face,r], c]
-		end
-		proj[face] = dot(facenorm(tmp), view_vec)
-	end
-	return proj
-end
-
 # ---------------------------------------------------------------------------------------------------
 function dataset_init(API::Ptr{Nothing}, ptr, actual_family::Vector{<:Integer})::Ptr{GMT_MATRIX}
 # Create empty Matrix container, associate it with julia data matrix, and use as GMT input.

src/imshow.jl

@@ -70,7 +70,9 @@ function imshow(arg1, x::AbstractVector{Float64}=Float64[], y::AbstractVector{Fl
 		isa(arg1, Matrix{<:Real}) && (arg1 = mat2ds(arg1))
 		ginfo = isa(arg1, GMTdataset) ? arg1.bbox : arg1[1].ds_bbox
 		CTRL.limits[1:4] = ginfo[1:4];		CTRL.limits[7:10] = ginfo[1:4]
-		call_plot3 = ((isa(arg1, GMTdataset) && arg1.geom == Gdal.wkbLineStringZ) || (isa(arg1, Vector{<:GMTdataset}) && arg1[1].geom == Gdal.wkbLineStringZ)) ? true : false		# Should evolve into a fun that detects the several plot3d cases.
+		call_plot3 = ((isa(arg1, GMTdataset) && arg1.geom == Gdal.wkbLineStringZ) ||
+		              (isa(arg1, Vector{<:GMTdataset}) && arg1[1].geom == Gdal.wkbLineStringZ) ||
+					  isFV(arg1)) ? true : false		# Should evolve into a fun that detects the several plot3d cases.
 		!call_plot3 && (call_plot3 = isplot3(kw))
 		return (call_plot3) ? plot3d(arg1; show=see, kw...) : plot(arg1; show=see, kw...)
 	elseif (isa(arg1, GMTcpt))

src/psxy.jl

@@ -44,7 +44,13 @@ function common_plot_xyz(cmd0::String, arg1, caller::String, first::Bool, is3D::
 	end
 
 	if (occursin('3', caller) && !haskey(d, :p) && !haskey(d, :view) && !haskey(d, :perspective))
-		d[:p] = "200/30"		# Need this before parse_BJR() so MAP_FRAME_AXES can be guessed.
+		d[:p] = "200/30"			# Need this before parse_BJR() so MAP_FRAME_AXES can be guessed.
+	end
+	cmd, opt_p = parse_p(d, cmd)	# Parse this one (view angle) aside so we can use it to remove invisible faces (3D)
+
+	if (is3D && isFV(arg1))			# case of 3D faces
+		arg1 = deal_faceverts(arg1, d, opt_p)
+		!haskey(d, :aspect3) && (d[:aspect3] = "equal")			# Needs thinking
 	end
 
 	isa(arg1, GMTdataset) && (arg1 = with_xyvar(d, arg1))		# See if we have a column request based on column names
@@ -82,7 +88,7 @@ function common_plot_xyz(cmd0::String, arg1, caller::String, first::Bool, is3D::
 
 	cmd, opt_JZ = parse_JZ(d, cmd; O=O, is3D=is3D)
 	#(is3D && O && opt_JZ == "" && CTRL.pocket_J[3] != "") && (cmd *= CTRL.pocket_J[3])
-	cmd, = parse_common_opts(d, cmd, [:a :e :f :g :p :t :w :params], first)
+	cmd, = parse_common_opts(d, cmd, [:a :e :f :g :t :w :params], first)
 	cmd, opt_l = parse_l(d, cmd)		# Parse this one (legend) aside so we can use it in classic mode
 	cmd, opt_f = parse_f(d, cmd)		# Parse this one (-f) aside so we can check against D.attrib
 	cmd  = parse_these_opts(cmd, d, [[:D :shift :offset], [:I :intens], [:N :no_clip :noclip]])
@@ -333,6 +339,20 @@ function common_plot_xyz(cmd0::String, arg1, caller::String, first::Bool, is3D::
 	return r
 end
 
+# ---------------------------------------------------------------------------------------------------
+function deal_faceverts(arg1, d, opt_p)
+	# Deal with the situation where we are plotting 3D FV's
+	spl = split(opt_p[4:end], '/')
+	az = parse(Float64, spl[1])
+	elev = length(spl) > 1 ? parse(Float64, spl[2]) : 90.0
+	arg1, dotprod = visible_faces(arg1, [sind(az) * cosd(elev), cosd(az) * cosd(elev), sind(elev)])
+	if (is_in_dict(d, [:G :fill]) === nothing)		# If fill not set we use the dotprod and a gray CPT to set the fill
+		d[:Z] = dotprod
+		(is_in_dict(d, CPTaliases) === nothing) && (d[:C] = gmt("makecpt -T0/1 -C150,210"))	# Users may still set a CPT
+	end
+	return arg1
+end
+
 # ---------------------------------------------------------------------------------------------------
 function frame_opaque(cmd::Vector{String}, oB::String, oR::String, oJ::String, oJZ::String=""; bot::Bool=true)
 	# Transparency affects the frame too, which is bad. So, if we have a transparency request we
@@ -1340,3 +1360,30 @@ function zoom_reactangle(_cmd, isplot::Bool)
 	append!(_cmd, [ins])				# Add the inset call again
 	return _cmd
 end
+
+# ---------------------------------------------------------------------------------------------------
+function faces_normals_view(V::GMTdataset{Float64,2}, F::GMTdataset{Int,2}, view_vec)
+	# Compute the dot product between the view vector and the normal of each face
+	# Faces hose dot product is <= 0 are not visible
+	n_faces = size(F.data, 1)		# Number of segments or faces (polygons)
+	n_rows  = size(F.data, 2)		# Number of rows (vertices of the polygon)
+	tmp  = zeros(n_rows, 3)
+	proj = zeros(n_faces)
+	for face = 1:n_faces 			# Each row in F (a face) is a new data segment (a polygon)
+		for c = 1:3, r = 1:n_rows
+			tmp[r,c] = V.data[F.data[face,r], c]
+		end
+		proj[face] = dot(facenorm(tmp), view_vec)
+	end
+	return proj
+end
+
+# ---------------------------------------------------------------------------------------------------
+function visible_faces(FV::Vector{<:GMTdataset}, view_vec)
+	# Remove the faces that are not visible from the normal view_vec
+	V::GMTdataset{Float64,2}, F::GMTdataset{Int,2} = FV[1], FV[2]
+	proj = faces_normals_view(V, F, view_vec)
+	is_vis = (proj .> 0)	# Visible faces
+	F2 = mat2ds(F.data[is_vis, :])
+	return [V,F2], proj[is_vis]
+end

src/solids.jl

@@ -77,6 +77,56 @@ function octahedron(r=1.0)
 	return [mat2ds(V, geom=wkbPointZ), mat2ds(F)]
 end
 
+# ----------------------------------------------------------------------------
+"""
+    FV = dodecahedron(r=1.0)
+
+Creates an dodecahedron mesh with radius `r`. 
+"""
+function dodecahedron(r=1.0)
+
+	ϕ = Base.MathConstants.golden # (1.0+sqrt(5.0))/2.0, Golden ratio
+	s = r/sqrt(3.0)
+	t = ϕ*s    
+	w = (ϕ-1.0)*s
+
+	V = [ s   s   s		# The Vertices
+		  w 0.0   t
+		 -t  -w 0.0
+		  t   w 0.0
+		 -s   s  -s
+		 0.0 -t  -w
+		 -t   w 0.0
+		  s  -s   s
+		 -s   s   s
+		 -s  -s   s
+		  s  -s  -s
+		  w 0.0  -t
+		 -s  -s  -s
+		 0.0 -t   w
+		 0.0  t  -w
+		 -w 0.0   t
+		  t  -w 0.0
+		 -w 0.0  -t
+		  s   s  -s
+		 0.0  t   w]
+
+	F = [20  9 16  2  1			# The Faces
+		  2 16 10 14  8
+		 16  9  7  3 10
+		  7  9 20 15  5
+		 18 13  3  7  5
+		  3 13  6 14 10
+		  6 13 18 12 11
+		  6 11 17  8 14
+		 11 12 19  4 17
+		  1  2  8 17  4
+		  1  4 19 15 20
+		 12 18  5 15 19]
+    
+	return [mat2ds(V, geom=wkbPointZ), mat2ds(F)]
+end
+
 # ----------------------------------------------------------------------------
 """
     FV = tetrahedron(r=1.0)

test/test_solidos.jl

@@ -2,6 +2,7 @@
 
 	GMT.icosahedron()
 	GMT.octahedron()
+	GMT.dodecahedron()
 	GMT.tetrahedron()
 	GMT.cube()
 	GMT.geosphere(2)