# Image Based Calculations

## Overview

The “Image-based calculations” feature allows to evaluate free-text expressions like

```1.5 * (a - b) / c**2 + gt(a, 7.3) * log(d)
```

where `a`, `b`, `c` and `d` are variables denoting image buffers, `gt()` is a thresholding function. This feature can be used interactively (File menu → Automate) and in batch-mode (documentation pending).

A family of thresholding functions can also be used in the free-text expression: in the example above, `gt()` will first create a clone of image buffer a and then set all voxels greater than the threshold (here: 7.3) to one, all voxels below and equal to the threshold to zero. A reference of all functions we have implemented (so far) is given below.

You can use most of the mathematical standard functions e.g. `log`, `sqrt`, `sin`, etc. A reference of standard functions we have implemented (so far) is given below.

Simple arithmetic operations involving image buffers are voxel-based (similar to the conventions used in the IDL programming language), e.g. the image buffer resulting from the operation `a+b` is derived from adding corresponding voxels of the image buffers referred to by the variables `a` and `b`, resp. Multiplications with a constant factor will multiply all voxels of an image buffer with that constant, adding a constant value will add that value to all voxels. In case of divisions, if a zero-voxel is encountered in the denominator, we will assign a value of zero to the corresponding voxel in the resulting image. Exponentiation can be done with any real number, undefined values e.g. `(-1.0)**0.5` will be set to zero.

All image-based operations (so far) work on “native” voxel volume of images. If the reslicing parameters of the images are different, Vinci will refuse computation giving the error message: Images are computed in image reset position, so the involved images MUST be co-registered in image reset position! Your images have different reslicing parameters. In this case, we prohibit the computation now as it could lead to unexpected results (maybe even unperceived). If you really have to use a resliced image (e.g. co-registered by MMM) you can first call “reset link channel” on one image, then save the other image using “current offset and rotation” and re-load the saved image.

Also, we will not automatically resample “incompatible” images: only image volumes with identical geometries (dimension and pixel size) are valid in one calculation.

Image-based calculations are among Vinci’s most versatile and powerful features. Part of the flexibility is due to our VinciPy package, written in the Python programming language: it allows to remotely-control Vinci through XML-based VinciScript, also at the heart of the high-level interfaces to C/C++, PERL and IDL. The full source code is provided in the external/Vinci_Python directory. The implementation of this feature is rather compact as most of the work is left to the Python interpreter: we overload operators in an image class, interested readers should have a look at `VinciPy/Vinci_ImageCalc.py`. Clearly, adding further functions is a straight-forward task, please let us know if you have any suggestions.

## Requirements

Image-Based Calculation requires a Python interpreter (versions 2.4-2.6 will do). Python should be available on all modern Unix systems (including MacOS X) and no further installation is required. Users of MS Windows systems can get a free standard Windows installer from here: http://python.org/download (installation is just a “double-click”, the distribution’s footprint is small).

## Reference

### “Logical” Thresholding

#### `gt(a, t)`

Create a clone of imagebuffer a and set all voxels above threshold t to one, and all voxels below and equal to the treshold to zero.

#### `ge(a, t)`

Similar to `gt()` but tests for “greater than or equal”.

#### `lt(a, t)`

Create a clone of imagebuffer a and set the all voxels below threshold t to one, and all voxels above and equal to the treshold to zero.

#### `le(a, t)`

Similar to `lt()` but tests for “lesser than or equal”.

#### `eq(a, t)`

Create a clone of imagebuffer a and set all voxels equal (i.e. with relative difference less than 10**(-6)) to threshold t to one, and all other voxels to zero.

#### `neq(a, t)`

Create a clone of imagebuffer a and set all voxels equal (i.e. with relative difference less than 10**(-6)) to threshold t to zero, and all other voxels to one.

### General Thresholding

#### `th_u(a, t, r)`

Create a clone of imagebuffer a and set the “upper part”, i.e. all voxels above threshold t to r, keep all other voxels.

#### `th_ue(a, t, r)`

Similar to th_u[pper] () but tests for “greater than or equal”.

