{This archive contains modified version of NIH-Image, where the following}
{functions have been added:}
{                                                                           }
{- XYQuantile, which is an extension of the median filtering, also }
{  allowing max or min filtering. It works on 2 dim images.}
{- ZQuantile, which does the same filtering "vertically" on a stack,}
{- Crest Pathway, which can be used to automatically trace contours}
{- Assisted Tracking: a routine for semi-automated contour tracing.}
{                                                                           }
{XYQuantile   }
{---------- }
{  The XYQuantile routine allows you to effect a median filtering on a      }
{rectangle of chosen dimensions.  The user is asked to specify the values   }
{of  filter width (FW), filter height (FH), and quantile (Q).  For each     }
{pixel, the routine selects a rectangle of FWxFH pixels, sorts all the      }
{points in this rectangle,  and replaces the original pixel by the          }
{quantile. The median value is given by Q=((FWxFH) div 2) +1. For           }
{instance, in a 5x3 rectangle, the median is the 8th quantile of the        }
{sorted values.   }
{  To show the potenial of this filter, I added a "Noisy Image". Compare    }
{the filtering of this image with  the Reduce Noise function if NIH, and    }
{with a median filter with the following parameters:  Filter Width=5,       }
{Filter Height=1, Quantile=3.   }
{   You can also specify that the routine retains another value than the    }
{median, for instance the local  minimum (Quantile=1) or maximum            }
{(Quantile=FHxFW).                                                          }
{                                                                           }
{Macro call: UserCode(1, FilterW, FilterH, Quantile);}
{                                                                           }
{                                                                           }
{ZQuantile   }
{---------   }
{  The ZQuantile filter works exclusively on stacks. It perfomrs            }
{basically the same operation than the XYQuantile  procedure, but           }
{vertically. It gives similar results as the Average function of the        }
{Stack menu, but there  is a weaker influence of the extreme values on      }
{the median than on the average.   }
{                                                                           }
{  Note that both XYQuantile and ZQuantile filterings are slow to           }
{process, and that the speed is inversely propotional to the size of the    }
{filter.   }
{                                                                           }
{Macro call:  UserCode(2,Quantile,0,0);}
{                                                                           }
{                                                                           }
{Crest Pathway:}
{-------------- }
{  This function allows you to trace a crest in an image by reiterately }
{stepping from a point the neighbourg, following the darkest direction.}
{Use the "Crest/Track Options" command to optimize the path finding }
{parameters.}
{                                                                           }
{Macro call: }
{  UserCode(3,0,0,0);}
{or:}
{  UserCode('CrestPathway',0,0,0);}
{                                                                           }
{                                                                           }
{Assisted Tracking:}
{------------------}
{  This function facilitates the manual surrouding of structures with a }
{dark edge. The principle is that you follow roughly the border of your }
{structure with the mouse, and the routines follows your mouse by the }
{darkest pathway. It is in fact a variant of the Crest Pathway routine. }
{Use "Crest/Track Options" to optimize the tracking parameters. }
{                                                                           }
{Macro call: }
{  Usercode(4,0,0,0);}
{or:}
{  UserCode('AssistedTracking',0,0,0);}
{                                                                           }
{                                                                           }
{Crest/Track Options:}
{-------------------- }
{Allows the user to choose parameters for the "Crest Pathway" and }
{"Assisted Tracking" commands.}
{                                                                           }
{1) Maximal path length: indicates the number of steps the routine will }
{   reiterate. Note that the routine stops earlier if the pathway is cycling.  }
{   A value of 0 for the maximal pathway indicates that the routine only stops }
{   after having accomplished a cycling pathway. }
{   Default value=0. }
{                                                                           }
{2) Initial direction: allows the user to choose the direction to follow }
{   just after clicking. Direction is measured in degrees (0-360).}
{   Range:    0-360}
{   Default:  90 (towards the top of the image)}
{  }
{3) Maximal angle: this parameters indicates the maximal deviation that }
{   the path can take at each step. Values are rounded to the nearest multiple }
{   of 45 degrees.}
{   Accepted values : 45, 90 or 135}
{   Default : 45}
{                                                                           }
{4) Averaging distance: to avoid influence of noisy pixels, you can ask }
{   the routine to average the gray level over a chosen distance, and to }
{   choose the direction with the darkest average. }
{   Range : 1-15}
{   Default value : 3. }
{                                                                              }
{5) Assistance radius: this option is only valid for the "ASsisted Tracking"   }
{   command. It represents the maximal distance (in pixels) between the        }
{   mouse and the selected path.}
{   Range : 1-45}
{   Default value=15. In practice, the best results are obtained with a value}
{   of 3 to 5 times the averaging distance.}
{                                                                              }
{6) Tail Eating:}
{   Allows you to eliminate the part of the pathway outside the cycle          }
{   (if your first clicking was not exactly within the circle).                }
{   Default value : false (no tail eating).}
{                                                                              }
{                                                                              }
{Protocol for surrounding fluorescent cells:}
{-------------------------------------------}
{I use the following protocol to surround fluorescently labeled cells:}
{1) XYQuantile (5,5,3) for smoothing the image without affecting edges. }
{2) Trace Edge function of NIH-Image, but with the Option key helt down, }
{   to avoid thresholding. }
{3) Crest Path to automatically surround the cells and select the border.}
{   PŠrameter values are 0 for Max Path Lenght, 45 for Maximal Angle, 3 for}
{   Averaging Distance, and true for Eating Tail. Initial Direction is variable}
{   and Assistance radius does not matter.}
{4) Measure.}
{                                                                           }
{An example stack is joined to the archive: "Find cells stack".}
{                                                                           }
{If you find some use of or have some comment about these commands,   }
{please contact jvanheld@dbm.ulb.ac.be   }



