[ITK-users] watershed on surface meshs

[Grothausmann, Roman Dr.]

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. [phone]
***@mh-hannover.de
http://www.mh-hannover.de/anatomie.html

[Dženan Zukić]

Hi Roman,

I don't know of any implementation of mesh watershed segmentation. ITK does
not really have mesh processing capabilities, it is mostly limited to
conversion between mesh and image representations.

A colleague <http://www.cg.informatik.uni-siegen.de/en/chiosa-iurie> worked
on something related, maybe he can provide the source code if you ask.

And we are moving to discourse <https://discourse.itk.org/>. Re-posting
this question there might be useful.

Regards,
DÅŸenan

On Fri, Nov 24, 2017 at 6:26 AM, Grothausmann, Roman Dr. <

Post by Grothausmann, Roman Dr.
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
http://teacher.en.rmutt.ac.th/ktw/Resources/Full%20paper%20P
DF/Partitioning%203D%20surface%20meshes%20using%20watershed%
20segmentation.pdf
https://www.google.de/url?sa=t&rct=j&q=&esrc=s&source=web&cd
=4&cad=rja&uact=8&ved=0ahUKEwjD0by1lafWAhVUGsAKHZ2MAbUQFgg_M
AM&url=http%3A%2F%2Fciteseerx.ist.psu.edu%2Fviewdoc%2Fdownlo
ad%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. [phone]
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
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Richard Beare]

I don't recall an implementation of this anywhere. However you may be able
to hack something together reasonably quickly. Watersheds are built around
minimal paths, so depending on what complexity is available in the
mesh-based shortest path tools. In the worst cast you could create a copy
of the mesh for each seed region, compute the minimal distance to each non
seed vertex, then do a vertex-wise min, tracking which mesh/seed region has
the minimum, and that's your label. Not very efficient because you compute
the entire mesh distance each time, and no obvious way to do a watershed
line consistently, but perhaps that doesn't matter for testing the idea.

On Fri, Nov 24, 2017 at 10:26 PM, Grothausmann, Roman Dr. <

Post by Grothausmann, Roman Dr.
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
http://teacher.en.rmutt.ac.th/ktw/Resources/Full%20paper%20P
DF/Partitioning%203D%20surface%20meshes%20using%20watershed%
20segmentation.pdf
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. [phone]
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
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Richard Beare]

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. <

Post by Grothausmann, Roman Dr.
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
http://teacher.en.rmutt.ac.th/ktw/Resources/Full%20paper%20P
DF/Partitioning%203D%20surface%20meshes%20using%20watershed%
20segmentation.pdf
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. [phone]
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
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Matt McCormick]

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

Post by Richard Beare
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.

Post by Grothausmann, Roman Dr.
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
http://teacher.en.rmutt.ac.th/ktw/Resources/Full%20paper%20PDF/Partitioning%203D%20surface%20meshes%20using%20watershed%20segmentation.pdf
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. [phone]
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
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
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
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
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/list