unit EDM;

{Euclidian distance maps, ultimate eroded points, and watershed segmentation}

interface

	uses
		Types, Memory, QuickDraw, QuickDrawText, Packages, Menus, Events, Fonts, 
		Scrap, ToolUtils, Resources, Errors, Palettes, StandardFile, Windows,
		Controls, TextEdit, Files, Dialogs, TextUtils, Finder, MixedMode, Processes,
		globals, Utilities, Graphics, Analysis;


	procedure MakeEDM(item: integer);


implementation

	type
		Image16Data = packed array[0..maxImageSize] of integer;
		Image16Ptr = ^Image16Data;
		StartsArray = array[0..255] of LongInt;


function FindUtimatePoints(MaxEDM,item: LongInt; image2: ImageP; 
		LevelStart: StartsArray; xCoordinate, yCoordinate: Image16Ptr): boolean;
{
Finds peaks in the EDM that contain pixels equal to or greater than
all of their neighbors.
}
var
	x, y, level, rowsize, offset, i, count: LongInt;
	CoordOffset, xmax, ymax: LongInt;
	NextTick: LongInt;
	image: ImageP;
	SetPixel: boolean;
begin
	with info^ do begin
		image := ImageP(PicbaseAddr);
		BlockMove(image, image2, PixMapSize);
		rowsize := BytesPerRow;
		xmax := PixelsPerLine - 1;
		ymax := nLines - 1;
		NextTick := TickCount + 15;
		for level := maxEDM - 1 downto 1 do begin
			repeat
				count := 0;
				for i := 0 to Histogram[level] - 1 do begin
					CoordOffset := LevelStart[level] + i;
			      	x:= xCoordinate^[CoordOffset];
		  			y:= yCoordinate^[CoordOffset];
					offset := x + y * rowsize;
					if image^[offset] <> 255 then begin
						SetPixel := false;
						if (x > 0) and (y > 0) then
							if image^[offset - rowSize - 1] > level then
								SetPixel := true;
						if y > 0 then
							if image^[offset - rowSize] > level then
								SetPixel := true;
						if (x < xmax) and (y > 0) then
							if image^[offset - rowSize + 1] > level then
								SetPixel := true;
						if x < xmax then
							if image^[offset + 1] > level then
								SetPixel := true;
						if (x < xmax) and (y < ymax) then
							if image^[offset + rowSize + 1] > level then
								SetPixel := true;
						if y < ymax then
							if image^[offset + rowSize] > level then
								SetPixel := true;
						if (x > 0) and (y < ymax) then
							if image^[offset + rowSize - 1] > level then
								SetPixel := true;
						if x > 0 then
							if image^[offset - 1] > level then
								SetPixel := true;
						if SetPixel then begin
							image^[offset] := 255;
							count := count + 1;
						end;
					end; {if}
				end; {for i}
			until count = 0;
			if TickCount > NextTick then begin
				ShowAnimatedWatch;
				if CommandPeriod then begin
					beep;
					undo;
					WhatToUndo := NothingToUndo;
					FindUtimatePoints := false;
					exit(FindUtimatePoints);
				end;
				NextTick := TickCount + 15;
			end;
		end; {for level}
		if item = WatershedItem then begin
			for i := 0 to info^.PixMapSize - 1 do
				if (image^[i] > 0) and (image^[i] < 255) then
					image2^[i] := 255;
			BlockMove(image2, image, PixMapSize);
		end else begin
			for i := 0 to info^.PixMapSize - 1 do
				if image^[i] = 255 then
					image^[i] := 0;
		end;
		FindUtimatePoints := true;
	end; {with}
end;


