Fwd: Fwd: MCX simulation of SFDI

[Qianqian Fang]

---------- Forwarded message ---------
From: JR <jingxu...@gmail.com>
Date: Fri, Jan 8, 2021 at 12:55 AM
Subject: MCX simulation of SFDI
To: <mcx-...@googlegroups.com>

Hi Qianqian,

Happy new year!

I'm trying to simulate spatially modulated sources with MCXlab and compare the result with SFDI theory. However, there seems an underestimation of the decay rate of the relative reflectance for MCX result. The basic idea to generate the reflectance (Rd) in MCX is as follows:

1. create a homogeneous slab with the size of 100x30x20 mm^3;

2. use the 'fourier' mode of the MCX source to simulate the fluence with three phases (120 deg interval);

3. for each spatial frequency, demodulate the fluence to get the Rd (i.e., peak-peak value of the sinusoids).

I've attached the code of this part below. I was wondering if such comparison has been made previously? Any possible reason for this? I know that maybe the slab size is limited so there's potential partial volume effect and cannot be directly compared with SFDI theory, but it would be great to know if you have any comments.

Thank you,

Jingxuan

[Qianqian Fang]

hi Jingxuan,

I've forwarded your message to my student Rahul Ragunathan, who had done comparisons between MCX and diffusion with SFDI sources. He will follow up with your question shortly.

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/3dded8d7-0702-7dd6-c52e-3d77496ace8c%40neu.edu.

[JR]

Thank you very much!

Jingxuan

[Rahul Ragunathan]

Hello Jingxuan,

 

Sorry for jumping around messages, but this way I can post to the group as well as keep Dr .Fang in the loop.

 

So I found multiple issues and there is an additional clarification I want to add:

 

1. The optical properties in the comparison file that was sent had an error (musp=12 vs a more realistic musp=1.2).  I have attached the corrected plot to this email. Note I also include a white MC PB+Hankel Transform LUT derived from an external software (used in “appSFDI” ). This is the green crosses  in the attached plots. I also normalized to fx and included a separate plot. Note that all 3 plots share very similar trends and that MCX provides an intermediate output between the Diffusion approximation and the openSFDI forward model. This can be explained by the boundary effect of the slab vs semi-infinite( very large slab in practice) assumption and the fact that I ran lower photon counts in wide field just for testing (4e7).

 

2. As it pertains to your code, there is an error in your fx=0 term. Specifically, because the cfg.normflag output is set to 1, the phase shifts in a DC image that ordinarily equate to an image of [1 0.25 0.25] for phases are set to [1 1 1] instead. I suggest normalizing manually by the mean of your 3 patterns. Note that since your normalization is by this fx=0 term, that your entire output curve would change. I have included an image showing that your 3 DC images effectively have the same exact intensity.

 

3. Dr.Fang may have a better explanation for this, but I noticed that you were looking at flux.data at 1 voxel into your src domain rather than flux.dref at the 0 padded layer. Also, your 0 padded layer was misplaced at (:,:,1) even though your src term is shining at the opposite surface. This would also cause a deviation from the analytical approximation because you are effectively sinking your output fluence into your medium rather than looking at what is actually reflected on the source side.

 

I have attached my modified version of the script as well as the aforementioned images. Dr. Fang might have additional input regarding the third point in specific.

 

Please let me know if this is useful and if you have any additional questions!

 

Rahul

[JR]

Hi Rahul,

Thank you for your comments. I've gone through the code and all your comments. Here is something I'd like to discuss with you.

1. I was using mus = 12 and g=0.9, so the musp = 1.2, which is the same as your simulation. Is there a way to directly set musp in MCX?

2. Sorry that I didn't explain it well in the previous email. For the fx=0 term, I did demodulate it by the mean of 3 patterns. But thanks for pointing out the phase assignment issue, which should give me a better understanding of MCX. Do you mean cfg.isnormalized when you talk about cfg.normflag? I didn't find cfg.normflag in the ReadMe file. By the way, the note for cfg.isnormalized is

