Re-writings to allow FV's with multiple geometries.

src/gmt_main.jl

@@ -1090,35 +1090,43 @@ function dataset_init(API::Ptr{Nothing}, Darr::Vector{<:GMTdataset}, direction::
 end
 
 # ---------------------------------------------------------------------------------------------------
-function dataset_init_FV(API::Ptr{Nothing}, FV)::Ptr{GMT_MATRIX}
-	n_segs = size(FV.faces[1], 1)			# Number of segments or faces (polygons)
-	n_rows = size(FV.faces[1], 2)			# Number of rows (vertices of the polygon)
-	dim = [1, n_segs, n_rows, 3]			# [1, GMT_SEG+1, GMT_ROW+1, GMT_COL+1]
+function dataset_init_FV(API::Ptr{Nothing}, FV::GMTfv)::Ptr{GMT_MATRIX}
+	# FV is a GMTfv type and these types may be complicated when they hold bodies with more than one geometry
+	# (e.g. triangles and quadrangles) and/or more than one body. This function deals with those complixities.
+	what_faces = isempty(FV.faces_view) ? FV.faces : FV.faces_view	# faces_view is a processed version for the current view
+	n_segs_tot = sum(size.(what_faces, 1))			# Total number of segments or faces (polygons)
+	n_rows_tot = sum(size.(what_faces, 2))			# Total number of rows (vertices of the polygon)
+	dim = [1, n_segs_tot, n_rows_tot, 3]			# [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
+	DT = unsafe_load(unsafe_load(DS.table))			# GMT_DATATABLE
 
-	n_records = 0
-	tmp = zeros(n_rows, 3)
+	n_records, _seg = 0, 0
+	tmp = zeros(maximum(size.(what_faces,2)), 3)	# The maximum number of vertices in any polygon of all geometries
+
+	for geom = 1:numel(what_faces)					# If we have more than one type of geometries (e.g. triangles and quadrangles)
+		n_segs = size(what_faces[geom], 1)			# Number of segments or faces (polygons) in this geometry (what_faces[geom])
+		n_rows = size(what_faces[geom], 2)			# Number of rows (vertices of the polygon) in this geometry
+		for seg = 1:n_segs 							# Each row in a face is a new data segment (a polygon)
+			_seg += 1								# Counter that keeps track on the current number of polygon. 
+			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, 3, "", DSv)	# Ptr{GMT_DATASEGMENT}
+			Sb = unsafe_load(S)						# GMT_DATASEGMENT;		Sb.data -> Ptr{Ptr{Float64}}
+
+			for c = 1:3, r = 1:n_rows
+				tmp[r,c] = FV.verts[what_faces[geom][seg, r], c]
+			end
+			for col = 1:3							# Copy the data columns
+				unsafe_copyto!(unsafe_load(Sb.data, col), pointer(tmp[:,col]), n_rows)
+			end
 
-	for seg = 1:n_segs 							# Each row in a face is a new data segment (a polygon)
-		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, 3, "", DSv) # Ptr{GMT_DATASEGMENT}
-		Sb = unsafe_load(S)						# GMT_DATASEGMENT;		Sb.data -> Ptr{Ptr{Float64}}
-
-		for c = 1:3, r = 1:n_rows
-			tmp[r,c] = FV.verts[FV.faces[1][seg, r], c]
-		end
-		for col = 1:3							# Copy the data columns
-			unsafe_copyto!(unsafe_load(Sb.data, col), pointer(tmp[:,col]), n_rows)
+			n_records += n_rows						# Must manually keep track of totals
+			DS.type_ = GMT_READ_DATA
+			unsafe_store!(S, Sb)
+			unsafe_store!(DT.segment, S, _seg)
 		end
-
-		n_records += n_rows						# Must manually keep track of totals
-		DS.type_ = GMT_READ_DATA
-		unsafe_store!(S, Sb)
-		unsafe_store!(DT.segment, S, seg)
 	end
 	DT.n_records, DS.n_records = n_records, n_records	# They are equal because our GMT_DATASET has only one table
 	Dt = unsafe_load(DS.table)

src/gmt_types.jl

@@ -292,7 +292,7 @@ GMTdataset(data::Array{Float32,2}, text::String) =
 GMTdataset(data::Array{Float32,2}) =
 	GMTdataset(data, Float64[], Float64[], DictSvS(), String[], String[], "", String[], "", "", 0, 0)
 
