unit Background;
{************************************************************************}
{*     by Michael Castle and Janice Keller                                                                                                    *}
{*     University of Michigan Mental Health Research Institute (MHRI)                                                     *}
{*     e-mail address: mike.castle@med.umich.edu                                                                                  *}
{************************************************************************}
{*     Rolling ball and rolling disk algorithms inspired by Stanley Sternberg's article, "Biomedical        *}
{*  Image Processing", IEEE Computer, January 1983.  Discussions with Rork Kuick and Tom               *}
{*  Ford-Holevinski also contributed to the development of the algorithms and improvements in their   *}
{*  efficiency.                                                                                                                                                *}
{************************************************************************}

interface

	uses
		QuickDraw, Palettes, PrintTraps, globals, Utilities, Graphics, Camera, Filters;

	procedure DoBackgroundMenuEvent (MenuItem: integer);


implementation

	type
		IntRow = array[0..MaxLine] of Integer;
		BackSubKindType = (RollingHorizontalArc, RollingVerticalArc, RollingBothArcs, RollingBall);

	var
		ArcTrimPer: integer;                                  {trim off percentage of each side of the rolling ball patch}
		shrinkfactor: integer;                                 {shrink the image and ball by this factor before rolling ball}
		BackSubKind: BackSubKindType;                {which kind of background subtraction are we doing}
		IntPlotWidth: Boolean;
		Intplotwidthval: integer;
		BoundRect: rect;
		xxcenter, yycenter: integer;                                        {center of rectangular mask used in MinIn2DMask}
		xmaskmin, ymaskmin: integer;                                   {upper left corner of mask used in AvgIn2DMask}
		backgroundptr, ballptr, smallimageptr: ptr;              {ptrs to background, rolling ball, shrunk image memory}
		backgroundaddr, balladdr, smallimageaddr: longint;   {addrs of background, rolling ball shrunk image memory}
		patchwidth: integer;                                                      {x or y dimension of the rolling ball patch}
		leftroll, rightroll, toproll, bottomroll: integer;         {bounds of the shrunk image}
		Aborting: boolean;



	function NewPtrToClearBuf (blockSize: Size): Ptr;
    {This function will return a pointer of size specified and will}
    {clear the memory to zeros . This is done to create an empty bit}
    {map containing nothing but white bits . }

    {MOVE . L  ( SP ) + , D0  ; get Size variable from stack}
    {_NewPtr , clear           ; make pointer }
    {MOVE.L  A0 , ( SP )       ; return pointer }
    {MOVE.W D0, MemErr     ; set up MemErr }
	inline
		$201F, $A31E, $2E88, $31C0, $0220;



	procedure RollBall;
{*******************************************************************************}
{*     RollBall 'rolls' a filtering object over a (shrunken) image in order to find the image's smooth continuous    *}
{*  background.  For the purpose of explaining this algorithm, imagine that the 2D grayscale image has a third     *}
{*  (height) dimension defined by the intensity value (0-255) at every point in the image.  The center of the      *}
{*  filtering object, a patch from the top of a sphere having radius BallRadius, is moved along each scan line of     *}
{*  the image so that the patch is tangent to the image at one or more points with every other point on the patch    *}
{*  below the corresponding (x,y) point of the image.  Any point either on or below the patch during this process*}
{*  is considered part of the background.  Shrinking the image before running this procedure is advised due to      *}
{*  the fourth-degree complexity of the algorithm.  Care has been taken to avoid unnecessary operations (exp.      *}
{*  multiplication inside loops) in this code.                                                                                                                *}
{*******************************************************************************}
		var
			halfpatchwidth,                                {distance in x or y from patch center to any edge}
			ptsbelowlastpatch,                           {number of points we may ignore because they were below last patch}
			left, right, top, bottom,                   {}
			xpt, ypt,                                           {current (x,y) point in the shrunken image}
			xpt2, ypt2,                                      {current (x,y) point in the patch relative to upper left corner}
			xval, yval,                                        {location in ball in shrunken image coordinates}
			zdif,                                                  {difference in z (height) between point on ball and point on image}
			zmin,                                                {smallest zdif for ball patch with center at current point}
			zctr,                                                 {current height of the center of the sphere of which the patch is a part}
			zadd: integer;                                    {height of a point on patch relative to the xy-plane of the shrunken image}
			ballpt,                                               {index to chunk of memory storing the precomputed ball patch}
			imgpt,                                               {index to chunk of memory storing the shrunken image}
			backgrpt,                                          {index to chunk of memory storing the calculated background}
			ybackgrpt,                                        {displacement to current background scan line}
			ybackgrinc,                                      {distance in memory between two shrunken y-points in background}
			smallimagewidth: longint;               {length of a scan line in shrunken image}
			p1, p2: ptr;                                      {temporary pointers to background, ball, or small image}
	begin
		UpdateMeter(20, 'Finding Background...');
		left := 1;
		right := rightroll - leftroll - 1;
		top := 1;
		bottom := bottomroll - toproll - 1;
		smallimagewidth := right - left + 3;
		halfpatchwidth := patchwidth div 2;
		ybackgrinc := shrinkfactor * (BoundRect.right - BoundRect.left);  {real dist btwn 2 adjacent (dy=1) shrunk pts}
		zctr := 0;                                            {start z-center in the xy-plane}
		for ypt := top to (bottom + patchwidth) do begin
				for xpt := left to (right + patchwidth) do {while patch is tangent to edges or within image...}
					begin                                           {xpt is far right edge of ball patch}
{do we have to move the patch up or down to make it tangent to but not above image?...}
						zmin := 255;                              {highest could ever be 255}
						ballpt := balladdr;
						ypt2 := ypt - patchwidth;          {ypt2 is top edge of ball patch}
						imgpt := smallimageaddr + ypt2 * smallimagewidth + xpt - patchwidth;
						while ypt2 <= ypt do begin
								xpt2 := xpt - patchwidth;      {xpt2 is far left edge of ball patch}
								while xpt2 <= xpt do            {check every point on ball patch}
									begin                                   {only examine points on }
										if ((xpt2 >= left) and (xpt2 <= right) and (ypt2 >= top) and (ypt2 <= bottom)) then begin
												p1 := ptr(ballpt);
												p2 := ptr(imgpt);
												zdif := BAND(p2^, 255) - (zctr + BAND(p1^, 255));  {curve - circle points}
												if (zdif < zmin) then begin {keep most negative, since ball should always be below curve}
														zmin := zdif;
													end;
											end;  {if xpt2,ypt2}
										ballpt := ballpt + 1;          {step thru the ball patch memory}
										xpt2 := xpt2 + 1;
										imgpt := imgpt + 1;
									end;  {while xpt2 }
								ypt2 := ypt2 + 1;
								imgpt := imgpt - patchwidth - 1 + smallimagewidth;
							end;  {while ypt2}
						if (zmin <> 0) then
							zctr := zctr + zmin;                {move ball up or down if we find a new minimum}
						if (zmin < 0) then
							ptsbelowlastpatch := halfpatchwidth    {ignore left half of ball patch when dz < 0}
						else
							ptsbelowlastpatch := 0;
{now compare every point on ball with background,  and keep highest number}
						yval := ypt - patchwidth;
						ypt2 := 0;
						ballpt := balladdr;
						ybackgrpt := backgroundaddr + longint(yval - top + 1) * ybackgrinc;
						while ypt2 <= patchwidth do begin
								xval := xpt - patchwidth + ptsbelowlastpatch;
								xpt2 := ptsbelowlastpatch;
								ballpt := ballpt + ptsbelowlastpatch;
								backgrpt := ybackgrpt + (xval - left + 1) * shrinkfactor;
								while xpt2 <= patchwidth do begin     {for all the points in the ball patch}
										if ((xval >= left) and (xval <= right) and (yval >= top) and (yval <= bottom)) then begin
												p1 := ptr(ballpt);
												zadd := zctr + BAND(p1^, 255);
												p1 := ptr(backgrpt);
												if (zadd > BAND(p1^, 255)) then  {keep largest adjustment}
													p1^ := zadd;
											end;
										ballpt := ballpt + 1;
										xval := xval + 1;
										xpt2 := xpt2 + 1;
										backgrpt := backgrpt + shrinkfactor;     {move to next point in x}
									end;  {while xpt2}
								yval := yval + 1;
								ypt2 := ypt2 + 1;
								ybackgrpt := ybackgrpt + ybackgrinc;       {move to next point in y}
							end;  {while ypt2}
					end;  {for xpt }
				if ((ypt mod 5) = 0) or not FasterBackgroundSubtraction then begin
						UpdateMeter(20 + ((ypt - top) * 70) div (bottom + patchwidth - top), 'Finding Background...');
						if CommandPeriod then begin
								beep;
								Aborting := true;
								Exit(RollBall);
							end;
					end;
			end;  {for ypt}
	end;


	function MinIn2DMask {(xmaskmin,ymaskmin: integer)}
		: integer;
{*******************************************************************************}
{*     MinInMask finds the minimum pixel value in a shrinkfactor X shrinkfactor mask.                                           *}
{*******************************************************************************}
		var
			i, j,                                           {loop indices to step through mask}
			thispixel,                                  {value at current pixel in mask}
			min,                                          {temporary minimum value in mask}
			nextrowoffset: integer;             {distance in memory from end of mask in this row to beginning in next}
			paddr: longint;                           {address of current mask pixel}
			p: ptr;                                        {pointer to current pixel in mask}
	begin
		with info^ do begin
				min := 255;
				nextrowoffset := bytesperrow - shrinkfactor;
				paddr := ord4(PicBaseAddr) + longint(ymaskmin) * bytesperrow + xmaskmin;
				for j := 1 to shrinkfactor do begin
						for i := 1 to shrinkfactor do begin
								p := ptr(paddr);
								thispixel := BAND(p^, 255);
								if (thispixel < min) then
									min := thispixel;
								paddr := paddr + 1;
							end;     {for i}
						paddr := paddr + nextrowoffset;
					end;     {for j}
				MinIn2DMask := min;
			end; {with}
	end;


	procedure GetRollingBall;
{******************************************************************************}
{*     This procedure computes the location of each point on the rolling ball patch relative to the center of the     *}
{*  sphere containing it.  The patch is located in the top half of this sphere.  The vertical axis of the sphere         *}
{*  passes through the center of the patch.  The projection of the patch in the xy-plane below is a square.           *}
{******************************************************************************}
		var
			rsquare,                                                                         {rolling ball radius squared}
			xtrim,                                                                            {# of pixels trimmed off each end of ball to make patch}
			xval, yval,                                                                     {x,y-values on patch relative to center of rolling ball}
			smallballradius, diam,                                                  {radius and diameter of rolling ball}
			temp,                                                                             {value must be >=0 to take square root}
			halfpatchwidth: integer;                                                {distance in x or y from center of patch to any edge}
			i,                                                                                    {index into rolling ball patch memory}
			ballsize: Size;                                                                {size of rolling ball memory}
			p: ptr;                                                                            {pointer to current point on the ball patch}
	begin
		smallballradius := ballradius div shrinkfactor;           {operate on small-sized image with small-sized ball}
		if smallballradius < 1 then
			smallballradius := 1;
		rsquare := smallballradius * smallballradius;
		diam := smallballradius * 2;
		xtrim := (ArcTrimPer * diam) div 100;                      {only use a patch of the rolling ball}
		patchwidth := diam - xtrim - xtrim;
		halfpatchwidth := smallballradius - xtrim;                   {this is half the patch width}
		ballsize := longint(patchwidth + 1) * longint(patchwidth + 1);
		ballptr := NewPtrToClearBuf(ballsize);
		p := ballptr;
		for i := 0 to ballsize - 1 do begin
				xval := i mod (patchwidth + 1) - halfpatchwidth;
				yval := i div (patchwidth + 1) - halfpatchwidth;
				temp := rsquare - (xval * xval) - (yval * yval);
				if (temp >= 0) then
					p^ := round(sqrt(temp))
				else
					p^ := 0;
				p := ptr(ord4(p) + 1);
			end;
	end;


	procedure InterpolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
{******************************************************************************}
{*     This procedure uses bilinear interpolation to find the points in the full-scale background given the points *}
{*  from the shrunken image background.  Since the shrunken background is found from an image composed of    *}
{*  minima (over a sufficiently large mask), it is certain that no point in the full-scale interpolated                 *}
{*  background has a higher pixel value than the corresponding point in the original image.                                  *}
{******************************************************************************}
		var
			i, ii,                                                   {horizontal loop indices}
			j, jj,                                                  {vertical loop indices}
			hloc, vloc,                                          {position of current pixel in calculated background}
			vinc,                                                   {memory offset from current calculated pos to current interpolated pos}
			lastvalue, nextvalue: integer;           {calculated pixel values between which we are interpolating}
			p,                                                        {pointer to current interpolated pixel value}
			bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values between which we are interpolating}
	begin
		vloc := 0;
		with BoundRect do begin
				for j := 1 to bottomroll - toproll - 1 do begin     {interpolate to find background interior}
						hloc := 0;
						vloc := vloc + shrinkfactor;
						for i := 1 to rightroll - leftroll do begin
								hloc := hloc + shrinkfactor;
								bgnextptr := ptr(backgroundaddr + vloc * longint(right - left) + hloc);
								bglastptr := ptr(ord4(bgnextptr) - shrinkfactor);
								nextvalue := BAND(bgnextptr^, 255);
								lastvalue := BAND(bglastptr^, 255);
								for ii := 1 to shrinkfactor - 1 do begin     {interpolate horizontally}
										p := ptr(ord4(bgnextptr) - ii);
										p^ := lastvalue + (shrinkfactor - ii) * (nextvalue - lastvalue) div shrinkfactor;
									end;     {for ii}
								for ii := 0 to shrinkfactor - 1 do begin     {interpolate vertically}
										bglastptr := ptr(backgroundaddr + (vloc - shrinkfactor) * longint(right - left) + hloc - ii);
										bgnextptr := ptr(backgroundaddr + vloc * longint(right - left) + hloc - ii);
										lastvalue := BAND(bglastptr^, 255);
										nextvalue := BAND(bgnextptr^, 255);
										vinc := 0;
										for jj := 1 to shrinkfactor - 1 do begin
												vinc := vinc - right + left;
												p := ptr(ord4(bgnextptr) + vinc);
												p^ := lastvalue + (shrinkfactor - jj) * (nextvalue - lastvalue) div shrinkfactor;
											end;     {for jj}
									end;     {for ii}
							end;     {for i}
					end;     {for j}
			end;   {with boundrect}
	end;


	procedure ExtrapolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
{******************************************************************************}
{*     This procedure uses linear extrapolation to find pixel values on the top, left, right, and bottom edges of      *}
{*  the background.  First it finds the top and bottom edge points by extrapolating from the edges of the                *}
{*  calculated and interpolated background interior.  Then it uses the edge points on the new calculated,               *}
{*  interpolated, and extrapolated background to find all of the left and right edge points.  If extrapolation yields *}
{*  values below zero or above 255, then they are set to zero and 255 respectively.                                              *}
{******************************************************************************}
		var
			ii, jj,                                                 {horizontal and vertical loop indices}
			hloc, vloc,                                          {position of current pixel in calculated/interpolated background}
			edgeslope,                                          {difference of last two consecutive pixel values on an edge}
			pvalue,                                               {current extrapolated pixel value}
			lastvalue, nextvalue: integer;           {calculated pixel values from which we are extrapolating}
			p,                                                        {pointer to current extrapolated pixel value}
			bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values from which we are extrapolating}
	begin
		with BoundRect do begin
				for hloc := shrinkfactor to shrinkfactor * (rightroll - leftroll) - 1 do begin     {extrapolate on top and bottom}
						bglastptr := ptr(backgroundaddr + shrinkfactor * longint(right - left) + hloc);
						bgnextptr := ptr(backgroundaddr + (shrinkfactor + 1) * longint(right - left) + hloc);
						lastvalue := BAND(bglastptr^, 255);
						nextvalue := BAND(bgnextptr^, 255);
						edgeslope := nextvalue - lastvalue;
						p := bglastptr;
						pvalue := lastvalue;
						for jj := 1 to shrinkfactor do begin
								p := ptr(ord4(p) - right + left);
								pvalue := pvalue - edgeslope;
								if (pvalue < 0) then
									p^ := 0
								else if (pvalue > 255) then
									p^ := 255
								else
									p^ := pvalue;
							end;     {for jj}
						bglastptr := ptr(backgroundaddr + (shrinkfactor * (bottomroll - toproll - 1) - 1) * longint(right - left) + hloc);
						bgnextptr := ptr(backgroundaddr + shrinkfactor * (bottomroll - toproll - 1) * longint(right - left) + hloc);
						lastvalue := BAND(bglastptr^, 255);
						nextvalue := BAND(bgnextptr^, 255);
						edgeslope := nextvalue - lastvalue;
						p := bgnextptr;
						pvalue := nextvalue;
						for jj := 1 to (bottom - top - 1) - shrinkfactor * (bottomroll - toproll - 1) do begin
								p := ptr(ord4(p) + right - left);
								pvalue := pvalue + edgeslope;
								if (pvalue < 0) then
									p^ := 0
								else if (pvalue > 255) then
									p^ := 255
								else
									p^ := pvalue;
							end;     {for jj}
					end;     {for hloc}
				for vloc := top to bottom - 1 do begin     {extrapolate on left and right}
						bglastptr := ptr(backgroundaddr + (vloc - top) * longint(right - left) + shrinkfactor);
						bgnextptr := ptr(ord4(bglastptr) + 1);
						lastvalue := BAND(bglastptr^, 255);
						nextvalue := BAND(bgnextptr^, 255);
						edgeslope := nextvalue - lastvalue;
						p := bglastptr;
						pvalue := lastvalue;
						for ii := 1 to shrinkfactor do begin
								p := ptr(ord4(p) - 1);
								pvalue := pvalue - edgeslope;
								if (pvalue < 0) then
									p^ := 0
								else if (pvalue > 255) then
									p^ := 255
								else
									p^ := pvalue;
							end;     {for ii}
						bgnextptr := ptr(backgroundaddr + (vloc - top) * longint(right - left) + shrinkfactor * (rightroll - leftroll - 1) - 1);
						bglastptr := ptr(ord4(bgnextptr) - 1);
						lastvalue := BAND(bglastptr^, 255);
						nextvalue := BAND(bgnextptr^, 255);
						edgeslope := nextvalue - lastvalue;
						p := bgnextptr;
						pvalue := nextvalue;
						for ii := 1 to (right - left - 1) - shrinkfactor * (rightroll - leftroll - 1) + 1 do begin
								p := ptr(ord4(p) + 1);
								pvalue := pvalue + edgeslope;
								if (pvalue < 0) then
									p^ := 0
								else if (pvalue > 255) then
									p^ := 255
								else
									p^ := pvalue;
							end;     {for ii}
					end;     {for vloc}
			end;   {with BoundRect}
	end;


	procedure SubtractBackground2D;
{*****************************************************************************}
{*     This procedure subtracts each pixel from the calculated/interpolated/extrapolated background from the  *}
{*  corresponding pixel value in the original image.  The resulting image is stored in place of the original        *}
{*  image.  Any pixel subtractions with results below zero are given the value zero.                                           *}
{*****************************************************************************}
		var
			hloc, vloc,                                          {current pixel location in image and background}
			pvalue: integer;                                 {difference at current pixel location}
			offset,                                                 {offset in memory from beginning of original image to current scan line}
			backgrpt: LongInt;                              {offset to current point in background}
			p: ptr;                                                {temporary pointer to image or background points}
			Databand: Linetype;                           {current scan line in image}
			ControlKey: boolean;
	begin
		backgrpt := 0;
		ControlKey := ControlKeyDown;
		with Info^, BoundRect do begin
				for vloc := top to bottom - 1 do begin
						GetLine(0, vloc, pixelsperline, Databand);
						for hloc := left to right - 1 do begin
								p := ptr(backgroundaddr + backgrpt);
								pvalue := Databand[hloc] - BAND(p^, 255);
								if ControlKey then
									pvalue := BAND(p^, 255);
								if pvalue < 0 then
									Databand[hloc] := 0
								else
									Databand[hloc] := pvalue;
								backgrpt := backgrpt + 1;
							end;     {for}
						offset := LongInt(vloc) * BytesPerRow;
						p := ptr(ord4(PicBaseAddr) + offset);
						BlockMove(@Databand, p, pixelsperline);
					end;  {for}
			end;     {with}
	end;


	procedure Background2D;
{******************************************************************************}
{*     This procedure implements a rolling-ball algorithm for the removal of smooth continuous background       *}
{*  from a two-dimensional gel image.  It rolls the ball (actually a square patch on the top of a sphere) on a       *}
{*  low-resolution (by a factor of 'shrinkfactor' times) copy of the original image in order to increase speed     *}
{*  with little loss in accuracy.  It uses interpolation and extrapolation to blow the shrunk image to full size.     *}
{******************************************************************************}
		var
			tport: Grafptr;
			i,                                     {loop index for shrunk image memory}
			backgroundsize,              {size of the background memory}
			smallimagesize: Size;     {size of the shrunk image memory}
			p: ptr;                             {pointer to current pixel in shrunk image memory}
			table: FateTable;             {not used}
	begin
		ShowWatch;
		UpdateMeter(0, 'Building Rolling Ball...');
		GetPort(tPort);
		with Info^ do begin
				SetPort(GrafPtr(osPort));
				BoundRect := roiRect;
			end;
		GetRollingBall;                                                                  {precompute the rolling ball}
		UpdateMeter(3, 'Building Rolling Ball...');
		balladdr := ord4(ballptr);
		with BoundRect do begin
				leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
				rightroll := right div shrinkfactor - 1;                      {on which to roll ball}
				toproll := top div shrinkfactor;
				bottomroll := bottom div shrinkfactor - 1;
				backgroundsize := longint(right - left) * longint(bottom - top);
				backgroundptr := NewPtrToClearBuf(backgroundsize);
				Aborting := backgroundptr = nil;
				backgroundaddr := ord4(backgroundptr);
				smallimagesize := longint(rightroll - leftroll + 1) * longint(bottomroll - toproll + 1);
				smallimageptr := NewPtrToClearBuf(smallimagesize);
				Aborting := Aborting or (smallimageptr = nil);
				smallimageaddr := ord4(smallimageptr);
				if not aborting then begin
						UpdateMeter(6, 'Smoothing Image ');
						filter(unweightedAvg, 1, table);                                {smooth image before shrinking}
						UpdateMeter(10, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
						for i := 0 to smallimagesize - 1 do begin                {create a lower resolution image for ball-rolling}
								p := ptr(smallimageaddr + i);
								xmaskmin := left + shrinkfactor * (i mod (rightroll - leftroll + 1));
								ymaskmin := top + shrinkfactor * (i div (rightroll - leftroll + 1));
								p^ := MinIn2DMask;                            {each point in small image is minimum of its neighborhood}
							end;
						if not aborting then begin
								Undo;        {restore original unsmoothed image}
								RollBall;
							end;
					end
				else
					beep;
				if not Aborting then begin
						UpdateMeter(90, 'Interpolating Background...');
						InterpolateBackground2D;                              {interpolate to find background interior}
						UpdateMeter(95, 'Extrapolating Background...');
						ExtrapolateBackground2D;                             {extrapolate on top and bottom}
						UpdateMeter(98, 'Subtracting Background...');
						SubtractBackground2D;                                  {subtract background from original image}
						UpdateMeter(100, 'Subtracting Background...');
					end;
			end;   {with boundrect}
		DisposPtr(backgroundptr);                           {free up background, rolling ball, shrunk image memory}
		DisposPtr(ballptr);
		DisposPtr(smallimageptr);
		DisposeWindow(MeterWindow);
		MeterWindow := nil;
		SetPort(tPort);
	end;


	procedure RollArc (left, rightminusone, diam: integer; var background, ballpoints: IntRow; var Dataline: Linetype);
		var
			xpt, xpt2, xval, ydif, ymin, yctr, bpt, yadd: integer;
	begin
		for xpt := left to rightminusone do begin
				background[xpt] := -255;         {init background curve to minimum values}
			end;
		yctr := 0;                                   {start y-center at the x axis}
		for xpt := left to (rightminusone + diam - 1) do {while semicircle is tangent to edges or within curve...}
			begin                                       {xpt is far right edge of semi-circle}
{do we have to move the circle?...}
				ymin := 256;                          {highest could ever be 255}
				bpt := 0;
				xpt2 := xpt - diam;                {xpt2 is far left edge of semi-circle}
				while xpt2 <= xpt do            {check every point on semicircle}
					begin
						if ((xpt2 >= left) and (xpt2 <= rightminusone)) then begin  {only examine points on curve}
								ydif := dataline[xpt2] - (yctr + ballpoints[bpt]);                {curve minus circle points}
								if (ydif < ymin) then begin {keep most negative, since ball should always be below curve}
										ymin := ydif;
									end;
							end;  {if xpt2 }
						bpt := bpt + 1;
						xpt2 := xpt2 + 1;
					end;  {while xpt2 }
				if (ymin <> 256) then{if we found a new minimum...}
					yctr := yctr + ymin;   {move circle up or down}
{now compare every point on semi with background,  and keep highest number}
				xval := xpt - diam;
				xpt2 := 0;
				while xpt2 <= diam do begin  {for all the points in the semicircle}
						if ((xval >= left) and (xval <= rightminusone)) then begin
								yadd := yctr + ballpoints[xpt2];
								if (yadd > background[xval]) then  {keep largest adjustment}
									background[xval] := yadd;
							end;
						xval := xval + 1;
						xpt2 := xpt2 + 1;
					end;  {while xpt2}
			end;  {for xpt }
	end;


	function MinIn1DMask (var Databand: LineType; xcenter: integer): integer;
{*******************************************************************************}
{*     MinIn1DMask finds the minimum pixel value in a (2*shrinkfactor-1) mask about the point xcenter in the *}
{*  current line.  This code must run FAST because it gets called OFTEN!                                                                   *}
{*******************************************************************************}
		var
			i,                                                                              {index to pixels in the mask}
			temp: integer;                                                          {temporary minimum value}
	begin
		temp := Databand[xcenter - shrinkfactor + 1];
		for i := xcenter - shrinkfactor + 2 to xcenter + shrinkfactor - 1 do
			if (Databand[i] < temp) then
				temp := Databand[i];
		MinIn1DMask := temp;
	end;


{******************************************************************************}
{*     This procedure computes the location of each point on the rolling arc relative to the center of the circle     *}
{*  containing it.  The arc is located in the top half of this circle.  The vertical axis of the circle passes through  *}
{*  the midpoint of the arc.  The projection of the arc on the x-axis below is a line segment.                                 *}
{******************************************************************************}
	procedure GetRollingArc (var arcpoints: IntRow; var arcwidth: integer);
		var
			xpt,                                                                                 {x-point along arc}
			xval,                                                                               {x-point in arc array}
			rsquare,                                                                         {shrunken arc radius squared}
			xtrim,                                                                            {points to be trimmed off each end of arc}
			smallballradius,                                                            {radius of shrunken arc which actually rolls}
			diam: integer;                                                                 {diameter of shrunken arc's circle}
	begin
		smallballradius := ballradius div shrinkfactor;            { operate on small-sized image with small-sized ball}
		rsquare := smallballradius * smallballradius;
		for xpt := -smallballradius to smallballradius do        { find the ballpoints for arc based at  (x,y)=(0,0) }
			begin
				xval := xpt + smallballradius;                                     {offset, can't have negative index}
				arcpoints[xval] := round(sqrt(abs(rsquare - (xpt * xpt))));  {Ys are positive, top half of circle}
			end;
		diam := smallballradius * 2;
		xtrim := (ArcTrimPer * diam) div 100;                       {how many points to trim off each end}
		arcwidth := diam - xtrim - xtrim;
		for xpt := -smallballradius to smallballradius - xtrim - xtrim do begin
				xval := xpt + smallballradius;
				arcpoints[xval] := arcpoints[xval + xtrim];
			end;
		for xpt := smallballradius - xtrim - xtrim + 1 to smallballradius do begin
				xval := xpt + smallballradius;
				arcpoints[xval] := 0;
			end;
	end;


	procedure ExtrapolateBackground1D (var Backline, Dataline: LineType; background: IntRow; leftroll, rightroll: integer);
{******************************************************************************}
{*     This procedure uses linear extrapolation to find pixel values on the left and right edges of the current        *}
{*  line of the background.  It finds the edge points by extrapolating from the edges of the calculated and               *}
{*  interpolated background interior.  If extrapolation yields values below zero or above 255, then they are set *}
{*  to zero and 255 respectively.                                                                                                                               *}
{******************************************************************************}
		var
			i,                                                                             {index to edges of background array}
			hloc,                                                                       {}
			edgeslope: integer;                                                 {}
	begin
		with BoundRect do begin
				edgeslope := (background[leftroll + 1] - background[leftroll + 2]) div shrinkfactor;
				for i := left to shrinkfactor * (leftroll + 1) - 1 do begin     {extrapolate on left edge}
						hloc := shrinkfactor * (leftroll + 1) - 1 + left - i;
						if (Backline[hloc + 1] + edgeslope < 0) then
							Backline[hloc] := 0
						else if (Backline[hloc + 1] + edgeslope > Dataline[hloc]) then
							Backline[hloc] := Dataline[hloc]
						else
							Backline[hloc] := Backline[hloc + 1] + edgeslope;
					end;
				edgeslope := (background[rightroll] - background[rightroll - 1]) div shrinkfactor;
				for hloc := shrinkfactor * rightroll to right - 1 do begin     {extrapolate on right edge}
						if (Backline[hloc - 1] + edgeslope < 0) then
							Backline[hloc] := 0
						else if (Backline[hloc - 1] + edgeslope > Dataline[hloc]) then
							Backline[hloc] := Dataline[hloc]
						else
							Backline[hloc] := Backline[hloc - 1] + edgeslope;
					end;
			end;     {with}
	end;

	procedure Background1D;
		var
			tport: Grafptr;
			hloc, vloc, arcwidth, leftroll, rightroll, numpixels: integer;
			left, right, top, bottom: integer;                      {image bounds; ROTATED if RollingVerticalArc}
			i, j, maskwidth: integer;
			background, arcpoints: IntRow;
			offset: LongInt;
			p: ptr;
			Dataline, Backline, Smalldataline: Linetype;
			str: str255;
	begin
		ShowWatch;
		UpdateMeter(0, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
		GetPort(tPort);
		with Info^ do begin
				SetPort(GrafPtr(osPort));
				BoundRect := roiRect;
			end;
		GetRollingArc(arcpoints, arcwidth);
		maskwidth := shrinkfactor + shrinkfactor - 1;
		case BackSubKind of
			RollingHorizontalArc:  begin
					left := BoundRect.left;
					top := BoundRect.top;
					right := BoundRect.right;
					bottom := BoundRect.bottom;
					numpixels := Info^.pixelsperline;
					str := 'Rolling Disk Horizontally...';
				end;
			RollingVerticalArc:  begin
					left := BoundRect.top;
					top := BoundRect.left;
					right := BoundRect.bottom;
					bottom := BoundRect.right;
					numpixels := Info^.nlines;
					str := 'Rolling Disk Vertically...';
				end;
		end;     {case}
		leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
		rightroll := right div shrinkfactor - 1;                      {on which to roll arc}
		for vloc := top to bottom - 1 do begin  {for ROI}
				case BackSubKind of
					RollingHorizontalArc: 
						GetLine(0, vloc, numpixels, Dataline);
					RollingVerticalArc: 
						GetColumn(vloc, 0, numpixels, Dataline);
				end;     {case}
				for i := leftroll + 1 to rightroll do
					smalldataline[i] := MinIn1DMask(Dataline, shrinkfactor * i - 1);
				RollArc(leftroll + 1, rightroll, arcwidth, background, arcpoints, smalldataline);  {roll arc on one line}
				for i := leftroll + 1 to rightroll do begin           {interpolate to find interior background points}
						hloc := shrinkfactor * i - 1;
						Backline[hloc] := background[i];
						for j := 1 to shrinkfactor - 1 do
							Backline[hloc - j] := background[i - 1] + (shrinkfactor - j) * (background[i] - background[i - 1]) div shrinkfactor;
					end;
				ExtrapolateBackground1D(Backline, Dataline, background, leftroll, rightroll);
				for i := left to right - 1 do begin                                {subtract background from current scan line}
						Dataline[i] := Dataline[i] - Backline[i];
						if Dataline[i] < 0 then
							Dataline[i] := 0;
					end;
				case BackSubKind of
					RollingHorizontalArc: 
						with Info^ do begin
								offset := LongInt(vloc) * BytesPerRow;
								p := ptr(ord4(PicBaseAddr) + offset);
								BlockMove(@Dataline, p, numpixels);            {fast whole line write}
							end;
					RollingVerticalArc: 
						PutColumn(vloc, 0, numpixels, Dataline);         {slow whole column write}
				end;     {case}
				if ((vloc mod 8) = 0) and (vloc > 16) then begin
						UpdateMeter((LongInt(vloc - top) * 100) div (bottom - top - 1), str);
						if CommandPeriod then begin
								beep;
								Aborting := true;
								leave;
							end;
					end;
			end;
		UpdateMeter(100, str);
		DisposeWindow(MeterWindow);
		MeterWindow := nil;
		SetPort(tPort);
	end;

	procedure SetUpGel;
		var
			frame: rect;
			AutoSelectAll: boolean;
			p: ptr;
	begin
		if NotinBounds or NotRectangular then
			exit(SetUpGel);
		StopDigitizing;
		AutoSelectAll := not Info^.RoiShowing;
		if AutoSelectAll then
			SelectAll(false);
		SetupUndoFromClip;
		with info^ do begin
				frame := roiRect;
				if ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction) then
					ApplyLookupTable;
				changes := true;
			end;
		case BackSubKind of
			RollingHorizontalArc, RollingVerticalArc: 
				Background1D;    {--------------> call background subtract <-------------------}
			RollingBall: 
				Background2D;
			RollingBothArcs:  begin
					BackSubKind := RollingHorizontalArc;           {remove horizontal streaks}
					Background1D;
					BackSubKind := RollingVerticalArc;               {remove vertical streaks}
					if not aborting then
						Background1D;
					BackSubKind := RollingBothArcs;                   {leave BackSubKind as we found it}
				end;
		end;     {case}
		UpdatePicWindow;
		SetUpRoiRect;
		WhatToUndo := UndoFilter;
		Info^.changes := true;
		ShowWatch;
		if AutoSelectAll then
			KillRoi;
		if Aborting then begin
				Undo;
				WhatToUndo := NothingToUndo;
				UpdatePicWindow;
			end;
	end;


	procedure GetBallRadius;
		var
			SaveRadius: integer;
			canceled: boolean;
	begin
		SaveRadius := BallRadius;
		BallRadius := GetInt('Rolling BallRadius:', BallRadius, canceled);
		if (BallRadius < 1) or (BallRadius > 319) or canceled then begin
				BallRadius := SaveRadius;
				if not canceled then
					beep;
			end;
	end;


	procedure DoBackgroundMenuEvent (MenuItem: integer);
		var
			map_array: Ptr;
	begin
		MeterWindow := nil;
		Aborting := false;
		case MenuItem of
			HorizontalItem, VerticalItem:  begin
					if FasterBackgroundSubtraction then begin
							ArcTrimPer := 20;
							shrinkfactor := 4;
						end
					else begin
							ArcTrimPer := 10;
							shrinkfactor := 2;
						end;
					if MenuItem = HorizontalItem then
						BackSubKind := RollingHorizontalArc
					else
						BackSubKind := RollingVerticalArc;
					SetUpGel;
				end;
			Sub2DItem:  begin
					if FasterBackgroundSubtraction then begin
							if BallRadius > 15 then begin
									ArcTrimPer := 20;     {trim 40% in x and y}
									shrinkfactor := 8;
								end
							else begin
									ArcTrimPer := 16;  {trim 32% in x and y}
									shrinkfactor := 4;
								end
						end
					else begin  {faster not checked}
							if BallRadius > 15 then begin
									ArcTrimPer := 16;  {trim 32% in x and y}
									shrinkfactor := 4;
								end
							else begin
									ArcTrimPer := 12;   {trim 24% in x and y}
									ShrinkFactor := 2;
								end
						end;
					BackSubKind := RollingBall;
					SetUpGel;
				end;
			RemoveStreaksItem:  begin
					ArcTrimPer := 20;
					shrinkfactor := 4;
					BackSubKind := RollingBothArcs;
					SetUpGel;
				end;
			FasterItem: 
				FasterBackgroundSubtraction := not FasterBackgroundSubtraction;
			RadiusItem: 
				GetBallRadius;
		end; {case}
	end;


end.