[1]-normalize the output fluence to unitary source, 0-no reflection.

What does the "no reflection" mean is this case? I'm assuming that isnormalized=1 makes the output fluence as the ratio of the photon weight with respect to a unitary source. Could you explain why this is affecting the source pattern phase?

3. Good catch! I miss-placed the air pad on the opposite side of the volume. I'm re-running my pipeline with corrections. I'll update the result with you when it's ready.

BTW, could you explain the difference between (a) fluence.data at the boundary layer and (b) fluence.dref at the air pad? My understanding is that (a) is the photon weight collected in the boundary layer, which is still the medium. However, (b) is the photon weight collected outside of the volume. If so, what is the acceptance angle of the detector for (b)?

Thank you very much and these discussions will really help me understand MCX.

Jingxuan

[Rahul Ragunathan]

Hello Jingxuan,

 

1. I see your point for the first statement. What you are doing is what I do in m simulations as well. I just did not catch it.
2. Sorry! I did mean cfg.isnormalized. I accidentally used the name of my variable in my function instead. So essentially, for SFDI we are projecting a pattern in the form of (cos(2*pi*fx+alpha)+1)/2. For fx=0 and alpha=[0 2p/3 4pi/3], this simplifies to intensities of [1 0.25 0.25]. When cfg.isnormalized is set to 1, you are normalizing your source patterns to a unitary source which means you are actually obtaining the results of [1 1 1] for your 3 phases. This is why in my attached image for the 3 DC terms, the fluence output was exactly the same despite the phase change ( we would expect [1 0.25 0.25] in terms of the fluences). You can simply run one full DC simulation at phase=0 and divide the output by 4 for your other 2 images.
3. I am unsure of the acceptance angle of photons present in the boundary layer. Dr. Fang, would likely be able to answer this better.

 

Thank you and I look forward to looking at the results!

 

Rahul

[JR]

Hi Rahul,

Thank you for the explanation. Here are my new results with dref and corrected air padding. Although the MCX result is still not matching the diffusion theory, I believe the problem is the normalization. Below is my analysis.

What I calculated for the MCX R_d is actually the demodulated reflectance M_ac. The unit is the default unit of flux.dref from MCX and the quantity is around 1e4 (see the unNormalized figure). However, in the diffusion model, R_d is a ratio normalized by its source power (amplitude of the sinusoidal). I naturally took the assumption that the source power of MCX is the same for each fx if I set cfg.isnormalized=1. That's why I normalize the M_ac by the DC reflectance from MCX and compare them with the diffusion theory (see the Normalized figure). So my questions are:

1. What is the source power when I set cfg.isnormalized=1 and what unit is it? According to your previous email, a pattern of (cos(2*pi*fx+alpha)+1)/2 is projected. So the source amplitude is 1? If so, why is the demodulated M_ac is around 1e4? In your previous email, you showed the R_d of MCX. Could you explain how you convert the M_ac to R_d?

2. What is the source power when I set cfg.isnormalized=0?

Thank you very much!

Jingxuan

[Rahul Ragunathan]

Hello Jingxuan,

 

I think a more direct way of addressing your comment is to list the normalization steps I used in my code. Note that I always keep my voxel size at 1 mm, and there likely is an area/volume scaling factor that would have to be accounted for in the script. Finding where that factor should be incorporated should be straightforward ( I had it working at one point). So for 1 mm voxel size:

 

1. I run all patterns with cfg.isnormalized set to 0,
2. I scale each image by the following PRIOR to demodulation for each of the 3 phases. Note the AreaBox term is constant and is equal to the source area ( aka number of voxels. I suspect with voxel size change , this should be in mm^2/mm^3 but please check!). croppedimages=flux.dref of your air layer in this case!
3.  

AreaScaling=1/AreaBox ; % Scaling for Area Projected

scalingfactor=numphotons*AreaScaling;

I1=croppedimages(:,:,1)./scalingfactor;

I2=croppedimages(:,:,2)./scalingfactor;

