SVMC also available within mcx environment?

[Bastiaan]

Dear Qianqian,

I'm using python to set-up (building json input files etc) and execute all my simulations. I'm very interested in using the SVMC algorithm. However, correct me if I'm wrong, it is currently only available within the MCXLAB environment. 

Will you also make the algorithm available within the mcx environment?

Cheers,

Bastiaan

[Qianqian Fang]

hi Bastiaan,

in a way, yes, it was already supported, but you will have to do some extra work to "repack" the 8-byte voxel records generated by mcxsvmc.m into 2x 4-byte segments, which is the voxel format used internally by mcx's GPU kernel.

I just committed two patches to avoid this manual step, so you can directly use the volume created by mcxsvmc.m (cast to uint8)

https://github.com/fangq/mcx/commit/de9850c1bf32750f03de061ef4a86a3e0f64a80f
https://github.com/fangq/mcx/commit/4010d99f85084fdfc834f1c61e592d1065953d26

here are the two needed steps:

1. create the SVMC formatted volume (8-byte per voxel, see Fig. 2g of the svmc paper) using mcxsvmc.m, cast the volume to uint8(), and write it to a .bin file (use fopen/fwrite/fclose in matlab), then set the .bin file as the volume input file in the JSON input, under Domain.VolumeFile

2. when running mcx, you need to tell mcx you are using the svmc 8-byte volume format by attaching -K svmc or  --mediabyte svmc, or replace svmc by 97, where identifier 97 is defined here

https://github.com/fangq/mcx/blob/master/src/mcx_const.h#L62

I just recompiled mcx nightly builds after these two patches, let me know if this works for you.

Qianqian

Cheers,

Bastiaan

--
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/bafdfd32-f1db-4681-9ed0-cac21f4eeb6fn%40googlegroups.com.

[Qianqian Fang]

hi Bastiaan,

I made additional commits and testing and now svmc is supported in the command line.

I also added an example

https://github.com/fangq/mcx/tree/master/example/svmcsphshells

let me know if you are able to run this example.

Qianqian

[Bastiaan]

Hi Qianqian,

Thanks a lot! I need some time to implement this on my machine but I will keep you posted. 

Cheers,

Bastiaan

[Bastiaan]

I'm able to run this example, thanks!

[Bastiaan]

Hi Qianqian, 

I'm heaving some problems with the nightly build version that I've downloaded to use the split-voxel algorithm with json input files. 

1) No photons are detected, although a detector is included that detects photons when running a the voxel based algorithm

2) when using a pencil beam, no photons are detected. 

Below the matlab scripts that I use to generate the input files and run the simulation. I think I followed your instructions above correctly however the program is not behaving as I am expecting. If you have time, could you comment on this? 

Thanks in advance!

Cheers,

Bastiaan

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Matlab program to scan a point source over a surface that includes

% a bloodvessel. Two types of simulations will be done:

% a)Voxel based

% b)Split-voxel based

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

clear all;

addpath ('C:\Sources\mcx-win-x86_64-nightlybuild_SVMC_25may\mcx\utils')

addpath('...iso2mesh-2018-allinone\iso2mesh')

mcx_source = 'C:\Sources\mcx-win-x86_64-v2020\mcx\bin\mcx.exe'

svmcx_source = 'C:\Sources\mcx-win-x86_64-nightlybuild_SVMC_25may\mcx\bin\mcx.exe'

FileVMC = '....\Simulations\VMC'

FileSVMC = '....\Simulations\SVMC'

shapefileVMC = strcat(FileVMC,'\shape.bin')

shapefileSVMC = strcat(FileSVMC,'\shape.bin')

%% common MC setup

cfg.nphoton=1e8;

cfg.seed=randi([1 2^31-1], 1, 1); %random seed

% pencil beam light source

%cfg.srctype = 'pencil'

cfg.srctype = 'disk'

cfg.srcparam1=[1 0 0]

cfg.srcdir=[0 0 1];

cfg.issrcfrom0=1;

% optical properties (tissue-like multi-layered media)

cfg.prop=[0.0 0.0 1.0 1.0;   % background (air, void)

          1 1.0 .9 1.37;  % tag1

          1 1.0 .9 1.5;]  % tag2

%          1.0 1.0 .9 1.37;  % tag3

%          1.0 1.0 .9 1.37]; % tag4

      

% time-domain simulation parameters

cfg.tstart=0;

cfg.tend=5e-9;

cfg.tstep=5e-9;

% enable boundary reflection/refraction

cfg.isreflect=1;

% spatial resolution

cfg.unitinmm=1;

% output fluence

cfg.outputtype='fluence';

% GPU settings

cfg.gpuid=1;

cfg.autopilot=1;

cfg.detpos=[30 30 1 30];

cfg_mcx = cfg

%% prepare mcx input volume

dim=60;

[xi,yi,zi]=ndgrid(1:dim,1:dim,1:dim);

dist=(xi-30.5).^2+(zi-8.5).^2;

mcxvol=zeros(size(xi));

mcxvol(2:59,2:59,2:59)=1;

mcxvol(dist<25)=2;

figure(1);

%subplot(131);

