unit User;

{
This module is the best place to put user additions to NIH Image. Uncomment
the call to InitUser in Image.p to activate the User menu. Edit the "User"
menu resource with ResEdit to customize the User menu.
}


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, Filters, Analysis;


	procedure InitUser;
	procedure DoUserCommand1;
	procedure DoUserCommand2;
	procedure DoUserMenuEvent (MenuItem: integer);
	procedure OldUserMacroCode (CodeNumber: integer; Param1, Param2, Param3: extended);
	procedure UserMacroCode (str: str255; Param1, Param2, Param3: extended);


implementation

{User global variables go here.}
	var
		color, MinSpacing: integer;
		SaveInfo: InfoPtr;
		PeakRadius, Peakedness: extended;


	procedure InitUser;
	begin
		UserMenuH := GetMenu(UserMenu);
		InsertMenu(UserMenuH, 0);
		DrawMenuBar;
{Additional user initialization code goes here.}
	end;


	procedure DrawDot (row, column, RowOffset, ColumnOffset: integer; big: boolean);
		var
			h, v: integer;
	begin
		if big then begin
				for h := -1 to 1 do
					for v := -1 to 1 do
						PutPixel(column * 16 + ColumnOffset * 4 + h + 16, row * 16 + RowOffset * 4 + v + 16, color)
			end
		else
			PutPixel(column * 16 + ColumnOffset * 4 + 16, row * 16 + RowOffset * 4 + 16, color);
	end;

	procedure DrawNeighborhood (i, row, column: integer);

	begin
		DrawDot(row, column, 0, 0, BitAnd(i, 1) = 1);
		DrawDot(row, column, 0, 1, BitAnd(i, 2) = 2);
		DrawDot(row, column, 0, 2, BitAnd(i, 4) = 4);
		DrawDot(row, column, 1, 2, BitAnd(i, 8) = 8);
		DrawDot(row, column, 2, 2, BitAnd(i, 16) = 16);
		DrawDot(row, column, 2, 1, BitAnd(i, 32) = 32);
		DrawDot(row, column, 2, 0, BitAnd(i, 64) = 64);
		DrawDot(row, column, 1, 0, BitAnd(i, 128) = 128);
		DrawDot(row, column, 1, 1, true);
	end;


	procedure SetColor (i: integer);
{Color neighborhoods to show which ones would be removed on the first pass(150), second pass(100),}
{or either pass(200) when using the Zhang and Suen thinning algorithm(CACM, Mar. 1984,236-239).}
		var
			p2, p3, p4, p5, p6, p7, p8, p9, A, B: integer;
	begin
		p2 := bsr(band(i, 2), 1);
		p3 := bsr(band(i, 4), 2);
		p4 := bsr(band(i, 8), 3);
		p5 := bsr(band(i, 16), 4);
		p6 := bsr(band(i, 32), 5);
		p7 := bsr(band(i, 64), 6);
		p8 := bsr(band(i, 128), 7);
		p9 := band(i, 1);
		A := 0;
		if (p2 = 0) and (p3 = 1) then
			A := A + 1;
		if (p3 = 0) and (p4 = 1) then
			A := A + 1;
		if (p4 = 0) and (p5 = 1) then
			A := A + 1;
		if (p5 = 0) and (p6 = 1) then
			A := A + 1;
		if (p6 = 0) and (p7 = 1) then
			A := A + 1;
		if (p7 = 0) and (p8 = 1) then
			A := A + 1;
		if (p8 = 0) and (p9 = 1) then
			A := A + 1;
		if (p9 = 0) and (p2 = 1) then
			A := A + 1;
		B := p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;
		color := 255;
		if A = 1 then
			if (B >= 2) and (B <= 6) then begin
					if ((p2 * p4 * p6 = 0) and (p4 * p6 * p8 = 0)) and ((p2 * p4 * p8 = 0) and (p2 * p6 * p8 = 0)) then
						color := 200
					else if (p2 * p4 * p6 = 0) and (p4 * p6 * p8 = 0) then
						color := 150
					else if (p2 * p4 * p8 = 0) and (p2 * p6 * p8 = 0) then
						color := 100;
				end;
	end;


	procedure DoUserCommand1;
{Generates a table showing all possible 3x3 neighborhoods. This table is used}
{ for making up the "fate table" used by the Skeletonize command and the Wand tool.}
		var
			row, column, index: integer;
	begin
		row := 0;
		column := 0;
		if NewPicWindow('Fate Table', 600, 200) then
			for index := 0 to 255 do begin
					SetColor(index);
					DrawNeighborhood(index, row, column);
					column := column + 1;
					if column = 32 then begin
							row := row + 1;
							column := 0;
						end;
				end;
	end;


	function isPeak (x, y, minValue: LongInt): boolean;
		var
			delta, angle, dx, dy: extended;
			v, i, v2, maxv2, x2, y2, v2count, nSamples: integer;
			sample: LineType;
			minlower, count, nLower, maxCount: integer;
			PeakFound: boolean;
			mask: rect;
	begin
		isPeak := false;
		v := MyGetPixel(x, y);
		if v < minValue then
			exit(isPeak);
		if v <= MyGetPixel(x + 1, y) then
			exit(isPeak);
		if v <= MyGetPixel(x + 1, y + 1) then
			exit(isPeak);
		if v <= MyGetPixel(x, y + 1) then
			exit(isPeak);
		if v <= MyGetPixel(x - 1, y + 1) then
			exit(isPeak);
		if v < MyGetPixel(x - 1, y) then
			exit(isPeak);
		if (v < MyGetPixel(x - 1, y - 1)) then
			exit(isPeak);
		if v < MyGetPixel(x, y - 1) then
			exit(isPeak);
		if v < MyGetPixel(x + 1, y - 1) then
			exit(isPeak);
		nSamples := round(4 * PeakRadius);
		delta := 2.0 * pi / nsamples;
		angle := 0.0;
		maxv2 := round((1.0 - Peakedness) * v);
		for i := 1 to nSamples do begin
				dx := PeakRadius * cos(angle);
				dy := PeakRadius * sin(angle);
				sample[i] := round(GetInterpolatedPixel(x + dx, y + dy));
				angle := angle + delta;
			end;
		minLower := round(0.677 * nsamples);
		PeakFound := false;
		count := 0;
		i := 1;
		nLower := 0;
		maxCount := nSamples + minLower;
		repeat
			if sample[i] <= maxv2 then
				nLower := nLower + 1
			else
				nLower := 0;
			PeakFound := nLower >= minLower;
			i := i + 1;
			if i > nSamples then
				i := 1;
			count := count + 1;
		until PeakFound or (count = maxCount);
		if PeakFound then begin
				info := SaveInfo;
				with info^ do begin
						SetRect(RoiRect, x - MinSpacing + 1, y - MinSpacing + 1, x + MinSpacing, y + MinSpacing);
						with RoiRect do begin
								if left < 0 then
									left := 0;
								if top < 0 then
									top := 0;
								if right > PicRect.right then
									right := PicRect.right;
								if bottom > PicRect.bottom then
									bottom := PicRect.bottom;
							end;
						GetRectHistogram;
						PeakFound := histogram[0] = 0;
					end; {with}
				Info := UndoInfo;
			end;
		isPeak := PeakFound;
	end;


	procedure FindPeaks (minValue, PeakRadiusP, PeakednessP: extended);
		var
			x, y, i, iMinValue: integer;
			AutoSelectAll: boolean;
			srect, mask: rect;
			count: LongInt;
			t: FateTable;
	begin
		if NotRectangular or NotInBounds or NoUndo then
			exit(FindPeaks);
		iMinValue := round(minValue);
		if iMinValue < 0 then
			iMinValue := 0;
		if iMinValue > 255 then
			iMinValue := 255;
		PeakRadius := PeakRadiusP;
		if PeakRadius = 0.0 then
			PeakRadius := 6.0;
		if PeakRadius < 1.0 then
			PeakRadius := 1.0;
		if PeakRadius > 50.0 then
			PeakRadius := 50.0;
		MinSpacing := round(PeakRadius) - 1;
		if MinSpacing < 1 then
			MinSpacing := 1;
		if MinSpacing > 4 then
			MinSpacing := 4;
		Peakedness := PeakednessP;
		if Peakedness = 0.0 then
			Peakedness := 0.2;
		if Peakedness < 0.05 then
			Peakedness := 0.05;
		if Peakedness > 0.95 then
			Peakedness := 0.95;
		AutoSelectAll := not Info^.RoiShowing;
		if AutoSelectAll then
			SelectAll(true);
		ShowWatch;
		SetupUndo;
		WhatToUndo := UndoEdit;
		SetupUndoInfoRec;
		SaveInfo := Info;
		srect := info^.roiRect;
		KillRoi;
		ChangeValues(0, 0, 1);
		info := UndoInfo;
		count := 0;
		with srect do
			for y := top to bottom - 1 do begin
					if CommandPeriod then begin
							beep;
							Info := SaveInfo;
							leave;
						end;
					for x := left to right - 1 do
						if isPeak(x, y, iMinValue) then begin
								count := count + 1;
								Info := SaveInfo;
								PutPixel(x, y, 0);
{PutPixel(x - 1, y, 0);}
{PutPixel(x - 1, y - 1, 0);}
{PutPixel(x, y - 1, 0);}
								SetRect(mask, x - 1, y - 1, x + 1, y + 1);
								UpdateScreen(mask);
								Info := UndoInfo;
								if count < MaxMeasurements then begin
										User1^[count] := x;
										User2^[count] := y;
									end;
								if (y mod 50) = 0 then ShowMessage(concat(long2str(y), '  ', long2str(count)));
							end;
				end;
		Info := SaveInfo;
		if count < MaxMeasurements then begin
				UnsavedResults := false;
				ResetCounter;
				for i := 1 to count do begin
						ClearResults(i);
						xcenter^[i] := User1^[i];
						ycenter^[i] := User2^[i];
					end;
				mCount := count;
				UpdateList;
				ShowInfo;
			end
		else
			PutError('"Max Measurements" is too small.');
		ShowMessage(concat('Count=', long2str(count), crStr, 'Threshold=', long2str(iMinValue)));
	end;



	procedure ComputeBirefringence (scale, offset: extended);
{This an example of how to do image math using a UserCode macro routine.}
{It executes the following formula}

      {SQRT ( ( I1 - I2 ) ^ 2 + ( I3 - I4 ) ^ 2 ) / ( I1 + I2 - I3 + I4 ) ,}

