Add a 'extrude' function to vertically extrude shapes, and a cylinder…

… FV body. (#1564)

* Add a cross2 version that deals with the cross product of 2D vectors.

* Update the fv2fv() for the new GMTfv members.

* Improve a docstring.

* Export cylinder and extrude

* Add the 'bfculling' member to GMTfv type.

* When FV.bfculling is set to false do not delete the (maybe) invisible faces.

* Add a 'extrude' function to vertically extrude shapes, and a cylinder FV body.

src/GMT.jl

@@ -184,9 +184,9 @@ export
 
 	lazinfo, lazread, lazwrite, lasread, laswrite,
 
-	cube, dodecahedron, icosahedron, sphere, octahedron, tetrahedron, torus, replicant,
+	cube, cylinder, dodecahedron, icosahedron, sphere, octahedron, tetrahedron, torus, replicant,
 
-	df2ds, ds2df, ODE2ds, fv2fv, surf2fv,
+	df2ds, ds2df, extrude, fv2fv, surf2fv, ODE2ds,
 	@?, @dir
 
 include("common_docs.jl")

src/gmt_types.jl

@@ -178,7 +178,7 @@ Base.size(C::GMTcpt) = size(C.range, 1)
 Base.isempty(C::GMTcpt) = (size(C) == 0)
 
 """
-mutable struct GMTps
+    mutable struct GMTps
 
 The fields of this struct are:
 - `postscript::String`:         Actual PS plot (text string)
@@ -292,12 +292,29 @@ GMTdataset(data::Array{Float32,2}, text::String) =
 GMTdataset(data::Array{Float32,2}) =
 	GMTdataset(data, Float64[], Float64[], DictSvS(), String[], String[], "", String[], "", "", 0, 0)
 
+"""
+    mutable struct GMTfv{T<:AbstractFloat} <: AbstractMatrix{T}
+
+The GMTfv struct is used to store a triangulated mesh.
+
+The fields of this struct are:
+- `verts::AbstractMatrix{T}`:             Mx3 Matrix with the data vertices
+- `faces`:
+- `faces_view`:
+- `bbox`:
+- `zscale`:
+- `bfculling`:
+- `proj4`:
+- `wkt`:
+- `eps`:
+"""
 Base.@kwdef mutable struct GMTfv{T<:AbstractFloat} <: AbstractMatrix{T}
 	verts::AbstractMatrix{T}=Matrix{Float64}(undef,0,0)
 	faces::Vector{<:AbstractMatrix{<:Integer}}=Vector{Matrix{Int}}(undef,0)
 	faces_view::Vector{Matrix{Int}}=Vector{Matrix{Int}}(undef,0)
 	bbox::Vector{Float64}=zeros(6)
 	zscale::Float64=1.0
+	bfculling::Bool=true
 	proj4::String=""
 	wkt::String=""
 	epsg::Int=0

src/proj_utils.jl

@@ -212,15 +212,18 @@ end
 # -------------------------------------------------------------------------------------------------
 """
     circgeo(lon, lat; radius=X, proj="", s_srs="", epsg=0, dataset=false, unit=:m, np=120, shape="")
-or
+Or
 
     circgeo(lonlat; radius=X, proj="", s_srs="", epsg=0, dataset=false, unit=:m, np=120, shape="")
 
-Args:
+Compute a geographical circle (and other shapes) in geographical or in projected coordinates.
 
+### Args
 - `lonlat`:   - longitude, latitude (degrees). If a Mx2 matrix, returns as many segments as number of rows.
                 Use this to compute multiple shapes at different positions. In this case output type is
-				always a vector of GMTdatasets.
+                always a vector of GMTdatasets.
+
+### Kwargs
 - `radius`:   - The circle radius in meters (but see `unit`) or circumscribing circle for the other shapes
 - `proj` or `s_srs`:  - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKT
 - `epsg`:     - Alternative way of specifying the projection [Default is WGS84]
@@ -233,9 +236,11 @@ Args:
 ### Returns
 - circ - A Mx2 matrix or GMTdataset with the circle coordinates
 
-## Example: Compute circle about the (0,0) point with a radius of 50 km
-
-    c = circgeo([0.,0], radius=50, unit=:k)
+### Example:
+Compute a circle about the (0,0) point with a radius of 50 km
+```julia
+    c = circgeo([0. 0], radius=50, unit=:k)
+```
 """
 circgeo(lon::Real, lat::Real; radius::Real=0., proj::String="", s_srs::String="", epsg::Integer=0, dataset::Bool=false, unit=:m, np::Int=120, shape="") =
 	circgeo([lon lat]; radius=radius, proj=proj, s_srs=s_srs, epsg=epsg, dataset=dataset, unit=unit, np=np, shape=shape)

src/psxy.jl

@@ -392,8 +392,8 @@ function deal_faceverts(arg1::GMTfv, d; del::Bool=true)::GMTfv
 	azim, elev = get_numeric_view()
 	arg1, dotprod = sort_visible_faces(arg1, azim, elev; del=del)		# Sort & kill (or not) 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
+		is_in_dict(d, [:Z :level :levels]) === nothing && (d[:Z] = abs.(dotprod))
+		(is_in_dict(d, CPTaliases) === nothing) && (d[:C] = gmt("makecpt -T0/1 -C140,220"))	# Users may still set a CPT
 	end
 	return arg1
 end
@@ -1450,16 +1450,21 @@ that the furthest faces are drawn first and hence do not hide the others. NOTE:
 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).
 
+### Args
 - `FV`: The Faces-Vertices dataset.
 - `azim`: Azimuth angle in degrees. Positive clock-wise from North.
 - `elev`: Elevation angle in degrees above horizontal plane.
-- `del`: Boolean to control whether to delete invisible faces. True by default
+
+### Kwargs
+- `del`: Boolean to control whether to delete invisible faces. True by default. But this can be
+  overwriten by the value of the ``bfculling`` member of the FV object.
 """
 function sort_visible_faces(FV::GMTfv, azim, elev; del::Bool=true)::Tuple{GMTfv, Vector{Float64}}
 	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]
 	projs = Float64[]
 
