{Prepare irregularly Cropped 16 bit Stack}
{Locate folder containing the raw images}
{Open mark image named 'CropMask'}
{Create a stack of the proper size.}
{Open images in sequence and squeeze into stack.}

{sqzmask(in,mask,out)}

var
 vectPid, mrkvPid, maskPid, markPid, newPid, playPid, x16Pid, b16Pid: Integer;
 xlate16uto8Pid, ymin, ymax: integer;
 kx, ky, kw, kh: integer;
 kernelPid: integer;

macro 'Create global variable save image';
begin
  setNewSize(32,32);
  MakeNewWindow('Global Variables');
end;

macro 'Save global variables';
var
  savePid, globalPid: integer;
begin
  savePid := pidNumber;
  SelectWindow('Global Variables');
  globalPid := pidNumber;
  PutPixVec32s(globalPid,1,vectPid);
  PutPixVec32s(globalPid,2,mrkvPid);
  PutPixVec32s(globalPid,3,maskPid);
  PutPixVec32s(globalPid,4,markPid);
  PutPixVec32s(globalPid,5,newPid);
  PutPixVec32s(globalPid,6,playPid);
  PutPixVec32s(globalPid,7,x16Pid);
  PutPixVec32s(globalPid,8,xlate16uto8Pid);
  PutPixVec32s(globalPid,9,ymin);
  PutPixVec32s(globalPid,10,ymax);
  PutPixVec32s(globalPid,11,kx);
  PutPixVec32s(globalPid,12,ky);
  PutPixVec32s(globalPid,13,kw);
  PutPixVec32s(globalPid,14,kh);
  PutPixVec32s(globalPid,15,kernelPid);
  PutPixVec32s(globalPid,16,b16Pid);
  SelectPic(savePid);
end;

macro 'Restore global variables';
var
  savePid, globalPid: integer;
begin
  savePid := pidNumber;
  SelectWindow('Global Variables');
  globalPid := pidNumber;
  vectPid := GetPixVec32s(globalPid,1);
  mrkvPid := GetPixVec32s(globalPid,2);
  maskPid := GetPixVec32s(globalPid,3);
  markPid := GetPixVec32s(globalPid,4);
  newPid := GetPixVec32s(globalPid,5);
  playPid := GetPixVec32s(globalPid,6);
  x16Pid := GetPixVec32s(globalPid,7);
  xlate16uto8Pid := GetPixVec32s(globalPid,8);
  ymin := GetPixVec32s(globalPid,9);
  ymax := GetPixVec32s(globalPid,10);
  kx := GetPixVec32s(globalPid,11);
  ky := GetPixVec32s(globalPid,12);
  kw := GetPixVec32s(globalPid,13);
  kh := GetPixVec32s(globalPid,14);
  kernelPid := GetPixVec32s(globalPid,15);
  b16Pid := GetPixVec32s(globalPid,16);
  SelectPic(savePid);
end;

macro 'Front image is mask image';
begin
 maskPid := pidNumber;
end;

macro 'Front image is mark image';
begin
 markPid := pidNumber;
end;

macro 'Front image is 16 bit cropped stack';
begin
 vectPid := pidNumber;
end;

macro 'Front image is 16 to 8 conversion table';
begin
 xlate16uto8Pid := pidNumber;
end;

macro 'Front image is 8 bit display buffer';
begin
 playPid := pidNumber;
end;

macro 'Front image is 16 bit slice buffer';
begin
 x16Pid := pidNumber;
end;

macro 'Front image is kernel';
begin
 kernelPid := PidNumber;
 kx := getPixel(0, 0);
 ky := getPixel(0, 1);
 kw := getPixel(0, 2);
 kh := getPixel(0, 3);
end;

macro 'Select data image';
begin
  SelectPic(MarkupPicNumber('data'));
end;

macro 'Select data16 image';
begin
  SelectPic(MarkupPicNumber('data16'));
end;

macro 'Select marks image';
begin
  SelectPic(MarkupPicNumber('mark'));
end;

macro 'Select bkg image';
begin
  SelectPic(MarkupPicNumber('bkg'));
end;

procedure checkPids;
 var
  left, top, width, height: integer;
