Gnome Crystal

Gnome Crystal(gcrystal)というソフトを使ってみました。今ではthe Gnome Chemistry Utilsの一部になっているそうです。Debian etchでパッケージ提供されないかな…
◆ビルドの記録

  336  apt-get install libglade
  337  apt-get install libglade2-0
  338  apt-get install libglade2-0-dev
  339  apt-get install libglade
  340  apt-cache search libglade
  341  apt-get install libglade2-dev
  342  apt-cache search libgnomeprint
  343  apt-get install libgnomeprintui2.2-0
  344  apt-get install libgnomeprintui2.2-dev
  345  apt-get install libgnomeui
  346  apt-cache search libgnomeui
  347  apt-cache search libgnomeui-dev
  348  apt-get install libgnomeui-dev
  349  apt-cache search gcu
  350  apt-get install libgcu0
  351  apt-get install libgcu-dev
  352  apt-get install libgcu-doc
  353  apt-cache search libpng
  354  apt-get install libpng12-dev
libgnome-devもかな？

$ ./configure --prefix=/home/bronyraur/usr
$ make
$ make install

◆ダイヤモンド構造の結晶の画像

FreeFem++ ex02.edp

//
// １次元の微分方程式（単振動）
// d^2u/dx^2 + u = 0
//
// xを時間、yは一様、uを位置とみたてて問題を解く
//

mesh Th = square(100,1,[0+2*pi*x, y]);  // 100x1 elements

plot(Th, wait=true, ps="ex02.mesh.eps");

//fespace Vh(Th, P1);    //X軸に対して値が対称にならない？
fespace Vh(Th, P2);
Vh u, v;

problem ex02(u, v)
  = int2d(Th)( -dx(u)*dx(v) ) + int2d(Th)( u*v )
    + on(4, u=1);  //x=0でu=1、du/dx=0,  x=2*piでdu/dx=0
//    + int1d(Th,1)(0*v)  // y=-1の境界条件 du/dy = 0 (指定しなくてもOK)
//    - int1d(Th,3)(0*v); // y=-1の境界条件 du/dy = 0 (指定しなくてもOK)
ex02;

plot(u, wait=1, fill=1, value=1, ps="ex02.value.eps");

// gnuplot用に結果を出力する
{
  ofstream ff("ex02.gnuplot");
  for(int i=0; i<Th.nt; i++){
    ff << Th[i][0].x << " " << Th[i][0].y << " " << u[][Vh(i,0)] << endl;
  }
}

FreeFem++ ex01.edp

メッシュ

計算結果

計算結果をgnuplotで表示

// ex01.edp
//
// １次元の微分方程式（物体の落下）
// d^2u/dx^2 + 9.8 = 0
//
// xを時間、yは一様、uを位置とみたてて問題を解く
//
mesh Th = square(100,1,[0+10*x, y]);  // 100x1 elements 0<=x<=10, 0<=y<=1

plot(Th, wait=true, ps="ex01.mesh.eps");

//fespace Vh(Th, P1);  //X軸に対して値が対称にならない？
fespace Vh(Th, P2);
Vh u, vv;

real v = 30;
real g = -9.8;
problem ex01(u, vv)
  = int2d(Th)(-dx(u)*dx(vv)) + int2d(Th)(-g*vv)
    +on(4, u=0) - int1d(Th,4)(v*vv)  //初速30, x=0でu=0
    + int1d(Th,2)((g*x+v)*vv);       //x=10での速さを指定
ex01;

plot(u, wait=1, fill=1, value=1, ps="ex01.value.eps");

// gnuplot用に結果を出力する
{
  ofstream ff("ex01.gnuplot");
  for(int i=0; i<Th.nt; i++){
    ff << Th[i][0].x << " " << Th[i][0].y << " " << u[][Vh(i,0)] << endl;
  }
}

vtk - ascviewer.py [01]

vtkを使ってDICOMファイルを表示するpythonスクリプトascviewer.pyを作成しました。
画像はこれを使って表示したものです。

使い方は、

'o': DICOMディレクトリを指定してロードする。
'v': visualizationのパラメータを設定する。
'u': 表示している画像を'image.png'というファイルで保存する。
SPACE, SHIFT+SPACE: Axial, Saggital, Coronalビューポートではスライスを増減する。3Dビューポートでは、ExtractVOIのSampleRateを増減する。

※3Dビューポートでは、通常のInteractorの操作ができます。

#! /usr/bin/env python
# -*- coding:utf-8 -*-
#
# ascviewer.py
#
# ボリュームデータから、それを横切る平面の画像とisosurfaceを作成する
# $Date: 2007-07-02 11:42:51 $
#

import vtk
import Tkinter
import tkFileDialog


