Friday, June 19, 2015

VTK 5 and FFMPEG shenanigans on Windows

Over the years, trying to get VTK 5 to work with FFMPEG has always been tricky. Thankfully, VTK 6 has no such problems, but for those that have legacy code that use VTK 5, this is still a problem.

In this post I will cover some necessary tweaks, steps and links to getting FFMPEG IO working with your project that uses VTK 5, such as the 'Animate' plugin for sMILX and SMILI.

Firstly, I have found VTK 5.8 to be far more stable and workable than 5.10.1. In my scientific visualisation software SMILI, certain features don't work at all with the same code that runs fine on VTK 6 and 5.8.

Secondly, it seems VTK 5 has issues with Visual Studio 2013 (required to support Windows 8), but a patched version of VTK 5.10.1 available on GitHub works OK, though I had to fix a few round brackets here and here.

For version 1.2 of FFMPEG, this is what I recently did to get VTK 5.8 working with FFMPEG. The best link for FFMPEG deprecation changes are found here. Summarised, I had to fix vtkFFMPEGWriter.cxx as the following:
  1. av_set_parameters  - comment out
  2. URL_WRONLY to AVIO_FLAG_WRITE
  3. url_fopen to avio_open and url_fclose to avio_close
  4. av_write_header to avformat_write_header, last arg to NULL
No changes are required for VTK 6.1.0 at the time of writing.

Hope that helps
Cheers Shakes - L3mming

Friday, June 5, 2015

Mesh/Surface Processing with Command-line SMILI tools

I recently had to do some processing for a prostate segmentation scheme I've been working on and thought some uses of SMILI tools would be helpful to others. Majority of surface processing tasks can be achieved with the milxModelApp in SMILI, a easy to use open source medical imaging library and viewer.

My tasks involved involved:
  1. copying scalars from an atlas mesh to lots of other meshes
  2. using scalars to clip the meshes
  3. computing the surface distance of the meshes to their corresponding labelled images that were created by clinical experts
  4. computing the scalar statistics across the population of meshes (my meshes already have point to point correspondence) 
Luckily, the milxModelApp application provides most of these without needing to script for multiple surfaces. Furthermore, the Qt framework is not required for these applications as they are VTK/ITK only applications (i.e. built on milxSMILI, the GUI-independent sub-library of SMILI).The following examples are for Version 0.998 (SMILI Release Beta 1).

Each one of the above tasks can be achieved as:
  1. The --scalarcopy  option in milxModelApp provides a way to do this. I did:
    milxModelApp --scalarcopy focus_atlases/focus_bladder_atlas.vtk results/bladder/asm_bladder_*.vtk -p results_scalars/bladder/
    Copy scalars from atlas mesh (focus_bladder_atlas.vtk) to other meshes in directory (results/bladder/asm_bladder_*.vtk) and output to another directory with same names.
  2. The --clip  option in milxModelApp provides a way to do this. I did:
    milxModelApp --clip 1 results_scalars/bladder/asm_bladder_*.vtk -p results_clipped/bladder/
    Clip meshes in directory keeping only parts with value of 1 and output to another directory with same names.
  3. We need to write a script for the milxHausdorffDistance application that provides this using distance maps in a straight forward way. This script is now in the repository as scripts/ batch_hausdorff.py. An example usage of the application:
    milxHausdorffDistance results_clipped/bone/asm_bone_4.vtk -p Hausdorff/bone/bone__ -o Hausdorff/bone/bone__004.vtk --label manuals_renamed/bone/manual_case_bone_004.nii.gz --labelvalue 1 -c 004
    The distance map of the manual is created and the distances stored on the mesh thus giving the surface distances to the manual.
  4. The --scalarstats option in milxModelApp provides this. I did:
    milxModelApp --scalarstats Hausdorff/bladder/bladder__*.vtk -o bladder_stats.vtk
    Compute scalars stats (mean, variance etc. per point) of meshes in directory and output single mesh (of first mesh) with stats as arrays in mesh. You can then use sMILX viewer to view the meshes and change the loaded arrays via the Right Click->Show->Load Array menu.
