👈🏿 Not following column-major layout

ZenanH · ZenanH · commit 3cfb136ebeb4 · 2025-08-10T22:08:27.000+02:00
diff --git a/example/2d.jl b/example/2d.jl
@@ -1,7 +1,7 @@
 using FastPointQuery
 using WGLMakie
 
-points = rand(2, 100)
+points = rand(100, 2)
 poly = get_polygon(points, ratio=0.1)
 stl_file = joinpath(@__DIR__, "2d_hole.stl")
 stl_model = readSTL2D(stl_file)
@@ -14,9 +14,9 @@ write_polygon(poly, joinpath(@__DIR__, "polygon.geojson"))
 write_polygon(joinpath(@__DIR__, "polygon.geojson"), poly)
 poly2 = read_polygon(joinpath(@__DIR__, "polygon.geojson"))
 
-pts1 = get_pts(poly, 0.1)
-pts2 = get_pts(poly, 0.1; fill=false)
-pts3 = get_pts(stl_model, 0.01, fill=false)
+pts1 = FastPointQuery._get_pts(poly, 0.1, true, true)
+pts2 = FastPointQuery._get_pts(poly, 0.1, false, false)
+pts3 = FastPointQuery._get_pts(stl_model, 0.01, false, true)
 
 let
     fig = Figure()
diff --git a/example/3d.jl b/example/3d.jl
@@ -4,7 +4,7 @@ using WGLMakie
 stl_path1 = get_resource().dragon_phoenix
 stl_path2 = get_resource().wheel
 
-points = rand(3, 10000)
+points = rand(10000, 3)
 stl_model1 = readSTL3D(stl_path1)
 stl_model2 = readSTL3D(stl_path2)
 
@@ -13,7 +13,7 @@ pip_query(stl_model2, points; nsamples=5)
 
 let
     h = 1
-    pts = FastPointQuery._get_pts_voxel(stl_model1, h)
+    pts = FastPointQuery._get_pts_voxel(stl_model1, h, true)
     fig = Figure()
     ax = LScene(fig[1, 1])
     scatter!(ax, pts, markersize=1)
@@ -22,7 +22,7 @@ end
 
 let
     h = 5
-    pts = FastPointQuery._get_pts_voxel(stl_model2, h; fill=true)
+    pts = FastPointQuery._get_pts_voxel(stl_model2, h, true)
     fig = Figure()
     ax = LScene(fig[1, 1])
     scatter!(ax, pts, markersize=1)
diff --git a/example/fileio.jl b/example/fileio.jl
@@ -1,7 +1,7 @@
 using FastPointQuery
 using WGLMakie
 
-points = rand(2, 100)
+points = rand(100, 2)
 saveply(points, joinpath(@__DIR__, "pointcloud.ply"))
 readply(joinpath(@__DIR__, "pointcloud.ply"), xy=true)
 
diff --git a/example/getpts.jl b/example/getpts.jl
@@ -4,7 +4,7 @@ using WGLMakie
 stl_model = readSTL3D(joinpath(@__DIR__, "wheel.stl"))
 
 h = 3
-pts = FastPointQuery._get_pts_ray(stl_model, h, ϵ="FP32")
+pts = FastPointQuery._get_pts_ray(stl_model, h, true)
 
 saveply(joinpath(@__DIR__, "wheel.ply"), pts)
 
diff --git a/src/FastPointQuery.jl b/src/FastPointQuery.jl
@@ -23,7 +23,6 @@ const pyfun = PythonCall.pybuiltins
 const res_dir = joinpath(@__DIR__, "../example")
 
 function __init__()
-    @info "initializing environment..."
     try # import Python modules
         PythonCall.pycopy!(np        , PythonCall.pyimport("numpy"       ))
         PythonCall.pycopy!(shapely   , PythonCall.pyimport("shapely"     ))
diff --git a/src/discretization/_polygon.jl b/src/discretization/_polygon.jl
@@ -30,9 +30,9 @@ function _get_pts(polygon::QueryPolygon, h::Real, fill::Bool, edge::Bool)
         ur = pts_cen + np.array([ dx,  dy])
         dl = pts_cen + np.array([-dx, -dy])
         dr = pts_cen + np.array([ dx, -dy])
