Skip to content

Commit

Permalink
seperate Mie scattering functions from mcx_utils.h
Browse files Browse the repository at this point in the history
  • Loading branch information
@ShijieYan
ShijieYan committed Jan 27, 2022
1 parent 9b74e4b commit c350c67
Show file tree
Hide file tree
Showing 5 changed files with 372 additions and 291 deletions.
2 changes: 1 addition & 1 deletion src/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ CPPOPT=-g -Wall -std=c99 #-DNO_LZMA # -DUSE_OS_TIMER
OBJSUFFIX=.o
EXESUFFIX=

FILES=mcx_core mcx_utils mcx_shapes tictoc mcextreme mcx_bench cjson/cJSON ubj/ubjw
FILES=mcx_core mcx_utils mcx_shapes tictoc mcextreme mcx_bench mcx_mie cjson/cJSON ubj/ubjw

ifeq ($(findstring CYGWIN,$(PLATFORM)), CYGWIN)
CC=nvcc
Expand Down
319 changes: 319 additions & 0 deletions src/mcx_mie.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,319 @@
/***************************************************************************//**
** \mainpage Monte Carlo eXtreme - GPU accelerated Monte Carlo Photon Migration
**
** \author Qianqian Fang <q.fang at neu.edu>
** \copyright Qianqian Fang, 2009-2021
**
** \section sref Reference:
** \li \c (\b Fang2009) Qianqian Fang and David A. Boas,
** <a href="http://www.opticsinfobase.org/abstract.cfm?uri=oe-17-22-20178">
** "Monte Carlo Simulation of Photon Migration in 3D Turbid Media Accelerated
** by Graphics Processing Units,"</a> Optics Express, 17(22) 20178-20190 (2009).
** \li \c (\b Yu2018) Leiming Yu, Fanny Nina-Paravecino, David Kaeli, and Qianqian Fang,
** "Scalable and massively parallel Monte Carlo photon transport
** simulations for heterogeneous computing platforms," J. Biomed. Optics,
** 23(1), 010504, 2018. https://doi.org/10.1117/1.JBO.23.1.010504
** \li \c (\b Yan2020) Shijie Yan and Qianqian Fang* (2020), "Hybrid mesh and voxel
** based Monte Carlo algorithm for accurate and efficient photon transport
** modeling in complex bio-tissues," Biomed. Opt. Express, 11(11)
** pp. 6262-6270. https://doi.org/10.1364/BOE.409468
**
** \section slicense License
** GPL v3, see LICENSE.txt for details
*******************************************************************************/

/***************************************************************************//**
\file mcx_utils.c
@brief Mie scattering parameters handling for polarized light simulations
*******************************************************************************/

#include <math.h>
#include <stdlib.h>
#include "mcx_utils.h"
#include "mcx_mie.h"
#include "mcx_const.h"

/**
* @brief Precompute scattering parameters based on Mie theory [bohren and huffman]
*
* For each combination of sphere and background medium, compute the scattering
* efficiency and scattering Mueller matrix w.r.t. different scattering angles.
*
* @param[in] x: sphere particle size parameters
* @param[in] m: complex relative refractive index
* @param[in] mu: precomputed cosine of sampled scattering angles
* @param[out] smatrix: scattering Mueller matrix
* @param[out] qsca: scattering efficiency
*/

