Can't use pattern and detector simultaneously

[Giovanna Tramontin]

Hello Professor Fang, 

I have a problem with mmclab. I am trying to simulate a slab with an inclusion and, apparently, the functionality of the code depends on the status of cfg.gpuid. 

If -1, them mmc is used to run the code and the patterns works just fine (as you can see in the plot) but the detector does not, meaning that there is no output for detp. 

If gpuid is 1, them openCL is used to run the code and the detector works just fine, but the patterns are not illuminated into the sample. 

I need to use both the patterns and the detector. But, I only have this problem for cfg.srctype = 'planar' or 'pattern'. With Fourier, I can use both ate the same time with OpenCl. Can you help me? 

Here is the code, you can try changing cfg.gpuid to -1 and to 1 and look at the outputs. 

%
clear cfg;
clear all;

%% Geometry    
[node1, face1] = meshabox([0 0 0],[65 65 15], 2, 2);
[node2, face2] = meshacylinder([20 0+1e-1 12], [20 65-1e-1 12], 2, 2, 2);  
[node, face] = mergemesh(node1, face1, node2, face2);
[cfg.node, cfg.elem] = surf2mesh(node, face, [0 0 0], [65 65 15], 1, 2, [1,1,1;21,1,13], [], 0);


%% Source
cfg.srcpos = [10 30 -1];
cfg.srcdir = [0 0 1];  
cfg.srctype = 'pattern';
   
%Generate patterns:
val(:,:,1) = [1 1; 1 1];
val(:,:,2) = [1 1; 0 0];
val(:,:,3) = [1 0; 1 0];
val(:,:,4) = [1 0; 0 1];
M=4;

cfg.srcpattern = val;
cfg.srcparam1=[20 0 0 2];
cfg.srcparam2=[0 20 0 2];

srcdef=struct('srctype',cfg.srctype,'srcpos',cfg.srcpos,'srcdir',cfg.srcdir,...
'srcparam1',cfg.srcparam1,'srcparam2',cfg.srcparam2);

[cfg.node,cfg.elem] = mmcaddsrc(cfg.node,cfg.elem,...
mmcsrcdomain(srcdef,[min(cfg.node);max(cfg.node)]));

%Detector

cfg.detpos = [10.0, 30.0, 16.0, 0.0];
cfg.detdir = [0 0 -1];
cfg.detparam1 = [20.0, 0.0, 0.0, 20];
cfg.detparam2 = [0.0, 20.0, 0.0, 20];
detpos = cfg.detpos(1:3);
detdef = struct('srctype', 'planar', 'srcpos', detpos, 'srcdir', cfg.detdir, ...
         'srcparam1', cfg.detparam1, 'srcparam2', cfg.detparam2);
                   
[cfg.node, cfg.elem] = mmcadddet(cfg.node, cfg.elem, detdef);

cfg.elemprop=cfg.elem(:,5);
cfg.elem=cfg.elem(:,1:4);

cfg.e0 = tsearchn(cfg.node,cfg.elem,cfg.srcpos);

%General config
cfg.unitinmm = 1;
cfg.nphoton = 1e6;
cfg.srcpos = [10 30 -1];
cfg.srcdir = [0 0 1];
cfg.tsart = 0;                
cfg.tend = 5e-9;            
cfg.tstep = 5e-10;            
cfg.debuglevel = 'TP';
cfg.method='elem';
cfg.prop = [0, 0, 1, 1; 0.1, 0.001, 0.89, 1.37; 1, 0.01, 0.89, 1.37];

%If 1 use OpenCL, if -1, classical MMC  - THIS CHANGES THE OUTPUT
cfg.gpuid = -1;
   
%% Run simulation
cfg.issaveexit = 1; %1 2 and 3 changes the output in detp. 1 and 3 are the same and 2 is different
cfg.outputtype = 'fluence';
[output, detp] = mmclab(cfg);
   
data = output.data;
cwdata = squeeze(sum(data,3));
   