Note that the milxDeformableModel app is an application derived from milxModelApp and additionally provides options for voxelising surfaces and applying image orientation to surfaces.

Hope this helps
Cheers Shakes - L3mming

Saturday, February 28, 2015

Compiling VTK 5.8.0 on Windows 64 with Visual Studio 2013

This is a quick post to cover the steps necessary to compile VTK 5.8.0, the most stable VTK 5 version in my experience, on Windows x64 with VS 2013 (to support Windows 8.1). The information is on the internet but scattered about so I thought I'd collect it just in case its useful for someone.

There are a few things that are incompatible with the new compiler. Firstly, the vtkOStreamWrapper error. This is covered in great detail here. Basically replace the offending line with:
//VTKOSTREAM_OPERATOR(ostream&);
vtkOStreamWrapper& vtkOStreamWrapper::operator << (ostream& a) {
    this->ostr << (void *)&a;
    return *this;
}
Then there are the ifstream->read() errors. Replace the
if ( this->IFile->read(result, 80) == 0)
with the
if ( this->IFile->read(result, 80).fail())
Lastly, there are the make_pair errors, solved here. You need to replace the
this->Map->insert(vtkstd::make_pair< vtkVariant, vtkVariant >(from, to));
with
this->Map->insert(vtkstd::make_pair(fromto)); // Addendum
in the vtkMapArrayValues.cxx file. You might need to add the include:
#ifdef _WINDOWS
  #include  // Addendum for vs11 to find 'greater'
#endif
to vtkAdjacencyMatrixToEdgeTable.cxx. I also had to add the above include and 
#include 
to vtkAdjacencyMatrixToEdgeTable.cxx, vtkNormalizeMatrixVectors.cxx, vtkPairwiseExtractHistogram2D.cxx, vtkParallelCoordinatesRepresentation.cxx, vtkChartXY.cxx, vtkControlPointsItem.cxx, to get rid of std max and min errors, 

HTHCheers Shakes - L3mming

Thursday, February 12, 2015

HP 4000b Bluetooth Mouse and Ubuntu

This quick post cover howto get the HP 4000b Bluetooth mouse working flawlessly with Ubuntu 14.10 + (and maybe even older versions). These bluetooth mice are useful for ultra-portable devices that have no or few USB ports, like my Dell XPS 13, keeping the port free for other uses.

The HP 4000b mouse is an affordable bluetooth mouse that does work out of the box with Ubuntu but drops out after use for a few minutes or so. This post covers the solution to this problem that worked for me, and this post helpful for TLP users.

In summary, add the following to your rc.local file (the one I used was /etc/rc.local):
# Prevents the Bluetooth USB card from getting reset which disconnects the mouse
BTUSB_DEV="8087:07dc"
BTUSB_BINDING="$(lsusb -d "$BTUSB_DEV" |
    cut -f 1 -d : |
    sed -e 's,Bus ,,' -e 's, Device ,/,' |
    xargs -I {} udevadm info -q path -n /dev/bus/usb/{} |
    xargs basename)"


echo "Disabling autosuspend for Bluetooth USB Soundcard: $BTUSB_BINDING..."
echo -1 > "/sys/bus/usb/devices/$BTUSB_BINDING/power/autosuspend_delay_ms"
And in the TLP config (/etc/default/tlp), add the USB_BLACKLIST="8087:07da". The ID for my device was 8087:07dc as in the rc.local file. These IDs need to match and can be found using lsusb.

HTH
Cheers Shakes - L3mming

Sunday, January 25, 2015

Visualizing 32-bit Integer PGMs (eg. from FTL) or other ASCII Image Formats

In this post I will post some Python code for visualising ASCII image formats such as PGMs (esp. from the Finite Transform Library (FTL)). These formats are commonly used because of their simplicity in implementation.