-        return PythonCall.PyArray(np.vstack([ul, ur, dl, dr]).T, copy=false)
+        return PythonCall.PyArray(np.vstack([ul, ur, dl, dr]), copy=false)
     else
-        return PythonCall.PyArray(pts_cen.T, copy=false)
+        return PythonCall.PyArray(pts_cen, copy=false)
     end
 end
 
diff --git a/src/discretization/_polyhedron.jl b/src/discretization/_polyhedron.jl
@@ -8,10 +8,10 @@
 function _get_pts_voxel(stl_model::STLInfo3D, h::Real, fill::Bool)
     # inputs check
     h > 0 || error("h must be a positive number")
-    step = h * 0.25
-    offsets = np.array([-1  1  1; -1  1 -1; 1  1  1; 1  1 -1
-                        -1 -1  1; -1 -1 -1; 1 -1  1; 1 -1 -1]) * step
-    py_all = pyslice(nothing, nothing, nothing)
+    # step = h * 0.25
+    # offsets = np.array([-1  1  1; -1  1 -1; 1  1  1; 1  1 -1
+    #                     -1 -1  1; -1 -1 -1; 1 -1  1; 1 -1 -1]) * step
+    # py_all = pyslice(nothing, nothing, nothing)
 
     # convert to trimesh
     vertices = stl_model.py_vertices
@@ -23,17 +23,18 @@ function _get_pts_voxel(stl_model::STLInfo3D, h::Real, fill::Bool)
     pts_cen = filled_vox.points
 
     if fill # filling mode
-        expanded = pts_cen[py_all, np.newaxis, py_all] + offsets[np.newaxis, py_all, py_all]
-        return PyArray(expanded.reshape(-1, 3))
+        # expanded = pts_cen[py_all, np.newaxis, py_all] + offsets[np.newaxis, py_all, py_all]
+        # PyArray(expanded.reshape(-1, 3))
+        return filling_pts(PyArray(pts_cen), h)
     else
         return PyArray(pts_cen)
     end
 end
 
-function _get_pts_ray(stl_model::STLInfo3D, h::Real, fill::Bool, ϵ::String)
+function _get_pts_ray(stl_model::STLInfo3D, h::Real, fill::Bool)
     # inputs check
     h > 0 || error("h must be a positive number")
-    points, pts2 = prepareprojection(stl_model, h, ϵ)
+    points, pts2 = prepareprojection(stl_model, h)
     @info "generating pts at h = $h"
     edges = projectionlist(stl_model, points)
     pts_cen = fill_particles(edges, pts2, h)
@@ -44,19 +45,19 @@ function _get_pts_ray(stl_model::STLInfo3D, h::Real, fill::Bool, ϵ::String)
     end
 end
 
-function prepareprojection(stl_model::STLInfo3D, h::Real, ϵ)
-    T = ϵ == "FP32" ? Float32 : Float64
+function prepareprojection(stl_model::STLInfo3D, h::Real)
+    T = Float64
     pts2 = meshbuilder(stl_model.vmin[1]-1.5h : h : stl_model.vmax[1]+1.5h,
-                       stl_model.vmin[2]-1.5h : h : stl_model.vmax[2]+1.5h, ϵ=ϵ)
+                       stl_model.vmin[2]-1.5h : h : stl_model.vmax[2]+1.5h)
     vzlimit = T(stl_model.vmin[3] - 10h)
-    points = vcat(pts2, vzlimit .* ones(T, 1, size(pts2, 2)))
+    points = hcat(pts2, vzlimit .* ones(T, size(pts2, 1)))
     return points, pts2
 end
 
 function projectionlist(stl_model::STLInfo3D, points::AbstractArray{T}) where T
     # inputs check
-    pts = points'; n, m = size(pts)
-    m == 3 || error("points must be a 3xN array")
+    pts = points; n, m = size(pts)
+    m == 3 || error("points must be a Nx3 array")
 
     mesh = stl_model.mesh
     pts = np.asarray(pts, dtype=np.float32)