void Mie(double x, double _Complex m, const double *mu, float4 *smatrix, double *qsca, double *g){
double _Complex *D,*s1,*s2;
double _Complex z1=0.+I*0.;
double _Complex an,bn,bnm1,anm1;
double *pi0,*pi1,*tau;
double _Complex xi,xi0,xi1;
double psi0,psi1;
double alpha,beta,factor;
long n,k,nstop,sign;

if(x<=0.0) MCX_ERROR(-6,"sphere size must be positive");
if(x>20000.0) MCX_ERROR(-6,"spheres with x>20000 are not validated");

if((creal(m)==0.0 && x<0.1) || (creal(m)>0.0 && cabs(m)*x<0.1)){
small_Mie(x,m,mu,smatrix,qsca,g);
return;
}

nstop=floor(x+4.05*pow(x,0.33333)+2.0);

s1=(double _Complex *)calloc(NANGLES,sizeof(double _Complex));
s2=(double _Complex *)calloc(NANGLES,sizeof(double _Complex));
pi0=(double *)calloc(NANGLES,sizeof(double));
tau=(double *)calloc(NANGLES,sizeof(double));
pi1=(double *)calloc(NANGLES,sizeof(double));
for(int i=0;i<NANGLES;i++) pi1[i]=1.0;

if(creal(m)>0.0){
double _Complex z=x*m;
D=(double _Complex *)calloc(nstop+1,sizeof(double _Complex));
if(fabs(cimag(m)*x)<((13.78*creal(m)-10.8)*creal(m)+3.9))
Dn_up(z,nstop,D);
else
Dn_down(z,nstop,D);
}

psi0=sin(x);
psi1=psi0/x-cos(x);
xi0=psi0+I*cos(x);
xi1=psi1+I*(cos(x)/x+sin(x));
*qsca=0.0;
*g=0.0;
sign=1;
anm1=0.0+I*0.0;
bnm1=0.0+I*0.0;

for(n=1;n<=nstop;n++){
if(creal(m)==0.0){
an=(n*psi1/x-psi0)/(n/x*xi1-xi0);
bn=psi1/xi1;
}else if(cimag(m)==0.0){
z1=creal(D[n])/creal(m)+n/x+I*cimag(z1);
an=(creal(z1)*psi1-psi0)/(creal(z1)*xi1-xi0);

z1=creal(D[n])*creal(m)+n/x+I*cimag(z1);
bn=(creal(z1)*psi1-psi0)/(creal(z1)*xi1-xi0);
}else{
z1=D[n]/m;
z1+=n/x;
an=(creal(z1)*psi1-psi0+I*cimag(z1)*psi1)/(z1*xi1-xi0);

z1=D[n]*m;
z1+=n/x;
bn=(creal(z1)*psi1-psi0+I*cimag(z1)*psi1)/(z1*xi1-xi0);
}

for(k=0;k<NANGLES;k++){
factor=(2.0*n+1.0)/(n+1.0)/n;
tau[k]=n*mu[k]*pi1[k]-(n+1)*pi0[k];
alpha=factor*pi1[k];
beta=factor*tau[k];
s1[k]+=alpha*creal(an)+beta*creal(bn)+I*(alpha*cimag(an)+beta*cimag(bn));
s2[k]+=alpha*creal(bn)+beta*creal(an)+I*(alpha*cimag(bn)+beta*cimag(an));
}

for(k=0;k<NANGLES;k++){
factor=pi1[k];
pi1[k]=((2.0*n+1.0)*mu[k]*pi1[k]-(n+1.0)*pi0[k])/n;
pi0[k]=factor;
}

factor=2.0*n+1.0;
*g+=(n-1.0/n)*(creal(anm1)*creal(an)+cimag(anm1)*cimag(an)+creal(bnm1)*creal(bn)+cimag(bnm1)*cimag(bn));
*g+=factor/n/(n+1.0)*(creal(an)*creal(bn)+cimag(an)*cimag(bn));
*qsca+=factor*(cabs(an)*cabs(an)+cabs(bn)*cabs(bn));
sign*=-1;

factor=(2.0*n+1.0)/x;
xi=factor*xi1-xi0;
xi0=xi1;
xi1=xi;

psi0=psi1;
psi1=creal(xi1);

anm1=an;
bnm1=bn;
}
/* compute scattering efficiency and smatrix */
(*qsca)*=2.0/(x*x);
(*g)*=4.0/(*qsca)/(x*x);
for(int i=0;i<NANGLES;i++){
smatrix[i].x=0.5*cabs(s2[i])*cabs(s2[i])+0.5*cabs(s1[i])*cabs(s1[i]);
smatrix[i].y=0.5*cabs(s2[i])*cabs(s2[i])-0.5*cabs(s1[i])*cabs(s1[i]);
smatrix[i].z=creal(conj(s1[i])*s2[i]);
smatrix[i].w=cimag(conj(s1[i])*s2[i]);
}

if(creal(m)>0.0) free(D);
free(s1);
free(s2);
free(pi0);
free(pi1);
free(tau);
}

