Default origin in MCXLAB and Iso2Mesh

[Hirvi Pauliina]

Dear Dr. Fang and MCX-/Iso2Mesh-users,

In MATLAB's intrinsic spatial coordinate system, the extreme outer corner of voxel (1,1,1) (the corner of the diagonal direction of the first voxel) is at (0.5,0.5,0.5), and the center of the voxel is at (1,1,1); see https://se.mathworks.com/help/images/image-coordinate-systems.html.

I would like to present the Iso2Mesh mesh nodes, and the MCXLAB src/det positions and the photon exit positions with respect to this MATLAB origin (for visualization). I wonder, if the following empirical interpretations are correct (sorry, if I missed some relevant documentation):

1) Does Iso2Mesh by default assume that the extreme corner of voxel (1,1,1) is at (0,0,0), so you have to translate the node coordinates by +[0.5,0.5,0.5] to present them with respect to the MATLAB origin? See the issue I posted previously for details: https://github.com/fangq/iso2mesh/issues/76

2) Do the MCXLAB photon exit positions also assume that the extreme corner of voxel (1,1,1) is at (0,0,0), so you have to translate the positions by +[0.5,0.5,0.5] to present them with respect to the MATLAB origin? This is despite setting "cfg.issrcfrom0=0" ,or independent of cfg.issrcfrom0?

3) When giving the src/det positions to MCXLAB with "cfg.issrcfrom0=0", we assume that the extreme corner of voxel (1,1,1) is at (1,1,1); see FAQ 8 in https://mcx.space/wiki/index.cgi?Doc/FAQ. Thus, we have to translate these positions by [-0.5, -0.5, -0.5] to present them with respect to the MATALB origin? Furthermore, if we select the positions wrt the MATLAB origin, we have to translate the coordinates by +[0.5,0.5,0.5] before giving them to MCXLAB with cfg.issrcfrom0=0, or by  [-0.5, -0.5, -0.5] with cfg.issrcfrom0=1?

4) Are the exit directions of photons the directions that they come towards the surface, or is refraction at the tissue-air interface accounted for with the Snell's law?

Attached is a short script and the corresponding figure for testing points 2) and 3). Since knowing the origin is essential when setting the detector location and radius, and when checking them afterwards from the spread of exit positions, I would greatly appreciate any help! Many thanks in advance.

Best regards,

Pauliina Hirvi

[Qianqian Fang]

hi Paullina,

sorry for my late reply to your iso2mesh ticket.

the iso2mesh mesh coordinates were adjusted so that the lower-bottom corner of the first voxel of a 3D volume is treated as [0,0,0].

the offset you observed was a result of using imagesc - which displays the center of the first voxel as [1,1] instead of [0.5,0.5] as in the iso2mesh mesh coordinates.

you can try these commands

>> [no,fc]=binsurface(box);
>> min(no)

ans =

     1     1     1

>> max(no)

ans =



     5     5     5







>> [no,el]=v2s(box,0.5,0.2);
extracting surfaces from a volume ...
region 1 centroid :    5.333333 5.666667 4.000000

processing threshold level 1...
Warning: You are meshing the surface with sub-pixel size. If this is not your your intent, please check if you set
...

done.
Final number of points: 1622
surface mesh generation is complete
>> min(no)

ans =

    0.9999    0.9999    0.9999

>> max(no)

ans =

    5.0001    5.0001    5.0001







>> [no,el,~,~]=vol2surf(box,1:1:size(box,1),1:1:size(box,2),1:1:size(box,3),1,1);
extracting surfaces from a volume ...
region 1 centroid :    5.333333 5.666667 4.000000
...
done.
Final number of points: 189
surface mesh generation is complete
>> min(no)

ans =

    1.4998    1.4998    1.4998

>> max(no)

ans =

    4.5002    4.5002    4.5002

you can see that in all 3 cases, the center of the output mesh is around [3,3,3], which is the center of a 6x6x6 cube of your input volume "box".

the reason that your vol2surf() created mesh shrinks by 0.5 voxel in each side is because you use a threshold of 1 instead of 0.5 - using 0.5, the interpolation between binary 0-1 value should split right at the middle.

