我想将 DM 中“Fit Image” plattlet 提供的函数（特别是 2D 多项式拟合到给定的输入图像并减去它）集成到脚本中，以自动化整个图像处理流程。但是，我找不到任何关于如何做到这一点的描述。
如果有人知道它，或者对此有一定的文档，我们将不胜感激。
问问题
269 次
1 回答
3
尚未正式支持/记录用于拟合的脚本功能。
但是，您可以使用以下示例来查看命令的工作原理：
命令
Boolean FitGaussian(Image* data, Image* errors, double* N, double* mu, double* sigma, double* chiSqr, double conv_cond)
ImageRef PlotGaussian(Image* data, double N, double mu, double sigma)
Boolean FitLorentzian(Image* data, Image* errors, double* I, double* x0, double* gamma, double* chiSqr, double conv_cond)
ImageRef PlotLorentzian(Image* data, double I, double x0, double gamma)
Boolean FitPolynomial(Image* data, Image* errors, Image* pars, Image* parsToFit, double* chiSqr, double conv_cond)
ImageRef PlotPolynomial(Image* data, Image* pars)
Boolean FitGaussian2D(Image* data, Image* errors, Image* pars, Image* parsToFit, double* chiSqr, double conv_cond)
ImageRef PlotGaussian2D(Image* data, Image* pars)
Boolean FitPolynomial2D(Image* data, Image* errors, Image* pars, Image* parsToFit, double* chiSqr, double conv_cond)
ImageRef PlotPolynomial2D(Image* data, Image* pars)
Boolean FitFormula(dm_string formulaStr, Image* data, Image* errors, Image* pars, Image* parsToFit, double* chiSqr, double conv_cond)
ImageRef PlotFormula(dm_string formulaStr, Image* data, Image* pars)
示例 1，一维公式拟合
// create the input image:
Image input := NewImage("formula test", 2, 100)
input = 500.5 - icol*11.1 + icol*icol*0.11
// add some random noise:
input += (random()-0.5)*sqrt(abs(input))
// create image with error data (not required)
Image errors := input.ImageClone()
errors = tert(input > 1, sqrt(input), 1)
// setup fit:
Image pars := NewImage("pars", 2, 3)
Image parsToFit := NewImage("pars to fit", 2, 3)
pars = 10; // starting values
parsToFit = 1;
Number chiSqr = 1e6
Number conv_cond = 0.00001
Result("\n starting pars = {")
Number xSize = pars.ImageGetDimensionSize(0)
Number i = 0
for (i = 0; i < xSize; i++)
{
Result(GetPixel(pars, i, 0))
if (i < (xSize-1)) Result(", ")
}
Result("}")
// fit:
String formulaStr = "p0 + p1*x + p2*x**2"
Number ok = FitFormula(formulaStr, input, errors, pars, parsToFit, chiSqr, conv_cond)
Result("\n results pars = {")
for (i = 0; i < xSize; i++)
{
Result(GetPixel(pars, i, 0))
if (i < (xSize-1)) Result(", ")
}
Result("}")
Result(", chiSqr ="+ chiSqr)
// plot results of fit:
Image plot := PlotFormula(formulaStr, input, pars)
// compare the plot and original data:
Image compare := NewImage("Compare Fit", 2, 100, 3)
compare[icol, 0] = input // original data
compare[icol, 1] = plot // fit function
compare[icol, 2] = input - plot // residuals
ImageDocument linePlotDoc = CreateImageDocument("Test Fitting")
ImageDisplay linePlotDsp = linePlotDoc.ImageDocumentAddImageDisplay(compare, 3)
linePlotDoc.ImageDocumentShow()
示例 2，2D 高斯拟合
// $BACKGROUND$
// create data image
Image img := NewImage("Gaussian2D", 2, 200, 200)
Image true_pars := NewImage("Gaussian2D Pars", 2, 6)
// true parameters
true_pars[0,0] = 1000 // height of gaussian
true_pars[1,0] = 60 // center in x
true_pars[2,0] = 50 // width in x
true_pars[3,0] = 40 // center in y
true_pars[4,0] = 80 // width in y
true_pars[5,0] = 0.7 // rotation in radians
Image data := PlotGaussian2D(img, true_pars)
data += (gaussianrandom())*sqrt(abs(data)) //add noise
ShowImage(data)
Image errors := data.ImageClone()
errors = tert(abs(data) > 1, sqrt(abs(data)), 1)
// starting parameters of fit
Image pars := NewImage("Gaussian2D Pars", 2, 6)
pars = 100
pars[0,0] = max(data) // estimate normalization from peak of data
pars[5,0] = 0 // 100 radians doesn't make sense
Image parsToFit := NewImage("tmp", 2, 6)
parsToFit = 1
Number chiSqr = 1e6
Number conv_cond = 0.00001
Result("\n starting pars = {")
Number xSize = pars.ImageGetDimensionSize(0)
Number i = 0
for (i = 0; i < xSize; i++)
{
Result(GetPixel(pars, i, 0))
if (i < (xSize-1)) Result(", ")
}
Result("}")
// fit
Number ok = FitGaussian2D(data, errors, pars, parsToFit, chiSqr, conv_cond)
if (chiSqr > 2)
ok = FitGaussian2D(data, errors, pars, parsToFit, chiSqr, conv_cond)
Image parDif = 100.0*(pars - true_pars)/true_pars
Result("\n results pars (% dif from true)= {")
for (i = 0; i < xSize; i++)
{
Result(GetPixel(parDif, i, 0))
if (i < (xSize-1)) Result(", ")
}
Result("}")
Result(", chiSqr ="+ chiSqr)
// show residuals
Image residuals := PlotGaussian2D(img, pars)
residuals = data - residuals
ShowImage(residuals)
示例 3，2D 多项式拟合
// $BACKGROUND$
// The number of parameters are defined by the order,
// nPar = (order+1)*(order+2)/2. For example, a
// 3rd order poly will have (3+1)*(3+2)/2 = 10 parameters:
//
// x^0 x^1 x^2 x^3
// --------------------------------
// y^0 | p0 p1 p2 p3
// y^1 | p4 p5 p6 -- (the -- terms are higher)
// y^2 | p7 p8 -- -- (order so are ignored )
// y^3 | p9 -- -- --
//
//
// i.e. f(x,y|p) = p0 + p1*x + p2*x^2 + p3*x^3 + p4*y + p5*x*y
// + p6*x^2*y + p7*y^2 + p8*x*y^2 + p9*y^3
Number xImgSize = 512
Number yImgSize = 512 // create data image
Image img := NewImage("Poly2D", 2, xImgSize, yImgSize)
Image pars_true := NewImage("Poly2D Pars", 2, 3, 3)
// true parameters
pars_true[0,0] = 100
pars_true[1,0] = 60
pars_true[2,0] = -0.05
pars_true[0,1] = 70
pars_true[1,1] = 0.01
pars_true[0,2] = -0.1
Image data := PlotPolynomial2D(img, pars_true)
data += (gaussianrandom())*sqrt(abs(data)) //add noise
ShowImage(data)
Image errors := data.ImageClone()
errors = tert(abs(data) > 1, sqrt(abs(data)), 1)
// starting parameters of fit
Image pars := NewImage("Poly2D Pars", 2, 3, 3)
pars = 10
Image parsToFit := NewImage("tmp", 2, 3, 3)
parsToFit = 1
Number chiSqr = 1e6
Number conv_cond = 0.00001
Result("\n starting pars = {")
Number xSize = pars.ImageGetDimensionSize(0)
Number ySize = pars.ImageGetDimensionSize(1)
Number i, j
for (j = 0; j < ySize; j++)
{
if (j > 0) Result(", ")
Result("{")
for (i = 0; i < xSize; i++)
{
if (i > 0) Result(", ")
if ((i+j) > 2)
Result("-")
else
Result(GetPixel(pars, i, j))
}
Result("}")
}
Result("}")
// fit
Number startTicks = GetHighResTickCount()
Number ok = FitPolynomial2D(data, errors, pars, parsToFit, chiSqr, conv_cond)
Number endTicks = GetHighResTickCount()
Number secs = CalcHighResSecondsBetween(startTicks, endTicks)
Image parDif = 100*(pars - pars_true)/pars_true
Result("\n results pars (% diff from true) = {")
for (j = 0; j < ySize; j++)
{
if (j > 0) Result(", ")
Result("{")
for (i = 0; i < xSize; i++)
{
if (i > 0) Result(", ")
if ((i+j) > 2)
Result("-")
else
Result(GetPixel(parDif, i, j))
}
Result("}")
}
Result("}")
Result(", chiSqr = "+ chiSqr)
Result(", Fit Time (s) = " + secs)
// show residuals
Image residuals := PlotPolynomial2D(img, pars)
residuals = data - residuals
ShowImage(residuals)
于 2016-04-28T11:42:59.973 回答