[vtkusers] Median of a ROI

[Dean Inglis]

Hi Sergi, Michael,

see tcl script below on how to
do region of interest analysis with
vtkImageTracerWidget, vtkImageStencil and 
vtkImageAccumulate. It is not too difficult
to extrapolate this to using vtkSplineWidget.

Dean



package require vtk
package require vtkinteraction

  vtkVolume16Reader v16
  v16 SetDataDimensions 64 64
  v16 SetDataByteOrderToLittleEndian
  v16 SetImageRange 1 93
  v16 SetDataSpacing 3.2 3.2 1.5
  v16 SetFilePrefix "$VTK_DATA_ROOT/Data/headsq/quarter"
  v16 ReleaseDataFlagOn
  v16 SetDataMask 0x7fff

  vtkExtractVOI extract
  extract SetVOI 31 31 0 63 0 92
  extract SetSampleRate 1 1 1
  extract SetInput [v16 GetOutput]
  extract ReleaseDataFlagOff
  extract Update

  set range [[extract GetOutput] GetScalarRange]
  set min [lindex $range 0]
  set max [lindex $range 1]
  set diff [expr $max - $min]
  set avg [expr 0.5*($max + $min)]
  set bounds [[extract GetOutput] GetBounds]

  vtkPoints points
  points SetNumberOfPoints 5

  set x [lindex $bounds 0]
  set y0 [expr 0.25*[lindex $bounds 3]]
  set y1 [expr 0.75*[lindex $bounds 3]]
  set z0 [expr 0.25*[lindex $bounds 5]]
  set z1 [expr 0.75*[lindex $bounds 5]]

  points SetPoint 0 $x $y0 $z0
  points SetPoint 1 $x $y1 $z0
  points SetPoint 2 $x $y1 $z1
  points SetPoint 3 $x $y0 $z1
  points SetPoint 4 $x $y0 $z0

  vtkCellArray cells
  cells InsertNextCell 5

  for {set i 0} {$i < 5} {incr i} {
    cells InsertCellPoint $i
    }

  vtkPolyData pathPoly
  pathPoly SetPoints points
  pathPoly SetLines cells

  vtkLinearExtrusionFilter extrude
  extrude SetScaleFactor 1
  extrude SetExtrusionTypeToNormalExtrusion
  extrude SetVector 1 0 0
  extrude SetInput pathPoly

  vtkPolyDataToImageStencil dataToStencil
  dataToStencil SetInput [extrude GetOutput]

  vtkImageStencil stencil
  stencil ReleaseDataFlagOff
  stencil SetInput [extract GetOutput]
  stencil SetStencil [dataToStencil GetOutput]
  stencil ReverseStencilOff
  stencil SetBackgroundValue $min

  vtkImageMapToWindowLevelColors imageMapper
  imageMapper SetOutputFormatToLuminance
  imageMapper SetInput [stencil GetOutput]
  imageMapper SetLevel $avg
  imageMapper SetWindow $diff

  vtkImageActor imageActor
  imageActor SetInput [imageMapper GetOutput]
  imageActor SetDisplayExtent 31 31 0 63 0 92
  imageActor InterpolateOff

  vtkImageAccumulate accum
  accum SetInput [extract GetOutput]
  accum ReverseStencilOff
  accum SetStencil [stencil GetStencil]
  accum SetComponentSpacing 1 0 0
  accum SetComponentExtent 0 $diff 0 0 0 0
  accum SetComponentOrigin $min 0 0
  accum ReleaseDataFlagOff


  vtkTextActor textActor
  textActor ScaledTextOff

  textActor SetInput "testing..."

  [textActor GetPositionCoordinate] SetCoordinateSystemToNormalizedDisplay
  [textActor GetPositionCoordinate] SetValue 0.01 0.99

  set tprop [textActor GetTextProperty]
  $tprop SetFontSize 16
  $tprop SetFontFamilyToArial
  $tprop BoldOff
  $tprop ItalicOff
  $tprop ShadowOff
  $tprop SetLineSpacing 0.9
  $tprop SetJustificationToLeft
  $tprop SetVerticalJustificationToTop
  $tprop SetColor 1 1 0

  vtkRenderer ren
  ren SetBackground 0.4 0.4 0.5

  vtkRenderWindow renWin
  renWin AddRenderer ren
  renWin SetSize 800 800

  vtkInteractorStyleImage interactor

  vtkRenderWindowInteractor iren
  iren SetInteractorStyle interactor
  iren SetRenderWindow renWin

