ImageJ

Plugin to read AIM ISQ RSQ files

ImageJ is good for Viewing, Converting, and Custom Analyses on Scanco images after installing the following plugin and importing AIM, ISQ, and RSQ files.

AIM Reader v1.02, reads uncompressed (short) AIM files v020, sets scale and calibration
ISQ Reader v1.01, reads ISQ files, sets scale
RSQ Reader v1.01, reads RSQ files, sets scale

You could Download the ImageJ_Plugins_SCANCO_Files.zip from BoydDrive:\\04 - Resources\04 - Software\ImageJ_Plugin.

Return to Contents or Return to Main Page

Generate TIFF from AIM

Refer to ImageJ SOP - generate TIFF from AIM.

Return to Contents or Return to Main Page

Amira

Script to read AIM files

Download Script to write Amira Tcl command to open AIM files.

Return to Contents or Return to Main Page

R

JAGS Statistics

Analysis of Bayesian models using Markov Chain Monte Carlo.

Overview

JAGS is a replacement for WinBUGS that can be run from within R. It is cross-platform and will run on any system (Windows, OS-X, Linux).

Installation

JAGS

You can get get JAGS from here: http://www-fis.iarc.fr/~martyn/software/jags/

R

You can get R from here: http://www.r-project.org/

runjags

You will also need the runjags package from R. This is the webpage for it: http://cran.stat.sfu.ca/web/packages/runjags/index.html . There is a very good PDF manual there, with examples. (The help is also available from within R).

Note that there is an alternate package rjags, which I haven't tried. I'm not sure which is better.

Installation of runjags is very easy. From within R, go to the menu bar and select Packages&Data -> Package Installer. Click Get List, then enter runjags in the search field. Select it and click Install Selected. You will also need to install the package coda.

Now go to Packages&Data ->Package Manager, scroll down and find coda and runjags and make sure that they are both selected. Note that if you click on runjags here will will get it's on-line help, which can be useful.

Reading a spreadsheet into R

Refer to http://cran.r-project.org/doc/manuals/R-intro.html#The-read_002etable_0028_0029-function . There is also a whole import/export manual for R at http://cran.r-project.org/doc/manuals/R-data.html

The first step is to export your spreadsheet as a text file from your spreadsheet program. You might want to copy only the relevant columns into a new spreadsheet.

If you are using Excel, I believe that you can export as a text file, although I am not certain of this.

I was using Apple's Numbers. In Numbers, select the spreadsheet or the columns that you want to Export, then start up TextEdit. Make sure that the new blank text document is plain text (Format -> Make Plain Test), then Paste in the spreadsheet. This resulted in a text file in exactly the right format for importing into R, like this

cage	fnbmd10
20	0.72
20.1	0.74
20.1	1.02
20.2	0.97
20.2	0.84
20.4	0.69
20.4	0.98
20.5	
20.6	0.82

It may happen that there are some missing data in your columns (as in the above example). This command (from the command line) will copy all lines from file fnbmd10.data to file fnbmd10_trimmed.data that have exactly two fields, which is what we want:

awk '(NF==2) {print $0}' fnbmd10.data > fnbmd10_trimmed.data

If you have a column for sex, and you want to split it into two files, one for each sex, you can do this:

awk '($2==0) {print $0}' input.txt > output_w.txt
awk '($2==1) {print $0}' input.txt > output_m.txt

In this example, the sex is the second field (hence $2). You will need to manually add the column headings again after doing this.

Once your text file looks correct, you can load it into R with the following command (from within R obviously):

dataframe <- read.table("fnbmd10_trimmed.data", header=TRUE)

I recommend at this point using some of R's graphical capabilities to plot and examine the data, just to make sure it seems to be reasonable.

plot (dataframe$cage, dataframe$fnbmd10)

Tip: Setting the Working Directory in R

setwd("/Users/eric/Documents/change_point")

Tip: Converting line endings

