How To's

HowTo Library

In addition to this HowTo library, much information is also available at the Scanco Medical FAQ.

• users@scanco.ch
• basser5dorf

Evaluation

HowTo: Understand the morphological parameters

BV/TV  -> VOX-BV/TV
BS/BV  -> TRI-BS/BV
Tb.Th  -> DT-Tb.Th
Tb.Sp  -> DT-Tb.Sp
Tb.N   -> DT-Tb.N
ConnD  -> Conn-Dens
SMI    -> TRI-SMI
DA     -> TRI-DA

Usually I'd use the names on the left for the plots, but these correspond to the names on the right of the arrow in the Scanco output files.

Note, if you use the TRI-Tb.Th rather than the DT-Tb.Th*, then you are using the derived (plate-model) morphological parameters. These are less accurate than the direct parameters. They are calculated by triangulating the surface (marching cubes) to get BS, and then calculate thickness Th = 2 * BV / BS, etc...

Mean Values

mean 1: mean of bone AND marrow --> 'apparent density'
mean 2: mean of what was considered bone only (i.e. what was segmented) 
'bone tissue density'

Note, there are five possible categories for the mean Mu. They are assigned by setting -region_number in the /voxgobj_scanco_param function. The region number is zero-based, so setting it to 0 corresponds to mean1, 1 corresponds to mean2, etc. So, although mean1 and mean2 are typically apparent and tissue density, respectively, this can be adjusted in the IPL script by setting the region_number. Mean3 might be useful for storing cortical bone density, for example.

/voxgobj_scanco_param      Voxel or Mu Scanco Parameters
  -input                     Gray scale or Segmented AIM object
  -gobj_filename             GOBJ: input or Automatic
  -peel_iter                 Peel Iterations or Automatic
  -region_number             Region number for Mu only

! -region_number             0   ! mean1
! -region_number             1   ! mean2
! -region_number             2   ! mean3
! -region_number             3   ! mean4
! -region_number             4   ! mean5

What you should also know: because of partial voluming effects, for Mean2, we are DISCARDING the 2 surface voxels of EVERY trabeculae - thus, to get meaningful Mean2, you have to do your scan with the chosen resolution high enough, such that the mean trabeculae is depicted by at the very least 5 pixels, better actually more than 6 pixels thickness. (You can do this check at the end in the XLS sheet: divide DT-Tb.Th* by el-size-mm.x, and you need to get numbers bigger than 6 or 7 or 10 or so.

For parameters such as ConnD where no GOBJ is taken as input, the TV is determined from the processing log stored in an AIM that has previously had a GOBJ applied (stored as the relative region of interest; ROI). So, it is important to apply a GOBJ prior to ConnD measurements to ensure the TV is correctly determined (and not taken as the whole AIM).

Trabecular Number (TbN) The epsilon values can be set to account for differences in scanner resolution, and hence 'surface roughness' for individual pixel suppression.

For microct, for good resolution: if there is a pixel 'sticking out' from a flat surface we do not assign that pixel with a thickness value of 1, but we let it 'get eaten' by the bigger sphere that sits in the center of the trabeculae. So this is nice for 'noise' (0.9/1.5).

For XT resolution, where also many 'good' features are just one or two pixels wide, I liked a lower epsilon value better (0.9/0.5).

HowTo: Include cortical thickness in the standard evaluation

Use the IPL command /dt_cortex_obj_param to measure cortical thickness using direct methods.

When extracting the results from the database using uct_list, specify the UCT_3D_LIST_W_CTH.TXT template in directory DISK0:[MICROCT_COMMMON.UIL].

$ uct_list
UCT_LIST.EXE Version V5.1
Template-File [UCT_3D_LIST_SHEET]: DISK0:[MICROCT_COMMON.UIL]UCT_3D_LIST_W_CTH.TXT

UCT_3D_LIST_W_CTH.TXT is now the default template as of March 2009.

HowTo: Evaluate a measurement multiple times and write out the results to a file

The basic procedure is to 1) clear the database of existing results, 2) do your evaluation, 3) extract the results, 4) paste results into excel template. Repeat as necessary. You must submit all jobs to a single queue!!!

1) Start ipl and clear the database of existing results.

$ IPL_BATCH

/db_clear_eval
-input "IPL_AIM

..
$

2) Insert your IPL evaluation.

3) Extract the results.

$@com:extract_results.com P1 P2 where P1 is the aim file (including path) and P2 is the result file (including path)

Example: $@com:extract_results.com DISK5:[BONELAB.PROJECTS.BUIE]MYAIM.AIM DISK5:[BONELAB.PROJECTS.BUIE]MYRESULTS.TXT

Note: P2 should be the same for all samples in your project so that results are compiled into a single file.

4) Paste results into excel template.

Download the template here: Excel Extract Results Template

Open your results text file and copy and paste the contents into cell A2. There's no need to use paste special since the sheet is protected.

HowTo: Run the Scanco cut and match function (on Percy)

DATA DISK: 
$ upat_list
mode 1, list mode.
Select patients and site, etc.
Select the appropriate template file, e.g. UE:UPAT_TEMPL_CUTMATCH_AUTO_PATTB_CT.TXT
Output to a file, then run that file.

PROJECT DISK: 
If you have moved the files to the project directory and renamed them, you can generate a script from:
$ @UE:SEARCH_PROJECTS_CUTMATCH.COM
-P1 = filename with device and directory
-P2 = site (RR, RL, TR, TL)
-P3 = output filename

For example, if we run:
$ @UE:SEARCH_PROJECTS_CUTMATCH.COM DISK3:[MICROCT.PROJECTS.GRACE]GRACE_01_00MTH_RR.AIM RR MYFILE.COM

The script will search DISK3:[MICROCT.PROJECTS.GRACE] for files in the form G*RR*.AIM with the same
length as GRACE_01_00MTH_RR, then check sample numbers to determine if the file is a follow-up scan.
Upat_list is called to obtain position and dimension data for matching slices.  The output file 
DISK3:[MICROCT.PROJECTS.GRACE]MYFILE.COM is created to submit and run the cut and match.  

Before running, open the file and add $ exit at the end.
$ @DISK3:[MICROCT.PROJECTS.GRACE]MYFILE.COM


