unit fft;

{
This file contains a Pascal language implementation of the 
Fast Hartley Transform algorithm which is covered under
United States Patent Number 4,646,256.

This code may therefore be freely used and distributed only
under the following conditions:

		1)  This header is included in every copy of the code; and
		2)  The code is used for noncommercial research purposes only.

 Firms using this code for commercial purposes will be infringing a United
 States patent and should contact the

					Office of Technology Licensing
					Stanford University
					857 Serra Street, 2nd Floor
					Stanford, CA   94305-6225
					(415) 723 0651

This implementation is based on Pascal code contibuted
by Arlo Reeves (arlo@apple.com).
}

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 doFFT(fftKind: fftType);
	procedure DisplayPowerSpectrum(x: rImagePtr);
	function   isPowerOf2 (x: LongInt): boolean;
	procedure SetFFTWindowName(doInverse: boolean);
	procedure RedisplayPowerSpectrum;
	procedure doSwapQuadrants;
	function isFFT: boolean;
	procedure ShowFFTValues (hloc, vloc, ivalue: LongInt);


implementation

	const
		UpdateTicks = 10; {Show progress 6 times/sec}
		
	var
		AbortFFT: boolean;
		C, S: rLinePtr;


	function log2 (x: LongInt): LongInt;
		var
			count: LongInt;
	begin
		count := 15;
		while not BTST(x, count) do
			count := count - 1;
		log2 := count;
	end;
	

	function BitRevX (x, bitlen: LongInt): LongInt;
		var
			i: LongInt;
			temp: longint;
	begin
		temp := 0;
		for i := 0 to bitlen do
			if BTST(x, i) then
				BSET(temp, bitlen - i - 1);
		BitRevX := LoWord(temp);
	end;
	

	procedure BitRevRArr (x: rLinePtr; bitlen, maxN: LongInt);
		var
			i: LongInt;
			tempArr: rLineType;
	begin
		for i := 0 to maxN - 1 do
			tempArr[i] := x^[BitRevX(i, bitlen)];
		BlockMove(@tempArr, x, maxN * SizeOf(real));
	end;
	

	procedure transposeR (x: rImagePtr; maxN: LongInt);
		var
			r, c: LongInt;
			rTemp: real;
	begin
		for r := 0 to maxN - 1 do
			for c := r to maxN - 1 do
				if r <> c then
					begin
						rTemp := x^[ r * maxN + c];
						x^[ r * maxN + c] := x^[c * maxN + r];
						x^[c * maxN + r] := rTemp;
					end;
	end;
	

	procedure FHTps (r, maxN: LongInt; inMat: rImagePtr;  var outArr: rLineType);
{ Power Spectrum of one row from 2D Hartley Transform }
		var
			c, base: LongInt;
	begin
		base := r * maxN;
		for c := 0 to maxN - 1 do
			outArr[c] := (sqr(inMat^[base + c]) + sqr(inMat^[((maxN - r) mod maxN) * maxN + (maxN - c) mod maxN])) / 2;
	end;


	procedure dfht3 (x: rLinePtr; inverse: boolean; maxN: LongInt);
	{ an optimized real FHT }
		var
			i, stage, gpNum, gpIndex, gpSize, numGps, Nlog2: LongInt;
			bfNum, numBfs: LongInt;
			Ad0, Ad1, Ad2, Ad3, Ad4, CSAd: LongInt;
			rt1, rt2, rt3, rt4: real;
	begin
		Nlog2 := log2(maxN);
		BitRevRArr(x, Nlog2, maxN);			{ bitReverse the input array }

		gpSize := 2;						{ first & second stages - do radix 4 butterflies once thru }
		numGps := maxN div 4;
		for gpNum := 0 to numGps - 1 do
			begin
				Ad1 := gpNum * 4;
				Ad2 := Ad1 + 1;
				Ad3 := Ad1 + gpSize;
				Ad4 := Ad2 + gpSize;
				rt1 := x^[Ad1] + x^[Ad2];			{ a + b }
				rt2 := x^[Ad1] - x^[Ad2];			{ a - b }
				rt3 := x^[Ad3] + x^[Ad4];			{ c + d }
				rt4 := x^[Ad3] - x^[Ad4];			{ c - d }
				x^[Ad1] := rt1 + rt3;				{ a + b + (c + d ) }
				x^[Ad2] := rt2 + rt4;				{ a - b + (c - d ) }
				x^[Ad3] := rt1 - rt3;				{ a + b - (c + d ) }
				x^[Ad4] := rt2 - rt4;				{ a - b - (c - d ) }
			end;

		if Nlog2 > 2 then
			begin						{ third + stages computed here }
				gpSize := 4;
				numBfs := 2;
				numGps := numGps div 2;
				for stage := 2 to Nlog2 - 1 do
					begin
						for gpNum := 0 to numGps - 1 do
							begin
								Ad0 := gpNum * gpSize * 2;
								Ad1 := Ad0;							{ 1st butterfly is different from others - no mults needed }
								Ad2 := Ad1 + gpSize;
								Ad3 := Ad1 + gpSize div 2;
								Ad4 := Ad3 + gpSize;
								rt1 := x^[Ad1];
								x^[Ad1] := x^[Ad1] + x^[Ad2];
								x^[Ad2] := rt1 - x^[Ad2];
								rt1 := x^[Ad3];
								x^[Ad3] := x^[Ad3] + x^[Ad4];
								x^[Ad4] := rt1 - x^[Ad4];
								for bfNum := 1 to numBfs - 1 do
									begin		{ subsequent BF's dealt with together }
										Ad1 := bfNum + Ad0;
										Ad2 := Ad1 + gpSize;
										Ad3 := gpSize - bfNum + Ad0;
										Ad4 := Ad3 + gpSize;

										CSAd := bfNum * numGps;
										rt1 := x^[Ad2] * C^[CSAd] + x^[Ad4] * S^[CSAd];
										rt2 := x^[Ad4] * C^[CSAd] - x^[Ad2] * S^[CSAd];

										x^[Ad2] := x^[Ad1] - rt1;
										x^[Ad1] := x^[Ad1] + rt1;
										x^[Ad4] := x^[Ad3] + rt2;
										x^[Ad3] := x^[Ad3] - rt2;

									end; 		{ for bfNum := 0 to ... }
							end;			{ for gpNum := 0 to ... }
						gpSize := gpSize * 2;
						numBfs := numBfs * 2;
						numGps := numGps div 2;
					end;
			end;		{ if }

		if inverse then
			for i := 0 to maxN - 1 do
				x^[i] := x^[i] / maxN;
	end;


	procedure rc2Dfht (x: rImagePtr; inverse: boolean; maxN: LongInt);
{ Row-column 2D FHT }
		var
			row, col, mRow, mCol, NextTicks: LongInt;
			A, B, C, D, E: real;
			theRow: rLinePtr;
	begin
		NextTicks := TickCount + UpdateTicks;
		for row := 0 to maxN - 1 do begin
			theRow := rLinePtr(ord4(x) + row * maxN * SizeOf(real));
			dfht3(theRow, inverse, maxN);
			if TickCount > nextTicks then begin
				UpdateMeter(round(row/maxN * 50.0), 'Computing FHT');
				nextTicks := TickCount + updateTicks;
				if CommandPeriod then begin
						beep;
						AbortFFT := true;
						exit(rc2Dfht)
					end;
			end;
		end;
		transposeR(x, maxN);
		for row := 0 to maxN - 1 do begin
			theRow := rLinePtr(ord4(x) + row * maxN * SizeOf(real));
			dfht3(theRow, inverse, maxN);
			if TickCount > nextTicks then begin
				UpdateMeter(round(row/maxN * 50.0) + 50, 'Computing FHT');
				nextTicks := TickCount + updateTicks;
				if CommandPeriod then begin
						beep;
						AbortFFT := true;
						exit(rc2Dfht)
					end;
			end;
		end;
		transposeR(x, maxN);
		for row := 0 to maxN div 2 do 			{ Now calculate actual Hartley transform }
			for col := 0 to maxN div 2 do
				begin
					mRow := (maxN - row) mod maxN;
					mCol := (maxN - col) mod maxN;
					A := x^[row * maxN + col];						{ see Bracewell, 'Fast 2D Hartley Transf.' IEEE Procs. 9/86 }
					B := x^[mRow * maxN + col];
					C := x^[row * maxN + mCol];
					D := x^[mRow * maxN + mCol];
					E := ((A + D) - (B + C)) / 2;
					x^[row * maxN + col] := A - E;
					x^[mRow * maxN + col] := B + E;
					x^[row * maxN + mCol] := C + E;
					x^[mRow * maxN + mCol] := D - E;
				end;
		UpdateMeter(-1, 'Hide Meter');
	end;
	
	
procedure SwapQuadrants;
{Swap quadrants 1 and 3 and quadrants 2 and 4 so the}
{power spectrum origin is at the center of the image.}
{2 1}
{3 4}
var
	row, col, halfMaxN, tmp, maxN: LongInt;
	line1, line2: LineType;
begin
	maxN := info^.PixelsPerLine;
	halfMaxN := maxN div 2;
	for row:= 0 to halfMaxN -1 do begin
		GetLine(0, row, maxN, line1);
		GetLine(0, row + halfMaxN, maxN, line2);
		for col := 0 to halfMaxN -1 do begin
			tmp := line1[col];
			line1[col] := line1[col + halfMaxN];
			line1[col + halfMaxN] := tmp;
			tmp := line2[col];
			line2[col] := line2[col + halfMaxN];
			line2[col + halfMaxN] := tmp;
		end;
		PutLine(0, row, maxN, line2);
		PutLine(0, row + halfMaxN, maxN, line1);
	end;
end;


	
	procedure DisplayRealImage(x: rImagePtr);
		var
			row, col, i, base, maxN: LongInt;
			r, min, max, scale: real;
			line: lineType;
	begin
		maxN := info^.PixelsPerLine;
		min := 1e99;
		max := -1e99;
		i := 1;
		for row := 0 to maxN - 1 do begin
			base := row * maxN;
			for col := 0 to maxN - 1 do begin
					r := x^[base + col];
					if r < min then min := r;
					if r > max then max := r;
			end;
		end;
		scale := 253.0 / (max - min);
		for row := 0 to maxN - 1 do begin
				base := row * maxN;
				for col := 0 to maxN - 1 do begin
						r := x^[base + col];
						line[col] := round((r - min) * scale) + 1;
				end;
				PutLine(0, row, maxN, line);
			end;
	end;


	procedure DisplayPSUsingBuffer(fht, ps: rImagePtr; var rLine: rLineType);
	var
		row, col, base, nextTicks, maxN: LongInt;
		r, min, max, scale: real;
		line: lineType;
		
	begin
		maxN := info^.PixelsPerLine;
		nextTicks := TickCount + updateTicks;
		min := 1e99;
		max := -1e99;
		for row := 0 to maxN - 1 do begin
			FHTps (row, maxN, fht, rLine);
			base := row * maxN;
			for col := 0 to maxN - 1 do begin
					r := rLine[col];
					if r < min then min := r;
					if r > max then max := r;
					ps^[base + col] := r;
			end;
			if TickCount > nextTicks then begin
				UpdateMeter(round(row/maxN * 80.0), 'Computing Power Spectrum');
				nextTicks := TickCount + updateTicks;
				if CommandPeriod then begin
						beep;
						AbortFFT := true;
						exit(DisplayPSUsingBuffer)
					end;
			end;
		end;
		if min < 1.0 then
			min := 0.0
		else
			min := ln(min);
		max := ln(max);
		scale := 253.0 / (max - min);
		for row := 0 to maxN - 1 do begin
				base := row * maxN;
				for col := 0 to maxN - 1 do begin
						r := ps^[base + col];
						if r < 1.0 then
							r := 0.0
						else
							r := ln(r);
						line[col] := round((r - min) * scale) + 1;
				end;
				PutLine(0, row, maxN, line);
				if TickCount > nextTicks then begin
					UpdateMeter(round(row/maxN * 20.0 ) + 80, 'Computing Power Spectrum');
					nextTicks := TickCount + updateTicks;
					if CommandPeriod then begin
							beep;
							AbortFFT := true;
							exit(DisplayPSUsingBuffer)
						end;
				end;
			end;
		SwapQuadrants;
		UpdateMeter(-1, 'Hide Meter');
	end;	
	
	
	procedure DisplayPowerSpectrum(fht: rImagePtr);
	var
		row, col, nextTicks, maxN: LongInt;
		r, min, max, scale: real;
		line: lineType;
		rLine: rLineType;
		tempH: handle;
		ps: rImagePtr;
	begin
		maxN := info^.PixelsPerLine;
		tempH := GetBigHandle(maxN * maxN * SizeOf(real));
		if tempH <> nil then begin
			hlock(tempH);
			ps := rImagePtr(tempH^);
			DisplayPSUsingBuffer(fht, ps, rLine);
			hunlock(tempH);
			DisposeHandle(tempH);
			exit(DisplayPowerSpectrum);
		end;
		min := 1e99;
		max := -1e99;
		nextTicks := TickCount + updateTicks;
		for row := 0 to maxN - 1 do begin
			FHTps (row, maxN, fht, rLine);
			for col := 0 to maxN - 1 do begin
					r := rLine[col];
					if r < min then min := r;
					if r > max then max := r;
			end;
			if TickCount > nextTicks then begin
				UpdateMeter(round(row/maxN * 35.0), 'Computing Power Spectrum');
				nextTicks := TickCount + updateTicks;
				if CommandPeriod then begin
						beep;
						AbortFFT := true;
						exit(DisplayPowerSpectrum)
					end;
			end;
		end;
		if min < 1.0 then
			min := 0.0
		else
			min := ln(min);
		max := ln(max);
		scale := 253.0 / (max - min);
		for row := 0 to maxN - 1 do begin
				FHTps (row, maxN, fht, rLine);
				for col := 0 to maxN - 1 do begin
						r := rLine[col];
						if r < 1.0 then
							r := 0.0
						else
							r := ln(r);
						line[col] := round((r - min) * scale) + 1;
				end;
				PutLine(0, row, maxN, line);
				if TickCount > nextTicks then begin
					UpdateMeter(round(row/maxN * 65.0 ) + 35, 'Computing Power Spectrum');
					nextTicks := TickCount + updateTicks;
					if CommandPeriod then begin
							beep;
							AbortFFT := true;
							exit(DisplayPowerSpectrum)
						end;
				end;
			end;
		SwapQuadrants;
		UpdateMeter(-1, 'Hide Meter');
	end;


procedure ConvertToReal;
var
	row, col, i, sum, base: LongInt;
	width, height: LongInt;
	mean, pixelCount: real;
	line: LineType;
	DataP: rImagePtr;
begin
	with info^ do begin
		width := pixelsPerLine;
		height := nLines;
		hlock(DataH);
		DataP := rImagePtr(DataH^);
	end;
	UpdateMeter(0, 'Converting to Real');
	{
	GetRectHistogram;
	sum := 0;
	pixelCount := width * height;
	for i := 0 to 255 do
		sum := sum + histogram[i] * i;
	mean := sum / pixelCount;
	}
	for row:= 0 to height - 1 do begin
		GetLine(0, row, width, line);
		base := row * width;
		for col := 0 to width - 1 do
			DataP^[base + col] := line[col] {- mean};
	end;
	hunlock(info^.DataH);
end;


function isPowerOf2 (x: LongInt): boolean;
var
	i, sum: integer;
begin
sum := 0;
x := abs(x);
for i := 0 to 15 do
	sum := sum + ord(BTST(x, i));
IsPowerOf2 := (sum <= 1);
end;


function PowerOf2Size: boolean;
var
	width, height: LongInt;

	procedure fail;
	begin
 		PutError('A square, power of two size image or selection (128x128, 256x256, etc.) required.');
 		AbortMacro;
 		PowerOf2Size := false;
 		exit(PowerOf2Size);
	end;

begin
	with info^ do begin
		if info = noInfo then
			fail;
		width := pixelsPerLine;
		height := nLines;
		if RoiShowing and (roiType = rectRoi) then with roiRect do begin
			width := right - left;
			height := bottom - top;
		end;
		if not isPowerOf2(width) or (width <> height) then
			fail;
		if width < 4 then begin
				PutMessage('Sequence Length must be >= 4.');
				PowerOf2Size := true;
				exit(PowerOf2Size);
			end;
		PowerOf2Size := true;
	end;
end;


function MakeSinCosTables(maxN: LongInt): boolean;
var
	i: LongInt;
	theta, dTheta: real;
begin
	C := rLinePtr(NewPtr(SizeOf(rLineType)));
	S := rLinePtr(NewPtr(SizeOf(rLineType)));
	if (C = nil) or (S = nil) then begin
		MakeSinCosTables := false;
		PutError('Out of Memory');
		exit(MakeSinCosTables);
	end;
	theta := 0;
	dTheta := 2 * pi / maxN;
	for i := 0 to maxN div 4 - 1 do
		begin
			C^[i] := cos(theta);
			S^[i] := sin(theta);
			theta := theta + dTheta;
		end;
		MakeSinCosTables := true;
end;


function MakeRealImage: boolean;
var
	TempH: handle;
	maxN: LongInt;
begin
	maxN := info^.PixelsPerLine;
	tempH := GetBigHandle(maxN * maxN * SizeOf(real));
	if TempH = nil then begin
		PutMemoryAlert;
		MakeRealImage := false;
		exit(MakeRealImage);
	end;
	if not Duplicate(StringOf('FFT ', nPics + 1:1), false) then begin
		DisposeHandle(TempH);
		exit(MakeRealImage);
	end;
	UpdatePicWindow;
	info^.DataH := tempH;
	ConvertToReal;
	UpdateTitleBar;
	MakeRealImage := true;
end;


procedure SetFFTWindowName(doInverse: boolean);
begin
	with info^ do begin
		if doInverse then
			title := StringOf('Inverse FFT ', picNum:1)
		else
			title := StringOf('FFT ', picNum:1);
		UpdateWindowsMenuItem;
		UpdateTitleBar;
	end;
end;


procedure ApplyFilter(rData: rImagePtr);
var
	row, col, width, height, base, i: LongInt;
	line: LineType;
	passMode: boolean;
	t:FateTable;
begin
	SwapQuadrants;
	with info^ do begin
		width := pixelsPerLine;
		height := nLines;
	end;
	for row:= 0 to height - 1 do begin
		GetLine(0, row, width, line);
		base := row * width;
		for col := 0 to width - 1 do
			rData^[base + col] := line[col]/255.0 * rData^[base + col];
	end;
end;


procedure doMasking(rData: rImagePtr);
var
	row, col, width, height, base, i: LongInt;
	line: LineType;
	passMode: boolean;
	t:FateTable;
begin
	GetRectHistogram;
	if (histogram[0] = 0) and (histogram[255] = 0) then
		exit(doMasking);
	UpdateMeter(0, 'Masking');
	passMode := histogram[255] <> 0;
	if passMode then
		ChangeValues(0,254,0)
	else
		ChangeValues(1,255,255);
	for i := 1 to 3 do
		Filter(UnweightedAvg, 0, t);
	UpdatePicWindow;
	ApplyFilter(rData);
end;


	procedure doFFT(fftKind: fftType);
	var
		startTicks, maxN: LongInt;
		trect: rect;
		RealData: rImagePtr;
		doInverse: boolean;
	begin
			doInverse := fftKind <> ForewardFFT;
			if not PowerOf2Size then
				exit(doFFT);
			startTicks := tickCount;
			if info^.DataH = nil then begin
				if doInverse then begin
					PutError('A real image is required to do an inverse transform.');
 					AbortMacro;
					exit(doFFT);
				end;
				if not MakeRealImage then begin
 					AbortMacro;
					exit(doFFT);
				end
			end else begin
				KillRoi;
				SetFFTWindowName(doInverse);
			end;
			hlock(info^.DataH);
			RealData := rImagePtr(info^.DataH^);
			ShowWatch;
			maxN := info^.PixelsPerLine;
			if not MakeSinCosTables(maxN) then
				exit(doFFT);
			AbortFFT := false;
			ShowMessage(CmdPeriodToStop);
			if doInverse then begin
				if fftKind = InverseFFTWithMask then 
					doMasking(RealData)
				else if fftKind = InverseFFTWithFilter then
					ApplyFilter(RealData);
				rc2DFHT(RealData, true, maxN);
				if not AbortFFT then
					DisplayRealImage(RealData);
			end else begin
				rc2DFHT(RealData, false, maxN);
				if not AbortFFT then
					DisplayPowerSpectrum(RealData);
			end;
			if AbortFFT then
				UpdateMeter(-1, 'Hide');
			hunlock(info^.dataH);
			UpdatePicWindow;
			SetRect(trect, 0, 0, maxN, maxN);
			ShowTime(startTicks, trect, '');
			UpdateWindowsMenuItem;
			DisposePtr(ptr(C));
			DisposePtr(ptr(S));
	end;
	

function isFFT: boolean;
begin
	isFFT := false;
	with info^ do
		if DataH <> nil then
			if pos('FFT', title) = 1 then
				isFFT := true;
end;


procedure RedisplayPowerSpectrum;
	var
		rData: rImagePtr;
	begin
			if info = noInfo then
				exit(RedisplayPowerSpectrum);
			KillRoi;
			if not PowerOf2Size then
				exit(RedisplayPowerSpectrum);
			if not isFFT then begin
					PutError('Real frequency domain image required.');
					exit(RedisplayPowerSpectrum);
				end;
			hlock(info^.DataH);
			rData := rImagePtr(info^.DataH^);
			DisplayPowerSpectrum(rData);
			hunlock(info^.dataH);
			UpdatePicWindow;
	end;


	Procedure doSwapQuadrants;
	begin
		if info = noInfo then
			exit(doSwapQuadrants);
		KillRoi;
		if not PowerOf2Size then
			exit(doSwapQuadrants);
		SetupUndo;
		WhatToUndo := UndoEdit;
		SwapQuadrants;
		UpdatePicWindow;
	end;
	
	
	function arctan2 (x, y: extended): extended;
{ returns angle in the correct quadrant }
	begin
		if x = 0 then
			x := 1E-30; { Could be improved }
		if x > 0 then
			if y >= 0 then
				arctan2 := arctan(y / x)
			else
				arctan2 := arctan(y / x) + 2 * pi
		else
			arctan2 := arctan(y / x) + pi;
	end;
	

	procedure ShowFFTValues (hloc, vloc, ivalue: LongInt);
		var
			tPort: GrafPtr;
			hstart, vstart: integer;
			r, theta, center: extended;
	begin
		with info^ do
			begin
				hstart := InfoHStart;
				vstart := InfoVStart;
				GetPort(tPort);
				SetPort(InfoWindow);
				TextSize(9);
				TextFont(Monaco);
				TextMode(SrcCopy);
				if hloc < 0 then
					hloc := -hloc;
				center := pixelsPerLine div 2;
				r := sqrt(sqr(hloc - center) + sqr(vloc - center));
				theta := arctan2(hloc - center, center - vloc);
				theta := theta * 180 / pi;
				MoveTo(xValueLoc, vstart);
				if SpatiallyCalibrated then begin
						DrawReal(pixelsPerLine / r / xScale, 6, 2);
						DrawString(xUnit);
						DrawString('/c ');
						DrawString('(');
						DrawReal(hloc - center, 4, 0);
						DrawString(')');
					end else begin
						DrawReal(pixelsPerLine / r, 6, 2);
						DrawString('p/c  ');
						DrawString('(');
						DrawReal(hloc - center, 4, 0);
						DrawString(')');
					end;
				DrawString('    ');
				vloc := PicRect.bottom - vloc - 1;
				if vloc < 0 then
					vloc := -vloc;
				MoveTo(yValueLoc, vstart + 10);
				DrawReal(theta, 6, 2);
				TextMode(srcOr);
				DrawString('¡    ');
				TextMode(srcCopy);
				DrawString('(');
				DrawReal(vloc - center + 1, 4, 0);
				DrawString(')');
				DrawString('    ');
				MoveTo(zValueLoc, vstart + 20);
				if fit <> uncalibrated then
					begin
						DrawReal(cvalue[ivalue], 6, 2);
						DrawString(' (');
						DrawLong(ivalue);
						DrawString(')');
					end
				else
					DrawLong(ivalue);
				DrawString('    ');
				SetPort(tPort);
			end;
	end;


end. {fft Unit}