In case it's required, this can be done from a command prompt in OS-X or Linux:

cat input.txt | tr '\r' '\n' > output.txt

Model file

Here is the JAGS model I used, findcp.model . I suggest thinking carefully about what distributions, with what values, are most appropriate for b.intercept, b.slope and tau as this will affect how the solution converges.

Obtain slopes for females #
25 to 45 for the BMD value. #
Determine the cut point from 25 to 45 . #

var

 n,           # Read from data file; number of subjects
 age[n],      # Read from data file; age of subject
 y[n],        # Read from data file; BMD of subject
 b.cutoff,    # change point age
 b.intercept, # constant value of BMD function
 b.slope,     # slope of BMD function
 mu[n],       # expected mean BMD of subject
 tau,         # tau of the normal distribution of BMD
 sigma;       # std. dev. of the normal distribution of BMD

model {

   for (i in 1:n) {
       mu[i] <- b.intercept + b.slope * step(age[i]-b.cutoff) * (age[i]-b.cutoff)
       y[i] ~ dnorm(mu[i],tau)
   }
   b.cutoff ~ dunif(25,55)
   b.intercept ~ dnorm(.9,6)
   b.slope ~ dnorm(0,9)
   tau ~ dgamma(1.0E-3,1.0E-3)
   sigma <- 1.0/sqrt(tau)

} </source>

Tip: Using a BUGS style model file instead of a JAGS style model file

JAGS and BUGS model files are very similar, however, the BUGS model file allows more type of things in the file, such as initialization values, that are not allowed in a JAGS model file (in JAGS these go in a separate file).

runjags can read BUGS-style files directly, without modification, with the R command autorun.jagsfile. This is instead of autorun.jags. See Below.

Running a model in R

This is pretty much just following the example from the runjags documentation for autorun.jags .

We determine the number of data points and stick our data in an R list format:

n <- length(dataframe$cage)
jagsdata <- list (age=dataframe$cage, y=dataframe$fnbmd10, n=n)

We read in the model as a test string. This worked for me, but I suspect that there is an easier way.

modelraw <- readBin("findcp.model", raw(), 9000)
model <- rawToChar(modelraw)

The last argument to readBin is the number of bytes to read. Just make sure that this is bigger than the file size and you will get the whole file.

Now we run the model:

results <- autorun.jags(model=model, monitor=c("b.cutoff", "b.slope", "b.intercept"), data=jagsdata)

You can analyze the traceplots of the results to assess convergence:

results$trace

To get a summary of the monitored variables:

results$summary

This will look something like this:

<source lang=text> Iterations = 1:21741 Thinning interval = 2 Number of chains = 2 Sample size per chain = 10871

1. Empirical mean and standard deviation for each variable,

  plus standard error of the mean:

                Mean        SD  Naive SE Time-series SE

b.cutoff 43.152034 5.9029519 4.003e-02 1.799e-01 b.slope -0.003935 0.0006817 4.624e-06 1.866e-05 b.intercept 0.832309 0.0102962 6.983e-05 2.305e-04

2. Quantiles for each variable:

                2.5%       25%       50%       75%     97.5%

b.cutoff 31.641838 39.016150 42.733250 47.586850 53.975075 b.slope -0.005512 -0.004329 -0.003827 -0.003447 -0.002851 b.intercept 0.813198 0.825117 0.832013 0.839125 0.853183 </source>

Summary of places to find Documentation

The OpenBUGS manual: http://www.openbugs.info/Manuals/Manual.html .
The JAGS manual : http://www-fis.iarc.fr/~martyn/software/jags/jags_user_manual.pdf .
The manual for runjags : http://cran.stat.sfu.ca/web/packages/runjags/runjags.pdf . Note that runjags help is also available on-line in R. For example, entering help(autorun.jags) will bring up the help on that command.

OBSOLETE: JAGS run script

This is not required when running JAGS from within R, which seems to be the best way to run JAGS.

