[Homer-users] Modified Beer Lambert Law Code [Matlab]

[thuppert]

	A while ago, there were a number of people who requested a MBLL code.  I'm
sorry it took so long to get back to you on this.  I was in Italy for the
Human Brain Mapping conference for almost a month and I just got back last
week.  I am just getting caught up now on emails etc.

There really isn't any short code that I can give you to do the MBLL from
HOMER.  Most of HOMER is book-keeping to get the variables into a structure
that the program uses.  I have a code that would do the MBLL if the data was
in the proper "HOMER" structure, but I don't think that is very useful for
general use.

Instead, below you will find a sample script which calculates HbO & HbR from
dOD data using the MBLL.  This code makes use of a file called
"GetExtinctions.m" (which is used in HOMER).  I placed this file on an FTP
site.
https://www.nmr.mgh.harvard.edu/facility/filedrop/showgroup/2762/1/world



Ted Huppert, M.Sc.

PhD student-Harvard Univ.
Dept of Biophysics
Photon Migration Imaging lab
Mass General Hospital/CNY

Tele: (617)726-9338

thuppert at nmr.mgh.harvard.edu


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here is an example of the MBLL:
% This can be pasted into the Matlab command line.
% MAKE SURE YOU DOWNLOAD THE GetExtinctions.m FILE FIRST
% AND THAT THIS FILE IS IN THE MATLAB PATH

%Sample Data
I_w1=rand(10,1)+1;
I_w2=rand(10,1)+1;
ppf=60; dpf=6; L=3;


dOD_w1 = -log(I_w1/mean(I_w1));
dOD_w2 = -log(I_w2/mean(I_w2));

E = GetExtinctions([830 690]);  %E = [e at 830-HbO e at 830-HbR e at 830-Lipid
e at 830-H2O e at 830-AA;
                                %     e at 690-HbO e at 690-HbR e at 690-Lipid
e at 690-H2O e at 690-AA];

E=E(1:2,1:2);   %Only keep HbO/HbR parts

dOD_w1_L = dOD_w1 * ppf/(dpf*L);        %Pre-divide by pathlength (dpf,pvc
etc)
dOD_w2_L = dOD_w2 * ppf/(dpf*L);

dOD_L = [dOD_w1_L dOD_w2_L];            %Concatinate the 2(or more)
wavelengths

%I put pathlength into dOD_L so that I can preform one matrix inversion
rather than                                          %one per #measurements.
You could do inv(E*L) instead.

Einv = inv(E'*E)*E';            %Linear inversion operator (for 2 or more
wavelengths)

HbO_HbR = Einv * dOD_L';
 %Solve for HbO and HbR (This is the least-squares solution for unlimited #
of wavelengths)

HbO = HbO_HbR(1,:);
HbR = HbO_HbR(2,:);