########################################################################
#ビューポートにフィットさせるための拡大率を計算する
########################################################################
def get_magnification_factor(renderer):
  global renderWindow, renderer_A, renderer_S, renderer_C
  global spacing, dimensions
  #print "get_magnification_factor"
  if renderer is renderer_A:
    image_w = dimensions[0]
    image_h = dimensions[1]
  elif renderer is renderer_S:
    image_w = dimensions[2]
    image_h = dimensions[1]
  elif renderer is renderer_C:
    image_w = dimensions[0]
    image_h = dimensions[2]
  else: return [0.0, 0.0]
#  print "image_w, image_h =", image_w, image_h
  renderer_viewport_rect = renderer.GetViewport()
  window_w, window_h = renderWindow.GetSize()
  renderer_viewport_w = renderer_viewport_rect[2] - renderer_viewport_rect [0]
  renderer_viewport_h = renderer_viewport_rect[3] - renderer_viewport_rect [1]
  renderer_w = window_w * renderer_viewport_w
  renderer_h = window_h * renderer_viewport_h
  magnification_factor = [float(renderer_w)/image_w, float(renderer_h)/image_h]
#  print "renderer_viewport_rect =", renderer_viewport_rect
#  print "window_w, window_h, renderer_w, renderer_h =", \
#         window_w, window_h, renderer_w, renderer_h
#  print "magnification_factor =", magnification_factor
  return magnification_factor
########################################################################

########################################################################
#DICOMファイルをロードする
########################################################################
def open_dicom():
  # ディレクトリの指定
  root = Tkinter.Tk()  #空Tkウィンドウのポップアップ表示を避ける
  root.withdraw()
  dirname = tkFileDialog.askdirectory()
  root.destroy()
#  print 'get_dicom_dirname(), dirname =', dirname
  if dirname == ():
    return
  # dicom source objectの作成
  global sourceProcessObject, renderWindow
  global dimensions, spacing, slice_index
  global color_window, color_level, contour_level
  global imageReslice_S
  dicomImageReader = vtk.vtkDICOMImageReader()
  dicomImageReader.SetDirectoryName(dirname)
  dicomImageReader.Update()
  spacing = dicomImageReader.GetDataSpacing()
  dimensions = dicomImageReader.GetOutput().GetDimensions()
  slice_index = [int(dimension/2) for dimension in dimensions]
  color_window  = 400
  color_level   = 128
  contour_level = 0.0
  sourceProcessObject = dicomImageReader
  print "dimensions =", dimensions
  print "spacing =", spacing
  initialize_visualization_pipeline()
  renderWindow.Render()
########################################################################

########################################################################
#可視化のパラメータを設定する
########################################################################
def configure_visualization_parameters():
  def apply():
    global color_level, color_window
    print 'apply'
    color_window = scale_color_window.get()
    color_level  = scale_color_level.get()
    print 'color_window =',
    print color_window
    print 'color_level =',
    print color_level
    global contour_level, voi_sample_rate
    contour_level = scale_contour_level.get()
    voi_sample_rate = scale_voi_sample_rate.get()
    print 'contour_level =',
    print contour_level
    print 'voi_sample_rate =',
    print voi_sample_rate
    tk.destroy()
    update_color()
    update_3D()
  tk = Tkinter.Tk()
  tk.title("Configure Visualization Parameters")
  #
  scale_color_window = Tkinter.Scale(tk, label='Color Window', length=300,
                         orient='horizontal', from_=-1024, to=4096)
  scale_color_window.pack()
  scale_color_window.set(color_window)
  scale_color_level = Tkinter.Scale(tk, label='Color Level', length=300,
                         orient='horizontal', from_=-1024, to=4096)
  scale_color_level.pack()
  scale_color_level.set(color_level)
  #
  separator = Tkinter.Frame(tk, width=300, height=2,
                borderwidth=2, relief='sunken')
  separator.pack(pady=10)
  #
  scale_voi_sample_rate = Tkinter.Scale(tk, label='ExtractVOI SampleRate',
                            length=300, orient='horizontal', from_=1, to=16)
  scale_voi_sample_rate.set(voi_sample_rate)
  scale_voi_sample_rate.pack()
  scale_contour_level = Tkinter.Scale(tk, label='Contour Level', length=300,
                          orient='horizontal', from_=-1024, to=4096)
  scale_contour_level.set(contour_level)
  scale_contour_level.pack()
  #
  separator = Tkinter.Frame(tk, width=300, height=2,
                borderwidth=2, relief='sunken')
  separator.pack(pady=10)
  #
  button_apply = Tkinter.Button(tk, text='Apply', command=apply)
  button_apply.pack(side='right')
  tk.mainloop()
########################################################################