For multiple files, do something like the following.
$ @UE:SEARCH_PROJECTS_CUTMATCH.COM DISK3:[MICROCT.PROJECTS.GRACE]GRACE_01_00MTH_RR.AIM RR MYFILE.COM
$ @UE:SEARCH_PROJECTS_CUTMATCH.COM DISK3:[MICROCT.PROJECTS.GRACE]GRACE_02_00MTH_RR.AIM RR MYFILE.COM
Open MYFILE.COM and add $ exit at the end.
$ @MYFILE.COM

NOTES:  The script assumes all follow-ups are provided.  If follow-up scans are missing, it will be 
necessary to modify the output file so that the position and dimension data matches the files or 
update the SEARCH_PROJECTS_CUTMATCH.COM script.

HowTo: Use UCT_EVALUATION tool version 6

$ run um:uct_k_v6

The new evaluation tool v6 should be used in place of the older v5. Please see the Documents section of this wiki for a manual on its use.

Known bugs:

When you are producing your own EVAL scripts, and want to call your own ipl file or such:

The 'EVAL:' keyword in the scripts has the filenames of the uct_evaluation_v6.com, then the ipl_v6_xxxx.com, then additional filenames (such as the _prs.com to print and the 3dresults_append.com)

Now, if you are in BASIC mode in the modify window, and you change any of the 3 filename fields, the last of the additional filenames is 'dropped'... ooops...

So please go to Advanced mode and check the long 'EVAL:' keyword that all wished files are still in the list, otherwise add the last again (most of the time, that's the useful _3dresult_append.com script...)

Will correct that of course in the next version.

HowTo: Start a 3D evaulation from the terminal

Warning: Do not change uct_evaulation_v4.com and ipl_fixed_batch_v4.com in .eval, but only work on your own file (i.e., eval_steve). First, an explanation of how the normal evaluation runs by hitting the 'Start 3d' button.

$ scanco_submit_command ue:uct_evaluation_v4.com/para=(xxx,yyy,zzz,…)

There are 8 parameters:

• P1: number of print outs of our results sheet
• P2: the ISQ filename
• P3-5: start x (pos x), y, z
• P6-8: dim x, y, z

The parameters p3-p8 come from the white bounding box. The problem is that when you start from the terminal you do not have these. This problem can be overcome by using the isq_to_aim_gobjbox.exe and then give the relevant parameters. You can $ run um: isq_to_aim_gobjbox.exe to see what it's asking for. By using gobjbox you do not need p3-p8. Therefore you can start the evaluation by giving:

$ scanco_submit_command ue:uct_evaluation_special_josh.com –
/para=(1,dk0:[microct.data.0000275.0000245]c0000332.ISQ)

Note that you can make uct_evaluation_special_josh.com such that there is no need for P1 (change p2 to p1 in the script, etc). [Alternatively to using the gobjbox you could use the normal isq_to_aim, but hard code the dim and pos.] Then call whatever IPL files you want to call. Look into fixed_batch_v4 to see how we turn on the database always on the first line of every IPL part, and how we calculate parameters.

HowTo: Associate the evaluation database with an evaluation

When /db_scanco_activate is used in the IPL files, the correct association with the evaluation database is through the Index Patient and Index Measurement in the AIM header.

Index Patient                                    3095
Index Measurement                                5962
!-------------------------------------------------------------------------------

Using IPL

$ ipl
$ read in xxx.aim
$ examine in log
$ db true           {activates database}
$ dt_object         {at end "writing results to database for measurement 1234"}

HowTo: Submit a job to the queue

$ submit/noprint /que=sys$slow ipl_voi_josh_m.com 
  - /log=sys$scratch:test.log

HowTo: Extract patient results on a slice-by-slice basis

There is an executable that was specially installed by Andres. Run it as user microct:

$ RUN DISK1:[CALIBRATION.UTIL]UPAT_EXPORT_MEASUREMENT.EXE

This gives you the .txt file in scratch containing the morphological parameters determined for each slice.

HowTo: Extract your results from the 3D database

Assuming you've used the SCANCO_SUBMIT_COMMAND to initialize your evaluations (either manually or through the uct_evaluation tool) then you can get all of your results stored into a text file that can be easily imported in Excel.

$ uct_list

The output text file will be stored in sys$scratch and you can ftp this to a PC so that it can be imported into excel and further analysis can be performed.

Importing into Excel requires setting a deliminator based on the space.

Customizing the uct_list output.

The default list template is set by the logical UCT_3D_LIST_SHEET which is normally set to point to UU:UCT_3D_LIST.TXT. See the templates by

$ dir uu:*.txt

Almost all the available keywords are in UCT_3D_LIST.TXT enclosed by \\ and //.

Only when you use the 'show directly on screen mode':

$ uct_list 1234

then you need to do this to override the default output:

$ def/user uct_3d_single_sheet  my-template_file.txt
$ uct_list 1234

Some example format files are:

UCT_3D_LIST.TXT
UCT_3D_SINGLE.TXT
UCT_3D_LIST_W_CRH.TXT
UCT_3D_MEANS.TXT       --> Use this to get the 5 mean values (i.e., calibration)

A summary of all the available keywords is available in UU:UCT_3D_SINGLE_SHEET.TXT

Name             :\\UCTMEASFILENAME//
Sample-Name      :\\SAMPLENAME//
Sample-Nr        :\\SAMPLENUMBER//
Meas-Nr          :\\UCTMEASNR//
Date of Birth    :\\SAMPLEDOB//
Sample-Remarks	 :\\SAMPLEREM//
Meas-Remarks     :\\UCTMEASREM//
Meas-NrofSlices  :\\UCTMEASNUMSLI//
Direct: -------------------------------------------------------------------
TV        : \\VOX_TV//
BV        : \\VOX_BV//
BV/TV     : \\VOX_BV_TV// 
Conn.Dens.: \\UCT_CONN//
SMI       : \\MC_SMI//
Tb.N*     : \\DT_N//
Tb.Th*    : \\DT_TH//
Tb.Sp*    : \\DT_SP//
Tb.1/N.SD*: \\DT_N_SD//
Tb.Th.SD* : \\DT_TH_SD//
Tb.Sp.SD* : \\DT_SP_SD//
Lin.Att.1 : \\UCT_MU0// +-\\UCT_MU0_SD//
Lin.Att.2 : \\UCT_MU1// +-\\UCT_MU1_SD//
Lin.Att.3 : \\UCT_MU2// +-\\UCT_MU2_SD//
Lin.Att.4 : \\UCT_MU3// +-\\UCT_MU3_SD//
Lin.Att.5 : \\UCT_MU4// +-\\UCT_MU4_SD//
TRI (Plate Model): --------------------------------------------------------
TV        : \\MC_TV//
BV        : \\MC_BV//
BS        : \\MC_BS//
BV/TV     : \\MC_BV_TV//
BS/BV     : \\MC_BS_BV//
Tb.Th     : \\MC_TH//
Tb.Sp     : \\MC_SP//
Tb.N      : \\MC_N//
Anisotropy:  --------------------------------------------------------------
|H1|      : \\MC_H1//  H1x: \\MC_H1X//  H1y: \\MC_H1Y//  H1z: \\MC_H1Z//
|H2|      : \\MC_H2//  H2x: \\MC_H2X//  H2y: \\MC_H2Y//  H2z: \\MC_H2Z//
|H3|      : \\MC_H3//  H3x: \\MC_H3X//  H3y: \\MC_H3Y//  H3z: \\MC_H3Z//
DA        : \\MC_DA//
Dim X	  : \\VOX_DX//
Dim Y	  : \\VOX_DY//
Dim Z	  : \\VOX_DZ//
Pos X	  : \\VOX_PX//
Pos Y	  : \\VOX_PY//
Pos Z	  : \\VOX_PZ//
El-Siz X  : \\VOX_ELSX//
El-Siz Y  : \\VOX_ELSY//
El-Siz Z  : \\VOX_ELSZ//  

HowTo: Tabulate your LOG file results from 3D evaluations

Only use this if you have not used the standard uct_evaluation tool to generate your 3D evaluations (in which case see the HowTo related to the 3D database). You can tabulate the LOG files:

$ search xxxx.log "ov "

Or whatever you want to search for can be stored in a text file

$ search/output=sys$scratch:ov.txt xxxx.log "ov "

To search more than one log file, and get the output into one single file:

$ fda/output=xxx.txt search *.log "ov "

The fda means 'file do all' and is a script made by Scanco to do one command (search) on a wildcard of files (*.log). Then put all the ov.txt, tv.txt, name.txt etc. files into excel. Josh has made a .com file that does a similar job as described above.

HowTo: Compile the evaluation logs into one file

Call the UE:UCT_EVALUATION_V4_3D_RESULT.COM with the same construction as $ eval_special. If you need help with this, see the HowTo on this Wiki site (HowTo: Assign special 3D evaluation to the UCT_EVALUATION tool).

Edit UE:UCT_EVALUATION_V4_3D_RESULT.COM file to call a specific ipl script as you wish.

When you run the evaluation, you will have three new files created:

• Two files in the measurement directory (on DK0!) (one with each eval 3D results in its own file, and one APPENDED_ALL.TXT file)
• Then in sys$scratch, you will have one SYS$SCRATCH:UCT_LIST_ALL_3D_RESULTS.LOG file with ALL 3D evaluations of all measurements.

Note: You may need to edit the SCANCO_SUBMIT_COMMAND to specify a single queue (i.e., /que=sceval$fast)

HowTo: Do cortical midshaft evaluation with moment of inertia calculation

This explains what happens when you type

 
$ midshaft       ! or
$ eval_midshaft 

Draw a rough contour around your slices: open the Contouring Sub-window, double-click on the big image, and it automatically finds the best contour. Then click on 'Selection' 'Forwards', then 'Iterate Forwards'.

Then 'cancel' the contouring sub-window, and start the 3D evaluation. Make sure you have already set the correct threshold, sigma, support in file

DISK2:[MICROCT.EVAL]IPL_FIXED_BATCH_MIDSHAFT.COM

When you type

$ show queue

you should find a job called UCT_EVALUATION_MIDSHAFT and a SYS$SCRATCH:*_MIDSHAFT_*.LOG

The IPL_FIXED_BATCH_MIDSHAFT.COM by default does not calculate the direct measurement results because you must have many more than 20 slices for these to be valid. You can uncomment these parts of the file if you contoured the whole stack (i.e., stack thickness is >> than cortex thickness).

TILDE: Note that you will see lots of tildes (~) on the result sheets. This is to tell you that, if the measurements were trabecular bone, those values would only be approximate. But now for the cortex, they are actually correct! (see below for more explanation about this) I would now use the tildes to check the 'excel' table that the operator has taken the correct script for the trab. or cort. evaluations for the respective measurements!

So what's the advantage of all this? Now, also VOX-TV is correct! So VOX-TV is the whole area, VOX-BV the doughnut, so TV - BV is the marrow area! And still, TRI-Tb.Th is the cortical thickness (TRI-TV is NOT correct, because the .gobj feature was turned off).

I say 'area' above: depending on how many slices one individual scan has, VOX-BV (and VOX-TV) might be different for different scans, so you will HAVE to convert the volumes to area again:

volume equals area times height
--> Bone Area = VOX-BV / height = BV / (dim.z * el_size_mm.z)

and same for Total area, of course :-)

