macro 'calculate pixelwise standard deviation of 16 bit stack';
var
    n, i, width, wpix, height: integer;
    inPid, sumXSqdPid, sqdSumXPid, XPid, StdDevPid, outPid: integer;
begin
    inPid := pidNumber;
    n := nSlices;
    getPicSize(width, height);
    wpix := width div 2;
    saveState;
    setBackground(0);
    setForeground(255);
    setNewSize(wpix * 4, height);
    makeNewWindow('sum X sqd');
    sumXSqdPid := pidNumber;
    makeNewWindow('sqd sum X');
    sqdSumXPid := pidNumber;
    makeNewWindow('X');
    XPid := pidNumber;
    makeNewWindow('StdDev real');
    StdDevPid := pidNumber;
    setNewSize(wpix * 2, height);
    makeNewWindow('Stdev 16 bit');
    outPid := pidNumber;
    restoreState;
    for i := 1 to n do begin
        choosePic(inPid);
        chooseSlice(i);
        cnvrt16uT32r(inPid,XPid);
        add32r(XPid,sqdSumXPid,sqdSumXPid);
        mpy32r(XPid,Xpid,XPid);
        add32r(Xpid,sumXSqdPid,sumXSqdPid);
    end;
    mpyK32r(n,sumXSqdPid,sumXSqdPid);
    mpy32r(sqdSumXPid,sqdSumXPid,sqdSumXPid);
    sub32r(sumXSqdPid,sqdSumXPid,StdDevPid);
    mpyK32r(1/n/(n-1),StdDevPid,StdDevPid);
    sqrt32r(StdDevPid,StdDevPid);
    cnvrt32rT16u(StdDevPid, outPid);
end;

macro '[m]make stack';
var
    i: integer;
    inPid: integer;
    v: real;
begin
    setNewSize(64,32);
    makeNewStack('input');
    inPid := pidNumber;
    makeRoi(2,1,60,30);
    setCounter(10);
    for i := 1 to 10 do begin
        if i <> 1 then addSlice;
        v := 100000*exp(-i)+1000;
        rY[i] := v;
        fill16u(v,inPid);
    end;
end;

macro '[d]dispose all';
begin
disposeAll;
end;

{procedure LeastSqr32r(sumX,n,sumXsqd: real; sumY,sumXY,m,b: pidNum);}
macro 'least squares fit to 16 bit stack';
var
    X, sumX, n, sumXsqd: real;
    YPid, sumYPid, sumXYPid, mPid, bPid: integer;
    inPid, outPid, n, i, width, wpix, height: integer;
begin
    inPid := pidNumber;
    n := nSlices;
    for i := 1 to n do rX[i] := i*i;
    getPicSize(width, height);
    wpix := width div 2;
    makeRoi(2,1,wpix*2-4,height-2);
    saveState;
    setBackground(0);
    setForeground(255);
    setNewSize(wpix * 2, height);
    makeNewStack('out');
    outPid := pidNumber;
    makeRoi(2,1,wpix*2-4,height-2);
    for i := 2 to n do addSlice;
    setNewSize(wpix * 4, height);
    makeNewWindow('Y');
    YPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2); {outside of ROI should stay zero}
    makeNewWindow('sum Y');
    sumYPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('sum XY');
    sumXYPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('m (slope)');
    mPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('b (Y intercept)');
    bPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    restoreState;
    sumX := 0.0;
    sumXsqd := 0.0;
    for i := 1 to n do begin
        X := rX[i]; {X coordinate is in the measurements array}
        sumX := sumX + X;
        sumXsqd := sumXsqd + X * X;
        choosePic(inPid);
        chooseSlice(i);
        cnvrt16uT32r(inPid,YPid);
        add32r(YPid,SumYPid,SumYPid);
        mpyK32r(X,YPid,YPid);
        add32r(YPid,sumXYPid,sumXYPid);
    end;
    LeastSqr32r(sumX, n, sumXsqd, sumYPid, sumXYPid, mPid, bPid);
    for i := 1 to n do begin
        X := rX[i];
        mpyK32r(X,mPid,YPid);
        add32r(bPid,YPid,YPid);
        choosePic(outPid);
        chooseSlice(i);
        cnvrt32rT16u(YPid,outPid);
    end;
end;

{ExpFtStep32r(BestResid,NextResid,BestBase,NextBase,DeltaBase: pidNum);}
macro 'exponential fit';
var
    X, sumX, n, sumXsqd: real;
    YPid, newYPid, sumYPid, sumXYPid, mPid, bPid: integer;
    inPid, outPid, i, width, wpix, height: integer;
    BestResidPid,NextResidPid,BestBasePid,NextBasePid,DeltaBasePid: integer;
    iter, last: integer;
