【问题标题】：Are there any description about using "fit image" function (DM FitTools) in script?是否有关于在脚本中使用“适合图像”功能（DM FitTools）的描述？
【发布时间】：2016-04-28 09:24:21
【问题描述】：

我想将 DM 中“Fit Image” plattlet 提供的功能（尤其是 2D 多项式拟合到给定的输入图像并减去它）集成到脚本中，以自动化整个图像处理流程。但是，我找不到任何关于如何做到这一点的描述。

如果有人知道它，或者对此有某些文档，我们将不胜感激。

【问题讨论】：

标签： 2d curve-fitting data-fitting dm-script

【解决方案1】：

尚未正式支持/记录用于拟合的脚本功能。

但是，您可以使用以下示例来了解命令的工作原理：

命令

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，二维高斯拟合

// $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，二维多项式拟合

// $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)

【讨论】：