@@ -88,7 +89,7 @@ function fill_particles(edges::AbstractArray, pts::AbstractMatrix{T}, h::Real) w
     @inbounds for i in eachindex(edges)
         v = edges[i]
         if !isempty(v) && iseven(length(v)) && any(!iszero, v)
-            sort!(v); xi, yi = pts[1, i], pts[2, i]
+            sort!(v); xi, yi = pts[i, 1], pts[i, 2]
             for j in 1:2:(length(v) - 1)
                 for z in v[j]:h:v[j + 1]
                     push!(x_all, xi)
@@ -98,9 +99,9 @@ function fill_particles(edges::AbstractArray, pts::AbstractMatrix{T}, h::Real) w
             end
         end
     end
-    xyz = Matrix{T}(undef, 3, length(x_all))
-    xyz[1, :] .= x_all
-    xyz[2, :] .= y_all
-    xyz[3, :] .= z_all
+    xyz = Matrix{T}(undef, length(x_all), 3)
+    xyz[:, 1] .= x_all
+    xyz[:, 2] .= y_all
+    xyz[:, 3] .= z_all
     return xyz
 end
diff --git a/src/fileio/asc.jl b/src/fileio/asc.jl
@@ -49,7 +49,7 @@ function readasc(filename::String)
 
             # 拼成 (N, 3) 的数组
             xyz = np.column_stack([xs, ys, zs])
-            return xyz.T
+            return xyz
     """ => py_tmp
     
     return py2ju(Array, py_tmp(filename, rasterio, np))
diff --git a/src/fileio/geojson.jl b/src/fileio/geojson.jl
@@ -43,7 +43,7 @@ Write the polygon to a GeoJSON file. The `polygon` is an instance of `QueryPolyg
 Example:
 ---
 ```julia
-polygon_xy = [0 0; 1 0; 1 1; 0 1]' # 2xN array
+polygon_xy = [0 0; 1 0; 1 1; 0 1] # Nx2 array
 poly = get_polygon(polygon_xy, ratio=1) # poly.polygon to visualize
 write_polygon(poly, "/path/to/polygon.geojson")
 # or
diff --git a/src/fileio/ply.jl b/src/fileio/ply.jl
@@ -17,26 +17,26 @@ Saves a point cloud to a PLY file.
 Example:
 ---
 ```julia
-points = rand(3, 1000) # 3D points
+points = rand(1000, 3) # 3D points
 saveply(points, "pointcloud.ply")
 # or
 saveply("pointcloud", points) # add .ply extension automatically
 
-points = rand(2, 1000) # 2D points
+points = rand(1000, 2) # 2D points
 saveply("pointcloud.ply", points)
 ```
 """
 function saveply(points::AbstractArray, filename::String)
     # input check
-    n, m = size(points)
-    n in [2, 3] || throw(ArgumentError("points must be a 2D array with 2 or 3 rows"))
+    m, n = size(points)
+    n in [2, 3] || throw(ArgumentError("points must be a 2D array with 2 or 3 columns"))
     if !endswith(filename, ".ply")
         filename *= ".ply"
     end
     if n == 2
-        pts = vcat(points, zeros(1, m))'
+        pts = hcat(points, zeros(m))
     else
-        pts = points'
+        pts = points
     end
 
     cloud = trimesh.points.PointCloud(pts)
@@ -65,9 +65,9 @@ pc = readply("pointcloud.ply", xy=true) # returns only x and y
 """
 function readply(filename::String; xy::Bool=false)
     isfile(filename) || throw(ArgumentError("file $(filename) does not exist"))
-    pc = py2ju(Array, trimesh.load(filename).vertices.T)
+    pc = py2ju(Array, trimesh.load(filename).vertices)
     if xy
-        return Array(pc[1:2, :])
+        return Array(pc[:, 1:2])
     else
         return pc
     end
diff --git a/src/fileio/tiff.jl b/src/fileio/tiff.jl
@@ -12,7 +12,7 @@ export readtiff
 Description:
 ---
 Reads a GeoTIFF file and returns the pixel coordinates and values as a 3D array
-with shape (3, N), where N is the number of pixels. `dst_crs` specifies the destination 
+with shape (N, 3), where N is the number of pixels. `dst_crs` specifies the destination 
 coordinate reference system (CRS) for the output coordinates. The default CRS is "EPSG:3857" 
 (Web Mercator).
 
@@ -63,7 +63,7 @@ function readtiff(filename::String; dst_crs::String="EPSG:3857")
             xs = np.array(xs)
             ys = np.array(ys)
             xyz = np.column_stack([xs, ys, zs])
-            return xyz.T
+            return xyz
     """ => py_tmp
 
     return py2ju(Array, py_tmp(filename, dst_crs, rasterio, np))