EDIT: The sMILX Viewer (v1.0 Alpha and above), part of the SMILI project for scientific visualisation now supports these PGMs. Just drag and drop the file into the viewer window to visualise them!

The ASCII PGM format for example has a simple header made up of a format string (made of two characters) such as 'P2' on the first line, followed by a comment line, the dimensions of the image (width and height) and the bit depth. Finally the data is next. An example is given below:
P2
# Generated PGM.
101 101
255
184 180 188 199 202 203 200 195 195 202 198 186 181 156 ........
Since most PGM and image viewers always expect a bit depth of 8-bit, the result isn't always shown correctly. I have written a simple Python module to do this type of reading and also utilise 32-bit PNGs as well. The full module can be found in this gist. The section for PGMs is as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def readPGM(name):
    '''
    Read a PGM image, where the PGM format is that of the FTL library for NTTs and FRTs.
    This is a slightly modified format where values are not limited to 8-bit values.
    
    Returns array of image and bit depth (image, depth)
    '''
    inFile = open(name,"r")
    
    #read header
    formatLine = inFile.readline()
    commentLine = inFile.readline()
    
    if not "P2" in formatLine:
        print "Error: PGM not in correct format (P2)"
    print "Comment:", commentLine
    
    width, height = [int(x) for x in inFile.readline().split()] # read dimensions
    print "PGM Size:", width, "x", height
    
    bitDepth = [int(x) for x in inFile.readline().split()] # read bit Depth
    
    imageList = []
    for line in inFile: # read remaining lines
        valueList = [int(x) for x in line.split()] #read integers on each line
        for value in valueList:
            imageList.append(value) #append as 1D list
#    print imageList
    #store as array
    image = np.array(imageList).reshape(height, width)
    
    return image, bitDepth

Then to load results from FTL, simply use a script as:


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# -*- coding: utf-8 -*-
"""
Test Read PGM member
Created on Sun Jan 25 21:52:38 2015

@author: shakes
"""
import imageio

image, depth = imageio.readPGM("lena.pgm")
frtSpace, depth = imageio.readPGM("frtSpace.pgm")

#Plot
import matplotlib.pyplot as plt

#plot images
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(16, 5))

plt.gray()

ax[0].imshow(image)
ax[0].axis('off')
ax[0].set_title('Image')
ax[1].imshow(frtSpace)
ax[1].axis('off')
ax[1].set_title('FRT of Lena')

plt.show()

The result is a Matplotlib plot of the images, both of Lena and her FRT space:
You can adapt it for your image format by just changing what/how the header is read. Feel free to fork the gist as I will put the latest version there. This module is part of a pure Python version of FTL coming soon, so watch this space..... better yet... Follow!

PS: A FTL plugin for my scientific visualisation software SMILI is on the way too. More soon in the next version.

Cheers Shakes - L3mming

Thursday, November 20, 2014

SMILI: A simple open-source framework for scientific visualisation

I'm pleased to announce the release of SMILI, a simple framework for scientific visualisation, under a BSD license on SourceForge and GitHub.


In this post, I thought I would cover some of the cool features and command-line utilities available for researchers and developers. The overall documentation can be found here. This post will hopefully be complimentary to this and the videos of SMILI already available on YouTube.