begin
 RequiresUser('Pixel16u', 1);
 if not pidExists(maskPid) then begin
   PutMessage('Please run macro "Front image is mask image"');
   Exit;
  end;
 if not pidExists(vectPid) then begin
   PutMessage('Please run macro "Front image is 16 bit cropped stack"');
   Exit;
  end;
 if not pidExists(xlate16uto8Pid) then begin
   PutMessage('Please run macro "Front image is 16 to 8 conversion table"');
   Exit;
  end;
 if not pidExists(playPid) then begin
   selectPic(maskPid);
   GetPicSize(width, height);
   SetNewSize(width, height);
   MakeNewWindow('Play 16 bit movie');
   playPid := pidNumber;
  end;
 if not pidExists(x16Pid) then begin
   selectPic(maskPid);
   GetPicSize(width, height);
   SetNewSize(2 * width, height);
   MakeNewWindow('16 bit data buffer');
   x16Pid := pidNumber;
  end;
 if not pidExists(b16Pid) then begin
   selectPic(maskPid);
   GetPicSize(width, height);
   SetNewSize(2 * width, height);
   MakeNewWindow('16 bit background');
   b16Pid := pidNumber;
  end;
 if not pidExists(markPid) then begin
   selectPic(maskPid);
   GetPicSize(width, height);
   SetNewSize(width, height);
   MakeNewWindow('mark');
   markPid := pidNumber;
  end;
end;

macro 'Show markup';
begin
  checkPids;
  ChoosePic(playPid);
  MarkupData8(picNumber);
  ChoosePic(x16Pid);
  MarkupData16(picNumber);
  ChoosePic(b16Pid);
  MarkupBkg(picNumber);
  ChoosePic(markPid);
  MarkupMarks(picNumber);
  MarkupShow;
end;

procedure setKernelSize;
begin
 kx := trunc(GetNumber('kernel center x', kx));
 ky := trunc(GetNumber('kernel center y', ky));
 kw := trunc(GetNumber('kernel width', kw));
 kh := trunc(GetNumber('kernel height', kh));
end;

macro 'set kernel size';
begin
 setKernelSize;
end;

procedure createKernel;
var
 x, y: integer;
 frontImage: integer;
 sum: real;
begin
 RequiresUser('getputpixel', 1);
 frontImage := pidNumber;
 if frontImage = kernelPid then
  frontImage := 0;
 if kernelPid <> 0 then begin
   x := kernelPid;
   kernelPid := 0;
   SelectPic(x);
   Dispose;
  end;
 if (kx = 0) or (ky = 0) or (kw = 0) or (kh = 0) then begin
   kx := 2;
   ky := 2;
   kw := 5;
   kh := 5;
  end;
 SaveState;
 SetNewSize(kw * 4, kh + 1);
 MakeNewWindow('kernel ', kx : 0, 'x ', ky : 0, 'y ', kw : 0, 'w ', kh : 0, 'h');
 kernelPid := pidNumber;
 RestoreState;
 putPixel(0, 0, kx);
 putPixel(1, 0, ky);
 putPixel(2, 0, kw);
 putPixel(3, 0, kh);
 sum := 0.0;
 for x := -kx to kw - kx - 1 do begin
   for y := -ky to kh - ky - 1 do begin
     sum := sum + exp(-sqrt(sqr(x) + sqr(y)));
    end;
  end;
 sum := 32000.0 / sum; {nearly maximum before overflow on 65535 pixel}
 for x := -kx to kw - kx - 1 do begin
   for y := -ky to kh - ky - 1 do begin
     putPixel32s(kernelPid, x + kx, y + ky + 1, sum * exp(-sqrt(sqr(x) + sqr(y))));
    end;
  end;
 sum := 0.0;
 for x := -kx to kw - kx - 1 do begin
   for y := -ky to kh - ky - 1 do begin
     sum := sum + getPixel32s(kernelPid, x + kx, y + ky + 1);
    end;
  end;
 showmessage('kernel sum = ', sum);
 SelectPic(kernelPid);
 MakeRoi(0, 1, kw * 4, kh);
 if frontImage <> 0 then
  SelectPic(frontImage);
end;

macro 'create Kernel';
begin
  createKernel;
end;

macro '[c] Make Irregularly Cropped 16 bit stack';
var
 i, n: integer;
 left, top, width, height: integer;
 origPid: integer;
begin
  {Find the folder}
 GetFolder;
  {Open the CropMask file}
 Open('CropMask');
  {Make sure name of front image is CropMask???}

 maskPid := pidNumber;
 ChoosePic(maskPid);

 Measure;
 GetROI(left, top, width, height);
 if width = 0 then
  GetPicSize(width, height);
 n := width * height - histogram[0];
 width := Sqrt(n) + 1;
 if width < 32 then
  width := 32;
 height := n div width + 1;
 SetNewSize(width * 2, height); {16 bits}
 MakeNewStack('packed Crop16');
 vectPid := pidNumber;
 SetNewSize(width,height);
 MakeNewStack('packed mark stack');
 mrkvPid := pidNumber;
 n := GetNumber('Number of frames', 1);
 for i := 1 to n do begin
   open('Frame', i : 4); {??? need import???}
   origPid := pidNumber;
   choosePic(vectPid);
   if i > nSlices then
    addSlice;
   chooseSlice(i);
   choosePic(mrkvPid);
   if i > nSlices then
     addSlice;
   chooseSlice(i);
   SqzMask16u(origPid, maskPid, vectPid);
   choosePic(origPid);
   dispose;
  end;
 selectPic(vectPid);
 selectSlice(1);