########################################################################
#color_window, color_levelを設定してrenderWindowを更新する
########################################################################
def update_color():
  global color_window, color_level
  global imageMapper_A, imageMapper_S, imageMapper_C
  global renderWindow
  for im in [imageMapper_A, imageMapper_S, imageMapper_C]:
    im.SetColorWindow(color_window)
    im.SetColorLevel(color_level)
  renderWindow.Render()
########################################################################

########################################################################
#contour_level, を設定してrenderWindowを更新する
########################################################################
def update_3D():
  global contour_level, voi_sample_rate
  global extractVOI_3, contourFilter_3, renderWindow
  contourFilter_3.SetValue(0, contour_level)
  extractVOI_3.SetSampleRate([voi_sample_rate for i in range(3)])
#  ,voi_sample_rate,voi_sample_rate)
  renderWindow.Render()
########################################################################

########################################################################
# この関数では、renderWindowで"uキー"が入力されたときの処理を記述します。
#
# この例では、renderWindowに表示されている画像を"image.png"という名前の
# ファイルに保存する処理を記述しています。
# この関数では、renderWindowで"uキー"が入力されたときの処理を記述します。
#
# この例では、renderWindowに表示されている画像を"image.png"という名前の
# ファイルに保存する処理を記述しています。
########################################################################
def cb_user_event(obj, event):
  print "cb_user_event"
  windowToImageFilter=vtk.vtkWindowToImageFilter()
  windowToImageFilter.SetInput(renderWindow)
  windowToImageFilter.Update()
  pngWriter=vtk.vtkPNGWriter()
  pngWriter.SetInput(windowToImageFilter.GetOutput())
  pngWriter.SetFileName("image.png")
  renderWindow.Render()
  pngWriter.Write()
########################################################################

########################################################################
# この関数では、renderWindowでキーが押されたときの処理を記述します。
#
########################################################################
def cb_key_press_event(obj, event):
  global renderer_3, renderer_A, renderer_S, renderer_C
  global imageReslice_A, imageReslice_S, imageReslice_C
  global dimensions, spacing, slice_index, contour_level
  #print "cb_key_press_event"
  #print type(obj), type(event)
  #print obj
  #print obj.GetClassName()
  pos = obj.GetEventPosition()
  renderer = obj.FindPokedRenderer(pos[0], pos[1])
  if renderer not in [renderer_3, renderer_A, renderer_S, renderer_C]:
    return

  #print "%d" % ord(obj.GetKeyCode())
  #print chr(ord(obj.GetKeyCode()))
  if obj.GetKeyCode() not in ' ov':
    #print "key code is not space"
    return

  if obj.GetKeyCode() == 'o':
    print "file open"
    open_dicom()
    return

  if obj.GetKeyCode() == 'v':
    print "configure_visualization_parameters"
    configure_visualization_parameters()
    return

  if renderer is renderer_3:
    print "renderer_3"
    d = [1, -1][obj.GetShiftKey()]
    m = [10, 1][obj.GetControlKey()]  # Controlキーを認識してくれない
    #print d, m, obj.GetControlKey()
    contour_level = contour_level + m*d
    contourFilter_3.SetValue(0, contour_level)
    #contourFilter3.Update()
  elif renderer is renderer_A:
    #print "renderer_A"
    slice_index[2] \
      = max(
          min(slice_index[2]+[1,-1][obj.GetShiftKey()], dimensions[2]-1),
          0)
    imageReslice_A.SetResliceAxesOrigin(
                     0,
             0,
             slice_index[2]*spacing[2])
  elif renderer is renderer_S:
    #print "renderer_S"
    slice_index[0] \
      = max(
          min(slice_index[0] + [1,-1][obj.GetShiftKey()], dimensions[0] - 1),
      0)
    imageReslice_S.SetResliceAxesOrigin(
                      slice_index[0]*spacing[0],
              0,
              0)
  elif renderer is renderer_C:
    #print "renderer_C"
    slice_index[1] \
      =max(
         min(slice_index[1] + [1,-1][obj.GetShiftKey()], dimensions[1] - 1),
         0)
    imageReslice_C.SetResliceAxesOrigin(
                      0,
                      slice_index[1]*spacing[1],
                      (dimensions[2]-1)*spacing[2])

  print "slice_index =", slice_index
  renderWindow.Render()
########################################################################