Here is the script I used to run JAGS, fnbmd10.jags. Of course, you could enter all these commands interactively into JAGS as well. You probably want to customize this. For example, I just guessed at the number of iterations. There are also several output options. Probably you want to use the coda command to save the monitor values in files that can be read by the coda package of R, and then do some sort of statistical analysis on the convergence.

Read in the model and the data

model in findcp.model data in fnbmd10.R

Required stuff

compile initialize

Do some iterations
(The adaptive part of the iterating needs to be done before monitors can be set)

update 500

Now set some monitors to track parameters over iterations

monitor b.cutoff, thin(50) monitor b.intercept, thin(50) monitor b.slope, thin(50)

Do lots of iterations

update 10000

Save parameters and monitors.
There are other kinds of outputs and files you can write; refer to the
JAGS documentation

parameters to parameters.out monitors to monitors.out

Bye

exit </source>

Return to Contents or Return to Main Page

ParaView

Getting ParaView With Numerics88 Extensions

A custom-compiled version of ParaView is available with extensions to open AIM files, and FAIM output files directly. This version is already installed on many lab computers. If you need a copy, go see Eric.

How-Tos

Display Slices of an AIM file

Open the AIM file. Make sure that Data on Cells is not selected. Click the Apply button in the Object Inspector, as shown.

Go to the Display tab of the Object Inspector, find the Style section, and set the Representation box to Slice.

Now you probably want to set the viewing direction along a particular axis. There are buttons on the toolbar for this purpose. Finally, for viewing slices, it is convenient to use parallel projection. Select View Settings from the Edit menu, and make sure that Use Parallel Projection is selected.

Now you can scroll through the slices with the Slice slider bar, which is located on the Display tab of the Object Inspector. You may want to click the Toggle Color Legend Visibility button on the toolbar in order to get a legend for the color scale.

Display an AIM file using the Threshold filter

Open the AIM file. Make sure that Data on Cells is selected, otherwise, the data will be interpolated in the threshold filter and this can result in artefacts (for example, it can make it appear that disconnected regions are connected, or vice versa). Click the Apply button.

Making sure that the AIM file is selected in the Pipeline Browser (it should be if you just opened it), choose a Threshold filter. You can do this from the icon on the toolbar, or Filters menu. Select suitable threshold values with the sliders. Click Apply. If your data file is large, you might not see anything yet, since ParaView may have chosen Outline as the best viewing mode. You can change this to Surface. This can either be done from the toolbar, where there is a drop-down list, or you can go to the Display tab of the Object Inspector.

Tips:

You can use multiple Threshold filters on the same data. This can be useful if your data consists of two different materials. In this case you can go to Display and set a different Solid Color for each Threshold filter.

Display an AIM file using the Contour filter

Contouring can provide smoother-looking surfaces than Thresholding. However, it is more computationally intensive.

Open the AIM file. Make sure that Data on Cells is not selected, but that Pad is. The reason for padding the image, which adds a layer of zero-valued points all around the image volume, is that this will result in outer surfaces being shown. Click Apply. Select a Contour filter, either from the toolbar or the Filter menu, and click Apply. Done!

Show exaggerated displacements of an FE model

Open the FE .out file. Make sure it is selected in the Pipeline Browser.
Go to the Filters menu and choose Common -> Warp By Vector
Make sure that the drop down box Vectors in Properties for the Warp By Vector filter has NODAL_DISPLACEMENTS selected. (It should already.)
Enter the desired Scale Factor and click Apply.

Show displacements only in the x,y directions

The example shows how to use the calculator to manipulate vector fields associated with data objects. Here we are eliminating the z component of the displacement. This may be useful when we want to investigate fine details of the x,y displacement without flattening the model due to excessive exaggerated z displacement.

Starting with an open FE model, select a Calculator filter.
Enter this expression, which subtracts the z-component from the vector field NODAL_DISPLACEMENTS:

NODAL_DISPLACEMENTS -  kHat * (kHat.NODAL_DISPLACEMENTS)

Click Apply.
Now apply a Warp By Vector filter as described previously.

Show the difference in displacements between two similar models

Method 1: Using a Programmable Filter

This example will demonstrate several things:

Using a filter with multiple inputs.
Using a programmable filter.
Using numpy arrays within a ParaView Programmable Filter. (See [1].)

Here are the steps:

Open both data files. In this example, one will be a FAIM output file, while the other will be in the form of a VTK Unstructured Grid file.
Select both files at once in the Pipeline Browser. This is done in OS-X by clicking on one, then holding down the Command key and clicking on the other.
Now select a Programmable Filter from the menu.
Enter this script:

<source lang=python> from paraview.vtk import dataset_adapter as DA

Get inputs/outputs as VTK objects

input0 = self.GetInputDataObject(0,0) input1 = self.GetInputDataObject(0,1) output = self.GetOutputDataObject(0)

We can get the array data from the "inputs" object, which is an alternative interface (more Pythonesque) to the same data.
This gives us numpy arrays.

x0 = inputs[0].PointData.GetArray("NODAL_DISPLACEMENTS")

Change the following to GetArray("NODAL_DISPLACEMENTS") if your second input is also a FAIM data file.

x1 = inputs[1].PointData.GetArray("displacement")

Numpy arrays can be operated on element-by-element with simple expressions. (Both x0 and x1 are Nx3 arrays.)
You could also do a loop over the elements, but it would be slower and more verbose.

xdiff = x1 - x0

Convert from a numpy array to a VTK array.
Call it "dispdiff", but you could name it anything.

xdiff_vtk = DA.numpyTovtkDataArray(xdiff, 'dispdiff')

The output is a copy of the first input, except we attach our new vector data to the points.

output.ShallowCopy(input0) output.GetPointData().AddArray(xdiff_vtk) </source>

Click OK
Change the Display of the programmable filter to "dispdiff", and click "Toggle Color Legend Visibility" button.

Method 2: Using a Python Calculator Filter

This is a little simpler than using a Programmable Filter. Instead of a Programmable Filer, choose a Python Calculator. In the expression field enter this line: <source lang=python> inputs[0].PointData.GetArray("NODAL_DISPLACEMENTS") - inputs[1].PointData.GetArray("displacement") </source> Set the Array Name if you want, then hit Apply.

Show differences between two images prior and after registration

Image differences

Load two images
Select a subregion with the same size in both images
Select two subregions (Ctrl-Click) and feed them to a Python Calculator. In the expression field enter this line:

<source lang=python> inputs[0].CellData.GetArray("AIMData") - inputs[1].CellData.GetArray("AIMData") </source> This assumes that the AIM image attaches pixel values to CellData.

Set 'Array Association = CellData'
Uncheck 'Copy Arrays'
Set the Array Name if you want, then hit Apply.
To visualize, choose Display->Style->Representation = Slice. Play with XY vs XZ vs YZ views and select a slice of interest. Modify the colormap and mapping of the grayscale values to have a better comparison.
Add a second '3D View' in a separate panel and repeat the above for the second image pair.
Add the image comparison as a 'custom filter' by selecting the python calculator you just created and choosing Tools->Custom Filter, give a name like 'Compare Images' and click Finish. Select any two images you want to compare and Choose Filter->Alphabetical->Compare Images - done :o)

Here is the end result:

Checkerboard

Load two images, unselect 'DataOnCells' so that pixal values are attached to PointData
Select a subregion with the same size in both images (Somehow we need 'subregion', otherwise it is not happy later)
Select two subregions (Ctrl-Click) and feed them to a Programmable Filter, the code is:

<source lang=python> from paraview.vtk import imaging as im

Get inputs/outputs as VTK objects

input0 = self.GetInputDataObject(0,0) input1 = self.GetInputDataObject(0,1) output = self.GetOutputDataObject(0)