unit User;


interface

	uses
		File1, NI_DAQ_MAC, DigitalIO, Camera, ImageTimers, QuickDraw, Palettes, PrintTraps, Globals, Utilities, Graphics, {Functions, }
		Filters, Analysis, Initialization;

	procedure InitUser;
	procedure DoUserMenuEvent (MenuItem: integer);
	procedure TestFindPeak (event: EventRecord);
	procedure OldUserMacroCode (CodeNumber: integer; Param1, Param2, Param3: extended);
	procedure UserMacroCode (str: str255; Param1, Param2, Param3: extended);

implementation

{User global variables go here.}
	const
		XYQuantileItem = 1;
		ZQuantileItem = 2;
		CrestPathwayItem = 3;
		AssistedTrackingItem = 4;
		CrestOptionsItem = 5;


{CrestPathway}
		MaxAveragingDist = 15;
		MaxAssistanceRadius = 45;
		MaxMaxAngle = 3;
		MaxMaxCount = 10000;
		DefAveragingDist = 3;
		DefMaxAngle = 1;
		DefInitialDir = 2;
		DefMaxCount = 0;


	var
		PeakRadius, Peakedness: extended;
		color, MinSpacing: integer;

{CrestPathway}
		AssistanceRadius, AveragingDist, MaxAngle, InitialDir, MaxCount: integer;
		TailEating, Canceled: boolean;



	procedure InitUser;
		var
			i: integer;
	begin
		UserMenuH := GetMenu(UserMenu);
		InsertMenu(UserMenuH, 0);
		DrawMenuBar;

{CrestPathway}
		AveragingDist := DefAveragingDist;
		AssistanceRadius := 5 * AveragingDist;
		MaxAngle := DefMaxAngle;
		InitialDir := DefInitialDir;
		MaxCount := DefMaxCount;
		TailEating := false;
	end;



