[ITK-users] watershed on surface meshs
Richard Beare
richard.beare at gmail.com
Fri Nov 24 22:16:42 EST 2017
Just remembered something that may be useful for testing. I wrote the
following for doing watersheds on geospatial graphs using R/igraph.
Certainly isn't efficient, but may be useful for testing if you can turn
your mesh+curvature into something like a text file.
This was set up for allocating points on maps to destinations based on
travel time. It performed a kind of smoothing to smooth non-sensical travel
time estimates, but the basic algorithm is there. It will take some
fiddling. You could probably do something equivalent with python, probably
doing better than this priorty queue, but using much the same approach with
igraph.
library(liqueueR)
## Mostly a copy of PriorityQueue from liqueueR
StablePriorityQueue <- setRefClass("StablePriorityQueue",
contains = "Queue",
fields = list(
count = "numeric",
entries = "numeric",
priorities = "numeric",
prioritise = "function"
),
methods = list(
sort_ = function() {
order = order(priorities, entries,
decreasing = TRUE, partial = size():1)
#
data <<- data[order]
priorities <<- priorities[order]
entries <<- entries[order]
},
push = function(item, priority = NULL)
{
'Inserts element into the queue,
reordering according to priority.'
callSuper(item)
#
if (is.null(priority)) priority =
prioritise(item)
#
priorities <<- c(priorities,
priority)
entries <<- c(entries, count)
count <<- count - 1
#
sort_()
},
pop = function(N = 1) {
# 'Removes and returns head of queue
(or raises error if queue is empty).'
if (N > size()) stop("insufficient
items in queue!")
priorities <<-
priorities[-seq_len(N)]
entries <<- entries[-seq_len(N)]
callSuper(N)
},
initialize = function(prioritise =
NULL, ...) {
'Creates a PriorityQueue object.'
callSuper(...)
#
## to enforce FIFO order
count <<- 0
if (is.null(prioritise))
.self$prioritise = function(x) 0
else
.self$prioritise = prioritise
#
.self
}
)
)
igraph.watershed <- function(Gr, labelfield, unlabelled, vertexid, alltimes)
{
Gres <- Gr
## Watershed, without marking boundary (Beucher)
## 1. find all marker nodes that have a background neighbour
lablist <- which(vertex_attr(Gr, labelfield) != unlabelled)
nlist <- ego(Gr, 1, lablist)
uu <- which(map_lgl(nlist, ~any(vertex_attr(Gr, labelfield,
.x)==unlabelled)))
## uu indexes the labelled vertexes
## need indexes into all vertexes
uu <- lablist[uu]
## create priority queue
qq <- StablePriorityQueue$new()
## Insert the boundary markers
kk <- lapply(uu, qq$push, priority=0)
all.labels <- get.vertex.attribute(Gr, labelfield)
all.ids <- get.vertex.attribute(Gr, vertexid)
alltimes <- subset(alltimes, from %in% all.ids, select=c("from",
"Hospital", "minutes"))
dd <- duplicated(alltimes[, c("from", "Hospital")])
alltimes <- alltimes[!dd,]
alltimes.wide <- spread_(alltimes, key="Hospital", value="minutes")
## Make the order the same as all.ids - so now we'll be able to index by
number
rownames(alltimes.wide) <- alltimes.wide$from
alltimes.wide <- alltimes.wide[all.ids, ]
cc <- 1:ncol(alltimes.wide)
names(cc) <- colnames(alltimes.wide)
while (qq$size() > 0) {
vid <- qq$pop()
## Get the neighbours
nb <- neighborhood(Gr, 1, vid, mode="all", mindist=1)[[1]]
nlabs <- all.labels[nb]
## Are any neighbours unknown?
ul <- nlabs==unlabelled
if (any(ul)) {
this.label <- all.labels[vid]
nbb <- nb[ul]
all.labels[nbb] <- this.label
this.label.idx <- cc[this.label]
priorities <- alltimes.wide[[this.label.idx]][nbb]
kk <- map2(nbb, priorities, ~qq$push(.x, priority = .y* -1))
}
}
return(data.frame(PlaceID=all.ids, Hospital=all.labels, stringsAsFactors
= FALSE))
}
On Fri, Nov 24, 2017 at 10:26 PM, Grothausmann, Roman Dr. <
grothausmann.roman at mh-hannover.de> wrote:
> Dear mailing list members,
>
>
> I need to separate a mesh at "curved corners" (see attached PNG, using the
> colored labels from a facet analysis do not suffice but go in the right
> direction). So my current thought is to run vtkCurvature to get a Gaussian
> curvature value per point/vertex and then try to separate regions of
> positive values around local maxima. Just thresholding the result of
> vtkCurvature does not fully separate each local max from neighboring ones,
> but to my understanding a surface watershed would. I found two publications
> by Mangan and Whitaker on this:
>
> Partitioning 3D surface meshes using watershed:
> http://teacher.en.rmutt.ac.th/ktw/Resources/Full%20paper%20P
> DF/Partitioning%203D%20surface%20meshes%20using%20watershed%
> 20segmentation.pdf
>
> Surface Segmentation Using Morphological Watersheds:
> https://www.google.de/url?sa=t&rct=j&q=&esrc=s&source=web&cd
> =4&cad=rja&uact=8&ved=0ahUKEwjD0by1lafWAhVUGsAKHZ2MAbUQFgg_
> MAM&url=http%3A%2F%2Fciteseerx.ist.psu.edu%2Fviewdoc%
> 2Fdownload%3Fdoi%3D10.1.1.464.2788%26rep%3Drep1%26type%
> 3Dpdf&usg=AFQjCNGX-p9-ElQFcpsUyBRO0pCjBKCmNg
>
> Does anybody know about an implementation for this in VTK/ITK or another
> open-source library? If not, would it be possible to transfer the ITK
> watershed implementation for voxel data to mesh data, e.g. to crate a
> VTKmorphWatershedFilter?
>
> Thanks for any help or hints.
> Roman
>
>
> --
> Dr. Roman Grothausmann
>
> Tomographie und Digitale Bildverarbeitung
> Tomography and Digital Image Analysis
>
> Medizinische Hochschule Hannover
> Institut für Funktionelle und Angewandte Anatomie
> OE 4120, Carl-Neuberg-Str. 1, 30625 Hannover, Deutschland
>
> Tel. +49 511 532-2900
> grothausmann.roman at mh-hannover.de
> http://www.mh-hannover.de/anatomie.html
>
> The ITK community is transitioning from this mailing list to
> discourse.itk.org. Please join us there!
> ________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/insight-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-users/attachments/20171125/cd6e3611/attachment.html>
More information about the Insight-users
mailing list