Create a checkerboard object

checker = im.vtkImageCheckerboard()

Set checkerboard inputs

checker.SetInput1(input0) checker.SetInput2(input1)

Set the number of tiles in each dimension

checker.SetNumberOfDivisions(4,4,2)

Copy to the output object

checker.Update() output.DeepCopy(checker.GetOutput()) </source>

Click apply and use slice for visualization.
Create a custom filter as described above. Note I don't know how you would expose the SetNumberOfDivisions(), that would be handy.

Here is the end result:

Wish List

This is a list of things that we would like to be able to do. They may be things that ParaView can inherently do already, but we haven't figured out how, or they may be things that will require the writing of a new plug in.

Save a subvolume of an AIM

Would like to be able to read an AIM, select a reduced volume, and write that out as another AIM.

Create an FE model interactively.

This would involve lots of smaller steps. For example, selecting sets of points for Boundary Conditions.

Registration within ParaView

There are some existing plugins. They might be useful, or posibly they could be further customized for our purposes. http://paraview.org/Wiki/User_Created_Plugins

World Peace button on the toolbar

Return to Contents or Return to Main Page

Writing Plugins for Paraview

Have a look at the paraview_plugins project in the Numerics88 source code. The plugins there are working examples.

Links

Here are some useful links.

http://paraview.org/Wiki/Plugin_HowTo

http://paraview.org/Wiki/ParaView:FAQ#How_do_I_add_my_own_vtk_filter.2Freader_to_ParaView.3F

http://public.kitware.com/pipermail/paraview/2003-September/000290.html

http://paraview.org/Wiki/ExtendingParaView

Tips

Link to Paraview's internal VTK

Your plugins, and any libraries that they use, must link to Paraview's internal version of VTK, and not some other VTK libraries. This is true even if both the internal and external versions of VTK were built from the same VTK sources! As a concrete example, paraview_plugins links to vtkski. The plugins will die with mysterious segmentation faults unless you build vtkski to link against the Paraview internal version of VTK. So for example, on my Linux system, cmake will automatically find vtk at /usr/lib/vtk-5.2 . However, I have to change the cmake variable VTK_DIR manually to /home/eric/source/ParaView3/build/VTK .

Implement RequestInformation for vtkImageData data types

This seems to be necessary for vtkImageData, although strangely not for vtkUnstructuredData.

Here is an example, where the extent, spacing, origin and data type are set in RequestInformation. This seems to be sufficient.

<source lang=c> int vtkAIMReaderPlugin::RequestInformation (

 vtkInformation       * vtkNotUsed( request ),
 vtkInformationVector** vtkNotUsed( inputVector ),
 vtkInformationVector * outputVector)

{

 // Make sure FileName is set
 if (!this->FileName)
 {
   vtkErrorMacro("A Filename must be specified");
   return 0;
 }
 vtkSmartPointer<vtkAIMReader> reader = vtkSmartPointer<vtkAIMReader>::New();
 reader->SetFileName (FileName);
 reader->UpdateInformation();
 vtkImageData* data = reader->GetOutput();
 vtkInformation* outInfo = outputVector->GetInformationObject(0);

 outInfo->Set (vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
               data->GetWholeExtent(), 6);
 outInfo->Set (vtkDataObject::SPACING(), data->GetSpacing(), 3);
 outInfo->Set (vtkDataObject::ORIGIN(), data->GetOrigin(), 3);
 vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_FLOAT, 1);
 return 1;

} </source>

You may need to convert to VTK_FLOAT

Paraview does not seem to handle byte (char) data well. In particular, there is sign/unsigned ambiguity with char types in paraview. If your source produces byte data, I recommend converting it in your plugin to float.

Scalar values must be named

If you don't assign a name to scalar values, paraview can't do anything with them.

Here is an example.

<source lang=c> int vtkAIMReaderPlugin::RequestData (vtkInformation*,

                                    vtkInformationVector**,
                                    vtkInformationVector* outputVector)

