Visualisation tools

visualisation/CMakeLists.txt

@@ -0,0 +1,37 @@
+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+
+project(GridViewer)
+
+list(APPEND CMAKE_PREFIX_PATH "/Users/peterboyle/QCD/vtk/VTK-9.4.2-install/")
+
+find_package(VTK COMPONENTS 
+  CommonColor
+  CommonCore
+  FiltersCore
+  FiltersModeling
+  IOImage
+  IOFFMPEG
+  InteractionStyle
+  InteractionWidgets
+  RenderingContextOpenGL2
+  RenderingCore
+  RenderingFreeType
+  RenderingGL2PSOpenGL2
+  RenderingOpenGL2
+)
+
+if (NOT VTK_FOUND)
+  message(FATAL_ERROR "GridViewer: Unable to find the VTK build folder.")
+endif()
+
+# Prevent a "command line is too long" failure in Windows.
+set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
+
+add_executable(FieldDensityAnimate MACOSX_BUNDLE FieldDensityAnimate.cxx )
+  target_link_libraries(FieldDensityAnimate PRIVATE ${VTK_LIBRARIES}
+)
+# vtk_module_autoinit is needed
+vtk_module_autoinit(
+  TARGETS FieldDensityAnimate
+  MODULES ${VTK_LIBRARIES}
+)

visualisation/FieldDensity.py