########################################################################
# この関数は、ModifiedEventが発生したときに呼び出されます。
# ここでは、Window のリサイズの処理に対応させています。
########################################################################
def cb_modified_event(obj, event):
  global imageResample_A, imageResample_S, imageResample_C
  global renderer_A, renderer_S, renderer_C
  if None in [renderer_A, renderer_S, renderer_C]:
    return
  #print obj
  #print event
  mf = get_magnification_factor(renderer_A)
  imageResample_A.SetAxisMagnificationFactor(0, mf[0])
  imageResample_A.SetAxisMagnificationFactor(1, mf[1])
  mf = get_magnification_factor(renderer_S)
  imageResample_S.SetAxisMagnificationFactor(0, mf[0])
  imageResample_S.SetAxisMagnificationFactor(1, mf[1])
  mf = get_magnification_factor(renderer_C)
  imageResample_C.SetAxisMagnificationFactor(0, mf[0])
  imageResample_C.SetAxisMagnificationFactor(1, mf[1])
  #renderWindow.Render()
########################################################################

########################################################################
# "Visualization Pipeline"を作成し、actorにmapperを設定します。
########################################################################
def initialize_visualization_pipeline():
  global vtk
  global processSourceObject
  global renderer_3, rednerer_A, renderer_S, renderer_C
  global extractVOI_3, contourFilter_3
  global imageReslice_A, imageResample_A, imageMapper_A
  global imageReslice_S, imageResample_S, imageMapper_S
  global imageReslice_C, imageResample_C, imageMapper_C
  global dimensions, spacing
  global slice_index
  global color_window, color_level, voi_sample_rate

  #3D (extract skin)
  extractVOI_3 = vtk.vtkExtractVOI()
  extractVOI_3.SetInput(sourceProcessObject.GetOutput())
  extractVOI_3.SetVOI(0,dimensions[0]-1, 0,dimensions[1]-1, 0,dimensions[2]-1)
  extractVOI_3.SetSampleRate([voi_sample_rate for i in range(3)])
  print extractVOI_3.GetSampleRate()
  contourFilter_3 = vtk.vtkContourFilter()
  contourFilter_3.SetInput(extractVOI_3.GetOutput())
  contourFilter_3.SetValue(0,contour_level)
  contourFilter_3.ComputeNormalsOff()
  contourFilter_3.ComputeGradientsOff()
  polyDataMapper_3=vtk.vtkPolyDataMapper()
  polyDataMapper_3.SetInput(contourFilter_3.GetOutput())
  polyDataMapper_3.ScalarVisibilityOff()
  actor_3=vtk.vtkActor()
  actor_3.SetMapper(polyDataMapper_3)
  # アクターをすべて取り除く
  actorCollections = renderer_3.GetActors()
#  print '3:actorCollections.GetNumberOfItems() =',
#  print actorCollections.GetNumberOfItems()
  for i in range(actorCollections.GetNumberOfItems()):
    renderer_3.RemoveActor(actorCollections.GetItemAsObject(i))
  actorCollections = renderer_3.GetActors()