{jvh}
	procedure CrestOptions;
		const
			AveragingDistID = 3;
			AssistanceRadiusID = 4;
			MaxAngleID = 5;
			InitialDirID = 6;
			MaxCountID = 7;
			TailEatingID = 8;
		var
			mylog: DialogPtr;
			item, i: integer;
			SaveAveragingDist, SaveAssistanceRadius, SaveMaxAngle, SaveInitialDir, SaveMaxCount: integer;
			SaveTailEating: boolean;

	begin
		InitCursor;
		SaveAveragingDist := AveragingDist;
		SaveAssistanceRadius := AssistanceRadius;
		SaveMaxAngle := MaxAngle;
		SaveInitialDir := InitialDir;
		SaveMaxCount := MaxCount;
		SaveTailEating := TailEating;
		MaxAngle := MaxAngle * 45;
		InitialDir := (540 - InitialDir * 45) mod 360;

		mylog := GetNewDialog(7000, nil, pointer(-1));
		SetDNum(MyLog, AveragingDistID, AveragingDist);
		SetDNum(MyLog, AssistanceRadiusID, AssistanceRadius);
		SetDNum(MyLog, MaxAngleID, MaxAngle);
		SetDNum(MyLog, InitialDirID, InitialDir);
		SetDNum(MyLog, MaxCountID, MaxCount);
		SetDNum(MyLog, MaxAngleID, MaxAngle);
		OutlineButton(MyLog, OK, 16);
		if TailEating then
			SetDialogItem(mylog, TailEatingID, 1);

		repeat
			ModalDialog(nil, item);
			case item of
				AveragingDistID: 
					begin
						AveragingDist := GetDNum(MyLog, AveragingDistID);
						SetDNum(MyLog, AveragingDistID, AveragingDist)
					end;
				AssistanceRadiusID: 
					begin
						AssistanceRadius := GetDNum(MyLog, AssistanceRadiusID);
						SetDNum(MyLog, AssistanceRadiusID, AssistanceRadius)
					end;
				MaxAngleID: 
					begin
						MaxAngle := GetDNum(MyLog, MaxAngleID);
						SetDNum(MyLog, MaxAngleID, MaxAngle)
					end;
				InitialDirID: 
					begin
						InitialDir := GetDNum(MyLog, InitialDirID);
						SetDNum(MyLog, InitialDirID, InitialDir)
					end;
				MaxCountID: 
					begin
						MaxCount := GetDNum(MyLog, MaxCountID);
						SetDNum(MyLog, MaxCountID, MaxCount)
					end;
				TailEatingID: 
					begin
						TailEating := not TailEating;
						SetDialogItem(mylog, TailEatingID, ord(TailEating));
					end;
				otherwise
					;
			end;
		until (item = ok) or (item = cancel);
		DisposDialog(mylog);
		if item = cancel then
			begin
				AveragingDist := SaveAveragingDist;
				AssistanceRadius := SaveAssistanceRadius;
				MaxAngle := SaveMaxAngle;
				InitialDir := SaveInitialDir;
				MaxCount := SaveMaxCount;
				TailEating := SaveTailEating;
			end
		else
			begin
				if AveragingDist > MaxAveragingDist then
					AveragingDist := MaxAveragingDist
				else if AveragingDist < 1 then
					AveragingDist := 1;

				if AssistanceRadius > MaxAssistanceRadius then
					AssistanceRadius := MaxAssistanceRadius
				else if AssistanceRadius < 1 then
					AssistanceRadius := 1;

				MaxAngle := round(MaxAngle / 45);
				if MaxAngle > MaxMaxAngle then
					MaxAngle := MaxMaxAngle
				else if MaxAngle < 1 then
					MaxAngle := 1;

				InitialDir := (round((540 - InitialDir) / 45)) mod 8;
				if InitialDir > 7 then
					InitialDir := 7
				else if InitialDir < 0 then
					InitialDir := 0;

				if MaxCount > MaxMaxCount then
					MaxCount := MaxMaxCount
				else if MaxCount < 0 then
					MaxCount := 0;
			end;
	end;




	procedure CrestPathway (AssistedTracking: boolean);
		const
			DeltaSignif = 0.001;
		var
			dx, dy: array[0..7] of integer;
			NewPicNum, OldPicNum, ExAequos, i, dir, OptimalDir, counter, dist, x, y: integer;
			CurrentDir, xdelta, ydelta, FirstPointInLoop, SaveBackGround, ImgH, ImgW: integer;
			Interrupted, MaxCountReached, OnABorder, Cycled: boolean;
			InitialAngle, MouseDist, avg, max: real;
			MousePt: point;
			Oldinfo, Newinfo: InfoPtr;
			TempRgn: RgnHandle;
			tPort: GrafPtr;

	begin
		if info^.PictureType = NullPicture then
			exit(CrestPathway);
		Oldinfo := Info;
		OldPicNum := Info^.PicNum;
		ImgH := Info^.nlines;
		ImgW := Info^.PixelsPerLine;

		SaveBackground := BackgroundIndex;
		BackgroundIndex := WhiteIndex;
		MaxCountReached := false;
		OnABorder := false;
		Cycled := false;
		counter := 1;
		FirstPointInLoop := 0;
		CurrentDir := InitialDir;
		Interrupted := false;
		dx[0] := -1;
		dy[0] := 0;
		dx[1] := -1;
		dy[1] := -1;
		dx[2] := 0;
		dy[2] := -1;
		dx[3] := 1;
		dy[3] := -1;
		dx[4] := 1;
		dy[4] := 0;
		dx[5] := 1;
		dy[5] := 1;
		dx[6] := 0;
		dy[6] := 1;
		dx[7] := -1;
		dy[7] := 1;

		if (AveragingDist < 1) or (AveragingDist > MaxAveragingDist) then
			CrestOptions;

		if (MaxAngle < 1) or (MaxAngle > MaxMaxAngle) then
			CrestOptions;

		if (InitialDir < 0) or (InitialDir > 7) then
			CrestOptions;

		if (MaxCount < 0) or (MaxCount > MaxMaxCount) then
			CrestOptions;

		if not NewPicWindow('Crest pathway', ImgW, ImgH) then
			begin
				PutMessage('Not enough memory to open destination image');
				BackgroundIndex := SaveBackground;
				exit(CrestPathway);
			end;
		NewInfo := Info;
		NewPicNum := Info^.PicNum;

		Info := OldInfo;
		SelectWindow(PicWindow[OldPicNum]);

		with Info^ do
			begin
{get the initial position}
				SetPort(wptr);
				PenNormal;
				PenPat(pat[PatIndex]);
				PenSize(1, 1);
				ShowMessage('Click on starting position');

				SetCursor(ToolCursor[SelectionTool]);
				RoiShowing := false;
				ActivateWindow;
				UpdatePicWindow;
				SetUpUndo;
				repeat
				until button;
				GetMouse(MousePt);
				x := MousePt.h;
				y := MousePt.v;

				if (x < 0) or (y < 0) or (x > ImgW - 1) or (y > ImgH - 1) then {point outside the image}
					begin
						Putmessage(concat('Click on starting porition inside ', title, ' image'));
						UpdatePicWindow;
						repeat
						until button;
						GetMouse(MousePt);
						x := MousePt.h;
						y := MousePt.v;
						if (x < 0) or (y < 0) or (x > ImgW - 1) or (y > ImgH - 1) then {point outside the image}
							exit(CrestPathway);
					end;

				xCoordinates^[counter] := x;
				yCoordinates^[counter] := y;
				nCoordinates := counter;
				MoveTo(x, y);