diff --git a/src/polygon.jl b/src/polygon.jl
@@ -20,7 +20,7 @@ edge of the polygon, you can set `edge=true` to check if it is considered as ins
 Example:
 ---
 ```julia
-polypon_xy = [0 0; 1 0; 1 1; 0 1]' # 2xN array
+polypon_xy = [0 0; 1 0; 1 1; 0 1] # Nx2 array
 poly = get_polygon(polygon_xy, ratio=1) # poly.polygon to visualize
 pip_query(poly, 0, 0) # check points (0, 0) and return false
 pip_query(poly, 0, 0, edge=true) # check points (0, 0) and return true
@@ -48,7 +48,7 @@ the points are on the edge of the polygon as well.
 Example:
 ---
 ```julia
-polygon_xy = [0 0; 1 0; 1 1; 0 1]' # 2xN array
+polygon_xy = [0 0; 1 0; 1 1; 0 1] # Nx2 array
 poly = get_polygon(polygon_xy, ratio=1) # poly.polygon to visualize
 pip_query(polygon, [0 0; 0.5 0.5; 1 1], edge=true) # returns [true, true, true]
 ```
@@ -60,9 +60,9 @@ function pip_query(
 )
     # inputs check for points
     n, m = size(points)
-    m ≥ 1 || error("points should have at least 1 point")
-    n == 2 || error("points should be a 2xN array")
-    py_points = shapely.points(points')
+    n ≥ 1 || error("points should have at least 1 point")
+    m == 2 || error("points should be a Nx2 array")
+    py_points = shapely.points(points)
     
     # check if the particle is inside the polygon
     println("\e[1;36m[start]:\e[0m querying points inside the polygon")
@@ -82,14 +82,14 @@ Description:
 ---
 Determine whether a set of points is inside a 3D mesh defined by a 'STLInfo'. The input 
 `stl_model` should be initiated by function `readSTL2D(stl_file)`, and `points` should be 
-a 2xN array where each column represents a point (x, y). If `edge` is set to true, it checks 
+a Nx2 array where each column represents a point (x, y). If `edge` is set to true, it checks 
 if the points are on the edge of the mesh as well.
 
 Example:
 ---
 ```julia
 stl_model = readSTL2D("path/to/your/file.stl")
-points = [0 0; 0.5 0.5; 1 1]'
+points = [0 0; 0.5 0.5; 1 1]
 pip_query(stl_model, points, edge=true)
 ```
 """
@@ -98,8 +98,8 @@ function pip_query(
     points   ::AbstractMatrix;
     edge     ::Bool=false
 )
-    size(points, 1) == 2 || error("points must be a 2xN array")
-    py_points = shapely.points(points')
+    size(points, 2) == 2 || error("points must be a Nx2 array")
+    py_points = shapely.points(points)
     triangle_coords_2d = stl_model.py_vertices[stl_model.py_triangles]
     tris2d = shapely.polygons(triangle_coords_2d)
     region = shapely.unary_union(tris2d)
diff --git a/src/polyhedron.jl b/src/polyhedron.jl
@@ -14,10 +14,10 @@ function pip_query(
 )
     # inputs check for points
     n, m = size(points)
-    n == 3 || error("points must be a 3xN array")
-    m >= 1 || error("at least 1 points are required")
+    m == 3 || error("points must be a Nx3 array")
+    n >= 1 || error("at least 1 points are required")
     (nsamples > 0 && isodd(nsamples)) || error("nsamples must be a positive odd integer")
-    py_points = np.array(points')
+    py_points = np.array(points)
 
     # query by using open3d
     vertices = stl_model.py_vertices
diff --git a/src/utils.jl b/src/utils.jl
