0

I have a georeferenced mesh produced by photogrammetry software such as Agisoft Photoscan. Within this mesh, I have some areas of interest, which I can define with 2D vector masks (on the x and y plane) in the form of a shapefile or a SpatialPolygonsDataFrame (R).

Can anyone reccomend a method to clip mesh files (typically in .obj format) using a mask?

This is like an extract-by-mask workflow, but on a mesh as opposed to a typical raster.

A typical mesh in this case looks something like this: enter image description here

Ideally I wish to use R or python.

  • Are these triangular mesh objects? Can you point us to a sample? – Spacedman Jun 29 '17 at 21:02
  • Yes they are triangular. Please find a small sample .obj and .mtl here - dropbox.com/s/3yxpkieaxon1qsu/samplemesh.zip?dl=0 – JPD Jun 30 '17 at 7:29
  • Is there a way to read these files into R or python? What output do you want - a mesh file of the same form but only with triangles that are fully or partially in your area of interest? – Spacedman Jun 30 '17 at 9:20
  • Regarding whether they can be read in, that's part of what I wanted to find out. The output would be triangles that fall within the area of interest, and ultimately the total area coverage of the triangles in the area of interest. – JPD Jun 30 '17 at 9:39
  • I've written a function to read the triangles. Not sure what the lines beginning "vt " are (6114 of them in your data). There's 2038 triangles from 1115 vertices, so there's 3 per triangle but they aren't the same coordinates as the vertex ("v ") lines... – Spacedman Jun 30 '17 at 10:48
2

Here's a morning's work:

library(sp)

readmesh <- function(file){
    all_lines = read.table(file, fill=TRUE, stringsAsFactors=FALSE)
    v = all_lines[all_lines$V1=="v",-1]
    names(v)=c("x","y","z")
    v$x = as.numeric(v$x)
    v$y = as.numeric(v$y)
    v$z = as.numeric(v$z)
    coordinates(v)=~x+y

    vt = all_lines[all_lines$V1=="vt",2:3]
    names(vt)=c("Z1","Z2")

    f =  all_lines[all_lines$V1=="f",-1]
    f = cbind(f, do.call(rbind,strsplit(f$V2,"/")))
    f = cbind(f, do.call(rbind,strsplit(f$V3,"/")))
    f = cbind(f, do.call(rbind,strsplit(f$V4,"/")))

    f = f[,c(4,6,8)]
    names(f) = c("V1","V2","V3")
    fn = function(x){as.numeric(as.character(x))}
    f$V1 = fn(f$V1)
    f$V2 = fn(f$V2)
    f$V3 = fn(f$V3)
    list(v=v, vt=vt, f=as.matrix(f))
}

plot_tri <- function(m, n, gf=lines, ...){
    pts = coordinates(m$v[m$f[n,],])
    gf(pts2tri(pts),...)
}

pts2tri <- function(pts3){
    pts3[c(1,2,3,1),]
}

tri_area <- function(p){
    0.5 * det(rbind(p[,"x"],p[,"y"],1))
}

tri_areas <- function(m, intri){
    tris = (1:nrow(m$f))[intri]
    sapply(tris,
           function(i){
               tri_area(coordinates(m$v[m$f[i,],]))
           }
           )
}

tris_in <- function(m, spoly){
    inout = which(!is.na(over(m$v,spoly)))
    intri = apply(m$f,1, function(r){ all(r %in% inout)})
    intri
}

tri_area3 <- function(x,y,z){
    x1 = x[1]; x2=x[2]; x3=x[3]
    y1 = y[1]; y2=y[2]; y3=y[3]
    z1 = z[1]; z2=z[2]; z3=z[3]
    0.5* sqrt(
        ((x2*y1) - (x3*y1) - (x1*y2) + (x3*y2) + (x1*y3) - (x2*y3))^2 +
        ((x2*z1) - (x3*z1) - (x1*z2) + (x3*z2) + (x1*z3) - (x2*z3))^2 +
        ((y2*z1) - (y3*z1) - (y1*z2) + (y3*z2) + (y1*z3) - (y2*z3))^2
        )
}

tri_areas3 <- function(m,intri){
    tris = (1:nrow(m$f))[intri]
    sapply(tris,
           function(i){
               pts = m$v[m$f[i,],]
               cpts = coordinates(pts)
               tri_area3(cpts[,1], cpts[,2], pts$z)
           }
           )
}

Usage - read in mesh object file and plot the vertex points:

> m = readmesh("./samplemesh.obj")
> plot(m$v)

Now click to create a polygon for cropping, convert to a spatial polygon object:

> p = locator(type="l") # mouse-2 to stop
> p = SpatialPolygons(list(Polygons(list(Polygon(p)),"P")))

Compute indexes of triangles within the polygon:

> intri = tris_in(m, p)

Draw those triangles:

> for(i in (1:nrow(m$f))[intri]){plot_tri(m,i,col="blue")}

Sum of those triangle areas, should be a bit less than the area of the polygon:

> sum(tri_areas(m,intri))
[1] 90.68329
> p@polygons[[1]]@area
[1] 121.5584

enter image description here

Edit

For 3d triangle areas, use the extra functions above and then you get (for a different polygon than the 2d example previously):

> sum(tri_areas(m,intri))
[1] 137.995
> sum(tri_areas3(m,intri))
[1] 142.4022
> p@polygons[[1]]@area
[1] 163.1959

So the 3d area is slightly larger - but that's because the height doesn't have as much of a range as the x and y coordinates and so doesn't have a great effect on the areas. Multiply the z heights by ten and it makes a great difference!

And now in THREE-D

If you want to visualise the mesh in 3d:

library(rgl)
mesh2tri <- function(m, ...){
    v = c(t(m$f))
    pts = m$v[v,]
    xyz = cbind(coordinates(pts),pts$z)
    triangles3d(xyz, ...)
}

then do:

mesh2tri(m,col="red")
axis3d("x"); axis3d("y");axis3d("z")

gives:

enter image description here

Invoice is in the post :)

  • Wow - that's some great work. Thanks very much. One question - given that the mesh is in 3D space, the sum area of triangles is likely to be greater than the area of the 2D mask. I'm slowly working through your code to understand it...does it calculate the areas from the triangles in 3D space? – JPD Jun 30 '17 at 12:48
  • No, its just the 2d area as projected onto the x-y coordinates. I'd have to look up the formula for the surface area of a triangle from points in 3d, which I guess is a projection to a plane and then the 2d formula in projected coordinates. You need to make sure your x, y, and z coords are all the same units as well. – Spacedman Jun 30 '17 at 13:11
  • Yep, the mesh is projected (EPSG:32643 in this case), and the clipping shapefile can be in that projection too. – JPD Jun 30 '17 at 13:31
  • Found the formula for triangle area in 3d, tested it, stuck it into my code as an edit. – Spacedman Jun 30 '17 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.