%% prepare subplots
nCol = round(sqrt(M));
nRow = round(M/nCol);
figure(10),figure(11)
for i = 1:M
     figure(10)
     subplot(nRow,nCol,i); title(['output pattern ', num2str(i)]);
     qmeshcut(cfg.elem(cfg.elemprop>0,1:4),cfg.node,(cwdata(i,:)'),'z=14','linestyle','none');
     figure(11)
     subplot(nRow,nCol,i); title(['input pattern ', num2str(i)]);
     plotmesh([cfg.node,cwdata(i,1:length(cfg.node))'],cfg.elem(cfg.elemprop>0,1:4),'z>0','linestyle','none');
     shading interp
     axis on
%     view([0,0,-1])
end


      

[Qianqian Fang]

hi Giovanna,

unfortunately, photon-sharing -i.e. simulating multiple pattern sources simultaneously- was only implemented in the CPU version of mmc, but not the GPU/OpenCL version.

I just created a ticket to add this to the TODO list. please keep an eye on this ticket for future updates.

https://github.com/fangq/mmc/issues/86

Qianqian

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/f57b0ede-50ad-423e-bc16-e93f4bd3fa2bn%40googlegroups.com.

[Giovanna Tramontin]

Thank you for your quickly reply! 

But the patter is working fine with gpuid = -1. The problem happens when I put the wide detector, because them I receive an empty output for detp after the simulation.  

Is the wide detector implemented for gpu and cpu? 

 

Giovanna

[Qianqian Fang]

hi Giovanna,

thanks for the report. I was able to reproduce the issue - when running on the CPU (cfg.gpuid=-1), but the detp outputs are empty. I also noticed that not just dept, seeds output were also empty if requested. both outputs were non-empty if using GPU

I created the below bug report

https://github.com/fangq/mmc/issues/87

and was also able to fix it with this one-line change (as part of a jumble commit in which I added support for saving photon trajectories in mmc last night)

https://github.com/fangq/mmc/commit/ebefe98577bee6b9507fb695c4763ad67245c61d#diff-76885cc6201a298c5ef20e311fb6c88a7fc2cb3745c84d02830524656cae9d95R310

please download the github or nightly build from this folder and test again

https://mcx.space/nightly/

let me know if you see any other problem - please also be aware that photon sharing is still not supported in the GPU.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/f770efa2-4a09-44ba-a1a7-0c61907a894an%40googlegroups.com.

[Giovanna Tramontin]

Hii Professor Fang, 

Thank you for fixing this issue! 

I am facing another issue now, which is that:

- MMC. I send 4 patterns as in the photon_sharing example with cfg.gpu = -1 and issavexit=1.
- I have a planar detector slightly outside the diffusive slab.

When I run the measurement in the detp structure I get many data such as weights and exit positions, but they seem to be as if the source is a single uniform planar source. Weights are 1 and no way to discriminate the output due to the 4 patterns.

If I put issavexit=2 I get in detp only a field Data which is correctly shaped as the detector I defined with also time dimension. The only problem is that, as before, the output is due to a planar uniform source. No weights of the photon derived from the different patterns…

How can I solve this issue? Is there a way to get information about the weights of the photon per pattern?

Best wishes,

Giovanna

[Qianqian Fang]

hi Giovanna,

I debugged this problem and noticed that the CPU code does not return either the launched photon weight (if single pattern) or the 1D index of the pattern pixel at which the photon was launched (if photon sharing, like we did in mcxlab). As a result, the output detp.w0 field is always 1 in all pattern sources.

I just made the following update to mmc, and now detp.w0 should return the launched pattern weight when using a single pattern

https://github.com/fangq/mmc/commit/da6d04e9f1f2a4890ac4708b00ca66b97363805c

or return the 1D integer index (starting from 0) of the cfg.srcpattern[] 2D array (assumed to be column-major order, like in matlab - that is index 0 is pixel [1,1], index 1 is pixel [2,1], etc) - in the photon-sharing mode. this is now consistent with the output of mcxlab in photon sharing mode.

using the returned index, you should be able to look up the launched values (as cfg.srcpattern(detp.w0(i)+1) or shifting by k*Nx*Ny for the (k+1)-th pattern),

when used in combination with wide-field det (i.e. setting issavedet=2), I believe at least single-pattern case should work.

let me know if this is the case after you update your copy to the latest version.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/04b539d6-2edd-4890-b6ae-b3a615a43479n%40googlegroups.com.