[feat] complete per-voxel mua/mus/g/n format support, close #203

fangq · Dec 29, 2023 · 59f18a7 · 59f18a7
1 parent 894ed6e
commit 59f18a7
Show file tree

Hide file tree

Showing 18 changed files with 119 additions and 31 deletions.
diff --git a/example/muamus/mua_mus_g_n.json b/example/muamus/mua_mus_g_n.json
@@ -0,0 +1,45 @@
+{
+	"Session":{
+		"ID":"mua_mus_g_n",
+		"DoAutoThread":1,
+		"Photons":1000000
+	},
+	"Forward":{
+		"T0":0,
+		"T1":5e-09,
+		"Dt":5e-09
+	},
+	"Optode":{
+		"Source":{
+			"Pos":[30,30,1],
+			"Dir":[0,0,1]
+		}
+	},
+	"Domain":{
+		"OriginType":0,
+		"Media":[
+			{
+				"mua":0,
+				"mus":0,
+				"g":1,
+				"n":1
+			},
+			{
+				"mua":0.005,
+				"mus":1,
+				"g":0,
+				"n":1.37
+			}
+		],
+		"MediaFormat":"asgn_float",
+		"Dim":[60,40,20],
+		"VolumeFile":""
+	},
+	"Shapes":{
+		"_ArrayType_":"single",
+		"_ArraySize_":[60,40,20,4],
+		"_ArrayZipType_":"zlib",
+		"_ArrayZipSize_":[1,192000],
+		"_ArrayZipData_":"eAHtx11r13UcAND13HCWS7OYONmKVmzlFhmr/b+/5WhKKyZMK2IDozSKTObFqgWb02GpSatJ5HLBzDBGWGBjIkZBNhbNapSbdTGkQlgLMSlcNYh6GZ+Lc+7Oof1f1OTlddXm/e+2lo9qD7m7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u3uo33/7jbmewzvd3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d096IcKGnPFzW+4u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7t70C+8uC334dVvu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7e9CPne3PpeF33d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3dPei7J47mTj75gbu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uwd9buS7XHPhkLu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uwf9n0PncjOffuLu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u5B/0jT5al904i7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7sH/Z1fFab8om/c3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d096K+tW5b2jU66u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7t70P92rCKVtZ1xd3d3d3d3d3d3d3d3d3d3d3d3d3d3d3f3oB+tujcNl067u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7t70B8cXJ3qx393d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3f3oN9a8nCa6Pjb3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3dg76574l05PpL7nN3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d/eYry5sTe2brnB3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d/egX7SrI608ke/u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u5B/2LXjpRfdI27u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7sH/VO5PWm89Tp3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3cP+nV/9aZ9o4vd3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d096OuG+tL64iXu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u4e9JWtA6msbZm7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7sHfXHF++n82E3u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u4e9AXTh9Nw6a3u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u4e9HMHh1Jne4W7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7sH/a/rj6f68Up3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3cP+tNLPk/zy1a4u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7t70I+c/jJNdNzj7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uQf/x3m9T/6nk7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uQX9gzWTaUF7n7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uQd8zbyqVb1/l7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uQd8x+kv644cGd3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3D/pnu2fS8eVr3N3d3d3d3d3d3d3d3d3d3d3d3d3d3d3dPegfq72Qtr+81t3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3YN+9dxsaph61N3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3YP+7qP/psK7Wtzd3d3d3d3d3d3d3d3d3d3d3d3d3d3d3T3oe49dmv24+3F3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3cP+l03XJkd+GmDu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7B31XW372dPUz7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uHvQvfF+QVfY85+7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7kG/uWpBNnt2i7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uwf9xp6F2We5593d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3T3oW84tzl7Z+5K7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7sH/doHi7LGmU53d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3cP+obBpdmrA9vc3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d096FdeVZKdH+t2d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3f3oK/eeHPWNLvD3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3dg375ibJsuHSnu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7B/0tJeVZUeNud3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3D/qlW+/IOtv3uLu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7e9AvmqrKfn7vNXd3d3d3d3d3d3d3d3d3d3d3d3d3d3d3dw/6eTUrsvrx193d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3T3oL+urzgbnet3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3T3o/7lYk80ve9Pd3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d2D/sK62mxL01vu7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u4e9NNH6rKJjj53d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3cP+jMLVmXVg/vd3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d096Cc3P5D1n3rH3d3d3d3d3d3d3d3d3d3d3d3d3d3d3d3dg/7rkw9leXkD7u7u7u7u7u7u7u7u7u7u7u7u7u7u7u7uHvT/AeowpPg="
+	}
+}
diff --git a/example/muamus/run_muamusgn.bat b/example/muamus/run_muamusgn.bat
@@ -0,0 +1 @@
+..\..\bin\mcx.exe -f mua_mus_g_n.json -D P %*
diff --git a/example/muamus/run_muamusgn.sh b/example/muamus/run_muamusgn.sh
@@ -0,0 +1,2 @@
+#!/bin/sh
+../../bin/mcx -f mua_mus_g_n.json -D P $@
diff --git a/mcxcloud/frontend/index.html b/mcxcloud/frontend/index.html
@@ -1566,6 +1566,8 @@ <h3>Backend</h3>
           "enum": [
             "byte",
             "short",
+            "integer",
+            "asgn_float",
             "svmc",
             "mixlabel",
             "labelplus",

diff --git a/mcxlab/examples/demo_continuous_mua_mus.m b/mcxlab/examples/demo_continuous_mua_mus.m
@@ -9,7 +9,7 @@
 
 % only clear cfg to avoid accidentally clearing other useful data
 clear cfg cfgs
-cfg.nphoton=1e8;
+cfg.nphoton=1e7;
 cfg.vol=uint8(ones(60,60,60));
 cfg.srcpos=[30 30 1];
 cfg.srcdir=[0 0 1];
@@ -60,3 +60,19 @@
 colormap(jet);
 title('fluence in continuously varying media (1/mm^2)')
 view([-25.5,21.2]);
+
+%% Approach #3: a 4 x Nx x Ny x Nz float32 array defines [mua,mus,g,n] per voxel
+
+cfg.vol(3,:,:,:)=0;
+cfg.vol(4,:,:,:)=1.37;
+
+flux=mcxlab(cfg);
+
+mcxplotvol(squeeze(double(cfg.vol(2,:,:,:))))
+title('continuously varying scattering coeff \mu_s (1/mm)')
+view([-25.5,21.2]);
+
+mcxplotvol(log10(double(flux.data)))
+colormap(jet);
+title('fluence in continuously varying media (1/mm^2)')
+view([-25.5,21.2]);
diff --git a/mcxlab/mcxlab.m b/mcxlab/mcxlab.m
@@ -51,6 +51,7 @@
 %                                         the optical properties are first read based on the label, then the i-th property is replaced by the property value
 %                        3 x Nx x Ny x Nz int16/uint16 array:  "mixlabel"/98 format: {[int16 label1][uint16 label2][0-65535 label1 percentage]}
 %                        4 x Nx x Ny x Nz uint8 array: "asgn_byte"/103 format: mua/mus/g/n gray-scale (0-255) interpolating between prop(2,:) and prop(3,:)
+%                        4 x Nx x Ny x Nz single array: "asgn_float"/96 format: per-voxel mua/mus/g/n floating point values (internally convert to half-precision)
 %                        2 x Nx x Ny x Nz uint16 array: "muamus_short"/104 format: mua/mus gray-scale (0-65535) interpolating between prop(2,:) and prop(3,:)
 %                        8 x Nx x Ny x Nz uint8 array: "svmc"/97 format: split-voxel MC (SVMC) hybrid domain, can be created by mcxsvmc.m
 %                        Example: <demo_continuous_mua_mus.m>. If voxel-based media are used, partial-path/momentum outputs are disabled

diff --git a/mcxstudio/mcxgui.pas b/mcxstudio/mcxgui.pas
@@ -2203,6 +2203,7 @@ procedure TfmMCX.sgConfigSelectEditor(Sender: TObject; aCol, aRow: Integer;
          TCustomComboBox(Editor).Items.Add('mua_float - per voxel mua in float');
          TCustomComboBox(Editor).Items.Add('muamus_half - per voxel mua/mus in half-float');
          TCustomComboBox(Editor).Items.Add('asgn_byte - per voxel mua/mus/g/n grayscale in byte');
+         TCustomComboBox(Editor).Items.Add('asgn_float - per voxel mua/mus/g/n float values');
          TCustomComboBox(Editor).Items.Add('muamus_short - per voxel mua/mus grayscale in short');
      end;
    end;

diff --git a/pmcx/README.md b/pmcx/README.md
@@ -4,7 +4,7 @@
 
 - Copyright: (C) Matin Raayai Ardakani (2022-2023) <raayaiardakani.m at northeastern.edu>, Qianqian Fang (2019-2023) <q.fang at neu.edu>, Fan-Yu Yen (2023) <yen.f at northeastern.edu>
 - License: GNU Public License V3 or later
-- Version: 0.2.7
+- Version: 0.2.8
 - URL: https://pypi.org/project/pmcx/
 - Github: https://github.com/fangq/mcx
 

diff --git a/pmcx/pmcx/__init__.py b/pmcx/pmcx/__init__.py
@@ -49,7 +49,7 @@
 # from .files import loadmc2, loadmch, load, save
 from .bench import bench
 
-__version__ = "0.2.7"
+__version__ = "0.2.8"
 
 __all__ = (
     "gpuinfo",

diff --git a/pmcx/setup.py b/pmcx/setup.py
@@ -123,7 +123,7 @@ def build_extension(self, ext):
 setup(
     name="pmcx",
     packages=["pmcx"],
-    version="0.2.7",
+    version="0.2.8",
     requires=["numpy"],
     license="GPLv3+",
     author="Matin Raayai Ardakani, Qianqian Fang, Fan-Yu Yen",
@@ -133,7 +133,7 @@ def build_extension(self, ext):
     long_description_content_type="text/markdown",
     maintainer="Qianqian Fang",
     url="https://github.com/fangq/mcx",
-    download_url="http://mcx.space",
+    download_url="https://mcx.space",
     keywords=[
         "Monte Carlo simulation",
         "Biophotonics",

diff --git a/schema/mcxinput.json b/schema/mcxinput.json
@@ -989,6 +989,8 @@
           "enum": [
             "byte",
             "short",
+            "integer",
+            "asgn_float",
             "svmc",
             "mixlabel",
             "labelplus",

diff --git a/src/mcx_core.cu b/src/mcx_core.cu
@@ -594,6 +594,11 @@ __device__ void updateproperty(Medium* prop, unsigned int& mediaid, RandType t[R
         prop->mus = fabsf(__half2float(val.h[1]));
         prop->n = gproperty[!(mediaid & MED_MASK) == 0].w;
     } else if (gcfg->mediaformat == MEDIA_ASGN_F2H) { //< [h3][h2][h1][h0]: h3/h2/h1/h0: half-prec n/g/mus/mua for every voxel
+        if (idx1d == OUTSIDE_VOLUME_MIN || idx1d == OUTSIDE_VOLUME_MAX) {
+            *((float4*)(prop)) = gproperty[0]; // out-of-bounds
+            return;
+        }
+
         union {
             unsigned int i[2];
 #if ! defined(__CUDACC_VER_MAJOR__) || __CUDACC_VER_MAJOR__ >= 9
@@ -602,6 +607,7 @@ __device__ void updateproperty(Medium* prop, unsigned int& mediaid, RandType t[R
             half h[4];
 #endif
         } val;
+
         val.i[0] = mediaid & MED_MASK;
         val.i[1] = media[idx1d + gcfg->dimlen.z];
         prop->mua = fabsf(__half2float(val.h[0]));

diff --git a/src/mcx_utils.c b/src/mcx_utils.c
@@ -1546,10 +1546,13 @@ void mcx_preprocess(Config* cfg) {
 
     if (cfg->isrowmajor) {
         /*from here on, the array is always col-major*/
-        if (cfg->mediabyte == MEDIA_2LABEL_SPLIT || cfg->mediabyte == MEDIA_ASGN_F2H) {
-            mcx_convertrow2col64((size_t**) & (cfg->vol), &(cfg->dim));
+        if (cfg->mediabyte == MEDIA_2LABEL_SPLIT) {
+            mcx_convertrow2col64((size_t*) (cfg->vol), &(cfg->dim));
+        } else if (cfg->mediabyte == MEDIA_ASGN_F2H) {
+            mcx_convertrow2col(cfg->vol, &(cfg->dim));
+            mcx_convertrow2col(cfg->vol + cfg->dim.x * cfg->dim.y * cfg->dim.z, &(cfg->dim));
         } else {
-            mcx_convertrow2col(&(cfg->vol), &(cfg->dim));
+            mcx_convertrow2col(cfg->vol, &(cfg->dim));
         }
 
         cfg->isrowmajor = 0;
@@ -3747,12 +3750,12 @@ void mcx_loadseedfile(Config* cfg) {
  * @param[in] dim: the dimensions of the 3D array
  */
 
-void  mcx_convertrow2col(unsigned int** vol, uint3* dim) {
+void  mcx_convertrow2col(unsigned int* vol, uint3* dim) {
     uint x, y, z;
     unsigned int dimxy, dimyz;
     unsigned int* newvol = NULL;
 
-    if (*vol == NULL || dim->x == 0 || dim->y == 0 || dim->z == 0) {
+    if (vol == NULL || dim->x == 0 || dim->y == 0 || dim->z == 0) {
         return;
     }
 
@@ -3763,11 +3766,11 @@ void  mcx_convertrow2col(unsigned int** vol, uint3* dim) {
     for (x = 0; x < dim->x; x++)
         for (y = 0; y < dim->y; y++)
             for (z = 0; z < dim->z; z++) {
-                newvol[z * dimxy + y * dim->x + x] = (*vol)[x * dimyz + y * dim->z + z];
+                newvol[z * dimxy + y * dim->x + x] = vol[x * dimyz + y * dim->z + z];
             }
 
-    free(*vol);
-    *vol = newvol;
+    memcpy(vol, newvol, sizeof(unsigned int) * dim->x * dim->y * dim->z);
+    free(newvol);
 }
 
 /**
@@ -3777,12 +3780,12 @@ void  mcx_convertrow2col(unsigned int** vol, uint3* dim) {
  * @param[in] dim: the dimensions of the 3D array
  */
 
-void  mcx_convertrow2col64(size_t** vol, uint3* dim) {
+void  mcx_convertrow2col64(size_t* vol, uint3* dim) {
     uint x, y, z;
     size_t dimxy, dimyz;
     size_t* newvol = NULL;
 
-    if (*vol == NULL || dim->x == 0 || dim->y == 0 || dim->z == 0) {
+    if (vol == NULL || dim->x == 0 || dim->y == 0 || dim->z == 0) {
         return;
     }
 
@@ -3793,11 +3796,11 @@ void  mcx_convertrow2col64(size_t** vol, uint3* dim) {
     for (x = 0; x < dim->x; x++)
         for (y = 0; y < dim->y; y++)
             for (z = 0; z < dim->z; z++) {
-                newvol[z * dimxy + y * dim->x + x] = (*vol)[x * dimyz + y * dim->z + z];
+                newvol[z * dimxy + y * dim->x + x] = (vol)[x * dimyz + y * dim->z + z];
             }
 
-    free(*vol);
-    *vol = newvol;
+    memcpy(vol, newvol, sizeof(sizeof(size_t) * dim->x * dim->y * dim->z));
+    free(newvol);
 }
 
 /**
@@ -3826,8 +3829,8 @@ void  mcx_convertcol2row(unsigned int** vol, uint3* dim) {
                 newvol[x * dimyz + y * dim->z + z] = (*vol)[z * dimxy + y * dim->x + x];
             }
 
-    free(*vol);
-    *vol = newvol;
+    memcpy(*vol, newvol, sizeof(unsigned int) * dim->x * dim->y * dim->z);
+    free(newvol);
 }
 
 /**

diff --git a/src/mcx_utils.h b/src/mcx_utils.h
@@ -290,8 +290,8 @@ int  mcx_remap(char* opt);
 void mcx_maskdet(Config* cfg);
 void mcx_dumpmask(Config* cfg);
 void mcx_version(Config* cfg);
-void mcx_convertrow2col(unsigned int** vol, uint3* dim);
-void mcx_convertrow2col64(size_t** vol, uint3* dim);
+void mcx_convertrow2col(unsigned int* vol, uint3* dim);
+void mcx_convertrow2col64(size_t* vol, uint3* dim);
 void mcx_convertcol2row(unsigned int** vol, uint3* dim);
 void mcx_convertcol2row4d(unsigned int** vol, uint4* dim);
 int  mcx_loadjson(cJSON* root, Config* cfg);

diff --git a/src/mcxlab.cpp b/src/mcxlab.cpp
@@ -683,6 +683,10 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
                 float f2h[2];
                 int offset = (cfg->mediabyte == MEDIA_ASGN_F2H);
 
+                if (cfg->mediabyte == MEDIA_ASGN_F2H) {
+                    cfg->vol = static_cast<unsigned int*>(realloc(cfg->vol, dimxyz * 2 * sizeof(unsigned int)));
+                }
+
                 for (i = 0; i < dimxyz; i++) {
                     f2h[0] = val[i << (1 + offset)] * cfg->unitinmm;       // mua
                     f2h[1] = val[(i << (1 + offset)) + 1] * cfg->unitinmm; // mus

diff --git a/src/pmcx.cpp b/src/pmcx.cpp
@@ -303,25 +303,28 @@ void parseVolume(const py::dict& user_cfg, Config& mcx_config) {
                     auto val = (float*) buffer.ptr;
 
                     float f2h[2];
-                    int offset = (cfg->mediabyte == MEDIA_ASGN_F2H);
+                    int offset = (mcx_config.mediabyte == MEDIA_ASGN_F2H);
+                    if(mcx_config.mediabyte == MEDIA_ASGN_F2H) {
+                        mcx_config.vol = static_cast<unsigned int*>(realloc(mcx_config.vol, dim_xyz * 2 * sizeof(unsigned int)));
+                    }
 
                     for (i = 0; i < dim_xyz; i++) {
                         f2h[0] = val[i << (1 + offset)] * mcx_config.unitinmm;        // mua
-                        f2h[1] = val[(i << (1 + offset)) + 1] * mcx_config->unitinmm; // mus
+                        f2h[1] = val[(i << (1 + offset)) + 1] * mcx_config.unitinmm; // mus
 
                         if (f2h[0] != f2h[0]
                                 || f2h[1] != f2h[1]) { /*if one of mua/mus is nan in continuous medium, convert to 0-voxel*/
                             mcx_config.vol[i] = 0;
                             continue;
                         }
 
-                        if (cfg->mediabyte == MEDIA_ASGN_F2H) {
-                            mcx_config->vol[i] = mcx_float2half2(f2h);
+                        if (mcx_config.mediabyte == MEDIA_ASGN_F2H) {
+                            mcx_config.vol[i] = mcx_float2half2(f2h);
                             f2h[0] = val[(i << 2) + 2];   // g
                             f2h[1] = val[(i << 2) + 3];   // n
-                            mcx_config->vol[i + dim_xyz] = mcx_float2half2(f2h);
+                            mcx_config.vol[i + dim_xyz] = mcx_float2half2(f2h);
                         } else {
-                            mcx_config->vol[i] = mcx_float2half2(f2h);
+                            mcx_config.vol[i] = mcx_float2half2(f2h);
                         }
                     }
 

diff --git a/utils/json2mcx.m b/utils/json2mcx.m
@@ -87,11 +87,11 @@
             mediaclass='uint8';
             if(isfield(json.Domain,'MediaFormat'))
                 idx=find(ismember({'byte','short','integer','muamus_float',...
-                    'mua_float','muamus_half','asgn_byte','muamus_short','svmc'},...
+                    'mua_float','muamus_half','asgn_byte','muamus_short','svmc','asgn_float'},...
                     lower(json.Domain.MediaFormat)));
                 if(idx)
-                    typebyte=[1,2,4,8,4,4,4,4,8];
-                    typenames={'uint8','uint16','uint32','single','single','uint16','uint8','uint16','uint8'};
+                    typebyte=[1,2,4,8,4,4,4,4,8,16];
+                    typenames={'uint8','uint16','uint32','single','single','uint16','uint8','uint16','uint8','single'};
                     bytelen=typebyte(idx);
                     mediaclass=typenames{idx};
                 else

diff --git a/utils/mcx2json.m b/utils/mcx2json.m
@@ -107,6 +107,8 @@ function mcx2json(cfg,filestub)
                     Domain.MediaFormat='mua_float';
                 elseif(size(cfg.vol,1)==2)
                     Domain.MediaFormat='muamus_float';
+                elseif(size(cfg.vol,1)==4)
+                    Domain.MediaFormat='asgn_float';
                 end
             end
         otherwise