unit Filters;

{}

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, File1, File2, Analysis, Camera, Lut, realUtils, Text;


	procedure ApplyLookupTable;
	procedure MakeBinary;
	procedure Filter (ftype: FilterType; pass: integer; var table: FateTable);
	procedure PhotoMode;
	function AllSameSize: boolean;
	procedure EnhanceContrast;
	procedure EqualizeHistogram;
	procedure Convolve (name: str255; RefNum: integer);
	procedure ConvolveUsingText;
	procedure PlotSurface;
	procedure MakeSkeleton;
	procedure DoErosion;
	procedure DoDilation;
	procedure DoOpening;
	procedure DoClosing;
	procedure SetBinaryCount;
	procedure SetIterations;
	procedure ChangeValues (v1, v2, v3: integer);
	procedure DoPropagate (MenuItem: integer);
	procedure NewPlotSurface;
	procedure AutoThreshold;
	procedure AutoDensitySlice;
	procedure ScaleAndRotate;
	procedure DoRankFilter;
	procedure DoShadowFilter;
	procedure MakeEDM(item: integer);


implementation

	const
		MaxW = 4000;

	type
		ktype = array[0..MaxW] of integer;
		SortArray = array[1..9] of integer;
		StartsArray = array[0..255] of LongInt;

	var
		PixelsRemoved: LongInt;
		DebugCount: integer;


{$PUSH}
{$D-}

	function DoApplyTableDialogBox: boolean;
		const
			Button1 = 3;
			Button2 = 4;
			Button3 = 5;
			Button4 = 6;
		var
			mylog: DialogPtr;
			item: integer;
			SaveA, SaveB: boolean;

		procedure SetButtons;
		begin
			SetDlogItem(mylog, Button1, ord(ThresholdToForeground));
			SetDlogItem(mylog, Button2, ord(not ThresholdToForeground));
			SetDlogItem(mylog, Button3, ord(NonThresholdToBackground));
			SetDlogItem(mylog, Button4, ord(not NonThresholdToBackground));
		end;

	begin
		InitCursor;
		SaveA := ThresholdToForeground;
		SaveB := NonThresholdToBackground;
		mylog := GetNewDialog(40, nil, pointer(-1));
		SetButtons;
		OutlineButton(MyLog, ok, 16);
		repeat
			ModalDialog(nil, item);
			if (item = Button1) or (item = button2) then begin
					ThresholdToForeground := not ThresholdToForeground;
					SetButtons;
				end;
			if (item = Button3) or (item = button4) then begin
					NonThresholdToBackground := not NonThresholdToBackground;
					SetButtons;
				end;
		until (item = ok) or (item = cancel);
		DisposeDialog(mylog);
		if item = cancel then begin
				ThresholdToForeground := SaveA;
				NonThresholdToBackground := SaveB;
				DoApplyTableDialogBox := false
			end
		else
			DoApplyTableDialogBox := true;
	end;


	procedure ApplyLookupTable;
		var
			table: LookupTable;
			ConvertingColorPic, GrayScaleImage: boolean;
			i: integer;
	begin
		with info^ do begin
				GrayScaleImage := (LUTMode = Grayscale) or (LUTMode = CustomGrayscale);
				ConvertingColorPic := not GrayScaleImage and not DensitySlicing;
				if ConvertingColorPic then
					KillRoi;
				if DensitySlicing and (not macro) then begin
						if not DoApplyTableDialogBox then
							exit(ApplyLookupTable);
					end;
				if thresholding then
					BinaryPic := true;
				GetLookupTable(table);
				if GrayscaleImage or ConvertingColorPic then
					ResetGrayMap;
				ApplyTable(table);
				if ConvertingColorPic then
					WhatToUndo := NothingToUndo;
				if fit <> uncalibrated then begin
						fit := uncalibrated;
						for i := 0 to 255 do
							cvalue[i] := i;
					end;
			end; {with}
	end;


	procedure MakeBinary;
		var
			table: LookupTable;
			SaveBackground, SaveForeground, i: integer;
	begin
		with info^ do begin
				if DensitySlicing then begin
						ThresholdToForeground := true;
						NonThresholdToBackground := true;
						SaveBackground := BackgroundIndex;
						SaveForeground := ForegroundIndex;
						BackgroundIndex := WhiteIndex;
						ForegroundIndex := BlackIndex;
						GetLookupTable(table);
						ResetGrayMap;
						ApplyTable(table);
						BackgroundIndex := SaveBackground;
						ForegroundIndex := SaveForeground;
						BinaryPic := true;
					end
				else if Thresholding then begin
						for i := 0 to 255 do
							if i < ColorStart then
								table[i] := WhiteIndex
							else
								table[i] := BlackIndex;
						ResetGrayMap;
						ApplyTable(table);
						BinaryPic := true;
					end
				else
					PutError('Sorry, but you must be thresholding or density slicing to use Make Binary.');
			end;
	end;


	function FindMedian (var a: SortArray): integer;
  {Finds the 5th largest of 9 values}
{$IFC PowerPC}
		var
			i, j, mj, max: integer;
	begin
		for i := 1 to 4 do begin
				max := 0;
				mj := 1;
				for j := 1 to 9 do
					if a[j] > max then begin
							max := a[j];
							mj := j;
						end;
				a[mj] := 0;
			end;
		max := 0;
		for j := 1 to 9 do
			if a[j] > max then
				max := a[j];
		FindMedian := max;
	end;
	
{$ELSEC}
   {In-line code contributed by Edward J. Huff(huff@mcclbo.med.nyu.edu).}
   {Assember source with comments and a test program are available by anonymous}
   {ftp from zippy.nimh.nih.gov, in the /pub/nih-image/documents directory.}
	inline
		$205F, $48E7, $1F00, $4C98, $00FF, $B041, $6502, $C340,{}
		$B443, $6502, $C742, $B243, $6504, $C540, $C741, $B845,{}
		$6502, $CB44, $BC47, $6502, $CF46, $BA47, $6504, $CD44,{}
		$CF45, $B245, $6508, $CF43, $CD42, $CB41, $C940, $3E10,{}
		$BC47, $6502, $CF46, $BA47, $6504, $CD44, $CF45, $B245,{}
		$6508, $CF43, $CD42, $CB41, $C940, $B246, $6534, $B242,{}
		$6514, $B244, $6504, $3001, $6062, $B644, $6504, $3004,{}
		$605A, $3003, $6056, $B444, $650C, $B445, $6504, $3005,{}
		$604A, $3002, $6046, $B644, $6504, $3004, $603E, $3003,{}
		$603A, $B645, $6504, $C942, $CB43, $B846, $651C, $B644,{}
		$650C, $B444, $6504, $3002, $6022, $3004, $601E, $B646,{}
		$6504, $3003, $6016, $3006, $6012, $B646, $6508, $B446,{}
		$65F4, $3002, $6006, $B644, $65E0, $3003, $4CDF, $00F8,{}
		$3E80;
{$ENDC}



	procedure PhotoMode;
{Erases the screen to the background color and then redraws}
{the contents of the active image window . }
		var
			tPort: GrafPtr;
			event: EventRecord;
			WinRect: rect;
			SaveVisRgn: rgnHandle;
	begin
		with info^ do begin
				KillRoi;
				if OptionKeyWasDown then begin {Move window up to top of screen.}
						GetWindowRect(wptr, WinRect);
						MoveWindow(wptr, WinRect.left, 0, false);
					end;
				with wptr^ do begin
						SaveVisRgn := visRgn;
						visRgn := NewRgn;
						RectRgn(visRgn, qd.ScreenBits.Bounds);
					end;
				FlushEvents(EveryEvent, 0);
				GetPort(tPort);
				EraseScreen;
				UpdatePicWindow;
				repeat
				until button;
				with wptr^ do begin
						DisposeRgn(visRgn);
						visRgn := SaveVisRgn;
					end;
				RestoreScreen;
				SetPort(tPort);
				FlushEvents(EveryEvent, 0);
				if OptionKeyWasDown then begin
						MoveWindow(wptr, WinRect.left, WinRect.top, false);
					end;
			end;
	end;


	function AllSameSize: boolean;
{Returns true if all currently open Images have the same dimensions.}
		var
			i: integer;
			SameSize: Boolean;
			TempInfo: InfoPtr;
	begin
		if nPics = 0 then begin
				AllSameSize := false;
				exit(AllSameSize);
			end;
		SameSize := true;
		for i := 1 to nPics do begin
				TempInfo := pointer(WindowPeek(PicWindow[i])^.RefCon);
				SameSize := SameSize and EqualRect(Info^.PicRect, TempInfo^.PicRect);
			end;
		AllSameSize := SameSize;
	end;


	procedure EnhanceContrast;
		var
			AutoSelectAll: boolean;
			min, max, i, threshold: integer;
			found, SaveRedirectFlag: boolean;
			sum: LongInt;
	begin
		with info^ do
			if LUTMode = ColorLUT then begin
					PutError('Sorry, but you can not contrast enhance true color images.');
					exit(EnhanceContrast)
				end;
		if NotInBounds or (ClipBuf = nil) then
			exit(EnhanceContrast);
		StopDigitizing;
		AutoSelectAll := not Info^.RoiShowing;
		if AutoSelectAll then
			SelectAll(false);
		SaveRedirectFlag := RedirectSampling;
		RedirectSampling := false;
		if info^.RoiType = RectRoi then
			GetRectHistogram
		else
			GetHistogram;
		RedirectSampling := SaveRedirectFlag;
		sum := 0;
		for i := 0 to 255 do
			sum := sum + histogram[i];
		threshold := sum div 5000;
		i := -1;
		repeat
			i := i + 1;
			found := histogram[i] > threshold;
		until found or (i = 255);
		min := i;
		i := 256;
		repeat
			i := i - 1;
			found := histogram[i] > threshold;
		until found or (i = 0);
		max := i;
		if max > min then
			with info^ do begin
					SetupLutUndo;
					if isGrayScaleLUT then
						LUTMode := grayscale;
					ColorStart := min;
					ColorEnd := max;
					DrawMap;
					UpdateLUT;
					changes := true;
					IdentityFunction := false;
				end;
		if AutoSelectAll then
			KillRoi;
		if ScreenDepth<>8 then UpdatePicWindow;
	end;


	procedure EqualizeHistogram;
		var
			AutoSelectAll, SaveRedirectFlag: boolean;
			i, sum, v: integer;
			isum: LongInt;
			ScaleFactor: extended;
	begin
		with info^ do
			if (LUTMode <> GrayScale) and (LutMode <> CustomGrayscale) then begin
					PutError('Sorry, but you can only do histogram equalization on grayscale images.');
					exit(EqualizeHistogram)
				end;
		if NotInBounds or (ClipBuf = nil) then
			exit(EqualizeHistogram);
		StopDigitizing;
		AutoSelectAll := not Info^.RoiShowing;
		if AutoSelectAll then
			SelectAll(false);
		SaveRedirectFlag := RedirectSampling;
		RedirectSampling := false;
		if info^.RoiType = RectRoi then
			GetRectHistogram
		else
			GetHistogram;
		RedirectSampling := SaveRedirectFlag;
		FindThresholdingMode;
		ComputeResults;
		isum := 0;
		for i := 0 to 255 do
			isum := isum + histogram[i];
		ScaleFactor := 255.0 / isum;
		sum := 0;
		with info^ do begin
				SetupLutUndo;
				for i := 255 downto 0 do
					with cTable[i].rgb do begin
							sum := round(sum + histogram[i] * ScaleFactor);
							if sum > 255 then
								sum := 255;
							v := sum * 256;
							red := v;
							green := v;
							blue := v;
						end;
				LoadLUT(cTable);
				LUTMode := CustomGrayscale;
				SetupPseudocolor;
				changes := true;
				DrawMap;
				IdentityFunction := false;
			end; {with info}
		if AutoSelectAll then
			KillRoi;
		if ScreenDepth<>8 then UpdatePicWindow;
	end;


	procedure GetKernel (var kernel: ktype; var n: integer; var name: str255; RefNum: integer);
		var
			rLine: RealLine;
			i, count, nValues, nRows: integer;
	begin
		count := 0;
		nRows := 0;
		InitTextInput(name, RefNum);
		while not TextEof and (nRows <= 63) do begin
				GetLineFromText(rLine, nValues);
				if count <> 0 then
					nRows := nRows + 1;
				if nRows = 1 then
					n := nValues;
				for i := 1 to nValues do begin
						count := count + 1;
						kernel[count - 1] := round(rLine[i]);
					end;
			end;
		if count <> (n * n) then
			n := 0;
	end;


	procedure DoOnePixel (nLess1, BytesPerLine: integer; corner: LongInt; var sum: LongInt; var kernel: ktype);
{$IFC PowerPC}
		 var
			row, column, k: integer;
			pp: ptr;
	begin
		k := 0;
		sum := 0;
		for row := 0 to nless1 do begin
				corner := corner + BytesPerLine;
				pp := ptr(corner);
				for column := 0 to nless1 do begin
						sum := sum + band(pp^,255) * kernel[k];
						k := k + 1;
						pp := ptr(ord(pp) + 1);
					end;
			end;
	end;
{$ELSEC}

