diff --git a/ChangeLog.txt b/ChangeLog.txt index 39e0b5b2..4de5e4c2 100644 --- a/ChangeLog.txt +++ b/ChangeLog.txt @@ -2,7 +2,7 @@ Major updates are marked with a "*" -== MCX v2023.11 (Interstellar Ion - 2.2), FangQ == +== MCX v2024.1 (Interstellar Ion - 2.2), FangQ == 2023-11-07 [a710ab3] allow converting integer cfg.vol to json 2023-10-31 [d6c64e4] [test] fix rng test after make double diff --git a/NOTICE.txt b/NOTICE.txt index d4350088..923c1a23 100644 --- a/NOTICE.txt +++ b/NOTICE.txt @@ -1,6 +1,6 @@ ============================================================= Monte Carlo eXtreme (MCX) Suite - version 2023.11 + version 2024.1 ============================================================= Copyright (c) 2009-2023 Qianqian Fang diff --git a/README.md b/README.md index eb57d10b..3a386727 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ Monte Carlo eXtreme (MCX) - CUDA Edition - Author: Qianqian Fang (q.fang at neu.edu) - License: GNU General Public License version 3 (GPLv3) -- Version: 2.2.pre (v2023.11, Interstellar Ion) +- Version: 2.2.pre (v2024.1, Interstellar Ion) - Website: ![Mex and Binaries](https://github.com/fangq/mcx/actions/workflows/build_all.yml/badge.svg)\ diff --git a/README.txt b/README.txt index ca3b090e..ab8fbf7a 100644 --- a/README.txt +++ b/README.txt @@ -5,7 +5,7 @@ *Author: Qianqian Fang *License: GNU General Public License version 3 (GPLv3) -*Version: 2.2.pre (v2023.11, Interstellar Ion) +*Version: 2.2.pre (v2024.1, Interstellar Ion) *Website: http://mcx.space --------------------------------------------------------------------- diff --git a/example/multisrc/multisrc.json b/example/multisrc/multisrc.json new file mode 100644 index 00000000..e2388a4c --- /dev/null +++ b/example/multisrc/multisrc.json @@ -0,0 +1,84 @@ +{ + "Session":{ + "ID":"multisrc", + "DoAutoThread":1, + "Photons":10000000 + }, + "Forward":{ + "T0":0, + "T1":5e-09, + "Dt":5e-09 + }, + "Optode":{ + "Source":{ + "Pos":[ + [50,30,0,1], + [44.14213562373095,44.14213562373095,0,1], + [30,50,0,1], + [15.85786437626905,44.14213562373095,0,1], + [10,30.00000000000001,0,1], + [15.85786437626905,15.85786437626905,0,1], + [29.99999999999999,10,0,1], + [44.14213562373094,15.85786437626904,0,1], + [25,25,0,4] + ], + "Dir":[ + [-0.5,0,0.8660254037844386,20], + [-0.3535533905932738,-0.3535533905932737,0.8660254037844386,20], + [-1.110223024625157e-16,-0.5,0.8660254037844386,20], + [0.3535533905932737,-0.3535533905932738,0.8660254037844386,20], + [0.5,-1.387778780781446e-16,0.8660254037844386,20], + [0.3535533905932739,0.3535533905932736,0.8660254037844386,20], + [2.220446049250313e-16,0.5,0.8660254037844386,20], + [-0.3535533905932736,0.353553390593274,0.8660254037844386,20], + [0,0,1,0] + ], + "Param1":[ + [-8,3,0,0], + [-7.778174593052023,-3.535533905932737,0,0], + [-3.000000000000002,-8,0,0], + [3.535533905932736,-7.778174593052024,0,0], + [8,-3.000000000000004,0,0], + [7.778174593052025,3.535533905932735,0,0], + [3.000000000000004,7.999999999999999,0,0], + [-3.535533905932733,7.778174593052025,0,0], + [10,0,0,0] + ], + "Param2":[ + [3,3,0,0], + [4.440892098500626e-16,4.242640687119286,0,0], + [-3,3.000000000000001,0,0], + [-4.242640687119286,1.332267629550188e-15,0,0], + [-3.000000000000001,-2.999999999999999,0,0], + [-1.77635683940025e-15,-4.242640687119286,0,0], + [2.999999999999999,-3.000000000000002,0,0], + [4.242640687119286,-2.664535259100376e-15,0,0], + [0,10,0,0] + ], + "Type":"planar" + } + }, + "Domain":{ + "OriginType":0, + "Media":[ + { + "mua":0, + "mus":0, + "g":1, + "n":1 + }, + { + "mua":0.005, + "mus":0.2, + "g":0, + "n":1.37 + } + ], + "MediaFormat":"byte", + "Dim":[60,60,40], + "VolumeFile":"" + }, + "Shapes":[ + {"Grid": {"Tag":1, "Size":[60,60,40]}} + ] +} diff --git a/example/multisrc/run_multisrc.bat b/example/multisrc/run_multisrc.bat new file mode 100644 index 00000000..5be8ca8f --- /dev/null +++ b/example/multisrc/run_multisrc.bat @@ -0,0 +1 @@ +..\..\bin\mcx.exe -f multisrc.json -D P %* diff --git a/example/multisrc/run_multisrc.sh b/example/multisrc/run_multisrc.sh new file mode 100755 index 00000000..fea1dfe5 --- /dev/null +++ b/example/multisrc/run_multisrc.sh @@ -0,0 +1,2 @@ +#!/bin/sh +../../bin/mcx -f multisrc.json -D P $@ diff --git a/mcxlab/NEWS b/mcxlab/NEWS index 39e0b5b2..4de5e4c2 100644 --- a/mcxlab/NEWS +++ b/mcxlab/NEWS @@ -2,7 +2,7 @@ Major updates are marked with a "*" -== MCX v2023.11 (Interstellar Ion - 2.2), FangQ == +== MCX v2024.1 (Interstellar Ion - 2.2), FangQ == 2023-11-07 [a710ab3] allow converting integer cfg.vol to json 2023-10-31 [d6c64e4] [test] fix rng test after make double diff --git a/mcxlab/README.txt b/mcxlab/README.txt index 02a54e63..7273df56 100644 --- a/mcxlab/README.txt +++ b/mcxlab/README.txt @@ -2,7 +2,7 @@ Author: Qianqian Fang License: GNU General Public License version 3 (GPLv3) -Version: this package is part of Monte Carlo eXtreme (MCX) v2023.11 +Version: this package is part of Monte Carlo eXtreme (MCX) v2024.1 diff --git a/mcxlab/examples/demo_multsrc.m b/mcxlab/examples/demo_multsrc.m new file mode 100644 index 00000000..a0b9548d --- /dev/null +++ b/mcxlab/examples/demo_multsrc.m @@ -0,0 +1,82 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang +% +% In this example, we show how to simultaneously simulate multiple sources +% of the same type. +% +% Specifically, we create a complex source made of an array of 9 planar +% sources - where 8 of the sources form a circular array by rotating one +% of the sources around the center of the bottom-face center; the 9th +% source is manually placed at the center of the array. +% +% One can individually control the settings for each of the sources, +% including the position (srcpos), launch direction (srcdir), and +% additional parameters (srcparam1/srcparam2). + +% By default, each source entry launches an equal fraction of photons, +% determined by cfg.nphoton/size(cfg.srcpos(1,:)). However, one can set the +% 4-th element of the srcpos() for each sourc to create different weights +% for each source. Here, we set the central source 4x of the intensity +% compared to other circular-array sources (weight of 1). This allows mcx +% to allocate 4/(8+4) = 1/3 of the total photons (cfg.nphoton) to the central +% source while each other source launches 1/12 of the total photon packets. +% +% Worth to mention that we also set the 4th element of cfg.srcdir, i.e. the +% focal length, to create a convergent beam for each planar source. +% +% This file is part of Monte Carlo eXtreme (MCX) URL:http://mcx.sf.net +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% only clear cfg to avoid accidentally clearing other useful data +clear cfg cfgs + +cfg.nphoton=1e7; +cfg.vol=uint8(ones(60,60,40)); +cfg.srctype='planar'; + +Rs=20; +delta=pi/4; + +ang=0:delta:(2*pi-delta); +rotmat2d=[cos(delta), -sin(delta); sin(delta), cos(delta)]; +offset2d=[30 30]; + +% first create an array of planar sources by rotating around origin + +% the 4th element of srcpos sets the initial weight +cfg.srcpos=repmat([Rs, 0, 0, 1], length(ang), 1); + +% the 4th element of srcdir sets the focal length (positive: convergent) +cfg.srcdir=repmat([-0.5, 0, sqrt(3)/2, 20], length(ang), 1); +cfg.srcparam1=repmat([-8,3,0,0], length(ang), 1); +cfg.srcparam2=repmat([3,3,0,0], length(ang), 1); + +for i=2:length(ang) + cfg.srcpos(i,1:2)=(rotmat2d*cfg.srcpos(i-1,1:2)')'; + cfg.srcdir(i,1:2)=(rotmat2d*cfg.srcdir(i-1,1:2)')'; + cfg.srcparam1(i,1:2)=(rotmat2d*cfg.srcparam1(i-1,1:2)')'; + cfg.srcparam2(i,1:2)=(rotmat2d*cfg.srcparam2(i-1,1:2)')'; +end + +% translate the circular array to the desired offset +cfg.srcpos(:,1)=cfg.srcpos(:,1)+offset2d(1); +cfg.srcpos(:,2)=cfg.srcpos(:,2)+offset2d(2); + +% manually add the 9th source to the center + +cfg.srcpos(end+1, :)=[offset2d(1)-Rs/4,offset2d(2)-Rs/4,0,4]; % initial weight is 4 +cfg.srcdir(end+1, :)=[0, 0 1 0]; +cfg.srcparam1(end+1, :)=[Rs/2,0,0,0]; +cfg.srcparam2(end+1, :)=[0,Rs/2,0,0]; + +cfg.gpuid=1; +% cfg.gpuid='11'; % use two GPUs together +cfg.autopilot=1; +cfg.prop=[0 0 1 1;0.005 0.2 0 1.37]; +cfg.tstart=0; +cfg.tend=5e-9; +cfg.tstep=5e-9; + +flux=mcxlab(cfg); + +mcxplotvol(log10(flux.data)); diff --git a/mcxlab/mcxlab.m b/mcxlab/mcxlab.m index 510a2834..62c61aac 100644 --- a/mcxlab/mcxlab.m +++ b/mcxlab/mcxlab.m @@ -78,6 +78,11 @@ % of the srcdir direction; if the focal length is -inf, the launch % angle will be computed based on the Lambertian (cosine) distribution. % +% Starting v2024, cfg.{srcpos,srcdir,srcparam1,srcparam2} accept multiple sources, +% with each source corresponding to a single row of the array. For all 4 components, +% srcpos/srcdir support 3 or 4 columns, and srcparam1/srcparam2 support 4 columns. +% If any of the 4 compnents present, they should have matching row number. +% %== MC simulation settings == % cfg.seed: seed for the random number generator (integer) [0] % if set to a uint8 array, the binary data in each column is used diff --git a/mcxstudio/README.txt b/mcxstudio/README.txt index 1ed2e0d0..d16922ba 100644 --- a/mcxstudio/README.txt +++ b/mcxstudio/README.txt @@ -4,7 +4,7 @@ Author: Qianqian Fang License: GNU General Public License version 3 (GPLv3) -Version: 1.2.pre (v2023.11) +Version: 1.2.pre (v2024.1) Website: http://mcx.space --------------------------------------------------------------------- diff --git a/pmcx/README.md b/pmcx/README.md index d9a5a8da..add70604 100644 --- a/pmcx/README.md +++ b/pmcx/README.md @@ -4,7 +4,7 @@ - Copyright: (C) Matin Raayai Ardakani (2022-2023) , Qianqian Fang (2019-2023) , Fan-Yu Yen (2023) - License: GNU Public License V3 or later -- Version: 0.2.8 +- Version: 0.2.9 - URL: https://pypi.org/project/pmcx/ - Github: https://github.com/fangq/mcx diff --git a/pmcx/pmcx/__init__.py b/pmcx/pmcx/__init__.py index fcf58df6..44c1aa38 100644 --- 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.8" +__version__ = "0.2.9" __all__ = ( "gpuinfo", diff --git a/pmcx/setup.py b/pmcx/setup.py index c885b80c..df050c39 100644 --- 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.8", + version="0.2.9", requires=["numpy"], license="GPLv3+", author="Matin Raayai Ardakani, Qianqian Fang, Fan-Yu Yen", diff --git a/src/Makefile b/src/Makefile index 9ebd8f73..f7f86d60 100644 --- a/src/Makefile +++ b/src/Makefile @@ -74,7 +74,7 @@ NVCCOMP=$(OMP) CUDA_STATIC=--cudart static -Xcompiler "-static-libgcc -static-libstdc++" CFLAGS+=-std=c99 -CPPFLAGS+=-g -Wall #-DNO_LZMA # -DUSE_OS_TIMER +CPPFLAGS+=-g -Wall -pedantic #-DNO_LZMA # -DUSE_OS_TIMER OBJSUFFIX=.o EXESUFFIX= diff --git a/src/mcx_const.h b/src/mcx_const.h index 8b099be7..ff22311c 100644 --- a/src/mcx_const.h +++ b/src/mcx_const.h @@ -31,7 +31,7 @@ #ifndef _MCEXTREME_CONSTANT_H #define _MCEXTREME_CONSTANT_H -#define MCX_VERSION "v2023.11" +#define MCX_VERSION "v2024.1" #define MCX_VERSION_MAJOR 2 #define MCX_VERSION_MINOR 2 diff --git a/src/mcx_core.cu b/src/mcx_core.cu index 9f19b4a7..fb38f06b 100644 --- a/src/mcx_core.cu +++ b/src/mcx_core.cu @@ -1038,6 +1038,15 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* float photontof[], MCXsp* nuvox) { *w0 = 1.f; //< reuse to count for launchattempt int canfocus = 1; //< non-zero: focusable, zero: not focusable + MCXSrc* launchsrc = &(gcfg->src); + + if (gcfg->extrasrclen) { + int srcid = rand_uniform01(t) * (gcfg->extrasrclen + 1); + + if (srcid) { + launchsrc = (MCXSrc*)(gproperty + gcfg->maxmedia + 1 + gcfg->detnum + ((srcid - 1) << 2)); + } + } /** * Early termination of simulation when the detphoton buffer is filled if issavedet is set to 3 @@ -1144,17 +1153,17 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* * Attempt to launch a new photon until success */ do { - *((float4*)p) = gcfg->ps; - *((float4*)v) = gcfg->c0; + *((float4*)p) = launchsrc->pos; + *((float4*)v) = launchsrc->dir; *((float4*)f) = float4(0.f, 0.f, gcfg->minaccumtime, f->ndone); - *idx1d = gcfg->idx1dorig; - *mediaid = gcfg->mediaidorig; + *idx1d = *((uint*)&launchsrc->param2.z); /**< pre-computed 1D index of the photon at launch for pencil/isotropic beams */ + *mediaid = *((uint*)&launchsrc->param2.w); /**< pre-computed media index of the photon at launch for pencil/isotropic beams */ if (gcfg->issaveseed) { copystate(t, photonseed); } - *rv = float3(gcfg->ps.x, gcfg->ps.y, gcfg->ps.z); //< reuse as the origin of the src, needed for focusable sources + *rv = float3(launchsrc->pos.x, launchsrc->pos.y, launchsrc->pos.z); //< reuse as the origin of the src, needed for focusable sources if (issvmc) { nuvox->sv.issplit = 0; //< initialize the tissue type indicator under SVMC mode @@ -1183,23 +1192,23 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* if (gcfg->srctype == MCX_SRC_PATTERN3D) { rz = rand_uniform01(t); - *((float4*)p) = float4(p->x + rx * gcfg->srcparam1.x, - p->y + ry * gcfg->srcparam1.y, - p->z + rz * gcfg->srcparam1.z, + *((float4*)p) = float4(p->x + rx * launchsrc->param1.x, + p->y + ry * launchsrc->param1.y, + p->z + rz * launchsrc->param1.z, p->w); } else { - *((float4*)p) = float4(p->x + rx * gcfg->srcparam1.x + ry * gcfg->srcparam2.x, - p->y + rx * gcfg->srcparam1.y + ry * gcfg->srcparam2.y, - p->z + rx * gcfg->srcparam1.z + ry * gcfg->srcparam2.z, + *((float4*)p) = float4(p->x + rx * launchsrc->param1.x + ry * launchsrc->param2.x, + p->y + rx * launchsrc->param1.y + ry * launchsrc->param2.y, + p->z + rx * launchsrc->param1.z + ry * launchsrc->param2.z, p->w); } if (gcfg->srctype == MCX_SRC_PATTERN) { // need to prevent rx/ry=1 here if (gcfg->srcnum <= 1) { - p->w = gcfg->ps.w * srcpattern[(int)(ry * JUST_BELOW_ONE * gcfg->srcparam2.w) * (int)(gcfg->srcparam1.w) + (int)(rx * JUST_BELOW_ONE * gcfg->srcparam1.w)]; + p->w = launchsrc->pos.w * srcpattern[(int)(ry * JUST_BELOW_ONE * launchsrc->param2.w) * (int)(launchsrc->param1.w) + (int)(rx * JUST_BELOW_ONE * launchsrc->param1.w)]; ppath[3] = p->w; } else { - *((uint*)(ppath + 2)) = ((int)(ry * JUST_BELOW_ONE * gcfg->srcparam2.w) * (int)(gcfg->srcparam1.w) + (int)(rx * JUST_BELOW_ONE * gcfg->srcparam1.w)); + *((uint*)(ppath + 2)) = ((int)(ry * JUST_BELOW_ONE * launchsrc->param2.w) * (int)(launchsrc->param1.w) + (int)(rx * JUST_BELOW_ONE * launchsrc->param1.w)); for (int i = 0; i < gcfg->srcnum; i++) { ppath[i + 3] = srcpattern[(*((uint*)(ppath + 2))) * gcfg->srcnum + i]; @@ -1209,12 +1218,12 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* } } else if (gcfg->srctype == MCX_SRC_PATTERN3D) { if (gcfg->srcnum <= 1) { - p->w = gcfg->ps.w * srcpattern[(int)(rz * JUST_BELOW_ONE * gcfg->srcparam1.z) * (int)(gcfg->srcparam1.y) * (int)(gcfg->srcparam1.x) + - (int)(ry * JUST_BELOW_ONE * gcfg->srcparam1.y) * (int)(gcfg->srcparam1.x) + (int)(rx * JUST_BELOW_ONE * gcfg->srcparam1.x)]; + p->w = launchsrc->pos.w * srcpattern[(int)(rz * JUST_BELOW_ONE * launchsrc->param1.z) * (int)(launchsrc->param1.y) * (int)(launchsrc->param1.x) + + (int)(ry * JUST_BELOW_ONE * launchsrc->param1.y) * (int)(launchsrc->param1.x) + (int)(rx * JUST_BELOW_ONE * launchsrc->param1.x)]; ppath[3] = p->w; } else { - *((uint*)(ppath + 2)) = ((int)(rz * JUST_BELOW_ONE * gcfg->srcparam1.z) * (int)(gcfg->srcparam1.y) * (int)(gcfg->srcparam1.x) + - (int)(ry * JUST_BELOW_ONE * gcfg->srcparam1.y) * (int)(gcfg->srcparam1.x) + (int)(rx * JUST_BELOW_ONE * gcfg->srcparam1.x)); + *((uint*)(ppath + 2)) = ((int)(rz * JUST_BELOW_ONE * launchsrc->param1.z) * (int)(launchsrc->param1.y) * (int)(launchsrc->param1.x) + + (int)(ry * JUST_BELOW_ONE * launchsrc->param1.y) * (int)(launchsrc->param1.x) + (int)(rx * JUST_BELOW_ONE * launchsrc->param1.x)); for (int i = 0; i < gcfg->srcnum; i++) { ppath[i + 3] = srcpattern[(*((uint*)(ppath + 2))) * gcfg->srcnum + i]; @@ -1223,12 +1232,12 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* p->w = 1.f; } } else if (gcfg->srctype == MCX_SRC_FOURIER) - p->w = gcfg->ps.w * (cosf((floorf(gcfg->srcparam1.w) * rx + floorf(gcfg->srcparam2.w) * ry - + gcfg->srcparam1.w - floorf(gcfg->srcparam1.w)) * TWO_PI) * (1.f - gcfg->srcparam2.w + floorf(gcfg->srcparam2.w)) + 1.f) * 0.5f; //between 0 and 1 + p->w = launchsrc->pos.w * (cosf((floorf(launchsrc->param1.w) * rx + floorf(launchsrc->param2.w) * ry + + launchsrc->param1.w - floorf(launchsrc->param1.w)) * TWO_PI) * (1.f - launchsrc->param2.w + floorf(launchsrc->param2.w)) + 1.f) * 0.5f; //between 0 and 1 else if (gcfg->srctype == MCX_SRC_PENCILARRAY) { - p->x = gcfg->ps.x + floorf(rx * gcfg->srcparam1.w) * gcfg->srcparam1.x / (gcfg->srcparam1.w - 1.f) + floorf(ry * gcfg->srcparam2.w) * gcfg->srcparam2.x / (gcfg->srcparam2.w - 1.f); - p->y = gcfg->ps.y + floorf(rx * gcfg->srcparam1.w) * gcfg->srcparam1.y / (gcfg->srcparam1.w - 1.f) + floorf(ry * gcfg->srcparam2.w) * gcfg->srcparam2.y / (gcfg->srcparam2.w - 1.f); - p->z = gcfg->ps.z + floorf(rx * gcfg->srcparam1.w) * gcfg->srcparam1.z / (gcfg->srcparam1.w - 1.f) + floorf(ry * gcfg->srcparam2.w) * gcfg->srcparam2.z / (gcfg->srcparam2.w - 1.f); + p->x = launchsrc->pos.x + floorf(rx * launchsrc->param1.w) * launchsrc->param1.x / (launchsrc->param1.w - 1.f) + floorf(ry * launchsrc->param2.w) * launchsrc->param2.x / (launchsrc->param2.w - 1.f); + p->y = launchsrc->pos.y + floorf(rx * launchsrc->param1.w) * launchsrc->param1.y / (launchsrc->param1.w - 1.f) + floorf(ry * launchsrc->param2.w) * launchsrc->param2.y / (launchsrc->param2.w - 1.f); + p->z = launchsrc->pos.z + floorf(rx * launchsrc->param1.w) * launchsrc->param1.z / (launchsrc->param1.w - 1.f) + floorf(ry * launchsrc->param2.w) * launchsrc->param2.z / (launchsrc->param2.w - 1.f); } *idx1d = (int(floorf(p->z)) * gcfg->dimlen.y + int(floorf(p->y)) * gcfg->dimlen.x + int(floorf(p->x))); @@ -1239,9 +1248,9 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* *mediaid = media[*idx1d]; } - *rv = float3(rv->x + (gcfg->srcparam1.x + gcfg->srcparam2.x) * 0.5f, - rv->y + (gcfg->srcparam1.y + gcfg->srcparam2.y) * 0.5f, - rv->z + (gcfg->srcparam1.z + gcfg->srcparam2.z) * 0.5f); + *rv = float3(rv->x + (launchsrc->param1.x + launchsrc->param2.x) * 0.5f, + rv->y + (launchsrc->param1.y + launchsrc->param2.y) * 0.5f, + rv->z + (launchsrc->param1.z + launchsrc->param2.z) * 0.5f); break; } @@ -1249,21 +1258,21 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* case (MCX_SRC_FOURIERX2D): { // [v1x][v1y][v1z][|v2|]; [kx][ky][phi0][M], unit(v0) x unit(v1)=unit(v2) float rx = rand_uniform01(t); float ry = rand_uniform01(t); - float4 v2 = gcfg->srcparam1; + float4 v2 = launchsrc->param1; // calculate v2 based on v2=|v2| * unit(v0) x unit(v1) - v2.w *= rsqrtf(gcfg->srcparam1.x * gcfg->srcparam1.x + gcfg->srcparam1.y * gcfg->srcparam1.y + gcfg->srcparam1.z * gcfg->srcparam1.z); - v2.x = v2.w * (gcfg->c0.y * gcfg->srcparam1.z - gcfg->c0.z * gcfg->srcparam1.y); - v2.y = v2.w * (gcfg->c0.z * gcfg->srcparam1.x - gcfg->c0.x * gcfg->srcparam1.z); - v2.z = v2.w * (gcfg->c0.x * gcfg->srcparam1.y - gcfg->c0.y * gcfg->srcparam1.x); - *((float4*)p) = float4(p->x + rx * gcfg->srcparam1.x + ry * v2.x, - p->y + rx * gcfg->srcparam1.y + ry * v2.y, - p->z + rx * gcfg->srcparam1.z + ry * v2.z, + v2.w *= rsqrtf(launchsrc->param1.x * launchsrc->param1.x + launchsrc->param1.y * launchsrc->param1.y + launchsrc->param1.z * launchsrc->param1.z); + v2.x = v2.w * (launchsrc->dir.y * launchsrc->param1.z - launchsrc->dir.z * launchsrc->param1.y); + v2.y = v2.w * (launchsrc->dir.z * launchsrc->param1.x - launchsrc->dir.x * launchsrc->param1.z); + v2.z = v2.w * (launchsrc->dir.x * launchsrc->param1.y - launchsrc->dir.y * launchsrc->param1.x); + *((float4*)p) = float4(p->x + rx * launchsrc->param1.x + ry * v2.x, + p->y + rx * launchsrc->param1.y + ry * v2.y, + p->z + rx * launchsrc->param1.z + ry * v2.z, p->w); if (gcfg->srctype == MCX_SRC_FOURIERX2D) { - p->w = gcfg->ps.w * (sinf((gcfg->srcparam2.x * rx + gcfg->srcparam2.z) * TWO_PI) * sinf((gcfg->srcparam2.y * ry + gcfg->srcparam2.w) * TWO_PI) + 1.f) * 0.5f; //between 0 and 1 + p->w = launchsrc->pos.w * (sinf((launchsrc->param2.x * rx + launchsrc->param2.z) * TWO_PI) * sinf((launchsrc->param2.y * ry + launchsrc->param2.w) * TWO_PI) + 1.f) * 0.5f; //between 0 and 1 } else { - p->w = gcfg->ps.w * (cosf((gcfg->srcparam2.x * rx + gcfg->srcparam2.y * ry + gcfg->srcparam2.z) * TWO_PI) * (1.f - gcfg->srcparam2.w) + 1.f) * 0.5f; //between 0 and 1 + p->w = launchsrc->pos.w * (cosf((launchsrc->param2.x * rx + launchsrc->param2.y * ry + launchsrc->param2.z) * TWO_PI) * (1.f - launchsrc->param2.w) + 1.f) * 0.5f; //between 0 and 1 } *idx1d = (int(floorf(p->z)) * gcfg->dimlen.y + int(floorf(p->y)) * gcfg->dimlen.x + int(floorf(p->x))); @@ -1274,9 +1283,9 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* *mediaid = media[*idx1d]; } - *rv = float3(rv->x + (gcfg->srcparam1.x + v2.x) * 0.5f, - rv->y + (gcfg->srcparam1.y + v2.y) * 0.5f, - rv->z + (gcfg->srcparam1.z + v2.z) * 0.5f); + *rv = float3(rv->x + (launchsrc->param1.x + v2.x) * 0.5f, + rv->y + (launchsrc->param1.y + v2.y) * 0.5f, + rv->z + (launchsrc->param1.z + v2.z) * 0.5f); break; } @@ -1287,8 +1296,8 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* // http://mathworld.wolfram.com/DiskPointPicking.html float phi, sphi, cphi, r; - if (gcfg->srctype == MCX_SRC_RING && (gcfg->srcparam1.z > 0.f || gcfg->srcparam1.w > 0.f)) { - phi = fabsf(gcfg->srcparam1.z - gcfg->srcparam1.w) * rand_uniform01(t) + fminf(gcfg->srcparam1.z, gcfg->srcparam1.w); + if (gcfg->srctype == MCX_SRC_RING && (launchsrc->param1.z > 0.f || launchsrc->param1.w > 0.f)) { + phi = fabsf(launchsrc->param1.z - launchsrc->param1.w) * rand_uniform01(t) + fminf(launchsrc->param1.z, launchsrc->param1.w); } else { phi = TWO_PI * rand_uniform01(t); } @@ -1296,12 +1305,12 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* sincosf(phi, &sphi, &cphi); if (gcfg->srctype == MCX_SRC_DISK || gcfg->srctype == MCX_SRC_RING) { - r = sqrtf(rand_uniform01(t) * fabsf(gcfg->srcparam1.x * gcfg->srcparam1.x - gcfg->srcparam1.y * gcfg->srcparam1.y) + gcfg->srcparam1.y * gcfg->srcparam1.y); - } else if (fabsf(gcfg->c0.w) < 1e-5f || fabsf(gcfg->srcparam1.y) < 1e-5f) { - r = sqrtf(0.5f * rand_next_scatlen(t)) * gcfg->srcparam1.x; + r = sqrtf(rand_uniform01(t) * fabsf(launchsrc->param1.x * launchsrc->param1.x - launchsrc->param1.y * launchsrc->param1.y) + launchsrc->param1.y * launchsrc->param1.y); + } else if (fabsf(launchsrc->dir.w) < 1e-5f || fabsf(launchsrc->param1.y) < 1e-5f) { + r = sqrtf(0.5f * rand_next_scatlen(t)) * launchsrc->param1.x; } else { - r = gcfg->srcparam1.x * gcfg->srcparam1.x * M_PI / gcfg->srcparam1.y; //Rayleigh range - r = sqrtf(0.5f * rand_next_scatlen(t) * (1.f + (gcfg->c0.w * gcfg->c0.w / (r * r)))) * gcfg->srcparam1.x; + r = launchsrc->param1.x * launchsrc->param1.x * M_PI / launchsrc->param1.y; //Rayleigh range + r = sqrtf(0.5f * rand_next_scatlen(t) * (1.f + (launchsrc->dir.w * launchsrc->dir.w / (r * r)))) * launchsrc->param1.x; } if ( v->z > -1.f + EPS && v->z < 1.f - EPS ) { @@ -1341,8 +1350,8 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* sincosf(ang, &sphi, &cphi); if (gcfg->srctype == MCX_SRC_CONE) { // a solid-angle section of a uniform sphere - ang = cosf(gcfg->srcparam1.x); - ang = (gcfg->srcparam1.y > 0.f) ? rand_uniform01(t) * gcfg->srcparam1.x : acos(rand_uniform01(t) * (1.0 - ang) + ang); //sine distribution + ang = cosf(launchsrc->param1.x); + ang = (launchsrc->param1.y > 0.f) ? rand_uniform01(t) * launchsrc->param1.x : acos(rand_uniform01(t) * (1.0 - ang) + ang); //sine distribution } else { if (gcfg->srctype == MCX_SRC_ISOTROPIC) { // uniform sphere ang = acosf(2.f * rand_uniform01(t) - 1.f); //sine distribution @@ -1361,7 +1370,7 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* float ang, stheta, ctheta, sphi, cphi; ang = TWO_PI * rand_uniform01(t); //next arimuth angle sincosf(ang, &sphi, &cphi); - ang = sqrtf(2.f * rand_next_scatlen(t)) * (1.f - 2.f * rand_uniform01(t)) * gcfg->srcparam1.x; + ang = sqrtf(2.f * rand_next_scatlen(t)) * (1.f - 2.f * rand_uniform01(t)) * launchsrc->param1.x; sincosf(ang, &stheta, &ctheta); rotatevector(v, stheta, ctheta, sphi, cphi); canfocus = 0; @@ -1373,15 +1382,15 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* float phi = TWO_PI * rand_uniform01(t); sincosf(phi, &sphi, &cphi); - float r = sqrtf(0.5f * rand_next_scatlen(t)) * gcfg->srcparam1.x; + float r = sqrtf(0.5f * rand_next_scatlen(t)) * launchsrc->param1.x; /** parameter to generate photon path from coordinates at focus (depends on focal distance and rayleigh range) */ - float tt = -gcfg->srcparam1.y / gcfg->srcparam1.z; - float l = rsqrtf(r * r + gcfg->srcparam1.z * gcfg->srcparam1.z); + float tt = -launchsrc->param1.y / launchsrc->param1.z; + float l = rsqrtf(r * r + launchsrc->param1.z * launchsrc->param1.z); /** if beam direction is along +z or -z direction */ float3 pd = float3(r * (cphi - tt * sphi), r * (sphi + tt * cphi), 0.f); // position displacement from srcpos - float3 v0 = float3(-r * sphi * l, r * cphi * l, gcfg->srcparam1.z * l); // photon dir. w.r.t the beam dir. v + float3 v0 = float3(-r * sphi * l, r * cphi * l, launchsrc->param1.z * l); // photon dir. w.r.t the beam dir. v /** if beam dir. is not +z or -z, compute photon position and direction after rotation */ if ( v->z > -1.f + EPS && v->z < 1.f - EPS ) { @@ -1422,23 +1431,23 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* case (MCX_SRC_LINE): // uniformally emitting line source, emitting cylindrically case (MCX_SRC_SLIT): { // a line source emitting only along a specified direction, like a light sheet float r = rand_uniform01(t); - *((float4*)p) = float4(p->x + r * gcfg->srcparam1.x, - p->y + r * gcfg->srcparam1.y, - p->z + r * gcfg->srcparam1.z, + *((float4*)p) = float4(p->x + r * launchsrc->param1.x, + p->y + r * launchsrc->param1.y, + p->z + r * launchsrc->param1.z, p->w); if (gcfg->srctype == MCX_SRC_LINE) { float sphi, cphi; - r = rsqrtf(gcfg->srcparam1.x * gcfg->srcparam1.x + gcfg->srcparam1.y * gcfg->srcparam1.y + gcfg->srcparam1.z * gcfg->srcparam1.z); - *((float4*)v) = float4(gcfg->srcparam1.x * r, gcfg->srcparam1.y * r, gcfg->srcparam1.z * r, v->nscat); + r = rsqrtf(launchsrc->param1.x * launchsrc->param1.x + launchsrc->param1.y * launchsrc->param1.y + launchsrc->param1.z * launchsrc->param1.z); + *((float4*)v) = float4(launchsrc->param1.x * r, launchsrc->param1.y * r, launchsrc->param1.z * r, v->nscat); r = TWO_PI * rand_uniform01(t); // phi sincosf(r, &sphi, &cphi); // y=sin(phi), x=cos(phi) rotatevector(v, 1.f, 0.f, sphi, cphi); } - *rv = float3(rv->x + (gcfg->srcparam1.x) * 0.5f, - rv->y + (gcfg->srcparam1.y) * 0.5f, - rv->z + (gcfg->srcparam1.z) * 0.5f); + *rv = float3(rv->x + (launchsrc->param1.x) * 0.5f, + rv->y + (launchsrc->param1.y) * 0.5f, + rv->z + (launchsrc->param1.z) * 0.5f); canfocus = (gcfg->srctype == MCX_SRC_SLIT); break; } @@ -1455,7 +1464,7 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* float ang, stheta, ctheta, sphi, cphi; - if (gcfg->c0.w > 0.f) { // if focal-length > 0, no interpolation, just read the angleinvcdf value + if (launchsrc->dir.w > 0.f) { // if focal-length > 0, no interpolation, just read the angleinvcdf value ang = fminf(rand_uniform01(t) * gcfg->nangle, gcfg->nangle - EPS); cphi = ((float*)(sharedmem))[(int)(ang) + gcfg->nphaselen]; } else { // odd number length, interpolate between neigboring values @@ -1470,8 +1479,8 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* ang = TWO_PI * rand_uniform01(t); //next arimuth angle sincosf(ang, &sphi, &cphi); - if (gcfg->c0.w < 1.5f && gcfg->c0.w >= 0.f) { - *((float4*)v) = gcfg->c0; + if (launchsrc->dir.w < 1.5f && launchsrc->dir.w >= 0.f) { + *((float4*)v) = launchsrc->dir; } rotatevector(v, stheta, ctheta, sphi, cphi); @@ -1480,25 +1489,25 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* * If beam focus is set, determine the incident angle */ - if (isnan(gcfg->c0.w)) { // isotropic if focal length is nan + if (isnan(launchsrc->dir.w)) { // isotropic if focal length is nan float ang, stheta, ctheta, sphi, cphi; ang = TWO_PI * rand_uniform01(t); //next arimuth angle sincosf(ang, &sphi, &cphi); ang = acosf(2.f * rand_uniform01(t) - 1.f); //sine distribution sincosf(ang, &stheta, &ctheta); rotatevector(v, stheta, ctheta, sphi, cphi); - } else if (gcfg->c0.w < 0.f && isinf(gcfg->c0.w)) { // lambertian (cosine distribution) if focal length is -inf + } else if (launchsrc->dir.w < 0.f && isinf(launchsrc->dir.w)) { // lambertian (cosine distribution) if focal length is -inf float ang, stheta, ctheta, sphi, cphi; ang = TWO_PI * rand_uniform01(t); //next arimuth angle sincosf(ang, &sphi, &cphi); stheta = sqrtf(rand_uniform01(t)); ctheta = sqrtf(1.f - stheta * stheta); rotatevector(v, stheta, ctheta, sphi, cphi); - } else if (gcfg->c0.w != 0.f) { - float Rn2 = (gcfg->c0.w > 0.f) - (gcfg->c0.w < 0.f); - rv->x += gcfg->c0.w * v->x; - rv->y += gcfg->c0.w * v->y; - rv->z += gcfg->c0.w * v->z; + } else if (launchsrc->dir.w != 0.f) { + float Rn2 = (launchsrc->dir.w > 0.f) - (launchsrc->dir.w < 0.f); + rv->x += launchsrc->dir.w * v->x; + rv->y += launchsrc->dir.w * v->y; + rv->z += launchsrc->dir.w * v->z; v->x = Rn2 * (rv->x - p->x); v->y = Rn2 * (rv->y - p->y); v->z = Rn2 * (rv->z - p->z); @@ -1651,7 +1660,7 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[], uint idx1d, idx1dold; //< linear index to the current voxel in the media array - uint mediaid = gcfg->mediaidorig; + uint mediaid = *((uint*)(&gcfg->src.param2.w)); uint mediaidold = 0; int isdet = 0; float n1; //< reflection var @@ -2671,10 +2680,10 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { /** \c param - constants to be used in the GPU, copied to GPU as \c gcfg, stored in the constant memory */ MCXParam param = {cfg->steps, minstep, 0, 0, cfg->tend, R_C0* cfg->unitinmm, (uint)cfg->issave2pt, (uint)cfg->isreflect, (uint)cfg->isrefint, (uint)cfg->issavedet, 1.f / cfg->tstep, - p0, c0, s0, maxidx, uint4(0, 0, 0, 0), cp0, cp1, uint2(0, 0), cfg->minenergy, + p0, c0, cfg->srcparam1, cfg->srcparam2, cfg->extrasrclen, s0, maxidx, uint4(0, 0, 0, 0), cp0, cp1, uint2(0, 0), cfg->minenergy, cfg->sradius* cfg->sradius, minstep* R_C0* cfg->unitinmm, cfg->srctype, - cfg->srcparam1, cfg->srcparam2, cfg->voidtime, cfg->maxdetphoton, - cfg->medianum - 1, cfg->detnum, cfg->polmedianum, cfg->maxgate, 0, 0, ABS(cfg->sradius + 2.f) < EPS /*isatomic*/, + cfg->voidtime, cfg->maxdetphoton, + cfg->medianum - 1, cfg->detnum, cfg->polmedianum, cfg->maxgate, ABS(cfg->sradius + 2.f) < EPS /*isatomic*/, (uint)cfg->maxvoidstep, cfg->issaveseed > 0, (uint)cfg->issaveref, cfg->isspecular > 0, cfg->maxdetphoton * hostdetreclen, cfg->seed, (uint)cfg->outputtype, 0, 0, cfg->faststep, cfg->debuglevel, cfg->savedetflag, hostdetreclen, partialdata, w0offset, cfg->mediabyte, @@ -2744,7 +2753,7 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { /** If domain is 2D, the initial photon launch direction in the singleton dimension is expected to be 0 */ if (is2d) { - float* vec = &(param.c0.x); + float* vec = &(param.src.dir.x); if (ABS(vec[is2d - 1]) > EPS) { mcx_error(-1, "input domain is 2D, the initial direction can not have non-zero value in the singular dimension", __FILE__, __LINE__); @@ -3051,17 +3060,6 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { param.dimlen = dimlen; param.cachebox = cachebox; - /** - * Additional constants to avoid repeated computation inside GPU - */ - if (p0.x < 0.f || p0.y < 0.f || p0.z < 0.f || p0.x >= cfg->dim.x || p0.y >= cfg->dim.y || p0.z >= cfg->dim.z) { - param.idx1dorig = 0; - param.mediaidorig = 0; - } else { - param.idx1dorig = (int(floorf(p0.z)) * dimlen.y + int(floorf(p0.y)) * dimlen.x + int(floorf(p0.x))); - param.mediaidorig = (cfg->vol[param.idx1dorig] & MED_MASK); - } - memcpy(&(param.bc), cfg->bc, 12); Vvox = cfg->steps.x * cfg->steps.y * cfg->steps.z; /*Vvox: voxel volume in mm^3*/ @@ -3128,6 +3126,10 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { CUDA_ASSERT(cudaMemcpyToSymbol(gproperty, cfg->prop, cfg->medianum * sizeof(Medium), 0, cudaMemcpyHostToDevice)); CUDA_ASSERT(cudaMemcpyToSymbol(gproperty, cfg->detpos, cfg->detnum * sizeof(float4), cfg->medianum * sizeof(Medium), cudaMemcpyHostToDevice)); + if (cfg->srcdata) { + CUDA_ASSERT(cudaMemcpyToSymbol(gproperty, cfg->srcdata, cfg->extrasrclen * 4 * sizeof(float4), cfg->medianum * sizeof(Medium) + cfg->detnum * sizeof(float4), cudaMemcpyHostToDevice)); + } + MCX_FPRINTF(cfg->flog, "init complete : %d ms\n", GetTimeMillis() - tic); /** diff --git a/src/mcx_core.h b/src/mcx_core.h index 847474f6..b23be05d 100644 --- a/src/mcx_core.h +++ b/src/mcx_core.h @@ -116,6 +116,13 @@ typedef union __align__(16) GProperty { float f[4]; } Gprop; +typedef struct __align__(16) MCXSource { + float4 pos; /**< initial position vector, for pencil beam */ + float4 dir; /**< initial directon vector, for pencil beam */ + float4 param1; /**< source parameters set 1 */ + float4 param2; /**< source parameters set 2 */ +} MCXSrc; + typedef unsigned char uchar; /** @@ -136,8 +143,8 @@ typedef struct __align__(16) KernelParams { unsigned int dorefint; /**< flag if mcx perform reflection calculations at internal boundaries */ unsigned int savedet; /**< flag if mcx outputs detected photon partial length data */ float Rtstep; /**< reciprocal of the step size */ - float4 ps; /**< initial position vector, for pencil beam */ - float4 c0; /**< initial directon vector, for pencil beam */ + MCXSrc src; /**< additional source data, including pos, dir, param1, param2 */ + unsigned int extrasrclen; /**< number of additional sources */ float4 s0; /**< initial stokes parameters, for polarized photon simulation */ float3 maxidx; /**< maximum index in x/y/z directions for out-of-bound tests */ uint4 dimlen; /**< maximum index used to convert x/y/z to 1D array index */ @@ -148,16 +155,12 @@ typedef struct __align__(16) KernelParams { float skipradius2; /**< square of the radius within which the data is cached (obsolete) */ float minaccumtime; /**< time steps for tMCimg like weight accummulation (obsolete) */ int srctype; /**< type of the source */ - float4 srcparam1; /**< source parameters set 1 */ - float4 srcparam2; /**< source parameters set 2 */ int voidtime; /**< flag if the time-of-flight in the background is counted */ unsigned int maxdetphoton; /**< max number of detected photons */ unsigned int maxmedia; /**< max number of media labels */ unsigned int detnum; /**< max number of detectors */ unsigned int maxpolmedia; /**< max number of media labels for polarized light */ unsigned int maxgate; /**< max number of time gates */ - unsigned int idx1dorig; /**< pre-computed 1D index of the photon at launch for pencil/isotropic beams */ - unsigned int mediaidorig; /**< pre-computed media index of the photon at launch for pencil/isotropic beams */ unsigned int isatomic; /**< whether atomic operations are used */ unsigned int maxvoidstep; /**< max steps that photon can travel in the background before entering non-zero voxels */ unsigned int issaveseed; /**< flag if one need to save the detected photon seeds for replay */ diff --git a/src/mcx_utils.c b/src/mcx_utils.c index b230a4bc..a556e1f8 100644 --- a/src/mcx_utils.c +++ b/src/mcx_utils.c @@ -332,6 +332,8 @@ void mcx_initcfg(Config* cfg) { cfg->invcdf = NULL; cfg->nangle = 0; cfg->angleinvcdf = NULL; + cfg->extrasrclen = 0; + cfg->srcdata = NULL; memset(cfg->jsonfile, 0, MAX_PATH_LENGTH); memset(cfg->bc, 0, 13); memset(&(cfg->srcparam1), 0, sizeof(float4)); @@ -454,6 +456,10 @@ void mcx_clearcfg(Config* cfg) { free(cfg->angleinvcdf); } + if (cfg->srcdata) { + free(cfg->srcdata); + } + mcx_initcfg(cfg); } @@ -1161,7 +1167,7 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c */ void mcx_printlog(Config* cfg, char* str) { - if (cfg->flog > 0) { /*stdout is 1*/ + if (cfg->flog) { /*stdout is 1*/ MCX_FPRINTF(cfg->flog, "%s\n", str); } } @@ -1599,6 +1605,43 @@ void mcx_preprocess(Config* cfg) { } } + if (cfg->srcpos.w == 0.f) { + cfg->srcpos.w = 1.f; + } + + if (cfg->extrasrclen) { + for (int i = 0; i < cfg->extrasrclen; i++) { + if (cfg->srcdata[i].srcpos.w == 0.f) { + cfg->srcdata[i].srcpos.w = 1.f; + } + } + } + + // for all point sources, precompute launch voxel index and media value and store those in srcparam2.z/.w internally + if (cfg->srctype <= MCX_SRC_CONE || cfg->srctype == MCX_SRC_ARCSINE || cfg->srctype == MCX_SRC_ZGAUSSIAN) { + if (cfg->srcpos.x < 0.f || cfg->srcpos.y < 0.f || cfg->srcpos.z < 0.f || cfg->srcpos.x >= cfg->dim.x || cfg->srcpos.y >= cfg->dim.y || cfg->srcpos.z >= cfg->dim.z) { + *((uint*)&cfg->srcparam2.z) = 0; + *((uint*)&cfg->srcparam2.w) = 0; + } else { + uint idx1dorig = ((int)(floorf(cfg->srcpos.z)) * (cfg->dim.y * cfg->dim.x) + (int)(floorf(cfg->srcpos.y)) * cfg->dim.x + (int)(floorf(cfg->srcpos.x))); + *((uint*)&cfg->srcparam2.z) = idx1dorig; + *((uint*)&cfg->srcparam2.w) = (cfg->vol[idx1dorig] & MED_MASK); + } + + if (cfg->extrasrclen) { + for (int i = 0; i < cfg->extrasrclen; i++) { + if (cfg->srcdata[i].srcpos.x < 0.f || cfg->srcdata[i].srcpos.y < 0.f || cfg->srcdata[i].srcpos.z < 0.f || cfg->srcdata[i].srcpos.x >= cfg->dim.x || cfg->srcdata[i].srcpos.y >= cfg->dim.y || cfg->srcdata[i].srcpos.z >= cfg->dim.z) { + *((uint*)&cfg->srcdata[i].srcparam2.z) = 0; + *((uint*)&cfg->srcdata[i].srcparam2.w) = 0; + } else { + uint idx1dorig = ((int)(floorf(cfg->srcdata[i].srcpos.z)) * (cfg->dim.y * cfg->dim.x) + (int)(floorf(cfg->srcdata[i].srcpos.y)) * cfg->dim.x + (int)(floorf(cfg->srcdata[i].srcpos.x))); + *((uint*)&cfg->srcdata[i].srcparam2.z) = idx1dorig; + *((uint*)&cfg->srcdata[i].srcparam2.w) = (cfg->vol[idx1dorig] & MED_MASK); + } + } + } + } + if (cfg->vol) { unsigned int dimxyz = cfg->dim.x * cfg->dim.y * cfg->dim.z; @@ -1850,7 +1893,7 @@ void mcx_loadconfig(FILE* in, Config* cfg) { if (cfg->seed == 0) { MCX_ASSERT(fscanf(in, "%d", &(cfg->seed) ) == 1); } else { - MCX_ASSERT(fscanf(in, "%d", &itmp ) == 1); + MCX_ASSERT(fscanf(in, "%u", &itmp ) == 1); } comm = fgets(comment, MAX_PATH_LENGTH, in); @@ -1862,7 +1905,7 @@ void mcx_loadconfig(FILE* in, Config* cfg) { MCX_ASSERT(fscanf(in, "%f %f %f", &(cfg->srcpos.x), &(cfg->srcpos.y), &(cfg->srcpos.z) ) == 3); comm = fgets(comment, MAX_PATH_LENGTH, in); - if (cfg->issrcfrom0 == 0 && comm != NULL && sscanf(comm, "%d", &itmp) == 1) { + if (cfg->issrcfrom0 == 0 && comm != NULL && sscanf(comm, "%u", &itmp) == 1) { cfg->issrcfrom0 = itmp; } @@ -1923,21 +1966,21 @@ void mcx_loadconfig(FILE* in, Config* cfg) { fprintf(stdout, "%s\nPlease specify the x voxel size (in mm), x dimension, min and max x-index [1.0 100 1 100]:\n\t", filename); } - MCX_ASSERT(fscanf(in, "%f %d %d %d", &(cfg->steps.x), &(cfg->dim.x), &(cfg->crop0.x), &(cfg->crop1.x)) == 4); + MCX_ASSERT(fscanf(in, "%f %u %u %u", &(cfg->steps.x), &(cfg->dim.x), &(cfg->crop0.x), &(cfg->crop1.x)) == 4); comm = fgets(comment, MAX_PATH_LENGTH, in); if (in == stdin) fprintf(stdout, "%f %d %d %d\nPlease specify the y voxel size (in mm), y dimension, min and max y-index [1.0 100 1 100]:\n\t", cfg->steps.x, cfg->dim.x, cfg->crop0.x, cfg->crop1.x); - MCX_ASSERT(fscanf(in, "%f %d %d %d", &(cfg->steps.y), &(cfg->dim.y), &(cfg->crop0.y), &(cfg->crop1.y)) == 4); + MCX_ASSERT(fscanf(in, "%f %u %u %u", &(cfg->steps.y), &(cfg->dim.y), &(cfg->crop0.y), &(cfg->crop1.y)) == 4); comm = fgets(comment, MAX_PATH_LENGTH, in); if (in == stdin) fprintf(stdout, "%f %d %d %d\nPlease specify the z voxel size (in mm), z dimension, min and max z-index [1.0 100 1 100]:\n\t", cfg->steps.y, cfg->dim.y, cfg->crop0.y, cfg->crop1.y); - MCX_ASSERT(fscanf(in, "%f %d %d %d", &(cfg->steps.z), &(cfg->dim.z), &(cfg->crop0.z), &(cfg->crop1.z)) == 4); + MCX_ASSERT(fscanf(in, "%f %u %u %u", &(cfg->steps.z), &(cfg->dim.z), &(cfg->crop0.z), &(cfg->crop1.z)) == 4); comm = fgets(comment, MAX_PATH_LENGTH, in); if (cfg->steps.x != cfg->steps.y || cfg->steps.y != cfg->steps.z) { @@ -1976,7 +2019,7 @@ void mcx_loadconfig(FILE* in, Config* cfg) { fprintf(stdout, "%f %d %d %d\nPlease specify the total types of media:\n\t", cfg->steps.z, cfg->dim.z, cfg->crop0.z, cfg->crop1.z); - MCX_ASSERT(fscanf(in, "%d", &(cfg->medianum)) == 1); + MCX_ASSERT(fscanf(in, "%u", &(cfg->medianum)) == 1); cfg->medianum++; comm = fgets(comment, MAX_PATH_LENGTH, in); @@ -2007,15 +2050,15 @@ void mcx_loadconfig(FILE* in, Config* cfg) { fprintf(stdout, "Please specify the total number of detectors and fiber diameter (in grid unit):\n\t"); } - MCX_ASSERT(fscanf(in, "%d %f", &(cfg->detnum), &(cfg->detradius)) == 2); + MCX_ASSERT(fscanf(in, "%u %f", &(cfg->detnum), &(cfg->detradius)) == 2); comm = fgets(comment, MAX_PATH_LENGTH, in); if (in == stdin) { fprintf(stdout, "%d %f\n", cfg->detnum, cfg->detradius); } - if (cfg->medianum + cfg->detnum > MAX_PROP_AND_DETECTORS) { - MCX_ERROR(-4, "input media types plus detector number exceeds the maximum total (4000)"); + if (cfg->medianum + cfg->detnum + (cfg->extrasrclen << 2) > MAX_PROP_AND_DETECTORS) { + MCX_ERROR(-4, "input media types plus detector number plus additional sources exceeds the maximum total (4000)"); } cfg->detpos = (float4*)malloc(sizeof(float4) * cfg->detnum); @@ -2492,22 +2535,77 @@ int mcx_loadjson(cJSON* root, Config* cfg) { if (src) { subitem = FIND_JSON_OBJ("Pos", "Optode.Source.Pos", src); - if (subitem) { + if (subitem && cJSON_IsArray(subitem)) { + cfg->extrasrclen = 0; + + if (cJSON_IsArray(subitem->child)) { + if (cfg->srcdata && cfg->extrasrclen != cJSON_GetArraySize(subitem) - 1) { + MCX_ERROR(-1, "Length of sub-elements of Pos/Dir/Param1/Param2 must match"); + } + + if (cfg->srcdata == NULL) { + cfg->extrasrclen = cJSON_GetArraySize(subitem) - 1; + } + + subitem = subitem->child; + } + cfg->srcpos.x = subitem->child->valuedouble; cfg->srcpos.y = subitem->child->next->valuedouble; cfg->srcpos.z = subitem->child->next->next->valuedouble; + + if (subitem->child->next->next->next) { + cfg->srcpos.w = subitem->child->next->next->next->valuedouble; + } + + if (cfg->extrasrclen > 0) { + int count = 0; + + if (cfg->srcdata == NULL) { + cfg->srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), cfg->extrasrclen); + } + + while (subitem->next && count < cfg->extrasrclen) { + + subitem = subitem->next; + cfg->srcdata[count].srcpos.x = subitem->child->valuedouble - (!cfg->issrcfrom0); + cfg->srcdata[count].srcpos.y = subitem->child->next->valuedouble - (!cfg->issrcfrom0); + cfg->srcdata[count].srcpos.z = subitem->child->next->next->valuedouble - (!cfg->issrcfrom0); + + if (subitem->child->next->next->next) { + cfg->srcdata[count].srcpos.w = subitem->child->next->next->next->valuedouble; + } + + count++; + } + } } - cfg->srcpos.w = FIND_JSON_KEY("Weight", "Optode.Source.Weight", src, 1.f, valuedouble); + if (FIND_JSON_OBJ("Weight", "Optode.Source.Weight", src)) { + cfg->srcpos.w = FIND_JSON_KEY("Weight", "Optode.Source.Weight", src, cfg->srcpos.w, valuedouble); + } subitem = FIND_JSON_OBJ("Dir", "Optode.Source.Dir", src); - if (subitem) { + if (subitem && cJSON_IsArray(subitem)) { + if (cJSON_IsArray(subitem->child)) { + if (cfg->srcdata && cfg->extrasrclen != cJSON_GetArraySize(subitem) - 1) { + MCX_ERROR(-1, "Length of sub-elements of Pos/Dir/Param1/Param2 must match"); + } + + if (cfg->srcdata == NULL) { + cfg->extrasrclen = cJSON_GetArraySize(subitem) - 1; + } + + subitem = subitem->child; + } + cfg->srcdir.x = subitem->child->valuedouble; cfg->srcdir.y = subitem->child->next->valuedouble; cfg->srcdir.z = subitem->child->next->next->valuedouble; if (subitem->child->next->next->next) { + if (cJSON_IsString(subitem->child->next->next->next)) { if (strcmp(subitem->child->next->next->next->valuestring, "_NaN_") == 0) { cfg->srcdir.w = NAN; @@ -2520,6 +2618,38 @@ int mcx_loadjson(cJSON* root, Config* cfg) { cfg->srcdir.w = subitem->child->next->next->next->valuedouble; } } + + if (cfg->extrasrclen > 0) { + int count = 0; + + if (cfg->srcdata == NULL) { + cfg->srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), cfg->extrasrclen); + } + + while (subitem->next && count < cfg->extrasrclen) { + + subitem = subitem->next; + cfg->srcdata[count].srcdir.x = subitem->child->valuedouble; + cfg->srcdata[count].srcdir.y = subitem->child->next->valuedouble; + cfg->srcdata[count].srcdir.z = subitem->child->next->next->valuedouble; + + if (subitem->child->next->next->next) { + if (cJSON_IsString(subitem->child->next->next->next)) { + if (strcmp(subitem->child->next->next->next->valuestring, "_NaN_") == 0) { + cfg->srcdata[count].srcdir.w = NAN; + } else if (strcmp(subitem->child->next->next->next->valuestring, "_Inf_") == 0) { + cfg->srcdata[count].srcdir.w = INFINITY; + } else if (strcmp(subitem->child->next->next->next->valuestring, "-_Inf_") == 0) { + cfg->srcdata[count].srcdir.w = -INFINITY; + } + } else { + cfg->srcdata[count].srcdir.w = subitem->child->next->next->next->valuedouble; + } + } + + count++; + } + } } subitem = FIND_JSON_OBJ("IQUV", "Optode.Source.IQUV", src); @@ -2540,7 +2670,19 @@ int mcx_loadjson(cJSON* root, Config* cfg) { cfg->srctype = mcx_keylookup((char*)FIND_JSON_KEY("Type", "Optode.Source.Type", src, "pencil", valuestring), srctypeid); subitem = FIND_JSON_OBJ("Param1", "Optode.Source.Param1", src); - if (subitem) { + if (subitem && cJSON_IsArray(subitem)) { + if (cJSON_IsArray(subitem->child)) { + if (cfg->srcdata && cfg->extrasrclen != cJSON_GetArraySize(subitem) - 1) { + MCX_ERROR(-1, "Length of sub-elements of Pos/Dir/Param1/Param2 must match"); + } + + if (cfg->srcdata == NULL) { + cfg->extrasrclen = cJSON_GetArraySize(subitem) - 1; + } + + subitem = subitem->child; + } + cfg->srcparam1.x = subitem->child->valuedouble; if (subitem->child->next) { @@ -2554,11 +2696,51 @@ int mcx_loadjson(cJSON* root, Config* cfg) { } } } + + if (cfg->extrasrclen > 0) { + int count = 0; + + if (cfg->srcdata == NULL) { + cfg->srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), cfg->extrasrclen); + } + + while (subitem->next && count < cfg->extrasrclen) { + + subitem = subitem->next; + cfg->srcdata[count].srcparam1.x = subitem->child->valuedouble; + + if (subitem->child->next) { + cfg->srcdata[count].srcparam1.y = subitem->child->next->valuedouble; + + if (subitem->child->next->next) { + cfg->srcdata[count].srcparam1.z = subitem->child->next->next->valuedouble; + + if (subitem->child->next->next->next) { + cfg->srcdata[count].srcparam1.w = subitem->child->next->next->next->valuedouble; + } + } + } + + count++; + } + } } subitem = FIND_JSON_OBJ("Param2", "Optode.Source.Param2", src); - if (subitem) { + if (subitem && cJSON_IsArray(subitem)) { + if (cJSON_IsArray(subitem->child)) { + if (cfg->srcdata && cfg->extrasrclen != cJSON_GetArraySize(subitem) - 1) { + MCX_ERROR(-1, "Length of sub-elements of Pos/Dir/param2/Param2 must match"); + } + + if (cfg->srcdata == NULL) { + cfg->extrasrclen = cJSON_GetArraySize(subitem) - 1; + } + + subitem = subitem->child; + } + cfg->srcparam2.x = subitem->child->valuedouble; if (subitem->child->next) { @@ -2572,6 +2754,34 @@ int mcx_loadjson(cJSON* root, Config* cfg) { } } } + + if (cfg->extrasrclen > 0) { + int count = 0; + + if (cfg->srcdata == NULL) { + cfg->srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), cfg->extrasrclen); + } + + while (subitem->next && count < cfg->extrasrclen) { + + subitem = subitem->next; + cfg->srcdata[count].srcparam2.x = subitem->child->valuedouble; + + if (subitem->child->next) { + cfg->srcdata[count].srcparam2.y = subitem->child->next->valuedouble; + + if (subitem->child->next->next) { + cfg->srcdata[count].srcparam2.z = subitem->child->next->next->valuedouble; + + if (subitem->child->next->next->next) { + cfg->srcdata[count].srcparam2.w = subitem->child->next->next->next->valuedouble; + } + } + } + + count++; + } + } } cfg->omega = FIND_JSON_KEY("Frequency", "Optode.Source.Frequency", src, 0.f, valuedouble); @@ -2714,8 +2924,8 @@ int mcx_loadjson(cJSON* root, Config* cfg) { } } - if (cfg->medianum + cfg->detnum > MAX_PROP_AND_DETECTORS) { - MCX_ERROR(-4, "input media types plus detector number exceeds the maximum total (4000)"); + if (cfg->medianum + cfg->detnum + (cfg->extrasrclen << 2) > MAX_PROP_AND_DETECTORS) { + MCX_ERROR(-4, "input media types plus detector number plus additional sources exceeds the maximum total (4000)"); } if (Session) { @@ -3447,6 +3657,14 @@ void mcx_validatecfg(Config* cfg, float* detps, int dimdetps[2], int seedbyte) { cfg->srcpos.x--; cfg->srcpos.y--; cfg->srcpos.z--; /*convert to C index, grid center*/ + + if (cfg->extrasrclen) { + for (int i = 0; i < cfg->extrasrclen; i++) { + cfg->srcdata[i].srcpos.x--; + cfg->srcdata[i].srcpos.y--; + cfg->srcdata[i].srcpos.z--; + } + } } if (cfg->tstep == 0.f) { @@ -4924,7 +5142,7 @@ int mcx_lookupindex(char* key, const char* index) { void mcx_version(Config* cfg) { const char ver[] = "$Rev:: $ " MCX_VERSION; - int v = 0; + uint v = 0; sscanf(ver, "$Rev::%x", &v); MCX_FPRINTF(cfg->flog, "MCX Revision %x\n", v); exit(0); diff --git a/src/mcx_utils.h b/src/mcx_utils.h index b613a8d8..8ade1903 100644 --- a/src/mcx_utils.h +++ b/src/mcx_utils.h @@ -94,6 +94,14 @@ typedef struct MCXPolarizeMedium { float nmed; /**< background medium refrative index */ } POLMedium; +typedef struct MCXExtraSource { + float4 srcpos; /**< initial position vector + initial weight */ + float4 srcdir; /**< initial directon vector + focal length */ + float4 srcparam1; /**< source parameters set 1 */ + float4 srcparam2; /**< source parameters set 2 */ +} ExtraSrc; + + /** * Header data structure in .mch/.mct files to store detected photon data * This header has a total of 256 bytes @@ -263,6 +271,8 @@ typedef struct MCXConfig { float* invcdf; /**< equal-space sampled inversion of CDF(cos(theta)) for the phase function of the zenith angle */ unsigned int nangle; /**< number of samples for inverse-cdf of launch angle, will be added by 2 to include -1 and 1 on the two ends */ float* angleinvcdf; /**< equal-space sampled inversion of CDF(cos(theta)) for the phase function of the zenith angle of photon launch */ + unsigned int extrasrclen; /**< length of additional sources */ + ExtraSrc* srcdata; /**< buffer to store multiple source input data */ } Config; #ifdef __cplusplus diff --git a/src/mcxlab.cpp b/src/mcxlab.cpp index 99086742..44dda611 100644 --- a/src/mcxlab.cpp +++ b/src/mcxlab.cpp @@ -564,15 +564,151 @@ 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_VEC34_FIELD(cfg, srcpos) - GET_VEC34_FIELD(cfg, srcdir) GET_VEC3_FIELD(cfg, steps) GET_VEC3_FIELD(cfg, crop0) GET_VEC3_FIELD(cfg, crop1) - GET_VEC4_FIELD(cfg, srcparam1) - GET_VEC4_FIELD(cfg, srcparam2) GET_VEC4_FIELD(cfg, srciquv) - else if (strcmp(name, "vol") == 0) { + else if (strcmp(name, "srcpos") == 0) { + arraydim = mxGetDimensions(item); + + if (arraydim[0] == 0 || arraydim[1] < 3 || arraydim[1] > 4) { + mexErrMsgTxt("the 'srcpos' field must have 3 or 4 columns (x,y,z,w0)"); + } + + double* val = mxGetPr(item); + + for (i = 0; i < arraydim[1]; i++) { + ((float*)(&cfg->srcpos.x))[i] = val[i * arraydim[0]]; + } + + printf("mcx.srcpos=[%g %g %g %g];\n", cfg->srcpos.x, cfg->srcpos.y, cfg->srcpos.z, cfg->srcpos.w); + + if (arraydim[0] == 1 && cfg->extrasrclen == 0) { + return; + } + + if (cfg->extrasrclen && cfg->extrasrclen != arraydim[0] - 1) { + mexErrMsgTxt("Length of sub-elements of srcpos/srcdir/srcparam1/srcparam2 must match"); + } else { + cfg->extrasrclen = arraydim[0] - 1; + } + + if (cfg->srcdata == NULL) { + cfg->srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), cfg->extrasrclen); + } + + for (j = 0; j < arraydim[1]; j++) + for (i = 0; i < cfg->extrasrclen; i++) { + ((float*)(&cfg->srcdata[i].srcpos.x))[j] = val[j * arraydim[0] + i + 1]; + } + + printf("mcx.extrasrclen=%d;\n", cfg->extrasrclen); + } else if (strcmp(name, "srcdir") == 0) { + arraydim = mxGetDimensions(item); + + if (arraydim[0] == 0 || arraydim[1] < 3 || arraydim[1] > 4) { + mexErrMsgTxt("the 'srcdir' field must have 3 or 4 columns (vx,vy,vz,focallength)"); + } + + double* val = mxGetPr(item); + + for (i = 0; i < arraydim[1]; i++) { + ((float*)(&cfg->srcdir.x))[i] = val[i * arraydim[0]]; + } + + printf("mcx.srcdir=[%g %g %g %g];\n", cfg->srcdir.x, cfg->srcdir.y, cfg->srcdir.z, cfg->srcdir.w); + + if (arraydim[0] == 1 && cfg->extrasrclen == 0) { + return; + } + + if (cfg->extrasrclen && cfg->extrasrclen != arraydim[0] - 1) { + mexErrMsgTxt("Length of sub-elements of srcdir/srcdir/srcparam1/srcparam2 must match"); + } else { + cfg->extrasrclen = arraydim[0] - 1; + } + + if (cfg->srcdata == NULL) { + cfg->srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), cfg->extrasrclen); + } + + for (j = 0; j < arraydim[1]; j++) + for (i = 0; i < cfg->extrasrclen; i++) { + ((float*)(&cfg->srcdata[i].srcdir.x))[j] = val[j * arraydim[0] + i + 1]; + } + + printf("mcx.extrasrclen=%d;\n", cfg->extrasrclen); + } else if (strcmp(name, "srcparam1") == 0) { + arraydim = mxGetDimensions(item); + + if (arraydim[0] == 0 || arraydim[1] != 4) { + mexErrMsgTxt("the 'srcparam1' field must have 4 columns"); + } + + double* val = mxGetPr(item); + + for (i = 0; i < arraydim[1]; i++) { + ((float*)(&cfg->srcparam1.x))[i] = val[i * arraydim[0]]; + } + + printf("mcx.srcparam1=[%g %g %g %g];\n", cfg->srcparam1.x, cfg->srcparam1.y, cfg->srcparam1.z, cfg->srcparam1.w); + + if (arraydim[0] == 1 && cfg->extrasrclen == 0) { + return; + } + + if (cfg->extrasrclen && cfg->extrasrclen != arraydim[0] - 1) { + mexErrMsgTxt("Length of sub-elements of srcparam1/srcparam1/srcparam1/srcparam2 must match"); + } else { + cfg->extrasrclen = arraydim[0] - 1; + } + + if (cfg->srcdata == NULL) { + cfg->srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), cfg->extrasrclen); + } + + for (j = 0; j < arraydim[1]; j++) + for (i = 0; i < cfg->extrasrclen; i++) { + ((float*)(&cfg->srcdata[i].srcparam1.x))[j] = val[j * arraydim[0] + i + 1]; + } + + printf("mcx.extrasrclen=%d;\n", cfg->extrasrclen); + } else if (strcmp(name, "srcparam2") == 0) { + arraydim = mxGetDimensions(item); + + if (arraydim[0] == 0 || arraydim[1] != 4) { + mexErrMsgTxt("the 'srcparam2' field must have 4 columns"); + } + + double* val = mxGetPr(item); + + for (i = 0; i < arraydim[1]; i++) { + ((float*)(&cfg->srcparam2.x))[i] = val[i * arraydim[0]]; + } + + printf("mcx.srcparam2=[%g %g %g %g];\n", cfg->srcparam2.x, cfg->srcparam2.y, cfg->srcparam2.z, cfg->srcparam2.w); + + if (arraydim[0] == 1 && cfg->extrasrclen == 0) { + return; + } + + if (cfg->extrasrclen && cfg->extrasrclen != arraydim[0] - 1) { + mexErrMsgTxt("Length of sub-elements of srcparam2/srcparam2/srcparam2/srcparam2 must match"); + } else { + cfg->extrasrclen = arraydim[0] - 1; + } + + if (cfg->srcdata == NULL) { + cfg->srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), cfg->extrasrclen); + } + + for (j = 0; j < arraydim[1]; j++) + for (i = 0; i < cfg->extrasrclen; i++) { + ((float*)(&cfg->srcdata[i].srcparam2.x))[j] = val[j * arraydim[0] + i + 1]; + } + + printf("mcx.extrasrclen=%d;\n", cfg->extrasrclen); + } else if (strcmp(name, "vol") == 0) { dimtype dimxyz; cfg->mediabyte = 0; arraydim = mxGetDimensions(item); diff --git a/src/pmcx.cpp b/src/pmcx.cpp index 0eb73ae0..96174a15 100644 --- a/src/pmcx.cpp +++ b/src/pmcx.cpp @@ -304,7 +304,8 @@ void parseVolume(const py::dict& user_cfg, Config& mcx_config) { float f2h[2]; int offset = (mcx_config.mediabyte == MEDIA_ASGN_F2H); - if(mcx_config.mediabyte == MEDIA_ASGN_F2H) { + + if (mcx_config.mediabyte == MEDIA_ASGN_F2H) { mcx_config.vol = static_cast(realloc(mcx_config.vol, dim_xyz * 2 * sizeof(unsigned int))); } @@ -428,16 +429,200 @@ 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_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); GET_VEC3_FIELD(user_cfg, mcx_config, crop1, uint); - GET_VEC4_FIELD(user_cfg, mcx_config, srcparam1, float); - GET_VEC4_FIELD(user_cfg, mcx_config, srcparam2, float); GET_VEC4_FIELD(user_cfg, mcx_config, srciquv, float); parseVolume(user_cfg, mcx_config); + if (user_cfg.contains("srcpos")) { + auto f_style_volume = py::array_t < float, py::array::f_style | py::array::forcecast >::ensure(user_cfg["srcpos"]); + + if (!f_style_volume) { + throw py::value_error("Invalid srcpos field value"); + } + + auto buffer_info = f_style_volume.request(); + + size_t arraydim[2] = {0}; + + if (buffer_info.shape.size() == 1) { + arraydim[0] = 1; + arraydim[1] = buffer_info.shape.at(0); + } else if (buffer_info.shape.size() > 1) { + arraydim[0] = buffer_info.shape.at(0); + arraydim[1] = buffer_info.shape.at(1); + } + + if (arraydim[0] == 0 || arraydim[1] < 3 || arraydim[1] > 4) { + throw py::value_error("the 'srcpos' field must have 3 or 4 columns (x,y,z,w0)"); + } + + auto val = static_cast(buffer_info.ptr); + + for (int i = 0; i < arraydim[1]; i++) { + ((float*)(&mcx_config.srcpos.x))[i] = val[i * arraydim[0]]; + } + + if (arraydim[0] > 1 || mcx_config.extrasrclen > 0) { + if (mcx_config.extrasrclen && mcx_config.extrasrclen != arraydim[0] - 1) { + throw py::value_error("Length of sub-elements of srcpos/srcdir/srcparam1/srcparam2 must match"); + } else { + mcx_config.extrasrclen = arraydim[0] - 1; + } + + if (mcx_config.srcdata == NULL) { + mcx_config.srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), mcx_config.extrasrclen); + } + + for (int j = 0; j < arraydim[1]; j++) + for (int i = 0; i < mcx_config.extrasrclen; i++) { + ((float*)(&mcx_config.srcdata[i].srcpos.x))[j] = val[j * arraydim[0] + i + 1]; + } + } + } + + if (user_cfg.contains("srcdir")) { + auto f_style_volume = py::array_t < float, py::array::f_style | py::array::forcecast >::ensure(user_cfg["srcdir"]); + + if (!f_style_volume) { + throw py::value_error("Invalid srcdir field value"); + } + + auto buffer_info = f_style_volume.request(); + + size_t arraydim[2] = {0}; + + if (buffer_info.shape.size() == 1) { + arraydim[0] = 1; + arraydim[1] = buffer_info.shape.at(0); + } else if (buffer_info.shape.size() > 1) { + arraydim[0] = buffer_info.shape.at(0); + arraydim[1] = buffer_info.shape.at(1); + } + + if (arraydim[0] == 0 || arraydim[1] < 3 || arraydim[1] > 4) { + throw py::value_error("the 'srcdir' field must have 3 or 4 columns (vx,vy,vz,focallength)"); + } + + auto val = static_cast(buffer_info.ptr); + + for (int i = 0; i < arraydim[1]; i++) { + ((float*)(&mcx_config.srcdir.x))[i] = val[i * arraydim[0]]; + } + + if (arraydim[0] > 1 || mcx_config.extrasrclen > 0) { + if (mcx_config.extrasrclen && mcx_config.extrasrclen != arraydim[0] - 1) { + throw py::value_error("Length of sub-elements of srcpos/srcdir/srcparam1/srcparam2 must match"); + } else { + mcx_config.extrasrclen = arraydim[0] - 1; + } + + if (mcx_config.srcdata == NULL) { + mcx_config.srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), mcx_config.extrasrclen); + } + + for (int j = 0; j < arraydim[1]; j++) + for (int i = 0; i < mcx_config.extrasrclen; i++) { + ((float*)(&mcx_config.srcdata[i].srcdir.x))[j] = val[j * arraydim[0] + i + 1]; + } + } + } + + if (user_cfg.contains("srcparam1")) { + auto f_style_volume = py::array_t < float, py::array::f_style | py::array::forcecast >::ensure(user_cfg["srcparam1"]); + + if (!f_style_volume) { + throw py::value_error("Invalid srcparam1 field value"); + } + + auto buffer_info = f_style_volume.request(); + + size_t arraydim[2] = {0}; + + if (buffer_info.shape.size() == 1) { + arraydim[0] = 1; + arraydim[1] = buffer_info.shape.at(0); + } else if (buffer_info.shape.size() > 1) { + arraydim[0] = buffer_info.shape.at(0); + arraydim[1] = buffer_info.shape.at(1); + } + + if (arraydim[0] == 0 || arraydim[1] != 4) { + throw py::value_error("the 'srcparam1' field must have 4 columns"); + } + + auto val = static_cast(buffer_info.ptr); + + for (int i = 0; i < arraydim[1]; i++) { + ((float*)(&mcx_config.srcparam1.x))[i] = val[i * arraydim[0]]; + } + + if (arraydim[0] > 1 || mcx_config.extrasrclen > 0) { + if (mcx_config.extrasrclen && mcx_config.extrasrclen != arraydim[0] - 1) { + throw py::value_error("Length of sub-elements of srcpos/srcdir/srcparam1/srcparam2 must match"); + } else { + mcx_config.extrasrclen = arraydim[0] - 1; + } + + if (mcx_config.srcdata == NULL) { + mcx_config.srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), mcx_config.extrasrclen); + } + + for (int j = 0; j < arraydim[1]; j++) + for (int i = 0; i < mcx_config.extrasrclen; i++) { + ((float*)(&mcx_config.srcdata[i].srcparam1.x))[j] = val[j * arraydim[0] + i + 1]; + } + } + } + + if (user_cfg.contains("srcparam2")) { + auto f_style_volume = py::array_t < float, py::array::f_style | py::array::forcecast >::ensure(user_cfg["srcparam2"]); + + if (!f_style_volume) { + throw py::value_error("Invalid srcparam2 field value"); + } + + auto buffer_info = f_style_volume.request(); + + size_t arraydim[2] = {0}; + + if (buffer_info.shape.size() == 1) { + arraydim[0] = 1; + arraydim[1] = buffer_info.shape.at(0); + } else if (buffer_info.shape.size() > 1) { + arraydim[0] = buffer_info.shape.at(0); + arraydim[1] = buffer_info.shape.at(1); + } + + if (arraydim[0] == 0 || arraydim[1] != 4) { + throw py::value_error("the 'srcparam2' field must have 4 columns"); + } + + auto val = static_cast(buffer_info.ptr); + + for (int i = 0; i < arraydim[1]; i++) { + ((float*)(&mcx_config.srcparam2.x))[i] = val[i * arraydim[0]]; + } + + if (arraydim[0] > 1 || mcx_config.extrasrclen > 0) { + if (mcx_config.extrasrclen && mcx_config.extrasrclen != arraydim[0] - 1) { + throw py::value_error("Length of sub-elements of srcpos/srcdir/srcparam1/srcparam2 must match"); + } else { + mcx_config.extrasrclen = arraydim[0] - 1; + } + + if (mcx_config.srcdata == NULL) { + mcx_config.srcdata = (ExtraSrc*)calloc(sizeof(ExtraSrc), mcx_config.extrasrclen); + } + + for (int j = 0; j < arraydim[1]; j++) + for (int i = 0; i < mcx_config.extrasrclen; i++) { + ((float*)(&mcx_config.srcdata[i].srcparam2.x))[j] = val[j * arraydim[0] + i + 1]; + } + } + } + if (user_cfg.contains("detpos")) { auto f_style_volume = py::array_t < float, py::array::f_style | py::array::forcecast >::ensure(user_cfg["detpos"]);