unit User;

{This module is a good place to put user additions to NIH Image. You will need }
{to uncomment the call to InitUser in Image.p.}


interface

	uses
		Memory, QuickDraw, Packages, Menus, Events, Fonts, Scrap, ToolUtils,Resources, Errors, Palettes, 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 < 10 then
			iMinValue := 10;
		if iMinValue > 150 then
			iMinValue := 150;
		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 Think 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.