+	!FV.bfculling && (del = false)		# Do not delete if bfculling is set to false (for example if FV is not closed)
 	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)
@@ -1483,7 +1488,7 @@ function sort_visible_faces(FV::GMTfv, azim, elev; del::Bool=true)::Tuple{GMTfv,
 		(k == 1) ? (FV.faces_view = [data]) : append!(FV.faces_view, [data])
 		projs = (k == 1) ? _projs[ind] : append!(projs, _projs[ind])
 	end
-	sum(size.(FV.faces_view, 1)) < 3 * sum(size.(FV.faces, 1)) &&
+	sum(size.(FV.faces_view, 1)) < sum(size.(FV.faces, 1) / 3) &&
 		@warn("More than 2/3 of the faces found invisible. This often indicates that the Z and X,Y units are not the same. Consider using the `zscale` field of the `FV` input.")
 
 	return FV, projs

src/solids.jl

@@ -3,7 +3,7 @@
 """
     FV = icosahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0))
 
-Creates an icosahedron mesh with radius `r`. 
+Create 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)`.
@@ -59,7 +59,7 @@ end
 """
     FV = octahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0))
 
-Creates an octahedron mesh with radius `r`. 
+Create 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)`.
@@ -94,7 +94,7 @@ end
 """
     FV = dodecahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0))
 
-Creates an dodecahedron mesh with radius `r`. 
+Create 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)`.
@@ -151,7 +151,7 @@ end
 """
     FV = tetrahedron(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0))
 
-Creates a tetrahedron mesh with radius `r`. 
+Create 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)`.
@@ -183,7 +183,7 @@ end
 """
     FV = cube(r=1.0; radius=1.0, origin=(0.0, 0.0, 0.0))
 