{find pathway}
				repeat
					Info := OldInfo;
					counter := counter + 1;
					max := -1;
					for i := (CurrentDir - MaxAngle) to (CurrentDir + MaxAngle) do
						begin
							dir := (i + 8) mod 8;
							avg := 0;
							for dist := 1 to AveragingDist do
								avg := avg + myGetPixel(x + dist * dx[dir], y + dist * dy[dir]);
							avg := avg / AveragingDist;
							if (avg > max) or (abs(avg - max) <= DeltaSignif * avg) then
								begin
									if (abs(avg - max) <= DeltaSignif * avg) then
										ExAequos := ExAequos + 1
									else
										ExAequos := 0;
									max := avg;
									OptimalDir := dir;
								end;
						end;	{for i := (CurrentDir - MaxAngle) to (CurrentDir + MaxAngle)}
					if ExAequos > 1 then
						OptimalDir := (OptimalDir - (ExAequos div 2) + 8) mod 8;

{next position}
					x := x + dx[OptimalDir];
					y := y + dy[OptimalDir];
					if x < 0 then
						x := 0
					else if x > ImgW - 1 then
						x := ImgW - 1;
					if y < 0 then
						y := 0
					else if y > ImgH - 1 then
						y := ImgH - 1;
					xCoordinates^[counter] := x;
					yCoordinates^[counter] := y;
					nCoordinates := counter;
					LineTo(x, y);

{Assisted tracking}
					if AssistedTracking then
						begin
							repeat
								GetMouse(MousePt);
								xdelta := MousePt.h - x;
								ydelta := MousePt.v - y;
								if abs(xdelta) > MaxAssistanceRadius then
									if xdelta > 0 then
										xdelta := MaxAssistanceRadius
									else
										xdelta := -MaxAssistanceRadius;
								if abs(ydelta) > MaxAssistanceRadius then
									if ydelta > 0 then
										ydelta := MaxAssistanceRadius
									else
										ydelta := -MaxAssistanceRadius;
								MouseDist := sqrt(sqr(xdelta) + sqr(ydelta));
							until (MouseDist > AssistanceRadius) or not button;
							if not button then
								Interrupted := true
							else
								begin
									InitialAngle := abs(180 / 3.1416 * arctan(ydelta / xdelta));
									if (InitialAngle < 23) then
										if xdelta > 0 then
											CurrentDir := 4
										else
											CurrentDir := 0
									else if InitialAngle > 66 then
										if ydelta > 0 then
											CurrentDir := 6
										else
											CurrentDir := 2
									else if ydelta > 0 then
										if xdelta > 0 then
											CurrentDir := 5
										else
											CurrentDir := 7
									else if xdelta > 0 then
										CurrentDir := 3
									else
										CurrentDir := 1;
								end;
						end
					else
						CurrentDir := OptimalDir;

					if (MaxCount > 0) and (counter >= MaxCount) then
						MaxCountReached := true;

					if ((x < 1) or (y < 1) or (x > ImgW - 2) or (y > ImgH - 2)) then
						OnABorder := true;

		{next pixel already in the path ?}
					Info := NewInfo;
					if (MyGetPixel(x, y) <> 0) and (counter > 10) and not OnABorder then
						begin
							Cycled := true;
							Info := NewInfo;
							UpdatePicWindow;
						end
					else
						PutPixel(x, y, 255);

				until MaxCountReached or Cycled or (OnABorder and not AssistedTracking) or Interrupted;
			end; {with Info^}

		if AssistedTracking and Cycled then
			begin
				ShowMessage('Cycle completed');
				repeat
				until not button;
			end;

		if Cycled and TailEating then {removing points outside of the loop}
			begin
{find first point in the loop}
				Info := NewInfo;
				repeat
					FirstPointInLoop := FirstPointInLoop + 1;
				until ((xCoordinates^[FirstPointInLoop] = x) and (yCoordinates^[FirstPointInLoop] = y)) or (FirstPointInLoop > counter);