@@ -0,0 +1,285 @@
+#!/usr/bin/env python
+
+# noinspection PyUnresolvedReferences
+import math
+import vtk
+import vtkmodules.vtkInteractionStyle
+# noinspection PyUnresolvedReferences
+import vtkmodules.vtkRenderingOpenGL2
+from vtkmodules.vtkCommonColor import vtkNamedColors
+from vtkmodules.vtkCommonCore import (
+    VTK_VERSION_NUMBER,
+    vtkVersion
+)
+from vtkmodules.vtkCommonCore import VTK_DOUBLE
+from vtkmodules.vtkCommonDataModel import vtkImageData
+from vtkmodules.vtkFiltersCore import (
+    vtkMarchingCubes,
+    vtkStripper
+)
+from vtkmodules.vtkFiltersModeling import vtkOutlineFilter
+from vtkmodules.vtkIOImage import (
+    vtkMetaImageReader,
+    vtkJPEGWriter,
+    vtkPNGWriter
+)
+from vtkmodules.vtkRenderingCore import (
+    vtkActor,
+    vtkCamera,
+    vtkPolyDataMapper,
+    vtkProperty,
+    vtkRenderWindow,
+    vtkRenderWindowInteractor,
+    vtkRenderer,
+    vtkWindowToImageFilter
+)
+
+
+class vtkTimerCallback():
+    def __init__(self, steps, imageData, iren):
+        self.timer_count = 0
+        self.steps = steps
+        self.imageData = imageData
+        self.iren = iren
+        self.timerId = None
+        self.step = 0
+
+    def execute(self, obj, event):
+
+        print(self.timer_count)
+
+        dims = self.imageData.GetDimensions()
+
+        t=self.step/10.0
+
+        z0 = 2
+        y0 = 4
+        x0 = 4
+        z1 = 14
+        y1 = 12
+        x1 = 12
+        for z in range(dims[2]):
+            for y in range(dims[1]):
+                for x in range(dims[0]):
+                    self.imageData.SetScalarComponentFromDouble(x, y, z, 0,
+                                                                math.sin(t)*math.exp(-0.25*((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0)))
+                                                                - math.cos(t)*math.exp(-0.25*((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1))))
+
+        self.imageData.Modified()
+
+        iren = obj
+        iren.GetRenderWindow().Render()
+        self.timer_count += 1
+        self.step += 1
+
+        if self.step >= self.steps :
+            iren.DestroyTimer(self.timerId)
+
+def WriteImage(fileName, renWin):
+    '''
+    '''
+
+    import os
+
+    if fileName:
+        # Select the writer to use.
+        path, ext = os.path.splitext(fileName)
+        ext = ext.lower()
+        if not ext:
+            ext = '.png'
+            fileName = fileName + ext
+        elif ext == '.jpg':
+            writer = vtkJPEGWriter()
+        else:
+            writer = vtkPNGWriter()
+
+        windowto_image_filter = vtkWindowToImageFilter()
+        windowto_image_filter.SetInput(renWin)
+        windowto_image_filter.SetScale(1)  # image quality
+        windowto_image_filter.SetInputBufferTypeToRGBA()
+
+        writer.SetFileName(fileName)
+        writer.SetInputConnection(windowto_image_filter.GetOutputPort())
+        writer.Write()
+    else:
+        raise RuntimeError('Need a filename.')
+
+
+def main():
+
+    colors = vtkNamedColors()
+
+    file_name = get_program_parameters()
+
+    colors.SetColor('InstantonColor', [240, 184, 160, 255])
+    colors.SetColor('BackfaceColor', [255, 229, 200, 255])
+    colors.SetColor('BkgColor', [51, 77, 102, 255])
+
+    # Create the renderer, the render window, and the interactor. The renderer
+    # draws into the render window, the interactor enables mouse- and
+    # keyboard-based interaction with the data within the render window.
+    #
+    a_renderer = vtkRenderer()
+    ren_win = vtkRenderWindow()
+    ren_win.AddRenderer(a_renderer)
+
+    iren = vtkRenderWindowInteractor()
+    iren.SetRenderWindow(ren_win)
+
+    # 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
+    # uses 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.)
+    imageData = vtkImageData()
+    imageData.SetDimensions(16, 16, 16)
+    imageData.AllocateScalars(VTK_DOUBLE, 1)
+
+    dims = imageData.GetDimensions()
+
+    # Fill every entry of the image data with '2.0'
+    # Set the input data
+    for z in range(dims[2]):
+        z0 = dims[2]/2
+        for y in range(dims[1]):
+            y0 = dims[1]/2
+            for x in range(dims[0]):
+                x0 = dims[0]/2
+                imageData.SetScalarComponentFromDouble(x, y, z, 0, math.exp(-0.25*((x-x0)*(x-x0)+(y-y0)*(y-y0)+z*z)) -  math.exp(-0.25*((x-x0)*(x-x0)+y*y+(z-z0)*(z-z0))))
+
+    instanton_extractor = vtkMarchingCubes()
+    instanton_extractor.SetInputData(imageData)
+    instanton_extractor.SetValue(0, 0.1)
+    
+    instanton_stripper = vtkStripper()
+    instanton_stripper.SetInputConnection(instanton_extractor.GetOutputPort())
+    
+    instanton_mapper = vtkPolyDataMapper()
+    instanton_mapper.SetInputConnection(instanton_stripper.GetOutputPort())
+    instanton_mapper.ScalarVisibilityOff()
+    
+    instanton = vtkActor()
+    instanton.SetMapper(instanton_mapper)
+    instanton.GetProperty().SetDiffuseColor(colors.GetColor3d('InstantonColor'))
+    instanton.GetProperty().SetSpecular(0.3)
+    instanton.GetProperty().SetSpecularPower(20)
+    instanton.GetProperty().SetOpacity(0.5)
+
+    # The triangle stripper is used to create triangle strips from the
+    # isosurface these render much faster on may systems.
+    antiinstanton_extractor = vtkMarchingCubes()
+    antiinstanton_extractor.SetInputData(imageData)
+    antiinstanton_extractor.SetValue(0, -0.1)
+    
+    antiinstanton_stripper = vtkStripper()
+    antiinstanton_stripper.SetInputConnection(antiinstanton_extractor.GetOutputPort())
+    
+    antiinstanton_mapper = vtkPolyDataMapper()
+    antiinstanton_mapper.SetInputConnection(antiinstanton_stripper.GetOutputPort())
+    antiinstanton_mapper.ScalarVisibilityOff()
+    
+    antiinstanton = vtkActor()
+    antiinstanton.SetMapper(antiinstanton_mapper)
+    antiinstanton.GetProperty().SetDiffuseColor(colors.GetColor3d('Ivory'))
+
+    # An outline provides box around the data.
+    outline_data = vtkOutlineFilter()
+    outline_data.SetInputData(imageData)
+    
+    map_outline = vtkPolyDataMapper()
+    map_outline.SetInputConnection(outline_data.GetOutputPort())
+    
+    outline = vtkActor()
+    outline.SetMapper(map_outline)
+    outline.GetProperty().SetColor(colors.GetColor3d('Black'))
+    
+    # It is convenient to create an initial view of the data. The FocalPoint
+    # and Position form a vector direction. Later on (ResetCamera() method)
+    # this vector is used to position the camera to look at the data in
+    # this direction.
+    a_camera = vtkCamera()
+    a_camera.SetViewUp(0, 0, -1)
+    a_camera.SetPosition(0, -100, 0)
+    a_camera.SetFocalPoint(0, 0, 0)
+    a_camera.ComputeViewPlaneNormal()
+    a_camera.Azimuth(30.0)
+    a_camera.Elevation(30.0)
+
+    # Actors are added to the renderer. An initial camera view is created.
+    # The Dolly() method moves the camera towards the FocalPoint,
+    # thereby enlarging the image.
+    a_renderer.AddActor(outline)
+    a_renderer.AddActor(instanton)
+    a_renderer.AddActor(antiinstanton)
+    a_renderer.SetActiveCamera(a_camera)
+    a_renderer.ResetCamera()
+    a_camera.Dolly(1.0)
+
+    # Set a background color for the renderer and set the size of the
+    # render window (expressed in pixels).
+    a_renderer.SetBackground(colors.GetColor3d('BkgColor'))
+    ren_win.SetSize(1024, 1024)
+    ren_win.SetWindowName('ExpoDemo')
+
+    # Note that when camera movement occurs (as it does in the Dolly()
+    # method), the clipping planes often need adjusting. Clipping planes
+    # consist of two planes: near and far along the view direction. The
+    # near plane clips out objects in front of the plane the far plane
+    # clips out objects behind the plane. This way only what is drawn
+    # between the planes is actually rendered.
+    a_renderer.ResetCameraClippingRange()
+
+    # write image
+    #    WriteImage('exp.jpg',ren_win)
+
+    # Sign up to receive TimerEvent
+    cb = vtkTimerCallback(200, imageData, iren)
+    iren.AddObserver('TimerEvent', cb.execute)
+    cb.timerId = iren.CreateRepeatingTimer(50)
+
+    # start the interaction and timer
+    ren_win.Render()
+    
+    # Initialize the event loop and then start it.
+    iren.Initialize()
+    iren.Start()
+
+
+def get_program_parameters():
+    import argparse
+    description = 'Simple lattice volumetric demo'
+    epilogue = '''
+    Derived from VTK/Examples/Cxx/Medical2.cxx
+    '''
+    parser = argparse.ArgumentParser(description=description, epilog=epilogue,
+                                     formatter_class=argparse.RawDescriptionHelpFormatter)
+    parser.add_argument('filename', help='FieldDensity.py')
+    args = parser.parse_args()
+    return args.filename
+
+
+def vtk_version_ok(major, minor, build):
+    """
+    Check the VTK version.
+
+    :param major: Major version.
+    :param minor: Minor version.
+    :param build: Build version.
+    :return: True if the requested VTK version is greater or equal to the actual VTK version.
+    """
+    needed_version = 10000000000 * int(major) + 100000000 * int(minor) + int(build)
+    try:
+        vtk_version_number = VTK_VERSION_NUMBER
+    except AttributeError:  # as error:
+        ver = vtkVersion()
+        vtk_version_number = 10000000000 * ver.GetVTKMajorVersion() + 100000000 * ver.GetVTKMinorVersion() \
+                             + ver.GetVTKBuildVersion()
+    if vtk_version_number >= needed_version:
+        return True
+    else:
+        return False
+
+
+if __name__ == '__main__':
+    main()