end;

macro '[m] play cropped 16 bit movie';
var
 i, n, fps: integer;
begin
 CheckPids;
 selectPic(playPid);
 choosePic(vectPid);
 n := nSlices;
 KillDelay(1);
 fps := GetNumber('Frames per second', 20);
 while not button do begin
   StartElapsed(1);
   for i := 1 to n do begin
     StartDelay(1, 1 / fps);
     choosePic(vectPid);
     chooseSlice(i);
     choosePic(mrkvPid);
     chooseSlice(i);
     CnvMsk16uTo8(vectPid, xlate16uto8Pid, maskPid, playPid);
     {ExpMask8u(mrkvPid,maskPid,markPid);}
     selectPic(playPid);
     {selectPic(markPid);--way too slow}
     WaitDelay(1);
    end;
   ShowMessage(concat(n / MeasureElapsed(1) : 0 : 2, ' frames per second\hold mouse button to stop.'));
  end;
  ExpMask8u(mrkvPid,maskPid,markPid);
  SelectPic(markPid);
end;

macro 'Copy mark changes into mark stack';
var
  n: integer;
begin
 CheckPids;
 choosePic(vectPid);
 n := sliceNumber;
 choosePic(mrkvPid);
 SelectSlice(n);
 SqzMask8u(markPid,maskPid,mrkvPid);
end;

macro '[n] next slice';
var
  n: integer;
begin
 CheckPids;
 choosePic(vectPid);
 if sliceNumber < nSlices then
  SelectSlice(sliceNumber + 1);
 n := sliceNumber;
 choosePic(mrkvPid);
 SelectSlice(n);
 CnvMsk16uTo8(vectPid, xlate16uto8Pid, maskPid, playPid);
 ExpMask16u(vectPid, maskPid, x16Pid);
 ExpMask8u(mrkvPid,maskPid,markPid);
 choosePic(markPid);
 SelectAll;
 selectPic(playPid);
end;

macro '[b] back to previous slice';
var
  n: integer;
begin
 CheckPids;
 choosePic(vectPid);
 if sliceNumber > 1 then
  selectSlice(sliceNumber - 1);
 n := sliceNumber;
 choosePic(mrkvPid);
 SelectSlice(n);
 CnvMsk16uTo8(vectPid, xlate16uto8Pid, maskPid, playPid);
 ExpMask16u(vectPid, maskPid, x16Pid);
 ExpMask8u(mrkvPid,maskPid,markPid);
 choosePic(markPid);
 SelectAll;
 selectPic(playPid);
end;

macro 'convert front 8 bit image to 16 bits';
var
 i: integer;
 inPid, xlatePid, outPid: integer;
 width, height: integer;
begin
 SaveState;
 inPid := pidNumber;
 GetPicSize(width, height);
 SetNewSize(512, 1);
 SetBackgroundColor(0);
 MakeNewWindow('8 to 16 conversion table');
 xlatePid := pidNumber;
 for i := 0 to 255 do
  putpixel(i * 2, 0, i);
 MakeRoi(0, 0, 512, 1);
 SetNewSize(width * 2, height);
 MakeNewWindow('16 bit image from ', inPid : 0);
 outPid := pidNumber;
 RestoreState;
 Cnvrt8to16u(inPid, xlatePid, outPid);
end;

macro 'set 16 to 8 conversion parameters';
begin
 ymin := GetNumber('minimum 8 bit pixel value', 1);
 ymax := GetNumber('maximum 8 bit pixel value', 254);
end;

macro 'make 16 to 8 conversion table';
var
 i, j, x, y: integer;
 width, height: integer;
 inPid, min, max: integer;
begin
 RequiresUser('pixel16u', 1);
 if ymax = 0 then begin
   ymin := GetNumber('minimum 8 bit pixel value', 1);
   ymax := GetNumber('maximum 8 bit pixel value', 254);
  end;
 SaveState;
 inPid := pidNumber;
 minMax16u(inPid, min, max);
 SetNewSize(256, 256);
 SetBackgroundColor(0);
 MakeNewWindow('16 to 8 conversion table');
 xlate16uto8Pid := pidNumber;
 linLUT16uto8(xlate16uto8Pid, min, max, ymin, ymax);
end;