{blank all points outside the loop (eats the tail)}
				if (FirstPointInLoop < Counter) and (FirstPointInLoop > 1) then
					begin
						for i := 0 to FirstPointInLoop - 1 do
							PutPixel(xCoordinates^[i], yCoordinates^[i], 0);
						for i := 0 to Counter - FirstPointInLoop + 1 do
							begin
								xCoordinates^[i] := xCoordinates^[i + FirstPointInLoop - 1];
								yCoordinates^[i] := yCoordinates^[i + FirstPointInLoop - 1];
							end;
						nCoordinates := Counter - FirstPointInLoop + 1;
					end;
			end;

		if Cycled then
			if MaxCount = 0 then
				ShowMessage(concat('Cycle of ', Long2str(Counter - FirstPointInLoop), ' steps'))
			else
				ShowMessage(concat('Cycled after ', Long2str(Counter), ' steps'))
		else if OnABorder then
			ShowMessage(concat('Reached a border after ', Long2str(Counter), ' steps'))
		else if MaxCountReached then
			ShowMessage(concat('Maximum number of steps (', Long2str(Counter), ') reached'))
		else
			ShowMessage('Completed');
		BackgroundIndex := SaveBackground;

{Close the Crest Pathway window}
{    Info := NewInfo;}
{    DoClose;}
{define a ROI along the pathway}
		Info := OldInfo;
		MakeOutline(FreeHandROI);
		UpdatePicWindow;

	end;			{CrestPathway}