# Set up the image tracer widget
#
vtkImageTracerWidget itw
#
# Set the tolerance for capturing last handle when near first handle
# to form closed paths.
#
  itw SetCaptureRadius 1.5
  [itw GetGlyphSource] SetColor 1 0 0
#
# Set the size of the glyph handle
#
  [itw GetGlyphSource] SetScale 3.0
#
# Set the initial rotation of the glyph if desired.  The default glyph
# set internally by the widget is a '+' so rotating 45 deg. gives a 'x'
#
  set bounds [imageActor GetBounds]
  set pos [lindex $bounds 0]

  [itw GetGlyphSource] SetRotationAngle 45.0
  [itw GetGlyphSource] Modified
  itw ProjectToPlaneOn
  itw SetProjectionNormalToXAxes
  itw SetProjectionPosition $pos
  itw SetViewProp imageActor
  itw SetInput [imageMapper GetOutput]
  itw SetInteractor iren
  itw PlaceWidget
#
# When the underlying vtkDataSet is a vtkImageData, the widget can be
# forced to snap to either nearest pixel points, or pixel centers.  Here
# it is turned off.
#
  itw SnapToImageOff
#
# Automatically form closed paths.
#
  itw AutoCloseOn
  itw AddObserver EndInteractionEvent AdjustROI
  itw On
  itw InitializeHandles points

  ren AddViewProp imageActor  

  iren AddObserver UserEvent {wm deiconify .vtkInteract}
  renWin Render

  set cam [ren GetActiveCamera]
  $cam SetViewUp 0 1 0
  $cam Azimuth 270
  $cam Roll 270
  $cam Dolly 1.7
  ren ResetCameraClippingRange

  renWin Render

  ren AddViewProp textActor

  accum Update

  set aN [accum GetVoxelCount]
  set amin [accum GetMin]
  set amax [accum GetMax]
  set amean [accum GetMean]
  set astd  [accum GetStandardDeviation]

  set atext {}
  append atext "N: "
  append atext [string toupper $aN]
  append atext ", \["
  append atext [string toupper [lindex $amin 0] ]
  append atext " -> "
  append atext [string toupper [lindex $amax 0] ]
  append atext "\] "
  append atext [string toupper [lindex $amean 0] ]
  append atext " +- "
  append atext [string toupper [lindex $astd 0] ]
  append atext " HU"

  textActor SetInput $atext

  renWin Render


# Prevent the tk window from showing up then start the event loop.
  wm withdraw .


proc AdjustROI { } {

  set npts [ itw GetNumberOfHandles ]

  set closed [itw IsClosed]
  
  if { $closed } {
    itw GetPath pathPoly 
    extrude SetInput pathPoly 
    imageMapper SetInput [stencil GetOutput]

  stencil UpdateWholeExtent
  accum UpdateWholeExtent


  set aN [accum GetVoxelCount]
  set amin [accum GetMin]
  set amax [accum GetMax]
  set amean [accum GetMean]
  set astd  [accum GetStandardDeviation]

  set atext {}
  append atext "N: "
  append atext [string toupper $aN]
  append atext ", \["
  append atext [string toupper [lindex $amin 0] ]
  append atext " -> "
  append atext [string toupper [lindex $amax 0] ]
  append atext "\] "
  append atext [string toupper [lindex $amean 0] ]
  append atext " +- "
  append atext [string toupper [lindex $astd 0] ]
  append atext " HU"

  textActor SetInput $atext

  textActor Modified


  } else {
    textActor SetInput "no ROI"
    imageMapper SetInput [extract GetOutput]
  }
}