/**
* @brief Precompute scattering parameters for small particles
* @param[in] x: sphere particle size parameters
* @param[in] m: complex relative refractive index
* @param[in] mu: precomputed cosine of sampled scattering angles
* @param[out] smatrix: scattering Mueller matrix
* @param[out] qsca: scattering efficiency
*/

void small_Mie(double x, double _Complex m, const double *mu, float4 *smatrix, double *qsca, double *g){
double _Complex ahat1,ahat2,bhat1;
double _Complex z0,m2,m4;
double x2,x3,x4;

m2=m*m;
m4=m2*m2;
x2=x*x;
x3=x2*x;
x4=x2*x2;
z0=-cimag(m2)+I*(creal(m2)-1.0);
{
double _Complex z1,z2,z3,z4,D;

if(creal(m)==0.0){
z3=0.0+I*(2.0/3.0*(1.0-0.2*x2));
D=1.0-0.5*x2+I*(2.0/3.0*x3);
}else{
z1=2.0/3.0*z0;
z2=1.0-0.1*x2+(4.0*creal(m2)+5.0)*x4/1400.0+I*(4.0*x4*cimag(m2)/1400.0);
z3=z1*z2;

z4=x3*(1.0-0.1*x2)*z1;
D=2.0+creal(m2)+(1-0.7*creal(m2))*x2+(8*creal(m4)-385*creal(m2)+350.0)/1400*x4+creal(z4)+
I*((-0.7*cimag(m2))*x2+(8*cimag(m4)-385*cimag(m2))/1400*x4+cimag(z4));
}
ahat1=z3/D;
}
{
double _Complex z2,z6,z7;
if(creal(m)==0.0){
bhat1=(0.0+I*(-(1.0-0.1*x2)/3.0))/(1+0.5*x2+I*(-x3/3.0));
}else{
z2=x2/45.0*z0;
z6=1.0+(2.0*creal(m2)-5.0)*x2/70.0+I*(cimag(m2)*x2/35.0);
z7=1.0-(2.0*creal(m2)-5.0)*x2/30.0+I*(-cimag(m2)*x2/15.0);
bhat1=z2*(z6/z7);
}
}
{
double _Complex z3,z8;

if(creal(m)==0.0){
ahat2=0.0+I*x2/30.0;
}else{
z3=(1.0-x2/14)*x2/15.0*z0;
z8=2.0*creal(m2)+3.0-(creal(m2)/7.0-0.5)*x2+
I*(2.0*cimag(m2)-cimag(m2)/7.0*x2);
ahat2=z3/z8;
}
}
{
double T;

T=cabs(ahat1)*cabs(ahat1)+cabs(bhat1)*cabs(bhat1)+5.0/3.0*cabs(ahat2)*cabs(ahat2);
*qsca=6.0*x4*T;
*g=(creal(ahat1)*(creal(ahat2)+creal(bhat1))+cimag(ahat1)*(cimag(ahat2)+cimag(bhat1)))/T;
}
{
double muj,angle;
double _Complex s1,s2;

x3*=1.5;
ahat1*=x3;
bhat1*=x3;
ahat2*=x3*5.0/3.0;
for(int j=0;j<NANGLES;j++){
muj=mu[j];
angle=2*muj*muj-1;
s1=ahat1+(bhat1+ahat2)*muj;
s2=bhat1+(ahat1+ahat2)*angle;

smatrix[j].x=0.5*cabs(s2)*cabs(s2)+0.5*cabs(s1)*cabs(s1);
smatrix[j].y=0.5*cabs(s2)*cabs(s2)-0.5*cabs(s1)*cabs(s1);
smatrix[j].z=creal(conj(s1)*s2);
smatrix[j].w=cimag(conj(s1)*s2);
}
}
}

/**
* @brief
* @param z
* @param n
*/