begin
    inPid := pidNumber;
    n := nSlices;
    for i := 1 to n do rX[i] := i;
    getPicSize(width, height);
    wpix := width div 2;
    makeRoi(2,1,wpix*2-4,height-2);
    saveState;
    setBackground(0);
    setForeground(255);
    setNewSize(wpix * 2, height);
    makeNewStack('out');
    outPid := pidNumber;
    makeRoi(2,1,wpix*2-4,height-2);
    for i := 2 to n do addSlice;
    setNewSize(wpix * 4, height);
    makeNewWindow('Y');
    YPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2); {outside of ROI should stay zero}
    makeNewWindow('new Y = exp(m * x + b) + base');
    newYPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2); {outside of ROI should stay zero}
    makeNewWindow('sum Y');
    sumYPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('sum XY');
    sumXYPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('m (slope)');
    mPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('b (Y intercept)');
    bPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('Best Residuals');
    BestResidPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('Next Residuals');
    NextResidPid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('Best Baseline');
    BestBasePid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('Next Baseline');
    NextBasePid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    makeNewWindow('Delta Baseline');
    DeltaBasePid := pidNumber;
    makeRoi(4,1,wpix*4-8,height-2);
    restoreState;
    fill32r(exp(ln(10)*20),BestResidPid);
    fill32r(0,BestBasePid);
    fill32r(0,NextBasePid);
    fill32r(exp(ln(10)*20),DeltaBasePid);
  last := 15;
  for iter := 1 to last do begin
    showMessage('iter=',iter);
    sumX := 0.0;
    sumXsqd := 0.0;
    fill32r(0,SumYPid);
    fill32r(0,SumXYPid);
    { add32r(BestBasePid,DeltaBasePid,NextBasePid);--already done by expftstep}
    if iter = last then begin
        addK32r(0,BestBasePid, NextBasePid);
        {this causes the last iteration to recalculate the best parameters}
    end;
    for i := 1 to n do begin
        showMessage('iter=',iter,'\i=',i);
        X := rX[i]; {X coordinate is in the measurements array}
        sumX := sumX + X;
        sumXsqd := sumXsqd + X * X;
        choosePic(inPid);
        chooseSlice(i);
        cnvrt16uT32r(inPid,YPid);
        if iter = 1 then begin
            sml32r(YPid,DeltaBasePid,DeltaBasePid);
        end;
        sub32r(YPid,NextBasePid,YPid);
        lrgK32r(0.001,YPid,YPid);
        rUser2[i] := GetPixel32r(YPid,10,10);
        ln32r(YPid,YPid);
        rUser1[i] := GetPixel32r(YPid,10,10);
        add32r(YPid,SumYPid,SumYPid);
        mpyK32r(X,YPid,YPid);
        add32r(YPid,SumXYPid,SumXYPid);
    end;
    LeastSqr32r(SumX, n, SumXsqd, SumYPid, SumXYPid, mPid, bPid);
    fill32r(0.0,nextResidPid);
    for i := 1 to n do begin
        showMessage('iter=',iter,'\resid i=',i);
        X := rX[i];
        choosePic(inPid);
        chooseSlice(i);
        cnvrt16uT32r(inPid,YPid);
        mpyK32r(X,mPid,newYPid);
        add32r(bPid,newYPid,newYPid);
        exp32r(newYPid,newYPid);
        add32r(newYPid,nextBasePid,newYPid);
        if iter = last then begin
            choosePic(outPid);
            chooseSlice(i);
            cnvrt32rT16u(newYPid,outPid);
        end;
        sub32r(newYPid,YPid,newYPid);
        mpy32r(newYPid,newYPid,newYPid);
        add32r(newYPid,nextResidPid,nextResidPid);
    end;
    if iter = 1 then begin
        addK32r(0,NextResidPid,BestResidPid);
        addK32r(0,NextBasePid,BestBasePid);
        mpyK32r(0.5,DeltaBasePid,DeltaBasePid);
        add32r(BestBasePid,DeltaBasePid,NextBasePid);
    end else if iter <> last then begin
        {ExpFtStep32r is:  for all pixels do begin
             if nextResid < bestResid then begin
                   bestResid := nextResid;
                   bestBase := nextBase;
             end else begin
                   deltaBase := -1/2*deltaBase;
             end;
             nextBase := bestBase + deltaBase;
         end}
        ExpFtStep32r(
            BestResidPid, NextResidPid,
            BestBasePid, NextBasePid, DeltaBasePid);
    end;
  end;
end;

macro '[p]get real pixel from front image';
var x,y: integer;
begin
   x := 10;
   y := 10;
{   x := getnumber('x',x);
   y := getnumber('y',y);}
   showmessage('pid',pidNumber,'x',x,'y',y,'getpixel32r',getpixel32r(pidNumber,x,x):1:4);
end;

macro 'delete slices 1-30';
var i: integer;
begin
   for i := 1 to 30 do deleteslice;
end;

macro 'delete slices to end';
var i: integer;
begin
    i := sliceNumber;
   while i <= nSlices do begin
     selectSlice(i);
     deleteslice;
    end;
end;

macro 'convert front real image to 8 bits and report scale';
var amin,amax: real;
    dup, out, width, height: integer;
begin
   duplicate('deleteme');
   dup := pidNumber;
   minmax32r(dup,amin,amax);
   subk32r(amin,dup,dup);
   mpyK32r(256/(amax-amin),dup,dup);
   getPicSize(width,height);
   width := width div 4;
   saveState;
   setnewsize(width, height);
    makenewwindow('scaled ',amin:0:4,' to ',amax:0:4);
    out := pidNumber;
   restoreState;
    cnvrt32rTo8(dup,out);
    choosePic(dup);
    dispose;
    selectPic(out);
end;