From a467813709359d73f41e00929a216988cd06c68e Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 21 May 2026 15:55:57 -0700 Subject: [PATCH 1/3] Added missing dMC detector processing to both load_results_script.m and loadMCResults.m, and corrected normalization for dMCdROfFxdMua and dMCdROfFxdMus detectors. --- matlab/post_processing/loadMCResults.m | 140 +++++++++++++++++-- matlab/post_processing/load_results_script.m | 60 +++++++- 2 files changed, 184 insertions(+), 16 deletions(-) diff --git a/matlab/post_processing/loadMCResults.m b/matlab/post_processing/loadMCResults.m index 98cdfe5b6..c0f2e0228 100644 --- a/matlab/post_processing/loadMCResults.m +++ b/matlab/post_processing/loadMCResults.m @@ -815,7 +815,46 @@ - imag(FluenceOfFxAndZ.Mean) .* imag(FluenceOfFxAndZ.Mean)) / json.N); end results{di}.FluenceOfFxAndZ = FluenceOfFxAndZ; - case 'RadianceOfRhoAndZAndAngle' + case 'RadianceOfRhoAtZ' + RadianceOfRhoAtZ.Name = detector.Name; + RadianceOfRhoAtZ.ZDepth = detector.ZDepth; + RadianceOfRhoAtZ.ZDirection = detector.ZDirection; + tempRho = detector.Rho; + RadianceOfRhoAtZ.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); + RadianceOfRhoAtZ.Rho_Midpoints = (RadianceOfRhoAtZ.Rho(1:end-1) + RadianceOfRhoAtZ.Rho(2:end))/2; + RadianceOfRhoAtZ.Mean = readBinaryData([datadir slash detector.Name],length(RadianceOfRhoAtZ.Rho)-1); + if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) + RadianceOfRhoAtZ.SecondMoment = readBinaryData([datadir slash detector.Name '_2'], ... + length(RadianceOfRhoAtZ.Rho)-1); + RadianceOfRhoAtZ.Stdev = sqrt((RadianceOfRhoAtZ.SecondMoment - (RadianceOfRhoAtZ.Mean .* RadianceOfRhoAtZ.Mean)) / json.N); + end + results{di}.RadianceOfRhoAtZ = RadianceOfRhoAtZ; + case 'RadianceOfRhoAndTimeAndMaxDepthAtZ' + RadianceOfRhoAndTimeAndMaxDepthAtZ.Name = detector.Name; + RadianceOfRhoAndTimeAndMaxDepthAtZ.ZDepth = detector.ZDepth; + RadianceOfRhoAndTimeAndMaxDepthAtZ.ZDirection = detector.ZDirection; + tempRho = detector.Rho; + tempTime = detector.Time; + tempMaxDepth = detector.MaxDepth; + RadianceOfRhoAndTimeAndMaxDepthAtZ.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); + RadianceOfRhoAndTimeAndMaxDepthAtZ.Rho_Midpoints = (RadianceOfRhoAndTimeAndMaxDepthAtZ.Rho(1:end-1) + RadianceOfRhoAndTimeAndMaxDepthAtZ.Rho(2:end))/2; + RadianceOfRhoAndTimeAndMaxDepthAtZ.Time = linspace((tempTime.Start), (tempTime.Stop), (tempTime.Count)); + RadianceOfRhoAndTimeAndMaxDepthAtZ.Time_Midpoints = (RadianceOfRhoAndTimeAndMaxDepthAtZ.Time(1:end-1) + RadianceOfRhoAndTimeAndMaxDepthAtZ.Time(2:end))/2; + RadianceOfRhoAndTimeAndMaxDepthAtZ.MaxDepth = linspace((tempMaxDepth.Start), (tempMaxDepth.Stop), (tempMaxDepth.Count)); + RadianceOfRhoAndTimeAndMaxDepthAtZ.MaxDepth_Midpoints = (RadianceOfRhoAndTimeAndMaxDepthAtZ.MaxDepth(1:end-1) + RadianceOfRhoAndTimeAndMaxDepthAtZ.MaxDepth(2:end))/2; + RadianceOfRhoAndTimeAndMaxDepthAtZ.Mean = readBinaryData([datadir slash detector.Name], ... + (length(RadianceOfRhoAndTimeAndMaxDepthAtZ.MaxDepth)-1)*(length(RadianceOfRhoAndTimeAndMaxDepthAtZ.Time)-1)*(length(RadianceOfRhoAndTimeAndMaxDepthAtZ.Rho)-1)); % read column major json binary + RadianceOfRhoAndTimeAndMaxDepthAtZ.Mean = reshape(RadianceOfRhoAndTimeAndMaxDepthAtZ.Mean, ... + [length(RadianceOfRhoAndTimeAndMaxDepthAtZ.MaxDepth)-1,length(RadianceOfRhoAndTimeAndMaxDepthAtZ.Time)-1,length(RadianceOfRhoAndTimeAndMaxDepthAtZ.Rho)-1]); % read column major json binary + if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) + RadianceOfRhoAndTimeAndMaxDepthAtZ.SecondMoment = readBinaryData([datadir slash detector.Name '_2'], ... + (length(RadianceOfRhoAndTimeAndMaxDepthAtZ.MaxDepth)-1)*(length(RadianceOfRhoAndTimeAndMaxDepthAtZ.Time)-1)*(length(RadianceOfRhoAndTimeAndMaxDepthAtZ.Rho)-1)); + RadianceOfRhoAndTimeAndMaxDepthAtZ.SecondMoment = reshape(RadianceOfRhoAndTimeAndMaxDepthAtZ.Mean, ... + [length(RadianceOfRhoAndTimeAndMaxDepthAtZ.MaxDepth)-1,length(RadianceOfRhoAndTimeAndMaxDepthAtZ.Time)-1,length(RadianceOfRhoAndTimeAndMaxDepthAtZ.Rho)-1]); + RadianceOfRhoAndTimeAndMaxDepthAtZ.Stdev = sqrt((RadianceOfRhoAndTimeAndMaxDepthAtZ.SecondMoment - (RadianceOfRhoAndTimeAndMaxDepthAtZ.Mean .* RadianceOfRhoAndTimeAndMaxDepthAtZ.Mean)) / json.N); + end + results{di}.RadianceOfRhoAndTimeAndMaxDepthAtZ = RadianceOfRhoAndTimeAndMaxDepthAtZ; + case 'RadianceOfRhoAndZAndAngle' RadianceOfRhoAndZAndAngle.Name = detector.Name; tempRho = detector.Rho; tempZ = detector.Z; @@ -1566,17 +1605,17 @@ dMCROfRhodMua.Stdev = sqrt((dMCdROfRhodMua.SecondMoment - (dMCdROfRhodMua.Mean .* dMCdROfRhodMua.Mean)) / (databaseInputJson.N)); end results{di}.dMCdROfRhodMua = dMCdROfRhodMua; - case 'dMCdROfRhodMus' - dMCdROfRhodMus.Name = detector.Name; + case 'dMCdROfRhoAndTimedMus' + dMCdROfRhoAndTimedMus.Name = detector.Name; tempRho = detector.Rho; - dMCdROfRhodMus.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); - dMCdROfRhodMus.Rho_Midpoints = (dMCdROfRhodMus.Rho(1:end-1) + dMCdROfRhodMus.Rho(2:end))/2; - dMCdROfRhodMus.Mean = readBinaryData([datadir slash detector.Name],length(dMCdROfRhodMus.Rho)-1); + dMCdROfRhoAndTimedMus.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); + dMCdROfRhoAndTimedMus.Rho_Midpoints = (dMCdROfRhoAndTimedMus.Rho(1:end-1) + dMCdROfRhoAndTimedMus.Rho(2:end))/2; + dMCdROfRhoAndTimedMus.Mean = readBinaryData([datadir slash detector.Name],length(dMCdROfRhoAndTimedMus.Rho)-1); if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) - dMCdROfRhodMus.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],length(dMCdROfRhodMus.Rho)-1); - dMCROfRhodMus.Stdev = sqrt((dMCdROfRhodMus.SecondMoment - (dMCdROfRhodMus.Mean .* dMCdROfRhodMus.Mean)) / (databaseInputJson.N)); + dMCdROfRhoAndTimedMus.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],length(dMCdROfRhoAndTimedMus.Rho)-1); + dMCROfRhoAndTimedMus.Stdev = sqrt((dMCdROfRhoAndTimedMus.SecondMoment - (dMCdROfRhoAndTimedMus.Mean .* dMCdROfRhoAndTimedMus.Mean)) / (databaseInputJson.N)); end - results{di}.dMCdROfRhodMus = dMCdROfRhodMus; + results{di}.dMCdROfRhoAndTimedMus = dMCdROfRhoAndTimedMus; case 'pMCROfRhoAndTime' pMCROfRhoAndTime.Name = detector.Name; tempRho = detector.Rho; @@ -1655,6 +1694,34 @@ pMCROfRhoAndTimeAndSubregionRecessed.ROfRhoStdev = sqrt((pMCROfRhoAndTimeAndSubregionRecessed.ROfRhoSecondMoment - (pMCROfRhoAndTimeAndSubregionRecessed.ROfRho .* pMCROfRhoAndTimeAndSubregionRecessed.ROfRho)) / (databaseInputJson.N)); end results{di}.pMCROfRhoAndTimeAndSubregionRecessed = pMCROfRhoAndTimeAndSubregionRecessed; + case 'dMCdROfRhoAndTimedMua' + dMCdROfRhoAndTimedMua.Name = detector.Name; + tempRho = detector.Rho; + dMCdROfRhoAndTimedMua.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); + dMCdROfRhoAndTimedMua.Rho_Midpoints = (dMCdROfRhoAndTimedMua.Rho(1:end-1) + dMCdROfRhoAndTimedMua.Rho(2:end))/2; + tempTime = detector.Time; + dMCdROfRhoAndTimedMua.Time = linspace((tempTime.Start), (tempTime.Stop), (tempTime.Count)); + dMCdROfRhoAndTimedMua.Time_Midpoints = (dMCdROfRhoAndTimedMua.Time(1:end-1) + dMCdROfRhoAndTimedMua.Time(2:end))/2; + dMCdROfRhoAndTimedMua.Mean = readBinaryData([datadir slash detector.Name],[length(dMCdROfRhoAndTimedMua.Time)-1,length(dMCdROfRhoAndTimedMua.Rho)-1]); + if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) + dMCdROfRhoAndTimedMua.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],[length(dMCdROfRhoAndTimedMua.Time)-1,length(dMCdROfRhoAndTimedMua.Rho)-1]); + dMCROfRhoAndTimedMua.Stdev = sqrt((dMCdROfRhoAndTimedMua.SecondMoment - (dMCdROfRhoAndTimedMua.Mean .* dMCdROfRhoAndTimedMua.Mean)) / (databaseInputJson.N)); + end + results{di}.dMCdROfRhoAndTimedMua = dMCdROfRhoAndTimedMua; + case 'dMCdROfRhoAndTimedMus' + dMCdROfRhoAndTimedMus.Name = detector.Name; + tempRho = detector.Rho; + dMCdROfRhoAndTimedMus.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); + dMCdROfRhoAndTimedMus.Rho_Midpoints = (dMCdROfRhoAndTimedMus.Rho(1:end-1) + dMCdROfRhoAndTimedMus.Rho(2:end))/2; + tempTime = detector.Time; + dMCdROfRhoAndTimedMus.Time = linspace((tempTime.Start), (tempTime.Stop), (tempTime.Count)); + dMCdROfRhoAndTimedMus.Time_Midpoints = (dMCdROfRhoAndTimedMus.Time(1:end-1) + dMCdROfRhoAndTimedMus.Time(2:end))/2; + dMCdROfRhoAndTimedMus.Mean = readBinaryData([datadir slash detector.Name],[length(dMCdROfRhoAndTimedMus.Time)-1,length(dMCdROfRhoAndTimedMus.Rho)-1]); + if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) + dMCdROfRhoAndTimedMus.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],[length(dMCdROfRhoAndTimedMus.Time)-1,length(dMCdROfRhoAndTimedMus.Rho)-1]); + dMCROfRhoAndTimedMus.Stdev = sqrt((dMCdROfRhoAndTimedMus.SecondMoment - (dMCdROfRhoAndTimedMus.Mean .* dMCdROfRhoAndTimedMus.Mean)) / (databaseInputJson.N)); + end + results{di}.dMCdROfRhoAndTimedMus = dMCdROfRhoAndTimedMus; case 'pMCROfXAndY' pMCROfXAndY.Name = detector.Name; tempX = detector.X; @@ -1731,7 +1798,9 @@ pMCROfFx.Fx = linspace((tempFx.Start), (tempFx.Stop), (tempFx.Count)); pMCROfFx.Fx_Midpoints = pMCROfFx.Fx; tempData = readBinaryData([datadir slash detector.Name],2*length(pMCROfFx.Fx)); - pMCROfFx.Mean = tempData(1:2:end) + 1i*tempData(2:2:end); + pMCROfFx.Mean = tempData(1:2:end) + 1i*tempData(2:2:end); + pMCROfFx.Amplitude = abs(pMCROfFx.Mean); + pMCROfFx.Phase = -angle(pMCROfFx.Mean); if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) tempData = readBinaryData([datadir slash detector.Name '_2'],2*length(pMCROfFx.Fx)); pMCROfFx.SecondMoment = tempData(1:2:end) + 1i*tempData(2:2:end); @@ -1739,6 +1808,57 @@ imag(pMCROfFx.Mean) .* imag(pMCROfFx.Mean)) / databaseInputJson.N); end results{di}.pMCROfFx = pMCROfFx; + case 'pMCROfFxAndTime' + pMCROfFxAndTime.Name = detector.Name; + tempFx = detector.Fx; + pMCROfFxAndTime.Fx = linspace((tempFx.Start), (tempFx.Stop), (tempFx.Count)); + pMCROfFxAndTime.Fx_Midpoints = pMCROfFxAndTime.Fx; + tempTime = detector.Time; + pMCROfFxAndTime.Time = linspace((tempTime.Start), (tempTime.Stop), (TempTime.Count)); + pMCROfFxAndTime.Time_Midpoints = pMCROfFxAndTime; + tempData = readBinaryData([datadir slash detector.Name],2*length(pMCROfFxAndTime.Fx)); + pMCROfFxAndTime.Mean = tempData(1:2:end) + 1i*tempData(2:2:end); + pMCROfFxAndTime.Amplitude = abs(pMCROfFxAndTime.Mean); + pMCROfFxAndTime.Phase = -angle(pMCROfFxAndTime.Mean); + if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) + tempData = readBinaryData([datadir slash detector.Name '_2'],2*length(pMCROfFxAndTime.Fx)); + pMCROfFxAndTime.SecondMoment = tempData(1:2:end) + 1i*tempData(2:2:end); + pMCROfFxAndTime.Stdev = sqrt((pMCROfFxAndTime.SecondMoment + real(pMCROfFxAndTime.Mean) .* real(pMCROfFxAndTime.Mean) + ... + imag(pMCROfFxAndTime.Mean) .* imag(pMCROfFxAndTime.Mean)) / databaseInputJson.N); + end + results{di}.pMCROfFxAndTime = pMCROfFxAndTime; + case 'dMCdROfFxdMua' + dMCdROfFxdMua.Name = detector.Name; + tempFx = detector.Fx; + dMCdROfFxdMua.Fx = linspace((tempFx.Start), (tempFx.Stop), (tempFx.Count)); + dMCdROfFxdMua.Fx_Midpoints = dMCdROfFxdMua.Fx; + tempData = readBinaryData([datadir slash detector.Name],2*length(dMCdROfFxdMua.Fx)); + dMCdROfFxdMua.Mean = tempData(1:2:end) + 1i*tempData(2:2:end); + dMCdROfFxdMua.Amplitude = abs(dMCdROfFxdMua.Mean); + dMCdROfFxdMua.Phase = -angle(dMCdROfFxdMua.Mean); + if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) + tempData = readBinaryData([datadir slash detector.Name '_2'],2*length(dMCdROfFxdMua.Fx)); + dMCdROfFxdMua.SecondMoment = tempData(1:2:end) + 1i*tempData(2:2:end); + dMCdROfFxdMua.Stdev = sqrt((dMCdROfFxdMua.SecondMoment + real(dMCdROfFxdMua.Mean) .* real(dMCdROfFxdMua.Mean) + ... + imag(dMCdROfFxdMua.Mean) .* imag(dMCdROfFxdMua.Mean)) / databaseInputJson.N); + end + results{di}.dMCdROfFxdMua = dMCdROfFxdMua; + case 'dMCdROfFxdMus' + dMCdROfFxdMus.Name = detector.Name; + tempFx = detector.Fx; + dMCdROfFxdMus.Fx = linspace((tempFx.Start), (tempFx.Stop), (tempFx.Count)); + dMCdROfFxdMus.Fx_Midpoints = dMCdROfFxdMus.Fx; + tempData = readBinaryData([datadir slash detector.Name],2*length(dMCdROfFxdMus.Fx)); + dMCdROfFxdMus.Mean = tempData(1:2:end) + 1i*tempData(2:2:end); + dMCdROfFxdMus.Amplitude = abs(dMCdROfFxdMus.Mean); + dMCdROfFxdMus.Phase = -angle(dMCdROfFxdMus.Mean); + if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) + tempData = readBinaryData([datadir slash detector.Name '_2'],2*length(dMCdROfFxdMus.Fx)); + dMCdROfFxdMus.SecondMoment = tempData(1:2:end) + 1i*tempData(2:2:end); + dMCdROfFxdMus.Stdev = sqrt((dMCdROfFxdMus.SecondMoment + real(dMCdROfFxdMus.Mean) .* real(dMCdROfFxdMus.Mean) + ... + imag(dMCdROfFxdMus.Mean) .* imag(dMCdROfFxdMus.Mean)) / databaseInputJson.N); + end + results{di}.dMCdROfFxdMus = dMCdROfFxdMus; case 'pMCTOfRho' pMCTOfRho.Name = detector.Name; tempRho = detector.Rho; diff --git a/matlab/post_processing/load_results_script.m b/matlab/post_processing/load_results_script.m index 8bc65da03..3e48052d6 100644 --- a/matlab/post_processing/load_results_script.m +++ b/matlab/post_processing/load_results_script.m @@ -73,7 +73,14 @@ show.pMCROfXAndYAndTimeAndSubregion = 1; show.pMCROfXAndYAndTimeAndSubregionRecessed = 1; show.pMCROfFx = 1; +show.pMCROfFxAndTime = 1; show.pMCTOfRho = 1; +show.dMCdROfRhodMua = 1; +show.dMCdROfRhodMus = 1; +show.dMCdROfRhoAndTimedMua = 1; +show.dMCdROfRhoAndTimedMus = 1; +show.dMCdROfFxdMua = 1; +show.dMCdROfFxdMus = 1; show.ReflectedMTOfRhoAndSubregionHist = 1; show.ReflectedMTOfXAndYAndSubregionHist = 1; show.TransmittedMTOfRhoAndSubregionHist = 1; @@ -1057,13 +1064,13 @@ rhodelta = results{di}.pMCROfRho.Rho(2)-results{di}.pMCROfRho.Rho(1); rhonorm = 2 * pi * results{di}.pMCROfRho.Rho_Midpoints * rhodelta; disp(['Total reflectance captured by pMCROfRho detector: ' num2str(sum(results{di}.pMCROfRho.Mean.*rhonorm'))]); - end + end if isfield(results{di}, 'pMCROfRhoRecessed') && show.pMCROfRhoRecessed figname = sprintf('log10(%s)',results{di}.pMCROfRhoRecessed.Name); figure; plot(results{di}.pMCROfRhoRecessed.Rho_Midpoints, log10(results{di}.pMCROfRhoRecessed.Mean)); title(figname); set(gcf,'Name', figname); xlabel('\rho [mm]'); ylabel('pMC R(\rho) [mm^-^2]'); rhodelta = results{di}.pMCROfRhoRecessed.Rho(2)-results{di}.pMCROfRhoRecessed.Rho(1); rhonorm = 2 * pi * results{di}.pMCROfRhoRecessed.Rho_Midpoints * rhodelta; disp(['Total reflectance captured by pMCROfRhoRecessed detector: ' num2str(sum(results{di}.pMCROfRhoRecessed.Mean.*rhonorm'))]); - end + end if isfield(results{di}, 'pMCROfRhoAndTime') && show.pMCROfRhoAndTime figname = sprintf('log10(%s)',results{di}.pMCROfRhoAndTime.Name); figure; imagesc(results{di}.pMCROfRhoAndTime.Rho_Midpoints, results{di}.pMCROfRhoAndTime.Time_Midpoints,log10(results{di}.pMCROfRhoAndTime.Mean)); colorbar; title(figname); set(gcf,'Name', figname);ylabel('time [ns]'); xlabel('\rho [mm]'); numtimes = length(results{di}.pMCROfRhoAndTime.Time)-1; @@ -1071,7 +1078,7 @@ rhodelta = results{di}.pMCROfRhoAndTime.Rho(2)-results{di}.pMCROfRhoAndTime.Rho(1); rhonorm = 2 * pi * results{di}.pMCROfRhoAndTime.Rho_Midpoints * rhodelta; disp(['Total reflectance captured by pMCROfRhoAndTime detector: ' num2str(sum(sum(timedelta*results{di}.pMCROfRhoAndTime.Mean.*repmat(rhonorm,[numtimes,1]))))]); - end + end if isfield(results{di}, 'pMCROfRhoAndTimeRecessed') && show.pMCROfRhoAndTimeRecessed figname = sprintf('log10(%s)',results{di}.pMCROfRhoAndTimeRecessed.Name); figure; imagesc(results{di}.pMCROfRhoAndTimeRecessed.Rho_Midpoints, results{di}.pMCROfRhoAndTimeRecessed.Time_Midpoints,log10(results{di}.pMCROfRhoAndTimeRecessed.Mean)); colorbar; title(figname); set(gcf,'Name', figname);ylabel('time [ns]'); xlabel('\rho [mm]'); numtimes = length(results{di}.pMCROfRhoAndTimeRecessed.Time)-1; @@ -1102,6 +1109,35 @@ rhonorm = 2 * pi * results{di}.pMCROfRhoAndTimeSubregionRecessed.Rho_Midpoints * rhodelta; disp(['Total reflectance captured by ROfRhoAndTimeAndSubregionRecessed detector: ' num2str(sum(sum(timedelta*results{di}.pMCROfRhoAndTimeAndSubregionRecessed.Mean.*repmat(rhonorm,[numtimes,1]))))]); end + if isfield(results{di}, 'dMCdROfRhodMua') && show.dMCdROfRhodMua + figname = sprintf('log10(%s)',results{di}.dMCdROfRhodMua.Name); figure; plot(results{di}.dMCdROfRhodMua.Rho_Midpoints, log10(results{di}.dMCdROfRhodMua.Mean)); title(figname); set(gcf,'Name', figname); xlabel('\rho [mm]'); ylabel('dMC dR(\rho)/dMua [mm^-^2]'); + rhodelta = results{di}.dMCdROfRhodMua.Rho(2)-results{di}.dMCdROfRhodMua.Rho(1); + rhonorm = 2 * pi * results{di}.dMCdROfRhodMua.Rho_Midpoints * rhodelta; + disp(['Total reflectance captured by dMCdROfRhodMua detector: ' num2str(sum(results{di}.dMCdROfRhodMua.Mean.*rhonorm'))]); + end + if isfield(results{di}, 'dMCdROfRhodMus') && show.dMCdROfRhodMus + figname = sprintf('log10(%s)',results{di}.dMCdROfRhodMus.Name); figure; plot(results{di}.dMCdROfRhodMus.Rho_Midpoints, log10(results{di}.dMCdROfRhodMus.Mean)); title(figname); set(gcf,'Name', figname); xlabel('\rho [mm]'); ylabel('dMC dR(\rho)/dMus [mm^-^2]'); + rhodelta = results{di}.dMCdROfRhodMus.Rho(2)-results{di}.dMCdROfRhodMus.Rho(1); + rhonorm = 2 * pi * results{di}.dMCdROfRhodMus.Rho_Midpoints * rhodelta; + disp(['Total reflectance captured by dMCdROfRhodMus detector: ' num2str(sum(results{di}.dMCdROfRhodMus.Mean.*rhonorm'))]); + end + if isfield(results{di}, 'dMCdROfRhoAndTimedMua') && show.dMCdROfRhoAndTimedMua + figname = sprintf('log10(%s)',results{di}.dMCdROfRhoAndTimedMua.Name); figure; imagesc(results{di}.dMCdROfRhoAndTimedMua.Rho_Midpoints, results{di}.dMCdROfRhoAndTimedMua.Time_Midpoints,log10(results{di}.dMCdROfRhoAndTimedMua.Mean)); colorbar; title(figname); set(gcf,'Name', figname);ylabel('time [ns]'); xlabel('\rho [mm]'); + numtimes = length(results{di}.dMCdROfRhoAndTimedMua.Time)-1; + timedelta = results{di}.dMCdROfRhoAndTimedMua.Time(2)-results{di}.dMCdROfRhoAndTimedMua.Time(1); + rhodelta = results{di}.dMCdROfRhoAndTimedMua.Rho(2)-results{di}.dMCdROfRhoAndTimedMua.Rho(1); + rhonorm = 2 * pi * results{di}.dMCdROfRhoAndTimedMua.Rho_Midpoints * rhodelta; + disp(['Total reflectance captured by dMCdROfRhoAndTimedMua detector: ' num2str(sum(sum(timedelta*results{di}.dMCdROfRhoAndTimedMua.Mean.*repmat(rhonorm,[numtimes,1]))))]); + end + if isfield(results{di}, 'dMCdROfRhoAndTimedMus') && show.dMCdROfRhoAndTimedMus + figname = sprintf('log10(%s)',results{di}.dMCdROfRhoAndTimedMus.Name); figure; imagesc(results{di}.dMCdROfRhoAndTimedMus.Rho_Midpoints, results{di}.dMCdROfRhoAndTimedMus.Time_Midpoints,log10(results{di}.dMCdROfRhoAndTimedMus.Mean)); colorbar; title(figname); set(gcf,'Name', figname);ylabel('time [ns]'); xlabel('\rho [mm]'); + numtimes = length(results{di}.dMCdROfRhoAndTimedMus.Time)-1; + timedelta = results{di}.dMCdROfRhoAndTimedMus.Time(2)-results{di}.dMCdROfRhoAndTimedMus.Time(1); + rhodelta = results{di}.dMCdROfRhoAndTimedMus.Rho(2)-results{di}.dMCdROfRhoAndTimedMus.Rho(1); + rhonorm = 2 * pi * results{di}.dMCdROfRhoAndTimedMus.Rho_Midpoints * rhodelta; + disp(['Total reflectance captured by dMCdROfRhoAndTimedMus detector: ' num2str(sum(sum(timedelta*results{di}.dMCdROfRhoAndTimedMus.Mean.*repmat(rhonorm,[numtimes,1]))))]); + end + if isfield(results{di}, 'pMCROfXAndY') && show.pMCROfXAndY figname = sprintf('log10(%s)',results{di}.pMCROfXAndY.Name); figure; imagesc(results{di}.pMCROfXAndY.X_Midpoints, results{di}.pMCROfXAndY.Y_Midpoints,log10(results{di}.pMCROfXAndY.Mean)); colorbar; title(figname); set(gcf,'Name', figname);ylabel('y [mm]'); xlabel('x [mm]'); xdelta = results{di}.ROfXAndYAndTime.X(2)-results{di}.ROfXAndYAndTime.X(1); @@ -1143,10 +1179,22 @@ end if isfield(results{di}, 'pMCROfFx') && show.pMCROfFx figname = sprintf('%s - Amplitude',results{di}.pMCROfFx.Name);figure;plot(results{di}.pMCROfFx.Fx_Midpoints, abs(results{di}.pMCROfFx.Mean));title(figname);set(gcf,'Name', figname);xlabel('f_x [/mm]');ylabel('R(f_x) [unitless]'); - Fxdelta = results{di}.pMCROfFx.Fx(2)-results{di}.pMCROfFx.Fx(1); - Fxnorm = 2 * pi * (results{di}.pMCROfFx.Fx_Midpoints * Fxdelta); - disp(['Total reflectance captured by ROfFx detector: ' num2str(sum(results{di}.pMCROfFx.Mean.*Fxnorm'))]); + disp(['Total reflectance captured by pMCROfFx detector: ' num2str(results{di}.pMCROfFx.Amplitude(1))]); end + if isfield(results{di}, 'dMCdROfFxdMua') && show.dMCdROfFxdMua + figname = sprintf('%s - Amplitude',results{di}.dMCdROfFxdMua.Name);figure;plot(results{di}.dMCdROfFxdMua.Fx_Midpoints, abs(results{di}.dMCdROfFxdMua.Mean));title(figname);set(gcf,'Name', figname);xlabel('f_x [/mm]');ylabel('dR(f_x)/dMua [unitless]'); + disp(['Total reflectance captured by dMCdROfFxdMua detector: ' num2str(results{di}.dMCdROfFxdMua.Amplitude(1))]); + end + if isfield(results{di}, 'dMCdROfFxdMus') && show.dMCdROfFxdMus + figname = sprintf('%s - Amplitude',results{di}.dMCdROfFxdMus.Name);figure;plot(results{di}.dMCdROfFxdMus.Fx_Midpoints, abs(results{di}.dMCdROfFxdMus.Mean));title(figname);set(gcf,'Name', figname);xlabel('f_x [/mm]');ylabel('dR(f_x)/dMua [unitless]'); + disp(['Total reflectance captured by dMCdROfFxdMus detector: ' num2str(results{di}.dMCdROfFxdMus.Amplitude(1))]); + end + if isfield(results{di}, 'pMCROfFxAndTime') && show.pMCROfFxAndTime + figname = sprintf('%s - Amplitude',results{di}.pMCROfFxAndTime.Name);figure;imagesc(results{di}.pMCROfFxAndTime.Fx_Midpoints, results{di}.pMCROfFxAndTime.Time_Midpoints, abs(results{di}.pMCROfFxAndTime.Mean));title(figname);set(gcf,'Name', figname);xlabel('f_x [/mm]');ylabel('time [ns]'),title('R(f_x) [unitless]'); + timedelta = results{di}.pMCROfFx.Time(2)-results{di}.pMCROfFx.Time(1); + disp(['Total reflectance captured by pMCROfFxAndTime detector: ' num2str(timedelta*results{di}.pMCROfFxAndTime.Amplitude(1))]); + end + if isfield(results{di}, 'pMCTOfRho') && show.pMCTOfRho figname = sprintf('log10(%s)',results{di}.pMCTOfRho.Name); figure; plot(results{di}.pMCTOfRho.Rho_Midpoints, log10(results{di}.pMCTOfRho.Mean)); title(figname); set(gcf,'Name', figname); xlabel('\rho [mm]'); ylabel('pMC T(\rho) [mm^-^2]'); disp(['Total reflectance captured by pMCTOfRho detector: ' num2str(sum(results{di}.pMCTOfRho.Mean(:)))]); From b6ed9ef3dfa217f7c0c739ca7d1a697894a9e2b3 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 21 May 2026 16:07:22 -0700 Subject: [PATCH 2/3] Editing snafu. Should be correct now. --- matlab/post_processing/loadMCResults.m | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/matlab/post_processing/loadMCResults.m b/matlab/post_processing/loadMCResults.m index c0f2e0228..26bd5ca8c 100644 --- a/matlab/post_processing/loadMCResults.m +++ b/matlab/post_processing/loadMCResults.m @@ -1602,20 +1602,20 @@ dMCdROfRhodMua.Mean = readBinaryData([datadir slash detector.Name],length(dMCdROfRhodMua.Rho)-1); if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) dMCdROfRhodMua.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],length(dMCdROfRhodMua.Rho)-1); - dMCROfRhodMua.Stdev = sqrt((dMCdROfRhodMua.SecondMoment - (dMCdROfRhodMua.Mean .* dMCdROfRhodMua.Mean)) / (databaseInputJson.N)); + dMCdROfRhodMua.Stdev = sqrt((dMCdROfRhodMua.SecondMoment - (dMCdROfRhodMua.Mean .* dMCdROfRhodMua.Mean)) / (databaseInputJson.N)); end results{di}.dMCdROfRhodMua = dMCdROfRhodMua; - case 'dMCdROfRhoAndTimedMus' - dMCdROfRhoAndTimedMus.Name = detector.Name; + case 'dMCdROfRhodMus' + dMCdROfRhodMus.Name = detector.Name; tempRho = detector.Rho; - dMCdROfRhoAndTimedMus.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); - dMCdROfRhoAndTimedMus.Rho_Midpoints = (dMCdROfRhoAndTimedMus.Rho(1:end-1) + dMCdROfRhoAndTimedMus.Rho(2:end))/2; - dMCdROfRhoAndTimedMus.Mean = readBinaryData([datadir slash detector.Name],length(dMCdROfRhoAndTimedMus.Rho)-1); + dMCdROfRhodMus.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); + dMCdROfRhodMus.Rho_Midpoints = (dMCdROfRhoAndTimedMus.Rho(1:end-1) + dMCdROfRhodMus.Rho(2:end))/2; + dMCdROfRhodMus.Mean = readBinaryData([datadir slash detector.Name],length(dMCdROfRhodMus.Rho)-1); if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) - dMCdROfRhoAndTimedMus.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],length(dMCdROfRhoAndTimedMus.Rho)-1); - dMCROfRhoAndTimedMus.Stdev = sqrt((dMCdROfRhoAndTimedMus.SecondMoment - (dMCdROfRhoAndTimedMus.Mean .* dMCdROfRhoAndTimedMus.Mean)) / (databaseInputJson.N)); + dMCdROfRhodMus.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],length(dMCdROfRhodMus.Rho)-1); + dMCdROfRhodMus.Stdev = sqrt((dMCdROfRhodMus.SecondMoment - (dMCdROfRhodMus.Mean .* dMCdROfRhodMus.Mean)) / (databaseInputJson.N)); end - results{di}.dMCdROfRhoAndTimedMus = dMCdROfRhoAndTimedMus; + results{di}.dMCdROfRhodMus = dMCdROfRhodMus; case 'pMCROfRhoAndTime' pMCROfRhoAndTime.Name = detector.Name; tempRho = detector.Rho; @@ -1705,7 +1705,7 @@ dMCdROfRhoAndTimedMua.Mean = readBinaryData([datadir slash detector.Name],[length(dMCdROfRhoAndTimedMua.Time)-1,length(dMCdROfRhoAndTimedMua.Rho)-1]); if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) dMCdROfRhoAndTimedMua.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],[length(dMCdROfRhoAndTimedMua.Time)-1,length(dMCdROfRhoAndTimedMua.Rho)-1]); - dMCROfRhoAndTimedMua.Stdev = sqrt((dMCdROfRhoAndTimedMua.SecondMoment - (dMCdROfRhoAndTimedMua.Mean .* dMCdROfRhoAndTimedMua.Mean)) / (databaseInputJson.N)); + dMCdROfRhoAndTimedMua.Stdev = sqrt((dMCdROfRhoAndTimedMua.SecondMoment - (dMCdROfRhoAndTimedMua.Mean .* dMCdROfRhoAndTimedMua.Mean)) / (databaseInputJson.N)); end results{di}.dMCdROfRhoAndTimedMua = dMCdROfRhoAndTimedMua; case 'dMCdROfRhoAndTimedMus' @@ -1719,7 +1719,7 @@ dMCdROfRhoAndTimedMus.Mean = readBinaryData([datadir slash detector.Name],[length(dMCdROfRhoAndTimedMus.Time)-1,length(dMCdROfRhoAndTimedMus.Rho)-1]); if(detector.TallySecondMoment && exist([datadir slash detector.Name '_2'],'file')) dMCdROfRhoAndTimedMus.SecondMoment = readBinaryData([datadir slash detector.Name '_2'],[length(dMCdROfRhoAndTimedMus.Time)-1,length(dMCdROfRhoAndTimedMus.Rho)-1]); - dMCROfRhoAndTimedMus.Stdev = sqrt((dMCdROfRhoAndTimedMus.SecondMoment - (dMCdROfRhoAndTimedMus.Mean .* dMCdROfRhoAndTimedMus.Mean)) / (databaseInputJson.N)); + dMCdROfRhoAndTimedMus.Stdev = sqrt((dMCdROfRhoAndTimedMus.SecondMoment - (dMCdROfRhoAndTimedMus.Mean .* dMCdROfRhoAndTimedMus.Mean)) / (databaseInputJson.N)); end results{di}.dMCdROfRhoAndTimedMus = dMCdROfRhoAndTimedMus; case 'pMCROfXAndY' From 4f1f337da79b010e8ac8b355bd3343e766d4c5a1 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 21 May 2026 16:20:03 -0700 Subject: [PATCH 3/3] Indentation fixes. --- matlab/post_processing/loadMCResults.m | 10 +++++----- matlab/post_processing/load_results_script.m | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/matlab/post_processing/loadMCResults.m b/matlab/post_processing/loadMCResults.m index 26bd5ca8c..2ddf68894 100644 --- a/matlab/post_processing/loadMCResults.m +++ b/matlab/post_processing/loadMCResults.m @@ -1565,7 +1565,7 @@ ReflectedTimeOfRhoAndSubregionHist.Stdev = sqrt((ReflectedTimeOfRhoAndSubregionHist.SecondMoment - (ReflectedTimeOfRhoAndSubregionHist.Mean .* ReflectedTimeOfRhoAndSubregionHist.Mean)) / (N)); end results{di}.ReflectedTimeOfRhoAndSubregionHist = ReflectedTimeOfRhoAndSubregionHist; - case 'pMCATotal' + case 'pMCATotal' pMCAtotal.Name = detector.Name; pMCATotal_txt = readAndParseJson([datadir slash detector.Name '.txt']); pMCATotal.Mean = pMCATotal_txt.Mean; @@ -1594,7 +1594,7 @@ pMCROfRhoRecessed.Stdev = sqrt((pMCROfRhoRecessed.SecondMoment - (pMCROfRhoRecessed.Mean .* pMCROfRhoRecessed.Mean)) / (databaseInputJson.N)); end results{di}.pMCROfRhoRecessed = pMCROfRhoRecessed; - case 'dMCdROfRhodMua' + case 'dMCdROfRhodMua' dMCdROfRhodMua.Name = detector.Name; tempRho = detector.Rho; dMCdROfRhodMua.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); @@ -1605,7 +1605,7 @@ dMCdROfRhodMua.Stdev = sqrt((dMCdROfRhodMua.SecondMoment - (dMCdROfRhodMua.Mean .* dMCdROfRhodMua.Mean)) / (databaseInputJson.N)); end results{di}.dMCdROfRhodMua = dMCdROfRhodMua; - case 'dMCdROfRhodMus' + case 'dMCdROfRhodMus' dMCdROfRhodMus.Name = detector.Name; tempRho = detector.Rho; dMCdROfRhodMus.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); @@ -1708,7 +1708,7 @@ dMCdROfRhoAndTimedMua.Stdev = sqrt((dMCdROfRhoAndTimedMua.SecondMoment - (dMCdROfRhoAndTimedMua.Mean .* dMCdROfRhoAndTimedMua.Mean)) / (databaseInputJson.N)); end results{di}.dMCdROfRhoAndTimedMua = dMCdROfRhoAndTimedMua; - case 'dMCdROfRhoAndTimedMus' + case 'dMCdROfRhoAndTimedMus' dMCdROfRhoAndTimedMus.Name = detector.Name; tempRho = detector.Rho; dMCdROfRhoAndTimedMus.Rho = linspace((tempRho.Start), (tempRho.Stop), (tempRho.Count)); @@ -1722,7 +1722,7 @@ dMCdROfRhoAndTimedMus.Stdev = sqrt((dMCdROfRhoAndTimedMus.SecondMoment - (dMCdROfRhoAndTimedMus.Mean .* dMCdROfRhoAndTimedMus.Mean)) / (databaseInputJson.N)); end results{di}.dMCdROfRhoAndTimedMus = dMCdROfRhoAndTimedMus; - case 'pMCROfXAndY' + case 'pMCROfXAndY' pMCROfXAndY.Name = detector.Name; tempX = detector.X; tempY = detector.Y; diff --git a/matlab/post_processing/load_results_script.m b/matlab/post_processing/load_results_script.m index 3e48052d6..5a748d061 100644 --- a/matlab/post_processing/load_results_script.m +++ b/matlab/post_processing/load_results_script.m @@ -76,7 +76,7 @@ show.pMCROfFxAndTime = 1; show.pMCTOfRho = 1; show.dMCdROfRhodMua = 1; -show.dMCdROfRhodMus = 1; +show.dMCdROfRhodMus = 1; show.dMCdROfRhoAndTimedMua = 1; show.dMCdROfRhoAndTimedMus = 1; show.dMCdROfFxdMua = 1; @@ -1129,7 +1129,7 @@ rhonorm = 2 * pi * results{di}.dMCdROfRhoAndTimedMua.Rho_Midpoints * rhodelta; disp(['Total reflectance captured by dMCdROfRhoAndTimedMua detector: ' num2str(sum(sum(timedelta*results{di}.dMCdROfRhoAndTimedMua.Mean.*repmat(rhonorm,[numtimes,1]))))]); end - if isfield(results{di}, 'dMCdROfRhoAndTimedMus') && show.dMCdROfRhoAndTimedMus + if isfield(results{di}, 'dMCdROfRhoAndTimedMus') && show.dMCdROfRhoAndTimedMus figname = sprintf('log10(%s)',results{di}.dMCdROfRhoAndTimedMus.Name); figure; imagesc(results{di}.dMCdROfRhoAndTimedMus.Rho_Midpoints, results{di}.dMCdROfRhoAndTimedMus.Time_Midpoints,log10(results{di}.dMCdROfRhoAndTimedMus.Mean)); colorbar; title(figname); set(gcf,'Name', figname);ylabel('time [ns]'); xlabel('\rho [mm]'); numtimes = length(results{di}.dMCdROfRhoAndTimedMus.Time)-1; timedelta = results{di}.dMCdROfRhoAndTimedMus.Time(2)-results{di}.dMCdROfRhoAndTimedMus.Time(1);