3D simulations / photons initial location

[akipulkkinen]

In 3D codes, the photons location is generate on a triangle with an intent of having photon uniformly distributed it:

double w0 = UnifOpen(), w1 = UnifOpen(), w2 = UnifOpen(); phot->pos[0] = (w0 r(BH(ib, 0), 0) + w1 r(BH(ib, 1), 0) + w2 r(BH(ib, 2), 0)) / (w0 + w1 + w2); phot->pos[1] = (w0 r(BH(ib, 0), 1) + w1 r(BH(ib, 1), 1) + w2 r(BH(ib, 2), 1)) / (w0 + w1 + w2); phot->pos[2] = (w0 r(BH(ib, 0), 2) + w1 r(BH(ib, 1), 2) + w2 * r(BH(ib, 2), 2)) / (w0 + w1 + w2);

however, this does not produce uniform distribution. Correct form would be

double r1 = UnifOpen(), r2 = UnifOpen(); phot->pos[0] = (1.0 - sqrt(r1)) r(BH(ib, 0), 0) + sqrt(r1) (1.0 - r2) r(BH(ib, 1), 0) + sqrt(r1) r2) r(BH(ib, 2), 0); phot->pos[1] = (1.0 - sqrt(r1)) r(BH(ib, 0), 1) + sqrt(r1) (1.0 - r2) r(BH(ib, 1), 1) + sqrt(r1) r2) r(BH(ib, 2), 1); phot->pos[2] = (1.0 - sqrt(r1)) r(BH(ib, 0), 2) + sqrt(r1) (1.0 - r2) r(BH(ib, 1), 2) + sqrt(r1) r2) * r(BH(ib, 2), 2);

see, e.g. https://www.cs.princeton.edu/~funk/tog02.pdf and eq. (1)

[aaleino]

The source code now includes this fix.

[Marisadat]

hi dr.alexi

I used your last source code for 3D simulation of skin tissue and got my results. Will the simulation results in the modified code change significantly?

Thanks a lot

[akipulkkinen]

The significance of the issue depends on the size of the problem domain, light source definition, and the properties of the medium.

In situations where the problem is optically diffusive and exitance or fluence are studied far from the source the impact can be low.

My recommendation is to assess the severity by running simulation with the old code and the fixed code and compare the results.

Visually the difference between the two approaches for photons initial location generation can be compared in MATLAB (code below, figure attached).

Visuals show that the old approach generates non-uniform distribution of photons initial location. This translates to non-uniform irradiance profile (same as in the figure) on any boundary element acting as a light source.

MATLAB code to generate the attached figure follows:

% Triangle
p1 = [0, 0];
p2 = [1, 0];
p3 = [0, 1];

% Number of photons and discretization of histogram
N = 100000000;
dh = 0.001;
M = ceil(1 / dh);

% Old / incorrect approach 
r1 = rand(N, 1);
r2 = rand(N, 1);
r3 = rand(N, 1);
r = (r1 * p1 + r2 * p2 + r3 * p3) ./ repmat(r1 + r2 + r3, [1, 2]);

% Its histogram
ix = min(1 + floor(r(:, 1) / dh), M);
iy = min(1 + floor(r(:, 2) / dh), M);
C = sparse(sub2ind([M, M], iy, ix), 1, 1, M^2, 1);
C = reshape(full(C), [M, M]);
C1 = C;

% New / correct approach
r1 = rand(N, 1);
r2 = rand(N, 1);
xi = 1 - sqrt(r1);
eta = r2 .* (1 - xi);
r = (1 - xi - eta) * p1 + xi * p2 + eta * p3;

% Its histogram
ix = min(1 + floor(r(:, 1) / dh), M);
iy = min(1 + floor(r(:, 2) / dh), M);
C = sparse(sub2ind([M, M], iy, ix), 1, 1, M^2, 1);
C = reshape(full(C), [M, M]);
C2 = C;

% Normalize histograms into pdfs
C1 = C1 / sum(dh^2 * C1(:));
C2 = C2 / sum(dh^2 * C2(:));

% Compare visually
xvec =  linspace(0, 1, M);
yvec =  linspace(0, 1, M);
figure;
subplot(121); 
imagesc(xvec, yvec, C1); 
axis xy equal tight; 
colorbar; 
title('OLD/INCORRDECT APPROACH');
subplot(122);
imagesc(xvec, yvec, C2); 
axis xy equal tight; 
colorbar;
title('NEW/CORRDECT APPROACH');

[Marisadat]

Thank you for your response. I checked my simulation with the new code and compared the results. My problem is a 2mm2mm0.82mm cube with a small cylinder in it and the medium is skin. My light source is a circular section in upper boundary. I put two detectors in either side of the light source as reflectance detectors and one in lower boundary as transmittance detector. The difference in light exitance between the two codes for the reflective detectors is in fourth digit and for the transittance detector is in fifth digit. Also, for some elements this difference is positive and for others it is negative, so it is not possible to say with certainty in which state more light has reached the detectors. I have taken the results in different cylinder modes and different wavelengths with the previous code and I think because they are close to the new code results I can consider the same. Am I right?