procedure DoWatershedSegmentation(MaxEDM: LongInt; image2: ImageP; 
          LevelStart: StartsArray; xCoordinate, yCoordinate: Image16Ptr;
          var iterations: LongInt);
{
Assumes the current image contains an EDM and that the peaks in the EDM
(the ultimate eroded points) have been set to 255. The EDM is dilated from
these peaks, starting at the highest peak and working down. At each level,
the dilation is constrained to pixels whose values are at that level and
also constrained to prevent features from merging.
}
var
	rowSize, NextTicks: LongInt;
	level, x, y, offset, xmax, ymax: LongInt;
	count, FirstCount: LongInt;
	table: FateTable;
	image: ImageP;

	procedure CheckAbort;
	begin
		UpdatePicWindow;
		if CommandPeriod then begin
			beep;
			undo;
			WhatToUndo := NothingToUndo;
			exit(DoWatershedSegmentation);
		end;
		NextTicks := TickCount + 30;
		if OptionKeyDown then
			wait(30);
	end;
	
	procedure MakeFateTable;
	{Add pixel on 1st pass if bit 0 is set, on 2nd pass if bit 1 is}
	{set, on 3rd pass if bit 2 is set, and on 4th pass if bit 3 is set.}
	{E.g. '4' = add on 3rd pass, '3' = add on either 1st or 2nd pass,}
	{'f' = add on any pass. There is a routine in 'user.p' that draws}
	{all 256 neighborhoods.}
	const
		s999 = '01234567890123456789012345678901';
		
		s000 = '0044004480cc80cc0000000080cc80cc';
		s032 = '1000000080ff80ff1000000090ff90ff';
		s064 = '00000000000000000000000000000000';
		s096 = '1000000090ff90ff1000000090ff90ff';
		s128 = '2266006600ff00ff0000000000ff00ff';
		s160 = '13ff00ffffffffff33ff00ffffffffff';
		s192 = '2266006600ff00ff0000000000ff00ff';
		s224 = '33ff00ffffffffff33ff00fffffffff';
	var
		s: str255;
		c: char;
		i: LongInt;
	begin
		s := concat(s000, s032, s064, s096, s128, s160, s192, s224);
		for i := 0 to 254 do begin
			c := s[i + 1];
			if c <= '9' then
				table[i] := ord(c) - ord('0')
			else
				table[i] := 10 + ord(c) - ord('a')
		end;
		table[255] := 15; {f}
	end;

	procedure ProcessLevel(level, pass: LongInt);
	var
		index, x, y, i, CoordOffset, offset: LongInt;
	begin
		BlockMove(image, image2, info^.PixMapSize);
		for i := 0 to Histogram[level] - 1 do begin
			CoordOffset := LevelStart[level] + i;
	      	x:= xCoordinate^[CoordOffset];
  			y:= yCoordinate^[CoordOffset];
			offset := x + y * rowsize;
			if image2^[offset] <> 255 then begin
				index := 0;
				if (x > 0) and (y > 0) then
					if image2^[offset - rowSize - 1] = 255 then
						index := bor(index, 1);
				if y > 0 then
					if image2^[offset - rowSize] = 255 then
						index := bor(index, 2);
				if (x < xmax) and (y > 0) then
					if image2^[offset - rowSize + 1] = 255 then
						index := bor(index, 4);
				if x < xmax then
					if image2^[offset + 1] = 255 then
						index := bor(index, 8);
				if (x < xmax) and (y < ymax) then
					if image2^[offset + rowSize + 1] = 255 then
						index := bor(index, 16);
				if y < ymax then
					if image2^[offset + rowSize] = 255 then
						index := bor(index, 32);
				if (x > 0) and (y < ymax) then
					if image2^[offset + rowSize - 1] = 255 then
						index := bor(index, 64);
				if x > 0 then
					if image2^[offset - 1] = 255 then
						index := bor(index, 128);
				case pass of
					1: if band(table[index], 1) = 1 then begin
							image^[offset] := 255;
							count := count + 1;
						end;
					2:  if band(table[index], 2) = 2 then begin
							image^[offset] := 255;
							count := count + 1;
						end;
					3:  if band(table[index], 4) = 4 then begin
							image^[offset] := 255;
							count := count + 1;
						end;
					4: if band(table[index], 8) = 8 then begin
							image^[offset] := 255;
							count := count + 1;
						end;
				end; {case}
			end; {if}
		end; {for i}
	end; {ProcessLevel}

	procedure PostProcess;
	var
		i: LongInt;
	begin
		for i := 0 to info^.PixMapSize - 1 do
			if image^[i] < 255 then
				image^[i] := 0
	end;
	
begin
	with info^ do begin
		rowSize := BytesPerRow;
		MakeFateTable;
		image := ImageP(PicbaseAddr);
		xmax := PixelsPerLine - 1;
		ymax := nLines - 1;
		NextTicks := TickCount;
		iterations := 0;
		for level := maxEDM - 1 downto 1 do begin
			repeat
				count := 0;
				ProcessLevel(level, 1);
				ProcessLevel(level, 3);
				ProcessLevel(level, 2);
				ProcessLevel(level, 4);
				iterations := iterations + 1;
			until count = 0;
		if TickCount > NextTicks then
			CheckAbort;
  	end; {for level}
	PostProcess;
	UpdatePicWindow;
  end;
