initial support of negative source and negative-patterns

some test still fails, but feature is mostly working, need more work

src/mcx_utils.c

@@ -314,6 +314,7 @@ void mcx_initcfg(Config* cfg) {
     cfg->detectedcount = 0;
     cfg->runtime = 0;
     cfg->faststep = 0;
+    cfg->srcpos.w = 1.f;
     cfg->srcdir.w = 0.f;
     cfg->issaveref = 0;
     cfg->isspecular = 0;
@@ -2480,6 +2481,8 @@ int mcx_loadjson(cJSON* root, Config* cfg) {
                 cfg->srcpos.z = subitem->child->next->next->valuedouble;
             }
 
+            cfg->srcpos.w = FIND_JSON_KEY("Weight", "Optode.Source.Weight", src, 1.f, valuedouble);
+
             subitem = FIND_JSON_OBJ("Dir", "Optode.Source.Dir", src);
 
             if (subitem) {

src/mcxlab.cpp

@@ -392,15 +392,29 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
                 fieldlen = fielddim[0] * fielddim[1] * fielddim[2] * fielddim[3] * fielddim[4] * fielddim[5];
 
                 if (cfg.issaveref) {
+                    int highdim = fielddim[3] * fielddim[4] * fielddim[5];
+                    int voxellen = cfg.dim.x * cfg.dim.y * cfg.dim.z;
                     float* dref = (float*)malloc(fieldlen * sizeof(float));
                     memcpy(dref, cfg.exportfield, fieldlen * sizeof(float));
 
-                    for (int i = 0; i < fieldlen; i++) {
-                        if (dref[i] < 0.f) {
-                            dref[i] = -dref[i];
-                            cfg.exportfield[i] = 0.f;
-                        } else {
-                            dref[i] = 0.f;
+                    for (int ix = 0; ix < cfg.dim.x; ix++) {
+                        for (int iy = 0; iy < cfg.dim.y; iy++) {
+                            for (int iz = 0; iz < cfg.dim.z; iz++) {
+                                int voxelid = iz * (cfg.dim.x * cfg.dim.y) + iy * cfg.dim.x + ix;
+
+                                if (cfg.vol[voxelid]) {
+                                    for (int gate = 0; gate < highdim; gate++)
+                                        for (int srcid = 0; srcid < cfg.srcnum; srcid++) {
+                                            dref[(gate * voxellen + voxelid) * cfg.srcnum + srcid] = 0.f;
+                                        }
+                                } else {
+                                    for (int gate = 0; gate < highdim; gate++)
+                                        for (int srcid = 0; srcid < cfg.srcnum; srcid++) {
+                                            dref[(gate * voxellen + voxelid) * cfg.srcnum + srcid] = -dref[(gate * voxellen + voxelid) * cfg.srcnum + srcid];
+                                            cfg.exportfield[(gate * voxellen + voxelid) * cfg.srcnum + srcid] = 0.f;
+                                        }
+                                }
+                            }
                         }
                     }
 
@@ -551,7 +565,7 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
     GET_ONE_FIELD(cfg, omega)
     GET_ONE_FIELD(cfg, issave2pt)
     GET_ONE_FIELD(cfg, lambda)
-    GET_VEC3_FIELD(cfg, srcpos)
+    GET_VEC34_FIELD(cfg, srcpos)
     GET_VEC34_FIELD(cfg, srcdir)
     GET_VEC3_FIELD(cfg, steps)
     GET_VEC3_FIELD(cfg, crop0)

src/pmcx.cpp

@@ -458,7 +458,7 @@ void parse_config(const py::dict& user_cfg, Config& mcx_config) {
     GET_SCALAR_FIELD(user_cfg, mcx_config, srcnum, py::int_);
     GET_SCALAR_FIELD(user_cfg, mcx_config, omega, py::float_);
     GET_SCALAR_FIELD(user_cfg, mcx_config, lambda, py::float_);
-    GET_VEC3_FIELD(user_cfg, mcx_config, srcpos, float);
+    GET_VEC34_FIELD(user_cfg, mcx_config, srcpos, float);
     GET_VEC34_FIELD(user_cfg, mcx_config, srcdir, float);
     GET_VEC3_FIELD(user_cfg, mcx_config, steps, float);
     GET_VEC3_FIELD(user_cfg, mcx_config, crop0, uint);
@@ -1037,15 +1037,29 @@ py::dict pmcx_interface(const py::dict& user_cfg) {
             auto dref_array = py::array_t<float, py::array::f_style>(array_dims);
 
             if (mcx_config.issaveref) {
+                int highdim = fielddim[3] * fielddim[4] * fielddim[5];
+                int voxellen = cfg.dim.x * cfg.dim.y * cfg.dim.z;
                 auto* dref = static_cast<float*>(dref_array.mutable_data());
                 memcpy(dref, mcx_config.exportfield, field_len * sizeof(float));
 
-                for (int i = 0; i < field_len; i++) {
-                    if (dref[i] < 0.f) {
-                        dref[i] = -dref[i];
-                        mcx_config.exportfield[i] = 0.f;
-                    } else {
-                        dref[i] = 0.f;
+                for (int ix = 0; ix < mcx_config.dim.x; ix++) {
+                    for (int iy = 0; iy < mcx_config.dim.y; iy++) {
+                        for (int iz = 0; iz < mcx_config.dim.z; iz++) {
+                            int voxelid = iz * (mcx_config.dim.x * mcx_config.dim.y) + iy * mcx_config.dim.x + ix;
+
+                            if (mcx_config.vol[voxelid]) {
+                                for (int gate = 0; gate < highdim; gate++)
+                                    for (int srcid = 0; srcid < mcx_config.srcnum; srcid++) {
+                                        dref[(gate * voxellen + voxelid) * mcx_config.srcnum + srcid] = 0.f;
+                                    }
+                            } else {
+                                for (int gate = 0; gate < highdim; gate++)
+                                    for (int srcid = 0; srcid < mcx_config.srcnum; srcid++) {
+                                        dref[(gate * voxellen + voxelid) * mcx_config.srcnum + srcid] = -dref[(gate * voxellen + voxelid) * mcx_config.srcnum + srcid];
+                                        mcx_config.exportfield[(gate * voxellen + voxelid) * mcx_config.srcnum + srcid] = 0.f;
+                                    }
+                            }
+                        }
                     }
                 }