-Creates a cube mesh with radius `r`. 
+Create 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)`.
@@ -217,9 +217,11 @@ end
 
 # ----------------------------------------------------------------------------
 """
-    FV = sphere(r=1; n=1, radius=1.0, center=(0.0, 0.0, 0.0))
+    FV = sphere(r=1; radius=1.0, n=1, center=(0.0, 0.0, 0.0))
 
-Generate a geodesic sphere triangulation based on the number of refinement iterations `n`
+Create a triangulated geodesic sphere. 
+
+Generates a geodesic sphere triangulation based on the number of refinement iterations `n`
 and the radius `r`. Geodesic spheres (aka Buckminster-Fuller spheres) are triangulations
 of a sphere that have near uniform edge lenghts.  The algorithm starts with a regular
 icosahedron. Next this icosahedron is refined `n` times, while nodes are pushed to a sphere
@@ -258,9 +260,9 @@ end
 
     V, F = subTriSplit(V, F, n=1)
 
-Splits the triangulation defined by the faces F, and the vertices V, n times. Each triangle is
-linearly split into 4 triangles with each iterations.
+Split the triangulation defined by the faces F, and the vertices V, n times.
 
+Each triangle is linearly split into 4 triangles with each iterations.
 First mode ingests a two elements vector of GMTdataset where first contains the vertices and the second
 the indices that define the faces and returns a same type. The second mode expects the vertices and faces
 as two separate arrays (but it also accepts GMTdatasets) and returns two matrices with the vertices and faces.
@@ -313,11 +315,11 @@ end
 """
     FV = torus(; r=2.0, R=5.0, center=(0.0, 0.0, 0.0), nx=100, ny=50) -> GMTfv
 
-Creates a torus mesh with radius `r`. 
+Create a torus mesh with radius `r`. 
 
 - `r`: the inner radius of the torus.
 - `R`: the outer radius of the torus.
-- `center`: A tuple of three numbers defining the origin of the body. Default is `(0.0, 0.0, 0.0)`.
+- `center`: A 3-element array or tuple of three numbers defining the origin of the body. Default is `(0.0, 0.0, 0.0)`.
 - `nx`: the number of vertices in the xx direction.	
 - `ny`: the number of vertices in the yy direction.
 """
@@ -329,3 +331,88 @@ function torus(; r=2.0, R=5.0, center=(0.0, 0.0, 0.0), nx=100, ny=50)::GMTfv
 	z = [r*sin(v) + center[3] for u in Θ, v in ϕ]
 	surf2fv(x, y, z)
 end
+
+# ----------------------------------------------------------------------------
+"""
+    FV = extrude(shape::Matrix{<:AbstractFloat}, h; base=0.0, closed=true) -> GMTfv
+
+Create an extruded 2D/3D shape.
+
+### Args
+- `shape`: The shape to extrude.
+- `h`: The height of the extrusion. Same units as in `shape`.
+
+### Kwargs
+- `base`: The base height of the 2D shape to extrude. Default is 0. Ignored if the shape is a 3D polygon.
+- `closed`: If true (the default), close the shapre at top and bottom.
+"""
+function extrude(shape::Matrix{<:AbstractFloat}, h; base=0.0, closed=true)::GMTfv
+	np = size(shape, 1)
+	if (size(shape, 2) < 3)			# 2D polygon
+		V = [shape fill(convert(eltype(shape), base), np); shape fill(convert(eltype(shape), h), np)]
+	else
+		V = [shape; shape .+ [0 0 convert(eltype(shape), h)]]
+	end
+
+	if (h * facenorm(shape, normalize=false)[3] > 0)  v1, v2, f2, f4 = np:-1:1,     np+1:2np, 1, np
+	else                                              v1, v2, f2, f4 = 2np:-1:np+1, 1:np,     np, 1
+	end
+
+	if (size(shape, 1) == 5 && shape[1,:] ≈ shape[end,:])		# When shape is a quad
+		F = (closed == 1) ? [[reshape([v1;], 1, :); reshape([v2;], 1, :); zeros(Int, np-1, 5)]] : [zeros(Int, np-1, 5)]
+		k = (closed == 1) ? 2 : 0		# Starting index of the faces vector that will contain the vertical faces
+		for n = 1:np-1					# Create vertical faces
+			k += 1 
+			F[1][k, 1], F[1][k, 2], F[1][k, 3], F[1][k, 4], F[1][k, 5] = n, n+f2, n+np+1, n+f4, n
+		end
+	else
+		if (closed == 1)  F = [[reshape([v1;], 1, :); reshape([v2;], 1, :)], zeros(Int, np-1, 5)]
+		else              F = [zeros(Int, np-1, 5)]
+		end
+		iF = (closed == 1) ? 2 : 1		# Index of the faces vector that will contain the vertical faces
+		for n = 1:np-1					# Create vertical faces
+			F[iF][n, 1], F[iF][n, 2], F[iF][n, 3], F[iF][n, 4], F[iF][n, 5] = n, n+f2, n+np+1, n+f4, n
+		end
+	end
+	fv2fv(F, V; bfculling=(closed == 1))
+end
+
+# ----------------------------------------------------------------------------
+extrude(shape::GMTdataset, h; base=0.0, closed=true)::GMTfv = extrude(shape.data, h; base=base, closed=closed)
+
+# ----------------------------------------------------------------------------
+"""
+    FV = cylinder(r, h; base=0.0, center=(0.0, 0.0, 0.0), geog=false, unit="m", np=36) -> GMTfv
+
+Create a cylinder with radius `r` and height `h`.
+
+### Args
+- `r`: The radius of the cylinder. For geographical cylinders, the default is meters. But see `unit` below.
+- `h`: The height of the cylinder. It should be in the same unit as `r`.
+
+### Kwargs
+- `base`: The base height of the cylinder. Default is 0.
+- `center`: A 3-element array or tuple of three numbers defining the origin of the body. Default is `(0.0, 0.0, 0.0)`.
+- `closed`: If true (the default), close the cylinder at top and bottom.
+- `geog`: If true, create a cylinder in geographical coordinates.
+- `unit`: For geographical cylinders only.If radius is not in meters use one of `unit=:km`, or `unit=:Nautical` or `unit=:Miles`
+- `np`: The number of vertices in the circle. Default is 36.
+
+Return a Faces-Vertices dataset.
+
+### Example
+```julia
+	FV = cylinder(50, 100)
+	viz(FV)
+```
+"""
+function cylinder(r, h; base=0.0, center=(0.0, 0.0, 0.0), closed=true, geog=false, unit="m", np=36)::GMTfv
+	h0 = (base != 0.0) ? base : length(center) == 3 ? center[3] : 0.0
+	if (geog == 1)
+		xy = circgeo(center[1], center[2]; radius=r, unit=unit)
+	else
+		t = linspace(0, 2pi, np)
+		xy = [(center[1] .+ r * cos.(t)) (center[2] .+ r * sin.(t))]
+	end
+	extrude(xy, h; base=h0, closed=closed)
+end

src/utils.jl

@@ -1073,7 +1073,7 @@ end
 """
     overlap, value = rect_overlap(xc_1, yc_1, xc_2, yc_2, width1, height1, width2, height2) -> Bool, Float64
 
-Checks if two rectangles, aligned with the axes, overlap.
+Check if two rectangles, aligned with the axes, overlap.
 
 - `xc_1, yc_1`: Center of the first rectangle
 - `xc_2, yc_2`: Center of the second rectangle
@@ -1099,7 +1099,7 @@ end
 """
     n = facenorm(M::Matrix{<:Real}; normalize=true, zfact=1.0)
 
-Calculates the normal vector of a polygon with vertices in `M`.
+Calculate the normal vector of a polygon with vertices in `M`.
 
 - `normalize`: By default, it returns the unit vector. If `false`, it returns the non-normalized 
   vector that represents the area of the polygon.
@@ -1115,9 +1115,9 @@ function facenorm(M::Matrix{<:Real}; zfact=1.0, normalize=true)
 	p1 = M[1,:]
 	last = (p1 == M[end,:]) ? size(M,1)-1 : size(M,1)
 	if (zfact == 1.0)
-		c = cross(M[last,:], p1)			# First edge
+		c = cross2(M[last,:], p1)			# First edge
 		for k = 1:last-1					# Loop from first to end-1
-			c += cross(M[k,:], M[k+1,:])	# Sum the rest of edges
+			c += cross2(M[k,:], M[k+1,:])	# Sum the rest of edges
 		end
 	else
 		p2 = M[last,:]
@@ -1132,6 +1132,17 @@ function facenorm(M::Matrix{<:Real}; zfact=1.0, normalize=true)
 	return normalize ? c / norm(c) : c	# The cross product gives us 2 * A(rea)
 end
 
+# Like the cross() function from LinearAlgebra but deals also with 2D vectors, though it returns a 3D vector
+function cross2(a::AbstractVector, b::AbstractVector)
+	(length(a) == length(b) == 2) && return [0, 0, a[1]*b[2] - a[2]*b[1]]
+	if !(length(a) == length(b) == 3)
+		throw(DimensionMismatch("cross product is only defined for vectors of length 3"))
+	end
+	a1, a2, a3 = a
+	b1, b2, b3 = b
+	[a2*b3-a3*b2, a3*b1-a1*b3, a1*b2-a2*b1]
+end
+
 # ---------------------------------------------------------------------------------------------------
 """
     A, B, C, D = eq_plane(azim, elev, dist)

src/utils_types.jl

@@ -701,13 +701,15 @@ Create a FacesVerices object from a matrix of faces indices and another matrix o
 - `wkt`:  A WKT SRS.
 - `epsg`: Same as `proj` but using an EPSG code
 """
-function fv2fv(F::Vector{<:AbstractMatrix{<:Integer}}, V; hdr="", comment=String[], proj="", proj4="", wkt="", epsg=0)::GMTfv
+function fv2fv(F::Vector{<:AbstractMatrix{<:Integer}}, V; zscale=1.0, bfculling=true, proj="", proj4="", wkt="", epsg=0)::GMTfv
 	(isempty(proj4) && !isempty(proj)) && (proj4 = proj)	# Allow both proj4 or proj keywords
 	bbox = extrema(V, dims=1)
-	GMTfv(verts=V, faces=F, bbox=[bbox[1][1], bbox[1][2], bbox[2][1], bbox[2][2], bbox[3][1], bbox[3][2]], proj4=proj4, wkt=wkt, epsg=epsg)
+	GMTfv(verts=V, faces=F, bbox=[bbox[1][1], bbox[1][2], bbox[2][1], bbox[2][2], bbox[3][1], bbox[3][2]],
+	      zscale=zscale, bfculling=bfculling, proj4=proj4, wkt=wkt, epsg=epsg)
 end
 
-fv2fv(F::Matrix{<:Integer}, V; proj="", proj4="", wkt="", epsg=0) = fv2fv([F], V; proj=proj, proj4=proj4, wkt=wkt, epsg=epsg)
+fv2fv(F::Matrix{<:Integer}, V; zscale=1.0, bfculling=true, proj="", proj4="", wkt="", epsg=0) =
+	fv2fv([F], V; zscale=zscale, bfculling=bfculling, proj=proj, proj4=proj4, wkt=wkt, epsg=epsg)
 
 """
 When using Meshing.jl we can use the output of the ``isosurface`` function, "verts, faces" as input to this function.
@@ -727,10 +729,11 @@ FV = fv2fv(fcs, vts)
 viz(FV, cmap=makecpt(T="0/1", cmap="darkgreen,lightgreen"))
 ```
 """
-function fv2fv(F::Vector{Tuple{Int, Int, Int}}, V::Vector{Tuple{Float64, Float64, Float64}}; proj="", proj4="", wkt="", epsg=0)::GMTfv
-	verts=reshape(reinterpret(Float64, V), (3,:))'
-	faces=[reshape(reinterpret(Int, F), (3,:))']
-	fv2fv(faces, verts; proj=proj, proj4=proj4, wkt=wkt, epsg=epsg)
+function fv2fv(F::Vector{Tuple{Int, Int, Int}}, V::Vector{Tuple{Float64, Float64, Float64}};
+               zscale=1.0, bfculling=true, proj="", proj4="", wkt="", epsg=0)::GMTfv
+	verts = reshape(reinterpret(Float64, V), (3,:))'
+	faces = [reshape(reinterpret(Int, F), (3,:))']
+	fv2fv(faces, verts; zscale=zscale, bfculling=bfculling, proj=proj, proj4=proj4, wkt=wkt, epsg=epsg)
 end
 
 # ---------------------------------------------------------------------------------------------------