end;


procedure SmoothEDM(image2: ImageP);
var
	x, y, rowsize, offset, sum: LongInt;
	xmax, ymax: LongInt;
	image: ImageP;
begin
	with info^ do begin
		image := ImageP(PicbaseAddr);
		BlockMove(image, image2, PixMapSize);
		rowsize := BytesPerRow;
		xmax := PixelsPerLine - 1;
		ymax := nLines - 1;
		for y := 0 to nLines - 1 do begin
			for x := 0 to PixelsPerLine - 1 do begin
				offset := x + y * rowsize;
				if image2^[offset] <> 1 then begin
					sum := image2^[offset] * 2;
					if (x > 0) and (y > 0) then
						sum := sum + image2^[offset - rowSize - 1];
					if y > 0 then
						sum := sum + image2^[offset - rowSize];
					if (x < xmax) and (y > 0) then
						sum := sum + image2^[offset - rowSize + 1];
					if x < xmax then
						sum := sum + image2^[offset + 1];
					if (x < xmax) and (y < ymax) then
						sum := sum + image2^[offset + rowSize + 1];
					if y < ymax then
						sum := sum + image2^[offset + rowSize];
					if (x > 0) and (y < ymax) then
						sum := sum + image2^[offset + rowSize - 1];
					if x > 0 then
						sum := sum + image2^[offset - 1];
					image^[offset] := sum div 10;
					end;
			end; {for x}
		end; {for y}
	end; {with}
end;


procedure MakeEDM(item: integer);
{
Converts a binary image into a grayscale Euclidian distance
map (EDM). Each foreground (black) pixel in the binary image
is assigned a value equal to its distance from the nearest 
background (white) pixel. Uses the two pass EDM algorithm
from the "Image Pricessing Handbook" by John Russ.  
 }
const
	one = 41;
	sqrt2 = 58; {~41*sqrt(2)}
	sqrt5 = 92; {~41*sqrt(5)}