#  print '3:actorCollections.GetNumberOfItems() =',
#  print actorCollections.GetNumberOfItems()
  #
  renderer_3.AddActor(actor_3)
  renderer_3.ResetCamera()

  #Axial
  #  vtkImageReslice -> vtkImageResample -> vtkImageMapper -> vtkActor2D
  imageReslice_A = vtk.vtkImageReslice()
  imageReslice_A.SetInput(sourceProcessObject.GetOutput())
  imageReslice_A.SetOutputDimensionality(2)
  imageReslice_A.SetResliceAxesOrigin(0,0,slice_index[2]*spacing[2])
  imageReslice_A.SetResliceAxesDirectionCosines(
                   spacing[0],0,0,
                   0,spacing[1],0,
                   0,0,spacing[2])
  imageReslice_A.InterpolateOn()
  imageReslice_A.SetOutputSpacing(1,1,1)   # spacingを強制的に1にする
  #
  mf = get_magnification_factor(renderer_A)
  imageResample_A = vtk.vtkImageResample()
  imageResample_A.SetInput(imageReslice_A.GetOutput())
  imageResample_A.SetDimensionality(2)
  imageResample_A.InterpolateOn()
  imageResample_A.SetInterpolationModeToCubic()
  imageResample_A.SetAxisMagnificationFactor(0, mf[0])
  imageResample_A.SetAxisMagnificationFactor(1, mf[1])
  #
  imageMapper_A = vtk.vtkImageMapper()
  imageMapper_A.SetInput(imageResample_A.GetOutput())
  imageMapper_A.SetColorWindow(color_window)
  imageMapper_A.SetColorLevel(color_level)
  # アクターをすべて取り除く
  actor2DCollections = renderer_A.GetActors2D()
  for i in range(actor2DCollections.GetNumberOfItems()):
    renderer_A.RemoveActor(actor2DCollections.GetItemAsObject(i))
  actor2DCollections = renderer_A.GetActors2D()
  #
  actor2D_A = vtk.vtkActor2D()
  actor2D_A.SetMapper(imageMapper_A)
  renderer_A.AddActor(actor2D_A)
  #
  textActor_A = vtk.vtkTextActor()
  textActor_A.ScaledTextOff()
  textActor_A.SetPosition(0,0)
  textActor_A.SetInput("Axial")
  textActor_A.GetTextProperty().SetFontSize(12)
  textActor_A.GetTextProperty().SetFontFamilyToArial()
  textActor_A.GetTextProperty().BoldOn()
  textActor_A.GetTextProperty().ShadowOn()
  textActor_A.GetTextProperty().SetColor(1,1,0)
  renderer_A.AddActor2D(textActor_A)


  #Saggital
  #  vtkImageReslice -> vtkImageResample -> vtkImageMapper -> vtkActor2D
  #
  # 断面画像の切り出し
  # ※2D画像の実験はココ(Saggital画像)で行う
  imageReslice_S = vtk.vtkImageReslice()
  imageReslice_S.SetInput(sourceProcessObject.GetOutput())
  imageReslice_S.SetOutputDimensionality(2)
  imageReslice_S.SetResliceAxesOrigin(slice_index[0]*spacing[0], 0, 0)
  imageReslice_S.SetResliceAxesDirectionCosines(
  #                 0,0,1,
  #                 0,1,0,
  #                 1,0,0)
                    0,0,spacing[2],
                    0,spacing[1],0,
                    spacing[0],0,0)
  imageReslice_S.InterpolateOn()
  imageReslice_S.SetOutputSpacing(1,1,1)  # spacingを強制的に1にする

  #
  # 断面画像の大きさをViewportの大きさにリサイズする
  mf = get_magnification_factor(renderer_S)
  #print "renderer_S: mf =", mf
  imageResample_S = vtk.vtkImageResample()
  imageResample_S.SetInput(imageReslice_S.GetOutput())
  imageResample_S.SetDimensionality(2)
  imageResample_S.InterpolateOn()
  imageResample_S.SetInterpolationModeToCubic()
  imageResample_S.SetAxisMagnificationFactor(0, mf[0])
  imageResample_S.SetAxisMagnificationFactor(1, mf[1])
  #print "imageResample_S.GetOutput().GetDimensions() =",\
  #       imageResample_S.GetOutput().GetDimensions()
  #print "imageResample_S.GetOutput().GetSpacing() =",\
  #       imageResample_S.GetOutput().GetSpacing()
  #print "imageResample_S.GetOutputExtent() =",\
  #       imageResample_S.GetOutputExtent()

  imageMapper_S = vtk.vtkImageMapper()
  imageMapper_S.SetInput(imageResample_S.GetOutput())
  imageMapper_S.SetColorWindow(color_window)
  imageMapper_S.SetColorLevel(color_level)

  # アクターをすべて取り除く
  actor2DCollections = renderer_S.GetActors2D()
  #print 'S:actor2DCollections.GetNumberOfItems() =',
  #print actor2DCollections.GetNumberOfItems()
  for i in range(actor2DCollections.GetNumberOfItems()):
    renderer_S.RemoveActor(actor2DCollections.GetItemAsObject(i))
  actor2DCollections = renderer_S.GetActors2D()
  #print 'S:actor2DCollections.GetNumberOfItems() =',
  #print actor2DCollections.GetNumberOfItems()
  #
  actor2D_S = vtk.vtkActor2D()
  actor2D_S.SetMapper(imageMapper_S)
  renderer_S.AddActor(actor2D_S)

  #  /usr/share/vtk/Annotation/Python/TestText.pyをもとに…
  textActor_S = vtk.vtkTextActor()
  textActor_S.ScaledTextOff()
  #textActor_S.SetDisplayPosition(200,20)
  textActor_S.SetPosition(0,0)
  textActor_S.SetInput("Saggital")
  #textActor_S.GetPosition2Coordinate().SetCoordinateSystemToNormalizedViewport()
  #textActor_S.GetPosition2Coordinate().SetValue(0.6,0.1)
  textActor_S.GetTextProperty().SetFontSize(12)
  textActor_S.GetTextProperty().SetFontFamilyToArial()
  #textActor_S.GetTextProperty().SetJustificationToCentered()
  textActor_S.GetTextProperty().BoldOn()
  #textActor_S.GetTextProperty().ItalicOn()
  textActor_S.GetTextProperty().ShadowOn()
  textActor_S.GetTextProperty().SetColor(1,1,0)
  renderer_S.AddActor2D(textActor_S)

  #Coronal
  #  vtkImageReslice -> vtkImageResample -> vtkImageMapper -> vtkActor2D
  imageReslice_C = vtk.vtkImageReslice()
  imageReslice_C.SetInput(sourceProcessObject.GetOutput())
  imageReslice_C.SetOutputDimensionality(2)
  imageReslice_C.SetResliceAxesOrigin(
                   0,
                   slice_index[1]*spacing[1],
                   (dimensions[2]-1)*spacing[2])
  imageReslice_C.SetResliceAxesDirectionCosines(
                   spacing[0],0,0,
                   0,0,-spacing[2],
                   0,-spacing[1],0)
  imageReslice_C.InterpolateOn()
  imageReslice_C.SetOutputSpacing(1,1,1)  # spacingを強制的に1にする
  #
  mf = get_magnification_factor(renderer_C)
  imageResample_C = vtk.vtkImageResample()
  imageResample_C.SetInput(imageReslice_C.GetOutput())
  imageResample_C.SetDimensionality(2)
  imageResample_C.InterpolateOn()
  imageResample_C.SetInterpolationModeToCubic()
  imageResample_C.SetAxisMagnificationFactor(0, mf[0])
  imageResample_C.SetAxisMagnificationFactor(1, mf[1])
  #
  imageMapper_C = vtk.vtkImageMapper()
  imageMapper_C.SetInput(imageResample_C.GetOutput())
  imageMapper_C.SetColorWindow(color_window)
  imageMapper_C.SetColorLevel(color_level)
  # アクターをすべて取り除く
  actor2DCollections = renderer_C.GetActors2D()
  for i in range(actor2DCollections.GetNumberOfItems()):
    renderer_C.RemoveActor(actor2DCollections.GetItemAsObject(i))
  actor2DCollections = renderer_C.GetActors2D()
  #
  actor2D_C = vtk.vtkActor2D()
  actor2D_C.SetMapper(imageMapper_C)
  renderer_C.AddActor(actor2D_C)
  #
  textActor_C = vtk.vtkTextActor()
  textActor_C.ScaledTextOff()
  textActor_C.SetPosition(0,0)
  textActor_C.SetInput("Coronal")
  textActor_C.GetTextProperty().SetFontSize(12)
  textActor_C.GetTextProperty().SetFontFamilyToArial()
  textActor_C.GetTextProperty().BoldOn()
  textActor_C.GetTextProperty().ShadowOn()
  textActor_C.GetTextProperty().SetColor(1,1,0)
  renderer_C.AddActor2D(textActor_C)