regardless - this is rather an iso2mesh convention, and most iso2mesh functions have been somewhat adjusted to follow this assumption. This assumption not only is used when converting a volume to a mesh, but also the reverse - rasterizing a surface to a volume, as shown by the match of contour lines in this example

https://github.com/fangq/iso2mesh/blob/master/sample/demo_surf2vol_ex1.m

this is different from mcx's issrcfrom0 flag - mcx uses a voxelated grid, has only two possible options - either issrcfrom0=0 by treating the bottom-lower corner of the 1st voxel as [1,1,1], or, when issrcfrom0=1 by setting that corner as [0,0,0]. This setting is unrelated to iso2mesh mesh/volume mapping. The first option was rather confusing, but we had it in mcx from the beginning just to be compatible with the tMCimg input file that David Boas developed.

over time, in many of my demo scripts, I have set issrcfrom0 to be 1. it is my intent to set it default to 1 instead of 0 in the future release to eliminate this confusion.

let me know if my comments above make sense to you, and if I did not address any part of your questions.

Qianqian

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/AS8P192MB132096995D57EE74E791FAA3E92B2%40AS8P192MB1320.EURP192.PROD.OUTLOOK.COM.

[Hirvi Pauliina]

Dear Dr. Fang,

Many thanks for your quick response and the very useful information on Iso2Mesh!

My questions 2–4 were actually more related to MCXLAB and matching the outputs with the MATLAB (image) origin. I try to repeat the questions here in a shorter form, if you could still provide some comments on them:

2) Do the photon exit positions assume that the extreme corner of voxel (1,1,1) is at (0,0,0), so you have to translate the positions by +[0.5,0.5,0.5] to present them with respect to the MATLAB origin (independent of cfg.issrcfrom0)? The originally attached figure shows "misplacement" in MATLAB image frame without this operation.

3) When giving the src/det positions to MCXLAB with "cfg.issrcfrom0=0", we assume that the extreme corner of voxel (1,1,1) is at (1,1,1), thus we have to translate these positions by [-0.5, -0.5, -0.5] to present them with respect to the MATALB origin?

4) Are the exit directions of photons the directions that they come towards the surface, or is refraction at the tissue-air interface accounted for with the Snell's law?

Many thanks once again for your kind help,
Pauliina

---

From: Qianqian Fang <fan...@northeastern.edu>
Sent: Tuesday, March 12, 2024 10:28 PM
To: mcx-...@googlegroups.com <mcx-...@googlegroups.com>; Hirvi Pauliina <pauliin...@aalto.fi>
Subject: Re: [mcx-users] Default origin in MCXLAB and Iso2Mesh

 

[Qianqian Fang]

hi Pauliina,

sorry that I did not read it carefully - see my comments below

On 3/12/24 17:52, Hirvi Pauliina wrote:

Dear Dr. Fang,

Many thanks for your quick response and the very useful information on Iso2Mesh!

My questions 2–4 were actually more related to MCXLAB and matching the outputs with the MATLAB (image) origin. I try to repeat the questions here in a shorter form, if you could still provide some comments on them:

2) Do the photon exit positions assume that the extreme corner of voxel (1,1,1) is at (0,0,0), so you have to translate the positions by +[0.5,0.5,0.5] to present them with respect to the MATLAB origin (independent of cfg.issrcfrom0)? The originally attached figure shows "misplacement" in MATLAB image frame without this operation.

the photon exit positions are computed inside the GPU kernel - in the GPU, all positions assume the extreme corner of first voxel is [0,0,0].

also, there is no [0.5,0.5,0.5] - mcx either subtract 1 from srcpos/detpos, or subtract nothing. there is no half-voxel shift.

you can run `grep issrcfrom0 *` inside mcx/src, and you can see that the only place issrcfrom0 is used is at the handling of user input:

if issrcfrom0 is 0, srcpos.x/y/z and detpos.x/y/z subtract 1

https://github.com/fangq/mcx/blob/875a35786fde240ca296ba962e3dab1033b5e001/src/mcx_utils.c#L1925-L1937

https://github.com/fangq/mcx/blob/875a35786fde240ca296ba962e3dab1033b5e001/src/mcx_utils.c#L2589-L2591

https://github.com/fangq/mcx/blob/875a35786fde240ca296ba962e3dab1033b5e001/src/mcx_utils.c#L2682-L2686