{jvh}
	procedure XYQuantile (FilterW, FilterH, Quantile: integer);
		const
			MaxFilterDim = 11;
			SqMaxFilterDim = 121; {Square of MaxFilterDim}
		const
			PixelsPerUpdate = 5000;
			DefFilterW = 5;
			DefFilterH = 5;
		type
			LineArray = array[1..MaxFilterDim] of LineType;
			LineArrayPtr = ^LineArray;
		var
			pt: point;
			t: FateTable;
			LinesPerUpdate, min, max, maxj, minj, row, width: integer;
			Mark, NewMark, Xmargin, Ymargin, c, i, j, FilterSize: integer;
			offset, StartTicks: LongInt;
			N: packed array[1..SqMaxFilterDim] of UnsignedByte;
			MaskRect, frame, tRect: rect;
			AutoSelectAll, UseMask: boolean;
			LineBuff: LineArrayPtr;
			str: str255;
			p1, p2, lptr: ptr;
	begin
		StartTicks := TickCount;
		if info^.PictureType = NullPicture then
			exit(XYQuantile);
		lptr := NewPtr(SizeOf(LineArray));
		if lptr = nil then
			begin
				PutMemoryAlert;
				exit(XYQuantile);
			end;
		LineBuff := pointer(lptr);
		DisableDensitySlice;
		if (FilterW < 1) or (FilterW > MaxFilterDim) then
			begin
				FilterW := GetInt(concat('X Dimension of Filter (1-', Long2str(MaxFilterDim), ')'), DefFilterW, Canceled);
				if Canceled or (FilterW > MaxFilterDim) or (FilterW < 1) then
					exit(XYQuantile);
			end;
		if (FilterH < 1) or (FilterH > MaxFilterDim) then
			begin
				FilterH := GetInt(concat('Y Dimension of Filter (1-', Long2str(MaxFilterDim), ')'), DefFilterH, Canceled);
				if Canceled or (FilterH > MaxFilterDim) or (FilterH < 1) then
					exit(XYQuantile);
			end;
		FilterSize := FilterW * FilterH;
		if FilterSize < 2 then
			exit(XYQuantile);
		if (Quantile < 1) or (Quantile > FilterSize) then
			begin
				Quantile := GetInt(concat('Quantile (1-', Long2str(FilterSize), ')'), (FilterSize div 2) + 1, Canceled);
				if Canceled or (Quantile > FilterSize) or (Quantile < 1) then
					exit(XYQuantile);
			end;
		if NotinBounds then
			exit(XYQuantile);
{the Original Reduce Noise function is Faster then this one}
		if (FilterH = 3) and (FilterW = 3) and (Quantile = 5) then
			begin
				Filter(ReduceNoise, 0, t);
				exit(XYQuantile);
			end;
		AutoSelectAll := not Info^.RoiShowing;
		if AutoSelectAll then
			SelectAll(false);
		if info^.RoiType <> RectRoi then
			UseMask := SetupMask
		else
			UseMask := false;
		UpdatePicWindow;
		SetupUndoFromClip;
		WhatToUndo := UndoFilter;
		frame := info^.RoiRect;
		with frame, Info^ do
			begin
				changes := true;
				Xmargin := FilterW div 2;
				Ymargin := FilterH div 2;
				if left < Xmargin then
					left := left + Xmargin;
				if right > (PicRect.right - Xmargin) then
					right := right - Xmargin;
				if top < Ymargin then
					top := top + Ymargin;
				if bottom > (PicRect.bottom - Ymargin) then
					bottom := bottom - Ymargin;
				SetPort(wptr);
				PenNormal;
				PenPat(pat[PatIndex]);
				tRect := frame;
				OffscreenToScreenRect(tRect);
				FrameRect(tRect);
				width := right - left;
				Mark := RoiRect.top;
				LinesPerUpdate := (PixelsPerUpdate div width) div 3;
				str := concat(' Filter Size = (', Long2str(FilterW), ',', Long2str(FilterH), ')', cr, 'Quantile = ', Long2str(Quantile));
				ValuesMessage := Concat(str, cr, CmdPeriodToStop);
				ShowValues;
				ShowWatch;
				if Quantile >= (FilterSize div 2) + 1 then
					for row := top to bottom do
						begin
							for c := left to right do
								begin
									for i := 1 to FilterH do
										begin
											offset := Longint(row - Ymargin + i - 1) * BytesPerRow + c - Xmargin;
											p1 := ptr(ord4(PicBaseAddr) + offset);
											p2 := ptr(ord4(@N) + (i - 1) * FilterW);
											BlockMove(p1, p2, FilterW);
										end;
									for i := 1 to FilterSize - Quantile do
										begin
											max := 0;
											maxj := 1;
											for j := 1 to FilterSize do
												if N[j] > max then
													begin
														max := N[j];
														maxj := j;
													end;
											N[maxj] := 0;
										end;
									max := 0;
									for j := 1 to FilterSize do
										if N[j] > max then
											max := N[j];
									LineBuff^[Ymargin + 1, c - left] := max;
								end; {for c := left to right}
							if (row > (top + Ymargin)) then
								if UseMask then
									PutLineUsingMask(left, row - Ymargin, width, LineBuff^[1])
								else
									PutLine(left, row - Ymargin, width, LineBuff^[1]);
							for i := 1 to Ymargin do
								BlockMove(@LineBuff^[i + 1], @LineBuff^[i], width);
							if (row mod LinesPerUpdate) = 0 then
								begin
									pt.h := RoiRect.left;
									pt.v := row + 1;
									NewMark := pt.v;
									with RoiRect do
										SetRect(MaskRect, left, mark, right, NewMark);
									UpdateScreen(MaskRect);
									Mark := NewMark;
									if magnification > 1.0 then
										Mark := Mark - 1;
									if CommandPeriod then
										begin
											UpdatePicWindow;
											beep;
											if AutoSelectAll then
												KillRoi;
											exit(XYQuantile)
										end;
								end;
						end {for row:=...}
				else
					for row := top to bottom do
						begin
							for c := left to right do
								begin
									for i := 1 to FilterH do
										begin
											offset := Longint(row - Ymargin + i - 1) * BytesPerRow + c - Xmargin;
											p1 := ptr(ord4(PicBaseAddr) + offset);
											p2 := ptr(ord4(@N) + (i - 1) * FilterW);
											BlockMove(p1, p2, FilterW);
										end;
									for i := 1 to Quantile - 1 do
										begin
											min := 255;
											minj := 1;
											for j := 1 to FilterSize do
												if N[j] < min then
													begin
														min := N[j];
														minj := j;
													end;
											N[minj] := 255;
										end;
									min := 255;
									for j := 1 to FilterSize do
										if N[j] < min then
											min := N[j];
									LineBuff^[Ymargin + 1, c - left] := min;
								end; {for c := left to right}
							if (row > (top + Ymargin)) then
								if UseMask then
									PutLineUsingMask(left, row - Ymargin, width, LineBuff^[1])
								else
									PutLine(left, row - Ymargin, width, LineBuff^[1]);
							for i := 1 to Ymargin do
								BlockMove(@LineBuff^[i + 1], @LineBuff^[i], width);
							if (row mod LinesPerUpdate) = 0 then
								begin
									pt.h := RoiRect.left;
									pt.v := row + 1;
									NewMark := pt.v;
									with RoiRect do
										SetRect(MaskRect, left, mark, right, NewMark);
									UpdateScreen(MaskRect);
									Mark := NewMark;
									if magnification > 1.0 then
										Mark := Mark - 1;
									if CommandPeriod then
										begin
											UpdatePicWindow;
											beep;
											if AutoSelectAll then
												KillRoi;
											exit(XYQuantile)
										end;
								end;
							if CommandPeriod then
								begin
									UpdatePicWindow;
									beep;
									exit(XYQuantile);
								end;
						end; {for row:=...}

				for i := 2 to Ymargin + 1 do
					if UseMask then
						PutLineUsingMask(left, row - Ymargin + i - 1, width, LineBuff^[i])
					else
						PutLine(left, row - Ymargin + i - 1, width, LineBuff^[i]);

				ShowTime(StartTicks, frame, str);
			end; {with}
		UpdatePicWindow;
		SetupRoiRect;
		if AutoSelectAll then
			KillRoi;
		Dispose(LineBuff);
	end; {XYQuantile}