########################################################################

########################################################################
#ここから
########################################################################
renderer_3=vtk.vtkRenderer()
renderer_3.SetViewport(0.0, 0.5, 0.5, 1.0)
renderer_A=vtk.vtkRenderer()
renderer_A.SetViewport(0.5, 0.0, 1.0, 0.5)
renderer_A.SetBackground(0.2,0.0,0.0)
renderer_S=vtk.vtkRenderer()
renderer_S.SetViewport(0.0, 0.0, 0.5, 0.5)
renderer_S.SetBackground(0.0,0.2,0.0)
renderer_C=vtk.vtkRenderer()
renderer_C.SetViewport(0.5, 0.5, 1.0, 1.0)
renderer_C.SetBackground(0.0,0.0,0.2)
renderWindow=vtk.vtkRenderWindow()
renderWindowInteractor=vtk.vtkRenderWindowInteractor()
renderWindowInteractor.AddObserver("UserEvent", cb_user_event)
renderWindowInteractor.AddObserver("KeyPressEvent", cb_key_press_event)
renderWindowInteractor.AddObserver("ModifiedEvent", cb_modified_event)
renderer_A.AddObserver("KeyPressEvent", cb_key_press_event)
renderWindowInteractor.SetRenderWindow(renderWindow)
renderWindow.AddRenderer(renderer_3)
renderWindow.AddRenderer(renderer_A)
renderWindow.AddRenderer(renderer_S)
renderWindow.AddRenderer(renderer_C)
renderWindow.SetSize(400, 400)
########################################################################
# ここでは、"Visualization Pipeline"を作成し、actorにmapperを設定します。

