如何为两个人拟合 EELS 背景 window
How to fitting EELS background for two window
dm脚本提供了拟合指数背景的功能,例如:
image src := GetFrontImage()
number bEs = 340 // eV
number bEe = 390 // eV
image signal := EELSSubtractPowerLawBackground(src, bEs, bEe)
但是,对于微弱的信号,我们会想通过拟合两个window背景来提高S/N比率,例如:
有谁知道如何通过脚本来适配这种背景?
提前致谢。
不幸的是,不存在针对此的明确命令。
但是,可以通过模仿手动执行的操作来解决此问题:添加 ROI 以便显示。
下面的脚本可以满足您的需求,但它确实需要显示光谱并将移除其上的其他 ROI。
number b1Es = 400 // eV
number b1Ee = 510 // eV
number b2Es = 800 // eV
number b2Ee = 850 // eV
image src := GetFrontImage()
imageDisplay disp = src.ImageGetImageDisplay(0) // Script requires spectrum to be displayed
// Need to convert range values into channels by calibration
number scale, origin
string units
src.ImageGetDimensionCalibration(0,origin,scale,units,0)
if ( units != "eV" ) Throw( "Spectrum not calibrated in eV" )
number b1Es_ch = trunc( (b1Es - origin)/scale )
number b1Ee_ch = trunc( (b1Ee - origin)/scale )
number b2Es_ch = trunc( (b2Es - origin)/scale )
number b2Ee_ch = trunc( (b2Ee - origin)/scale )
// Remove any existing ROIs on the spectrum
while( disp.ImageDisplayCountROIS() )
disp.ImageDisplayDeleteROI( disp.ImageDisplayGetRoi(0))
ROI bkgr1 = NewROI()
bkgr1.ROISetRange(b1Es_ch,b1Ee_ch)
ROI sigr = NewROI()
sigr.ROISetRange(b1Ee_ch,b2Es_ch)
ROI bkgr2 = NewROI()
bkgr2.ROISetRange(b2Es_ch,b2Ee_ch)
// Add ROIs in order (like when adding by mouse)
disp.ImageDisplayAddROI(bkgr1) // BACKGROUND 1
disp.ImageDisplayAddROI(sigr) // SIGNAL
disp.ImageDisplayAddROI(bkgr2) // BACKGROUND 2 (and more...)
// Optionally: Change background model away from default
// ChooseMenuItem( "EELS", "Background Model", "Linear" )
image signal := src{"Edge"}.ImageClone() // copy edge-slice
// Remove ROIs again (cleanup)
while( disp.ImageDisplayCountROIS() )
disp.ImageDisplayDeleteROI( disp.ImageDisplayGetRoi(0))
// Show "signal"
signal.ShowImage()
感谢BmyGuest,但我想建立一个全自动程序。
幸运的是,我找到了一个简单的写下脚本的解决方案。
这是代码:
// reference:
// https://towardsdatascience.com/mathematics-for-machine-learning-linear-regression-least-square-regression-de09cf53757c
// key information:
// y = mx + b
// where m = sum ((xdata-xmean)*(ydata-ymean)) / sum ((xdata-xmean)**2)
// b = ymean - m*xmean
class FittingTwoWindowBkg : object {
FittingTwoWindowBkg(object self) {
Result("Script object 'FittingTwoWindowBkg' ["+self.ScriptObjectGetID().hex()+"] constructed.\n");
};
~FittingTwoWindowBkg(object self) {
Result("Script object 'FittingTwoWindowBkg' ["+self.ScriptObjectGetID().hex()+"] destructed.\n");
};
image TwoWindowBkgFitting(object self, image src, number Eng1, number Rng1, number Eng2, number Rng2) {
number Eorignal = src.ImageGetDimensionOrigin(2)
number Estep = src.ImageGetDimensionScale(2)
number sx, sy, sz
src.Get3DSize(sx, sy, sz)
image axis = src.ImageClone()
axis = iplane*Estep + Eorignal
axis.SetName("Energy axis")
// the power law bkg I = A E^-r can be simplified to
// the linear formula by : log(I) = -r*log(E) + log(A)
image data = src.ImageClone()
data = log(data)
data.SetName("log SI")
number ch1 = trunc((Eng1-Eorignal)/Estep)
number r1= trunc(Rng1/Estep)
number ch2 = trunc((Eng2-Eorignal)/Estep)
number r2= trunc(Rng2/Estep)
number nsz = r1+r2
// prepare for regression
image xdata := RealImage("xdata",4,sx,sy, nsz)
xdata.Slice3(0,0,0,0,sx,1,1,sy,1,2,r1,1) = axis.Slice3(0,0,ch1,0,sx,1,1,sy,1,2,r1,1)
xdata.Slice3(0,0,r1,0,sx,1,1,sy,1,2,r2,1) = axis.Slice3(0,0,ch2,0,sx,1,1,sy,1,2,r2,1)
image ydata := RealImage("ydata",4,sx,sy, nsz)
ydata.Slice3(0,0,0,0,sx,1,1,sy,1,2,r1,1) = data.Slice3(0,0,ch1,0,sx,1,1,sy,1,2,r1,1)
ydata.Slice3(0,0,r1,0,sx,1,1,sy,1,2,r2,1) = data.Slice3(0,0,ch2,0,sx,1,1,sy,1,2,r2,1)
image xmean := RealImage("xmean",4,sx,sy, nsz)
image ymean := RealImage("ymean",4,sx,sy, nsz)
for (number i=0; i<nsz; i++) {
// to repeat the mean value as a 3d stack
// the sum value of each pixel can be rapidly done by the Project function
xmean.Slice2(0,0,i,0,sx,1,1,sy,1) = Project(xdata,2) / nsz
ymean.Slice2(0,0,i,0,sx,1,1,sy,1) = Project(ydata,2) / nsz
};
// do regression
image m = Project((xdata-xmean)*(ydata-ymean), 2) / Project((xdata-xmean)**2, 2)
image b = ymean[icol,irow,0] - m*xmean[icol,irow,0]
// to apply y=mx+b to all the planes
image bkg := src.ImageClone()
for (number i=0; i<sz; i++) {
bkg.Slice2(0,0,i,0,sx,1,1,sy,1) = m*axis.Slice2(0,0,i,0,sx,1,1,sy,1)+b
};
// reconstruct to the power law bkg
bkg = exp(bkg)
return bkg
};
}
{
Object FTW = Alloc(FittingTwoWindowBkg)
number Eng1 = 1500 // eV, 1st Energy of fitting window
number Rng1 = 200 // eV, fitting range of 1st window
number Eng2 = 2600 // eV, 2nd Energy of fitting window
number Rng2 = 200 // eV, fitting range of 2nd window
image src := GetFrontImage()
image bkg := FTW.TwoWindowBkgFitting(src, Eng1, Rng1, Eng2, Rng2)
bkg.ShowImage()
};
dm脚本提供了拟合指数背景的功能,例如:
image src := GetFrontImage()
number bEs = 340 // eV
number bEe = 390 // eV
image signal := EELSSubtractPowerLawBackground(src, bEs, bEe)
但是,对于微弱的信号,我们会想通过拟合两个window背景来提高S/N比率,例如:
有谁知道如何通过脚本来适配这种背景? 提前致谢。
不幸的是,不存在针对此的明确命令。 但是,可以通过模仿手动执行的操作来解决此问题:添加 ROI 以便显示。
下面的脚本可以满足您的需求,但它确实需要显示光谱并将移除其上的其他 ROI。
number b1Es = 400 // eV
number b1Ee = 510 // eV
number b2Es = 800 // eV
number b2Ee = 850 // eV
image src := GetFrontImage()
imageDisplay disp = src.ImageGetImageDisplay(0) // Script requires spectrum to be displayed
// Need to convert range values into channels by calibration
number scale, origin
string units
src.ImageGetDimensionCalibration(0,origin,scale,units,0)
if ( units != "eV" ) Throw( "Spectrum not calibrated in eV" )
number b1Es_ch = trunc( (b1Es - origin)/scale )
number b1Ee_ch = trunc( (b1Ee - origin)/scale )
number b2Es_ch = trunc( (b2Es - origin)/scale )
number b2Ee_ch = trunc( (b2Ee - origin)/scale )
// Remove any existing ROIs on the spectrum
while( disp.ImageDisplayCountROIS() )
disp.ImageDisplayDeleteROI( disp.ImageDisplayGetRoi(0))
ROI bkgr1 = NewROI()
bkgr1.ROISetRange(b1Es_ch,b1Ee_ch)
ROI sigr = NewROI()
sigr.ROISetRange(b1Ee_ch,b2Es_ch)
ROI bkgr2 = NewROI()
bkgr2.ROISetRange(b2Es_ch,b2Ee_ch)
// Add ROIs in order (like when adding by mouse)
disp.ImageDisplayAddROI(bkgr1) // BACKGROUND 1
disp.ImageDisplayAddROI(sigr) // SIGNAL
disp.ImageDisplayAddROI(bkgr2) // BACKGROUND 2 (and more...)
// Optionally: Change background model away from default
// ChooseMenuItem( "EELS", "Background Model", "Linear" )
image signal := src{"Edge"}.ImageClone() // copy edge-slice
// Remove ROIs again (cleanup)
while( disp.ImageDisplayCountROIS() )
disp.ImageDisplayDeleteROI( disp.ImageDisplayGetRoi(0))
// Show "signal"
signal.ShowImage()
感谢BmyGuest,但我想建立一个全自动程序。 幸运的是,我找到了一个简单的写下脚本的解决方案。 这是代码:
// reference:
// https://towardsdatascience.com/mathematics-for-machine-learning-linear-regression-least-square-regression-de09cf53757c
// key information:
// y = mx + b
// where m = sum ((xdata-xmean)*(ydata-ymean)) / sum ((xdata-xmean)**2)
// b = ymean - m*xmean
class FittingTwoWindowBkg : object {
FittingTwoWindowBkg(object self) {
Result("Script object 'FittingTwoWindowBkg' ["+self.ScriptObjectGetID().hex()+"] constructed.\n");
};
~FittingTwoWindowBkg(object self) {
Result("Script object 'FittingTwoWindowBkg' ["+self.ScriptObjectGetID().hex()+"] destructed.\n");
};
image TwoWindowBkgFitting(object self, image src, number Eng1, number Rng1, number Eng2, number Rng2) {
number Eorignal = src.ImageGetDimensionOrigin(2)
number Estep = src.ImageGetDimensionScale(2)
number sx, sy, sz
src.Get3DSize(sx, sy, sz)
image axis = src.ImageClone()
axis = iplane*Estep + Eorignal
axis.SetName("Energy axis")
// the power law bkg I = A E^-r can be simplified to
// the linear formula by : log(I) = -r*log(E) + log(A)
image data = src.ImageClone()
data = log(data)
data.SetName("log SI")
number ch1 = trunc((Eng1-Eorignal)/Estep)
number r1= trunc(Rng1/Estep)
number ch2 = trunc((Eng2-Eorignal)/Estep)
number r2= trunc(Rng2/Estep)
number nsz = r1+r2
// prepare for regression
image xdata := RealImage("xdata",4,sx,sy, nsz)
xdata.Slice3(0,0,0,0,sx,1,1,sy,1,2,r1,1) = axis.Slice3(0,0,ch1,0,sx,1,1,sy,1,2,r1,1)
xdata.Slice3(0,0,r1,0,sx,1,1,sy,1,2,r2,1) = axis.Slice3(0,0,ch2,0,sx,1,1,sy,1,2,r2,1)
image ydata := RealImage("ydata",4,sx,sy, nsz)
ydata.Slice3(0,0,0,0,sx,1,1,sy,1,2,r1,1) = data.Slice3(0,0,ch1,0,sx,1,1,sy,1,2,r1,1)
ydata.Slice3(0,0,r1,0,sx,1,1,sy,1,2,r2,1) = data.Slice3(0,0,ch2,0,sx,1,1,sy,1,2,r2,1)
image xmean := RealImage("xmean",4,sx,sy, nsz)
image ymean := RealImage("ymean",4,sx,sy, nsz)
for (number i=0; i<nsz; i++) {
// to repeat the mean value as a 3d stack
// the sum value of each pixel can be rapidly done by the Project function
xmean.Slice2(0,0,i,0,sx,1,1,sy,1) = Project(xdata,2) / nsz
ymean.Slice2(0,0,i,0,sx,1,1,sy,1) = Project(ydata,2) / nsz
};
// do regression
image m = Project((xdata-xmean)*(ydata-ymean), 2) / Project((xdata-xmean)**2, 2)
image b = ymean[icol,irow,0] - m*xmean[icol,irow,0]
// to apply y=mx+b to all the planes
image bkg := src.ImageClone()
for (number i=0; i<sz; i++) {
bkg.Slice2(0,0,i,0,sx,1,1,sy,1) = m*axis.Slice2(0,0,i,0,sx,1,1,sy,1)+b
};
// reconstruct to the power law bkg
bkg = exp(bkg)
return bkg
};
}
{
Object FTW = Alloc(FittingTwoWindowBkg)
number Eng1 = 1500 // eV, 1st Energy of fitting window
number Rng1 = 200 // eV, fitting range of 1st window
number Eng2 = 2600 // eV, 2nd Energy of fitting window
number Rng2 = 200 // eV, fitting range of 2nd window
image src := GetFrontImage()
image bkg := FTW.TwoWindowBkgFitting(src, Eng1, Rng1, Eng2, Rng2)
bkg.ShowImage()
};