[ITK-users] watershed on surface meshs
Matt McCormick
matt.mccormick at kitware.com
Sat Dec 2 10:11:30 EST 2017
Hi Roman,
Another possible idea is to use ITK's mesh fast marching capabilities:
https://itk.org/ITKExamples/src/Filtering/FastMarching/ComputeGeodesicDistanceOnMesh/Documentation.html
with curvature as a speed function in same optimization approach used
to compute minimal paths on images:
http://www.insight-journal.org/browse/publication/213
HTH,
Matt
On Fri, Nov 24, 2017 at 10:16 PM, Richard Beare <richard.beare at gmail.com> wrote:
> 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%20PDF/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
>>
>
>
> 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
>
More information about the Insight-users
mailing list