# デフォルトのprocessSourceObujectを作っておく
# default source object
# Quadric definition. This is a type of implicit function. Here the
# coefficients to the equations are set.
quadric = vtk.vtkQuadric()
quadric.SetCoefficients(1,1,1, 0,0,0, 0,0,0, -100)
# The vtkSampleFunction uses the quadric function and evaluates function
# value over a regular lattice (i.e., a volume).
sampleFunction = vtk.vtkSampleFunction()
#sampleFunction.SetSampleDimensions(21,51,101)
sampleFunction.SetSampleDimensions(51,51,51)
#sampleFunction.SetSampleDimensions(21,21,21)
sampleFunction.SetModelBounds(0,50, 0,50, 0,50)
#print "sampleFunction.GetModelBounds() =",
#print  sampleFunction.GetModelBounds()
sampleFunction.SetImplicitFunction(quadric)
sampleFunction.ComputeNormalsOff()
sampleFunction.Update()
#print "sampleFunction.GetOutput().GetSpacing() =",
#print  sampleFunction.GetOutput().GetSpacing()
#print "sampleFunction.GetOutput().GetOrigin() =",
#print  sampleFunction.GetOutput().GetOrigin()
#print "sampleFunction.GetOutput().GetDimensions() =",
#print  sampleFunction.GetOutput().GetDimensions()
origin=sampleFunction.GetOutput().GetOrigin()
#今のところは未使用(0,0,0)として計算
spacing=sampleFunction.GetOutput().GetSpacing()
dimensions=sampleFunction.GetOutput().GetDimensions()
print "dimensions =", dimensions
print "spacing =", spacing
slice_index = [int(dimension/2) for dimension in dimensions]
color_window = 2000
color_level =  1000
contour_level = 0
voi_sample_rate = 2
sourceProcessObject = sampleFunction

if False:    # ./CTのDICOMファイルをロードする処理
  # dicom source object
  dicomImageReader = vtk.vtkDICOMImageReader()
  dicomImageReader.SetDirectoryName("CT")
  dicomImageReader.Update()
  spacing = dicomImageReader.GetDataSpacing()
  dimensions = dicomImageReader.GetOutput().GetDimensions()
  slice_index = [int(dimension/2) for dimension in dimensions]
  color_window  = 400
  color_level   = 128
  contour_level = 0.0
  sourceProcessObject = dicomImageReader
  print 'spacing =', spacing
  print 'dimensions =', dimensions
elif False:  # headsqデータのロード
  # headsq
  import vtk.util.misc
  volume16Reader = vtk.vtkVolume16Reader()
  volume16Reader.SetDataDimensions(64, 64)
  volume16Reader.SetDataByteOrderToLittleEndian()
  volume16Reader.SetFilePrefix(vtk.util.misc.vtkGetDataRoot()
                    + "/Data/headsq/quarter")
  volume16Reader.SetImageRange(1, 93)
  volume16Reader.SetDataSpacing(3.2, 3.2, 1.5)
  spacing = [3.2, 3.2, 1.5]
  dimensions = [64, 64, 93]
  slice_index = [int(dimension/2) for dimension in dimensions]
  color_window  = 2000
  color_level   = 500
  contour_level = 500
  sourceProcessObject = volume16Reader
  print spacing, dimensions

initialize_visualization_pipeline()
########################################################################

#
renderWindow.Render()
renderWindowInteractor.Start()

vtk - headsq05.py

とりあえず、DICOM sample image sets から取得したデータを使ってレンダリング。

vtk - vtk_skel.py

vtk の　python スケルトンスクリプトです。
ポリンゴンで作られた球体をレンダリングするウィンドウを表示します。
レンダーウィンドウ内では、マウスを使ってカメラの位置と方向を色々かえられます。
また、キーボード操作で下に示す操作ができます。

• 'e' : 終了
• 'w': ワイヤーフレーム表示
• 's': サーフェイス表示
• '3': ステレオ表示
• 'r': 表示位置リセット
• 'p': マウスカーソルの下のオブジェクトを選択
  
• 'u': ユーザメソッド（表示画像を'image.png'という名のファイルに保存）
• 'a': アクター中心に（モデル座標での）回転？（マウスカーソルの下にあるアクターが選択される）
• 'c': ワールド座標のTransformation操作？
  

下の画像は、このスクリプトで保存したものです。

#! /usr/bin/env python
# -*- coding:utf-8 -*-
#
# vtk_skel.py

import vtk

########################################################################
# この関数では、renderWindowで"uキー"が入力されたときの処理を記述します。
#
# この例では、renderWindowに表示されている画像を"image.png"という名前の
# ファイルに保存する処理を記述しています。
# この関数では、renderWindowで"uキー"が入力されたときの処理を記述します。
#
# この例では、renderWindowに表示されている画像を"image.png"という名前の
# ファイルに保存する処理を記述しています。
def userMethod(obj, arg):
  print "userMethod"
  windowToImageFilter=vtk.vtkWindowToImageFilter()
  windowToImageFilter.SetInput(renderWindow)
  windowToImageFilter.Update()
  pngWriter=vtk.vtkPNGWriter()
  pngWriter.SetInput(windowToImageFilter.GetOutput())
  pngWriter.SetFileName("image.png")
  renderWindow.Render()
  pngWriter.Write()
########################################################################

#
renderer=vtk.vtkRenderer()
renderWindow=vtk.vtkRenderWindow()
renderWindowInteractor=vtk.vtkRenderWindowInteractor()
renderWindowInteractor.AddObserver("UserEvent", userMethod)
renderWindow.AddRenderer(renderer)
renderWindowInteractor.SetRenderWindow(renderWindow)