{a0=^corner/^sum}
{a1=^kernel}
{a2=^pixels}

{d0=n-1}
{d1=BytesPerLine}
{d2=sum}
{d3=n-1(outer loop)}
{d4=n-1(inner loop)}
{d5=temp}

inline
	$4E56, $0000, {  link	a6,#0}
	$48E7, $FCE0,  {  movem.l	a0-a2/d0-d5,-(sp)}
	$4280,              {  clr.l	d0}
	$302E, $0012, {  move.w	18(a6),d0}
	$4281,              {  clr.l	d1}
	$322E, $0010, {  move.w	16(a6),d1}
	$206E, $000C, {  movea.l	12(a6),a0}
	$226E, $0004, {  movea.l	4(a6),a1}

	$4282,             {  clr.l	d2}
	$2600,             {  move.l	d0,d3}

	$D1C1,             {A adda.l	d1,a0}
	$2448,            {  move.l	a0,a2}
	$2800,            {  move.l	d0,d4}
	$4285,            {B clr.l	d5                   (2)}
	$1A1A,             {  move.b	(a2)+,d5    (6) }
	$CBD9,             {  muls	(a1)+,d5     (29!)}
	$D485,             {  add.l	d5,d2          (2)}
	$51CC, $FFF6, {  dbra	d4,B                (6)}
	$51CB, $FFEC, {  dbra	d3,A}

	$206E, $0008, {  move.l	8(a6),a0}
	$2082,              {  move.l	d2,(a0)}
	$4CDF, $073F, {  movem.l	(sp)+,a0-a2/d0-d5}
	$4E5E,              {  unlk	a6}
	$DEFC, $0010; {  add.w	#16,sp}
{$ENDC}



procedure DoConvolution (var kernel: ktype; n: integer);
	const
		skip = 7;
	var
		row, width, column, error: integer;
		margin, i, nless1: integer;
		frame, MaskRect, tRect: rect;
		AutoSelectAll, ScalingNeeded: boolean;
		SrcCenter, DstCenter, sum, max, offset, wsum, cscale, StartTicks: LongInt;
		value, MinResult, MaxResult: LongInt;
		p: ptr;
		str, str2: str255;
		ScaleFactor: extended;
		SaveRow:integer;
		NextUpdate: LongInt;
begin
	if NotinBounds or NotRectangular then
		exit(DoConvolution);
	StopDigitizing;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	SetupUndoFromClip;
	WhatToUndo := UndoFilter;
	frame := info^.RoiRect;
	with frame, Info^ do begin
			if ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction) then
				ApplyLookupTable;
			changes := true;
			margin := n div 2;
			if left < margin then
				left := left + margin;
			if right > (PicRect.right - margin) then
				right := right - margin;
			if top < margin then
				top := top + margin;
			if bottom > (PicRect.bottom - margin) then
				bottom := bottom - margin;
			SetPort(wptr);
			PenNormal;
			PenPat(AntPattern[PatIndex]);
			tRect := frame;
			OffscreenToScreenRect(tRect);
			FrameRect(tRect);
			width := right - left;
			max := n * n - 1;
			wsum := 0;
			for i := 0 to max do
				wsum := wsum + kernel[i];
			NumToString(n, str);
			NumToString(wsum, str2);
			InfoMessage := Concat(str, ' x ', str, ' kernel', crStr, 'sum = ', str2, crStr, crStr, CmdPeriodToStop);
			ShowInfo;
			if wsum <> 0 then
				cscale := wsum
			else
				cscale := 1;
			offset := -(n div 2) * BytesPerRow - BytesPerRow - n div 2;
			nless1 := n - 1;
			StartTicks := TickCount;
			str := '';
			SaveRow:=top;
			NextUpdate:=TickCount+6; {Update screen 10 times per second}
			if ScaleConvolutions then begin
					MinResult := MaxLongInt;
					MaxResult := -MaxLongInt;
					row := top;
					while row < bottom do begin
							SrcCenter := ord4(ClipBufInfo^.PicBaseAddr) + row * BytesPerRow + left;
							column := left;
							while column < (left + width) do begin
									DoOnePixel(nless1, BytesPerRow, SrcCenter + offset, sum, kernel);
									value := sum div cscale;
									if value < MinResult then
										MinResult := value;
									if value > MaxResult then
										MaxResult := value;
									SrcCenter := SrcCenter + skip;
									column := column + skip;
								end; {while column}
							row := row + skip;
						end; {while row...}
					ScalingNeeded := (MinResult < 0) or (MaxResult > 255);
					if ScalingNeeded then
						ScaleFactor := 253.0 / (MaxResult - MinResult)
					else
						ScaleFactor := 1.0;
					RealToString(ScaleFactor, 1, 4, str);
					str := concat('min=', long2str(MinResult), crStr, 'max=', long2str(MaxResult), crStr, 'scale factor= ', str);
					for row := top to bottom - 1 do begin
							SrcCenter := ord4(ClipBufInfo^.PicBaseAddr) + row * BytesPerRow + left;
							DstCenter := ord4(PicBaseAddr) + row * BytesPerRow + left;
							for column := left to left + width - 1 do begin
									DoOnePixel(nless1, BytesPerRow, SrcCenter + offset, sum, kernel);
									value := sum div cscale;
									if ScalingNeeded then begin
											if value < MinResult then
												value := MinResult;
											if value > MaxResult then
												value := MaxResult;
											value := round((value - MinResult) * ScaleFactor + 1);
										end;
									p := ptr(DstCenter);
									p^ := BAND(value, 255);
									SrcCenter := SrcCenter + 1;
									DstCenter := DstCenter + 1;
								end; {for column:=}
							if TickCount>=NextUpdate then begin
								SetRect(MaskRect, left, SaveRow, right, row + 1);
								UpdateScreen(MaskRect);
								SaveRow:=row+1;
								NextUpdate:=TickCount+6;
							end;
							if CommandPeriod then begin
									UpdatePicWindow;
									beep;
									exit(DoConvolution)
								end;
						end; {for row:=...}
				end  {Scale Convolutions}
			else
				for row := top to bottom - 1 do begin
						SrcCenter := ord4(ClipBufInfo^.PicBaseAddr) + row * BytesPerRow + left;
						DstCenter := ord4(PicBaseAddr) + row * BytesPerRow + left;
						for column := left to left + width - 1 do begin
								DoOnePixel(nless1, BytesPerRow, SrcCenter + offset, sum, kernel);
								value := sum div cscale;
								if value < MinResult then
									MinResult := value;
								if value > MaxResult then
									MaxResult := value;
								if value > 255 then
									value := 255;
								if value < 0 then
									value := 0;
								p := ptr(DstCenter);
								p^ := BAND(value, 255);
								SrcCenter := SrcCenter + 1;
								DstCenter := DstCenter + 1;
							end; {for column:=}
							if TickCount>=NextUpdate then begin
								SetRect(MaskRect, left, SaveRow, right, row + 1);
								UpdateScreen(MaskRect);
								SaveRow:=row+1;
								NextUpdate:=TickCount+6;
							end;
						if CommandPeriod then begin
								UpdatePicWindow;
								beep;
								exit(DoConvolution)
							end;
					end; {for row:=...}
			SetRect(MaskRect, left, SaveRow, right, bottom);
			UpdateScreen(MaskRect);
			ShowTime(StartTicks, frame, str);
		end; {with}
	SetupRoiRect;
	if AutoSelectAll then
		KillRoi;
end;


procedure Convolve (name: str255; RefNum: integer);
	var
		kernel: ktype;
		n, count: integer;
begin
	if noUndo then exit(Convolve);
	if name = '' then begin
			RefNum := 0;
			if not GetTextFile(name, RefNum) then
				exit(convolve)
			else
				KernelsRefNum := RefNum;
		end;
	DisableDensitySlice;
	GetKernel(kernel, n, name, RefNum);
	count := n * n;
	UpdatePicWindow;
	if (n >= 3) and (n <= 63) then
		DoConvolution(kernel, n)
	else
		PutError('Kernels must be n x n square matrices with 3 <= n <= 63.');
end;


procedure ConvolveUsingText;
	var
		f: integer;
		err: OSErr;
		count: LongInt;
begin
	if noUndo then exit(ConvolveUsingText);
	err := fsdelete('TempKernel', SystemRefNum);
	err := create('TempKernel', SystemRefNum, 'imag', 'TEXT');
	if err = NoErr then
		err := fsopen('TempKernel', SystemRefNum, f);
	if err <> NoErr then begin
			PutError('Unable to open temporary file.');
			exit(ConvolveUsingText);
		end;
	if TextInfo <> nil then
		with TextInfo^ do begin
				count := TextTE^^.TELength;
				err := fswrite(f, count, TextTE^^.hText^);
				err := fsclose(f);
				Convolve('TempKernel', SystemRefNum);
				err := fsdelete('TempKernel', SystemRefNum);
			end;
end;


procedure ApplyTableToLine (data: ptr; var table: LookupTable; width: LongInt);
{$IFC PowerPC}
	var
		line: LinePtr;
		i: integer;
begin
	line := LinePtr(data);
	for i := 0 to width - 1 do
		Line^[i] := table[Line^[i]];
end;
{$ELSEC}

{a0 = data}
{a1 = lookup table}
{d0 = width }
{d1 = pixel value}
inline
	$4E56, $0000, {  link a6,#0}
	$48E7, $C0C0, {  movem.l a0-a1/d0-d1,-(sp)}
	$206E, $000C, {  move.l 12(a6),a0}
	$226E, $0008, {  move.l 8(a6),a1}
	$202E, $0004, {  move.l 4(a6),d0}
	$5380,       {  subq.l #1,d0}
	$4281,       {  clr.l d1}
	$1210,       {L move.b (a0),d1}
	$10F1, $1000, {  move.b 0(a1,d1.w),(a0)+}
	$51C8, $FFF8, {  dbra d0,L}
	$4CDF, $0303, {  movem.l (sp)+,a0-a1/d0-d1}
	$4E5E,       {  unlk a6}
	$DEFC, $000C; {  add.w #12,sp}
{$ENDC}



function NewSurfacePlot: boolean;
	const
		WidthID = 5;
		HeightID = 6;
		TitleID = 8;
		WireframeID = 11;
		GrayscaleID = 12;
		DefaultName = 'Surface Plot';
	var
		mylog: DialogPtr;
		item: integer;
		SaveWidth, SaveHeight: integer;
		name: str255;
begin
	name := DefaultName;
	if not macro and not OptionKeyWasDown then begin
			InitCursor;
			SaveWidth := NewPicWidth;
			SaveHeight := NewPicHeight;
			mylog := GetNewDialog(190, nil, pointer(-1));
			SetDNum(MyLog, WidthID, NewPicWidth);
			SelectdialogItemText(MyLog, WidthID, 0, 32767);
			SetDNum(MyLog, HeightID, NewPicHeight);
			SetDString(MyLog, TitleID, DefaultName);
			SetDlogItem(MyLog, WireframeID, ord(WireframeSurfacePlots));
			SetDlogItem(MyLog, GrayscaleID, ord(GrayscaleSurfacePlots));
			OutlineButton(MyLog, ok, 16);
			repeat
				ModalDialog(nil, item);
				if item = WidthID then begin
						NewPicWidth := GetDNum(MyLog, WidthID);
						if (NewPicWidth < 0) or (NewPicWidth > MaxPicSize) then begin
								NewPicWidth := SaveWidth;
								SetDNum(MyLog, WidthID, NewPicWidth);
							end;
					end;
				if item = HeightID then begin
						NewPicHeight := GetDNum(MyLog, HeightID);
						if (NewPicHeight < 0) or (NewPicHeight > MaxPicSize) then begin
								NewPicHeight := SaveHeight;
								SetDNum(MyLog, HeightID, NewPicHeight);
							end;
					end;
				if item = WireframeID then begin
						WireframeSurfacePlots := not WireframeSurfacePlots;
						if not WireframeSurfacePlots then
							GrayscaleSurfacePlots := true;
						SetDlogItem(MyLog, WireframeID, ord(WireframeSurfacePlots));
						SetDlogItem(MyLog, GrayscaleID, ord(GrayscaleSurfacePlots));
					end;
				if item = GrayscaleID then begin
						GrayscaleSurfacePlots := not GrayscaleSurfacePlots;
						if not GrayscaleSurfacePlots then
							WireframeSurfacePlots := true;
						SetDlogItem(MyLog, WireframeID, ord(WireframeSurfacePlots));
						SetDlogItem(MyLog, GrayscaleID, ord(GrayscaleSurfacePlots));
					end;
			until (item = ok) or (item = cancel);
			if item = ok then
				name := GetDString(MyLog, TitleID);
			DisposeDialog(mylog);
			if NewPicWidth < 32 then
				NewPicWidth := 32;
			if NewPicHeight < 16 then
				NewPicHeight := 16;
			if item = cancel then begin
					NewPicWidth := SaveWidth;
					NewPicHeight := SaveHeight;
					NewSurfacePlot := false;
					exit(NewSurfacePlot);
				end;
		end; {if not macro}
	if not GrayscaleSurfacePlots then
		SetBackgroundColor(WhiteIndex);
	NewSurfacePlot := NewPicWindow(name, NewPicWidth, NewPicHeight);
	if GrayscaleSurfacePlots then
		SetBackgroundColor(WhiteIndex);
end;


procedure PlotSurface;
      {The code for grayscale/color surface plots was contributed by}
      {Norbert Vischer (norbert@mc.bio.uva.nl). Instead of erasing and framing}
      {a polygon, it  draws a section pattern into a buffer. The section pattern}
      {depends on calibration and LUT. 'CopyBits' is used to transfer the }
      {unframed , colored polygon into the surface plot window.}
	var
		hend, vend, h, v, DataWidth, DataHeight, i: LongInt;
		htemp, vtemp, ivalue: LongInt;
		skip, DataLeft, DataRight, DataTop, DataBottom: LongInt;
		hLoc, vLoc, hMin, hMax, vMin, vMax, MinIValue, MaxIValue: LongInt;
		hstart, vstart, dh, dv, hbase, vbase, vscale, nPlotLines, CalValue: extended;
		peak, MaxPeak, hinc, vinc, nLines, MinCalValue, MaxCalValue: extended;
		poly: PolyHandle;
		SaveInfo, PlotInfo: InfoPtr;
		aLine: LineType;
		MaskRect: rect;
		AutoSelectAll, ApplyLUT: boolean;
		table: LookupTable;
		StartTicks: LongInt;
		PolyRgn: RgnHandle;
		n, LineColor, pixelvalue: LongInt;
		sourcerect, destrect: rect;
		invers: boolean;
		penheight, dx, dy: LongInt;
		thisLevel, previousLevel: LongInt;
		SectionPatternInfo: Infoptr;
		FootPrintcolor, ignore: LongInt;
		SaveForeground, SaveBackground: LongInt;
		SaveGDevice: GDHandle;

	procedure FindVinc(var vstart, MaxPeak, vinc, hinc, dh, dv:extended); {ppc-bug}
	begin
		with PlotInfo^.PicRect do begin
				vstart := 5.0 + MaxPeak - dv * DataWidth;
				skip := round(DataHeight / ((bottom - vstart - 5.0) / vinc));
				if skip = 0 then
					skip := 1;
				nPlotLines := DataHeight / skip;
				vinc := (bottom - vstart - 5.0) / nPlotLines;
				vinc := vinc / 0.95;
				repeat
					vinc := vinc * 0.95;
					hinc := vinc / 2.0;
				until (5.0 + hinc * nPlotLines + dh * DataWidth) < right;
			end;
	end;

begin
	if NotRectangular or NotInBounds then
		exit(PlotSurface);
	StopDigitizing;
	DisableDensitySlice;
	SaveInfo := Info;
	SaveForeground := ForegroundIndex;
	SaveBackground := BackgroundIndex;
	if not NewSurfacePlot then begin
			KillRoi;
			exit(PlotSurface)
		end;
	PlotInfo := info;
	info := SaveInfo;
	AutoSelectAll := not Info^.RoiShowing;
	ShowWatch;
	if AutoSelectAll then
		SelectAll(true);
	if TooWide then
		exit(PlotSurface);
	with info^ do
		ApplyLUT := ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction);
	if ApplyLUT then
		GetLookupTable(table);
	Measure;
	UndoLastMeasurement(true);
	with results do begin
			MinIValue := MinIndex;
			MaxIValue := MaxIndex;
		end;
	if ApplyLut then begin
			MinIvalue := table[MinIValue];
			MaxIvalue := table[MaxIValue];
		end;
	MinCalValue := 10e100;
	MaxCalValue := -10e100;
	for i := MinIValue to MaxIValue do begin
			ivalue := i;
			if ApplyLUT then
				ivalue := table[ivalue];
			calValue := cvalue[i];
			if calValue < MinCalValue then
				MinCalValue := calValue;
			if calValue > MaxCalValue then
				MaxCalValue := calValue;
		end;
	WhatToUndo := NothingToUndo;
	with results do
		if (maxCValue - minCValue) <> 0.0 then
			vscale := (255.0 / (maxCValue - minCValue)) * 0.5
		else
			vscale := 0.5;
	with info^.RoiRect do begin
			DataLeft := left;
			DataRight := right;
			DataTop := top;
			DataBottom := bottom;
			DataWidth := DataRight - DataLeft;
			DataHeight := DataBottom - DataTop;
		end;
	dh := (0.65 * PlotInfo^.PicRect.right) / DataWidth;
	dv := -0.4 * dh;
	hstart := 5.0;
	vinc := 2.0;
	MaxPeak := (MaxCalValue - MinCalValue) * vscale * 0.5;
	FindVinc(vstart, MaxPeak, vinc, hinc, dh, dv); {First estimate}
	MaxPeak := MaxPeak * 2.0;
	hmin := DataRight + round(MaxPeak / dv);
	if hmin < 0 then
		hmin := 0;
	vmax := DataTop + round(MaxPeak / vinc);
	if vmax > DataBottom then
		vmax := DataBottom;
	MaxPeak := 0.0;
	vloc := DataTop;
	skip := 3;
	repeat
		hloc := hmin;
		repeat
			ivalue := MyGetPixel(hloc, vloc);
			if ApplyLUT then
				ivalue := table[ivalue];
			calValue := cvalue[ivalue];
			peak := (calValue - MinCalValue) * vscale + (DataRight - hloc) * dv - (vloc - DataTop) * vinc;
			if peak > MaxPeak then
				MaxPeak := peak;
			hloc := hloc + skip;
		until hloc > DataRight;
		vloc := vloc + skip;
	until vloc > vmax;
	FindVinc(vstart, MaxPeak, vinc, hinc, dh, dv);
	v := DataTop;
	SaveGDevice := GetGDevice;
	SetGDevice(osGDevice);
	SetPort(GrafPtr(PlotInfo^.osPort));
	if GrayscaleSurfacePlots then begin
			if not NewPicWindow('Section Pattern', NewPicWidth, NewPicHeight) then begin
					SetGDevice(SaveGDevice);
					exit(PlotSurface);
				end;
			SectionPatternInfo := Info;
			sourceRect := PlotInfo^.osPort^.portRect;
			SetPort(GrafPtr(SectionPatternInfo^.osPort));
			IndexToRgbForeColor(BlackIndex);
			IndexToRgbBackColor(WhiteIndex);
			PenNormal;
			previousLevel := sourcerect.bottom + 10; {add a safety band to the bottom}
			invers := cvalue[0] > cvalue[255];
			dx := sourcerect.right * 3 div 2;{>0.65}
			dy := dx * 2 div 5;{* 0.4}
			FootPrintcolor := -1;
			for n := 0 to 255 do begin{build inclined scale from bottom to top}
					if invers then
						pixelvalue := 255 - n
					else
						pixelvalue := n;
					if ApplyLUT then
						LineColor := table[pixelvalue]
					else
						LineColor := pixelvalue;
					thisLevel := round(sourcerect.bottom - vscale * (cvalue[LineColor] - MinCalValue));
					penheight := previousLevel - thisLevel;
					if penheight > 0 then begin{only if color was not dropped}
							MoveTo(0, thisLevel);
							IndexToRgbForeColor(pixelvalue);
							if FootPrintcolor = -1 then{catch the lowest color}
								FootPrintcolor := pixelvalue;
							pensize(1, penheight);
							line(dx, -dy);
							previousLevel := thisLevel;
						end;
				end;
			move(0, -10);
			pensize(1, 10);
			line(-dx, dy);{add a safety band to the top}
			PolyRgn := NewRgn;
		end; {if GrayscaleSurfacePlots}
	info := PlotInfo;
	SetGDevice(SaveGDevice);
	SelectWindow(PlotInfo^.wptr);
	SetGDevice(osGDevice);
	SetPort(GrafPtr(PlotInfo^.osPort));
	IndexToRgbForeColor(BlackIndex);
	IndexToRgbBackColor(WhiteIndex);
	PenNormal;
	UpdatePicWindow;
	StartTicks := TickCount;
	repeat
		hmax := 0;
		vmin := 9999;
		poly := OpenPoly;
		hbase := hstart;
		vbase := vstart;
		Info := SaveInfo;
		GetLine(DataLeft, v, DataWidth, aLine);
		info := PlotInfo;
		if ApplyLUT then
			ApplyTableToLine(@aLine, table, DataWidth);
		MoveTo(round(hbase), round(vbase - vscale * (cvalue[aLine[0]] - MinCalValue)));
		for i := 0 to DataWidth - 1 do begin
				hbase := hbase + dh;
				vbase := vbase + dv;
				hLoc := round(hbase);
				vLoc := round(vbase - vscale * (cvalue[aLine[i]] - MinCalValue));
				LineTo(hloc, vloc);
				if hloc > hmax then
					hmax := hloc;
				if vloc < vmin then
					vmin := vloc;
			end;
		LineTo(round(hbase), round(vbase));
		LineTo(round(hstart), round(vstart));
		LineTo(round(hstart), round(vstart - vscale * (cvalue[aLine[0]] - MinCalValue)));
		hmin := round(hstart);
		vmax := round(vstart);
		ClosePoly;
		if GrayscaleSurfacePlots then begin
				destrect := poly^^.polybbox;
                {now adjust size and position of sourcerect to get the right colors into the section : }
				sourcerect.right := destrect.right - destrect.left;
				sourcerect.top := sourcerect.bottom - (destrect.bottom - destrect.top);
				OpenRgn;
				FramePoly(poly);
				CloseRgn(PolyRgn);
				CopyBits(BitMapHandle(SectionPatternInfo^.osPort^.portPixMap)^^, BitMapHandle(info^.osPort^.portPixMap)^^, sourcerect, destrect, srcCopy, PolyRgn);
				if WireframeSurfacePlots then
					FramePoly(poly);
			end    {if GrayscaleSurfacePlots}
		else begin
				ErasePoly(poly);
				FramePoly(poly);
			end;
		KillPoly(poly);
		SetRect(MaskRect, hmin, vmin, hmax, vmax);
		UpdateScreen(MaskRect);
		hstart := hstart + hinc;
		vstart := vstart + vinc;
		v := v + skip;
	until (v >= DataBottom) or CommandPeriod;
	ShowTime(StartTicks, SaveInfo^.RoiRect, '');
	if CommandPeriod then
		beep;
	SetForegroundColor(SaveForeground);
	SetBackgroundColor(SaveBackground);
	info^.changes := true;
	if GrayscaleSurfacePlots then begin
			disposergn(PolyRgn);
			ignore := CloseAWindow(SectionPatternInfo^.wptr);
			info := PlotInfo;
		end;
	SetGDevice(SaveGDevice);
end;


{$POP}


function isqrt(x:LongInt):LongInt;
{The algorithm used by this integer square root function is described
in the document "square-root-code.txt" in the /pub/nih-image/documents
directory on zippy.nih.gov.} 
var
	r,nr,m:LongInt;
begin
    r:=0;
    m:=$40000000;
    repeat
    	nr := r + m;
        if (nr <= x) then begin
            x := x - nr;
            r := nr + m;
         end;
        r := bsr(r,1);
        m := bsr(m,2);
    until m = 0;
    if x > r then r:=r+1;
    isqrt:=r;
end;


procedure Filter (ftype: FilterType; pass: integer; var table: FateTable);
	var
		row, width, r1, r2, r3, c, value, error, sum, center: integer;
		tmp, SaveRow: integer;
		t1, t2, t3, t4: integer;
		MaskRect, frame, trect: rect;
		Filler1: integer;
		L1: LineType;
		Filler2: integer;
		L2: LineType;
		Filler3: integer;
		L3: LineType;
		Filler4: integer;
		result: LineType;
		a: SortArray;
		AutoSelectAll, UseMask, ProcessEdgePixels: boolean;
		L, T, R, B, index, code, FirstRow, LastRow: integer;
		StartTicks: LongInt;
		min, max, i, SaveBackground: integer;
		NextUpdate,tmp1,tmp2,tmp3:LongInt;
begin
	if NotinBounds then
		exit(Filter);
	StopDigitizing;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		with info^ do begin
				SelectAll(false);
				SetPort(wptr);
				PenNormal;
				PenPat(AntPattern[PatIndex]);
				FrameRect(wrect);
			end;
	if TooWide then
		exit(Filter);
	ShowWatch;
	if info^.RoiType <> RectRoi then
		UseMask := SetupMask
	else
		UseMask := false;
	if pass = 0 then begin
			SetupUndoFromClip;
			ShowMessage(CmdPeriodToStop);
			WhatToUndo := UndoFilter;
		end;
	frame := info^.RoiRect;
	StartTicks := TickCount;
	SaveBackground := BackgroundIndex;
	if ftype = MinRankFilter then
		BackgroundIndex := BlackIndex
	else
		BackgroundIndex := WhiteIndex;
	ProcessEdgePixels := ftype in [Erosion, Dilation, OutlineFilter, Skeletonize, Dither, MinRankFilter, MaxRankFilter];
	with frame, Info^ do begin
			changes := true;
			RoiShowing := false;
			width := right - left;
			if ProcessEdgePixels then begin
					FirstRow := top;
					LastRow := bottom - 1;
					L1[-1] := BackgroundIndex; {Out of bounds stores are intentional}
					L2[-1] := BackgroundIndex;
					L3[-1] := BackgroundIndex;
					if width < MaxLine then begin
							L1[width] := BackgroundIndex;
							L2[width] := BackgroundIndex;
							L3[width] := BackgroundIndex;
						end
				end
			else begin
					FirstRow := top + 1;
					LastRow := bottom - 2;
				end;
			GetLine(left, FirstRow - 1, width, L2);
			GetLine(left, FirstRow, width, L3);
			SaveRow := RoiRect.top;
			NextUpdate:=TickCount+6; {Update screen 10 times per second}
			for row := FirstRow to LastRow do begin
             {Move Convolution Window Down}
					BlockMove(@L2, @L1, width);
					BlockMove(@L3, @L2, width);
					GetLine(left, row + 1, width, L3);
                  {Process One Row}
					case ftype of
						FindEdges: 
							for c := 1 to width - 2 do begin
									tmp1 := L1[c - 1] + L1[c]*2 + L1[c + 1] - L3[c - 1] - L3[c]*2 - L3[c + 1];
									tmp2 := L1[c + 1] + L2[c + 1]*2 + L3[c + 1] - L1[c - 1] - L2[c - 1]*2 - L3[c - 1];
									tmp3:=isqrt(tmp1*tmp1+tmp2*tmp2);
									if tmp3>255 then tmp3:=255;
									result[c]:=tmp3;
								end;
						EdgeDetect: 
							for c := 1 to width - 2 do begin
									t1 := L1[c - 1] + L1[c] + L1[c + 1] - L3[c - 1] - L3[c] - L3[c + 1];
									t1 := abs(t1);
									t2 := L1[c + 1] + L2[c + 1] + L3[c + 1] - L1[c - 1] - L2[c - 1] - L3[c - 1];
									t2 := abs(t2);
									if t1 > t2 then
										tmp := t1
									else
										tmp := t2;
									if OptionKeyWasDown then begin
											if tmp > 255 then
												tmp := 255;
											if tmp < 0 then
												tmp := 0;
										end
									else if tmp > 35 then
										tmp := 255
									else
										tmp := 0;
									result[c] := tmp;
								end;
						ReduceNoise:  {median filter}
							for c := 0 to width - 1 do begin
									a[1] := L1[c - 1];
									a[2] := L1[c];
									a[3] := L1[c + 1];
									a[4] := L2[c - 1];
									a[5] := L2[c];
									a[6] := L2[c + 1];
									a[7] := L3[c - 1];
									a[8] := L3[c];
									a[9] := L3[c + 1];
									result[c] := FindMedian(a);
								end;
						MinRankFilter: 
							for c := 0 to width - 1 do begin
									min := 255;
									value := L1[c - 1];
									if value < min then
										min := value;
									value := L1[c];
									if value < min then
										min := value;
									value := L1[c + 1];
									if value < min then
										min := value;
									value := L2[c - 1];
									if value < min then
										min := value;
									value := L2[c];
									if value < min then
										min := value;
									value := L2[c + 1];
									if value < min then
										min := value;
									value := L3[c - 1];
									if value < min then
										min := value;
									value := L3[c];
									if value < min then
										min := value;
									value := L3[c + 1];
									if value < min then
										min := value;
									result[c] := min;
								end;
						MaxRankFilter: 
							for c := 0 to width - 1 do begin
									max := 0;
									value := L1[c - 1];
									if value > max then
										max := value;
									value := L1[c];
									if value > max then
										max := value;
									value := L1[c + 1];
									if value > max then
										max := value;
									value := L2[c - 1];
									if value > max then
										max := value;
									value := L2[c];
									if value > max then
										max := value;
									value := L2[c + 1];
									if value > max then
										max := value;
									value := L3[c - 1];
									if value > max then
										max := value;
									value := L3[c];
									if value > max then
										max := value;
									value := L3[c + 1];
									if value > max then
										max := value;
									result[c] := max;
								end;
						Dither:  {Floyd-Steinberg Algorithm}
							for c := 0 to width - 1 do begin
									value := L2[c];
									if value < 128 then begin
											result[c] := 0;
											error := -value;
										end
									else begin
											result[c] := 255;
											error := 255 - value
										end;
									tmp := L2[c + 1];              {A}
									tmp := tmp - (7 * error) div 16;
									if tmp < 0 then
										tmp := 0;
									if tmp > 255 then
										tmp := 255;
									L2[c + 1] := tmp;
									tmp := L3[c + 1];              {B}
									tmp := tmp - error div 16;
									if tmp < 0 then
										tmp := 0;
									if tmp > 255 then
										tmp := 255;
									L3[c + 1] := tmp;
									tmp := L3[c];              {C}
									tmp := tmp - (5 * error) div 16;
									if tmp < 0 then
										tmp := 0;
									if tmp > 255 then
										tmp := 255;
									L3[c] := tmp;
									tmp := L3[C - 1];                {D}
									tmp := tmp - (3 * error) div 16;
									if tmp < 0 then
										tmp := 0;
									if tmp > 255 then
										tmp := 255;
									L3[C - 1] := tmp;
								end;
						UnweightedAvg: 
							for c := 1 to width - 2 do begin
									tmp := (L1[C - 1] + L1[c] + L1[c + 1] + L2[c - 1] + L2[c] + L2[c + 1] + L3[c - 1] + L3[c] + L3[c + 1]) div 9;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						WeightedAvg: 
							for c := 1 to width - 2 do begin
									tmp := (L1[c - 1] + L1[c] + L1[c + 1] + L2[c - 1] + L2[c] * 4 + L2[c + 1] + L3[c - 1] + L3[c] + L3[c + 1]) div 12;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						fsharpen: 
							for c := 1 to width - 2 do begin
									if OptionKeyWasDown then
										tmp := L2[c] * 9 - L1[c - 1] - L1[c] - L1[c + 1] - L2[c - 1] - L2[c + 1] - L3[c - 1] - L3[c] - L3[c + 1]
									else begin
											tmp := L2[c] * 12 - L1[c - 1] - L1[c] - L1[c + 1] - L2[c - 1] - L2[c + 1] - L3[c - 1] - L3[c] - L3[c + 1];
											tmp := tmp div 4;
										end;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						SharpenMore: 
							for c := 1 to width - 2 do begin
										tmp := L2[c] * 9 - L1[c - 1] - L1[c] - L1[c + 1] - L2[c - 1] - L2[c + 1] - L3[c - 1] - L3[c] - L3[c + 1];
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						ShadowN: 
							for c := 1 to width - 2 do begin
									tmp:= L1[c-1] + L1[c]*2 + L1[c+1] + L2[c] - L3[c-1] - L3[c]*2 - L3[c+1];
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						ShadowNE: 
							for c := 1 to width - 2 do begin
									tmp:= L1[c] + L1[c+1]*2 - L2[c-1] + L2[c] + L2[c+1] - L3[c-1]*2 - L3[c];
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						ShadowE: 
							for c := 1 to width - 2 do begin
									tmp:= L1[c+1] -L1[c-1] - L2[c-1]*2 + L2[c] + L2[c+1]*2 - L3[c-1] + L3[c+1];
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						ShadowSE: 
							for c := 1 to width - 2 do begin
									tmp:= L3[c] + L3[c+1]*2 - L1[c-1]*2 - L1[c] - L2[c-1] + L2[c] + L2[c+1];
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						ShadowS: 
							for c := 1 to width - 2 do begin
									tmp:= L3[c-1] + L3[c]*2 + L3[c+1] - L1[c-1] - L1[c]*2 - L1[c+1] + L2[c];
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						ShadowSW: 
							for c := 1 to width - 2 do begin
									tmp:= L2[c-1] + L2[c] - L2[c+1] + L3[c-1]*2 + L3[c] - L1[c] - L1[c+1]*2;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						ShadowW: 
							for c := 1 to width - 2 do begin
									tmp:= L1[c-1] - L1[c+1] + L2[c-1]*2 + L2[c] - L2[c+1]*2 + L3[c-1] - L3[c+1];
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						ShadowNW: 
							for c := 1 to width - 2 do begin
									tmp:= L1[c-1]*2 + L1[c] + L2[c-1] + L2[c] - L2[c+1] - L3[c] - L3[c+1]*2;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						Erosion: 
							for c := 0 to width - 1 do begin
									center := L2[c];
									if center = BlackIndex then begin
											sum := L1[c - 1] + L1[c] + L1[c + 1] + L2[c - 1] + L2[c + 1] + L3[c - 1] + L3[c] + L3[c + 1];
											if (2040 - sum) >= BinaryThreshold then
												center := WhiteIndex;
										end;
									result[c] := center;
								end;
						Dilation: 
							for c := 0 to width - 1 do begin
									center := L2[c];
									if center = WhiteIndex then begin
											sum := L1[c - 1] + L1[c] + L1[c + 1] + L2[c - 1] + L2[c + 1] + L3[c - 1] + L3[c] + L3[c + 1];
											if sum >= BinaryThreshold then
												center := BlackIndex;
										end;
									result[c] := center;
								end;
						OutlineFilter: 
							for c := 0 to width - 1 do begin
									center := L2[c];
									if center = BlackIndex then begin
											if (L2[c - 1] = WhiteIndex) or (L1[c] = WhiteIndex) or (L2[c + 1] = WhiteIndex) or (L3[c] = WhiteIndex)
											or (L1[c-1] = WhiteIndex) or (L1[c+1] = WhiteIndex) or (L3[c-1] = WhiteIndex) or (L3[c+1] = WhiteIndex) then
												center := BlackIndex
											else
												center := WhiteIndex;
										end;
									result[c] := center;
								end;

						Skeletonize: 
							for c := 0 to width - 1 do begin
									center := L2[c];
									if center = BlackIndex then begin
											index := 0;											if L1[c - 1] = BlackIndex then
												index := bor(index, 1);
											if L1[c] = BlackIndex then
												index := bor(index, 2);
											if L1[c + 1] = BlackIndex then
												index := bor(index, 4);
											if L2[c + 1] = BlackIndex then
												index := bor(index, 8);
											if L3[c + 1] = BlackIndex then
												index := bor(index, 16);
											if L3[c] = BlackIndex then
												index := bor(index, 32);
											if L3[c - 1] = BlackIndex then
												index := bor(index, 64);
											if L2[c - 1] = BlackIndex then
												index := bor(index, 128);

											code := table[index];
											if odd(pass) then begin
													if (code = 2) or (code = 3) then begin
															center := WhiteIndex;
															PixelsRemoved := PixelsRemoved + 1;
														end;
												end
											else begin {even pass}
													if (code = 1) or (code = 3) then begin
															center := WhiteIndex;
															PixelsRemoved := PixelsRemoved + 1;
														end;
												end;
										end; {if}
									result[c] := center;
								end; {for}
					end; {case}
					if not ProcessEdgePixels then begin
							result[0] := L2[0];
							result[width - 1] := L2[width - 1];
						end;
					if UseMask then
						PutLineUsingMask(left, row, width, result)
					else
						PutLine(left, row, width, result);
					if TickCount>=NextUpdate then begin
							with RoiRect do
								SetRect(MaskRect, left, SaveRow, right, row+1);
							UpdateScreen(MaskRect);
							SaveRow:=row+1;
							if magnification > 1.0 then
								SaveRow := SaveRow - 1;
							NextUpdate:=TickCount+6;
							if CommandPeriod then begin
									UpdatePicWindow;
									beep;
									PixelsRemoved := 0;
									if AutoSelectAll then
										KillRoi;
									BackgroundIndex := SaveBackground;
									exit(filter)
								end;
						end;
				end; {for row:=...}
			trect := frame;
			InsetRect(trect, 1, 1);
			ShowTime(StartTicks, trect, '');
		end; {with}
	with frame do
		SetRect(MaskRect, left, SaveRow, right, bottom);
	UpdateScreen(MaskRect);
	SetupRoiRect;
	if AutoSelectAll then
		KillRoi;
	BackgroundIndex := SaveBackground;
end; {Filter}



procedure MakeSkeleton;
{This table-driven parallel thinning routine is based on an algorithm}
{by Zhang and Suen (CACM, March 1984, 236-239). There is}
{an entry in the table for each of the 256 possible 3x3 neighborhood}
{configurations. An entry of '1' means delete pixel on first pass, '2' means}
{delete pixel on second pass, and '3' means delete on either pass. There is a}
{routine in 'user.p' that draws all 256 neighborhoods.}
	const
		s999 = '01234567890123456789012345678901';
		s000 = '00010013003110130000000020203033';
		s032 = '00000000300000002000000020003022';
		s064 = '00000000000000000000000000000000';
		s096 = '20000000200020003000000030003020';
		s128 = '01310013000000010000000000000001';
		s160 = '31000000000000002000000000000000';
		s192 = '23130013000000010000000000000000';
		s224 = '2301000100000000330100002200200';
	var
		table: FateTable;
		s: str255;
		i, pass: integer;
begin
	s := concat(s000, s032, s064, s096, s128, s160, s192, s224);
	for i := 0 to 254 do
		table[i] := ord(s[i + 1]) - ord('0');
	table[255] := 0;
	pass := 0;
	repeat
		PixelsRemoved := 0;
		filter(skeletonize, pass, table);
		pass := pass + 1;
		if not CommandPeriod then
			filter(skeletonize, pass, table);
		pass := pass + 1;
	until (PixelsRemoved = 0) or CommandPeriod;
end;


procedure DoErosion;
	var
		i: integer;
		t: FateTable;
begin
	for i := 0 to BinaryIterations - 1 do begin
			filter(Erosion, i, t);
			if CommandPeriod then
				leave;
		end;
end;


procedure DoDilation;
	var
		i: integer;
		t: FateTable;
begin
	for i := 0 to BinaryIterations - 1 do begin
			filter(Dilation, i, t);
			if CommandPeriod then
				leave;
		end;
end;


procedure DoOpening;
	var
		i: integer;
		t: FateTable;
begin
	for i := 0 to BinaryIterations - 1 do begin
			filter(Erosion, i, t);
			if CommandPeriod then
				exit(DoOpening);
		end;
	for i := 0 to BinaryIterations - 1 do begin
			filter(Dilation, i + BinaryIterations, t);
			if CommandPeriod then
				exit(DoOpening);
		end;
end;

procedure DoClosing;
	var
		i: integer;
		t: FateTable;
begin
	for i := 0 to BinaryIterations - 1 do begin
			filter(Dilation, i, t);
			if CommandPeriod then
				exit(DoClosing);
		end;
	for i := 0 to BinaryIterations - 1 do begin
			filter(Erosion, i + BinaryIterations, t);
			if CommandPeriod then
				exit(DoClosing);
		end;
end;

procedure SetBinaryCount;
	var
		TempCount: integer;
		Canceled: boolean;
begin
	TempCount := GetInt('Neighborhood Pixel Count(1-8):', BinaryCount, Canceled);
	if Canceled then
		exit(SetBinaryCount);
	if (TempCount >= 1) and (TempCount <= 8) then begin
			BinaryCount := TempCount;
			BinaryThreshold := BinaryCount * 255
		end
	else
		beep;
end;

procedure SetIterations;
	var
		TempIterations: integer;
		Canceled: boolean;
begin
	TempIterations := GetInt('Number of Iterations:', BinaryIterations, Canceled);
	if Canceled then
		exit(SetIterations);
	if (TempIterations >= 1) and (TempIterations < 100) then
		BinaryIterations := TempIterations
	else
		beep;
end;


procedure ChangeValues (v1, v2, v3: integer);
      {Changes all pixels in the current selection with a value in the range v1 to v2 to a value of v3.}
	var
		i, value: integer;
		table: LookupTable;
begin
	for i := 0 to 255 do begin
			value := i;
			if (value >= v1) and (value <= v2) then
				value := v3;
			table[i] := value;
		end;
	ApplyTable(table);
end;


procedure DoPropagate (MenuItem: integer);
      {Copies the current Look-Up Table, spatial calibration, or density calibration to all open windows.}
	var
		TempInfo: InfoPtr;
		i: integer;

	procedure CopyLUTInfo;
	begin
		with info^ do begin
				TempInfo^.RedLUT := RedLUT;
				TempInfo^.GreenLUT := GreenLUT;
				TempInfo^.BlueLUT := BlueLUT;
				TempInfo^.ColorStart := ColorStart;
				TempInfo^.ColorEnd := ColorEnd;
				TempInfo^.nColors := nColors;
				TempInfo^.LutMode := LUTMode;
				TempInfo^.cTable := cTable;
				TempInfo^.FillColor1 := FillColor1;
				TempInfo^.FillColor2 := FillColor2;
				TempInfo^.FillColor1 := FillColor1;
				TempInfo^.SaveFill1 := SaveFill1;
				TempInfo^.SaveFill2 := SaveFill2;
			end;
	end;

	procedure CopySpatialCalibration;
		var
			SaveInfo: InfoPtr;
	begin
		with info^ do begin
				TempInfo^.xScale := xScale;
				TempInfo^.yScale := yScale;
				TempInfo^.PixelAspectRatio := PixelAspectRatio;
				TempInfo^.xUnit := xUnit;
				TempInfo^.changes := true;
				TempInfo^.SpatiallyCalibrated := SpatiallyCalibrated;
			end;
		SaveInfo := Info;
		Info := TempInfo;
		UpdateTitleBar;
		Info := SaveInfo;
	end;

	procedure CopyDensityCalibration;
		var
			SaveInfo: InfoPtr;
	begin
		with info^ do begin
				TempInfo^.ZeroClip := ZeroClip;
				TempInfo^.fit := fit;
				TempInfo^.nCoefficients := nCoefficients;
				TempInfo^.Coefficient := Coefficient;
				TempInfo^.UnitOfMeasure := UnitOfMeasure;
				TempInfo^.changes := true;
			end;
		SaveInfo := Info;
		Info := TempInfo;
		UpdateTitleBar;
		Info := SaveInfo;
	end;

begin
	for i := 1 to nPics do begin
			TempInfo := pointer(WindowPeek(PicWindow[i])^.RefCon);
			case MenuItem of
				1: 
					CopyLUTInfo;
				2: 
					CopySpatialCalibration;
				3: 
					CopyDensityCalibration;
			end; {case}
		end;
	WhatToUndo := NothingToUndo;
end;


procedure AutoThreshold;
       {Iterative thresholding technique, described originally by Ridler & Calvard in}
       {"PIcture Thresholding Using an Iterative Selection Method", IEEE transactions}
       { on Systems, Man and Cybernetics, August, 1978. }
	var
		AutoSelectAll, SaveRedirectFlag: boolean;
		index, MovingIndex, level: integer;
		tempSum1, tempSum2, tempSum3, tempSum4, result: extended;
begin
	AutoSelectAll := not info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	SaveRedirectFlag := RedirectSampling;
	RedirectSampling := false;
	if info^.RoiType = RectRoi then
		GetRectHistogram
	else
		GetHistogram;
	RedirectSampling := SaveRedirectFlag;
	OptionKeyWasDown := OptionKeyDown;
	if not OptionKeyWasDown then begin
              {Default is to set to these to null so erased areas won't be included in the threshold }
			Histogram[0] := 0;
			Histogram[255] := 0;
		end;
	with Results do begin {From ComputeResults}
			MinIndex := 0;
			while (histogram[MinIndex] = 0) and (MinIndex < 255) do
				MinIndex := MinIndex + 1;
			MaxIndex := 255;
			while (histogram[MaxIndex] = 0) and (MaxIndex > 0) do
				MaxIndex := MaxIndex - 1;
			if (MinIndex >= MaxIndex) then begin
					level := 128;
					ShowMessage(concat('Threshold=', Long2Str(level)));
					EnableThresholding(level);
					exit(AutoThreshold);
				end;
			MovingIndex := MinIndex;
			repeat
				tempSum1 := 0;
				tempSum2 := 0;
				tempSum3 := 0;
				tempSum4 := 0;
				for index := MinIndex to MovingIndex do begin
						tempSum1 := tempSum1 + index * Histogram[index];
						tempSum2 := tempSum2 + Histogram[index];
					end;
				for index := (MovingIndex + 1) to MaxIndex do begin
						tempSum3 := tempSum3 + index * Histogram[index];
						tempSum4 := tempSum4 + Histogram[index];
					end;
				Result := (tempSum1 / TempSum2 / 2.0) + (tempSum3 / tempSum4 / 2.0);
				MovingIndex := MovingIndex + 1;
			until ((MovingIndex + 1) > result) or (MovingIndex > (MaxIndex - 1));
			level := Round(result);
			EnableThresholding(level);
			ShowMessage(concat('Threshold=', Long2Str(level)));
		end; {with}
end;


procedure AutoDensitySlice;
	var
		AutoSelectAll: boolean;
		sigmak1k2, sigmax, nsum: extended;
		i, j, maxk1, maxk2, temp: integer;
		musubt, omegak1, omegak2, muk1, muk2: extended;
		part1, part2, part3: extended;
		intermed1, intermed2, intermed3: extended;
begin
	ResetGrayMap;
	AutoSelectAll := not info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	if info^.RoiType = RectRoi then
		GetRectHistogram
	else
		GetHistogram;
	maxk1 := 0;
	maxk2 := 0;
	musubt := 0.0;
	nsum := 0.0;
	for i := 1 to 254 do begin
			nsum := nsum + histogram[i];
		end;
	for i := 1 to 254 do begin
			musubt := musubt + (i * (histogram[i] / nsum));
		end;
	sigmak1k2 := 0.0;
	sigmax := 0.0;
	omegak1 := 0.0;
	muk1 := 0.0;
	for i := 1 to 253 do begin
			temp := i + 1;
			omegak2 := 0.0;
			muk2 := 0.0;
			omegak1 := omegak1 + (histogram[i] / nsum);
			muk1 := muk1 + (i * (histogram[i] / nsum));
			if omegak1 > 0.0 then begin
					for j := temp to 254 do begin
							omegak2 := omegak2 + (histogram[j] / nsum);
							muk2 := muk2 + (j * (histogram[j] / nsum));
							if omegak1 * omegak2 * (1.0 - omegak1 - omegak2) > 0.0 then begin
									part1 := ((omegak1 * muk2) - (omegak2 * muk1)) * ((omegak1 * muk2) - (omegak2 * muk1));
									intermed1 := omegak2 * omegak1;
									part2 := ((omegak1 * (musubt - muk2)) - (muk1 * (1 - omegak2))) * ((omegak1 * (musubt - muk2)) - (muk1 * (1 - omegak2)));
									intermed2 := omegak1 * (1 - omegak1 - omegak2);
									part3 := ((omegak2 * (musubt - muk1)) - (muk2 * (1 - omegak1))) * ((omegak2 * (musubt - muk1)) - (muk2 * (1 - omegak1)));
									intermed3 := omegak2 * (1 - omegak1 - omegak2);
									if intermed1 * intermed2 * intermed3 > 0.0 then begin
											sigmak1k2 := part1 / intermed1 + part2 / intermed2 + part3 / intermed3;
										end;
								end;
							if sigmak1k2 > sigmax then begin
									maxk1 := i;
									maxk2 := j;
									sigmax := sigmak1k2;
								end;
						end;
				end;
		end;
	SliceStart := maxk1;
	SliceEnd := maxk2;
end;

function GetScaleAndAngle: boolean;
	const
		hScaleID = 7;
		vScaleID = 8;
		AngleID = 9;
		NearestNeighborID = 10;
		BilinearID = 11;
		NewWindowID = 12;
	var
		mylog: DialogPtr;
		item, i: integer;
		vScaleUnchanged: boolean;
		str: str255;
begin
	vScaleUnchanged := true;
	InitCursor;
	mylog := GetNewDialog(50, nil, pointer(-1));
	SetDReal(MyLog, AngleID, rsAngle, 2);
	SetDReal(MyLog, hScaleID, rsHScale, 2);
	SelectdialogItemText(MyLog, hScaleID, 0, 32767);
	SetDReal(MyLog, vScaleID, rsVScale, 2);
	SetDlogItem(mylog, NewWindowID, ord(rsCreateNewWindow));
	SetDlogItem(mylog, BilinearID, ord(rsMethod = Bilinear));
	SetDlogItem(mylog, NearestNeighborID, ord(rsMethod = NearestNeighbor));
	repeat
		ModalDialog(nil, item);
		if item = AngleID then begin
				rsAngle := GetDReal(MyLog, AngleID);
				if rsAngle > 180.0 then
					rsAngle := 180.0;
				if rsAngle < -180.0 then
					rsAngle := -180.0;
			end;
		if item = hScaleID then begin
				str := GetDString(MyLog, hScaleID);
				rsHScale := StringToReal(str);
				if rsHScale = BadReal then
					rsHScale := 1.0;
				if vScaleUnchanged then begin
						rsVScale := rsHScale;
						SetDString(MyLog, vScaleID, str);
					end;
				if rsHScale < 0.05 then
					rsHScale := 0.05;
			end;
		if item = vScaleID then begin
				rsVScale := GetDReal(MyLog, vScaleID);
				if rsVScale = BadReal then
					rsVScale := 1.0;
				if rsVScale < 0.05 then
					rsVScale := 0.05;
				vScaleUnchanged := false;
			end;
		if item = NewWindowID then begin
				rsCreateNewWindow := not rsCreateNewWindow;
				SetDlogItem(mylog, NewWindowID, ord(rsCreateNewWindow));
			end;
		if (item = BilinearID) or (item = NearestNeighborID) then begin
				if item = BilinearID then
					rsMethod := Bilinear;
				if item = NearestNeighborID then
					rsMethod := NearestNeighbor;
				SetDlogItem(mylog, BilinearID, ord(rsMethod = Bilinear));
				SetDlogItem(mylog, NearestNeighborID, ord(rsMethod = NearestNeighbor));
			end;
	until (item = ok) or (item = cancel);
	DisposeDialog(mylog);
	GetScaleAndAngle := item <> cancel;
end;



procedure ScaleAndRotate;
	type
		EraseType = (Erase, DontErase);
	var
		value, DstWidth, DstHeight, hstart, vstart, hend, vend: integer;
		SaveWidth, SaveHeight: integer;
		hSrcCenter, vSrcCenter, hDstCenter, vDstCenter: integer;
		hbase, vbase, SrcWidth, SrcHeight: integer;
		SrcInfo, DstInfo, SaveInfo: InfoPtr;
		AutoSelectAll, UseNearestNeighbor, Rotate: boolean;
		MaskRect, SourceRect, DstRect: rect;
		StartTicks: LongInt;
		UseSameWindow: boolean;
		LinesPerUpdate: integer;

	procedure DoInterpolatedScaling;
    {Does interpolated scaling, but no rotation, using scaled integer arithmetic.}
		const
			CountsPerUpdate = 5;
		var
			SrcLeft, hloc, vloc, vbase, hbase, hrel: integer;
			LineCount, oldvloc, LastLine: integer;
			DstLine, SrcLine1, SrcLine2: LineType;
			MaskRect: rect;
			v, SrcTop: extended;
			h, hFraction, vFraction, UpperAverage, LowerAverage: LongInt;
			scale, scale2, hscale: LongInt;
	begin
		if (SrcWidth > MaxLine) or (DstWidth > MaxLine) then begin
			PutError(StringOf('Width must be less than ', MaxLine:1, ' pixels.'));
			exit(DoInterpolatedScaling);
		end;
		scale := 500;
		scale2 := scale * scale;
		hscale := round(rsHScale * scale);
		LastLine := SrcInfo^.PicRect.bottom - 1;
		with SourceRect do begin
				SrcLeft := left;
				SrcTop := top;
			end;
		with DstRect do begin
				oldvloc := top;
				LineCount := 0;
				for vloc := top to bottom - 1 do begin
						v := SrcTop + (vloc - top) / rsVScale;
						vbase := trunc(v);
						vFraction := round((v - vbase) * scale);
						Info := SrcInfo;
						GetLine(SrcLeft, vbase, SrcWidth, SrcLine1);
						SrcLine1[SrcWidth] := SrcLine1[SrcWidth - 1];
						if vbase <> LastLine then begin
								GetLine(SrcLeft, vbase + 1, SrcWidth, SrcLine2);
								SrcLine2[SrcWidth] := SrcLine2[SrcWidth - 1];
							end;
						for hloc := left to right - 1 do begin
								hrel := hloc - left;
								h := hrel * scale2 div hscale;
								hbase := hrel * scale div hscale;
								hFraction := h mod scale;
								LowerAverage := SrcLine1[hbase] + hFraction * (SrcLine1[hbase + 1] - SrcLine1[hbase]) div scale;
								UpperAverage := SrcLine2[hbase] + hFraction * (SrcLine2[hbase + 1] - SrcLine2[hbase]) div scale;
								DstLine[hrel] := (LowerAverage + vfraction * (UpperAverage - LowerAverage) div scale);
							end;
						Info := DstInfo;
						PutLine(left, vloc, DstWidth, DstLine);
						LineCount := LineCount + 1;
						if LineCount >= CountsPerUpdate then begin
								LineCount := 0;
								SetRect(MaskRect, left, oldvloc, right, vloc + 1);
								UpdateScreen(MaskRect);
								oldvloc := vloc;
							end;
						if CommandPeriod then begin
								beep;
								exit(DoInterpolatedScaling)
							end;
					end; {for vloc:=}
				SetRect(MaskRect, left, oldvloc, right, vloc + 1);
				UpdateScreen(MaskRect);
			end;
	end;


	procedure ScaleUsingCopyBits;
		var
			srcPort: cGrafPtr;
			SavePort: GrafPtr;
			MaskRect: rect;
			SaveGDevice: GDHandle;
	begin
		with DstInfo^ do begin
				SaveGDevice := GetGDevice;
				SetGDevice(osGDevice);
				GetPort(SavePort);
				SetPort(GrafPtr(osPort));
				IndexToRgbForeColor(BlackIndex);
				IndexToRgbBackColor(WhiteIndex);
				srcPort := SrcInfo^.osPort;
				CopyBits(BitMapHandle(srcPort^.portPixMap)^^, BitMapHandle(osPort^.PortPixMap)^^, SourceRect, DstRect, gCopyMode, nil);
				IndexToRgbForeColor(ForegroundIndex);
				IndexToRgbBackColor(BackgroundIndex);
				SetPort(SavePort);
				SetGDevice(SaveGDevice);
			end;
		if UseSameWindow then begin
				MaskRect := DstRect;
				UpdateScreen(MaskRect);
			end;
	end;

	procedure DoRotate;
		var
			vloc, hloc: integer;
			hRel, vRel: integer;
			DoScaling: boolean;
			scale, scale2: LongInt;
			AngleInRadians: extended;
			SinAngle, CosAngle, vSinAngle, vCosAngle: LongInt;
			hTemp, vTemp, h, v: LongInt;
			sHSrcCenter, sVSrcCenter: LongInt;
			hdiv, vdiv, hmod, vmod: integer;
			xbase, ybase: integer;
			LowerLeft, LowerRight, UpperLeft, UpperRight: integer;
			xfraction, yfraction, UpperAverage, LowerAverage: LongInt;
			offset, NextUpdate: LongInt;
			SaveVLoc:integer;
			MaskRect:rect;
	begin
		scale := 1000;
		scale2 := scale div 2;
		AngleInRadians := -((rsAngle + 270.0) / 360.0) * 2.0 * pi;
		SinAngle := round(sin(AngleInRadians) * scale);
		CosAngle := round(cos(AngleInRadians) * scale);
		sHSrcCenter := hSrcCenter * scale;
		sVSrcCenter := vSrcCenter * scale;
		DoScaling := (rsHScale <> 1.0) or (rsVScale <> 1.0);
		SaveVLoc:=vStart;
		NextUpdate:=TickCount+6; {Update screen 10 times per second}
		for vloc := vStart to vEnd do begin
				vrel := vloc - vDstCenter;
				vSinAngle := vrel * SinAngle;
				vCosAngle := vrel * CosAngle;
				for hloc := hStart to hEnd do begin
						hrel := hloc - hDstCenter;
						htemp := hrel * SinAngle + vCosAngle;
						vtemp := vSinAngle - hrel * CosAngle;
						if DoScaling then begin
								htemp := round(htemp / rsHScale);
								vtemp := round(vtemp / rsVScale);
							end;
						h := htemp + sHSrcCenter;
						v := vtemp + sVSrcCenter;
						info := SrcInfo;
						if UseNearestNeighbor then
							value := MyGetPixel((h + scale2) div scale, (v + scale2) div scale)
						else begin{Use bilinear interpolation}
								xbase := h div scale;
								ybase := v div scale;
								with info^ do
									if (xbase < 0) or (ybase < 0) or (xbase >= (PixelsPerLine - 1)) or (ybase >= (nlines - 1)) then
										value := 0
									else begin
											offset := ybase * BytesPerRow + xbase;
											LowerLeft := ImageP(PicBaseAddr)^[offset];
											LowerRight := ImageP(PicBaseAddr)^[offset + 1];
											UpperLeft := ImageP(PicBaseAddr)^[offset + BytesPerRow];
											UpperRight := ImageP(PicBaseAddr)^[offset + BytesPerRow + 1];
											xFraction := h mod scale;
											yFraction := v mod scale;
											UpperAverage := UpperLeft * scale + (xfraction * (UpperRight - UpperLeft));
											LowerAverage := LowerLeft * scale + (xfraction * (LowerRight - LowerLeft));
											value := ((LowerAverage + yfraction * ((UpperAverage - LowerAverage) + scale2) div scale) + scale2) div scale;
										end;
							end;
						Info := DstInfo;
						PutPixel(hloc, vloc, value);
					end; {for hloc:=}
				if TickCount>=NextUpdate then begin
						SetRect(MaskRect, hstart, SaveVLoc, hend, vloc + 1);
						UpdateScreen(MaskRect);
						ShowAnimatedWatch;
						SaveVLoc:=vloc+1;
						NextUpdate:=TickCount+6;
					end;
				if CommandPeriod then begin
						beep;
						{$ifc not PowerPC} {ppc-bug}
						KillRoi;
						exit(ScaleAndRotate)
						{$endc}
					end; {for hloc:=}
			end; {for vloc:=}
		SetRect(MaskRect, hstart, SaveVLoc, hend, vend + 1);
		UpdateScreen(MaskRect);
	end;


begin
	if NotRectangular or NotInBounds then
		exit(ScaleAndRotate);
	if not (macro and not rsInteractive) then
		if not GetScaleAndAngle then
			exit(ScaleAndRotate);
	UpdatePicWindow;
	UseSameWindow := not rsCreateNewWindow;
	if UseSameWindow then
		with info^ do
			if NoUndo then begin
					AbortMacro;
					exit(ScaleAndRotate)
				end;
	with info^ do
		UseNearestNeighbor := rsMethod = NearestNeighbor;
	DrawTools;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(true);
	ShowAnimatedWatch;
	if UseSameWindow then begin
			SetupUndo;
			WhatToUndo := UndoEdit;
			SetupUndoInfoRec;
			SrcInfo := UndoInfo;
			DstInfo := Info;
			if rsAngle = 0.0 then
				DoOperation(EraseOp);
		end
	else
		SrcInfo := info;
	with info^ do begin
			SourceRect := RoiRect;
			DstRect := RoiRect;
		end;
	with SourceRect do begin
			SrcWidth := right - left;
			SrcHeight := bottom - top;
			hSrcCenter := left + (SrcWidth div 2);
			vSrcCenter := top + (SrcHeight div 2);
			DstWidth := SrcWidth;
			DstHeight := SrcHeight;
		end;
	if UseSameWindow then
		with DstRect, info^ do begin
				if rsHScale <> 1.0 then begin
						DstWidth := round(SrcWidth * rsHScale);
						SaveWidth := DstWidth;
						left := left - (DstWidth - SrcWidth) div 2;
						if DstWidth > PicRect.right then
							DstWidth := PicRect.right;
						if left < 0 then
							left := 0;
						right := left + DstWidth;
						if DstWidth <> SaveWidth then begin
								SrcWidth := round(SrcWidth * (DstWidth / SaveWidth));
								SourceRect.left := hSrcCenter - SrcWidth div 2;
								SourceRect.right := SourceRect.left + SrcWidth;
							end;
					end;
				if rsVScale <> 1.0 then begin
						DstHeight := round(SrcHeight * rsVScale);
						SaveHeight := DstHeight;
						top := top - (DstHeight - SrcHeight) div 2;
						if DstHeight > PicRect.bottom then
							DstHeight := PicRect.bottom;
						if top < 0 then
							top := 0;
						bottom := top + DstHeight;
						if DstHeight <> SaveHeight then begin
								SrcHeight := round(SrcHeight * (DstHeight / SaveHeight));
								SourceRect.top := vSrcCenter - SrcHeight div 2;
								SourceRect.bottom := SourceRect.top + SrcHeight;
							end;
					end
			end {with}
	else begin
			DstWidth := round(SrcWidth * rsHScale);
			DstHeight := round(SrcHeight * rsVScale);
			if not NewPicWindow('Untitled', DstWidth, DstHeight) then begin
					KillRoi;
					exit(ScaleAndRotate)
				end;
			DstInfo := info;
			DstRect := info^.PicRect;
		end;
	with DstRect do begin
			hStart := left;
			vStart := top;
			hDstCenter := left + (DstWidth div 2);
			vDstCenter := top + (DstHeight div 2);
		end;
	hend := hstart + DstWidth - 1;
	vend := vstart + DstHeight - 1;
	rotate := rsAngle <> 0.0;
	if not macro then
		ShowMessage(CmdPeriodToStop);
	StartTicks := TickCount;
	if not rotate and (rsMethod = NearestNeighbor) then
		ScaleUsingCopyBits
	else if not rotate and not UseNearestNeighbor then
		DoInterpolatedScaling
	else
		DoRotate;
	if not macro then
		ShowTime(StartTicks, DstRect, '');
	KillRoi;
	with info^ do begin
			changes := true;
			if not UseSameWindow and (PixMapSize > UndoBufSize) then
				PutWarning;
			if SpatiallyCalibrated and (not UseSameWindow) then begin
					xScale := xScale * (DstWidth / SrcWidth);
					PixelAspectRatio := PixelAspectRatio * rsHScale / rsVScale;
					yScale := xScale / PixelAspectRatio;
				end;
		end;
	if not UseSameWindow and AutoSelectAll then begin
			SaveInfo := Info;
			Info := SrcInfo;
			KillRoi;
			Info := SaveInfo;
		end;
	if UseSameWindow then
		with NoInfo^ do begin
				roiType := RectRoi;
				RoiRect := DstRect;
				RectRgn(roiRgn, DstRect);
			end;
end;


function DoRankDialogBox: boolean;
	const
		MedianID = 5;
		MinID = 6;
		MaxID = 7;
		OpenID = 8;
		CloseID = 9;
		IterationsID = 11;
	var
		mylog: DialogPtr;
		item: integer;

	procedure SetButtons;
	begin
		SetDlogItem(mylog, MedianID, ord(RankFilter = MedianRank));
		SetDlogItem(mylog, MinID, ord(RankFilter = MinRank));
		SetDlogItem(mylog, MaxID, ord(RankFilter = MaxRank));
		SetDlogItem(mylog, OpenID, ord(RankFilter = OpeningRank));
		SetDlogItem(mylog, CloseID, ord(RankFilter = ClosingRank));
	end;

begin
	InitCursor;
	mylog := GetNewDialog(260, nil, pointer(-1));
	SetButtons;
	SetDNum(mylog, IterationsID, RankIterations);
	SelectdialogItemText(mylog, IterationsID, 0, 32767);
	repeat
		ModalDialog(nil, item);
		if (item >= MedianID) and (item <= CloseID) then begin
				if item = MedianID then
					RankFilter := MedianRank
				else if item = MinID then
					RankFilter := MinRank
				else if item = MaxID then
					RankFilter := MaxRank
				else if item = OpenID then
					RankFilter := OpeningRank
				else
					RankFilter := ClosingRank;
				SetButtons;
			end;
		if item = IterationsID then begin
				RankIterations := GetDNum(mylog, IterationsID);
				if RankIterations < 1 then
					RankIterations := 1;
				SetDNum(mylog, IterationsID, RankIterations);
			end;
	until (item = ok) or (item = cancel);
	DisposeDialog(mylog);
	DoRankDialogBox := item = ok;
end;


procedure DoRankFilter;
	var
		t: FateTable;
		i: integer;
begin
	if not macro then
		if not DoRankDialogBox then
			exit(DoRankFilter);
	case RankFilter of
		MedianRank: 
			for i := 0 to RankIterations - 1 do begin
					Filter(ReduceNoise, i, t);
					if CommandPeriod then
						leave;
				end;
		MinRank: 
			for i := 0 to RankIterations - 1 do begin
					Filter(MinRankFilter, i, t);
					if CommandPeriod then
						leave;
				end;
		MaxRank: 
			for i := 0 to RankIterations - 1 do begin
					Filter(MaxRankFilter, i, t);
					if CommandPeriod then
						leave;
				end;
		OpeningRank:  begin
				for i := 0 to RankIterations - 1 do begin
						Filter(MinRankFilter, i, t);
						if CommandPeriod then
							leave;
					end;
				for i := RankIterations to RankIterations * 2 - 1 do begin
						Filter(MaxRankFilter, i, t);
						if CommandPeriod then
							leave;
					end;
			end;
		ClosingRank:  begin
				for i := 0 to RankIterations - 1 do begin
						Filter(MaxRankFilter, i, t);
						if CommandPeriod then
							leave;
					end;
				for i := RankIterations to RankIterations * 2 - 1 do begin
						Filter(MinRankFilter, i, t);
						if CommandPeriod then
							leave;
					end;
			end
	end;
end;


procedure DoShadowFilter;
	var
		mylog: DialogPtr;
		item: integer;
		t: FateTable;
begin
	mylog := GetNewDialog(270, nil, pointer(-1));
	repeat
		ModalDialog(nil, item);
	until ((item >=1) and (item<=9));
	DisposeDialog(mylog);
	case item of
		1: Filter(ShadowN,  0, t);
		2: Filter(ShadowNE, 0, t);
		3: Filter(ShadowE,  0, t);
		4: Filter(ShadowSE, 0, t);
		5: Filter(ShadowS,  0, t);
		6: Filter(ShadowSW, 0, t);
		7: Filter(ShadowW,  0, t);
		8: Filter(ShadowNW, 0, t);
		otherwise;
	end;
end;


end.