https://github.com/fangq/mcx/blob/875a35786fde240ca296ba962e3dab1033b5e001/src/mcx_utils.c#L2947-L2951

from that point on, everything in the code assumes the origin is [0,0,0], including everything in the GPU kernels (and outputs)

3) When giving the src/det positions to MCXLAB with "cfg.issrcfrom0=0", we assume that the extreme corner of voxel (1,1,1) is at (1,1,1), thus we have to translate these positions by [-0.5, -0.5, -0.5] to present them with respect to the MATALB origin?

if you set cfg.srcpos/cfg.detpos in matlab, and set cfg.issrcfrom0=0, all it does is to trigger the -1 offset as shown in the above links (by calling mcx_validatecfg() in mcx_utils.c).

there is no 0.5 offset.

4) Are the exit directions of photons the directions that they come towards the surface, or is refraction at the tissue-air interface accounted for with the Snell's law?

if you use MCX releases after 2017, the detected existing directions considers refraction

however, before 2017, mcx did not perform that extra step, and this was fixed in the below commit and carried on since then

https://github.com/fangq/mcx/commit/c7895a6c5a70b82973f6aa9739cf0c6f780a0a53

let me know if this is helpful.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/AS8P192MB13205E8D19C1668C3E4F4916E92B2%40AS8P192MB1320.EURP192.PROD.OUTLOOK.COM.

[Hirvi Pauliina]

Dear Dr. Fang,

Many thanks for the careful explanations! Now I know how to perform the half-voxel shifts correctly in MATLAB when visualizing MCX inputs and outputs in MATLAB image frame. 

 

I'm also very grateful for the advice on how to get the Iso2Mesh surface meshes to follow the voxel boundaries exactly (no shrink). I guess you only need to specify the level sets/"isovalues" if you use "cgalsurf"? When using "cgalmesh" with v2s.m, v2m.m or vol2mesh.m, the "isovalues" input is not used, since eventually the function will call "[node,elem,face]=cgalv2m(vol,opt,maxvol)". Thus, it seems that "cgalmesh" will always output meshes that follow exterior and interior regions' voxel boundaries exactly?

One final question, if you may: How should one interpret the node indices, i.e., the fourth column values from "cgalmesh"? Why are they 1.5 in volume type 1, instead of 1 or some range?

Many thanks for all the help!

Best regards,

Pauliina

---

From: Qianqian Fang <fan...@northeastern.edu>
Sent: Wednesday, March 13, 2024 12:20 AM

[Qianqian Fang]

On 3/14/24 10:08, Hirvi Pauliina wrote:

Dear Dr. Fang,

Many thanks for the careful explanations! Now I know how to perform the half-voxel shifts correctly in MATLAB when visualizing MCX inputs and outputs in MATLAB image frame. 

 

I'm also very grateful for the advice on how to get the Iso2Mesh surface meshes to follow the voxel boundaries exactly (no shrink). I guess you only need to specify the level sets/"isovalues" if you use "cgalsurf"?

correct, cgalsurf determines the boundary using linear interpolation based upon the specified isovalue.

When using "cgalmesh" with v2s.m, v2m.m or vol2mesh.m, the "isovalues" input is not used, since eventually the function will call "[node,elem,face]=cgalv2m(vol,opt,maxvol)". Thus, it seems that "cgalmesh" will always output meshes that follow exterior and interior regions' voxel boundaries exactly?

for cgalmesh option, this can be achieved by adjusting the `opt.distbound` option in v2m - the exact definition of this parameter is denoted as "facet_distance" parameter

https://doc.cgal.org/latest/Mesh_3/index.html#:~:text=surface%20Delaunay%20balls.-,facet_distance,-.%20This%20parameter

my understanding is that this number limits how far away the triangles can deviate from the label based boundary. a small number means the boundary is more accurate, but requires lots of refinement, thus denser near the boundary.

One final question, if you may: How should one interpret the node indices, i.e., the fourth column values from "cgalmesh"? Why are they 1.5 in volume type 1, instead of 1 or some range?

to be honest, I never used the nodal labels returned by cgalmesh, I also wasn't able to find a documentation on its meaning, but it appears that the decimal numbers suggests the node is at the boundary of two labels.

I usually just use the element labels

Qianqian