#### `th_l(a, t, r)`

Create a clone of imagebuffer a and set the “lower part”, i.e. all voxels below threshold t to r, keep all other voxels.

#### `th_le(a, t, r)`

Similar to th_l[ower] () but tests for “lesser than or equal”.

#### `th_eq(a, t, r)`

Create a clone of imagebuffer a and set all voxels equal to threshold t to r, keep all other voxels.

#### `th_neq(a, t, r)`

Create a clone of imagebuffer a and set all voxels not equal to threshold t to r, keep all other voxels.

#### `th_u0(a, t)`

:= `th_u(a, t, 0)`

#### `th_ue0(a, t)`

:= `th_ue(a, t, 0)`

#### `th_u1(a, t)`

:= `th_u(a, t, 1)`

#### `th_ue1(a, t)`

:= `th_ue(a, t, 1)`

#### `th_l0(a, t)`

:= `th_l(a, t, 0)`

#### `th_le0(a, t)`

:= `th_le(a, t, 0)`

#### `th_l1(a, t)`

:= `th_l(a, t, 1)`

#### `th_le1(a, t)`

:= `th_le(a, t, 1)`

### Mathematical Standard Functions

In the following definitions, a denotes an image buffer.

#### `pow(a, t)`

computes `a**t` with any real number exponent t. pow is defined for all voxel values greater zero, for voxels values equal zero with exponents greater equal zero, and for voxel values less than zero with integer exponents; offending result voxels will be set to zero.

#### `a**t`

:= `pow(a, t)`

#### `sqrt(a)`

defined for voxel values greater or equal zero; offending result voxels will be set to zero.

#### `log(a)`

defined for voxel values greater zero; offending result voxels will be set to zero.

#### `exp(a)`

result voxels with under- and overflow will be set to zero.

#### `sin(a)`

defined for voxel values in the range [-263,+263]; offending result voxels will be set to zero. This is a technical limit set by the programming language C++, 263 is about 83.66 π.

#### `cos(a)`

defined for voxel values in the range [-263,+263]; offending result voxels will be set to zero. This is a technical limit set by the programming language C++, 263 is about 83.66 π.

#### `tan(a)`

defined for voxel values in the range [-263,+263]; offending result voxels will be set to zero. This is a technical limit set by the programming language C++, 263 is about 83.66 π.

#### `asin(a)`

defined for voxel values in the range [-1,+1]; offending result voxels will be set to zero.

#### `acos(a)`

defined for voxel values in the range [-1,+1]; offending result voxels will be set to zero.

#### `atan(a)`

results voxels will be in the range (-π/2,+π/2).

We could easily add other functions. Please, let us now if you need them.

### Multi-Frame Files

#### `sum(a)`

Python’s built-in sum function works now as the index function and the length function are defined for images to give a specific frame and the total number of frames. BEWARE: frame numbers for the index function always start with zero, they are transported to the first “native” frame index according to the image type: e.g. “1” for ECAT7 images and “0” for HRRT images.

#### `msum(a, iStart, iEnd)`

Sums all frames from frame <iStart> up to and including frame <iEnd> of the multi-frame file associated with image buffer a (in case of ECAT7 data) or the multi-frame series (in case of HRRT data). “Native” frame numbers are used: the first frame of an ECAT7 multi-frame file is “1” (usually), the first frame of a HRRT multi-frame series is “0”.

#### `msum(a, iStart)`

Sums all frames of the multi-frame file, or the multi-frame series starting with frame <iStart>.

#### `msum(a)`

Sums all frames of the multi-frame file, or the multi-frame series.

#### `mavg(a, iStart, iEnd)`

Similar to `msum(a, iStart, iEnd)` but scales the result with 1/#frames. Please note: this function is currently only available for ECAT7 multi-frame files.

#### `mavg(a, iStart)`

Similar to `msum(a, iStart)` but scales the result with 1/#frames. Please note: this function is currently only available for ECAT7 multi-frame files.

#### `mavg(a)`

Similar to `msum(a)` but scales the result with 1/#frames. Please note: this function is currently only available for ECAT7 multi-frame files.