home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / aijournl / ai_may89.arc / AIEDGE2.PAS < prev    next >
Pascal/Delphi Source File  |  1988-10-04  |  56KB  |  1,717 lines
  1. unit AiEdge2;
  2.  
  3. Interface
  4.  
  5. uses aiglob,bordunit,aimath,aiuser;
  6.  
  7. Var ProAddr : pointer;
  8.  
  9.  
  10. Function IsItForeground(x,y:word;backgroundvalue:byte;distance:integer):boolean;
  11. Function IsItBackground(x,y:word;backgroundvalue:byte;distance:integer):boolean;
  12. Function IsItForegroundv(x,y:word;backgroundvalue:byte;distance:integer):boolean;
  13. Function IsItBackgroundv(x,y:word;backgroundvalue:byte;distance:integer):boolean;
  14. Procedure Erosion(x1,y1,x2,y2 : word);
  15. Procedure Erosion2(x1,y1,x2,y2:word);
  16.  
  17. Function IntensityCheck(x,y:word;size:byte):boolean;
  18. Function SpotContrast(x,y:word;nucsize:byte;Var goodifsmall:boolean):boolean;
  19. Function ScanEdge(x1,y1,x2,y2:word):word;
  20. Function FindArea(x1,y1,x2,y2:word;Var _fore,_std:double):word;
  21. Procedure MakeDark(x1,y1,x2,y2:word);
  22.  
  23. Procedure SetAddress;
  24. Procedure FillIn(x1,y1,x2,y2:word;small:boolean;
  25.          nucsize:byte;observed:byte);
  26. Procedure HowMuchFore(x,y:word;size:byte;Var AmtFore,_stdev:double);
  27. Procedure EScan(x,y:word;Nucsize,Cv:byte;
  28.             Var da,db:byte);
  29. Function Mscan(x,y:word;size:byte;Var bstuff:double):byte;
  30.  
  31. Procedure LearnFromDeletion(Num:byte);
  32. Procedure LearnFromAddition(x,y:word;Nucsize,Width,Height:byte;Mval:double);
  33. Procedure HistoAnalysis(x,y:word;nucsize:byte;Var below:byte;
  34.     var ku,stout,rx,rx2:double;Var CytCond,Abscyt : boolean);
  35. Function ShellScan(x,y,nucsize:word;Uncertain:boolean;
  36.                     Var goodifnucleolus:boolean):boolean;
  37.  
  38. Implementation
  39.  
  40. {These two routines are used by the four routines immediately below them.}
  41.  
  42. function digit(x,y:word):byte;
  43. begin
  44.   digit := oldgrayvalue(x,y);
  45. end;
  46.  
  47. Procedure Setaddress;
  48. begin
  49.  proaddr := @digit;
  50. end;
  51.  
  52. {The following functions, given (x,y), scan DISTANCE pixels to the right or
  53.  left (depending on the sign) in order to determine how many consecutive
  54.  pixels are above or below the backgroundvalue.  Two functions scan
  55.  horizontally, the other two scan vertically.}
  56.  
  57. function isitbackground(x,y:word;backgroundvalue:byte;distance:integer):boolean;
  58. begin
  59. inline($8b/$4e/<distance/        {mov cx,distance}
  60.        $83/$f9/$00/              {cmp cx,00}
  61.        $74/$24/                  {jz dgd}
  62. {again}$8b/$46/<x/               {mov ax,x}
  63.        $03/$c1/                  {add ax,cx}
  64.        $51/                      {push cx}
  65.        $50/                      {push ax}
  66.        $ff/$76/<y/               {push y}
  67.        $ff/$16/proaddr/          {call digit}
  68.        $59/                      {pop cx}
  69.        $3a/$46/<backgroundvalue/ {cmp al,backgroundvalue}
  70.        $72/$05/                  {jb dchk}
  71. {nogd} $b3/$00/                  {mov bl,00}
  72.        $eb/$0e/                  {jmp done}
  73.        $90/
  74. {dchk} $83/$f9/$00/              {cmp cx,+00}
  75.        $7f/$03/                  {ja pos}
  76.  {neg} $41/                      {inc cx}
  77.        $eb/$da/                  {jmp again}
  78.  {pos} $49/                      {dec cx}
  79.        $eb/$d7/                  {jmp again}
  80.  {dgd} $b3/$01/                  {mov bl,01}
  81.  {done}$88/$5e/$ff);             {mov [bp-01],bl}
  82. end;
  83.  
  84. function isitforeground(x,y:word;backgroundvalue:byte;distance:integer):boolean;
  85. begin
  86. inline($8b/$4e/<distance/        {mov cx,distance}
  87.        $83/$f9/$00/              {cmp cx,00}
  88.        $74/$24/                  {jz dgd}
  89. {again}$8b/$46/<x/               {mov ax,x}
  90.        $03/$c1/                  {add ax,cx}
  91.        $51/                      {push cx}
  92.        $50/                      {push ax}
  93.        $ff/$76/<y/               {push y}
  94.        $ff/$16/proaddr/               {call digit}
  95.        $59/                      {pop cx}
  96.        $3a/$46/<backgroundvalue/ {cmp al,backgroundvalue}
  97.        $77/$05/                  {ja dchk}
  98. {nogd} $b3/$00/                  {mov bl,00}
  99.        $eb/$0e/                  {jmp done}
  100.        $90/
  101. {dchk} $83/$f9/$00/              {cmp cx,+00}
  102.        $7f/$03/                  {ja pos}
  103.  {neg} $41/                      {inc cx}
  104.        $eb/$da/                  {jmp again}
  105.  {pos} $49/                      {dec cx}
  106.        $eb/$d7/                  {jmp again}
  107.  {dgd} $b3/$01/                  {mov bl,01}
  108.  {done}$88/$5e/$ff);             {mov [bp-01],bl}
  109. end;
  110.  
  111.  
  112. function isitbackgroundv(x,y:word;backgroundvalue:byte;distance:integer):boolean;
  113. begin
  114. inline($8b/$4e/<distance/        {mov cx,distance}
  115.        $83/$f9/$00/              {cmp cx,00}
  116.        $74/$24/                  {jz dgd}
  117. {again}$8b/$46/<y/               {mov ax,y}
  118.        $03/$c1/                  {add ax,cx}
  119.        $51/                      {push cx}
  120.        $ff/$76/<x/               {push x}
  121.        $50/                      {push ax}
  122.        $ff/$16/proaddr/          {call digit}
  123.        $59/                      {pop cx}
  124.        $3a/$46/<backgroundvalue/ {cmp al,backgroundvalue}
  125.        $72/$05/                  {jb dchk}
  126. {nogd} $b3/$00/                  {mov bl,00}
  127.        $eb/$0e/                  {jmp done}
  128.        $90/
  129. {dchk} $83/$f9/$00/              {cmp cx,+00}
  130.        $7f/$03/                  {ja pos}
  131.  {neg} $41/                      {inc cx}
  132.        $eb/$da/                  {jmp again}
  133.  {pos} $49/                      {dec cx}
  134.        $eb/$d7/                  {jmp again}
  135.  {dgd} $b3/$01/                  {mov bl,01}
  136.  {done}$88/$5e/$ff);             {mov [bp-01],bl}
  137. end;
  138.  
  139. function isitforegroundv(x,y:word;backgroundvalue:byte;distance:integer):boolean;
  140. begin
  141. inline($8b/$4e/<distance/        {mov cx,distance}
  142.        $83/$f9/$00/              {cmp cx,00}
  143.        $74/$24/                  {jz dgd}
  144. {again}$8b/$46/<y/               {mov ax,y}
  145.        $03/$c1/                  {add ax,cx}
  146.        $51/                      {push cx}
  147.        $ff/$76/<x/               {push x}
  148.        $50/                      {push ax}
  149.        $ff/$16/proaddr/               {call digit}
  150.        $59/                      {pop cx}
  151.        $3a/$46/<backgroundvalue/ {cmp al,backgroundvalue}
  152.        $77/$05/                  {ja dchk}
  153. {nogd} $b3/$00/                  {mov bl,00}
  154.        $eb/$0e/                  {jmp done}
  155.        $90/
  156. {dchk} $83/$f9/$00/              {cmp cx,+00}
  157.        $7f/$03/                  {ja pos}
  158.  {neg} $41/                      {inc cx}
  159.        $eb/$da/                  {jmp again}
  160.  {pos} $49/                      {dec cx}
  161.        $eb/$d7/                  {jmp again}
  162.  {dgd} $b3/$01/                  {mov bl,01}
  163.  {done}$88/$5e/$ff);             {mov [bp-01],bl}
  164. end;
  165.  
  166. {Find lone RED pixels}
  167.  
  168. Function Erode1(x,y : word):boolean;
  169. begin
  170.   Erode1 := TRUE;
  171.   If (oldgrayvalue(x-1,y-1) and 1 = 1) or (oldgrayvalue(x,y-1) and 1 = 1) or
  172.      (oldgrayvalue(x+1,y-1) and 1 = 1) or (oldgrayvalue(x-1,y) and 1 = 1) or
  173.      (oldgrayvalue(x+1,y) and 1 = 1) or (oldgrayvalue(x-1,y+1) and 1 = 1) or
  174.      (oldgrayvalue(x,y+1) and 1 = 1) or (oldgrayvalue(x+1,y+1) and 1 = 1) then
  175.      Erode1 := FALSE;
  176. end; {end function erode1}
  177.  
  178.  
  179. {Erase single red dots}
  180.  
  181. Procedure Erosion(x1,y1,x2,y2 : word);
  182. Var
  183.     j,k : word;
  184. begin
  185.    newgrayvalue(1,1,oldgrayvalue(1,1));
  186.    For k := y1 to y2 do
  187.      for j := x1 to x2 do
  188.        If (oldgrayvalue(j,k) and 1 = 1) and Erode1(j,k) then
  189.      newgrayvalue(j,k,(oldgrayvalue(j,k) and $FE));
  190. end; {end procedure erosion}
  191.  
  192. {Erase all red dots}
  193.  
  194. Procedure Erosion2(x1,y1,x2,y2 : word);
  195. Var
  196.     j,k : word;
  197. begin
  198.    newgrayvalue(1,1,oldgrayvalue(1,1));
  199.    For k := y1 to y2 do
  200.      for j := x1 to x2 do
  201.        If (oldgrayvalue(j,k) and 1 = 1) then
  202.      newgrayvalue(j,k,(oldgrayvalue(j,k) and $FE));
  203. end; {end procedure erosion}
  204.  
  205. {This function scans within the region defined by (x1,y1,x2,y2) and
  206.  counts the number of RED marks to calculate the area.  In addition,
  207.  the routine calculates the average gray level and standard deviation
  208.  of the shaded region.}
  209.  
  210. Function FindArea(x1,y1,x2,y2:word;Var _fore,_std:double):word;
  211. Var j,k:word;
  212.     area : word;
  213.     gray1 : byte;
  214.     count : word;
  215.     imagdata : imagtype2;
  216. begin
  217.   area := 0;
  218.   count := 0;
  219.   For k := y1 to y2 do                {scan within box}
  220.     for j := x1 to x2 do
  221.     begin
  222.       gray1 := oldgrayvalue(j,k);     {get value}
  223.       If gray1 and 1 = 1 then         {is it RED?}
  224.       begin
  225.     area := area + 1;             {increment area count}
  226.     If gray1 > lowdiv then        {is it part of nucleolus?}
  227.     begin
  228.       count := count+1;           {If not then use value to calculate   }
  229.       imagdata[count] := gray1;   {Mean and StDev.  This helps to focus }
  230.     end;                          {only values describing the nucleus.  }
  231.       end;
  232.   end;
  233.   _fore := Mean(imagdata,1,count);
  234.   _std := stdev(imagdata,1,count,_fore);
  235.   FindArea := area;
  236. end;{end function findarea}
  237.  
  238. {In this procedure we fill in the object by alternating between different
  239.  erosion and dilating techniques.
  240.  SMALL describes whether the main program thinks it is a big nucleus or
  241.  a small one, and NUCSIZE is a value bigger than the largest Nucleus
  242.  diameter and is used as a maximum scanning distance.}
  243.  
  244. Procedure FillIn(x1,y1,x2,y2:word;small:boolean;
  245.          nucsize:byte;observed:byte);
  246. Var
  247.     j,k      : word;                      {general x,y counters}
  248.     gray1    : byte;                      {gray level value}
  249.     xa,ya    : word;                      {center coord}
  250.     i,q,r    : word;                      {common variables}
  251.     f        : integer;                   {used in erase routine}
  252.     Highest,                              {high gray value}
  253.     lowest   : byte;                      {low gray value}
  254.     debug    : boolean;
  255.     redcount :word;                       {used when counting red dots}
  256.     leg      : word;                      {largest diagonal from center
  257.                          to corner}
  258.     imagdata : imagtype2;                 {used to find backgnd}
  259.     count    : byte;
  260.     _mean,_f : double;
  261.     EraseMode: boolean;                   {used in erase routine}
  262.     diagdist : byte;                      {used in erase routine}
  263.     halfnuc,                              {size parameter of cell}
  264.     hnuc     : byte;
  265.     xhigh,                                {coordinates of brightest pixel}
  266.     yhigh    : word;
  267.     xpart,ypart,                          {width and height variables}
  268.     backgnd,                              {background threshold}
  269.     lowcount : byte;                      {amount of nucleolus}
  270.     adjust,
  271.     obs_adjust   : double;
  272.  
  273. begin
  274.   {..............................scan for values....................}
  275.  
  276.     nucsize := round(1.1*nucsize);                  {get a larger value     }
  277.     xa := (x1+x2) shr 1;                            {get center and diagonal}
  278.     ya := (y1+y2) shr 1;                            {   to the corner       }
  279.     leg := max(abs(xa-x1),abs(ya-y1));
  280.     leg := round (sqrt( sqr(leg+1) + sqr(leg+1) ));
  281.     diagdist := nucsize;
  282.     debug := true;
  283.     halfnuc := nucsize shr 1;                       {size up other variables}
  284.     hnuc := round(halfnuc/2);
  285.     If hnuc = 1 then
  286.        hnuc := 2;
  287.  
  288.     highest := 0;                                  {find highest and lowest}
  289.     lowest  := 255;                                {values within a sampled}
  290.     count := 0;                                    {region as well as coords}
  291.     lowcount := 0;                                 {of highest pixel value}
  292.     for k := ya-(halfnuc shr 1) to ya+(halfnuc shr 1) do
  293.     for j := xa-(halfnuc) to xa+(halfnuc) do
  294.     begin
  295.       gray1 := oldgrayvalue(j,k);
  296.       count := count+1;
  297.       imagdata[count] := gray1;
  298.       If gray1 < lowdiv then
  299.         lowcount := lowcount+1;
  300.       if (gray1 > highest) then
  301.       begin
  302.         highest := gray1;
  303.         xhigh := j;
  304.         yhigh := k;
  305.       end
  306.       else if gray1 < lowest then
  307.         lowest := gray1;
  308.     end;
  309.  
  310.     highest := 0;                                  {move to highest region  }
  311.     lowest  := 255;                                {and scan again for high }
  312.     count := 0;                                    {and low values.         }
  313.     for k := yhigh-(hnuc shr 1) to yhigh+(hnuc shr 1) do
  314.     for j := xhigh-(hnuc shr 1) to xhigh+(hnuc shr 1) do
  315.     begin
  316.       gray1 := oldgrayvalue(j,k);
  317.           If gray1 > lowdiv then
  318.           begin
  319.         count := count+1;
  320.         imagdata[count] := gray1;
  321.           end;
  322.       if (gray1 > highest) then
  323.         highest := gray1
  324.       else if gray1 < lowest then
  325.         lowest := gray1;
  326.     end;
  327.     _mean := mean(imagdata,1,count);               {compute a Mean gray level}
  328.  
  329. {...........................Determine background threshold.................}
  330.     adjust := 0;
  331.     obs_adjust := 1;
  332.  
  333.    If previous then
  334.      backgnd := round( (0.85 + (observed*0.01))*_mean)
  335.    else
  336.      backgnd := round(0.85*_mean);
  337.    lowdiv := 60;
  338.  
  339.    nucsize := round(nucsize/1.1);               {reset nucsize}
  340.    xpart := round(0.95*nucsize);                {Make width shorter and }
  341.    ypart := round(1.3*nucsize);                 {height longer since the }
  342.                  {spot-scanner will probably start finding values}
  343.                  {at the top of the nucleus.                     }
  344.  
  345. {........................................pass1...............................}
  346.     for k := ya-nucsize to ya+ypart do           {scan horizontally}
  347.       for j := xa-xpart to xa+xpart do           {If pixel is within bounds }
  348.       begin                                      {and is in line with 3 other}
  349.        gray1 := oldgrayvalue(j,k);               {pixels above backgnd value }
  350.        if  (gray1 < 1.005*highest) and           {then shade RED (OR low bit)}
  351.        (isitforeground(j,k,backgnd,4) or isitforeground(j,k,backgnd,-4))
  352.     then
  353.        newgrayvalue(j,k,oldgrayvalue(j,k) or 1);
  354.        end;
  355. {..............................pass2...................................}
  356.     for j := xa-xpart to xa+xpart do             {scan vertically}
  357.       for k := ya-nucsize to ya+ypart do
  358.       begin
  359.     gray1 := oldgrayvalue(j,k);
  360.     if (gray1 < 1.005*highest) and
  361.        (isitforegroundv(j,k,backgnd,4)  or isitforegroundv(j,k,backgnd,-4))
  362.     or (gray1 < lowdiv) then
  363.        newgrayvalue(j,k,gray1 or 1);
  364.       end;
  365.  
  366. {-------------------------filter little stuff-------------------------------}
  367.  
  368. for k := ya-ypart to ya+ypart do
  369.   for j := xa-xpart to xa+xpart do
  370.     if (oldgrayvalue(j,k) and 1 = 1) then            {matrix 3x3}
  371.     begin
  372.       q := 0;
  373.       for f := k-1 to k+1 do
  374.     for r := j-1 to j+1 do
  375.       if (oldgrayvalue(r,f) and 1 = 1) then
  376.         q := q + 1;
  377.       if q < 5 then
  378.     newgrayvalue(j,k,oldgrayvalue(j,k) and $FE);
  379.     end;
  380.             {determine if shaded region after first
  381.              filtering is within size limits          }
  382.  
  383. q := findarea(x1,y1,x2,y2,_f,_f);
  384.  
  385. If (q > minarea) and (q < 1.5*maxarea) then   {if area in limits}
  386. begin
  387. {-----------------------rebuild inside along both axis----------------------}
  388.    newgrayvalue(1,1,1);
  389.     for k := ya-nucsize to ya+nucsize do             {rebuild}
  390.     begin
  391.      for j := xa to xa+nucsize do
  392.        if (oldgrayvalue(j+1,k) and 1 = 1) then
  393.      newgrayvalue(j,k,(oldgrayvalue(j,k) or 1));
  394.      j := xa;
  395.      while (j >= xa-nucsize) do
  396.      begin
  397.        if (oldgrayvalue(j-1,k) and 1 = 1) then
  398.      newgrayvalue(j,k,(oldgrayvalue(j,k) or 1));
  399.        j := j - 1;
  400.      end;
  401.     end;{for k}
  402.  
  403.     newgrayvalue(1,1,1);
  404.  
  405.     for j := xa-nucsize to xa+nucsize do             {rebuild}
  406.     begin
  407.      for k := ya to ya+nucsize do
  408.        if (oldgrayvalue(j,k+1) and 1 = 1) then
  409.      newgrayvalue(j,k,(oldgrayvalue(j,k) or 1));
  410.      k := ya;
  411.      while (k >= ya-nucsize) do
  412.      begin
  413.        if (oldgrayvalue(j,k-1) and 1 = 1) then
  414.      newgrayvalue(j,k,(oldgrayvalue(j,k) or 1));
  415.        k := k - 1;
  416.      end;
  417.     end;{for j}
  418.  
  419. {These filters scan the UPPER and LOWER RIGHT and LEFT QUADRANTS.  The
  420.  filter starts scaning from the center of the box.  If a nucleus exists then
  421.  a round region should be shaded in the center. The cytoplasm, which is
  422.  darker, should not be shaded except for some lightly stained regions.  The
  423.  region between bordering cells may also be shaded because it is lighter.
  424.  The routine scans line by line outward from the center, counting the number
  425.  of unshaded regions.  If the gap is large enough then all pixels beyond that
  426.  point on the same line are erased.  In theory, this will erase
  427.  everything outside of the shaded nucleus.}
  428.  
  429.  
  430. {------------------filter regions not connected to center region------------}
  431.  
  432. {first four scan for HORIZONTAL gaps}
  433.  
  434. k := ya;                             {erase nocontinuos segments}
  435. while (k > ya-nucsize-2) do          {scan from center up}
  436. begin
  437.   redcount := 0;                     {no RED found yet}
  438.   EraseMode := FALSE;                {do not erase pixels yet}
  439.   for j := xa to xa+nucsize do       {scan from center to right, making this}
  440.   begin                              {an UPPER RIGHT QUADRANT scan.}
  441.     if EraseMode then                {If erase mode is set then erase RED}
  442.        NewGrayvalue(j,k,Oldgrayvalue(j,k) and $FE)
  443.     else if oldgrayvalue(j,k) and 1 <> 1 then  {else if NOT RED then up count}
  444.        redcount := redcount+1;
  445.     If Redcount > 3 then             {If less than three REDs have been found}
  446.        EraseMode := TRUE;            {we must be in a gap so start erasing}
  447.   end;
  448.   k := k-1;                          {move up one line}
  449. end;
  450. k := ya;                             {reset to vertical center}
  451. while (k > ya-nucsize-2) do          {scan UPPER LEFT QUADRANT}
  452. begin
  453.   redcount := 0;
  454.   J := xa;
  455.   erasemode := false;
  456.   while (j > xa-nucsize) do
  457.   begin
  458.     if EraseMode then
  459.        NewGrayvalue(j,k,Oldgrayvalue(j,k) and $FE)
  460.     else if oldgrayvalue(j,k) and 1 <> 1 then
  461.        redcount := redcount+1;
  462.     If Redcount > 3 then
  463.        EraseMode := TRUE;
  464.     j := j-1;
  465.   end;
  466.   k := k-1;
  467. end;
  468. for k := ya to ya+nucsize+2 do       {scan LOWER QUADRANTS}
  469. begin
  470.   redcount := 0;
  471.   erasemode := false;
  472.   for j := xa to xa+nucsize do
  473.   begin
  474.     If EraseMode then
  475.       NewGrayvalue(j,k,oldgrayvalue(j,k) and $FE)
  476.     else if oldgrayvalue(j,k) and 1 <> 1 then
  477.       redcount := redcount+1;
  478.     If Redcount > 3 then
  479.        EraseMode := TRUE;
  480.   end;
  481. end;
  482. for k := ya to ya+nucsize+2 do
  483. begin
  484.   redcount := 0;
  485.   erasemode := false;
  486.   j := xa;
  487.   while (j > xa-nucsize) do
  488.   begin
  489.     If EraseMode then
  490.       NewGrayvalue(j,k,oldgrayvalue(j,k) and $FE)
  491.     else if oldgrayvalue(j,k) and 1 <> 1 then
  492.       redcount := redcount+1;
  493.     If Redcount > 3 then
  494.        EraseMode := TRUE;
  495.     j := j-1;
  496.   end;
  497. end;
  498.  
  499. {these four scan for VERTICAL gaps}
  500.  
  501. j := xa;
  502. while (j > xa-nucsize-2) do
  503. begin
  504.   redcount := 0;
  505.   EraseMode := FALSE;
  506.   for k := ya to ya+nucsize do
  507.   begin
  508.     if EraseMode then
  509.        NewGrayvalue(j,k,Oldgrayvalue(j,k) and $FE)
  510.     else if oldgrayvalue(j,k) and 1 <> 1 then
  511.        redcount := redcount+1;
  512.     If Redcount > 3 then
  513.        EraseMode := TRUE;
  514.   end;
  515.   j := j-1;
  516. end;
  517. j := xa;
  518. while (j > xa-nucsize-2) do
  519. begin
  520.   redcount := 0;
  521.   k := ya;
  522.   erasemode := false;
  523.   while (k > ya-nucsize) do
  524.   begin
  525.     if EraseMode then
  526.        NewGrayvalue(j,k,Oldgrayvalue(j,k) and $FE)
  527.     else if oldgrayvalue(j,k) and 1 <> 1 then
  528.        redcount := redcount+1;
  529.     If Redcount > 3 then
  530.        EraseMode := TRUE;
  531.     k := k-1;
  532.   end;
  533.   j := j-1;
  534. end;
  535. for j := xa to xa+nucsize+2 do
  536. begin
  537.   redcount := 0;
  538.   erasemode := false;
  539.   for k := ya to ya+nucsize do
  540.   begin
  541.     If EraseMode then
  542.       NewGrayvalue(j,k,oldgrayvalue(j,k) and $FE)
  543.     else if oldgrayvalue(j,k) and 1 <> 1 then
  544.       redcount := redcount+1;
  545.     If Redcount > 3 then
  546.        EraseMode := TRUE;
  547.   end;
  548. end;
  549. for j := xa to xa+nucsize+2 do
  550. begin
  551.   redcount := 0;
  552.   erasemode := false;
  553.   k := ya;
  554.   while (k > ya-nucsize) do
  555.   begin
  556.     If EraseMode then
  557.       NewGrayvalue(j,k,oldgrayvalue(j,k) and $FE)
  558.     else if oldgrayvalue(j,k) and 1 <> 1 then
  559.       redcount := redcount+1;
  560.     If Redcount > 3 then
  561.        EraseMode := TRUE;
  562.     k := k-1;
  563.   end;
  564. end;
  565.  
  566. {After the filtering above some unfiltered regions may still exits.  This
  567.  pass filter attempts to remove these regions.  The filter is basically
  568.  the same except now there are only UPPER and LOWER QUADRANTS.  A whole
  569.  line is scanned.  If there are not enough red pixels on that line then
  570.  all parallel lines above are erased.}
  571.  
  572. {---------filter again: erase segments not fully connected to center------}
  573.  
  574. k := ya;                             {erase nocontinuos segments}
  575. EraseMode := FALSE;                  {first verticals}
  576. nucsize := round(nucsize*1.5);
  577. while (k > ya-nucsize-2) do          {scan up from center}
  578. begin
  579.   redcount := 0;
  580.   for j := xa-nucsize to xa+nucsize do     {scan entire horizontal line}
  581.     if EraseMode then
  582.        NewGrayvalue(j,k,Oldgrayvalue(j,k) and $FE)
  583.     else if oldgrayvalue(j,k) and 1 = 1 then
  584.        redcount := redcount+1;
  585.   If Redcount <= 3 then                     {If less than four REDs then }
  586.     EraseMode := TRUE;                      {erase all lines parallel.   }
  587.   k := k-1;
  588. end;
  589. EraseMode := FALSE;
  590. for k := ya to ya+nucsize+2 do
  591. begin
  592.   redcount := 0;
  593.   for j := xa-nucsize to xa+nucsize do
  594.     If EraseMode then
  595.       NewGrayvalue(j,k,oldgrayvalue(j,k) and $FE)
  596.     else if oldgrayvalue(j,k) and 1 = 1 then
  597.       redcount := redcount+1;
  598.     If Redcount <= 3 then
  599.        EraseMode := TRUE;
  600. end;
  601. j := xa;                             {now horizontals}
  602. EraseMode := FALSE;
  603. while (j > xa-nucsize-2) do
  604. begin
  605.   redcount := 0;
  606.   for k := ya-nucsize to ya+nucsize do
  607.     if EraseMode then
  608.        NewGrayvalue(j,k,Oldgrayvalue(j,k) and $FE)
  609.     else if oldgrayvalue(j,k) and 1 = 1 then
  610.        redcount := redcount+1;
  611.   If RedCount <= 3 then
  612.     EraseMode := TRUE;
  613.   j := j-1;
  614. end;
  615. EraseMode := FALSE;
  616. for j := xa to xa+nucsize+2 do
  617. begin
  618.   redcount := 0;
  619.   for k := ya-nucsize to ya+nucsize do
  620.     If EraseMode then
  621.       NewGrayvalue(j,k,oldgrayvalue(j,k) and $FE)
  622.     else if oldgrayvalue(j,k) and 1 = 1 then
  623.       Redcount := redcount+1;
  624.     If RedCount <= 3 then
  625.        EraseMode := TRUE;
  626. end;
  627.  
  628. {........use a simple matrix filter again to remove small spots of RED......}
  629.  
  630. for k := ya-nucsize to ya+nucsize do
  631.   for j := xa-nucsize to xa+nucsize do
  632.     if (oldgrayvalue(j,k) and 1 = 1) then            {matrix 3x3}
  633.     begin
  634.       q := 0;
  635.       for f := k-1 to k+1 do
  636.     for r := j-1 to j+1 do
  637.       if (oldgrayvalue(r,f) and 1 = 1) then
  638.         q := q + 1;
  639.       if q < 5 then
  640.     newgrayvalue(j,k,oldgrayvalue(j,k) and $FE);
  641.     end;
  642.  
  643. end;{end if findarea}
  644. end;{end procedure FillIn}
  645.  
  646. {This is the first Unit called and is part of the Spot-Scanner.  If looks
  647.  at three pixels and determines which is the brightest.  If it is above
  648.  the lowest allowable value and below the highest the routine looks to
  649.  see if there is a contasting region nearby.  If so then the value of
  650.  TRUE is returned.  This is an ON/OFF unit.}
  651.  
  652. Function IntensityCheck(x,y:word;Size:byte):boolean;
  653. Const Ratio = 1.05;                  {contrast ratio}
  654. Var j,k:word;
  655.     gray1,gray2 : byte;
  656.     high,low : word;
  657.     Rfactor,rf2  : double;
  658. Begin
  659.  high := 0;
  660.  low := 255;
  661.  IntensityCheck := false;            {get maximum brightness}
  662.  gray1 := Max(oldgrayvalue(x,y),oldgrayvalue(x-1,y));
  663.  gray1 := Max(gray1,oldgrayvalue(x+1,y));
  664.                      {check if within bounds}
  665.  If (gray1 > graystrike) then
  666.  begin
  667.    for j := x-size to x+size do
  668.    begin
  669.      gray2 := oldgrayvalue(j,y);
  670.      if gray2 > high then
  671.        high := gray2
  672.      else if gray2 < low then
  673.        low := gray2;
  674.     end;
  675.     Rfactor := High/(low+1);
  676.                      {check if horizontal contrast}
  677.     if ((Rfactor > ratio) and (gray1 > 0.9*high)) then
  678.     begin
  679.       low := 255;
  680.       high := 0;
  681.       for k := y-size to y+size do
  682.       begin
  683.     gray2 := oldgrayvalue(x,k);
  684.     if gray2 > high then
  685.       high := gray2
  686.     else if gray2 < low then
  687.       low := gray2;
  688.       end;
  689.       Rf2 := High/(low+1);
  690.                      {check if vertical contrast}
  691.       if ((Rf2 > ratio) and (gray1 > 0.9*high)) then
  692.       IntensityCheck := TRUE
  693.       else
  694.       IntensityCheck := FALSE;
  695.       end
  696.     else
  697.       IntensityCheck := FALSE;
  698.   end
  699.   else
  700.       IntensityCheck := FALSE;
  701. end;{end function IntensityCheck}
  702.  
  703. {This unit returns the %foreground and standard deviation of a pixel
  704.  sampling of the nucleus.  The SIZE of the sample is related to the
  705.  value of NUCSIZE.  The values returned are mainly for use when
  706.  Learning is required.}
  707.  
  708. Procedure HowMuchFore(x,y:word;size:byte;Var AmtFore,_stdev:double);
  709. Var
  710.     j,k      : word;
  711.     count    : byte;
  712.     mark     : byte;
  713.     gray1    : byte;
  714.     imagdata : imagtype2;
  715.     _mean    : double;
  716. begin
  717.     count := 0;
  718.     mark := 0;
  719.     for k := y-size to y+size do                  {sample region}
  720.       for j := x-size to x+size do
  721.       begin
  722.     gray1 := oldgrayvalue(j,k);               {get pixel value}
  723.     if ((gray1 > criticalvalue) or (gray1 < lowdiv)) then
  724.       count := count+1;                       {store if good}
  725.     if (gray1 > criticalvalue) then           {do not include nucleolus}
  726.     begin                                     { when calculating st. dev.}
  727.       mark := mark + 1;
  728.       imagdata[mark] := gray1;
  729.     end;
  730.       end;
  731.     _mean  := mean(imagdata,1,mark);
  732.     _stdev := stdev(imagdata,1,mark,_mean);
  733.     AmtFore := Count/(((2*size) + 1)*((2*size) + 1));
  734. end;
  735.  
  736. {This routine returns the average gray level of a sample and the
  737.  amount of nucleolus found.}
  738.  
  739. Function Mscan(x,y:word;size:byte;Var bstuff:double):byte;
  740. var
  741.     j,k          : word;
  742.     count,count2 : byte;
  743.     imagdata     : imagtype2;
  744.     gray1        : byte;
  745. begin
  746.   count := 0;
  747.   count2 := 0;
  748.   for k := y-size to y+size do              {sample region}
  749.     for j := x-size to x+size do
  750.     begin
  751.       gray1 := oldgrayvalue(j,k);
  752.       If gray1 > lowdiv then                {store values above nucleolus}
  753.       begin
  754.     count := count+1;
  755.     imagdata[count] := gray1;
  756.       end
  757.       else
  758.     count2 := count2+1;
  759.     end;
  760.     Bstuff := count2/(count2+count);        {number of nucleolus pixels}
  761.     Mscan := round(mean(imagdata,1,count)); {all values above nucleolus}
  762. end; {end procedure Mscan}
  763.  
  764.  
  765. {This procedure will scan top,bottom,left,and right cytoplasm values vs.
  766.  the nuclear gray level.  A critical ratio must be met.  In order to account
  767.  for size variation the scan begins at a maximum nuclear-radius and moves
  768.  inward if acceptable values are not found.  Then the values must fall with
  769.  a certain range in order to assure uniformity.  Finally, "random data" is
  770.  generated and compared with the limits.}
  771.  
  772. Function SpotContrast(x,y:word;nucsize:byte;Var goodifsmall:boolean):boolean;
  773. Const
  774.      cratio  = 1.3;            {critical upper nuc/cyt segment ratio}
  775.      uratio  = 1.03;
  776.      Ul      = 0.96;            {lower limit}
  777.      UsumMin = 3;            {minimum sum}
  778.      UsumMax = 9;             {maximum sum}
  779.      pzhigh  = 1.6;                   {random data thresholds}
  780.      pzlow   = 1.4;
  781.      pzszlow = 0.31;
  782.      pzszhigh = 0.36;
  783.      upzhigh = 0.36;
  784.      upzlow  = 0.30;
  785.      diffx  = 0.004;
  786.      upzszlow  = 6;
  787.      upzszhigh = 7.1;
  788. var j,k,
  789.     s,t          : word;
  790.     debug        : boolean;
  791.     nold         : byte;                   {minimum distance from nucleus}
  792.     notdone,                               {flag to check if routine is done}
  793.     continue     : boolean;                {flag to check ratio limits}
  794.     mean1,mean2,                           {nuclear and cytoplasm averages}
  795.     ratio        : double;                 {nuclear/cytoplasm ratio}
  796.     displ15,                               {displacements}
  797.     displ14,
  798.     displ13      : byte;
  799.     r1,r2,r3,r4,                           {individual ratios}
  800.     a,b,                                   {used with Uniformity_ratio}
  801.     uniform_ratio,                         {uniformity of ratios}
  802.     sumz,prodz   : double;                 {sums and products of ratios}
  803. begin
  804.   j := x;
  805.   k := y;
  806.   debug := false;
  807.   SpotContrast := FALSE;
  808.  
  809.   notdone := TRUE;
  810.   Nold    := 1+(nucsize shr 2);           {set smallest distance}
  811.   Mean2 := 0;                             {get nuclear sample value}
  812.   for s := x-1 to x+1 do
  813.      for t := y-1 to y+1 do
  814.       mean2:=oldgrayvalue(s,t)+Mean2;
  815.   Mean2 := Round(Mean2/3);
  816.  
  817.   If Mean2/3 > 0.98*graystrike then
  818.     goodifsmall := FALSE
  819.   else
  820.     goodifsmall := TRUE;
  821.  
  822.   If Mean2/3 > 0.95*graystrike then
  823.                       {scan for cytoplasm values}
  824.   WHILE (NotDone) DO                      {Repeat until good energy is     }
  825.   BEGIN                                   {achieved or NUCSIZE becomes too }
  826.                                           {small.                          }
  827.     displ15 := nucsize+3;                 {Displacement values }
  828.     displ14 := nucsize+2;
  829.     displ13 := nucsize+1;
  830.                                            {Sample of cytoplasm consists of
  831.                                            three points}
  832.     mean1:=oldgrayvalue(x-displ15,k)+oldgrayvalue(x-displ14,k)+
  833.       oldgrayvalue(x-displ13,k);
  834.  
  835.     ratio := Mean2/(Mean1+1);               {compute nuclear/cytoplasm ratio}
  836.     if debug then writeln('RATIO1: ',ratio);
  837.     r1:=ratio;
  838.     If (ratio > cratio) then              {If ratio is above threshold  }
  839.      continue := TRUE                     {then continue                }
  840.     else
  841.      continue := FALSE;
  842.     If continue then                      {get next cytoplasm value}
  843.     begin
  844.      Mean1:=oldgrayvalue(x+displ15,k)+oldgrayvalue(x+displ14,k)+
  845.         oldgrayvalue(x+displ13,k);
  846.  
  847.      ratio := Mean2/(mean1+1);
  848.      if debug then writeln('ratio2: ',ratio);
  849.      r2 := ratio;
  850.      If (ratio>cratio) then
  851.        continue := TRUE
  852.      else
  853.        continue := FALSE;
  854.     end;
  855.     If continue then
  856.     begin
  857.       Mean1:=oldgrayvalue(j,y+displ15)+oldgrayvalue(j,y+displ14)+
  858.         oldgrayvalue(j,y+displ13);
  859.  
  860.       ratio := Mean2/(mean1+1);
  861.       if debug then writeln('ratio3: ',ratio);
  862.     r3 := ratio;
  863.       If (ratio>cratio) then
  864.        continue := TRUE
  865.       else
  866.        continue := FALSE;
  867.     end;
  868.     If continue then
  869.     begin
  870.       Mean1:=oldgrayvalue(j,y-displ15)+oldgrayvalue(j,y-displ14)+
  871.         oldgrayvalue(j,y-displ13);
  872.  
  873.       ratio := Mean2/(mean1+1);
  874.       if debug then writeln('ratio4: ',ratio);
  875.       r4:=ratio;
  876.       If (ratio>cratio) then
  877.        continue := TRUE
  878.       else
  879.        continue := FALSE;
  880.     end;
  881. {    if continue then
  882.     begin
  883.       spotcontrast := TRUE;
  884.       notdone := false;
  885.     end;
  886.  }
  887.  
  888.   {If this point is reached then the individual ratios are ok. Now
  889.    generate random data to check if the relationships amoung these ratios
  890.    is compatable with the desired pattern.}
  891.  
  892.     If continue then
  893.     begin
  894.      A := MaxMinRatio(r1,r2);          {Uniform_ratio checks that the    }
  895.      B := MaxMinRatio(r3,r4);          {difference between the cytoplasm }
  896.      Uniform_ratio := A/B;             {gray levels on opposite sides is }
  897.      sumz := r1+r2+r3+r4;
  898.      prodz := r1*r2*r3*r4;
  899.      writeln('UNIFORM: ',uniform_ratio:5:3,' USUM: ',sumz:5:3,
  900.          'U*: ',prodz/sumz:5:5,'up: ',uniform_ratio*prodz/sumz:5:3);
  901.      If (Uniform_ratio > 0.8) and (Uniform_ratio < 3) and
  902.         (sumz < 10) and (Prodz/sumz < 2) then
  903.      begin
  904.        notdone := FALSE;
  905.        spotcontrast := TRUE;
  906.      end;
  907.    end;
  908.  
  909.     If Nucsize-1 > Nold then                     {Decrease distance from    }
  910.      Nucsize := nucsize-1                        {nucleus. If too small then}
  911.     else                                         {then end routine and pass }
  912.      NotDone := FALSE;                           {back FALSE.               }
  913.  
  914.  
  915.    END;
  916.  
  917. end;{end function SpotContrast}
  918.  
  919. {This routine will track around the edge of an object, where the boundary
  920.  is delimited by RED.  A box sets the limits on where the object is.  The
  921.  routine scans for the first RED pixel and start from there.}
  922.  
  923. Function ScanEdge(x1,y1,x2,y2:word):word;
  924. {                      7              This is the chain code. The numbers
  925.             6     0           represent eight orientations about
  926.           5    x    1         the center point.
  927.             4     2
  928.                3
  929. }
  930. Const
  931.     OffsetDir = 6;                    {Starting direction}
  932. var j,k : word;
  933.     x_old,y_old : word;
  934.     j_old,k_old : word;
  935.     ChainCode,
  936.     ChainStart   : byte;
  937.     foundfirst,
  938.     finished,
  939.     done        : boolean;
  940.     Perimeter   : word;
  941.                          {This subroutine is given the current (x,y)
  942.                           coordinates and chaincode.  It then
  943.                           calculates the new (x,y) coordinates to
  944.                           look for an edge.}
  945. Procedure GetPoint(Var x,y:word;ChainCode : byte);
  946. begin
  947.   Case ChainCode of
  948.   1:    x := x+1;                                {y unchanged}
  949.   2:  begin
  950.     x := x+1;
  951.     y := y+1;
  952.       end;
  953.   3:    y := y+1;
  954.   4:  begin
  955.     x := x-1;
  956.     y := y+1;
  957.       end;
  958.   5:    x := x-1;
  959.   6:  begin
  960.     x := x-1;
  961.     y := y-1;
  962.       end;
  963.   7:  y := y-1;
  964.   0:  begin
  965.     x := x+1;
  966.     y := y-1;
  967.       end;
  968.   end;{end case}
  969. end; {end procedure GetPoint}
  970.                               {This function transforms the chain code where
  971.                                the edge was found and determines how many
  972.                                chain codes from chain code '1' it is
  973.                                located going clockwise.}
  974. Function TransChain(ChainCode:byte):byte;
  975. Var
  976.     temp : byte;
  977. begin
  978.     temp := (7+ChainCode) mod 8;
  979.     TransChain := temp;
  980. end;{end function Transchain}
  981.  
  982. begin
  983.   Perimeter := 0;                                {perimeter is zero}
  984.   foundfirst := false;                           {look for first red}
  985.   finished := FALSE;
  986.   k := y1;
  987.   Repeat                                         {vertical values}
  988.     j := x1;
  989.     Repeat                                       {scan horizontally}
  990.       If oldgrayvalue(j,k) and 1 = 1 then
  991.     foundfirst := TRUE
  992.       else
  993.       j := j+1;
  994.     Until (j > x2) or FoundFirst;
  995.     If Not(FoundFirst) then
  996.       k := k+1;
  997.   Until (k > y2) or FoundFirst;
  998.  
  999.   If foundfirst then                                {did we find a RED?}
  1000.   begin
  1001.     x_old := j;                                     {Set to coordinates of}
  1002.     y_old := k;                                     {first RED pixel      }
  1003.     Perimeter := 1;
  1004.     chainCode := OffsetDir;                         {this is first direction}
  1005.                   {Within this Repeat loop we scan around the
  1006.                    entire object till we come back to the staring point}
  1007.     REPEAT                                          {scan whole object}
  1008.       Done := False;
  1009.       ChainStart := ChainCode;
  1010.       j_old := j;                                   {Save our position      }
  1011.       k_old := k;                                   {so we can look around  }
  1012.                                                     {in all eight directions}
  1013.                  {Within this loop we scan around a red point in
  1014.                   search of the next red (edge) point. If none are found
  1015.                   then there must be only one RED point and the routine is
  1016.                   done.}
  1017.       Repeat
  1018.        GetPoint(j,k,chaincode);                     {get point to scan}
  1019.        If oldgrayvalue(j,k) and 1 = 1 then          {is it RED?}
  1020.      done := TRUE
  1021.        else                                         {If not then             }
  1022.             ChainCode := (ChainCode+1) mod 8;          {look in next direction  }
  1023.        If Not(done) then                            {If we didn't find an    }
  1024.        begin                                        {edge reset center point.}
  1025.      j := j_old;
  1026.      k := k_old;
  1027.        end;
  1028.       Until done or (chaincode = chainstart); {then perimeter = 1}
  1029.       If (j = x_old) and (k = y_old) then            {did we go around object?}
  1030.      Finished := TRUE                            {if so then we are done}
  1031.       else
  1032.       begin                          {otherwise we                         }
  1033.      Perimeter := Perimeter+1;   {increment the perimeter and ROTATE   }
  1034.                                      {the chain code matrix around the edge}
  1035.          ChainCode := (OffsetDir + TransChain(ChainCode)) mod 8;
  1036.                   (*The formula above says:  We always start scanning in the
  1037.                     6 (OffsetDir) direction.  We simply figure how many
  1038.                     chain codes from chain code 1 we moved and add this to
  1039.                     OffsetDir.  Modular division by eight simply insures
  1040.                     that we only have eight chain codes.
  1041.         SUPPOSE that we are at (0,0) and the next edge is at (1,-1).  The
  1042.         chain code direction is 2. TransChain says this is 1 chain code
  1043.         away from chain code 1.  We want to start scanning for the next
  1044.         edge at chain code 6 RELATIVE to current point which is why 6 is
  1045.         added making    7
  1046.         seven the     6   0
  1047.         next chain   5  x  1      <--Starting axis for first point horizont.
  1048.         code.            .              New axis is diagonal relative to
  1049.                           . 7           the first.
  1050.                           6.  0
  1051.                          5  x          Basically, an algorithm was needed
  1052.                                        that would give us the first point
  1053.                                        to scan that was immediately after the
  1054.                                        imaginary line between the two x's in
  1055.                                        the clockwise direction.  If we rotate
  1056.                                        the line between the two x's and make
  1057.                                        it horizontal (1 chain code
  1058.                                        counterclockwise) we see that
  1059.                                        relative to the second x-point
  1060.                                        we are scanning at chain code 6.*)
  1061.       end;
  1062.     UNTIL finished;
  1063.     ScanEdge := Perimeter;
  1064.   end;
  1065. end;{end function scanedge}
  1066.  
  1067. {After all nuclei are found they are converted from being shaded
  1068.  RED to having a gray value of 20 (dark).}
  1069.  
  1070. Procedure MakeDark(x1,y1,x2,y2:word);
  1071. var j,k : word;
  1072. begin
  1073.   for k := y1 to y2 do
  1074.     for j := x1 to x2 do
  1075.       if oldgrayvalue(j,k) and 1 = 1 then
  1076.     newgrayvalue(j,k,20);
  1077. end; {end procedure MakeDark}
  1078.  
  1079. {This procedure scans the region and determines the distances from the
  1080.  center-point (x,y).  These distances are returned as the length and
  1081.  width of the nucleus (DA,DB).}
  1082.  
  1083. Procedure EScan(x,y:word;Nucsize,cv:byte;var da,db:byte);
  1084. Var
  1085.     j,k   : word;
  1086.     done  : boolean;
  1087.     temp,
  1088.     gray1,
  1089.     valx  : byte;
  1090. begin
  1091.   valx := round(0.95*cv);             {get threshold}
  1092.   done := FALSE;
  1093.   j := x;                             {scan along horizontal}
  1094.   While (j <= x+nucsize) and Not(done) do
  1095.   begin
  1096.     gray1 := oldgrayvalue(j,y);
  1097.     newgrayvalue(j,y,gray1 or 1);
  1098.     If (isitbackground(j,y,valx,3)) and (gray1 > lowdiv) then
  1099.       done := TRUE
  1100.     else
  1101.       j := j+1;
  1102.   end;
  1103.   temp := j-x;
  1104.   done := FALSE;
  1105.   j := x;
  1106.   While (j > x-nucsize) and Not(done) do
  1107.   begin
  1108.     gray1 := oldgrayvalue(j,y);
  1109.     newgrayvalue(j,y,gray1 or 1);
  1110.     If (isitbackground(j,y,valx,-3)) and (gray1 > lowdiv) then
  1111.       done := TRUE
  1112.     else
  1113.       j := j-1;
  1114.   end;
  1115.   da := temp + (x-j);
  1116.   done := FALSE;
  1117.   k := y;                             {scan along vertical}
  1118.   While (k <= y+nucsize) and Not(done) do
  1119.   begin
  1120.     gray1 := oldgrayvalue(x,k);
  1121.     newgrayvalue(x,k,gray1 or 1);
  1122.     If (isitbackgroundv(x,k,valx,3)) and (gray1 > lowdiv) then
  1123.       done := TRUE
  1124.     else
  1125.       k := k+1;
  1126.   end;
  1127.   temp := k-y;
  1128.   done := FALSE;
  1129.   k := y;
  1130.   While (k > y-nucsize) and Not(done) do
  1131.   begin
  1132.     gray1 := oldgrayvalue(x,k);
  1133.     newgrayvalue(x,k,gray1 or 1);
  1134.     If (isitbackgroundv(x,k,valx,-3)) and (gray1 > lowdiv) then
  1135.       done := TRUE
  1136.     else
  1137.       k := k-1;
  1138.   end;
  1139.   db := temp + (y-k);
  1140. end;{end procedure EScan}
  1141.  
  1142. {This is an Energy routine.  It sets up several concentric square shells
  1143.  around the point (x,y) and samples the pixel values.  It then computes
  1144.  a ratio with the center sample.  Based on the relationship of these
  1145.  ratios and the state of SEENBEFORE the function returns TRUE or FALSE.
  1146.  However, certain values require information from other units.
  1147.  GOODIFNUCLEOUS, if TRUE, says that the ratios are
  1148.  good only if there is a nucleolus in this nucleolus.  The presence of
  1149.  one is determined by other units so the value is returned.}
  1150.  
  1151. Function ShellScan(x,y,nucsize:word;Uncertain:boolean;
  1152.                     var goodifnucleolus:boolean):boolean;
  1153. var
  1154.     j,k,
  1155.     r                : word;
  1156.    rq                : double;
  1157.    rx                : array[1..4] of double;
  1158.    gray1             : byte;
  1159.    q1,q2,q3          : double;
  1160.    count             : word;
  1161.    sum               : double;
  1162.    s,w,
  1163.    lowcount          : byte;
  1164. begin
  1165.   r := 0;
  1166.   for k := y-1 to y+1 do                    {sample center}
  1167.     for j := x-1 to x+1 do
  1168.     begin
  1169.       gray1 := oldgrayvalue(j,k);
  1170.       If gray1 > lowdiv then
  1171.     r := r+gray1
  1172.       else
  1173.     r := r+graystrike;
  1174.     end;
  1175.   rq := r/9;                                 {center average}
  1176.   s := 1;
  1177.   For w := 1 to 4 do                         {Get four other samples}
  1178.   begin                                      { around the center    }
  1179.     count := 0;
  1180.     s := s+2;
  1181.     rx[w] := 0;
  1182.     for j := x-s to x+s do                   {get top and bottom}
  1183.     begin
  1184.       gray1 := oldgrayvalue(j,y-s);
  1185.       If gray1 > lowdiv then
  1186.       begin
  1187.     count := count+1;
  1188.     rx[w] := rx[w]+gray1;
  1189.       end;
  1190.       gray1 := oldgrayvalue(j,y+s);
  1191.       If gray1 > lowdiv then
  1192.       begin
  1193.     count := count+1;
  1194.     rx[w] := rx[w]+gray1;
  1195.       end;
  1196.     end;
  1197.     for k := y-(s-1) to y+(s-1) do           {get right and left sides}
  1198.     begin
  1199.       gray1 := oldgrayvalue(x-s,k);
  1200.       If gray1 > lowdiv then
  1201.       begin
  1202.     count := count+1;
  1203.     rx[w]:= rx[w]+gray1;
  1204.       end;
  1205.       gray1 := oldgrayvalue(x+s,k);
  1206.       If gray1 > lowdiv then
  1207.       begin
  1208.     count := count+1;
  1209.     rx[w] := rx[w]+gray1;
  1210.       end;
  1211.     end;
  1212.     rx[w] := rx[w]/count;
  1213.     rx[w] := rq/rx[w];                       {store value}
  1214.   end;{end w}
  1215.   lowcount := 0;
  1216.  
  1217.   q1 := rx[2]/rx[1];
  1218.   q2 := rx[3]/rx[2];
  1219.   q3 := rx[4]/rx[3];
  1220.   count := round(int(q1) + int(q2) + int(q3));
  1221.   sum := abs(rx[1]-rx[2])+abs(rx[2]-rx[3])+abs(rx[3]-rx[4]);
  1222.   writeln('BEGIN SHELL ROUTINE');
  1223.   writeln('r1: ',rx[1]:5:5,' r2: ',rx[2]:5:5,' r3: ',
  1224.       rx[3]:5:5,' r4: ',rx[4]:5:5);
  1225.  
  1226.  
  1227. If   ( (Uncertain and (rx[1] > 1)) or (rx[1] > 1.02) ) and
  1228.      (rx[1] < rx[2]) and  (rx[2] < rx[3]) and (rx[2] < 1.5) and
  1229.      (rx[3] - rx[4] < 0.03) and (rx[4] < 1.65) then
  1230. begin
  1231.       if rx[3] < rx[4] then
  1232.         goodifnucleolus := FALSE
  1233.       else
  1234.         goodifnucleolus := TRUE;
  1235.       ShellScan := TRUE;
  1236. end
  1237. else
  1238. begin
  1239.       ShellScan := FALSE;
  1240. end;
  1241.  
  1242. end; {end procedure ShellScan}
  1243.  
  1244. {As the procedure name says, this computes the center of gravity of
  1245.  a RED shaded nucleus.}
  1246.  
  1247. Procedure FindCenterGravity(x,y,nucsize:word;var xc,yc,w:word);
  1248. Var
  1249.    j,k   : word;
  1250.    gray1 : byte;
  1251.    xc1,
  1252.    yc1   : double;
  1253. begin
  1254.    w := 0;
  1255.    xc1 := 0;
  1256.    yc1 := 0;
  1257.    for k := y-nucsize to y+nucsize do
  1258.      for j := x-nucsize to x+nucsize do
  1259.      begin
  1260.     gray1 := oldgrayvalue(j,k);
  1261.     If gray1 and 1 = 1 then
  1262.     begin
  1263.       xc1 := xc1+j;
  1264.       yc1 := yc1+k;
  1265.       w  := w+1;
  1266.     end;
  1267.      end;
  1268.    If w > 0 then
  1269.    begin
  1270.      xc := round(xc1/w);
  1271.      yc := round(yc1/w);
  1272.    end
  1273.    else
  1274.    begin
  1275.      xc := round(xc1);
  1276.      yc := round(yc1);
  1277.    end;
  1278. end;{end procedure findcentergravity}
  1279.  
  1280. {This routine performs several unit functions.  The outputs of these
  1281.  units are evaluated in the main routine with the outputs of other
  1282.  units.  Thus, this is part of an energy routine.  The routine determines
  1283.  the amount of nucleolus, the kurtosis and st. dev. of the nucleus, two
  1284.  nuclear/cytoplasm ratios, and the st. dev. of the immediately surrounding
  1285.  cytoplasm.}
  1286.  
  1287. Procedure HistoAnalysis(x,y:word;nucsize:byte;var below:byte;
  1288.       var ku,stout,rx,rx2:double;Var CytCond,AbsCyt : boolean);
  1289. Const
  1290.     histratio = 0.94;
  1291.     histcrit  = 0.85;
  1292.     histdiff  = 0.015;
  1293.     histdiff2 = 0.019;
  1294. Var
  1295.     j,k       : word;
  1296.     gray1     : byte;
  1297.     imagdata,
  1298.     im2       : imagtype2;
  1299.     _mean,
  1300.     outval    : double;
  1301.     w,w2,
  1302.     xc,yc,
  1303.     ns,
  1304.     meanx     : word;
  1305.     meanq,
  1306.     c1,c2,
  1307.     c3,c4,
  1308.     tot,
  1309.     mn,stmn  : double;
  1310. begin                         {Get centergrav and area (w)}
  1311.    findcentergravity(x,y,nucsize,xc,yc,w);
  1312.    c1 := 0;
  1313.    c2 := 0;
  1314.    c3 := 0;
  1315.    c4 := 0;
  1316.    below := 0;
  1317.    ns := round(1.2*sqrt(w/pi));  {approximate radius slightly larger}
  1318.    w := 0;                       {than true radius.  Thus, we draw a}
  1319.    w2 := 0;                      {square around the nucleus.        }
  1320.    below := 0;
  1321.    for k := yc-ns to yc+ns do
  1322.      for j := xc-ns to xc+ns do
  1323.      begin
  1324.        gray1 := oldgrayvalue(j,k);
  1325.        if (gray1 and 1 <> 1) then   {Get data on surrounding cytoplasm}
  1326.        begin
  1327.       w := w+1;
  1328.       imagdata[w] := gray1;
  1329.        end
  1330.        else if (gray1 > lowdiv) then  {get data on nucleus}
  1331.        begin
  1332.       w2 := w2+1;
  1333.       im2[w2] := gray1;
  1334.        end
  1335.        else
  1336.      below := below+1;            {store data on nucleolus}
  1337.      end;
  1338.  
  1339.     for j := xc-ns to xc+ns do
  1340.     begin
  1341.       c1 := c1 + oldgrayvalue(j,yc-ns);
  1342.       c2 := c2 + oldgrayvalue(j,yc+ns);
  1343.     end;
  1344.     for k := yc-ns to yc+ns do
  1345.     begin
  1346.       c3 := c3 + oldgrayvalue(xc-ns,k);
  1347.       c4 := c4 + oldgrayvalue(xc+ns,k);
  1348.     end;
  1349.     tot := 1+ (2*ns);
  1350.     c1 := c1/tot;
  1351.     c2 := c2/tot;
  1352.     c3 := c3/tot;
  1353.     c4 := c4/tot;
  1354.  
  1355.  
  1356.     outval := round(mean(imagdata,1,w));   {get avg gray-level of cytoplasm}
  1357.     stout  := stdev(imagdata,1,w,outval);  {get st. dev. of cytoplasm}
  1358.     _mean  := mean(im2,1,w2);              {get mean gray-level of nucleus}
  1359.     ku := kurtosis(im2,w2,round(_mean));            {get ku of nucleus}
  1360.  
  1361.     c1 := outval/(1+c1);
  1362.     c2 := outval/(1+c2);
  1363.     c3 := outval/(1+c3);
  1364.     c4 := outval/(1+c4);
  1365.     If (c1 < histcrit) or (c2 < histcrit) or (c3 < histcrit) or
  1366.        (c4 < histcrit) then
  1367.        Abscyt := FALSE
  1368.     else
  1369.        Abscyt := TRUE;
  1370.     mn := (c1+c2+c3+c4)/4;
  1371.     stmn := sqr(c1-mn) + sqr(c2-mn) + sqr(c3-mn) + sqr(c4-mn);
  1372.     stmn := sqrt(stmn)/3;
  1373.     If (c1 > histratio) and (c2 > histratio) and
  1374.        (c3 > histratio) and (c4 > histratio) and
  1375.        (mn < 1.1) and (mn > 0.95) and
  1376.        (stmn < 0.1) then
  1377.        CytCond := TRUE
  1378.     else
  1379.        CytCond := FALSE;
  1380.  
  1381.     writeln('ratios       : ',c1:5:5,' ',c2:5:5,' ',c3:5:5,' ',c4:5:5);
  1382.     writeln('diffs:         ',abs(c1-c3):5:5,' ',abs(c2-c4):5:5);
  1383.     writeln('sums :         ',c1+c2+c3+c4:5:5);
  1384.     writeln('mean,st :      ',mn:5:5,' ',stmn:5:5);
  1385.     writeln('d2   : ',abs(c1-mn):5:5,' ',abs(c2-mn):5:5,' ',abs(c3-mn):5:5,
  1386.             ' ',abs(c4-mn):5:5);
  1387.  
  1388.   meanx := 0;
  1389.   for k := y-1 to y+1 do             {sampe center 9 pixels}
  1390.     for j := x-1 to x+1 do
  1391.     begin
  1392.       gray1 := oldgrayvalue(j,k);
  1393.       If gray1 > lowdiv then
  1394.     meanx := meanx+gray1
  1395.       else
  1396.     meanx := meanx+graystrike;
  1397.      end;
  1398.   meanq := meanx/9;                  {get center ratio}
  1399.   If outval > 0 then
  1400.   begin
  1401.     rx := Meanq/outval;              {ratio of sample/cytoplasm}
  1402.     rx2 := _mean/outval;             {ratio of whole nucleus/cytoplasm}
  1403.   end
  1404.   else
  1405.   begin
  1406.     rx := 0;
  1407.     rx2 := 0;
  1408.   end;
  1409. end;{end procedure Histoanalysis}
  1410.  
  1411. {--------------------------LEARNING ALGORITHMS-------------------------}
  1412.  
  1413. Procedure LearnFromDeletion(Num:byte);
  1414. Var i,
  1415.     graylow,
  1416.     grayhigh : byte;
  1417.     AreaLow,
  1418.     AreaHigh,
  1419.     DaDbmin  : word;
  1420.     BlackMin,
  1421.     BlackMax,
  1422.     stqmax,
  1423.     stqxmax,
  1424.     forxmin,
  1425.     forgndmin,
  1426.     MvalMax,
  1427.     MvalMin,
  1428.     ShapeInd,
  1429.     ShapeMax,
  1430.     ShapeMin,
  1431.     Kumax,Kumin,
  1432.     cytomax,hypmin,
  1433.     rx2max,rx2min,
  1434.     rx1min,rx1max   : double;
  1435. begin
  1436.   graylow := 255;
  1437.   grayhigh := 0;
  1438.   AreaLow  := 9999;
  1439.   AreaHigh := 0;
  1440.   BlackMin := 250;
  1441.   blackMax := 0;
  1442.   StqMax   := 0;
  1443.   StqxMax  := 0;
  1444.   Forgndmin := 250;
  1445.   ForxMin  := 250;
  1446.   MvalMax  := 0;
  1447.   MvalMin  := 250;
  1448.   DaDbMin  := 250;
  1449.   ShapeMax := 0;
  1450.   ShapeMin := 250;
  1451.   kumax    := 0;
  1452.   kumin    := 255;
  1453.   rx2min     := 255;
  1454.   rx1min     := 255;
  1455.   rx2max     := 0;
  1456.   rx1max     := 0;
  1457.   cytomax    := 0;
  1458.   hypmin     := 255;
  1459.  
  1460.   For i := 1 to cellcount do               {set values to compare with}
  1461.   If AiCells[i].good then
  1462.   With Aicells[i] do
  1463.    begin
  1464.      If gray > grayhigh then
  1465.        grayhigh := gray
  1466.      else if gray < graylow then
  1467.        graylow := gray;
  1468.      If area > areahigh then
  1469.        areahigh := area
  1470.      else if area < arealow then
  1471.        arealow := area;
  1472.      If (black > blackmax) then
  1473.        blackmax := black
  1474.      else if (black < blackmin) and (black <> 0) then
  1475.        blackmin := black;
  1476.      If _stdev > stqmax then
  1477.        stqmax := _stdev;
  1478.      If stdx > stqxmax then
  1479.        stqxmax := stdx;
  1480.      If Mval > MvalMax then
  1481.        MvalMax := Mval
  1482.      else if Mval < MvalMin then
  1483.        MvalMin := Mval;
  1484.      If (black = 0) and (dadb < dadbmin) then
  1485.        DaDbmin := dadb;
  1486.      If area/dadb < hypmin then
  1487.        hypmin := area/dadb;
  1488.      If forx < ForxMin then
  1489.        ForxMin := Forx;
  1490.      If foregnd < Forgndmin then
  1491.        forgndmin := foregnd;
  1492.      ShapeInd := perimeter*perimeter/(Area*12.56);
  1493.      If shapeInd > ShapeMax then
  1494.        ShapeMax := shapeInd
  1495.      else if shapeInd < shapeMin then
  1496.        Shapemin := shapeInd;
  1497.      If kux > kumax then
  1498.        kumax := kux
  1499.      else if kux < kumin then
  1500.        kumin := kux;
  1501.      If cytost > cytomax then
  1502.        cytomax := cytost;
  1503.      If rx2 < rx2min then
  1504.        rx2min := rx2
  1505.      else if rx2 > rx2max then
  1506.        rx2max := rx2;
  1507.      If rx1 < rx1min then
  1508.        rx1min := rx1
  1509.      else if rx1 > rx1max then
  1510.        rx1max := rx1;
  1511.    end;
  1512.  
  1513.    If BlackMin > 1 then
  1514.      blackMin := MinBlack;
  1515.    If BlackMax = 0 then
  1516.      blackMax := MaxBlack;
  1517.  
  1518.   With Aicells[num] do
  1519.   begin
  1520.     If (gray = graylow) then
  1521.       graystrike := round(1.01*graylow)
  1522.     else if (gray = grayhigh) then
  1523.       CriticalHigh := round(0.99*grayhigh);
  1524.     If (area = arealow) then
  1525.       MinArea := round(1.01*AreaLow)
  1526.     else if (area = areaHigh) then
  1527.       MaxArea := round(0.99*AreaHigh);
  1528.     If (Black = BlackMin) and (Black >= MinBlack) then
  1529.       MinBlack := (1.01)*BlackMin
  1530.     else if (black = blackMax) and (blackmax > blackmin) then
  1531.       MaxBlack := (0.99)*BlackMax;
  1532.     If (_stdev = stqmax) then
  1533.       _stqset := (0.99)*Stqmax;
  1534.     If (stdx = stqxmax) then
  1535.       _stqxset := (0.99)*stqxmax;
  1536.     If (Mval = MvalMin) then
  1537.       Mvalx := (1.01)*Mvalmin;
  1538.  
  1539.    If (black = 0) and (dadb = dadbmin) then
  1540.      DaDbx := dadbmin+2;
  1541.    If area/dadb < hypmin then
  1542.      lowhyp := 1.01*hypmin;
  1543.    If forx = forxmin then
  1544.       forxset := 1.01*ForxMin;
  1545.    If foregnd = forgndmin then
  1546.       forset := 1.01*forgndmin;
  1547.  
  1548.    shapeInd := perimeter*perimeter/(12.56*Area);
  1549.    If ShapeInd > ShapeMax then
  1550.      ShapeHigh := 0.99*ShapeInd
  1551.    else if shapeInd < shapeMin then
  1552.      ShapeLow := 1.01*ShapeInd;
  1553.  
  1554.    If kux < kumin then
  1555.      kulow := 1.01*kumin
  1556.    else if kux > kumax then
  1557.      kuhigh := 0.99*kumax;
  1558.    If cytoset > cytomax then
  1559.      cytoset := 0.99*cytomax;
  1560.  
  1561.    If rx2 < rx2min then
  1562.      rx2low := 1.01*rx2min
  1563.    else if rx2 > rx2max then
  1564.      rx2high := 0.99*rx2max;
  1565.    If rx1 < rx1min then
  1566.      rx1low := 1.01*rx1min
  1567.    else if rx1 > rx1max then
  1568.      rx1high := 0.99*rx1max;
  1569.   end;
  1570. end;{end procedure learnfromdeletion}
  1571.  
  1572. Procedure LearnFromAddition(x,y:word;Nucsize,Width,Height:byte;Mval:double);
  1573. Var i : byte;
  1574.     Grayx : byte;
  1575.     aq,pq : word;
  1576.     mvq      : double;
  1577.     ShapeInd,
  1578.     _forex,_stdx : double;
  1579.     Cmval : double;
  1580.     blackcomp : double;
  1581.     da,db : byte;
  1582.     s,t : word;
  1583.     Mhigh : byte;
  1584.     xm,ym : word;
  1585.     forecomp,_stdev : double;
  1586.     below    : byte;
  1587.     ku,stout,rxa,rx2a : double;
  1588.     cytcond,abscyt : boolean;
  1589.  
  1590. begin
  1591.  
  1592.    histoanalysis(x,y,20,below,ku,stout,rxa,rx2a,cytcond,abscyt);
  1593.  
  1594.    If rxa < rx2a then
  1595.    begin
  1596.      mvq := rx2a;
  1597.      rx2a := rxa;
  1598.      rxa  := mvq;
  1599.    end;
  1600.  
  1601.    If ku < kulow then
  1602.      kulow := 0.99*ku
  1603.    else if ku > kuhigh then
  1604.      kuhigh := 1.01*ku;
  1605.    If stout > cytoset then
  1606.      cytoset := 1.01*stout;
  1607.    If rx2a < rx2low then
  1608.      rx2low := 0.99*rx2a
  1609.    else if rx2a > rx2high then
  1610.      rx2high := 1.01*rx2a;
  1611.  
  1612.    If rxa > rx1high then
  1613.      rx1high := 1.01*rxa
  1614.    else if rxa < rx1low then
  1615.      rx1low := 0.99*rxa;
  1616.  
  1617.  
  1618.    Grayx := oldgrayvalue(x,y);
  1619.    mvq   := getgray(x,y,5);
  1620.    If (0.96*mvq) < criticalvalue then
  1621.    begin
  1622.      If 0.95*mvq < mvalx then
  1623.        mvalx := 0.95*mvq;
  1624.      criticalvalue := (criticalvalue + round(0.96*mvq)) shr 1;
  1625.      graystrike    := (graystrike + round(1.03*criticalvalue)) shr 1;
  1626.    end;
  1627.  
  1628.    if 0.98*Grayx > criticalhigh then
  1629.      CriticalHigh := round(1.01*criticalHigh);
  1630.  
  1631.    Aq := 1+findarea(x-(width shr 1),y-(height shr 1),
  1632.              x+(width shr 1),y+(height shr 1),_forex,_stdx);
  1633.    Pq := scanedge(x-(width shr 1),y-(height shr 1),
  1634.              x+(width shr 1),y+(height shr 1));
  1635.  
  1636.    escan(x,y,nucsize,round(mval),da,db);
  1637.  
  1638.  
  1639.    shapeInd := pq*pq/(12.56*Aq);
  1640.    If (shapeind > shapehigh) then
  1641.      ShapeHigh := 1.01*ShapeInd
  1642.    else if (shapeind < shapelow) then
  1643.      ShapeLow := 0.99*ShapeInd;
  1644.  
  1645.    If (Aq < Minarea) then
  1646.      MinArea := round(0.99*Aq)
  1647.    else if (Aq > MaxArea) then
  1648.      MaxArea := round(1.01*Aq);
  1649.  
  1650.    If (_forex < forxset) then
  1651.      ForxSet := 0.99*_Forex;
  1652.    If (_stdx > _stqxset) then
  1653.      _StqxSet := 1.01*_Stdx;
  1654.  
  1655.    If (da*db <> 0) and (aq/(da*db) > dadbq) then
  1656.      dadbq := aq/(da*db);
  1657.    If (da*db <> 0) and (aq/(da*db) < lowhyp) then
  1658.      lowhyp := aq/(da*db);
  1659.    Cmval := Mscan(x,y,round(nucsize/3),blackcomp);
  1660.  
  1661.    If (blackcomp < Minblack) and (blackcomp <> 0) then
  1662.        MinBlack := 0.99*Blackcomp
  1663.    else if (blackcomp > Maxblack) and (blackcomp < 1) then
  1664.        MaxBlack := 1.01*Blackcomp;
  1665.  
  1666.    If Cmval < Mvalx then
  1667.      Mvalx := 0.99*cMval;
  1668.    If (blackcomp = 0) and ((da*db)-1 > 0) and (da*db < dadbx) then
  1669.      DaDbx := round(da*db);
  1670.  
  1671.    Mhigh := 0;
  1672.    for t := y-2 to y+2 do
  1673.      for s := x-2 to x+2 do
  1674.      begin
  1675.        grayx := oldgrayvalue(s,t);
  1676.        if grayx > Mhigh then
  1677.        begin
  1678.      Mhigh := grayx;
  1679.      xm := s;
  1680.      ym := t;
  1681.        end;
  1682.      end;
  1683.      Howmuchfore(xm,ym,(nucsize shr 2)+1,forecomp,_stdev);
  1684.      If (forecomp < forset) then
  1685.        forset := 0.99*forecomp;
  1686.      If (_stdev > _stqset) then
  1687.        _StqSet := 1.01*_stdev;
  1688.  
  1689.    {-------enter data for printer report---------------}
  1690.    cellcount := cellcount+1;
  1691.    with aicells[cellcount] do
  1692.    begin
  1693.      area := aq;
  1694.      perimeter := pq;
  1695.      _area := aq;
  1696.      _perim := pq;
  1697.      good := true;
  1698.               { ..........set values..........}
  1699.      gray := grayx;
  1700.      black := MinBlack;
  1701.      foregnd := forset;
  1702.      _stdev := _stqset;
  1703.      _forex := forxset;
  1704.      _stdx  := _stqxset;
  1705.      dadb   := da*db;
  1706.      mval   := mvalx;
  1707.      xcoord := x;
  1708.      ycoord := y;
  1709.      rx1    := rxa;
  1710.      rx2    := rx2a;
  1711.      cytost := stout;
  1712.      kux    := ku;
  1713.    end;
  1714. end;  {end procedure learnfromaddition}
  1715.  
  1716. END.
  1717.