double _Complex Lentz_Dn(double _Complex z,long n){
double _Complex alpha_j1,alpha_j2,zinv,aj;
double _Complex alpha,ratio,runratio;

zinv=2.0/z;
alpha=(n+0.5)*zinv;
aj=(-n-1.5)*zinv;
alpha_j1=aj+1.0/alpha;
alpha_j2=aj;
ratio=alpha_j1/alpha_j2;
runratio=alpha*ratio;

do{
aj=zinv-aj;
alpha_j1=1.0/alpha_j1+aj;
alpha_j2=1.0/alpha_j2+aj;
ratio=alpha_j1/alpha_j2;
zinv=-zinv;
runratio*=ratio;
}while(fabs(cabs(ratio)-1.0)>1e-12);

return -n/z+runratio;
}

/**
* @brief
* @param z
* @param nstop
* @param D
*/

void Dn_up(double _Complex z, long nstop, double _Complex *D){
double _Complex zinv,k_over_z;
zinv=1.0/z;

D[0]=1.0/ctan(z);
for(long k=1;k<nstop;k++){
k_over_z=k*zinv;
D[k]=1.0/(k_over_z-D[k-1])-k_over_z;
}
}

/**
* @brief
* @param z
* @param nstop
* @param D
*/

void Dn_down(double _Complex z, long nstop, double _Complex *D){
double _Complex zinv,k_over_z;
zinv=1.0/z;

D[nstop-1]=Lentz_Dn(z,nstop);
for(long k=nstop-1;k>=1;k--){
k_over_z=k*zinv;
D[k-1]=k_over_z-1.0/(D[k]+k_over_z);
}
}
51 changes: 51 additions & 0 deletions src/mcx_mie.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
/***************************************************************************//**
** \mainpage Monte Carlo eXtreme - GPU accelerated Monte Carlo Photon Migration
**
** \author Qianqian Fang <q.fang at neu.edu>
** \copyright Qianqian Fang, 2009-2021
**
** \section sref Reference:
** \li \c (\b Fang2009) Qianqian Fang and David A. Boas,
** <a href="http://www.opticsinfobase.org/abstract.cfm?uri=oe-17-22-20178">
** "Monte Carlo Simulation of Photon Migration in 3D Turbid Media Accelerated
** by Graphics Processing Units,"</a> Optics Express, 17(22) 20178-20190 (2009).
** \li \c (\b Yu2018) Leiming Yu, Fanny Nina-Paravecino, David Kaeli, and Qianqian Fang,
** "Scalable and massively parallel Monte Carlo photon transport
** simulations for heterogeneous computing platforms," J. Biomed. Optics,
** 23(1), 010504, 2018. https://doi.org/10.1117/1.JBO.23.1.010504
** \li \c (\b Yan2020) Shijie Yan and Qianqian Fang* (2020), "Hybrid mesh and voxel
** based Monte Carlo algorithm for accurate and efficient photon transport
** modeling in complex bio-tissues," Biomed. Opt. Express, 11(11)
** pp. 6262-6270. https://doi.org/10.1364/BOE.409468
**
** \section slicense License
** GPL v3, see LICENSE.txt for details
*******************************************************************************/

/***************************************************************************//**
\file mcx_mie.h
@brief MCX Mie scattering functions header
*******************************************************************************/

#ifndef _MCEXTREME_MIE_H
#define _MCEXTREME_MIE_H

#include <complex.h>
#include <vector_types.h>

#ifdef __cplusplus
extern "C" {
#endif

void Mie(double x, double _Complex m, const double *mu, float4 *smatrix, double *qsca, double *g);
void small_Mie(double x, double _Complex m, const double *mu, float4 *smatrix, double *qsca, double *g);
double _Complex Lentz_Dn(double _Complex z,long n);
void Dn_up(double _Complex z, long nstop, double _Complex *D);
void Dn_down(double _Complex z, long nstop, double _Complex *D);

#ifdef __cplusplus
}
#endif

#endif
Loading

0 comments on commit c350c67

Please sign in to comment.