POROSITY: If you want porosity of the bone itself draw two contours: one around, counter-clockwise, on on the inside, clock-wise (thus a 'negative', excluding contour). Now, TV will be just that ring area, BV will be the bone within that ring area, and BV/TV will be 1 - porosity, thus BV/TV 94% means 6% porosity.

METAPHYSIS: Application of this script to metaphysis scan (with trab. in middle):

$ midshaft

Draw one contour around the cortex, and in the same slices, a NEGATIVE contour on the inside (CW contour) to exclude the trabecular bone. You may potentially work with the BP (break-point) system to be able to morph the negative contours, even though there is already an 'outer' contour around the bone in every slice.

Start 3D eval, read out the data-base with $ uct_list.

NOTE: You CANNOT give the .gobj a special name such as c0001234_cortex.gobj if you are using the automatic evaluation button on UCT_EVALUATION. Only after the eval (or midshaft eval) has finished, you may copy it to a more meaningful name, but all of the evaluations are just taking the STANDARD gobj name!

Afterwards (after renaming or copying the _seg.aim to _cortex.aim, and .gobj to _cortex.gobj), you may want to start the very normal eval with normal (inner) contours to evaluate the trabecular compartement.

Remember: the values in our 3D result database then are OVERWRITTEN.

If you have a large study, I would recommend to run the $ midshaft for ALL samples, $ uct_list to e.g. sys$scratch:mystudy_cortex.txt, rename/copy gobj and seg for all, and THEN run the trabecular evaluation for all samples. $ uct_list and you have two nice complete excel sheets.