########################################################################
# ここでは、"Visualization Pipeline"を作成し、actorにmapperを設定します。
# この例では、SphereSourceとPolyDataMapperのpipelineを作り、
# PolyDataMapperをActorに設定しています。
sphereSource=vtk.vtkSphereSource()
sphereSource.SetThetaResolution(100)
sphereSource.SetPhiResolution(10)
polyDataMapper=vtk.vtkPolyDataMapper()
polyDataMapper.SetInput(sphereSource.GetOutput())
#
actor=vtk.vtkActor()
actor.SetMapper(polyDataMapper)
renderer.AddActor(actor)
########################################################################

#
renderWindow.Render()
renderWindowInteractor.Start()

vtk - heart [01]

vtkを使ってハートマークを作ってみました。Heart Curve というサイトで作図方法がわかったのでできました。pythonスクリプトで書いたので、そのうちスクリプトを載せようと思います（画像を保存する機能を追加したら）。

ある平面で曲面が滑らかになっていないのですが、なぜでしょう…

下の画像はワイヤーフレーム表示したものです。

IFrIT and VTK [02]

VTKDataの中のheadsqデータを、IFrITのテキスト形式に変換しました。次のようなpythonスクリプトを作って。

#!/usr/bin/env python

#
# "/Data/headsq/quarter.*" -> ifrit (*.txt)
#
import vtk
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()

# The following reader is used to read a series of 2D slices (images)
# that compose the volume. The slice dimensions are set, and the
# pixel spacing. The data Endianness must also be specified. The reader
# usese the FilePrefix in combination with the slice number to construct
# filenames using the format FilePrefix.%d. (In this case the FilePrefix
# is the root name of the file: quarter.)
v16 = vtk.vtkVolume16Reader()
v16.SetDataDimensions(64, 64)
v16.SetDataByteOrderToLittleEndian()
v16.SetFilePrefix(VTK_DATA_ROOT + "/Data/headsq/quarter")
v16.SetImageRange(1, 93)
v16.SetDataSpacing(3.2, 3.2, 1.5)

#
#
v16.Update()
ds = v16.GetOutput()
nC = ds.GetNumberOfCells()
nP = ds.GetNumberOfPoints()
nDim = ds.GetDimensions()
print nDim[0], nDim[1], nDim[2]
for z in range(nDim[2]):
  for y in range(nDim[1]):
    for x in range(nDim[0]):
      print ds.GetScalarComponentAsFloat(x,y,z,0)

画像はそのテキスト形式データをIFrITで表示したものです。
縦方向のスケールがおかしいのは仕方ないか…

IFrIT and VTK [01]

VisualizationツールVTKとifritをためしてみました。

なかなか面白い。

本家をみてみると、最新バージョン(3.1.4)では、拡張ソースを使うと、VTKのファイルを読めるらしい。なのでソースを落してビルド（debianパッケージ版のソースのパッチが役に立った）。

で、VTKDataに入っているファイルを開けるようになったのですが、あまり面白いファイルがみあたらない。Medicalサンプルは、Data/headsq/quarterのバイナリファイルを読み込んでいるらしい。これをIFrITで試してみたくなったので、仕方無くVTKを使うことに…でも、pythonで使えるからまだ楽なのかな。pythonもあまり良くわかってないけど、こんな感じのプログラムを作って、

#!/usr/bin/env python

#
# "/Data/headsq/quarter.*" -> "headsq.vtk"
#
import vtk
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()

# The following reader is used to read a series of 2D slices (images)
# that compose the volume. The slice dimensions are set, and the
# pixel spacing. The data Endianness must also be specified. The reader
# usese the FilePrefix in combination with the slice number to construct
# filenames using the format FilePrefix.%d. (In this case the FilePrefix
# is the root name of the file: quarter.)
v16 = vtk.vtkVolume16Reader()
v16.SetDataDimensions(64, 64)
v16.SetDataByteOrderToLittleEndian()
v16.SetFilePrefix(VTK_DATA_ROOT + "/Data/headsq/quarter")
v16.SetImageRange(1, 93)
v16.SetDataSpacing(3.2, 3.2, 1.5)

#
w = vtk.vtkStructuredPointsWriter()
w.SetFileName("headsq.vtk")
w.SetInput(v16.GetOutput())
w.Write()

headsq.vtkをIFrITで読み込んでみました。

かなり不気味な画像。
しかし、これでも完全ではないみたい。サーフェイスレンダリングはいいけど、クロスセクションができないし、ボリュームレンダリングのときのレンジがいまいちだったりする。
なのでpython-vtkを使ってIFrITのデータ形式をつくり出すことを考えてます。