unit UMQuantile;

{Contributed by Jacques van Helden <jvanheld@dbm.ulb.ac.be>}

{User Macro "Quantile" package. }
{Macros which use these extensions should specify "requiresUser('Quantile',1)".}

{Code moved from User.p to UMQuantile.p by Edward J. Huff huff@mcclb0.med.nyu.edu}
{The original was "User.p (Quantile+Crest)" last modified Mon, Aug 23, 1993, 11:20 AM.}

{These instructions apply if you have received only UMQuantile.p and}
{a copy of UMacroDef.p and UMacroRun.p}

{You need the Think Pascal compiler to use this source code. }

{  Special for UMQuantile.p:}
{- Place the PutLineUsingMask procedure in the interfaces of the Functions  }
{    unit (Functions.p file)   }

{Installation instructions.  At the moment, the UMacroDef.p and UMacroRun.p}
{have not been accepted for distribution with standard Image.  If they are, then}
{these instructions will apply.  Otherwise, you will have to obtain a copy of image}
{with UMacroDef.p etc. already installed.  This will typically be available in}
{the /pub/nih-image/contrib directory of zippy.nimh.nih.gov and the file}
{name will include "UMX".}

{If you have no other user macro extension packages installed,}
{then just replace the default UMacroDef.p and UMacroRun.p with the enclosed}
{versions.  Otherwise, you must merge the changes, which should be found only}
{at well marked places.  Generally, this means simply adding lines, although you }
{may also need to add a comma here and there.  Changes should be found at}
{one place in UMacroDef.p and seven places in UMacroRun.p.}

{You will also have to use the "Project/Add File..." menu item to tell Think Pascal }
{to compile this file, and use the "Windows/Image.proj" item to get the project window, }
{and drag this file from the bottom of the list to somewhere between UMacroDef.p and }
{UMacroRun.p.}

{Finally, recompile and rebuild:  Use "Run/Build" (command B) and }
{"Project/Build Application...".}

{  In case of problems when installing the XYQuantile and ZQuantile routines in}
{  your version of  NIH Image, contact Jacques van Helden(jvanheld@dbm.ulb.ac.be)  }
{  or Edward J. Huff <huff@mcclb0.med.nyu.edu> for problems with UMacroDef.p etc.}

{  This file contains two filters based on the principle of median            }
{  filtering. Median filtering consists in replacing each pixel by the        }
{  medin value of its neighbours. It allows an edge-preserving noise          }
{  reduction.  The Reduce Noise function of NIH-Image consists in a median    }
{  filtering on a 3x3 square.   }
{}
{  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).   }
{     }
{  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.   }
{}
{    If you find some use of or have some comment about these filters,   }
{  please contact jvanheld@dbm.ulb.ac.be   }
{}


interface

	uses
		QuickDraw, Palettes, PrintTraps, Globals, Utilities, Graphics, Lut, Analysis, Filters, UMacroDef;

	procedure UMQuantileInit;
	procedure UMQuantileFinal;
	procedure UMQuantileAdd;
	procedure UMQuantileLookup (var uma: UserMacroArgs);
	procedure UMQuantileRun (var uma: UserMacroArgs);

implementation

	var
		Canceled: boolean;

{Called from procedure InitUserMacros in UMacroRun.p, }
{which is called from Image.p early in initialization.}
	procedure UMQuantileInit;
	begin
	end;


{Called from procedure FinalUserMacros in UMacroRun,p.}
{This will eventually be guaranteed to run prior to any exit, and}
{is intended for things which MUST be done prior to exit, like removing}
{timers from the system timer list.}
	procedure UMQuantileFinal;
	begin
	end;


{AddUMSym calls:}
{Add one call for each macro command, function, or string function}
{you wish to add to the macro language.}
{First argument is a string, case is ignored, truncated to SymbolSize characters.}
{Second argument must be one of UserCommandT, UserFuncT, or UserStrFuncT}
{Third argument is the UserCommandType item associated with the name.}


{Called from procedure AddUserMacros in UMacroRun.p.}
{This runs once each time macros are loaded from a file or a text window.}
	procedure UMQuantileAdd;
	begin
		AddUMSym('XYQuantile', UserCommandT, XYQuantileUC);
		AddUMSym('ZQuantile', UserCommandT, ZQuantileUC);
	end;

{Called from procedure LookupUserMacro in UMMacroRun.p}
{This runs every time the macro is executed, just prior to}
{parsing the arguments.}
	procedure UMQuantileLookup (var uma: UserMacroArgs);
	begin
		with uma do
			case UserMacroCommand of
				XYQuantileUC:  begin
						nArgs := 3;
						arg[1].atype := UMATinteger; {width}
						arg[2].atype := UMATinteger; {height}
						arg[3].atype := UMATinteger; {quantile}
					end;
				ZQuantileUC:  begin
						nArgs := 1;
						arg[1].atype := UMATinteger;{quantile}
					end;
				otherwise begin
						ErrorOccurred := true;
						str := 'UMQuantile.p LookupUserMacro';
					end;
			end;
	end;

{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));
				InfoMessage := Concat(str, cr, CmdPeriodToStop);
				ShowInfo;
				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
{- Place the PutLineUsingMask procedure in the interfaces of the Functions  }
{    unit (Functions.p file)   }
									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));
		InfoMessage := Concat(str, cr, CmdPeriodToStop);
		ShowInfo;
		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}


{Called from procedure DoUserMacro in UMacroRun.p .}
{Do not change uma.nArgs or uma.arg[i].argt here.}
{This runs once each time the macro is used, after parsing the arguments.}
	procedure UMQuantileRun (var uma: UserMacroArgs);
	begin
		with uma do
			case UserMacroCommand of
				XYQuantileUC: 
					XYQuantile(arg[1].ival, arg[2].ival, arg[3].ival);
				ZQuantileUC: 
					ZQuantile(arg[1].ival);
				otherwise begin
						ErrorOccurred := true;
						str := 'UMQuantile.p DoUserMacro';
					end;
			end;
	end;


end.