-Base.@kwdef struct GMTfv{T<:AbstractFloat} <: AbstractMatrix{T}
+Base.@kwdef mutable struct GMTfv{T<:AbstractFloat} <: AbstractMatrix{T}
 	verts::AbstractMatrix{T}=Matrix{Float64}(undef,0,0)
 	faces::Vector{<:AbstractMatrix{Int}}=Vector{<:AbstractMatrix{Int}}(undef,0)
 	faces_view::Vector{Matrix{Int}}=Vector{Matrix{Int}}(undef,0)

src/gmtreadwrite.jl

@@ -684,8 +684,14 @@ end
 Read a Wavefront .obj file and return the result in a FaceVertices object.
 """
 function read_obj(fname)
-	n_vert, count_v, count_f, vr, fr = 0, 0, 0, 0, 0
+	# This functions became convoluted because we want to be able to read .obj files that may have
+	# a different number of vertices per face. For example we may have a sphere made of quadrangles
+	# and triangles to close the poles. We deal with these cases by storing the faces matices in a
+	# vector that will have as many elements as there are different geometries. This way we workaround
+	# the limitation that matrices must have the same number of columns, but we add complexity.
+	count_v, count_f, vr, fr = 0, 0, 0, 0
 	first, has_slash = true, false
+	n_vert_vec = Int[]		# Vector to store the number of vertices per each face
 
 	# Do a first round to find out how many vertices and faces there are as well as if indices are simple or with slashes
 	fid = open(fname)
@@ -694,14 +700,32 @@ function read_obj(fname)
 		isempty(it) && continue
 		if (it[1] == 'v' && it[2] == ' ') count_v += 1
 		elseif (it[1] == 'f')
-			first && (has_slash = contains(it, '/'); n_vert = length(split(it))-1;	first = false)
+			first && (has_slash = contains(it, '/'); first = false)
+			push!(n_vert_vec, length(split(it))-1)		# Store the number of vertices of this face
 			count_f += 1
 		end
 	end
 	close(fid)		# Don't know how to rewind 'iter', so close file and open it again.
 
+	u = sort(unique(n_vert_vec))		# Number of faces geometries (tri, quad, etc) in growing order of vertices number
+	n_geoms = length(u)
+	off = minimum(n_vert_vec) - 1
+	n_vert_vec .-= off					# Shift indices so they start at one and we can easily compute a histogram
+	hst = zeros(Int, n_geoms)
+	@inbounds for k = 1:numel(n_vert_vec)  hst[n_vert_vec[k]] += 1  end		# Micro histogram
+	n_vert_vec .+= off
+	F = Vector{Matrix{Int}}(undef, n_geoms)
+	for k = 1:n_geoms
+		F[k] = fill(0, hst[k], n_vert_vec[findfirst(n_vert_vec .== u[k])])	# Initialize each matrix of faces (onre for each geometry type)
+	end
+
+	ind_F = ones(Int, count_f)			# Indices of faces
+	for k = 2:n_geoms
+		ind_F[n_vert_vec .== u[k]] .= k	# For each face, we store the index of the corresponding geometry. That is we know which geometry it belongs to
+	end
+	count_F_rows = zeros(Int, n_geoms)
+
 	V = Matrix{Float64}(undef, count_v, 3)
-	F = Matrix{Int}(undef, count_f, n_vert)
 	fid = open(fname)
 	iter = eachline(fid)
 	for it in iter
@@ -712,14 +736,16 @@ function read_obj(fname)
 			V[vr, 1], V[vr, 2], V[vr, 3] = parse(Float64, line[2]), parse(Float64, line[3]), parse(Float64, line[4])
 		elseif (line[1] == "f")
 			fr += 1
-			if (has_slash)		# Vertex with texture/normals coordinate indices
-				for k = 1:n_vert
+			i_t = ind_F[fr]
+			count_F_rows[i_t] += 1		# 
+			if (has_slash)				# Vertex with texture/normals coordinate indices
+				for k = 1:n_vert_vec[fr]
 					spli = split(line[k+1], "/")
-					F[fr, k] = parse(Int, spli[1])
+					F[i_t][count_F_rows[i_t], k] = parse(Int, spli[1])
 				end
-			else				# simple vertex indices
-				for k = 1:n_vert
-					F[fr, k] = parse(Int, line[k+1])
+			else						# simple vertex indices
+				for k = 1:n_vert_vec[fr]
+					F[i_t][count_F_rows[i_t], k] = parse(Int, line[k+1])
 				end
 			end
 		end

src/psxy.jl

@@ -1438,66 +1438,44 @@ end
     FV, projs = sort_visible_faces(FV, 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 first and hence do not hide the others.
+that the furthest faces are drawn first and hence do not hide the others. NOTE: to not change the original
+FV we store the resultant matix(ces) of faces in a FV0s fiels called "faces_view", which is a vector
+of matrices, one for each geometry (e.g. triangles, quadrangles, etc).
 
 - `FV`: The Faces-Vertices dataset.
 - `azim`: Azimuth angle in degrees. Positive clock-wise from North.
 - `elev`: Elevation angle in degrees above horizontal plane.
 """
