home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / mac / developm / source / image144.sit / Background.p < prev    next >
Encoding:
Text File  |  1992-01-17  |  39.9 KB  |  819 lines
  1. unit Background;
  2. {************************************************************************}
  3. {*     by Michael Castle and Janice Keller                                                                                                    *}
  4. {*     University of Michigan Mental Health Research Institute (MHRI)                                                     *}
  5. {*     e-mail address: mike_castle@ub.cc.umich.edu                                                                                   *}
  6. {************************************************************************}
  7. {*     Rolling ball and rolling disk algorithms inspired by Stanley Sternberg's article, "Biomedical        *}
  8. {*  Image Processing", IEEE Computer, January 1983.  Discussions with Rork Kuick and Tom               *}
  9. {*  Ford-Holevinski also contributed to the development of the algorithms and improvements in their   *}
  10. {*  efficiency.                                                                                                                                                *}
  11. {************************************************************************}
  12.  
  13. interface
  14.  
  15.     uses
  16.         QuickDraw, Palettes, PrintTraps, globals, Utilities, Graphics, Camera, Functions;
  17.  
  18.     procedure DoBackgroundMenuEvent (MenuItem: integer);
  19.  
  20.  
  21. implementation
  22.  
  23.     type
  24.         IntRow = array[0..MaxLine] of Integer;
  25.         BackSubKindType = (RollingHorizontalArc, RollingVerticalArc, RollingBothArcs, RollingBall);
  26.  
  27.     var
  28.         ArcTrimPer: integer;                                  {trim off percentage of each side of the rolling ball patch}
  29.         shrinkfactor: integer;                                 {shrink the image and ball by this factor before rolling ball}
  30.         BackSubKind: BackSubKindType;                {which kind of background subtraction are we doing}
  31.         IntPlotWidth: Boolean;
  32.         Intplotwidthval: integer;
  33.         BoundRect: rect;
  34.         xxcenter, yycenter: integer;                                        {center of rectangular mask used in MinIn2DMask}
  35.         backgroundptr, ballptr, smallimageptr: ptr;              {ptrs to background, rolling ball, shrunk image memory}
  36.         backgroundaddr, balladdr, smallimageaddr: longint;   {addrs of background, rolling ball shrunk image memory}
  37.         patchwidth: integer;                                                      {x or y dimension of the rolling ball patch}
  38.         leftroll, rightroll, toproll, bottomroll: integer;         {bounds of the shrunk image}
  39.         Aborting: boolean;
  40.  
  41.  
  42.  
  43.     function NewPtrToClearBuf (blockSize: Size): Ptr;
  44.     {This function will return a pointer of size specified and will}
  45.     {clear the memory to zeros . This is done to create an empty bit}
  46.     {map containing nothing but white bits . }
  47.  
  48.     {MOVE . L  ( SP ) + , D0  ; get Size variable from stack}
  49.     {_NewPtr , clear           ; make pointer }
  50.     {MOVE.L  A0 , ( SP )       ; return pointer }
  51.     {MOVE.W D0, MemErr     ; set up MemErr }
  52.     inline
  53.         $201F, $A31E, $2E88, $31C0, $0220;
  54.  
  55.  
  56.  
  57.     procedure RollBall;
  58. {*******************************************************************************}
  59. {*     RollBall 'rolls' a filtering object over a (shrunken) image in order to find the image's smooth continuous    *}
  60. {*  background.  For the purpose of explaining this algorithm, imagine that the 2D grayscale image has a third     *}
  61. {*  (height) dimension defined by the intensity value (0-255) at every point in the image.  The center of the      *}
  62. {*  filtering object, a patch from the top of a sphere having radius BallRadius, is moved along each scan line of     *}
  63. {*  the image so that the patch is tangent to the image at one or more points with every other point on the patch    *}
  64. {*  below the corresponding (x,y) point of the image.  Any point either on or below the patch during this process*}
  65. {*  is considered part of the background.  Shrinking the image before running this procedure is advised due to      *}
  66. {*  the fourth-degree complexity of the algorithm.  Care has been taken to avoid unnecessary operations (exp.      *}
  67. {*  multiplication inside loops) in this code.                                                                                                                *}
  68. {*******************************************************************************}
  69.         var
  70.             halfpatchwidth,                                {distance in x or y from patch center to any edge}
  71.             ptsbelowlastpatch,                           {number of points we may ignore because they were below last patch}
  72.             left, right, top, bottom,                   {}
  73.             xpt, ypt,                                           {current (x,y) point in the shrunken image}
  74.             xpt2, ypt2,                                      {current (x,y) point in the patch relative to upper left corner}
  75.             xval, yval,                                        {location in ball in shrunken image coordinates}
  76.             zdif,                                                  {difference in z (height) between point on ball and point on image}
  77.             zmin,                                                {smallest zdif for ball patch with center at current point}
  78.             zctr,                                                 {current height of the center of the sphere of which the patch is a part}
  79.             zadd: integer;                                    {height of a point on patch relative to the xy-plane of the shrunken image}
  80.             ballpt,                                               {index to chunk of memory storing the precomputed ball patch}
  81.             imgpt,                                               {index to chunk of memory storing the shrunken image}
  82.             backgrpt,                                          {index to chunk of memory storing the calculated background}
  83.             ybackgrpt,                                        {displacement to current background scan line}
  84.             ybackgrinc,                                      {distance in memory between two shrunken y-points in background}
  85.             smallimagewidth: longint;               {length of a scan line in shrunken image}
  86.             p1, p2: ptr;                                      {temporary pointers to background, ball, or small image}
  87.     begin
  88.         UpdateMeter(20, 'Finding Background...');
  89.         left := 1;
  90.         right := rightroll - leftroll - 1;
  91.         top := 1;
  92.         bottom := bottomroll - toproll - 1;
  93.         smallimagewidth := right - left + 3;
  94.         halfpatchwidth := patchwidth div 2;
  95.         ybackgrinc := shrinkfactor * (BoundRect.right - BoundRect.left);  {real dist btwn 2 adjacent (dy=1) shrunk pts}
  96.         zctr := 0;                                            {start z-center in the xy-plane}
  97.         for ypt := top to (bottom + patchwidth) do begin
  98.                 for xpt := left to (right + patchwidth) do {while patch is tangent to edges or within image...}
  99.                     begin                                           {xpt is far right edge of ball patch}
  100. {do we have to move the patch up or down to make it tangent to but not above image?...}
  101.                         zmin := 255;                              {highest could ever be 255}
  102.                         ballpt := balladdr;
  103.                         ypt2 := ypt - patchwidth;          {ypt2 is top edge of ball patch}
  104.                         imgpt := smallimageaddr + ypt2 * smallimagewidth + xpt - patchwidth;
  105.                         while ypt2 <= ypt do begin
  106.                                 xpt2 := xpt - patchwidth;      {xpt2 is far left edge of ball patch}
  107.                                 while xpt2 <= xpt do            {check every point on ball patch}
  108.                                     begin                                   {only examine points on }
  109.                                         if ((xpt2 >= left) and (xpt2 <= right) and (ypt2 >= top) and (ypt2 <= bottom)) then begin
  110.                                                 p1 := ptr(ballpt);
  111.                                                 p2 := ptr(imgpt);
  112.                                                 zdif := BAND(p2^, 255) - (zctr + BAND(p1^, 255));  {curve - circle points}
  113.                                                 if (zdif < zmin) then begin {keep most negative, since ball should always be below curve}
  114.                                                         zmin := zdif;
  115.                                                     end;
  116.                                             end;  {if xpt2,ypt2}
  117.                                         ballpt := ballpt + 1;          {step thru the ball patch memory}
  118.                                         xpt2 := xpt2 + 1;
  119.                                         imgpt := imgpt + 1;
  120.                                     end;  {while xpt2 }
  121.                                 ypt2 := ypt2 + 1;
  122.                                 imgpt := imgpt - patchwidth - 1 + smallimagewidth;
  123.                             end;  {while ypt2}
  124.                         if (zmin <> 0) then
  125.                             zctr := zctr + zmin;                {move ball up or down if we find a new minimum}
  126.                         if (zmin < 0) then
  127.                             ptsbelowlastpatch := halfpatchwidth    {ignore left half of ball patch when dz < 0}
  128.                         else
  129.                             ptsbelowlastpatch := 0;
  130. {now compare every point on ball with background,  and keep highest number}
  131.                         yval := ypt - patchwidth;
  132.                         ypt2 := 0;
  133.                         ballpt := balladdr;
  134.                         ybackgrpt := backgroundaddr + longint(yval - top + 1) * ybackgrinc;
  135.                         while ypt2 <= patchwidth do begin
  136.                                 xval := xpt - patchwidth + ptsbelowlastpatch;
  137.                                 xpt2 := ptsbelowlastpatch;
  138.                                 ballpt := ballpt + ptsbelowlastpatch;
  139.                                 backgrpt := ybackgrpt + (xval - left + 1) * shrinkfactor;
  140.                                 while xpt2 <= patchwidth do begin     {for all the points in the ball patch}
  141.                                         if ((xval >= left) and (xval <= right) and (yval >= top) and (yval <= bottom)) then begin
  142.                                                 p1 := ptr(ballpt);
  143.                                                 zadd := zctr + BAND(p1^, 255);
  144.                                                 p1 := ptr(backgrpt);
  145.                                                 if (zadd > BAND(p1^, 255)) then  {keep largest adjustment}
  146.                                                     p1^ := zadd;
  147.                                             end;
  148.                                         ballpt := ballpt + 1;
  149.                                         xval := xval + 1;
  150.                                         xpt2 := xpt2 + 1;
  151.                                         backgrpt := backgrpt + shrinkfactor;     {move to next point in x}
  152.                                     end;  {while xpt2}
  153.                                 yval := yval + 1;
  154.                                 ypt2 := ypt2 + 1;
  155.                                 ybackgrpt := ybackgrpt + ybackgrinc;       {move to next point in y}
  156.                             end;  {while ypt2}
  157.                     end;  {for xpt }
  158.                 if ((ypt mod 5) = 0) or not FasterBackgroundSubtraction then begin
  159.                         UpdateMeter(20 + ((ypt - top) * 70) div (bottom + patchwidth - top), 'Finding Background...');
  160.                         if CommandPeriod then begin
  161.                                 beep;
  162.                                 Aborting := true;
  163.                                 Exit(RollBall);
  164.                             end;
  165.                     end;
  166.             end;  {for ypt}
  167.     end;
  168.  
  169.  
  170.     function MinIn2DMask {(xcenter, ycenter: integer)}
  171.         : integer;
  172. {*******************************************************************************}
  173. {*     MinInMask finds the minimum pixel value in a (2*shrinkfactor-1) X (2*shrinkfactor-1) mask about the *}
  174. {*  point (xcenter,ycenter).  This code must run FAST because it gets called OFTEN!                                                *}
  175. {*******************************************************************************}
  176.         var
  177.             i, j,                                                               {loop indices to step through mask}
  178.             thispixel,                                                      {value at current pixel in mask}
  179.             temp,                                                            {temporary minimum value in mask}
  180.             xmaskmin, xmaskmax,                                {x-bounds of the mask}
  181.             ymaskmin, ymaskmax,                                {y-bounds of the mask}
  182.             bytesperrow,                                                {number of bytes in one image or ROI scan line}
  183.             nextrowoffset: integer;                                {distance in memory from end of mask in this row to beginning in next}
  184.             paddr: longint;                                              {address of current mask pixel}
  185.             p: ptr;                                                           {pointer to current pixel in mask}
  186.     begin
  187.         temp := 255;
  188.         bytesperrow := Info^.BytesPerRow;
  189.         xmaskmin := xxcenter - shrinkfactor div 2 + 1;
  190.         xmaskmax := xxcenter + shrinkfactor div 2;
  191.         ymaskmin := yycenter - shrinkfactor div 2 + 1;
  192.         ymaskmax := yycenter + shrinkfactor div 2;
  193.         nextrowoffset := bytesperrow - xmaskmax + xmaskmin - 1;
  194.         paddr := ord4(Info^.PicBaseAddr) + longint(ymaskmin) * bytesperrow + xmaskmin;
  195.         for j := ymaskmin to ymaskmax do begin
  196.                 for i := xmaskmin to xmaskmax do begin
  197.                         p := ptr(paddr);
  198.                         thispixel := BAND(p^, 255);
  199.                         if (thispixel < temp) then
  200.                             temp := thispixel;
  201.                         paddr := paddr + 1;
  202.                     end;     {for i}
  203.                 paddr := paddr + nextrowoffset;
  204.             end;     {for j}
  205.         MinIn2DMask := temp;
  206.     end;
  207.  
  208.  
  209.     procedure GetRollingBall;
  210. {******************************************************************************}
  211. {*     This procedure computes the location of each point on the rolling ball patch relative to the center of the     *}
  212. {*  sphere containing it.  The patch is located in the top half of this sphere.  The vertical axis of the sphere         *}
  213. {*  passes through the center of the patch.  The projection of the patch in the xy-plane below is a square.           *}
  214. {******************************************************************************}
  215.         var
  216.             rsquare,                                                                         {rolling ball radius squared}
  217.             xtrim,                                                                            {# of pixels trimmed off each end of ball to make patch}
  218.             xval, yval,                                                                     {x,y-values on patch relative to center of rolling ball}
  219.             smallballradius, diam,                                                  {radius and diameter of rolling ball}
  220.             temp,                                                                             {value must be >=0 to take square root}
  221.             halfpatchwidth: integer;                                                {distance in x or y from center of patch to any edge}
  222.             i,                                                                                    {index into rolling ball patch memory}
  223.             ballsize: Size;                                                                {size of rolling ball memory}
  224.             p: ptr;                                                                            {pointer to current point on the ball patch}
  225.     begin
  226.         smallballradius := ballradius div shrinkfactor;            {operate on small-sized image with small-sized ball}
  227.         rsquare := smallballradius * smallballradius;
  228.         diam := smallballradius * 2;
  229.         xtrim := (ArcTrimPer * diam) div 100;                      {only use a patch of the rolling ball}
  230.         patchwidth := diam - xtrim - xtrim;
  231.         halfpatchwidth := smallballradius - xtrim;                   {this is half the patch width}
  232.         ballsize := longint(patchwidth + 1) * longint(patchwidth + 1);
  233.         ballptr := NewPtrToClearBuf(ballsize);
  234.         p := ballptr;
  235.         for i := 0 to ballsize - 1 do begin
  236.                 xval := i mod (patchwidth + 1) - halfpatchwidth;
  237.                 yval := i div (patchwidth + 1) - halfpatchwidth;
  238.                 temp := rsquare - (xval * xval) - (yval * yval);
  239.                 if (temp >= 0) then
  240.                     p^ := round(sqrt(temp))
  241.                 else
  242.                     p^ := 0;
  243.                 p := ptr(ord4(p) + 1);
  244.             end;
  245.     end;
  246.  
  247.  
  248.     procedure InterpolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
  249. {******************************************************************************}
  250. {*     This procedure uses bilinear interpolation to find the points in the full-scale background given the points *}
  251. {*  from the shrunken image background.  Since the shrunken background is found from an image composed of    *}
  252. {*  minima (over a sufficiently large mask), it is certain that no point in the full-scale interpolated                 *}
  253. {*  background has a higher pixel value than the corresponding point in the original image.                                  *}
  254. {******************************************************************************}
  255.         var
  256.             i, ii,                                                   {horizontal loop indices}
  257.             j, jj,                                                  {vertical loop indices}
  258.             hloc, vloc,                                          {position of current pixel in calculated background}
  259.             vinc,                                                   {memory offset from current calculated pos to current interpolated pos}
  260.             lastvalue, nextvalue: integer;           {calculated pixel values between which we are interpolating}
  261.             p,                                                        {pointer to current interpolated pixel value}
  262.             bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values between which we are interpolating}
  263.     begin
  264.         vloc := 0;
  265.         with BoundRect do begin
  266.                 for j := 1 to bottomroll - toproll - 1 do begin     {interpolate to find background interior}
  267.                         hloc := 0;
  268.                         vloc := vloc + shrinkfactor;
  269.                         for i := 1 to rightroll - leftroll do begin
  270.                                 hloc := hloc + shrinkfactor;
  271.                                 bgnextptr := ptr(backgroundaddr + vloc * longint(right - left) + hloc);
  272.                                 bglastptr := ptr(ord4(bgnextptr) - shrinkfactor);
  273.                                 nextvalue := BAND(bgnextptr^, 255);
  274.                                 lastvalue := BAND(bglastptr^, 255);
  275.                                 for ii := 1 to shrinkfactor - 1 do begin     {interpolate horizontally}
  276.                                         p := ptr(ord4(bgnextptr) - ii);
  277.                                         p^ := lastvalue + (shrinkfactor - ii) * (nextvalue - lastvalue) div shrinkfactor;
  278.                                     end;     {for ii}
  279.                                 for ii := 0 to shrinkfactor - 1 do begin     {interpolate vertically}
  280.                                         bglastptr := ptr(backgroundaddr + (vloc - shrinkfactor) * longint(right - left) + hloc - ii);
  281.                                         bgnextptr := ptr(backgroundaddr + vloc * longint(right - left) + hloc - ii);
  282.                                         lastvalue := BAND(bglastptr^, 255);
  283.                                         nextvalue := BAND(bgnextptr^, 255);
  284.                                         vinc := 0;
  285.                                         for jj := 1 to shrinkfactor - 1 do begin
  286.                                                 vinc := vinc - right + left;
  287.                                                 p := ptr(ord4(bgnextptr) + vinc);
  288.                                                 p^ := lastvalue + (shrinkfactor - jj) * (nextvalue - lastvalue) div shrinkfactor;
  289.                                             end;     {for jj}
  290.                                     end;     {for ii}
  291.                             end;     {for i}
  292.                     end;     {for j}
  293.             end;   {with boundrect}
  294.     end;
  295.  
  296.     procedure ExtrapolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
  297. {******************************************************************************}
  298. {*     This procedure uses linear extrapolation to find pixel values on the top, left, right, and bottom edges of      *}
  299. {*  the background.  First it finds the top and bottom edge points by extrapolating from the edges of the                *}
  300. {*  calculated and interpolated background interior.  Then it uses the edge points on the new calculated,               *}
  301. {*  interpolated, and extrapolated background to find all of the left and right edge points.  If extrapolation yields *}
  302. {*  values below zero or above 255, then they are set to zero and 255 respectively.                                              *}
  303. {******************************************************************************}
  304.         var
  305.             ii, jj,                                                 {horizontal and vertical loop indices}
  306.             hloc, vloc,                                          {position of current pixel in calculated/interpolated background}
  307.             edgeslope,                                          {difference of last two consecutive pixel values on an edge}
  308.             pvalue,                                               {current extrapolated pixel value}
  309.             lastvalue, nextvalue: integer;           {calculated pixel values from which we are extrapolating}
  310.             p,                                                        {pointer to current extrapolated pixel value}
  311.             bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values from which we are extrapolating}
  312.     begin
  313.         with BoundRect do begin
  314.                 for hloc := shrinkfactor to shrinkfactor * (rightroll - leftroll) - 1 do begin     {extrapolate on top and bottom}
  315.                         bglastptr := ptr(backgroundaddr + shrinkfactor * longint(right - left) + hloc);
  316.                         bgnextptr := ptr(backgroundaddr + (shrinkfactor + 1) * longint(right - left) + hloc);
  317.                         lastvalue := BAND(bglastptr^, 255);
  318.                         nextvalue := BAND(bgnextptr^, 255);
  319.                         edgeslope := nextvalue - lastvalue;
  320.                         p := bglastptr;
  321.                         pvalue := lastvalue;
  322.                         for jj := 1 to shrinkfactor do begin
  323.                                 p := ptr(ord4(p) - right + left);
  324.                                 pvalue := pvalue - edgeslope;
  325.                                 if (pvalue < 0) then
  326.                                     p^ := 0
  327.                                 else if (pvalue > 255) then
  328.                                     p^ := 255
  329.                                 else
  330.                                     p^ := pvalue;
  331.                             end;     {for jj}
  332.                         bglastptr := ptr(backgroundaddr + (shrinkfactor * (bottomroll - toproll - 1) - 1) * longint(right - left) + hloc);
  333.                         bgnextptr := ptr(backgroundaddr + shrinkfactor * (bottomroll - toproll - 1) * longint(right - left) + hloc);
  334.                         lastvalue := BAND(bglastptr^, 255);
  335.                         nextvalue := BAND(bgnextptr^, 255);
  336.                         edgeslope := nextvalue - lastvalue;
  337.                         p := bgnextptr;
  338.                         pvalue := nextvalue;
  339.                         for jj := 1 to (bottom - top - 1) - shrinkfactor * (bottomroll - toproll - 1) do begin
  340.                                 p := ptr(ord4(p) + right - left);
  341.                                 pvalue := pvalue + edgeslope;
  342.                                 if (pvalue < 0) then
  343.                                     p^ := 0
  344.                                 else if (pvalue > 255) then
  345.                                     p^ := 255
  346.                                 else
  347.                                     p^ := pvalue;
  348.                             end;     {for jj}
  349.                     end;     {for hloc}
  350.                 for vloc := top to bottom - 1 do begin     {extrapolate on left and right}
  351.                         bglastptr := ptr(backgroundaddr + (vloc - top) * longint(right - left) + shrinkfactor);
  352.                         bgnextptr := ptr(ord4(bglastptr) + 1);
  353.                         lastvalue := BAND(bglastptr^, 255);
  354.                         nextvalue := BAND(bgnextptr^, 255);
  355.                         edgeslope := nextvalue - lastvalue;
  356.                         p := bglastptr;
  357.                         pvalue := lastvalue;
  358.                         for ii := 1 to shrinkfactor do begin
  359.                                 p := ptr(ord4(p) - 1);
  360.                                 pvalue := pvalue - edgeslope;
  361.                                 if (pvalue < 0) then
  362.                                     p^ := 0
  363.                                 else if (pvalue > 255) then
  364.                                     p^ := 255
  365.                                 else
  366.                                     p^ := pvalue;
  367.                             end;     {for ii}
  368.                         bgnextptr := ptr(backgroundaddr + (vloc - top) * longint(right - left) + shrinkfactor * (rightroll - leftroll - 1) - 1);
  369.                         bglastptr := ptr(ord4(bgnextptr) - 1);
  370.                         lastvalue := BAND(bglastptr^, 255);
  371.                         nextvalue := BAND(bgnextptr^, 255);
  372.                         edgeslope := nextvalue - lastvalue;
  373.                         p := bgnextptr;
  374.                         pvalue := nextvalue;
  375.                         for ii := 1 to (right - left - 1) - shrinkfactor * (rightroll - leftroll - 1) + 1 do begin
  376.                                 p := ptr(ord4(p) + 1);
  377.                                 pvalue := pvalue + edgeslope;
  378.                                 if (pvalue < 0) then
  379.                                     p^ := 0
  380.                                 else if (pvalue > 255) then
  381.                                     p^ := 255
  382.                                 else
  383.                                     p^ := pvalue;
  384.                             end;     {for ii}
  385.                     end;     {for vloc}
  386.             end;   {with BoundRect}
  387.     end;
  388.  
  389.  
  390.     procedure SubtractBackground2D;
  391. {*****************************************************************************}
  392. {*     This procedure subtracts each pixel from the calculated/interpolated/extrapolated background from the  *}
  393. {*  corresponding pixel value in the original image.  The resulting image is stored in place of the original        *}
  394. {*  image.  Any pixel subtractions with results below zero are given the value zero.                                           *}
  395. {*****************************************************************************}
  396.         var
  397.             hloc, vloc,                                          {current pixel location in image and background}
  398.             pvalue: integer;                                 {difference at current pixel location}
  399.             offset,                                                 {offset in memory from beginning of original image to current scan line}
  400.             backgrpt: LongInt;                              {offset to current point in background}
  401.             p: ptr;                                                {temporary pointer to image or background points}
  402.             Databand: Linetype;                           {current scan line in image}
  403.     begin
  404.         backgrpt := 0;
  405.         with Info^, BoundRect do begin
  406.                 for vloc := top to bottom - 1 do begin
  407.                         GetLine(0, vloc, pixelsperline, Databand);
  408.                         for hloc := left to right - 1 do begin
  409.                                 p := ptr(backgroundaddr + backgrpt);
  410.                                 pvalue := Databand[hloc] - BAND(p^, 255);
  411.                                 if pvalue < 0 then
  412.                                     Databand[hloc] := 0
  413.                                 else
  414.                                     Databand[hloc] := pvalue;
  415.                                 backgrpt := backgrpt + 1;
  416.                             end;     {for}
  417.                         offset := LongInt(vloc) * BytesPerRow;
  418.                         p := ptr(ord4(PicBaseAddr) + offset);
  419.                         BlockMove(@Databand, p, pixelsperline);
  420.                     end;  {for}
  421.             end;     {with}
  422.     end;
  423.  
  424.  
  425.     procedure Background2D;
  426. {******************************************************************************}
  427. {*     This procedure implements a rolling-ball algorithm for the removal of smooth continuous background       *}
  428. {*  from a two-dimensional gel image.  It rolls the ball (actually a square patch on the top of a sphere) on a       *}
  429. {*  low-resolution (by a factor of 'shrinkfactor' times) copy of the original image in order to increase speed     *}
  430. {*  with little loss in accuracy.  It uses interpolation and extrapolation to blow the shrunk image to full size.     *}
  431. {******************************************************************************}
  432.         var
  433.             tport: Grafptr;
  434.             i,                                                                                    {loop index for shrunk image memory}
  435.             backgroundsize,                                                             {size of the background memory}
  436.             smallimagesize: Size;                                                    {size of the shrunk image memory}
  437.             p: ptr;                                                                            {pointer to current pixel in shrunk image memory}
  438.     begin
  439.         ShowWatch;
  440.         UpdateMeter(0, 'Building Rolling Ball...');
  441.         GetPort(tPort);
  442.         with Info^ do begin
  443.                 SetPort(GrafPtr(osPort));
  444.                 BoundRect := roiRect;
  445.             end;
  446.         GetRollingBall;                                                                  {precompute the rolling ball}
  447.         UpdateMeter(3, 'Building Rolling Ball...');
  448.         balladdr := ord4(ballptr);
  449.         with BoundRect do begin
  450.                 leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
  451.                 rightroll := right div shrinkfactor - 1;                      {on which to roll ball}
  452.                 toproll := top div shrinkfactor;
  453.                 bottomroll := bottom div shrinkfactor - 1;
  454.                 backgroundsize := longint(right - left) * longint(bottom - top);
  455.                 backgroundptr := NewPtrToClearBuf(backgroundsize);
  456.                 Aborting := backgroundptr = nil;
  457.                 backgroundaddr := ord4(backgroundptr);
  458.                 smallimagesize := longint(rightroll - leftroll + 1) * longint(bottomroll - toproll + 1);
  459.                 smallimageptr := NewPtrToClearBuf(smallimagesize);
  460.                 Aborting := Aborting or (smallimageptr = nil);
  461.                 smallimageaddr := ord4(smallimageptr);
  462.                 if not aborting then begin
  463.                         UpdateMeter(10, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
  464.                         for i := 0 to smallimagesize - 1 do begin                {create a lower resolution image for ball-rolling}
  465.                                 p := ptr(smallimageaddr + i);
  466.                                 xxcenter := left + shrinkfactor * (i mod (rightroll - leftroll + 1));
  467.                                 yycenter := top + shrinkfactor * (i div (rightroll - leftroll + 1));
  468.                                 p^ := MinIn2DMask;                                   {each point in small image is min in its neighborhood}
  469.                             end;
  470.                         RollBall;
  471.                     end
  472.                 else
  473.                     beep;
  474.                 if not Aborting then begin
  475.                         UpdateMeter(90, 'Interpolating Background...');
  476.                         InterpolateBackground2D;                              {interpolate to find background interior}
  477.                         UpdateMeter(95, 'Extrapolating Background...');
  478.                         ExtrapolateBackground2D;                             {extrapolate on top and bottom}
  479.                         UpdateMeter(98, 'Subtracting Background...');
  480.                         SubtractBackground2D;                                  {subtract background from original image}
  481.                         UpdateMeter(100, 'Subtracting Background...');
  482.                     end;
  483.             end;   {with boundrect}
  484.         DisposPtr(backgroundptr);                           {free up background, rolling ball, shrunk image memory}
  485.         DisposPtr(ballptr);
  486.         DisposPtr(smallimageptr);
  487.         DisposeWindow(MeterWindow);
  488.         MeterWindow := nil;
  489.         SetPort(tPort);
  490.     end;
  491.  
  492.  
  493.     procedure RollArc (left, rightminusone, diam: integer; var background, ballpoints: IntRow; var Dataline: Linetype);
  494.         var
  495.             xpt, xpt2, xval, ydif, ymin, yctr, bpt, yadd: integer;
  496.     begin
  497.         for xpt := left to rightminusone do begin
  498.                 background[xpt] := -255;         {init background curve to minimum values}
  499.             end;
  500.         yctr := 0;                                   {start y-center at the x axis}
  501.         for xpt := left to (rightminusone + diam - 1) do {while semicircle is tangent to edges or within curve...}
  502.             begin                                       {xpt is far right edge of semi-circle}
  503. {do we have to move the circle?...}
  504.                 ymin := 256;                          {highest could ever be 255}
  505.                 bpt := 0;
  506.                 xpt2 := xpt - diam;                {xpt2 is far left edge of semi-circle}
  507.                 while xpt2 <= xpt do            {check every point on semicircle}
  508.                     begin
  509.                         if ((xpt2 >= left) and (xpt2 <= rightminusone)) then begin  {only examine points on curve}
  510.                                 ydif := dataline[xpt2] - (yctr + ballpoints[bpt]);                {curve minus circle points}
  511.                                 if (ydif < ymin) then begin {keep most negative, since ball should always be below curve}
  512.                                         ymin := ydif;
  513.                                     end;
  514.                             end;  {if xpt2 }
  515.                         bpt := bpt + 1;
  516.                         xpt2 := xpt2 + 1;
  517.                     end;  {while xpt2 }
  518.                 if (ymin <> 256) then{if we found a new minimum...}
  519.                     yctr := yctr + ymin;   {move circle up or down}
  520. {now compare every point on semi with background,  and keep highest number}
  521.                 xval := xpt - diam;
  522.                 xpt2 := 0;
  523.                 while xpt2 <= diam do begin  {for all the points in the semicircle}
  524.                         if ((xval >= left) and (xval <= rightminusone)) then begin
  525.                                 yadd := yctr + ballpoints[xpt2];
  526.                                 if (yadd > background[xval]) then  {keep largest adjustment}
  527.                                     background[xval] := yadd;
  528.                             end;
  529.                         xval := xval + 1;
  530.                         xpt2 := xpt2 + 1;
  531.                     end;  {while xpt2}
  532.             end;  {for xpt }
  533.     end;
  534.  
  535.  
  536.     function MinIn1DMask (var Databand: LineType; xcenter: integer): integer;
  537. {*******************************************************************************}
  538. {*     MinIn1DMask finds the minimum pixel value in a (2*shrinkfactor-1) mask about the point xcenter in the *}
  539. {*  current line.  This code must run FAST because it gets called OFTEN!                                                                   *}
  540. {*******************************************************************************}
  541.         var
  542.             i,                                                                              {index to pixels in the mask}
  543.             temp: integer;                                                          {temporary minimum value}
  544.     begin
  545.         temp := Databand[xcenter - shrinkfactor + 1];
  546.         for i := xcenter - shrinkfactor + 2 to xcenter + shrinkfactor - 1 do
  547.             if (Databand[i] < temp) then
  548.                 temp := Databand[i];
  549.         MinIn1DMask := temp;
  550.     end;
  551.  
  552.  
  553. {******************************************************************************}
  554. {*     This procedure computes the location of each point on the rolling arc relative to the center of the circle     *}
  555. {*  containing it.  The arc is located in the top half of this circle.  The vertical axis of the circle passes through  *}
  556. {*  the midpoint of the arc.  The projection of the arc on the x-axis below is a line segment.                                 *}
  557. {******************************************************************************}
  558.     procedure GetRollingArc (var arcpoints: IntRow; var arcwidth: integer);
  559.         var
  560.             xpt,                                                                                 {x-point along arc}
  561.             xval,                                                                               {x-point in arc array}
  562.             rsquare,                                                                         {shrunken arc radius squared}
  563.             xtrim,                                                                            {points to be trimmed off each end of arc}
  564.             smallballradius,                                                            {radius of shrunken arc which actually rolls}
  565.             diam: integer;                                                                 {diameter of shrunken arc's circle}
  566.     begin
  567.         smallballradius := ballradius div shrinkfactor;            { operate on small-sized image with small-sized ball}
  568.         rsquare := smallballradius * smallballradius;
  569.         for xpt := -smallballradius to smallballradius do        { find the ballpoints for arc based at  (x,y)=(0,0) }
  570.             begin
  571.                 xval := xpt + smallballradius;                                     {offset, can't have negative index}
  572.                 arcpoints[xval] := round(sqrt(abs(rsquare - (xpt * xpt))));  {Ys are positive, top half of circle}
  573.             end;
  574.         diam := smallballradius * 2;
  575.         xtrim := (ArcTrimPer * diam) div 100;                       {how many points to trim off each end}
  576.         arcwidth := diam - xtrim - xtrim;
  577.         for xpt := -smallballradius to smallballradius - xtrim - xtrim do begin
  578.                 xval := xpt + smallballradius;
  579.                 arcpoints[xval] := arcpoints[xval + xtrim];
  580.             end;
  581.         for xpt := smallballradius - xtrim - xtrim + 1 to smallballradius do begin
  582.                 xval := xpt + smallballradius;
  583.                 arcpoints[xval] := 0;
  584.             end;
  585.     end;
  586.  
  587.  
  588.     procedure ExtrapolateBackground1D (var Backline, Dataline: LineType; background: IntRow; leftroll, rightroll: integer);
  589. {******************************************************************************}
  590. {*     This procedure uses linear extrapolation to find pixel values on the left and right edges of the current        *}
  591. {*  line of the background.  It finds the edge points by extrapolating from the edges of the calculated and               *}
  592. {*  interpolated background interior.  If extrapolation yields values below zero or above 255, then they are set *}
  593. {*  to zero and 255 respectively.                                                                                                                               *}
  594. {******************************************************************************}
  595.         var
  596.             i,                                                                             {index to edges of background array}
  597.             hloc,                                                                       {}
  598.             edgeslope: integer;                                                 {}
  599.     begin
  600.         with BoundRect do begin
  601.                 edgeslope := (background[leftroll + 1] - background[leftroll + 2]) div shrinkfactor;
  602.                 for i := left to shrinkfactor * (leftroll + 1) - 1 do begin     {extrapolate on left edge}
  603.                         hloc := shrinkfactor * (leftroll + 1) - 1 + left - i;
  604.                         if (Backline[hloc + 1] + edgeslope < 0) then
  605.                             Backline[hloc] := 0
  606.                         else if (Backline[hloc + 1] + edgeslope > Dataline[hloc]) then
  607.                             Backline[hloc] := Dataline[hloc]
  608.                         else
  609.                             Backline[hloc] := Backline[hloc + 1] + edgeslope;
  610.                     end;
  611.                 edgeslope := (background[rightroll] - background[rightroll - 1]) div shrinkfactor;
  612.                 for hloc := shrinkfactor * rightroll to right - 1 do begin     {extrapolate on right edge}
  613.                         if (Backline[hloc - 1] + edgeslope < 0) then
  614.                             Backline[hloc] := 0
  615.                         else if (Backline[hloc - 1] + edgeslope > Dataline[hloc]) then
  616.                             Backline[hloc] := Dataline[hloc]
  617.                         else
  618.                             Backline[hloc] := Backline[hloc - 1] + edgeslope;
  619.                     end;
  620.             end;     {with}
  621.     end;
  622.  
  623.     procedure Background1D;
  624.         var
  625.             tport: Grafptr;
  626.             hloc, vloc, arcwidth, leftroll, rightroll, numpixels: integer;
  627.             left, right, top, bottom: integer;                      {image bounds; ROTATED if RollingVerticalArc}
  628.             i, j, maskwidth: integer;
  629.             background, arcpoints: IntRow;
  630.             offset: LongInt;
  631.             p: ptr;
  632.             Dataline, Backline, Smalldataline: Linetype;
  633.             str: str255;
  634.     begin
  635.         ShowWatch;
  636.         UpdateMeter(0, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
  637.         GetPort(tPort);
  638.         with Info^ do begin
  639.                 SetPort(GrafPtr(osPort));
  640.                 BoundRect := roiRect;
  641.             end;
  642.         GetRollingArc(arcpoints, arcwidth);
  643.         maskwidth := shrinkfactor + shrinkfactor - 1;
  644.         case BackSubKind of
  645.             RollingHorizontalArc:  begin
  646.                     left := BoundRect.left;
  647.                     top := BoundRect.top;
  648.                     right := BoundRect.right;
  649.                     bottom := BoundRect.bottom;
  650.                     numpixels := Info^.pixelsperline;
  651.                     str := 'Rolling Disk Horizontally...';
  652.                 end;
  653.             RollingVerticalArc:  begin
  654.                     left := BoundRect.top;
  655.                     top := BoundRect.left;
  656.                     right := BoundRect.bottom;
  657.                     bottom := BoundRect.right;
  658.                     numpixels := Info^.nlines;
  659.                     str := 'Rolling Disk Vertically...';
  660.                 end;
  661.         end;     {case}
  662.         leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
  663.         rightroll := right div shrinkfactor - 1;                      {on which to roll arc}
  664.         for vloc := top to bottom - 1 do begin  {for ROI}
  665.                 case BackSubKind of
  666.                     RollingHorizontalArc: 
  667.                         GetLine(0, vloc, numpixels, Dataline);
  668.                     RollingVerticalArc: 
  669.                         GetColumn(vloc, 0, numpixels, Dataline);
  670.                 end;     {case}
  671.                 for i := leftroll + 1 to rightroll do
  672.                     smalldataline[i] := MinIn1DMask(Dataline, shrinkfactor * i - 1);
  673.                 RollArc(leftroll + 1, rightroll, arcwidth, background, arcpoints, smalldataline);  {roll arc on one line}
  674.                 for i := leftroll + 1 to rightroll do begin           {interpolate to find interior background points}
  675.                         hloc := shrinkfactor * i - 1;
  676.                         Backline[hloc] := background[i];
  677.                         for j := 1 to shrinkfactor - 1 do
  678.                             Backline[hloc - j] := background[i - 1] + (shrinkfactor - j) * (background[i] - background[i - 1]) div shrinkfactor;
  679.                     end;
  680.                 ExtrapolateBackground1D(Backline, Dataline, background, leftroll, rightroll);
  681.                 for i := left to right - 1 do                                  {subtract background from current scan line}
  682.                     Dataline[i] := Dataline[i] - Backline[i];
  683.                 case BackSubKind of
  684.                     RollingHorizontalArc: 
  685.                         with Info^ do begin
  686.                                 offset := LongInt(vloc) * BytesPerRow;
  687.                                 p := ptr(ord4(PicBaseAddr) + offset);
  688.                                 BlockMove(@Dataline, p, numpixels);            {fast whole line write}
  689.                             end;
  690.                     RollingVerticalArc: 
  691.                         PutColumn(vloc, 0, numpixels, Dataline);         {slow whole column write}
  692.                 end;     {case}
  693.                 if ((vloc mod 8) = 0) and (vloc > 16) then begin
  694.                         UpdateMeter((LongInt(vloc - top) * 100) div (bottom - top - 1), str);
  695.                         if CommandPeriod then begin
  696.                                 beep;
  697.                                 Aborting := true;
  698.                                 leave;
  699.                             end;
  700.                     end;
  701.             end;
  702.         UpdateMeter(100, str);
  703.         DisposeWindow(MeterWindow);
  704.         MeterWindow := nil;
  705.         SetPort(tPort);
  706.     end;
  707.  
  708.     procedure SetUpGel;
  709.         var
  710.             frame: rect;
  711.             AutoSelectAll: boolean;
  712.             p: ptr;
  713.     begin
  714.         if NotinBounds or NotRectangular then
  715.             exit(SetUpGel);
  716.         StopDigitizing;
  717.         AutoSelectAll := not Info^.RoiShowing;
  718.         if AutoSelectAll then
  719.             SelectAll(false);
  720.         SetupUndoFromClip;
  721.         with info^ do begin
  722.                 frame := roiRect;
  723.                 if ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction) then
  724.                     ApplyLookupTable;
  725.                 changes := true;
  726.             end;
  727.         case BackSubKind of
  728.             RollingHorizontalArc, RollingVerticalArc: 
  729.                 Background1D;    {--------------> call background subtract <-------------------}
  730.             RollingBall: 
  731.                 Background2D;
  732.             RollingBothArcs:  begin
  733.                     BackSubKind := RollingHorizontalArc;           {remove horizontal streaks}
  734.                     Background1D;
  735.                     BackSubKind := RollingVerticalArc;               {remove vertical streaks}
  736.                     if not aborting then
  737.                         Background1D;
  738.                     BackSubKind := RollingBothArcs;                   {leave BackSubKind as we found it}
  739.                 end;
  740.         end;     {case}
  741.         UpdatePicWindow;
  742.         SetUpRoiRect;
  743.         WhatToUndo := UndoFilter;
  744.         Info^.changes := true;
  745.         ShowWatch;
  746.         if AutoSelectAll then
  747.             KillRoi;
  748.         if Aborting then begin
  749.                 Undo;
  750.                 WhatToUndo := NothingToUndo;
  751.                 UpdatePicWindow;
  752.             end;
  753.     end;
  754.  
  755.  
  756.     procedure GetBallRadius;
  757.         var
  758.             SaveRadius: integer;
  759.             canceled: boolean;
  760.     begin
  761.         SaveRadius := BallRadius;
  762.         BallRadius := GetInt('Rolling BallRadius:', BallRadius, canceled);
  763.         if (BallRadius < 5) or (BallRadius > 319) or canceled then begin
  764.                 BallRadius := SaveRadius;
  765.                 if not canceled then
  766.                     beep;
  767.             end;
  768.     end;
  769.  
  770.  
  771.     procedure DoBackgroundMenuEvent (MenuItem: integer);
  772.         var
  773.             map_array: Ptr;
  774.     begin
  775.         MeterWindow := nil;
  776.         Aborting := false;
  777.         case MenuItem of
  778.             HorizontalItem, VerticalItem:  begin
  779.                     if FasterBackgroundSubtraction then begin
  780.                             ArcTrimPer := 20;
  781.                             shrinkfactor := 4;
  782.                         end
  783.                     else begin
  784.                             ArcTrimPer := 10;
  785.                             shrinkfactor := 2;
  786.                         end;
  787.                     if MenuItem = HorizontalItem then
  788.                         BackSubKind := RollingHorizontalArc
  789.                     else
  790.                         BackSubKind := RollingVerticalArc;
  791.                     SetUpGel;
  792.                 end;
  793.             Sub2DItem:  begin
  794.                     if FasterBackgroundSubtraction then begin
  795.                             ArcTrimPer := 20;     {trim 40% in x and y}
  796.                             shrinkfactor := 8;
  797.                         end
  798.                     else begin
  799.                             ArcTrimPer := 16;     {trim 32% in x and y}
  800.                             shrinkfactor := 4;
  801.                         end;
  802.                     BackSubKind := RollingBall;
  803.                     SetUpGel;
  804.                 end;
  805.             RemoveStreaksItem:  begin
  806.                     ArcTrimPer := 20;
  807.                     shrinkfactor := 4;
  808.                     BackSubKind := RollingBothArcs;
  809.                     SetUpGel;
  810.                 end;
  811.             FasterItem: 
  812.                 FasterBackgroundSubtraction := not FasterBackgroundSubtraction;
  813.             RadiusItem: 
  814.                 GetBallRadius;
  815.         end; {case}
  816.     end;
  817.  
  818.  
  819. end.