var
	x, y, xmax, ymax, i, MaxEDM: LongInt;
	StartTicks, NextTicks, offset, rowsize: LongInt;
	iterations: LongInt;
	PointsOK: boolean;
	image16: Image16Ptr;
	image2: ImageP;
	str: str255;
	LevelStart: StartsArray;
	xCoordinate, yCoordinate: Image16Ptr;
	
	function ConvertToIntegers:boolean;
	var
		x, y, offset: LongInt;
	begin
		with info^ do begin
			image16 := Image16Ptr(NewPtr(rowsize * nLines * SizeOf(integer)));
			if image16 = nil then begin
				ConvertToIntegers := false;
				exit(ConvertToIntegers);
			end;
			for y := 0 to nLines - 1 do
			  for x := 0 to PixelsPerLine - 1 do begin
			    offset := x + y * rowsize;
					image16^[offset] := ImageP(PicBaseAddr)^[offset] * one;
			  end;
			ConvertToIntegers := true;
		end;
	end;

	procedure ConvertToBytes;
	var
		x, y, v, round, offset: LongInt;
	begin
		round := one div 2;
		MaxEDM := 0;
		with info^ do begin
			for y := 0 to nLines - 1 do
			  for x := 0 to PixelsPerLine - 1 do begin
			    offset := x + y * rowsize;
			    v := (image16^[offset] + round) div one;
			    if v > 255 then
						v := 255;
					if v > maxEDM then
						maxEDM := v;
					ImageP(PicBaseAddr)^[offset] := v;
			  end;
		end;
	end;

	procedure SetEdgeValue;
	var
		min, inc, v: LongInt;
		r1, r2, r3, r4, r5, offimage: LongInt;
	begin
		r1 := offset - rowsize - rowsize - 2;
		r2 := r1 + rowsize;
		r3 := r2 + rowsize;
		r4 := r3 + rowsize;
		r5 := r4 + rowsize;
		min := 32767;
	  	
		offimage := image16^[r3 + 2];
		
		if y < 2 then
			v := offImage + one
		else
			v := image16^[r2 + 2] + one;
		if v < min then
			min := v;
			
		if x < 2 then
			v := offImage + one
		else
			v := image16^[r3 + 1] + one;
		if v < min then
			min := v;
			
		if x > xmax then
			v := offImage + one
		else
			v := image16^[r3 + 3] + one;
		if v < min then
			min := v;
			
		if y > ymax then
			v := offImage + one
		else
			v := image16^[r4 + 2] + one;
		if v < min then
			min := v;
		
		if (x < 2) or (y < 2) then
			v := offImage + sqrt2
		else
			v := image16^[r2 + 1] + sqrt2;
		if v < min then
			min := v;
			
		if (x > xmax) or (y < 2) then
			v := offImage + sqrt2
		else
			v := image16^[r2 + 3] + sqrt2;
		if v < min then
			min := v;
			
		if (x < 2) or (y > ymax) then
			v := offImage + sqrt2
		else
			v := image16^[r4 + 1] + sqrt2;
		if v < min then
			min := v;
			
		if (x > xmax) or (y > ymax) then
			v := offImage + sqrt2
		else
			v := image16^[r4 + 3] + sqrt2;
		if v < min then
			min := v;
		
		if (x < 2) or (y < 2) then
			v := offImage + sqrt5
		else
			v := image16^[r1 + 1] + sqrt5;
		if v < min then
			min := v;
			
		if (x > xmax) or (y < 2) then
			v := offImage + sqrt5
		else
			v := image16^[r1 + 3] + sqrt5;
		if v < min then
			min := v;
			
		if (x > xmax) or (y < 2) then
			v := offImage + sqrt5
		else
			v := image16^[r2 + 4] + sqrt5;
		if v < min then
			min := v;
			
		if (x > xmax) or (y > ymax) then
			v := offImage + sqrt5
		else
			v := image16^[r4 + 4] + sqrt5;
		if v < min then
			min := v;
			
		if (x > xmax) or (y > ymax) then
			v := offImage + sqrt5
		else
			v := image16^[r5 + 3] + sqrt5;
		if v < min then
			min := v;
			
		if (x < 2) or (y > ymax) then
			v := offImage + sqrt5
		else
		v := image16^[r5 + 1] + sqrt5;
		if v < min then
			min := v;
			
		if (x < 2) or (y > ymax) then
			v := offImage + sqrt5
		else
			v := image16^[r4] + sqrt5;
		if v < min then
			min := v;
			
		if (x < 2) or (y < 2) then
			v := offImage + sqrt5
		else
			v := image16^[r2] + sqrt5;
		if v < min then
			min := v;
		
		image16^[offset] := min;
	end; {SetEdgeValue}

	procedure SetValue;
	var
		min, inc, v: LongInt;
		r1, r2, r3, r4, r5: LongInt;
	begin
		r1 := offset - rowsize - rowsize - 2;
		r2 := r1 + rowsize;
		r3 := r2 + rowsize;
		r4 := r3 + rowsize;
		r5 := r4 + rowsize;
		min := 32767;
	  	
		v := image16^[r2 + 2] + one;
		if v < min then
			min := v;
		v := image16^[r3 + 1] + one;
		if v < min then
			min := v;
		v := image16^[r3 + 3] + one;
		if v < min then
			min := v;
		v := image16^[r4 + 2] + one;
		if v < min then
			min := v;
		
		v := image16^[r2 + 1] + sqrt2;
		if v < min then
			min := v;
		v := image16^[r2 + 3] + sqrt2;
		if v < min then
			min := v;
		v := image16^[r4 + 1] + sqrt2;
		if v < min then
			min := v;
		v := image16^[r4 + 3] + sqrt2;
		if v < min then
			min := v;
		
		v := image16^[r1 + 1] + sqrt5;
		if v < min then
			min := v;
		v := image16^[r1 + 3] + sqrt5;
		if v < min then
			min := v;
		v := image16^[r2 + 4] + sqrt5;
		if v < min then
			min := v;
		v := image16^[r4 + 4] + sqrt5;
		if v < min then
			min := v;
		v := image16^[r5 + 3] + sqrt5;
		if v < min then
			min := v;
		v := image16^[r5 + 1] + sqrt5;
		if v < min then
			min := v;
		v := image16^[r4] + sqrt5;
		if v < min then
			min := v;
		v := image16^[r2] + sqrt5;
		if v < min then
			min := v;
		
		image16^[offset] := min;
	end; {SetValue}
	

	procedure CheckAbort;
	begin
		ShowAnimatedWatch;
		if CommandPeriod then begin
			beep;
			DisposePtr(ptr(image16));
			exit(MakeEDM);
		end;
		NextTicks := TickCount + 15;
	end;


	procedure MakeCoordinateArrays;
	{Generates the XY coordinate arrays that allow pixels at each level to
	be accesed directly without searching through the entire image. This
	method, suggested by Stein Roervik (steinr@kjemi.unit.no), greatly
	speeds up the watershed segmentation routine.}
	var
		i, x, y, rowsize, offset, v, ArraySize: LongInt;
		LevelOffset: array[0..255] of LongInt;
		image: ImageP;
	begin
		with info^ do begin
			RoiRect := PicRect;
			GetRectHistogram;
			ArraySize := 0;
			for i := 1 to maxEDM - 1 do
				ArraySize := ArraySize + histogram[i];
			xCoordinate := Image16Ptr(NewPtr(ArraySize * SizeOf(integer)));
			yCoordinate := Image16Ptr(NewPtr(ArraySize * SizeOf(integer)));
			if (xCoordinate = nil) or (yCoordinate = nil) then begin
				if xCoordinate <> nil then
					DisposePtr(ptr(xCoordinate));
				DisposePtr(ptr(image2));
				PutError('Out of memory');
				undo;
				WhatToUndo := NothingToUndo;
				exit(MakeEDM);
			end;
			image := ImageP(PicbaseAddr);
			offset := 0;
			for i := 0 to 255 do begin
				LevelStart[i] := offset;
				if (i >= 1) and (i <= (MaxEDM - 1)) then
					offset := offset + histogram[i];
			end;
			for I := 0 to 255 do
				LevelOffset[i] := 0;
			rowsize := BytesPerRow;
			for y := 0 to nLines - 1 do begin
				for x := 0 to PixelsPerLine - 1 do begin
					v := image^[x + y * rowsize];
					if (v >= 1) and (v <= (MaxEDM - 1)) then begin
						offset := LevelStart[v] + LevelOffset[v];
						xCoordinate^[offset] := x;
						yCoordinate^[offset] := y;
						LevelOffset[v] := LevelOffset[v] + 1;
					end;
				end; {for x}
			end; {for y}
		end; {with}
	end;