visualisation/FieldDensityAnimate.cxx

@@ -0,0 +1,490 @@
+// Derived from VTK/Examples/Cxx/Medical2.cxx
+// The example reads a volume dataset, extracts two isosurfaces that
+// represent the skin and bone, and then displays them.
+//
+// Modified heavily by Peter Boyle to display lattice field theory data as movies and compare multiple files
+
+#include <vtkActor.h>
+#include <vtkCamera.h>
+#include <vtkMetaImageReader.h>
+#include <vtkNamedColors.h>
+#include <vtkNew.h>
+#include <vtkOutlineFilter.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkProperty.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+#include <vtkStripper.h>
+#include <vtkImageData.h>
+#include <vtkVersion.h>
+#include <vtkCallbackCommand.h>
+#include <vtkTextActor.h>
+#include <vtkTextProperty.h>
+
+#define MPEG
+#ifdef MPEG
+#include <vtkFFMPEGWriter.h>
+#endif
+
+#include <vtkProperty2D.h>
+#include <vtkSliderWidget.h>
+#include <vtkSliderRepresentation2D.h>
+#include <vtkWindowToImageFilter.h>
+
+#include <array>
+#include <string>
+
+#include <Grid/Grid.h>
+
+#define USE_FLYING_EDGES
+#ifdef USE_FLYING_EDGES
+#include <vtkFlyingEdges3D.h>
+typedef vtkFlyingEdges3D isosurface;
+#else
+#include <vtkMarchingCubes.h>
+typedef vtkMarchingCubes isosurface;
+#endif
+
+int mpeg = 0 ;
+int xlate = 0 ;
+
+template <class T> void readFile(T& out, std::string const fname){
+#ifdef HAVE_LIME
+  Grid::emptyUserRecord record;
+  Grid::ScidacReader RD;
+  RD.open(fname);
+  RD.readScidacFieldRecord(out,record);
+  RD.close();
+#endif
+}
+using namespace Grid;
+
+class FrameUpdater : public vtkCallbackCommand
+{
+public:
+
+  FrameUpdater() {
+    TimerCount = 0;
+    xoff       = 0;
+    t          = 0;
+    imageData = nullptr;
+    grid_data = nullptr;
+    timerId = 0;
+    maxCount = -1;
+  }
+  
+  static FrameUpdater* New()
+  {
+    FrameUpdater* cb = new FrameUpdater;
+    cb->TimerCount = 0;
+    return cb;
+  }
+
+  virtual void Execute(vtkObject* caller, unsigned long eventId,void* vtkNotUsed(callData))
+  {
+    const int max=256;
+    char text_string[max];
+
+    if (this->TimerCount < this->maxCount) {
+
+      if (vtkCommand::TimerEvent == eventId)
+	{
+	  ++this->TimerCount;
+	  
+	  // Make a new frame
+	  auto latt_size = grid_data->Grid()->GlobalDimensions();
+	  for(int xx=0;xx<latt_size[0];xx++){
+	    for(int yy=0;yy<latt_size[1];yy++){
+	      for(int zz=0;zz<latt_size[2];zz++){
+		int x = (xx+xoff)%latt_size[0];
+		Coordinate site({x,yy,zz,t});
+		RealD value = real(peekSite(*grid_data,site));
+		imageData->SetScalarComponentFromDouble(xx,yy,zz,0,value);
+	  }}}
+
+	  if ( xlate ) { 
+	    xoff = (xoff + 1)%latt_size[0];
+	    if ( xoff== 0 ) t = (t+1)%latt_size[3];
+	  } else {
+	    t = (t+1)%latt_size[3];
+	    if ( t== 0 ) 	xoff = (xoff + 1)%latt_size[0];
+	  }
+
+	  snprintf(text_string,max,"T=%d",t);
+	  text->SetInput(text_string);
+      
+	  std::cout << this->TimerCount<<"/"<<maxCount<< " xoff "<<xoff<<" t "<<t  <<std::endl;
+	  imageData->Modified();
+
+	  vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
+	  iren->GetRenderWindow()->Render();
+	  
+	}
+    }
+    
+    if (this->TimerCount >= this->maxCount) {
+      vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
+      if (this->timerId > -1)
+      {
+        iren->DestroyTimer(this->timerId);
+      }
+    }
+  }
+  
+
+private:
+  int TimerCount;
+  int xoff;
+  int t;
+public:
+  Grid::LatticeComplexD * grid_data;
+  vtkImageData* imageData = nullptr;
+  vtkTextActor* text = nullptr;
+  vtkFFMPEGWriter *writer = nullptr;
+  int timerId ;
+  int maxCount ;
+  double rms;
+  isosurface * posExtractor;
+  isosurface * negExtractor;
+};
+class SliderCallback : public vtkCommand
+{
+public:
+    static SliderCallback* New()
+    {
+        return new SliderCallback;
+    }
+    virtual void Execute(vtkObject* caller, unsigned long eventId, void* callData)
+    {
+        vtkSliderWidget *sliderWidget = vtkSliderWidget::SafeDownCast(caller);
+        if (sliderWidget)
+        {
+	  contour = ((vtkSliderRepresentation *)sliderWidget->GetRepresentation())->GetValue();
+        }
+	for(int i=0;i<fu_list.size();i++){
+	  fu_list[i]->posExtractor->SetValue(0,  SliderCallback::contour*fu_list[i]->rms);
+	  fu_list[i]->negExtractor->SetValue(0, -SliderCallback::contour*fu_list[i]->rms);
+	  fu_list[i]->posExtractor->Modified();
+	  fu_list[i]->negExtractor->Modified();
+	}
+    }
+public:
+  static double contour;
+  std::vector<FrameUpdater *> fu_list;
+};
+
+
+double SliderCallback::contour;
+
+int main(int argc, char* argv[])
+{
+  using namespace Grid;
+
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+  GridCartesian    Grid(latt_size, simd_layout, mpi_layout);
+
+ 
+  std::cout << argc << " command Line arguments "<<std::endl;
+  for(int c=0;c<argc;c++) {
+    std::cout << " - "<<argv[c]<<std::endl;
+  }
+
+  std::vector<std::string> file_list;
+  double default_contour = 1.0;
+  std::string arg;
+#ifdef MPEG
+  if( GridCmdOptionExists(argv,argv+argc,"--mpeg") ){
+    mpeg = 1;
+  }
+#endif
+
+  if( GridCmdOptionExists(argv,argv+argc,"--xlate") ){
+    xlate = 1;
+  }
+
+  if( GridCmdOptionExists(argv,argv+argc,"--isosurface") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--isosurface");
+    GridCmdOptionFloat(arg,default_contour);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--file1") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--file1");
+    file_list.push_back(arg);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--file2") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--file2");
+    file_list.push_back(arg);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--file3") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--file3");
+    file_list.push_back(arg);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--file4") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--file4");
+    file_list.push_back(arg);
+  }
+  for(int c=0;c<file_list.size();c++) {
+    std::cout << " file: "<<file_list[c]<<std::endl;
+  }
+
+  // Common things:
+  vtkNew<vtkNamedColors> colors;
+  std::array<unsigned char, 4> posColor{{240, 184, 160, 255}};  colors->SetColor("posColor", posColor.data());
+  std::array<unsigned char, 4> bkg{{51, 77, 102, 255}};         colors->SetColor("BkgColor", bkg.data());
+
+  // Create the renderer, the render window, and the interactor. The renderer
+  // draws into the render window, the interactor enables mouse- and
+  // keyboard-based interaction with the data within the render window.
+  //
+  vtkNew<vtkRenderWindow> renWin;
+  vtkNew<vtkRenderWindowInteractor> iren;
+  iren->SetRenderWindow(renWin);
+
+  
+  std::vector<LatticeComplexD> data(file_list.size(),&Grid);
+  FieldMetaData header;
+
+  
+  int frameCount;
+  if ( mpeg ) frameCount = latt_size[3];
+  else        frameCount = latt_size[0] * latt_size[3];
+
+  std::vector<FrameUpdater *> fu_list;
+  for (int f=0;f<file_list.size();f++){
+
+    // It is convenient to create an initial view of the data. The FocalPoint
+    // and Position form a vector direction. Later on (ResetCamera() method)
+    // this vector is used to position the camera to look at the data in
+    // this direction.
+    vtkNew<vtkCamera> aCamera;
+    aCamera->SetViewUp(0, 0, -1);
+    aCamera->SetPosition(0, -1000, 0);
+    aCamera->SetFocalPoint(0, 0, 0);
+    aCamera->ComputeViewPlaneNormal();
+    aCamera->Azimuth(30.0);
+    aCamera->Elevation(30.0);
+
+    
+    vtkNew<vtkRenderer> aRenderer;
+    renWin->AddRenderer(aRenderer);
+    
+    double vol = data[f].Grid()->gSites();
+
+    std::cout << "Reading "<<file_list[f]<<std::endl;
+    readFile(data[f],file_list[f]);
+
+    auto nrm    = norm2(data[f]);
+    auto nrmbar = nrm/vol;
+    auto rms    = sqrt(nrmbar);
+
+    double contour = default_contour * rms; // default to 1 x RMS
+
+    // 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
+    // uses 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.)
+    vtkNew<vtkImageData> imageData;
+    imageData->SetDimensions(latt_size[0],latt_size[1],latt_size[2]);
+    imageData->AllocateScalars(VTK_DOUBLE, 1);
+    for(int xx=0;xx<latt_size[0];xx++){
+      for(int yy=0;yy<latt_size[1];yy++){
+	for(int zz=0;zz<latt_size[2];zz++){
+	  Coordinate site({xx,yy,zz,0});
+	  RealD value = real(peekSite(data[f],site));
+	  imageData->SetScalarComponentFromDouble(xx,yy,zz,0,value);
+    }}}
+
+    vtkNew<isosurface> posExtractor;
+    posExtractor->SetInputData(imageData);
+    posExtractor->SetValue(0, contour);
+  
+    vtkNew<vtkStripper> posStripper;
+    posStripper->SetInputConnection(posExtractor->GetOutputPort());
+
+    vtkNew<vtkPolyDataMapper> posMapper;
+    posMapper->SetInputConnection(posStripper->GetOutputPort());
+    posMapper->ScalarVisibilityOff();
+
+    vtkNew<vtkActor> pos;
+    pos->SetMapper(posMapper);
+    pos->GetProperty()->SetDiffuseColor(colors->GetColor3d("posColor").GetData());
+    pos->GetProperty()->SetSpecular(0.3);
+    pos->GetProperty()->SetSpecularPower(20);
+    pos->GetProperty()->SetOpacity(0.5);
+
+    // An isosurface, or contour value is set
+    // The triangle stripper is used to create triangle strips from the
+    // isosurface; these render much faster on may systems.
+    vtkNew<isosurface> negExtractor;
+    negExtractor->SetInputData(imageData);
+    negExtractor->SetValue(0, -contour);
+
+    vtkNew<vtkStripper> negStripper;
+    negStripper->SetInputConnection(negExtractor->GetOutputPort());
+
+    vtkNew<vtkPolyDataMapper> negMapper;
+    negMapper->SetInputConnection(negStripper->GetOutputPort());
+    negMapper->ScalarVisibilityOff();
+
+    vtkNew<vtkActor> neg;
+    neg->SetMapper(negMapper);
+    neg->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData());
+
+    // An outline provides context around the data.
+    vtkNew<vtkOutlineFilter> outlineData;
+    outlineData->SetInputData(imageData);
+
+    vtkNew<vtkPolyDataMapper> mapOutline;
+    mapOutline->SetInputConnection(outlineData->GetOutputPort());
+
+    vtkNew<vtkActor> outline;
+    outline->SetMapper(mapOutline);
+    outline->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());
+
+    vtkNew<vtkTextActor> Text;
+    Text->SetInput(file_list[f].c_str());
+    Text->SetPosition2(0,0);
+    Text->GetTextProperty()->SetFontSize(48);
+    Text->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
+
+    vtkNew<vtkTextActor> TextT;
+    TextT->SetInput("T=0");
+    TextT->SetPosition(0,.9*1025);
+    TextT->GetTextProperty()->SetFontSize(48);
+    TextT->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
+    
+  
+    // Actors are added to the renderer. An initial camera view is created.
+    // The Dolly() method moves the camera towards the FocalPoint,
+    // thereby enlarging the image.
+    aRenderer->AddActor(Text);
+    aRenderer->AddActor(TextT);
+    aRenderer->AddActor(outline);
+    aRenderer->AddActor(pos);
+    aRenderer->AddActor(neg);
+
+    // Sign up to receive TimerEvent
+    vtkNew<FrameUpdater> fu;
+    fu->imageData = imageData;
+    fu->grid_data = &data[f];
+    fu->text      = TextT;
+    fu->maxCount = frameCount;
+    fu->posExtractor = posExtractor;
+    fu->negExtractor = negExtractor;
+    fu->rms = rms;
+      
+    iren->AddObserver(vtkCommand::TimerEvent, fu);
+
+    aRenderer->SetActiveCamera(aCamera);
+    aRenderer->ResetCamera();
+    aRenderer->SetBackground(colors->GetColor3d("BkgColor").GetData());
+    aCamera->Dolly(1.0);
+
+    double nf = file_list.size();
+    std::cout << " Adding renderer " <<f<<" of "<<nf<<std::endl;
+    aRenderer->SetViewport((1.0/nf)*f, 0.0,(1.0/nf)*(f+1) , 1.0);
+
+    // Note that when camera movement occurs (as it does in the Dolly()
+    // method), the clipping planes often need adjusting. Clipping planes
+    // consist of two planes: near and far along the view direction. The
+    // near plane clips out objects in front of the plane; the far plane
+    // clips out objects behind the plane. This way only what is drawn
+    // between the planes is actually rendered.
+    aRenderer->ResetCameraClippingRange();
+    
+    fu_list.push_back(fu);
+  }
+
+
+  // Set a background color for the renderer and set the size of the
+  // render window (expressed in pixels).
+  // Initialize the event loop and then start it.
+  renWin->SetSize(1024*file_list.size(), 1024);
+  renWin->SetWindowName("FieldDensity");
+  renWin->Render();
+
+  iren->Initialize();
+
+  if ( mpeg ) {
+#ifdef MPEG
+    vtkWindowToImageFilter *imageFilter = vtkWindowToImageFilter::New();
+    imageFilter->SetInput( renWin );
+    imageFilter->SetInputBufferTypeToRGB();
+    
+    vtkFFMPEGWriter *writer = vtkFFMPEGWriter::New();
+    writer->SetFileName("movie.avi");
+    writer->SetRate(1);
+    writer->SetInputConnection(imageFilter->GetOutputPort());
+    writer->Start();
+
+    for(int i=0;i<fu_list[0]->maxCount;i++){
+      for(int f=0;f<fu_list.size();f++){
+	fu_list[f]->Execute(iren,vtkCommand::TimerEvent,nullptr);
+      }
+      imageFilter->Modified();
+      writer->Write();
+    }
+    writer->End();
+    writer->Delete();
+#else
+    assert(-1 && "MPEG support not compiled");
+#endif
+  } else { 
+  
+    // Add control of contour threshold
+    // Create a slider widget
+    vtkSmartPointer<vtkSliderRepresentation2D> sliderRep = vtkSmartPointer<vtkSliderRepresentation2D>::New();
+    sliderRep->SetMinimumValue(0.1);
+    sliderRep->SetMaximumValue(5.0);
+    sliderRep->SetValue(1.0);
+    sliderRep->SetTitleText("Fraction RMS");
+    // Set color properties:
+
+    // Change the color of the knob that slides
+    //  sliderRep->GetSliderProperty()->SetColor(colors->GetColor3d("Green").GetData());
+    sliderRep->GetTitleProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
+    sliderRep->GetLabelProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
+    sliderRep->GetSelectedProperty()->SetColor(colors->GetColor3d("DeepPink").GetData());
+
+    // Change the color of the bar
+    sliderRep->GetTubeProperty()->SetColor(colors->GetColor3d("MistyRose").GetData());
+    sliderRep->GetCapProperty()->SetColor(colors->GetColor3d("Yellow").GetData());
+    sliderRep->SetSliderLength(0.05);
+    sliderRep->SetSliderWidth(0.025);
+    sliderRep->SetEndCapLength(0.02);
+
+    double nf = file_list.size();
+    sliderRep->GetPoint1Coordinate()->SetCoordinateSystemToNormalizedDisplay();
+    sliderRep->GetPoint1Coordinate()->SetValue(0.1, 0.1);
+    sliderRep->GetPoint2Coordinate()->SetCoordinateSystemToNormalizedDisplay();
+    sliderRep->GetPoint2Coordinate()->SetValue(0.9/nf, 0.1);
+  
+    vtkSmartPointer<vtkSliderWidget> sliderWidget = vtkSmartPointer<vtkSliderWidget>::New();
+    sliderWidget->SetInteractor(iren);
+    sliderWidget->SetRepresentation(sliderRep);
+    sliderWidget->SetAnimationModeToAnimate();
+    sliderWidget->EnabledOn();
+  
+    // Create the slider callback
+    vtkSmartPointer<SliderCallback> slidercallback = vtkSmartPointer<SliderCallback>::New();
+    slidercallback->fu_list = fu_list;
+    sliderWidget->AddObserver(vtkCommand::InteractionEvent, slidercallback);
+
+    int timerId = iren->CreateRepeatingTimer(300);
+    std::cout << "timerId: " << timerId << std::endl;
+
+    // Start the interaction and timer
+    iren->Start();
+  }
+
+  Grid_finalize();
+
+  return EXIT_SUCCESS;
+}