NOTE ON PLATE MODEL: More explanations of why the plate model is correct

Assume a ring (cortex).

Area = thickness * 'middle line'.
'middle line' = average of outer and inner circumference
--> 'middle line' = 0.5 * 'bone surface'

--> thickness = bone volume * 2 / bone surface.

And this is nothing else but the 'plate model' formula for Tb.Th!

NORMAL TRABECULAR EVALUATION? And why not just run the 'normal trabecular eval'? Because in the 'normal' TRI calculation, we are EXCLUDING all surfaces that are created artificially by a gobj line that cuts through a trabecula.

However, for the midshaft, this means the whole outer (or even inner, if also inner contour is there) is SUPPRESSED, because it is touching the gobj line. So this is where that midshaft script comes in: it turns off that suppression.

MORE EXPLANATIONS: Why the direct DT-Tb.Th* is NOT correct for a scan with just 12 slices.

If we fill the structure with spheres, those spheres should be bounded by the width of the cortex, not just by the fact that we only have 12 slices (and even more so, we are always doing 'aim box' suppression with the spheres) if a sphere touches the boundary of the aim box, it is suppressed!

So this is quite crucial for Tb.Sp and Tb.N. To get meaningful Tb.Sp, you would need a LOT of slices. Also, then, you cannot have the inner contour - otherwise the negative inner contour will delete the Tb.Sp. spheres.

If you have enough slices, you may opt to turn the /dt_object on again in the midshaft script (second half). Then DT-Tb.Th* and TRI-Tb.Th just have a slightly different behaviour if there are little pores within the cortex:

TRI: the little pore just adds a little BS, thus the influence is not that big.
DT: the little pore cuts the thickness sphere in half, thus at that location, the thickness is halved.

MEAN1 and MEAN2: The script is also calculating Mean1 and Mean2 now.

Mean1: 'apparent' density, average over whole TV
Mean2: 'tissue' density: average over what was segmented as bone only.

These mean values will be in mg HA/ccm if you have the calibration option for bone density installed (and activated for those scans). Remember: the bone has to have some soft-tissue (or equivalent) on the inside, and preferably also on the outside then.

HowTo: Temporarily change the evaluation for any UCT_EVALUATION

In a terminal window, we need to redefine the logical that determines which 3D evaluation is associated with the UCT_EVALUATION tool.

For example, typing VOI runs the UCT_EVALUATION tool by executing UE:UCT_EVALUATION_VOI.COM

Type

$ show logical UCT_EVALUATION

and then you can redefine this to point to another evaluation script.

$ define UCT_EVALUATION UE:UCT_EVALUATION_CUTOUT.COM

This is only redefined in the particular terminal you are using, and it will reset to its default in a new terminal. It will also not work if you try to submit the job to a queue. It only works interactively.

HowTo: Assign special 3D evaluation to the UCT_EVALUATION tool

Warning: This is for advanced users ONLY. Do not attempt this without contacting Steve!

This will show how to set up $ eval_kyle: It is easy to set up multiple .COM files for the 3D evalulation tool, one such set-up is already prepared on your computer:

$ eval_special

will do something similar to $ voi where the normal contour program opens, but if you start the 3D eval, the special script um:uct_evaluation_special.com is called (associated with eval_special), which then in turn calls ipl_fixed_batch_special.com, where you can specify your own set of thresholds, etc.

So it is a daisy chain of files which we will need to follow to set up a special xxx_kyle.com script, and where ever we encounter _special.com, we copy and change it to _kyle.com.

So to set up $ eval_kyle you will need to do the following:

• 1. Edit disk1:[microct]microct_setup_user.com (double-click on file), please be careful here, it's very important! (and I hope you have done your 'weekly full backup' regularly, so you could always retrieve a good version of login.com...)

• 2. Then, towards to bottom and duplicate the line

$ eval_special :== @com:eval_var.com um:uct_evaluation_special.com and change duplicate to $ eval_kyle  :== @com:eval_var.com um:uct_evaluation_kyle.com

This will make this command eval_kyle available after you login freshly. Save and Exit.

• 3. Then go to disk1:[microct.eval] and copy the file uct_evaluation_special.com to the name uct_evaluation_kyle.com ('Commands', 'Copy ...', BUT CLICK ON THE THREE DOTS and give new output name.

If you do not see the eval_special: go to terminal, $ dir ue:*special*.com and take it from the microct_common directory)

• 4. Edit the _kyle.com file: at the top, change ipl_......_special.com to _kyle.com

• 5. Copy ipl_fixed_batch_special.com to ipl_fixed_batch_kyle.com

• 6. Edit ipl_fixed_batch_kyle.com to your liking. If you want $ eval_kyle to not perform the full evaluation, but only something akin to voi, then, in the steps above, copy uct_eval_voi.com to uct_evaluation_kyle.com, and ipl_voi.com to ipl_fixed_batch_kyle.com (thus you take the voi script as template for the _kyle script).

• 7. Then open a new terminal ('applications', 'decterm'), $ eval_kyle should work now!

Notes: When you start the 3D eval, then $ que will show you your batch job with your name on it now!

And now a few more explanations how the whole thing works: in all eval programs, hitting 'Start 3d' just submits the logical uct_evaluation, $ show log uct_evaluation shows you that by default, this points to ue:uct_evaluation_v4.com.

Those little voi and eval_steve scripts then just temporarily redefine this logical so that a different script is submitted when you 'start 3d eval'.

About the corelut eval: Then we start a whole different eval program: the .exe_corelut. So, if you set up a eval_steve_c (or whatever you want to call it): Produce a eval_var_c.com in the directory commands (maybe copy from microct_common to [microct.commands]), and there, on the last line, instead of starting $ uct_evaluation (the normal program), replace by $ evalc which starts the .exe_corelut - and because the eval_var.com has over-ridden the uct_evaluation logical just before, then you are starting different script upon starting 3d eval. NB: we have a logical AND a symbol called uct_evalution. Sorry for confusion. Symbol uct_evaluation starts normal eval program, so this is what you replace by symbol evalc now).

HowTo: Convert CT numbers to mg HA/ccm and or Houndsfield units (HU)

The two byte numbers (short) can be converted to mu values by first dividing by mu_scaling which is probably 4096 (uCT35, vivaCT40) or 8192 (XtremeCT).

Houndsfield units (HU)

Here is the conversion formula from mu (lin.att. in 1/cm, our original, 'physical' image numbers) to HU:

HU(pixel) = -1000 + mu(pixel) * 1000/mu(water)

so the mu(water) is the 'constant' that is needed for this conversion of any pixel value mu(pixel) to HU of that pixel.

divide the 2-byte numbers in our dicom/aim/isq files by 4096 to get mu. (4096 --> 'mu_scaling' in our headers). E.g. you see gray-scale number 8192 in our files --> mu = 2 [1/cm].

A reasonable result for the uCT35 is:

HU = (CT / 4096) * 1718.19 - 1000

and that means 1000/mu(water) = 1718.19 --> mu(water) = 0.58299 which is a reasonable number.

mg HA/ccm

mg HA/ccm = (AIM/mu_scaling) * slope + intercept
mg HA/ccm = (AIM/4096) * 295.494 - 189.695    ! vivaCT40
mg HA/ccm = (AIM/4096) * 256.51 - 172.99         ! uCT35
mg HA/ccm = (AIM/8192) * 1402.549 - 340.953  ! XtremeCT

These numbers are available in the AIM header, or the ISQ header ($ run ut:ctheader).

To look up slope, intercept and mu_scaling, check either the AIM header

$ aix C1234567.AIM

or

$ run ut:uct_update_calibration

and look at a record, and you'll find slope and intercept and mu-water.

vivaCT40
!---------------------------------------------------------------------- 
Mu_Scaling                    4096 
Calibration Data              55 kVp, 200mg HA wedge, 4096 scaling 
Calib. default unit type      2 (Density) 
Density: unit                 mg HA/ccm 
Density: slope                2.95493988e+02 
Density: intercept            -1.89695007e+02 
HU: mu water                  0.58200 
!----------------------------------------------------------------------

XtremeCT
!-------------------------------------------------------------------------------
Mu_Scaling                                       8192
Calibration Data              60kVp, BH: 200 mg HA/ccm, Scaling 8192            
Calib. default unit type      2 (Density)                                       
Density: unit                 mg HA/ccm                                         
Density: slope                         1.40254895e+03
Density: intercept                    -3.40953003e+02
HU: mu water                                  0.23900
!-------------------------------------------------------------------------------

HowTo: Create a list of registered measurements

• To generate a list of measurements that are registered in the data base, to determine where backups are stored, or if the data is still on disk:

$ write_info

• It asks you for a range of measurement numbers or sample names, and creates sys$scratch:info.txt

HowTo: View the scout views

Go to scout view directory

XtremeCT side:
$ cd DISK1:[MICROCT.DATA.SCOUTS]

VivaCT side:
$ cd ADISK1:[VIVACT.DATA.SCOUTS]

Use the program ximage

$ x asdf.SCV

HowTo: Rotate a dataset with IPL's turn3d function

The turn3d function can be used to perform a single rotation about any axis in 3d. The turnaxis_angle input is the angle between the x, y, and z axes and the line about which the rotation will be performed. The turnangle input is the amount to be rotated about the specified line. To use this function to perform three sequential rotations about the image's original x, y, and z axes requires compensation for the previous rotations.

Further instructions can be found here: turn3d.pdf

HowTo: Fix SAP colours when green

(using the standard gray colourmap), if changing color maps doesn't help.

• Menu point 'Norm, Default is Bone = 16000.

• Take max for 32000, or 'Adapt slice'.

Note: Do NOT take 'adapt vol' - it will be extremely slow! Every single mouse-click it re-calculates whole volume...

HowTo: Run SAP from the command line

$ sap -s myfile.aim

HowTo: Use SAP, CoreLUT or XV with XWin32 or Hummingbird Exceed

You need two things on your computer: (1) an X windows server (e.g., XWin32 or Hummingbird Exceed), and (2) an SSH client to run your X applications (e.g., Putty or SSH Secure Client).

For SAP, CoreLUT or XV to work with XWin32, the visual display must be set to 8-bit. The evaluation tool evalc can be used in place of uct_evaluation via the 8-bit display.

XWin32

In the XWin Configuration, choose the "Windows" tab, check the box for "Legacy PseudoColor Support". It may be necessary to restart Putty and Xwin32 before the settings take effect. If you want to switch back to 24-bit display, uncheck the "Legacy PseudoColor Support" and restart the programs.

Hummingbird Exceed

Under server settings, be sure to select "pseudocolor" display.

SSH Client

Make sure the X-tunneling is turned on so that X applications can be used with these SSH clients. Two SSH clients are recommended:

• putty
• SSH Secure Shell v3.2.9 http://www.ssh.com

HowTo: Fix Remove and Insert key function in Cutout or Evalc

You probably have caplocks activated. Deactivate it and the cut and paste functions of the remove and insert keys will be functional again.

Print outs and visualization

HowTo: Print out report sheets

The postscript report sheets can be generated and printed on the command line. For the XtremeCT system, use UPAT_PRINT_LIST.EXE which can be run (using the symbol):

$ UPAT_LIST

Look in the UU: directory to see options for the template files (UPAT*.TXT) and the postscript output files (*.PS).

• UU:UPAT_EVAL_FUP05_SCUSA.PS
• UU:UPAT_EVAL_BASELINE_SCUSA.PS

For other options, see the upat_list help:

$ UPAT_LIST -h

One particular option that may be useful, is to turn off slice matching:

$ UPAT_LIST -nomatch

For reports from non-patient measurements (i.e., the vivaCT system or non-patient scans on the XtremeCT), the template files are different (UCT*.TXT) and the postscript outputs are different (UCT*.PS). To print a result sheet use:

$ RUN UM:UCT_RAY_BATCH.EXE

If you want to create a postscript file, but not send it to the printer, set the number of copies to zero (interactive). Then, find the postscript file from the SYS$SCRATCH directory and transfer it to your PC. You must transfer the postscript file in ASCII mode. The files are typically called TEMP_3D.PS or UPAT_PAT_62_RR.PS. A program such as Adobe Photoshop can convert the postscript file to a different format.

See also our custom scripts:

$ @COM:PRINT_2D.COM
$ @COM:PRINT_3D.COM

Combining two volumes to create transparent image

First create your two start "seg.aim" images Start by running either the normal eval, or $ voi to get one _seg.aim file of one region

Rename it to a better name in the session manager, e.g. c0001234_all.aim

Then either by choosing different contours, or a different threshold, start the $ voi or the normal evaluation again, which will create another _seg.aim file.

Load them in the 3D viewer ('all files') to check that they look like expected.