{where I1 , I2 , I3 , I4  are the first four slices of the current stack.}
{The result in the fifth slice of the stack.}

		var
			i1, i2, i3, i4, i5: LineType;
			i, slice, row: integer;
			mask: rect;
			v, min, max: extended;
			minstr, maxstr: str255;
	begin
		with info^ do begin
				if StackInfo = nil then
					exit(ComputeBirefringence);
				if StackInfo^.nSlices <> 5 then
					exit(ComputeBirefringence);
				min := 1.0e12;
				max := -1.0e12;
				for row := 0 to nLines - 1 do begin
						SelectSlice(1);
						GetLine(0, row, PixelsPerLine, i1);
						SelectSlice(2);
						GetLine(0, row, PixelsPerLine, i2);
						SelectSlice(3);
						GetLine(0, row, PixelsPerLine, i3);
						SelectSlice(4);
						GetLine(0, row, PixelsPerLine, i4);
						for i := 0 to PixelsPerLine - 1 do begin
								v := sqrt(sqr(I1[i] - I2[i]) + sqr(I3[i] - I4[i])) / (I1[i] + I2[i] - I3[i] + I4[i]);
								if v < min then
									min := v;
								if v > max then
									max := v;
								if v > 255 then
									v := 255;
								if v < 0 then
									v := 0;
								v := v * scale + offset;
								i5[i] := round(v);
							end;
						SelectSlice(5);
						PutLine(0, row, PixelsPerLine, i5);
						SetRect(mask, 0, row, PixelsPerLine, row + 1);
						UpdateScreen(mask);
						if CommandPeriod then
							leave;
					end;
			end;
		RealToString(min, 1, 4, minstr);
		RealToString(max, 1, 4, maxstr);
		ShowMessage(concat('min=', minstr, crStr, 'max=', maxstr));
	end;


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


	procedure DoUserCommand2;
	begin
		ShowNoCodeMessage
	end;


	procedure DoUserMenuEvent (MenuItem: integer);
	begin
		case MenuItem of
			1: 
				DoUserCommand1;
			2: 
				DoUserCommand2;
		end;
	end;


	procedure OldUserMacroCode (CodeNumber: integer; Param1, Param2, Param3: extended);
  {Obsolete version kept for backward compatibilty.}
	begin
		case CodeNumber of
			1: 
				ShowNoCodeMessage;
			2: 
				ShowNoCodeMessage;
			3: 
				ShowNoCodeMessage;
			4: 
				ShowNoCodeMessage;
			5: 
				FindPeaks(param1, param2, param3);
			otherwise
				ShowNoCodeMessage;
		end;
	end;


	procedure UserMacroCode (str: str255; Param1, Param2, Param3: extended);
	begin
		MakeLowerCase(str);
		if pos('peaks', str) <> 0 then begin
				FindPeaks(param1, param2, param3);
				exit(UserMacroCode);
			end;
		if pos('birefringence', str) <> 0 then begin
				ComputeBirefringence(param1, param2);
				exit(UserMacroCode);
			end;
		ShowNoCodeMessage;
	end;


end.