-#=
-function sort_visible_faces(FV, azim, elev)::Tuple{Vector{GMTdataset{T, 2} where T<:Real}, Vector{Float64}}
-	V::GMTdataset{Float64,2}, F::GMTdataset{Int,2} = FV[1], FV[2]
-	cos_az, cos_el, sin_az, sin_el = cosd(azim), cosd(elev), sind(azim), sind(elev)
-	view_vec = [sin_az * cos_el, cos_az * cos_el, sin_el]
-
-	n_faces, n_verts = size(F.data, 1), size(F.data, 2)	# Number of faces (polygons) and vertices of the polygons
-	tmp = zeros(n_verts, 3)
-	isVisible = fill(false, n_faces)
-	projs = Float64[]
-	dists = NTuple{2,Float64}[]
-	for face = 1:n_faces
-		for c = 1:3, v = 1:n_verts						# Build the polygon from the FV collection
-			tmp[v,c] = V.data[F.data[face,v], c]
-		end
-		this_proj = dot(facenorm(tmp), view_vec)
-		if ((isVisible[face] = this_proj > 0))
-			cx, cy, cz = sum(tmp[:,1]), sum(tmp[:,2]), sum(tmp[:,3])	# Pseudo-centroids. Good enough for the sorting purpose
-			push!(dists, (cx * sin_az + cy * cos_az, cz * sin_el))
-			push!(projs, this_proj)
-		end
-	end
-	data = F.data[isVisible, :]
-	ind = sortperm(dists)
-	data = data[ind, :]
-	return [V, GMTdataset(data=data)], projs[ind]
-end
-=#
-
 function sort_visible_faces(FV::GMTfv, azim, elev)
 	cos_az, cos_el, sin_az, sin_el = cosd(azim), cosd(elev), sind(azim), sind(elev)
 	view_vec = [sin_az * cos_el, cos_az * cos_el, sin_el]
-
-	n_faces, n_verts = size(FV.faces[1], 1), size(FV.verts, 2)	# Number of faces (polygons) and vertices of the polygons
-	tmp = zeros(n_verts, 3)
-	isVisible = fill(false, n_faces)
 	projs = Float64[]
-	dists = NTuple{2,Float64}[]
-	for face = 1:n_faces
-		for c = 1:3, v = 1:n_verts						# Build the polygon from the FV collection
-			tmp[v,c] = FV.verts[FV.faces[1][face,v], c]
-		end
-		this_proj = dot(facenorm(tmp), view_vec)
-		if ((isVisible[face] = this_proj > 0))
-			cx, cy, cz = sum(tmp[:,1]), sum(tmp[:,2]), sum(tmp[:,3])	# Pseudo-centroids. Good enough for the sorting purpose
-			push!(dists, (cx * sin_az + cy * cos_az, cz * sin_el))
-			push!(projs, this_proj)
+
+	for k = 1:numel(FV.faces)			# Loop over number of face groups (we can have triangles, quads, etc)
+		n_faces, n_verts = size(FV.faces[k], 1), size(FV.verts, 2)	# Number of faces (polygons) and vertices of the polygons
+		tmp = zeros(n_verts, 3)
+		isVisible = fill(false, n_faces)
+		dists = NTuple{2,Float64}[]
+		_projs = Float64[]
+		for face = 1:n_faces
+			for c = 1:3, v = 1:n_verts						# Build the polygon from the FV collection
+				tmp[v,c] = FV.verts[FV.faces[k][face,v], c]
+			end
+			this_proj = dot(facenorm(tmp), view_vec)
+			if ((isVisible[face] = this_proj > 0))
+				cx, cy, cz = sum(tmp[:,1]), sum(tmp[:,2]), sum(tmp[:,3])	# Pseudo-centroids. Good enough for the sorting purpose
+				push!(dists, (cx * sin_az + cy * cos_az, cz * sin_el))
+				push!(_projs, this_proj)
+			end
 		end