Now we will go into IPL and concatenate those two volumes.

$ ipl
/read all xx_all.aim
/read trab xx_seg.aim
/transparent
  -input_solid               trab
  -input_transp              all
  -output                    out

Note: the solid volume is OVERWRITING the transparent input volume, so usually, it is the smaller (inner) of the two.

/write out xxx_transp.aim
/quit

Then load the _transp.aim file in 3D, go to 3D program FAQ to see how to change transparency levels in the Display options there.

See IPL FAQ or 3D program FAQ if you need help with these parts.

Step by Step Example: Showing your chosen trab VOI within whole bone

• $ voi, enter threshold in ipl_voi.com, start 3d eval WITHOUT contours
• WAIT UNTIL THE JOB IS FINISHED!
• Rename xxxx_seg.aim to xxxx_all.aim
• Start up normal eval program, draw trabecular contours, check threshold in ipl_fixed_batch_v4.com, start 3D eval.
• WAIT UNTIL THE JOB IS FINISHED!
• go into ipl and combine the two, see code above
• Voila, finished

Step by Step Example: Showing high-dense volume surronded by lower density volume

• Start normal 3D eval with HIGH threshold entered in ipl_fixed_batch_v4.com
• WAIT UNTIL THE JOB IS FINISHED!
• Rename xxxx_seg.aim to xxxx_high.aim
• If wished: read out 3d result database with $ uct_list
• Start normal 3D eval with LOW threshold entered in ipl_fixed_batch_v4.com
• WAIT UNTIL THE JOB IS FINISHED!
• go into ipl and combine the two, see code above
• Voila, finished

Make an animation

To calculate an Animation, two (separate) tools of the 3D program are used: the set-up of start and end point of the Animation, and the 'Auto-Print' tool that automatically creates the TIF files for every rendered image.

This is a 'step by step' guide to create an animation using TIFF files:

• Open your sample in the 3D program.
• Now, open 'options','Animation'
• Render your Image the way you want your animation to start. e.g. observer rotation -179 degrees. You may interrupt if you want. If you are satisfied with all the settings (rotation, elevation, transparency, colors etc.) the Click on 'Get Settings A'
• Then, adjust all the settings for the endpoint of your animation (with the mouse or the sliders) e.g. observer rotation +179, and start the rendering. If satisfied, click 'Get settings B'

Note: To test, put # of views to 5 or 10... You can also interrupt the build-ups for faster viewing.

• Click the START button again, wait till the build-up of the 3D image is finished.
• DO NOT HIT 'START CALCULATION' YET!
• First, you also need to set-up the 'Auto-Print' mode:
• Go to the 'Print' menu (file --> print), choose print to TIF and to File and enter THE BEGINNING OF A FILENAME, e.g. sys$scratch:test_ (the program will append 0000.TIF, 0001.TIF etc afterwards) Or enter (possibly with middle-mouse-click): "DK0:[MICROCT.DATA.00001234.00003456]ANIMATION_" to have the result go to the measurement directory.
• Click on the 'Auto-Print' box. Then 'ok'
• Adjust the number of views (e.g. 200), then click on 'Start Calculation', and it will render. This may well take a long time.
• Take those TIFF files to your PC, and combine the sequences in the program that you are using to convert the TIFF files to AVI or MOV files (like Windows Movie Maker 2, or Quick Time Pro, or Adobe Premiere etc.)

Notes/Hints:

• The order in the above list can also be reversed: You can first set up the 'Auto-Print', and then go into the 'Animation' menu, choose start and end settings, and hit 'Start Calculation'
• The positions A and B can be any: e.g. with subdim, you can first take everything (A), then slide it back to nothing (B).
• You can even slowly change the color of an object. Or slowly hiding objects: choose 'transparent' for object, but change from 0.0 (A) to 1.0 (B)
• If you want to create a movie with different sections, e.g. first rotate around, then slowly make bone transparent, then cut in: Do this in multiple sections with the 3D program, e.g. first create Auto-Print files 'rotation_0000 to rotation_0200', then 'transp_0000 to transp_0050', and then 'cut_0000 to cut_0100'. Make sure ALL the settings from the last view of the previous section match those of the first view of the next section so you get smoot transitions.
• The 'Auto-Print' box is also usable without the 'Animation' Menu: Activate it, Like above, with sys$scratch:test_ only the beginning of filename. 'ok' Then, whenever a 3D image is completely built-up (so when you are happy with a view and you do not interrupt the build-up), automatically a TIF file is generated, numbered sequentially.

Convert AIM to TIF series

Use the utility SEQ_TO_IMAGE located in the UT directory. A symbol has been set up to run this called TIF.

$ mcr ut:seq_to_image
Program: seq_to_image V1.1b
Input-File []: WISDOM_L_R05.AIM
File Format [TIF]: 
Default TIF Names: WISDOM_L_R05_xxxx.TIF

Output-File [WISDOM_L_R05]: 
!> Reading AimVersion020 
Start Slice [0]: 
End Slice [179]: 5
Number of images thereof (count) [6]: 
Scale Factor (Integer) [1]: 
global_pos (1=true, 0=false) [0]: 
Max value to norm to 255 [32767]: 
Slice: 0
Slice: 1
Slice: 2
Slice: 3
Slice: 4
Slice: 5
$ dir wisdom*

Directory DISK1:[MICROCT.SCRATCH]

WISDOM_L_R05.AIM;1  WISDOM_L_R05_0000.TIF;1                 WISDOM_L_R05_0001.TIF;1
WISDOM_L_R05_0002.TIF;1                 WISDOM_L_R05_0003.TIF;1
WISDOM_L_R05_0004.TIF;1                 WISDOM_L_R05_0005.TIF;1

Total of 7 files.
$ 

Writing Scripts

HowTo: Get information about a sample or measurement from database

The following will create a sys$scratch:info.txt file.

$ write_info

HowTo: Create an MSQ image and view it

This is based on marching cubes. See IPL_UPAT_EVAL.COM for an example as per below:

/read seg C0000214_SEG.AIM

/sup_divide                               ! careful of the order of inputs here
  -input                     seg
  -suppos_pixels_local       -1 -1 -1
  -supdim_numbers            1 2 1
  -subdim_pixels             -1 -1 -1
  -testoff_pixels            0 0 0

/sub_pick
  -input                     seg
  -output                    segsub
  -subpos_numbers            0 1 0