{

 vtkSmartPointer<vtkImageData> temp = vtkSmartPointer<vtkImageData>::New();

 // ACTUAL WORKING CODE TO FILL temp REMOVED FOR CLARITY
 // see vtkAIMReaderPlugin.cxx for details.
 //   .
 //   .
 //   .

 vtkImageData* out = vtkImageData::GetData (outputVector);
 out->ShallowCopy(temp);

 // Must have a name for the scalar data otherwise paraview can do nothing
 // with it.
 out->GetPointData()->GetScalars()->SetName("AIMData");

 return 1;

} </source>

Copy your data to the output, instead of connecting the output to a VTK filter

I don't understand why this is sometimes required, but using ShallowCopy seems be necessary in certain cases. Note there should in theory be no performance impact to using ShallowCopy (no new copy of the data will be made in memory).

Here is an example.

<source lang="c"> int vtkAIMReaderPlugin::RequestData (vtkInformation*,

                                    vtkInformationVector**,
                                    vtkInformationVector* outputVector)

{

 // Make sure FileName is set
 if (!this->FileName)
   {
   vtkErrorMacro("A Filename must be specified");
   return 0;
   }

 vtkSmartPointer<vtkAIMReader> reader = vtkSmartPointer<vtkAIMReader>::New();
 reader->SetFileName (FileName);
 reader->UpdateInformation();   // Not sure if this is redundant
 reader->Update();
 vtkImageData* data = reader->GetOutput();

 vtkImageData* out = vtkImageData::GetData (outputVector);
 if (!out)
 {
   return 0;
 }
 // Not sure why we need a temporary object for the caster, but without it,
 // an error occurs in paraview.
 vtkSmartPointer<vtkImageData> temp = vtkSmartPointer<vtkImageData>::New();
 // Convert to float, regardless of original type
 vtkSmartPointer<vtkImageCast> caster = vtkSmartPointer<vtkImageCast>::New();
 caster->SetInput (data);
 caster->SetOutput (temp);
 caster->SetOutputScalarType (VTK_FLOAT);
 caster->Update();

 out->ShallowCopy(temp);

 // Must have a name for the scalar data otherwise paraview can do nothing
 // with it.
 out->GetPointData()->GetScalars()->SetName("AIMData");

 return 1;

} </source>

This is the error that you will get if you fail to use ShallowCopy and just set the output of caster in the above example to 'out':

ERROR: In /home/eric/source/paraview-3.6.1/Servers/Common/vtkPVDataInformation.cxx, line 586
vtkPVDataInformation (0x2c391a0): Could not cast object to a known data set: vtkAlgorithmOutput

Return to Contents or Return to Main Page

Software Plugins and Codes

ImageJ

Plugin to read AIM ISQ RSQ files

Generate TIFF from AIM

Amira

Script to read AIM files

R

JAGS Statistics

Overview

Installation

JAGS

R

runjags

Reading a spreadsheet into R

Tip: Setting the Working Directory in R

Tip: Converting line endings

Model file

Tip: Using a BUGS style model file instead of a JAGS style model file

Running a model in R

Summary of places to find Documentation

OBSOLETE: JAGS run script

ParaView

Getting ParaView With Numerics88 Extensions

How-Tos

Display Slices of an AIM file

Display an AIM file using the Threshold filter

Display an AIM file using the Contour filter

Show exaggerated displacements of an FE model

Show displacements only in the x,y directions

Show the difference in displacements between two similar models

Method 1: Using a Programmable Filter

Method 2: Using a Python Calculator Filter

Show differences between two images prior and after registration

Image differences

Checkerboard

Wish List

Writing Plugins for Paraview

Links

Tips

Link to Paraview's internal VTK

Implement RequestInformation for vtkImageData data types

You may need to convert to VTK_FLOAT

Scalar values must be named

Copy your data to the output, instead of connecting the output to a VTK filter

Navigation menu

Search