imagesc(squeeze(mcxvol(:,30,:))')

axis equal

fileID = fopen(shapefileVMC, 'w');

fwrite(fileID, uint8(mcxvol));

fclose(fileID);

cfg_mcx.VolumeFile = shapefileVMC;

cfg_mcx.MediaFormat = 'byte';

cfg_mcx.dim = [dim,dim,dim];%to dimx, dimy, dimz

svmcvol=mcxsvmc(mcxvol,'smoothing',1 ); % 1: enable gaussian smoothing 0: otherwise

fprintf('SVMC preprocessing complete, ');

%display

[bmask,uniqidx,nc,nn,totalarea,iso] =mcxsvmc(mcxvol,'smoothing',1 );

figure(3)

plotmesh(iso.vertices,iso.faces,'facealpha',0.4,'facecolor','b','edgealpha',0.2)

%display

fileID = fopen(shapefileSVMC, 'w');

fwrite(fileID, uint8(svmcvol));

fclose(fileID);

cfg_svmc=cfg;

cfg_svmc.VolumeFile = shapefileSVMC;

cfg_svmc.MediaFormat = 'svmc';

cfg_svmc.dim = [dim,dim,dim];%to dimx, dimy, dimz

%loop over srcpos, save json and .bat file and run .bat file

for xpos = 30:31;

    cfg_svmc.srcpos=[xpos 30.5 1];

    cfg_mcx.srcpos=[xpos 30.5 1];

    

    jsonfileSVMC = strcat(FileSVMC,'\xpos_',strrep(sprintf('%.1f',xpos),'.','p'));

    jsonfileVMC = strcat(FileVMC,'\xpos_',strrep(sprintf('%.1f',xpos),'.','p'));

    mcx2json_bf(cfg_svmc, jsonfileSVMC);

    mcx2json_bf(cfg_mcx, jsonfileVMC);

    

    sv_batstring = strcat(svmcx_source, ' -f',32, strcat(jsonfileSVMC,'.json'),' -n 1e6 -w DSPXV -K svmc -F bnii -D P');

    v_batstring  = strcat(mcx_source, ' -f',32, strcat(jsonfileVMC,'.json'),' -n 1e6 -w DSPXV -F bnii -D P');

    system(sv_batstring);

    system(v_batstring);

end

and mcx2json_bf:

function mcx2json(cfg,filestub)

%

% Format:

%    mcx2json(cfg,filestub)

%

% Save MCXLAB simulation configuration to a JSON file for MCX binary

%

% Author: Qianqian Fang <q.fang at neu.edu>

%

% Input:

%    cfg: a struct defining the parameters associated with a simulation. 

%         Please run 'help mcxlab' or 'help mmclab' to see the details.

%         mcxpreview supports the cfg input for both mcxlab and mmclab.

%    filestub: the filestub is the name stub for all output files,including

%         filestub.json: the JSON input file

%         filestub_vol.bin: the volume file if cfg.vol is defined

%         filestub_shapes.json: the domain shape file if cfg.shapes is defined

%         filestub_pattern.bin: the domain shape file if cfg.pattern is defined

%

% Dependency:

%    this function depends on the savejson/saveubjson functions from the 

%    Iso2Mesh toolbox (http://iso2mesh.sf.net) or JSONlab toolbox 

%    (http://iso2mesh.sf.net/jsonlab)

%

% This function is part of Monte Carlo eXtreme (MCX) URL: http://mcx.space

%

% License: GNU General Public License version 3, please read LICENSE.txt for details

%

[fpath, fname, fext]=fileparts(filestub);

filestub=fullfile(fpath,fname);

%% define the optodes: sources and detectors

Optode.Source=struct();

Optode.Source=copycfg(cfg,'srcpos',Optode.Source,'Pos');

Optode.Source=copycfg(cfg,'srcdir',Optode.Source,'Dir');

Optode.Source=copycfg(cfg,'srcparam1',Optode.Source,'Param1');

Optode.Source=copycfg(cfg,'srcparam2',Optode.Source,'Param2');

Optode.Source=copycfg(cfg,'srctype',Optode.Source,'Type');

Optode.Source=copycfg(cfg,'srcnum',Optode.Source,'SrcNum');

if(isfield(cfg,'detpos') && ~isempty(cfg.detpos))

    Optode.Detector=struct();

    Optode.Detector=cell2struct(mat2cell(cfg.detpos, ones(1,size(cfg.detpos,1)),[3 1]), {'Pos','R'} ,2);

    if(length(Optode.Detector)==1)

        Optode.Detector={Optode.Detector};

    end

end

if(isfield(cfg,'srcpattern') && ~isempty(cfg.srcpattern))

    Optode.Source.Pattern.Nx=size(cfg.srcpattern,1);

    Optode.Source.Pattern.Ny=size(cfg.srcpattern,2);

    Optode.Source.Pattern.Nz=size(cfg.srcpattern,3);

    Optode.Source.Pattern.Data=[filestub '_pattern.bin'];

    fid=fopen(Optode.Source.Pattern.Data,'wb');

    fwrite(fid,cfg.srcpattern,'float32');

    fclose(fid);

end

%% define the domain and optical properties

Domain=struct();

Domain=copycfg(cfg,'issrcfrom0',Domain,'OriginType',0);

Domain=copycfg(cfg,'unitinmm',Domain,'LengthUnit');

Domain=copycfg(cfg,'MediaFormat',Domain,'MediaFormat');

Domain.Media=cell2struct(num2cell(cfg.prop), {'mua','mus','g','n'} ,2)';

if(isfield(cfg,'shapes') && ischar(cfg.shapes))

    Shapes=loadjson(cfg.shapes);

    Shapes=Shapes.Shapes;

end

    

if(isfield(cfg,'vol') && ~isempty(cfg.vol) && ~isfield(Domain,'VolumeFile'))

    

    switch(class(cfg.vol))

        case {'uint8','int8'}

            Domain.MediaFormat='byte';

            if(ndims(cfg.vol)==4 && size(cfg.vol,1)==4)

                Domain.MediaFormat='asgn_byte';

            elseif(ndims(cfg.vol)==4 && size(cfg.vol,1)==8)

              

                dim = size(cfg.vol);

                %cfg.vol=reshape(typecast(uint8(cfg.vol(:)),'uint64'),size(cfg.vol,2,3,4));%matlab 2019b notation

                cfg.vol=reshape(typecast(uint8(cfg.vol(:)),'uint64'),dim(2:end));

                Domain.MediaFormat='svmc';

            end

        case {'uint16','int16'}

            Domain.MediaFormat='short';

            if(ndims(cfg.vol)==4 && size(cfg.vol,1)==2)

                Domain.MediaFormat='muamus_short';

            end

        case {'uint32','int32'}

            Domain.MediaFormat='integer';

        case {'single','double'}

            if(isa(cfg.vol,'double'))

                cfg.vol=single(cfg.vol);

            end

            if(all(mod(cfg.vol(:),1) == 0))

                Domain.MediaFormat='integer';

            elseif(ndims(cfg.vol)==4)

                if(size(cfg.vol,1))==1

                    Domain.MediaFormat='mua_float';

                elseif(size(cfg.vol,1)==2)

                    Domain.MediaFormat='muamus_float';

                end

            end

        otherwise

            error('cfg.vol has format that is not supported');

    end

    Domain.Dim=size(cfg.vol);

    if(length(Domain.Dim)==4)

        Domain.Dim(1)=[];

    end

    if(exist('Shapes','var'))

        Domain.VolumeFile=[filestub '_vol.bin'];

        fid=fopen(Domain.VolumeFile,'wb');

        fwrite(fid,cfg.vol,class(cfg.vol));

        fclose(fid);

    else

        Domain.VolumeFile='';

        Shapes=cfg.vol;

    end

end

if(isfield(cfg,'VolumeFile')) 

    Domain.VolumeFile=cfg.VolumeFile;

    %Domain.MediaFormat='svmc';

    Domain.Dim=cfg.dim;

end

%% define the simulation session flags

Session.ID=filestub;

Session=copycfg(cfg,'isreflect',Session,'DoMismatch');

Session=copycfg(cfg,'issave2pt',Session,'DoSaveVolume');

Session=copycfg(cfg,'issavedet',Session,'DoPartialPath');

Session=copycfg(cfg,'issaveexit',Session,'DoSaveExit');

Session=copycfg(cfg,'issaveseed',Session,'DoSaveSeed');

Session=copycfg(cfg,'isnormalize',Session,'DoNormalize');

Session=copycfg(cfg,'outputformat',Session,'OutputFormat');

Session=copycfg(cfg,'outputtype',Session,'OutputType');

Session=copycfg(cfg,'debuglevel',Session,'Debug');

Session=copycfg(cfg,'autopilot',Session,'DoAutoThread');

if(isfield(cfg,'seed') && numel(cfg.seed)==1)

    Session.RNGSeed=cfg.seed;

end

Session=copycfg(cfg,'nphoton',Session,'Photons');

Session=copycfg(cfg,'rootpath',Session,'RootPath');

%% define the forward simulation settings

Forward.T0=cfg.tstart;

Forward.T1=cfg.tend;

Forward.Dt=cfg.tstep;

%% assemble the complete input, save to a JSON or UBJSON input file

mcxsession=struct('Session', Session, 'Forward', Forward, 'Optode',Optode, 'Domain', Domain);

if(exist('Shapes','var'))

    mcxsession.Shapes=Shapes;

end

if(strcmp(fext,'ubj'))

    saveubjson('',mcxsession,[filestub,'.ubj']);

else

    savejson('',mcxsession,'filename',[filestub,'.json'],'compression','zlib');

end

disp(Domain)

function outdata=copycfg(cfg,name,outroot,outfield,defaultval)

if(nargin>=5)

    outroot.(outfield)=defaultval;

end

if(isfield(cfg,name))

    outroot.(outfield)=cfg.(name);

end

outdata=outroot;

[Fang, Qianqian]

On 6/14/21 11:08 AM, Bastiaan wrote:

Hi Qianqian, 

I'm heaving some problems with the nightly build version that I've downloaded to use the split-voxel algorithm with json input files. 

1) No photons are detected, although a detector is included that detects photons when running a the voxel based algorithm

2) when using a pencil beam, no photons are detected. 

Below the matlab scripts that I use to generate the input files and run the simulation. I think I followed your instructions above correctly however the program is not behaving as I am expecting. If you have time, could you comment on this?

hi Bastiaan,

it looks like you used matlab to create the .json file and run mcx from the command line.

does this happen when you define the problem using mcxlab cfg data structure (instead of saving as json)? I am CCing my student Shijie, who implemented the svmc feature. From what I can tell, the detector processing in an SVMC volume was handled properly in our code.

if mcxlab does work, but .json does not, please send me your .json file so I can test (can not run your script because of local folders).

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/011e8982-050a-4029-9a98-f031e9cfb8dfn%40googlegroups.com.

[Bastiaan]

Hi Qianqian, 

With mcxlab I can run the examples and I can also my own simulations. With .json as input file in mcx I can but there are no photons detected with the svmc code and input file. 

I have placed the input files for the svmc mcx simulation here:

https://drive.google.com/drive/folders/19srt6LhEMUszH5AGKDzUBS3OKdCQL28D?usp=sharing

Thanks in advance for all your help!

Cheers,

Bastiaan

[Fang, Qianqian]

On 6/22/21 7:11 AM, Bastiaan wrote:

Hi Qianqian, 

With mcxlab I can run the examples and I can also my own simulations.

just to clarify - your simulation did return detected photons when you run it in mcxlab?

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/5b7dbd21-5f6c-49f3-8370-21728d43100dn%40googlegroups.com.

[Bastiaan Florijn]

Yes, that's right. 

Op di 22 jun. 2021 19:53 schreef Fang, Qianqian <q.f...@northeastern.edu>:

[Bastiaan]

Hi Qianqian, 

With the nighltybuild mcxlab version downloaded on the 23rd of June I have the same problem, no photons are detected. I will place the matlab script that I use in the folder I recently shared. 

Also, the for loop in my matlab program (where I call the mcxlab function) is continuously stopped because of  a error caused by the GPU:

'MCXLAB ERROR -46 in unit mcx_core.cu:2398: all CUDA-capable devices are busy or unavailable

Error from thread (0): all CUDA-capable devices are busy or unavailable

C++ Error: MCXLAB Terminated due to an exception!'

I don't have this problem when running simulations with mcx using .json input files. 

Cheers,

Bastiaan

[Bastiaan]

Hi Qianqian,

Do you also see similar problems with the split voxel nightly build version?

Cheers,

Bastiaan

[Fang, Qianqian]

On 7/1/21 3:41 AM, Bastiaan wrote:

Hi Qianqian,

Do you also see similar problems with the split voxel nightly build version?

hi Bastiaan,

I spent a bit time yesterday looking into this issue, and I was able to reproduce both problems - a memory error, and then a cuda device busy message after an CUDA error.

the memory error appears in both mcxlab and command line, despite sometimes mcxlab did not crash. by adding -d 0, I traced the the issue to invalid lower/upper labels in the svmc data (max label # is 2, but the labels read is over 120).

I plan to do more debugging today to hopefully come up with a fix. stay tuned.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/6af1e763-0800-473b-9bd4-29553e9b3962n%40googlegroups.com.

[Bastiaan]

Thanks!

[Bastiaan]

Hi Qianqian, 

Any updates on the split voxel algorithm? 

Cheers,

Bstiaan

[Qianqian Fang]

On 7/13/21 11:07 AM, Bastiaan wrote:

Hi Qianqian, 

Any updates on the split voxel algorithm?

hi Bstiaan,

I am CCing my student Shijie - we spent a couple of hours on zoom discussing this last weekend and had some confusing findings -

basically, running the script testSVMCbloodvesselMCXLAB.m in matlab as well as using the mcx2json_bf.m exported .json file in the command line results in no memory errors. However, running another script testSVMCbloodvesselMCXLAB2.m on one of my machines consistently raised a memory error in my user account, but not in Shijie's account - we are not entirely sure why this was the case.

on the other side, we do see that the svmc method does not detect photons, shijie is currently working on debugging this issue.

I am attaching my related files.

can you let me know what are the results running these two testSVMCbloodvessel scripts on your machine? do they have memory errors or not?

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/1a362a9d-0966-4013-b70c-bb7babc72d43n%40googlegroups.com.

[Bastiaan]

Hi Qianqian,

Sorry for my late response. 

I can run both files now without any memory errors. However, when setting the for loop to run from xpos = 28:40 I get a memory error in  tesSVMCbloodvesselMCXLAB2.m. Setting back the for loop to run from xpos = 28:28, also gives me a memory error again. 

Cheers,

Bastiaan

[Shijie Yan]

Hi Bstiaan,

Can you try setting the medium in the plane y = 1 and y = 60 to void/air (mcxvol(:,1,:)=0; mcxvol(:,60,:)=0;) and then rerun your simulation? Basically, doing so assures that your domain is enclosed by a layer of zero-labeled voxels.

Shijie

---

From: mcx-...@googlegroups.com <mcx-...@googlegroups.com> on behalf of Bastiaan <hcbfl...@gmail.com>
Sent: Tuesday, July 20, 2021 3:26 AM
To: mcx-users <mcx-...@googlegroups.com>
Subject: Re: [mcx-users] SVMC also available within mcx environment?

 

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/c9362616-4519-444e-a706-1d028460b49fn%40googlegroups.com.

[Shijie Yan]

Hi Bastiaan,

The photon detection issue under SVMC mode has been fixed (https://github.com/fangq/mcx/pull/120). Could you please download the latest nightly build MCXLAB and try running the attached MATLAB script? I also attached the final simulation results generated from my end.

Please let us know if you have any further questions.

Shijie

---

From: mcx-...@googlegroups.com <mcx-...@googlegroups.com> on behalf of Shijie Yan <yan....@northeastern.edu>
Sent: Friday, July 30, 2021 6:43 PM

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/PH0PR06MB71412766AC27EDDEEC505FAD9CEC9%40PH0PR06MB7141.namprd06.prod.outlook.com.

[Fang, Qianqian]

hi Bstiaan,

just want to let yo know Shijie had submitted a fix to the photon detection bug under svmc. please see this PR:

https://github.com/fangq/mcx/pull/120

let me know if this patch works for you. nightly-built should have all updated.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/PH0PR06MB71412766AC27EDDEEC505FAD9CEC9%40PH0PR06MB7141.namprd06.prod.outlook.com.

[Bastiaan]

Hi Qianqian and Shijie, 

Thanks for your updates! I will download the patches and test them on my pc. 

Cheers,

Bastiaan

[Bastiaan]

Hi Qianqian and Shijie, 

I have tested the split voxel fix but I still have problems with the detected photons. 

When I run the nightly build, that I have downloaded on the 2nd of September, I can a different amount of detected photons when using the split voxel algorithm (SVMC) or the voxelated algorithm (VMC). This difference is very large. For example, when simulating a volume consisting out of a single tissue with ua = .1, us = 1.0, g = 0.9 and n = 1.37 there are 11 photons detected when using the SVMC algorithm while there are 9250 photons detected when using the VMC algorithm.

Below the matlab code that I use to compare both algorithms. Maybe I made a mistake but as for now, I can't spot one.  Can you figure out what is going on? 

Cheers.

Bastiaan

clear all;

addpath ('C:\Sources\mcx-win-x86_64-nightlybuild_2sept2021\mcx\utils')

addpath('\...\Software\iso2mesh-2018-allinone\iso2mesh')

%mcx_source = 'C:\Sources\mcx-win-x86_64-v2020\mcx\bin\mcx.exe'

mcx_source ='C:\Sources\mcx-win-x86_64-nightlybuild_2sept2021\mcx\bin\mcx.exe'

svmcx_source = 'C:\Sources\mcx-win-x86_64-nightlybuild_2sept2021\mcx\bin\mcx.exe'

dd=datestr(now,'yyyymmdd');

tt = datestr(now,'HHMMSS');

FileVMC = '...\Simulations\VMC'

FileSVMC =   '...  \Simulations\SVMC'

FileVMC = strcat(FileVMC, '\',dd,'_',tt);

FileSVMC = strcat(FileSVMC, '\',dd,'_',tt);

mkdir(FileSVMC);

mkdir(FileVMC);

shapefileVMC = strcat(FileVMC,'\shape.bin')

shapefileSVMC = strcat(FileSVMC,'\shape.bin')

%% common MC setup

cfg.nphoton=1e8;

cfg.seed=randi([1 2^31-1], 1, 1); %random seed

% pencil beam light source

%cfg.srctype = 'pencil'

cfg.srctype = 'disk'

cfg.srcparam1=[2 0 0]

cfg.srcdir=[0 0 1];

cfg.issrcfrom0=1;

% optical properties (tissue-like multi-layered media)

cfg.prop=[0.0 0.0 1.0 1.37;   % background (air, void)

          .1 1.0 .9 1.37;]  % tag1

%          3 .5 .9 1.4;]; % tag2

%          1.0 1.0 .9 1.37;  % tag3

%          1.0 1.0 .9 1.37]; % tag4

      

% time-domain simulation parameters

cfg.tstart=0;

cfg.tend=5e-9;

cfg.tstep=5e-9;

% enable boundary reflection/refraction

cfg.isreflect=1;

% spatial resolution

cfg.unitinmm=1;

% output fluence

cfg.outputtype='energy';

% GPU settings

cfg.gpuid=1;

cfg.autopilot=1;

cfg.detpos=[30 30 1 30];

cfg_mcx = cfg

%% prepare mcx input volume

dim=60;

[xi,yi,zi]=ndgrid(1:dim,1:dim,1:dim);

dist=(xi-30.5).^2+(zi-8).^2;

mcxvol=zeros(size(xi));

mcxvol(2:59,2:59,2:59)=1;

%mcxvol(25:29,25:29,25:29)=2;

%mcxvol(dist<36)=2;

figure(1);

%subplot(131);

imagesc(squeeze(mcxvol(:,30,:))')

axis equal

fileID = fopen(shapefileVMC, 'w');

fwrite(fileID, uint8(mcxvol));

fclose(fileID);

cfg_mcx.VolumeFile = shapefileVMC;

cfg_mcx.MediaFormat = 'byte';

cfg_mcx.dim = [dim,dim,dim];%to dimx, dimy, dimz

svmcvol=mcxsvmc(mcxvol,'smoothing',1 ); % 1: enable gaussian smoothing 0: otherwise

fprintf('SVMC preprocessing complete, ');

%display

[bmask,uniqidx,nc,nn,totalarea,iso] =mcxsvmc(mcxvol,'smoothing',1 );

figure(3)

plotmesh(iso.vertices,iso.faces,'facealpha',0.4,'facecolor','b','edgealpha',0.2)

%display

fileID = fopen(shapefileSVMC, 'w');

fwrite(fileID, uint8(svmcvol));

fclose(fileID);

cfg_svmc=cfg;

cfg_svmc.VolumeFile = shapefileSVMC;

cfg_svmc.MediaFormat = 'svmc';

%cfg_svmc.VolumeFile = shapefileSVMC;

%cfg_svmc.MediaFormat = 'byte';

cfg_svmc.dim = [dim,dim,dim];%to dimx, dimy, dimz

%loop over srcpos, save json and .bat file and run .bat file

for xpos = 10:50;%:31;

    cfg_svmc.srcpos=[xpos 30.5 1];

    cfg_svmc.detpos = [xpos 30.5 1 1];

    

    cfg_mcx.srcpos=[xpos 30.5 1];

    cfg_mcx.detpos = [xpos 30.5 1 1];

    

    jsonfileSVMC = strcat(FileSVMC,'\xpos_',strrep(sprintf('%.1f',xpos),'.','p'));

    jsonfileVMC = strcat(FileVMC,'\xpos_',strrep(sprintf('%.1f',xpos),'.','p'));

    mcx2json_bf(cfg_svmc, jsonfileSVMC);%mcx2json_bf

    mcx2json_bf(cfg_mcx, jsonfileVMC);%mcx2json_bf

    

    sv_batstring = strcat(svmcx_source, ' -f',32, strcat(jsonfileSVMC,'.json'),' -n 1e6 -D M  -F bnii -w DSPXV -D P -K svmc');

    v_batstring  = strcat(mcx_source, ' -f',32, strcat(jsonfileVMC,'.json'),' -n 1e6 -w DSPXV -F bnii -D P');

    system(sv_batstring);

    system(v_batstring);

end

[Qianqian Fang]

hi Shijie, can you take a look?

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/8b0be6c0-93d4-4e65-88d7-50d52bc96056n%40googlegroups.com.

[Shijie Yan]

Hi Bastiaan,

We have fixed the photon detection issue for the SVMC mode (https://github.com/fangq/mcx/pull/123). Once it is merged, please download the nightly-build and try your script again. Thank you for reporting this.

Best,

Shijie

Fix photon detection issue in SVMC mode by ShijieYan · Pull Request #123 · fangq/mcx The issue was reported from this thread. github.com

---

From: Fang, Qianqian <q.f...@northeastern.edu>
Sent: Wednesday, September 15, 2021 9:56 AM
To: mcx-...@googlegroups.com <mcx-...@googlegroups.com>; Bastiaan <hcbfl...@gmail.com>
Cc: Shijie Yan <yan....@northeastern.edu>

[Bastiaan]

Thanks a lot!

When will it be merged so I can download the nightly-build?

Cheers,

Bastiaan

[Fang, Qianqian]

On 9/21/21 3:53 AM, Bastiaan wrote:

Thanks a lot!

When will it be merged so I can download the nightly-build?

already done and should be in the nightly built packages.

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/2ec4fa9b-5e69-484a-b983-85b20e463e30n%40googlegroups.com.

[Bastiaan]

Hi Qianqian and Shijie,

Do you know if someone is also working on a python version of the split voxel algorithm? A Python version of mcxsvmc().m? 

Kind regards,

Bastiaan

[Qianqian Fang]

hi Bastiaan,

I haven't tried it myself, but I assume the recent pybind-based pymcx module contributed by Matin should work with svmc, see

https://github.com/fangq/mcx/pull/149

the interface is nearly identical to mcxlab, you use numpy.ndarray or array.array to define all 3D arrays. for svmc, I assume it also accept 8-byte-per-voxel volumes

https://github.com/fangq/mcx/blob/master/src/pymcx.cpp#L123-L133

I haven't added documentation to this "official" python interface, but you can see the commands for compilation and running basic simulations in this tweet

CCing Matin in case you have feedback to this module.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/0432eed2-76d9-4a3f-8f54-d808ad15023dn%40googlegroups.com.

[Qianqian Fang]

sorry, forgot to attach the link to the tweet with screenshots

https://twitter.com/FangQ/status/1525903673095618560

https://pbs.twimg.com/media/FS0ZSMyXEAEMPdF?format=png&name=4096x4096

[Bastiaan Florijn]

Thanks! I will dive into this!

[Matin Raayai Ardakani]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA512

Bastiaan and Qianqian,
The functionality between mcxlab and pymcx should be the same. I don't
expect it to be bug-free because we haven't extensively tested every
feature in Pymcx, so I recommend running the same input arguments with
both mcxlab and pymcx first. let us know if you run into any issues so
that we can fix it.
Best,
Matin

> > > > > > > Fix photon detection issue in SVMC mode by ShijieYan ·
> > > > > > > Pull Request #123 · fangq/mcx
> > > > > > > The issue was reported from this thread.
> > > > > > > github.com
> > > > > > >

> > > > > > > > > > > > > > > > > > > > > > v2020\mcx\bin\mcx.exe'

> > > > > > > > > > > > > > > > > > > > > > svmcx_source =
> > > > > > > > > > > > > > > > > > > > > > 'C:\Sources\mcx-win-x86_64-

> > > > > > > > > > > > > > > > > > > > > > rrep(sprintf('%.1f',xpos),'
> > > > > > > > > > > > > > > > > > > > > > .','p'));
> > > > > > > > > > > > > > > > > > > > > > jsonfileVMC =
> > > > > > > > > > > > > > > > > > > > > > strcat(FileVMC,'\xpos_',str

> > > > > > > > > > > > > > > > > > > > > > https://groups.google.com/d/msgid/mcx-users/011e8982-050a-4029-9a98-f031e9cfb8dfn%40googlegroups.com
> > > > > > > > > > > > > > > > > > > > > > .

> > > > > > > > > > > > > > > > > > > > >
> > > > > > > > > > > > > > > > > > > > >
> > > > > > > > > > > > > > > > > > > >
> > > > > > > > > > > > > > > > > > > > --
> > > > > > > > > > > > > > > > > > > > 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/5b7dbd21-5f6c-49f3-8370-21728d43100dn%40googlegroups.com
> > > > > > > > > > > > > > > > > > > > .

> > > > > > > > > > > > > > > > > > >
> > > > > > > > > > > > > > > > > > >
> > > > > > > > > > > > > > > >
> > > > > > > > > > > > > > > > --
> > > > > > > > > > > > > > > > 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/6af1e763-0800-473b-9bd4-29553e9b3962n%40googlegroups.com
> > > > > > > > > > > > > > > > .

> > > > > > > > > > > > > > >
> > > > > > > > > > > > > > >
> > > > > > > > > > > > >
> > > > > > > > > > > > > --
> > > > > > > > > > > > > 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/1a362a9d-0966-4013-b70c-bb7babc72d43n%40googlegroups.com
> > > > > > > > > > > > > .

> > > > > > > > > > >
> > > > > > > > > > > --
> > > > > > > > > > > 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/c9362616-4519-444e-a706-1d028460b49fn%40googlegroups.com
> > > > > > > > > > > .

> > > > > > > > > > > --
> > > > > > > > > > > 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/PH0PR06MB71412766AC27EDDEEC505FAD9CEC9%40PH0PR06MB7141.namprd06.prod.outlook.com
> > > > > > > > > > > .

> > > > > > > > > >
> > > > > > > > > >
> > > > > > > >
> > > > > > > > --
> > > > > > > > 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/8b0be6c0-93d4-4e65-88d7-50d52bc96056n%40googlegroups.com
> > > > > > > > .

> > > > > >
> > > > > > --
> > > > > > 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/2ec4fa9b-5e69-484a-b983-85b20e463e30n%40googlegroups.com
> > > > > > .

> > > > >
> > > > >
> > > > >
> > > >
> > > > --
> > > > 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/0432eed2-76d9-4a3f-8f54-d808ad15023dn%40googlegroups.com
> > > > .
-----BEGIN PGP SIGNATURE-----

iQIzBAEBCgAdFiEEvVc+fJGkHWMaFP7CCWbW9yPrH6IFAmLHMagACgkQCWbW9yPr
H6KprQ/+OOjOzb86hPM5BVJb9jKAVD72ZyjDEEEMle9AWxct/NgyyFf3ipYqaINn
9m6KbRa4wui6wIAVgR5fyhJEPR/a/mSsAVY+7e4Kf98HR295O3kI2pddrJpD7FsC
wc3Q5h4Ch13CVQxLICF+nvuppjoPOqwfytbBTBPRgpU6nKhk9mjSrNP7jB2+Oywr
Un3iWG+QaaO1tkRCxYrkBGcTpT35gMfM3U4l9O62dcK52h9+ZMfN4NpZnAO7N4E/
kc1oQ9PRDu5TahRlijOHQ3pTvrhpmmRc+wTXsREUFD+nugAjzKVvQQBQ2yRIro6c
y01niKFn3cNsBFyG9Xo9Lh5lsG3M9LWB7n1kb0Gd3jXSbmh42IJ28Vuw3Lpjdpt5
23B8bml6ZvIZln+Mn5mtasE0xpPcUBc1e5XMAA5AUvvhigWXtlCTayOKXWVYyBBM
zGnd4F9UoJWpx5IsDlXu8L+9B3fiaQrsO6J9vvBgd8y/XkgRWf0qoz2hLQj2n2/H
SFMTPI4fLYOhcS2RRNoSi+UZ1JNVh7MNRuPcLA7i3r1oivWxYZ6YLOzJGCaoPyLH
q/eYm5BFI8E0ZZDDn3DxFbGD8m9x6SnjqehA4JSaLOYcAwX1CdI0Mjz9I+18YEZS
0Dl2UYzcxXy5QTuK9nzXUpwyVNkadGYHvh3SchhleWISHuED0BA=
=fu2J
-----END PGP SIGNATURE-----

[V Zoutenbier]

Dear Martin, 

Sorry if this post finds you twice- I couldn't find where I originally posted it. 

I'm trying to get MCX running on a windows machine using WLS2, which makes the GPU accessible to an ubuntu command prompt, inside my windows machine (We're trying to use the Split Voxel MCX in a python environment). 

When I run CMAKE, I find two errors: 

-- Could NOT find Python3 (missing: Python3_LIBRARIES Python3_INCLUDE_DIRS Development) (found version "3.9.5")
-- Could NOT find Matlab (missing: Matlab_INCLUDE_DIRS Matlab_MEX_LIBRARY Matlab_MEX_EXTENSION Matlab_ROOT_DIR Matlab_MX_LIBRARY) (found version "NOTFOUND")

The first I believe has to do with conflicting Python library installations.  When I run python3 --version it returns with v 3.8; but CMAKE finds v3.9.  I can troubleshoot this further.  

The second error is a bit more troubling; what should be installed from Matlab to get Python to work?  We are trying to minimize (remove) our MATLAB licenses, and would prefer not to have to install the full version.  I tried installing the MATLAB RTE only on my ubuntu PC, but it isn't finding the library/dirs it expects.  Is a full version of MATLAB required to run PyMCX? 

Thanks in advance for any guidence/help in getting PyMCX running. 

Best, 

Vincent

Op zaterdag 9 juli 2022 om 03:20:25 UTC+2 schreef Matin Raayai Ardakani:

[V Zoutenbier]

Hi Quianquian and Matin, 

This question lies somewhere between SVMC and PMCX. I've been working on it for about a week without luck.  I'm trying to implement SVMC in the recently released pmcx, but am running into issues with the size of the matrix fed into the run command.

In python, when I feed a volume processed by SVMC in MATLAB, I see: 

Exception has occurred: ValueError

Invalid volume dims for double volume.

File "\\PATH\test_cylinder.py", line 26, in <module> res=pmcx.run(mcfg) ValueError: Invalid volume dims for double volume.

Can you describe what kind of data type and size that output should be? I tried matching the output when processing in Matlab, which is 8x[dims] in size. Below I show exactly how I do this, but I have Python casting them directly from MATLAB as double variables. Casting the variable into uint8 to align with the comment above to write an svmc volume to a file from matlab as uint8 also did not work. I've tried saving as .bin from matlab and loading into python, as well as the more direct version outlined below by calling the matlab mcxsvmc.m function directly from python. Any breadcrumbs on this appreciated.

Vincent

Troubleshooting below ----

To avoid rewriting mcxsvmc.m into python, I am calling the Matlab mcxsvmc.m function on a volume made in python using the Matlab 2022b engine for python. I show this functionality in the code below.  I attach the .pickle files with variables "vol", which gets sent through mcxsvmc, and mcfg.pickle which contains the rest of my pmcx cfg commands so you can run the below commands.  I also give instructions for installing the matlab.engine further below in case you don't use it.

import pickle

import pmcx

import numpy as np

import matlab.engine

# Load unsmoothed volume

with open(r'\\tsn.tno.nl\RA-Data\SV\sv-070492\07_MCSimulaties\userWorkspaces\Zoutenbier\temp\volBeforeMCXSVMC.pickle', 'rb') as handle:

    vol = pickle.load(handle)

print('Loaded original volume')

# Start MATLAB

eng = matlab.engine.start_matlab()

print('Starting SVMC in matlab')

 # Call SVMC

svmcvol = eng.mcxsvmc(vol,'smoothing',1)

   

# Exit MATLAB

eng.quit()

cfg={}

print('Loading cfg for pmcx')

with open(r'\\tsn.tno.nl\RA-Data\SV\sv-070492\07_MCSimulaties\userWorkspaces\Zoutenbier\temp\mcfg.pickle', 'rb') as handle:

    cfg = pickle.load(handle)

cfg['mediabyte']="svmc"

# Reshape array from matlab connector

cfg['vol']=np.array(svmcvol).reshape(svmcvol.size, order='F')

#Run pmcx

pmcx.run(cfg)

print('hold')

Running this code works for me with the attached .pickles until the pmcx.run command.  I was able to get mcxsvmc working; I could pass a numpy array in, and transform the matlabarray type coming out also into a numpy array.  

I could also query the size of the volume created by the mcxsvmc.m process after reshaping, that's:

cfg['vol'].shape
(8, 100, 200, 100)

This seems correct to me.

I dont see what's wrong.

cfg dict passed into pmcx:

'vol':
array([[[[1., 1., 1., ..., 1., 1., 1.],
         [1., 1., 1., ..., 1., 1., 1.],
         [1., 1., 1., ..., 1., 1., 1.],
         ..., (TRUNCATED for display)
'dim': (100, 200, 100)
'prop': array([[24.4 ,  0.1 ,  0.9 ,  1.41],
       [24.4 ,  0.1 ,  0.9 ,  1.41],
       [24.4 ,  0.1 ,  0.9 ,  1.41],
       [24.4 ,  0.1 ,  0.9 ,  1.41],
       [ 0.  ,  0.  ,  1.  ,  1.38]])
'issrcfrom0': 1
'unitinmm': 0.005
'mediabyte': svmc'
'tstart': 0.0
'tend': 5e-09
'tstep': 5e-09
'srctype': 'line'
'srcpos': [50.0, 0.0, 50.0]
'srcdir': [0, 1, 0]
'srcparam1': [0.0, 200.0, 0.0, 0]
'srcparam2': [0, 0, 0, 0]
'nphoton': 500000000.0

'bc': '______000000'
'seed': 99999

Installation of matlab functionality in python: 

with matlab 2022b installed

python -m pip install matlabengine

 

I was able to prove functionality with:

import matlab.engine
eng = matlab.engine.start_matlab()
tf = eng.isprime(37)
print(tf)

this returns true for me (as it should!)

more info: <https://nl.mathworks.com/help/matlab/matlab_external/call-matlab-functions-from-python.html>)

Op donderdag 7 juli 2022 om 06:37:13 UTC+2 schreef q.fang:

[V Zoutenbier]

Sorry, I didn't attach the link for the .pickle download.  Available through my university file sharing service:

https://surfdrive.surf.nl/files/index.php/s/8EJe1GKBYxP7y6y

Op woensdag 29 maart 2023 om 19:53:41 UTC+2 schreef V Zoutenbier:

[Qianqian Fang]

hi Vincent

I don't have matlab.engine on my machine, so I can't test this. But one possible issue I can see is that you are passing on a c-ordererd (row-major) 3D numpy array to matlab (which is FORTRAN- or column-major order) to call mcxsvmc.m, I think this can create problems.

please recreate your vol raw data in the 'F' order and see if that can fix the issue

https://numpy.org/doc/stable/reference/generated/numpy.asfortranarray.html

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/6a5eaf2a-61e4-4764-a019-1f459b669825n%40googlegroups.com.