{jvh}
	procedure ZQuantile (Quantile: integer);
{ when 0 is given as parameter, a dialog box asks the}
{ user to enter the value of the quantile}
		const
			MaxWidth = 768;
			MaxDepth = 100;
		type
			ArrayBuff = array[1..MaxDepth] of LineType;
			SortingArrayPtr = ^ArrayBuff;
		var
			slices, sRow, aRow, s, i, j, k, SaveSlice, place: integer;
			min, minj, max, maxj, width, height, hstart, vStart: integer;
			iptr, p: Ptr;
			OldInfo, NewInfo: InfoPtr;
			aLine: LineType;
			mask, frame: rect;
			SortedLine: SortingArrayPtr;
			offset, StartTicks: Longint;
			AutoSelectAll: boolean;
			str: str255;
	begin
		with info^, info^.StackInfo^ do
			begin
				if CurrentSlice = 0 then
					begin
						PutMessage('This function is only available for stacks');
						Exit(ZQuantile);
					end;
				if nSlices < 2 then
					begin
						PutMessage('Z Quantile requires at least 2 slices.');
						macro := false;
						exit(ZQuantile);
					end;
			end;
		OldInfo := Info;
		with info^ do
			begin
				AutoSelectAll := not Info^.RoiShowing;
				if AutoSelectAll then
					SelectAll(true);
				with RoiRect do
					begin
						hStart := left;
						vStart := top;
						width := right - left;
						height := bottom - top;
					end;
				if width > MaxWidth then
					begin
						PutMessage(concat('Image can''t average selections wider than ', Long2str(MaxWidth), ' pixels.'));
						exit(ZQuantile);
					end;
				with StackInfo^ do
					begin
						slices := StackInfo^.nSlices;
						SaveSlice := CurrentSlice;
						Canceled := false;
						if (Quantile < 1) or (Quantile > slices) then
							begin
								Quantile := GetInt(concat('Quantile (1-', Long2str(slices), ')'), slices div 2 + 1, Canceled);
								if Canceled then
									exit(ZQuantile);
								if (Quantile > slices) or (Quantile < 1) then
									begin
										PutMessage(concat('Choose a number between 1 and ', Long2str(slices)));
										exit(ZQuantile);
									end;
							end;
					end;
			end;
		str := concat('Stack depth =', Long2str(slices), cr, 'Quantile = ', Long2str(Quantile));
		ValuesMessage := Concat(str, cr, CmdPeriodToStop);
		ShowValues;
		iptr := NewPtr(SizeOf(ArrayBuff));
		if iptr = nil then
			begin
				PutMemoryAlert;
				exit(ZQuantile);
			end;
		SortedLine := pointer(iptr);
		if not NewPicWindow(concat('Quantile ', Long2str(Quantile), ' of ', Long2str(slices)), width, height) then
			begin
				exit(ZQuantile);
				Dispose(SortedLine);
			end;
		NewInfo := Info;
		ShowWatch;
		StartTicks := TickCount;
		aRow := 0;

		if Quantile >= (slices div 2) + 1 then
			for sRow := vStart to vStart + height - 1 do
				begin
					info := OldInfo;
					for s := 1 to slices do
						begin
							SelectSlice(s);
							GetLine(hStart, sRow, width, SortedLine^[s]);
						end;
					for k := 0 to width - 1 do
						begin
							for i := 1 to (slices - Quantile) do
								begin
									max := 0;
									maxj := 1;
									for j := 1 to slices do
										if max < SortedLine^[j, k] then
											begin
												max := SortedLine^[j, k];
												maxj := j;
											end;
									SortedLine^[maxj, k] := 0;
								end;
							max := 0;
							for j := 1 to slices do
								if max < SortedLine^[j, k] then
									max := SortedLine^[j, k];
							aLine[k] := max;
						end;
					info := NewInfo;
					PutLine(0, aRow, width, aLine);
					SetRect(mask, 0, aRow, width, aRow + 1);
					aRow := aRow + 1;
					UpdateScreen(mask);
					if CommandPeriod then
						leave;
				end
		else
			for sRow := vStart to vStart + height - 1 do
				begin
					info := OldInfo;
					for s := 1 to slices do
						begin
							SelectSlice(s);
							GetLine(hStart, sRow, width, SortedLine^[s]);
						end;
					for k := 0 to width - 1 do
						begin
							for i := 1 to (Quantile - 1) do
								begin
									min := 255;
									minj := 1;
									for j := 1 to slices do
										if min > SortedLine^[j, k] then
											begin
												min := SortedLine^[j, k];
												minj := j;
											end;
									SortedLine^[minj, k] := 255;
								end;
							min := 255;
							for j := 1 to slices do
								if min > SortedLine^[j, k] then
									min := SortedLine^[j, k];
							aLine[k] := min;
						end;
					info := NewInfo;
					PutLine(0, aRow, width, aLine);
					SetRect(mask, 0, aRow, width, aRow + 1);
					aRow := aRow + 1;
					UpdateScreen(mask);
					if CommandPeriod then
						leave;
				end;
		NewInfo^.Changes := true;
		Dispose(SortedLine);
		info := OldInfo;
		frame := info^.RoiRect;
		SelectSlice(SaveSlice);
		ShowTime(StartTicks, frame, str);
		if AutoSelectAll then
			KillRoi;
	end; {ZQuantile}


	procedure DoUserMenuEvent (MenuItem: integer);
		var
			ResetSteps: integer;
	begin
		case MenuItem of
			XYQuantileItem: 
				begin
					XYQuantile(0, 0, 0);
				end;
			ZQuantileItem: 
				ZQuantile(0);
			CrestPathwayItem: 
				CrestPathway(False);
			AssistedTrackingItem: 
				CrestPathway(True);
			CrestOptionsItem: 
				CrestOptions;
			otherwise
				;
		end;
	end;

	procedure OldUserMacroCode (CodeNumber: integer; Param1, Param2, Param3: extended);
	begin
		case CodeNumber of
			1: 
				XYQuantile(round(Param1), round(Param2), round(Param3));
			2: 
				ZQuantile(round(Param1));
			3: 
				CrestPathway(false);
			4: 
				CrestPathway(true);
			5: 
				CrestOptions;
			otherwise
				;
		end;
	end;

	procedure ShowNoCodeMessage;
	begin
		PutMessage('Requires user written Think Pascal routine. ');
	end;


	procedure UserMacroCode (str: str255; Param1, Param2, Param3: extended);
	begin
		MakeLowerCase(str);
		if pos('assistedtracking', str) <> 0 then
			begin
				CrestPathway(true);
				exit(UserMacroCode);
			end;
		if pos('crestpathway', str) <> 0 then
			begin
				CrestPathway(false);
				exit(UserMacroCode);
			end;
		ShowNoCodeMessage;
	end;