macro 'convert front 16 bit image to 8 bits';
var
 i, j, x, y: integer;
 inPid, outPid: integer;
 width, height: integer;
 min, max: integer;
begin
 RequiresUser('pixel16u', 1);
 inPid := pidNumber;
 GetPicSize(width, height);
 SaveState;
 ChoosePic(xlate16uto8Pid);
 MakeRoi(0, 0, 256, 256);
 SetNewSize(width div 2, height);
 MakeNewWindow('8 bit image from ', inPid : 0);
 outPid := pidNumber;
 RestoreState;
 Cnvrt16uto8(inPid, xlate16uto8Pid, outPid);
end;

{mark image is from markup, indicates segmentation}
{mask image:  zero pixels in mask are invalid and cannot }
{contribute to the background calculation. (no data there)}
{The  "mask" image used with CnvMsk16uTo8}
{is opposite polarity from mask used with Convolve16u,}
{zero means no data for CnvMsk,}
{nonzero means no data for Convolve.}

macro 'calculate integrated intensities';
var
 flagPid, outPid, tmskPid : integer;
 width, height, wpix, hpix: integer;
begin
  SaveState;
 CheckPids;
 if not pidExists(kernelPid) then
  createKernel;
 ChoosePic(kernelPid);
 kx := getPixel(0, 0);
 ky := getPixel(1, 0);
 kw := getPixel(2, 0);
 kh := getPixel(3, 0);
 MakeRoi(0,1,kw*4,kh);
 ChoosePic(x16Pid);
 KillRoi;
 GetPicSize(width, height);
 SetNewSize(width, height);
 wpix := width div 2;
 hpix := height;
 MakeNewWindow('output image for ', x16Pid : 0);
 outPid := pidNumber;
 ChoosePic(markPid);
 KillRoi;
 GetPicSize(width, height);
  if wpix <> width then begin
    putmessage('mark width is incorrect');
    Exit;
  end;
  if hpix <> height then begin
    putmessage('mark height is incorrect');
    Exit;
  end;
 ChoosePic(maskPid);
 KillRoi;
 GetPicSize(width, height);
  if wpix <> width then begin
    putmessage('mask width is incorrect');
    Exit;
  end;
  if hpix <> height then begin
    putmessage('mask height is incorrect');
    Exit;
  end;
  SetNewSize(wpix, hpix);
  MakeNewWindow('flag image for ', x16Pid : 0);
  flagPid := PidNumber;
  SetNewSize(wpix+kw-1, hpix+kh-1);
  MakeNewWindow('temp mask image for ', x16Pid : 0);
  tMskPid := PidNumber;
  RestoreState;


  ChoosePic(markPid);
  SelectAll;
  Copy;


  {Flag image should be zero only in marked pixels}
  {These are the pixels that get an output value calculated.}
  {It is a waste of time to calculate values for other pixels.}
  ChoosePic(flagPid);
  KillRoi;
  Paste;
  KillRoi;
  ChangeValues(1,255,255);
  Invert;

  {Mask image indicates which pixels do not contribute }
  {to background.  This includes foreground pixels, and }
  {pixels which have no value (outside image).  These}
  {pixels have nonzero mask values.  Zero mask pixels}
  {will contribute to the background calculation.}

  ChoosePic(maskPid);
  SelectAll;
  Copy;
  KillRoi;
  ChoosePic(tMskPid);
  KillRoi;
  Paste;
  KillRoi;
  changeValues(1,255,255);
  Invert;
  ChoosePic(markPid);
  SelectAll;
  Copy;
  KillRoi;
  ChoosePic(tMskPid);
  KillRoi;
  SetForegroundColor(255);
  SetBackgroundColor(0);
  Paste;
  doOr; {doesn't work if foreground not 255}

  {Pixels outside of the input image must be nonzero}
  {in the mask image.  The pixels in question should }
  {already nonzero but change them just in case.}
  SetForegroundColor(255);
  MakeRoi(0,0,kx,hpix+kh);
  Fill;
  MakeRoi(kx+wpix,0,kw-kx-1,hpix+kh)
  Fill;
  MakeRoi(kx,0,wpix,ky);
  Fill;
  MakeRoi(kx,ky+hpix,wpix,kh-ky-1);
  Fill;
  {Mask image must have an ROI same size as image and}
  {with borders matching kernel, thus:}
  MakeRoi(kx,ky,wpix,hpix);

  {Finally, do the convolution to calculate the background}

  Convolve16u(flagPid,x16Pid,kernelPid,kx,ky,tMskPid,outPid); 
  {alphabetical order}

  ChoosePic(flagPid);
  {Dispose;}
  ChoosePic(tMskPid);
  {Dispose;}
end;