/tri
  -input                     segsub
  -output                    tri
  -gobj_filename             none
  -peel_iter                 0
  -ip_sigma                  0.8
  -ip_support                1
  -ip_threshold              40
  -nr_ave_iter               2
  -size_image                512 512
  -scale_image               0.800000
  -edges                     false
  -nr_views                  1             ! Important

/msq_from_aim
  -aim_name                  tri
  -msq_filename              C0000123_3D.MSQ

Then, use UCT_3D to view it. If you click on "animations" it will give you the option to read the MSQ file.

HowTo: Create your own IPL script with some tips

To make a script that would submit a job to the queue with the command as follows:

$ @um:ipl_eval_only.com 'aim_file' 'gobj_file' 'p1' 0

you need to create your own IPL file would be translated to the following command:

$ submit /noprint /que=sys$slow ue:IPL_VOI_JOSH_M.COM 
  - /log=sys$scratch:test.log /para=('aim_file','gobj_file','seg_file')

The IPL file would look something like this…

$   if p1 .EQS. "" 
$ 	THEN 
$	write sys$output "Give DK0:[xxx]C0001234_SEG.aim on the line."
$	write sys$output "Or: "
$	write sys$output "Type $ scanco_submit_command um:ipl_eval_autostart/para=(dk0:[xxx]c00123_seg.aim) "
$ 	exit
$   endif
$
$ gobj_file = p1-F$PARSE(p1,,,"VERSION")-"_SEG.AIM" + ".GOBJ"
$ aim_file = p1-F$PARSE(p1,,,"VERSION")-"_SEG.AIM" + ".AIM"
$
$!
$ @um:ipl_eval_only.com 'aim_file' 'gobj_file' 'p1' 0
$ exit

Some more handy IPL sequences are:

$    INPUT = F$SEARCH("UCT_MEASUREMENT_DATA:[*.*''P1']*.RSQ",1)

and then you just need to give meas. number as P1 (and maybe slice range as p2),

$ if p2 .eqs. "" then p2 := A

OR, some nice starts are:

$    IF P3 .EQS. ""
$    THEN
$      WRITE SYS$OUTPUT "P1: Measurement Number to duplicate"
$      WRITE SYS$OUTPUT "P2: New Sample Number from Database"
$      WRITE SYS$OUTPUT "P3: New Start Letter of ISQ file"
$      WRITE SYS$OUTPUT "P4: Optional: Operator Number"
$      WRITE SYS$OUTPUT "P5: Optional: Controlfile Number"
$      EXIT
$    ENDIF

where you say which parameters HAVE to be there (up to p3), and which ones are optional.

All together, the script may look like this:

$    IF P1 .EQS. "" 
$    THEN 
$      WRITE SYS$OUTPUT "P1: Measurement Number to reco" 
$      WRITE SYS$OUTPUT "P2: Optional: Slice number N or Range M:N. Auto:All" 
$      WRITE SYS$OUTPUT "P3: Optional: Queue. Auto: sys$slow" 
$      EXIT 
$    ENDIF 
$ if p2 .eqs. "" then p2 := A 
$ if p3 .eqs. "" then p3 := SYS$SLOW 
$ slices  = p2 
$ quename = p3 
$ INPUT = F$SEARCH("UCT_MEASUREMENT_DATA:[*.*''P1']*.RSQ",1) 
$ INPUT = INPUT - F$PARSE(INPUT,,,"VERSION") 
$ logfile = "sys$scratch:RECO_" + f$PARSE(input,,,"NAME") + ".LOG" 
$ submit/noprint/queue='quename' uct_reconstruction_file - 
    /para=('input','slices','logfile') 
$ exit

Note that uct_reconstruction_file is a logical that points to um:uct_reconstruction.com.

Logicals are preferred in scripts to make them more universal.

General

Wireless Network in Morphometrics Lab

The network is hidden so that we don't get any unwelcome guests or trouble from the internet Nazis in HiiTec.

SSID (network name): bhlab
WPA2 encryption key: hallgrimssonboyd

Errors and Fixes

HowTo: Align the whole body device to be in the field of view

It's important to have your whole sample in the field of view for accurate density measurements. Follow the instructions in this PDF to align the device.

HowTo: Fix "last slice" error problem when extracting AIM from ISQ

Sometimes the evaluation script (i.e., voi) converting the ISQ to an AIM fails because it tries to read one slice more than what is available in the ISQ data. The result of this failure is that the 3D analysis appears to fail and no .AIM or _SEG.AIM is created in the data directory. For example, it may have tried to read 261 slices, but there were only 260 slices available. The easiest work around is to make sure that the white bounding box in the uct_evaluation program does not include the last slice (largest Z dimension). Scanco is aware of this problem and will fix it in future software releases (V5.2).

HowTo: Fix sample data when results fail to print/list/show after evaluation

There are a number of reasons the results may not print:

Incorrect first measurement number

When sample information (such as patient name, DOB, remark) is changed DURING scanning, this can over write the 1st Scan measurement number, resulting in the evaluation sheet not being printed, and data not shown when running upat_list. For example, this happens when:

1) open sample program, enter name, hit save, leave program open

2) start scan (which would write correct record with 1st meas. in sample database)

3) in sample program: then only enter DOB, and hit 'save' again in the sample program --> it 'overwrites' the 1st meas. with an empty string again.

If you want to enter DOB after the scan has started, you HAVE to 'Find...' and re-load that patient first (so it correctly loads the 1st meas. entry also). It is best to just not use the sample management program whilst scanning is taking place.

However, if you do encounter this problem, use the following to fix (shown using patient 668 as an example. 1st scan measurement number was 2672):

 $ set host 0
 calibration
 xxxxx

 $ spb
 $ @[microct]microct_setup
 $ run [calibration.util]uct_correct_first
 $ r UCT_CORRECT_FIRST

 Patient Nr.: 668
 SampleNr.:          668
 First:            0

 2672

 Patient Nr.:

and continue with next non-working patient, give patient number 0 to exit.

$ logout

Evaluation completed (aim files produced) but results not written to database

If the first measurement number and DOB are correct... check the .log file of the evaluation... However, if there is no upat*.log file available for that scan, use:

upat_list -v=2 and you will see (for eg sample 594):

Found sample 594, CAMOS_0018
Found measurement 2536 made 19-MAR-2007 12:52:58.48
Measurement site RL found in list at index 0
Measurement 2536 within measurement period
ERROR finding evaluation for measurement 2536

