From 5fefd9d53df5b2edfe68211cd480422c8c9024c5 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 18 Oct 2024 19:18:56 +0100 Subject: [PATCH 1/5] Complicate the read_obj() function to allow reading .obj with more than one geometry. --- src/gmtreadwrite.jl | 44 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 35 insertions(+), 9 deletions(-) diff --git a/src/gmtreadwrite.jl b/src/gmtreadwrite.jl index 0a25e3390..13e4c23ce 100644 --- a/src/gmtreadwrite.jl +++ b/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 From cf251840b7eaf4559b7e178e649626a328b2848a Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 18 Oct 2024 19:19:49 +0100 Subject: [PATCH 2/5] Make the GMTfv struct mutable. --- src/gmt_types.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gmt_types.jl b/src/gmt_types.jl index 1b99d7456..b44c48ba0 100644 --- a/src/gmt_types.jl +++ b/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) From cd5a6f0cf46afc63d867d575fb611033398c42fe Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 18 Oct 2024 19:20:51 +0100 Subject: [PATCH 3/5] Rework the dataset_init_FV() to deal with FV types with multiple geometries. --- src/gmt_main.jl | 52 ++++++++++++++++++++++++++++--------------------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/src/gmt_main.jl b/src/gmt_main.jl index 3b990c380..5c338be8b 100644 --- a/src/gmt_main.jl +++ b/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) From f1b001a2b026abf6087ff14ee840ae6abbf63ad7 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 18 Oct 2024 19:22:22 +0100 Subject: [PATCH 4/5] Re-writings to allow FV's with multiple geometries. --- src/psxy.jl | 110 ++++++++++++++++++++++------------------------------ 1 file changed, 47 insertions(+), 63 deletions(-) diff --git a/src/psxy.jl b/src/psxy.jl index 5cd448f45..d534aa3d3 100644 --- a/src/psxy.jl +++ b/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) From aedbd1f6ca8979fd4264e84ad2878986600d50aa Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 18 Oct 2024 19:22:45 +0100 Subject: [PATCH 5/5] Delete repeated test. --- test/test_solidos.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_solidos.jl b/test/test_solidos.jl index 32e5d2954..6434e2ff4 100644 --- a/test/test_solidos.jl +++ b/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));