visualisation/README

@@ -0,0 +1,87 @@
+========================================
+Visualisation of Grid / SciDAC format density fields using VTK (visualisation toolkit). Peter Boyle, 2025.
+========================================
+
+Uses:
+
+https://vtk.org
+
+Files are, for example, those produced by
+     Grid/HMC/ComputeWilsonFlow.cc
+and
+     Grid/HMC/site_plaquette.cc
+
+========================================
+Prerequisites:
+========================================
+
+
+1) Install ffmpeg-7.0.2  (developer install, includes headers and libraries).
+   MacOS note: must install ffmpeg from source -- homebrew only installs binaries.
+
+   https://ffmpeg.org/download.html#releases
+
+
+   Note: the latest ffmpeg (7.1.1) broke software compatibility with VTK.
+
+
+2) Build and install VTK-9.4.2, with FFMEG support enabled.
+
+   This is particularly involved on MacOS, so documented here.
+
+   cd VTK-9.4.2
+   mkdir build
+   cd build
+   ccmake ..
+
+
+   Using cmake editor, set:
+
+   FFMPEG_DIR                       /usr/local      
+
+   toggle "advanced mode"
+
+
+   CMAKE_EXE_LINKER_FLAGS to:
+   
+   -framework Foundation -framework AudioToolbox -framework CoreAudio -liconv -lm -framework AVFoundation -framework CoreVideo -framework CoreMedia -framework CoreGraphics -framework AudioToolbox -framework OpenGL -framework OpenGL -framework VideoToolbox -framework CoreImage -framework AppKit -framework CoreFoundation -framework CoreServices -lz -lbz2 -Wl,-framework,CoreFoundation -Wl,-framework,Security -L/usr/local/lib -lavdevice -lavfilter -lavformat -lavcodec -lswresample -lswscale -lavutil
+
+   VTK really should make it easier to pick up the flags required for FFMPEG linkage, especially as they are very quirky on MacOS.
+
+========================================
+Grid:
+========================================
+
+3) Build and install a version of Grid
+
+4) Ensure "grid-config" is in your path.
+
+5) cd Grid/visualisation/
+   libs=`grid-config --libs`
+   ldflags=`grid-config --ldflags`
+   cxxflags=`grid-config --cxxflags`
+   cxx=`grid-config --cxx`
+
+   mkdir build
+   cd build
+
+   LDFLAGS="$ldflags $libs " cmake .. -DCMAKE_CXX_COMPILER=$cxx -DCMAKE_CXX_FLAGS=$cxxflags 
+
+   make
+   
+6) Invoke as:
+
+   FieldDensityAnimate --isosurface <float-from-0-to-5> --grid X.Y.Z.T --file1 SciDacDensityFile1 [--xlate] [--mpeg] 
+   FieldDensityAnimate --isosurface <float-from-0-to-5> --grid X.Y.Z.T --file1 SciDacDensityFile1 --file2 SciDacDensityFile2 [--xlate] [--mpeg] 
+
+==================================
+Extensions
+==================================
+
+7) Direct calling from Grid ?:
+
+   Not yet implemented, but could develop sufficient interface to write a Lattice scalar field into MPEG direct from running code.
+
+8) Example python code: FieldDensity.py . This is not interfaced to Grid.
+
+

visualisation/cmake-command

@@ -0,0 +1,9 @@
+libs=`grid-config --libs`
+ldflags=`grid-config --ldflags`
+cxxflags=`grid-config --cxxflags`
+cxx=`grid-config --cxx`
+
+mkdir build
+cd build
+
+LDFLAGS="$ldflags $libs " cmake .. -DCMAKE_CXX_COMPILER=$cxx -DCMAKE_CXX_FLAGS=$cxxflags