No valid evaluation found for measurement 2536

So it seems the 3D upat eval never wrote anything to the database for that scan. In this case, you can run:

$ @ue:ipl_upat_eval_again_quick.com

and let it run for the missing data. NOTE it will only need the aim files and gobj files that you still have in your 'models' directory.

Incorrect sequence of measurement numbers

Finally, if none of these options work... check the order of measurement numbers: especially if an import has occured, first check the first measurement number is correct (as above). Then run sys$Scratch:uct_correct_prev_next.exe and correct the order of the measurements: The database keeps a 'daisy' chain of meas. numbers, each meas. knows which one was the previous meas and the next meas.

0 means: no prev. meas (or no next meas) So corrected the order for you like this:

$ r UCT_CORRECT_prev_next

Measurement No. (0 to stop) [0]: 2821
MeasNo.:         2821 PatientNo:          723

Previous [0]:   <-- leave 0
Next     [0]: 2828  <-- give next meas number
Measurement No. (0 to stop) [2828]:
MeasNo.:         2828 PatientNo:          723

Previous [0]: 2821  <-- give prev. meas number
Next     [0]:  <-- leave 0 for very last meas.

Please note: the series of meas. numbers per patient is NOT separate per meas. site, thus radius and tibia meas. are linking to each other. Do not mess up!

File Transfer

HowTo: SC5019 (vivaCT) files to Sodor

$ ftp hal.kin.ucalgary.ca
220 Microsoft FTP Service
Connected to hal.kin.ucalgary.ca. 
Name (hal.kin.ucalgary.ca:microct): hpl-nt\skboyd
331 Password required for hpl-nt\skboyd.
Password: 
230 User hpl-nt\skboyd logged in.
FTP> 

Note that you can log in as either:

• hpl-nt\skboyd
• skboyd@kin.ucalgary.ca

On the HPL network, open the window file program and type: \\hal\FTP Home

Other options for file transfer:

$ sftp   ! haven't tried it much
$ scp    ! doesn't work well with openvms

HowTo: FTP to douglas and donald (only in KIN network)

There are FTP servers running on douglas and donald for file transfer WITHIN the KIN network only, e.g. percy or your PC. If you are connecting from home, establish a PPTP VPN connection to nt-hpl first to be virtually in the KIN network. Access from scanners at the hospital (VivaCT) is not possible at the moment. Use user guest and correpsonding password for login to either douglas.kin.ucalgary.ca or donald.kin.ucalgary.ca:

$ douglas.kin.ucalgary.ca
220 Welcome to Dougla's FTP Service
Name (douglas.kin.ucalgary.ca:ypauchard): guest
331 Password required for guest.
Password: 
230 Login successful.
Remote system type is UNIX.
Using binary mode to transfer files.
FTP> 

Scanco Finite Element Analysis

Running the FE code on the VMS system

Read the PDF file that is located in the "documents" directory for documentation and pick a problem number and output type. Make sure that your image is suitable for the slected analysis (see page 7 in the Scanco FE guide). You can use your segmented AIM from morphological analysis for FE analysis.

Edit the IPl file used in batch job mode DISK1:[MICROCT.EVAL]IPL_FE_VM.COM on percy to include the desired analysis, e.g.:

FE_PROBLEM :== 9 !combination of 1 - 6
FE_PAR_NO :== 42 ! von Mises equivalent stress

Start the evaluation as a batch job in queue SC3309_FE

$ fem DISK3:[MICROCT.PROJECTS.DEFROM.MODELS]SHAM_A12_W12_SEG.AIM

Look at the log file SYS$SCRATCH:SHAM_A12_W12_SEG_FE.LOG during execution to verify the process.

The job is submitted using:

SUBMIT/NOPRINT/NOTIF/LOG=SYS$SCRATCH:SHAM_A12_W12_SEG_FE.LOG/QUE=UCT_FE_QUEUE UE:IPL_FE_VM.COM/PARA=DISK3:[MICROCT.PROJECTS.DEFROM.MODELS]SHAM_A12_W12_SEG.AIM

You can copy the IPL_FE_VM.COM to your project directory and use the above command to submit it.

Alternatively, FE analysis can be run interactively

$ iplfe

Extracting results

Creates a total result file called .post file. Then, by default, it extracts one parameter (e.g. von mises stress), and produces stats file .postlist and a colorful .aim file (_fe_par_42.aim) which you can view with uct_3d.

If you then want other parameters extracted, or other colorful images, you only need to do a very quick operation in iplfe:

$ iplfe
/fe_post

give .post file and write out the resulting .aim file, and it produces also the .postlist file.

Scanco Image Registration

On percy, check out the DISK1:[MICROCT.EVAL]IPL_REGISTER*.COM scripts. You can also use it interactively by starting

$ iplreg
ipl> register
ipl> transform

Use the ipl help function for more info and read the PDF guide (in Documents).

SimilarityRegistration

Understanding rotation parameters (versors)

In SimilarityRegistration, former AppSimilarityRegistrationV3, used when registering images with Shape, the config file allows for setting inital rotations initial_transform_axis_x, initial_transform_axis_y, initial_transform_axis_z. These correspond to the versor components and are calculated as follows:

initial_transform_axis_x = n_x * sin(alpha/2)
initial_transform_axis_y = n_y * sin(alpha/2)
initial_transform_axis_z = n_z * sin(alpha/2)

where n=[n_x,n_y,n_z]^T is the unit vecotr representing the axis of rotation and alpha is the angle of rotation. Example 45 degree rotation around x-axis:

initial_transform_axis_x = 1 * sin(45/2) = 0.382
initial_transform_axis_y = 0 * sin(45/2) = 0
initial_transform_axis_z = 0 * sin(45/2) = 0

The rotations after registration final_transform_axis_x, final_transform_axis_y, final_transform_axis_z can be interpreted the same way. ITK provides the functions GetAxis() and GetAngle to obtain [n_x, n_y, n_z] and alpha respectively, which are calculated as follows:

v_1 = initial_transform_axis_x
v_2 = initial_transform_axis_y
v_3 = initial_transform_axis_z
magnitude = sqrt(v_1^2 + v_2^2 + v_3^2)
n_x = v_1/magnitude
n_y = v_2/magnitude
n_z = v_3/magnitude
alpha = 2 * atan(magnitude/1)

Statistics and Analysis

JAGS (analysis of Bayesian models using Markov Chain Monte Carlo)

See JAGS Usage.