Some of the additional features of SMILI are the ability to share the same processing and display capabilities of the GUI applications with their command-line counterparts. For example, you have the following command-line tools to assist you in your research:

  1. milxImageViewer and milxModelViewer - These are fast, no-nonsense viewers for n-D images and 3D polygonal data, such as triangulated surfaces. They don't load pesky plugins, they just allow you to view your data fast. The right click options are still available, especially the processing elements. However, any interaction between windows is no longer possible for obvious reasons. You can use sMILX for this.
  2. milxOverlay and milxAnimate - These applications allow you to take screenshots/movies of models and images together with pre-defined views (which can be created in sMILX). This is great for visually inspecting your results or for websites when there are a lot of them. There is a batch script in 'scripts/' directory to batch these applications over multiple threads for quick generation.
  3. milxImageApp and milxModelApp - These are the 'swiss-army knife' like applications for images and models respectively. With them you can simply process all the images/models you provide via the command-line in the same way with many of the algorithms available in sMILX. For example, to threshold labelled images (having Nifti format) within a certain range and storing result in 'auto', we could simply:
    milxImageApp --threshold --above 225 --below 180 *.nii.gz -p auto/auto_
  4. milxLabelVisualisation - Sometimes it is necessary to visualise labelled images using iso-surfacing or marching cubes or volume rendering with certain colour maps. This application provides a way of doing these things with off-screen rendering.
  5. milxAssistant - This application provides a simple web browser interface to explore the SMILI documentation in a similar fashion to Qt Assistant.

As updates occur, I will post notices on twitter and Google+. Major development updates will be posted here in this blog. Currently, a journal publication for SMILI is under review entitled:
"Focused Shape Visualisation via the Simple Medical Imaging Library Interface"
Visualisation and Computer Graphics, IEEE Transactions on, 2014, Submitted
Upon acceptance, more details, plugins and revision history will be released. If you have any feedback, feel free to post as a ticket, email the mailing-list or message me on SourceForge.

Cheers Shakes - L3mming

Monday, October 6, 2014

Robust Digital Image Reconstruction Example

In this post, we discuss how to employ the digital image reconstruction technique of Chandra et al. (2014):

Robust digital image reconstruction via the discrete Fourier slice theorem
S Chandra, N Normand, A Kingston, JP Guedon, I Svalbe
IEEE Sig. Proc. Lett. (2014)

using the FTL (implemented in C, available via LGPL license).

This method takes a sufficient set of discrete (rational angle) projections assuming the Dirac pixel model, i.e. digital image sampling where lines have said to have sampled a pixel iff the line goes through the centre of the pixel, and reconstructs them in O(nlogn), where n=N^2. Sufficiency is classified as those projections meeting the Katz criterion, i.e. basically all bins are sampled at least once and no unambiguous  solution (i.e. a ghost) cannot fit within the image. See also:

Fast Mojette transform for discrete tomography
SS Chandra, N Normand, A Kingston, J Gu├ędon, I Svalbe
arXiv preprint arXiv:1006.1965

Once you have FTL built, you should have four binaries for this method:

  • fmt_angles - Select the angle type and generate the rational angles for given n and N, the image and reconstruction (FFT space) sizes respectively.
  • mt - Compute the discrete (rational angle) projections of the image, also known as the Mojette transform (MT).
  • mt2frt - Convert the projections to those of the FRT/DRT, which are the inverse FFTs of the slices of the 2D FFT.
  • ifrt - To reconstruct the resulting FRT projections in O(nlogn), no relation to the n before.

To illustrate the process we given a tutorial here of the whole process.

1. First, we crop Lena image to a 128x128 image of Lena from the centre:

./crop lena512.pgm 128 128 0 lena128.pgm


2. We create the angle set, we choose the L1 minimal set since it has a nice symmetry:
./fmt_angles 128 256 1 mt_angles_128_in_256.txt
3. Next we compute the MT:
./mt lena128.pgm mt_angles_128_in_256.txt mt_lena128.pgm
    Note that if you already have projections, such as those of a sinogram, then see this Google Groups     discussion. You can find the publication by my colleague Andrew Kingston on how to do this here.
4. Convert the MT projections into FRT space
 ./mt2frt mt_lena128.pgm mt_angles_128_in_256.txt 128 256 1.0 frt_lena128.pgm
5. Invert the FRT projections in O(nlogn) using the discrete Fourier slice theorem
./ifrt frt_lena128.pgm recon_lena128.pgm

This gives our nxn result reconstructed and padded into the NxN space.

EDIT: See my other post about visualising these and FTL (or PGM files) results in general.

HTH
Cheers Shakes - L3mming