end.
{}
{}
{}
procedure CrestOptions;
begin
	AveragingDist := GetInt(concat('Averaging Distance (1-', Long2str(MaxAveragingDist), ')'), AveragingDist, Canceled);
	if Canceled or (AveragingDist > MaxAveragingDist) or (AveragingDist < 1) then
		exit(CrestOptions);

	AssistanceRadius := GetInt(concat('Assistance Radius (1-', Long2str(MaxAssistanceRadius), ')'), AssistanceRadius, Canceled);
	if Canceled or (AssistanceRadius > 3 * MaxAssistanceRadius) or (AssistanceRadius < 1) then
		exit(CrestOptions);

	MaxAngle := GetInt(concat('MaxAngle (1-', Long2str(MaxMaxAngle), ')'), MaxAngle, Canceled);
	if Canceled or (MaxAngle > 3) or (MaxAngle < 1) then
		exit(CrestOptions);

	InitialDir := GetInt(concat('Initial Direction (0-7)'), InitialDir, Canceled);
	if Canceled or (InitialDir > 7) or (InitialDir < 0) then
		exit(CrestOptions);

	MaxCount := GetInt(concat('Maximal Path Lenght (0-', Long2str(MaxMaxCount), ')'), MaxCount, Canceled);
	if Canceled or (MaxCount > MaxMaxCount) or (MaxCount < 0) then
		exit(CrestOptions);
end;
{}