I3=croppedimages(:,:,3)./scalingfactor;

4. Because our patterns are discretized, I then normalize by the mean of all 3 patterns ( all 3 phases for a given fx). Please see the following chunk of code.

Note that this is the mean of the normalized sinusoidal patterns of each phase ( 0 to 1). This should also account for fx=0, where your intensities of the pattern will be [1 0.25 0.25].  You will therefore divide by 0.5. This term is necessary as our patterns are discretized and we are simulating over a slab where a phase shift will not retain consistent src power over the entire surface. I found this to be the best match. Since you used the fourier option as your source, the mean of your pattern may need to be approximated over your surface independently.

 

%% EQN 20 Cuccia

term1=sqrt(2)/3;

term2=(I1-I2).^2;

term3=(I2-I3).^2;

term4=(I3-I1).^2;

MAC.reflectance_d=(term1*((term2+term3+term4).^(1/2)))./mean(meanpatterns);  % This should roughly be equal to R_d derived analytically or through white MC using a pencil beam convolution

 

Please let me know if you have any other questions.

 

Rahul

From: JR <jingxu...@gmail.com>
Sent: Thursday, January 14, 2021 4:58 PM
To: Rahul Ragunathan <raguna...@northeastern.edu>
Cc: Fang, Qianqian <q.f...@northeastern.edu>; mcx-...@googlegroups.com
Subject: Re: [mcx-users] Fwd: Fwd: MCX simulation of SFDI

 

Hi Rahul,

 

Thank you for the explanation. Here are my new results with dref and corrected air padding. Although the MCX result is still not matching the diffusion theory, I believe the problem is the normalization. Below is my analysis.

 

What I calculated for the MCX R_d is actually the demodulated reflectance M_ac. The unit is the default unit of flux.dref from MCX and the quantity is around 1e4 (see the unNormalized figure). However, in the diffusion model, R_d is a ratio normalized by its source power (amplitude of the sinusoidal). I naturally took the assumption that the source power of MCX is the same for each fx if I set cfg.isnormalized=1. That's why I normalize the M_ac by the DC reflectance from MCX and compare them with the diffusion theory (see the Normalized figure). So my questions are:

 

1. What is the source power when I set cfg.isnormalized=1 and what unit is it? According to your previous email, a pattern of (cos(2*pi*fx+alpha)+1)/2 is projected. So the source amplitude is 1? If so, why is the demodulated M_ac is around 1e4? In your previous email, you showed the R_d of MCX. Could you explain how you convert the M_ac to R_d?

 

2. What is the source power when I set cfg.isnormalized=0?

 

Thank you very much!

 

Jingxuan

 

 

 

 

On Wed, Jan 13, 2021 at 6:45 PM Rahul Ragunathan <raguna...@northeastern.edu> wrote:

Hello Jingxuan,

 

1.      I see your point for the first statement. What you are doing is what I do in m simulations as well. I just did not catch it.

2.      Sorry! I did mean cfg.isnormalized. I accidentally used the name of my variable in my function instead. So essentially, for SFDI we are projecting a pattern in the form of (cos(2*pi*fx+alpha)+1)/2. For fx=0 and alpha=[0 2p/3 4pi/3], this simplifies to intensities of [1 0.25 0.25]. When cfg.isnormalized is set to 1, you are normalizing your source patterns to a unitary source which means you are actually obtaining the results of [1 1 1] for your 3 phases. This is why in my attached image for the 3 DC terms, the fluence output was exactly the same despite the phase change ( we would expect [1 0.25 0.25] in terms of the fluences). You can simply run one full DC simulation at phase=0 and divide the output by 4 for your other 2 images.

3.      I am unsure of the acceptance angle of photons present in the boundary layer. Dr. Fang, would likely be able to answer this better.

[JR]

Hello Rahul,

I tried your scaling steps but eventually found that those scaling factors are constant for all the fx because of the illumination area setup in my simulation.