begin
	with info^ do begin
		ShowWatch;
		KillRoi;
		rowsize := BytesPerRow;
		xmax := PixelsPerLine - 3;
		ymax := nLines - 3;
		StartTicks := TickCount;
		NextTicks := TickCount;
		SetupUndo;
		WhatToUndo := UndoFilter;
		if not ConvertToIntegers then begin
			AbortMacro;
			PutError('Out of Memory');
			exit(MakeEDM);
		end;
		for y := 0 to nLines - 1 do begin
			for x := 0 to PixelsPerLine - 1 do begin
				offset := x + y * rowsize;
				if image16^[offset] > 0.0 then begin
					if (x < 2) or (x > xmax) or (y < 2) or (y > ymax) then
						SetEdgeValue
					else
						SetValue;
				end;
			end;
			if TickCount > NextTicks then
				CheckAbort;
		end;
		for y := nLines - 1 downto 0 do begin
			for x := PixelsPerLine - 1 downto 0 do begin
		    offset := x + y * rowsize;
				if image16^[offset] > 0.0 then begin
					if (x < 2) or (x > xmax) or (y < 2) or (y > ymax) then
						SetEdgeValue
					else
						SetValue;
				end;
			end;
			if TickCount > NextTicks then
				CheckAbort;
		end;
		ConvertToBytes;
		DisposePtr(ptr(image16));
		if (item = UltimateItem) or (item = WatershedItem) then begin
			image2 := ImageP(NewPtr(PixMapSize));
			if image2 = nil then
				exit(MakeEDM);
			if ((item = WatershedItem) and not OptionKeyWasDown) or ((item = UltimateItem) and OptionKeyWasDown) then
				SmoothEDM(image2);
			MakeCoordinateArrays;
			PointsOK := FindUtimatePoints(MaxEDM, item, image2, LevelStart, xCoordinate, yCoordinate);
		end;
		if (item = WatershedItem) and PointsOK then begin
			DoWatershedSegmentation(MaxEDM, image2, LevelStart, xCoordinate, yCoordinate, iterations);
			str := StringOf(MaxEDM - 1:1, ' levels', crStr, iterations/(MaxEDM - 1):1:2, ' iterations/level');
		end else
			str := '';
		if (item = UltimateItem) or (item = WatershedItem) then begin
			DisposePtr(ptr(xCoordinate));
			DisposePtr(ptr(yCoordinate));
			DisposePtr(ptr(image2));
		end;
		ShowTime(StartTicks, PicRect, str);
		changes := true;
		UpdatePicWindow;
	end; {with}
end;

end.