+		data = FV.faces[k][isVisible, :]
+		ind = sortperm(dists)
+		data = data[ind, :]
+		(k == 1) ? (FV.faces_view = [data]) : append!(FV.faces_view, [data])
+		projs = (k == 1) ? _projs[ind] : append!(projs, _projs[ind])
 	end
-	data = FV.faces[1][isVisible, :]
-	ind = sortperm(dists)
-	data = data[ind, :]
-	FV.faces[1] = data
-	return FV, projs[ind]
+
+	return FV, projs
 end
 
 # ---------------------------------------------------------------------------------------------------
@@ -1652,34 +1630,40 @@ end
 function replicant_worker(FV::GMTfv, xyz, azim, elev, cval, cpt, scales)::Vector{GMTdataset}
 	# This guy is the one who does the replicant work
 
-	FV, normals = sort_visible_faces(FV, azim, elev)	# Sort & kill invisible
+	FV, normals = sort_visible_faces(FV, azim, elev)	# Return a modified FV containing info about the sorted visible faces.
 
-	n_faces = size(FV.faces[1], 1)			# Number of faces of base body
-	n_rows = size(FV.faces[1], 2)			# Number of rows (vertices of the polygon)
+	n_faces_tot = sum(size.(FV.faces_view, 1))			# Total number of segments or faces (polygons)
 
 	# ---------- First we convert the FV into a vector of GMTdataset
-	D1 = Vector{GMTdataset}(undef, n_faces)
-	t = zeros(eltype(FV.verts), n_rows, 3)
-	for face = 1:n_faces					# Loop over faces
-		for c = 1:3, r = 1:n_rows
-			t[r,c] = FV.verts[FV.faces[1][face, r], c]
+	D1 = Vector{GMTdataset}(undef, n_faces_tot)
+	t = zeros(eltype(FV.verts), maximum(size.(FV.faces_view,2)), 3)	# The maximum number of vertices in any polygon of all geometries
+	count_face = 0
+
+	for geom = 1:numel(FV.faces_view)					# If we have more than one type of geometries (e.g. triangles and quadrangles)
+		n_faces = size(FV.faces_view[geom], 1)			# Number of segments or faces (polygons) in this geometry
+		n_rows  = size(FV.faces_view[geom], 2)			# Number of rows (vertices of the polygon) in this geometry
+		for face = 1:n_faces							# Loop over faces
+			count_face += 1								# Counter that keeps track on the current number of polygon. 
+			for c = 1:3, r = 1:n_rows
+				t[r,c] = FV.verts[FV.faces_view[geom][face, r], c]
+			end
+			D1[count_face] = GMTdataset(data=copy(t))
 		end
-		D1[face] = GMTdataset(data=copy(t))
 	end
 
-	P::Ptr{GMT_PALETTE} = palette_init(G_API[1], cpt)		# A pointer to a GMT 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)
+	D2 = Vector{GMTdataset}(undef, size(xyz, 1) * n_faces_tot)
 
 	# ---------- Now we do the replication
-	for k = 1:size(xyz, 1)					# Loop over number of positions. For each of these we have a new body
+	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
+		for face = 1:n_faces_tot						# 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 * scales[k] .+ xyz[k:k,1:3]), header=txt)
+			D2[(k-1)*n_faces_tot+face] = GMTdataset(data=(D1[face].data * scales[k] .+ xyz[k:k,1:3]), header=txt)
 		end
 	end
 	return set_dsBB(D2)

test/test_solidos.jl

@@ -5,7 +5,6 @@
 	GMT.dodecahedron()
 	GMT.tetrahedron()
 	GMT.cube()
-	GMT.sphere()
 
 	FV = sphere();
 	D  = GMT.replicant(FV, replicate=(centers=rand(10,3), scales=0.1));