1. This time, I also ran the forward simulation with cfg.isnormalized set to 0.

2. My source area is 100 mm*30 mm for all the spatial freq. With 0.5 mm voxel size, the grid size of the illumination area is 200*60. Here's how I define the source pattern.

    cfg.srcparam1 = [60, 0, 0, phase]; 
    cfg.srcparam2 = [0, 200, 0, kx];   % kx is number of cycles = fx*100

3. Following your scaling process, the ScaleFactor is actually constant for all fx.

    AreaBox = 200*60;  % area of illumination
    AreaScale=1/AreaBox ; % Scaling for Area Projected
    ScaleFactor=cfg.nphoton*AreaScale;

4. Since I was using a src area = 100mm, there are always integer numbers of cycles for the spatial frequency I used here. The src power will retain constant for different phase shifts. And the calculated means of all three phases are always 0.5 for those fx. I've attached a figure to show the sample dref (and therefore reflecting the src pattern) for fx=0.01 /mm. The phase shift should not change the total illumination power. Therefore, the scaling by mean(meanpatterns) does not change the relative R_d among different fx.

However, the absolute R_d does make more sense by these scaling steps (see attached). So the problem is still not solved in my case. I suspect that my assignment of srcparam may have issues. Since your mean(meanpatterns) is not constant for different fx, could you share the way you assign the src pattern?

Thank you,

Jingxuan

[Rahul Ragunathan]

 

I elected to use cfg.srctype=’pattern’ in order to utilize MCX’s photon sharing method and so that I could use other pattern types other than fourier using the same base code. The readme and help command should be useful but in general:

 

Cfg.srcpattern= patterns; à patterns is a 3D array of the form (pat idx, x,y) . I use the standard SFDI equation to solve over this grid using a 0.1 spacing. This can be different than your voxel spacing.)

Cfg.srcparam1=[height 0 0 size(patterns,2)] à height is src height in voxels

Cfg.srcparam2=[width 0 0 size(patterns,3)] à width is src width in voxels

Cfg.srcnum=size(patterns,1);

 

I compute the mean manually for each pattern when I create them. Yes, the mean is roughly 0.5 for all AC patterns as expected. The real purpose of the mean(mean(patterns))) is for the DC term as discussed before. I still feel that the issue is due to the change in voxel size not being accounted for. It seems that the R_d is consistently off what would be expected by a factor of 2 in your case which happens to be what the voxel size factor is ( 0.5 vs 1). Just curious, could you try to make your voxel size 0.25 mm^3 and tell me what you observe? Also is 1 mm^3 working fine?

 

Please let me know and hopefully what I have sent can be useful in solving this issue.

 

Thank you!

[Rahul Ragunathan]

Sorry about the double post.  Disregard my scaling factor comment as it does not appear to be a consistent 2x difference . Please do try the srcpattern change though.  I do not have much experience using the fourier srctype , so I wonder if this could potentially fix the problem.

 

Thank you!

Rahul

 

From: JR <jingxu...@gmail.com>
Sent: Friday, January 15, 2021 3:33 AM

[JR]

Thanks for the suggestion. I'll try the srcpattern. Just want to confirm: when you say "use the standard SFDI equation to solve over this grid using a 0.1 spacing", do you mean that

1. You calculate cfg.srcpattern (a 3D floating array) for each fx to create sinusoidal patterns with the strength of each voxel between 0 and 1?

AND

2. Such pattern is calculated on a volume with grid size of 0.1?

Thanks,

Jingxuan

[Christian Crane]

Hey Rahul!

I'm currently working with the SlabSFDI_RRmodified code in matlab. I just wanted to understand why your source position parameter is [0,0,23]./scale I understand scaling this based on voxel size. This may be a dumb question but why 23 of all values? Should this value scale with volume dimension like Lz etc..?

Thanks for any responses I appreciate it! I'll also be trying to understand and adapt the code further for my understanding and use in SFDI simulations and may end up posting further updates or files as time passes. 

~Christian

University of Maine Biomedical Engineering