unit Math;

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, RealUtils, Graphics, Camera, Filters, Lut, fft;


	procedure SetPasteMode (item: integer);
	procedure DoMouseDownInPasteControl (loc: point);
	procedure ShowPasteControl;
	procedure DrawPasteControl;
	procedure DoArithmetic (MenuItem: integer; constant: extended);
	function GetMathRoi(Src1PicNum, Src2PicNum:integer; var roi:rect):boolean;
	procedure DoMath (Src1PicNum, Src2PicNum: integer; DstInfo:InfoPtr; roi:rect);
	procedure DoPasteMath;
	procedure DoImageMath;
	function GetInfoPtr (PicN: integer): InfoPtr;


implementation

	const
		Src1Item = 7;
		Src2Item = 8;
		OpItem = 9;


	procedure DoPasteMath;
		const
			PixelsPerUpdate = 15000;
		var
			nrows, ncols, hSrcStart, vSrcStart, hDstStart, vDstStart: integer;
			SaveInfo: InfoPtr;
			h, v, vDst, PixelCount, offset: integer;
			Src, Dst: LineType;
			tmp, range, min, max, StartTicks: LongInt;
			x, xmax, xmin, xrange, xxscale: extended;
	begin
		if TooWide then
			exit(DoPasteMath);
		ShowWatch;
		OpPending := false;
		WhatToUndo := UndoPaste;
		KillRoi;
		with info^.RoiRect do begin
				ncols := right - left;
				nrows := bottom - top;
				hDstStart := left;
				vDstStart := top;
			end;
		with ClipBufInfo^.RoiRect do begin
				hSrcStart := left;
				vSrcStart := top;
			end;
		if hDstStart < 0 then begin
				offset := -hDstStart;
				hDstStart := 0;
				hSrcStart := hSrcStart + offset;
				ncols := ncols - offset;
			end;
		if vDstStart < 0 then begin
				offset := -vDstStart;
				vDstStart := 0;
				vSrcStart := vSrcStart + offset;
				nrows := nrows - offset;
			end;
		with info^.PicRect do begin
				if hDstStart + ncols > right then
					ncols := right - hDstStart;
				if vDstStart + nrows > bottom then
					nrows := bottom - vDstStart;
			end;
		SaveInfo := info;
		vDst := vDstStart;
		min := 999999;
		max := -999999;
		xmin := 999999.0;
		xmax := -999999.0;
		StartTicks := TickCount;
       {First pass to find result range}
		if ScaleArithmetic then begin
				for v := vSrcStart to vSrcStart + nRows - 1 do begin
						Info := ClipBufInfo;
						GetLine(hSrcStart, v, nCols, Src);
						Info := SaveInfo;
						GetLine(hDstStart, vDst, nCols, Dst);
						case CurrentOp of
							AddOp:  begin
									for h := 0 to nCols - 1 do begin
											tmp := Src[h] + Dst[h];
											if tmp > max then
												max := tmp;
											if tmp < Min then
												min := tmp;
										end;
								end;
							SubtractOp:  begin
									for h := 0 to nCols - 1 do begin
											tmp := Dst[h] - Src[h];
											if tmp > max then
												max := tmp;
											if tmp < Min then
												min := tmp;
										end;
								end;
							MultiplyOp:  begin
									for h := 0 to nCols - 1 do begin
											tmp := Dst[h];
											tmp := tmp * Src[h];
											if tmp > max then
												max := tmp;
											if tmp < min then
												min := tmp;
										end;
								end;
							DivideOp:  begin
									for h := 0 to nCols - 1 do begin
											tmp := Src[h];
											if tmp = 0 then
												tmp := 1;
											x := Dst[h] / tmp;
											if x > xmax then begin
													xmax := x;
												end;
											if x < xmin then
												xmin := x;
										end;
								end;
						end;
						vDst := vDst + 1;
					end;
				vDst := vDstStart;
				if CurrentOp = DivideOp then begin
						xrange := xmax - xmin;
						if xrange <> 0.0 then
							xxscale := 256.0 / xrange
						else
							xxscale := 1;
					end
				else
					range := max - min;
			end; {if ScaleArithmetic=true}
		PixelCount := 0;
       {Second pass to do arithmetic and scaling}
		for v := vSrcStart to vSrcStart + nRows - 1 do begin
				Info := ClipBufInfo;
				GetLine(hSrcStart, v, nCols, Src);
				Info := SaveInfo;
				GetLine(hDstStart, vDst, nCols, Dst);
				case CurrentOp of
					AddOp: 
						if ScaleArithmetic then
							for h := 0 to nCols - 1 do begin
									tmp := Dst[h] + Src[h] - min;
									if range <> 0 then
										tmp := tmp * 256 div range
									else
										tmp := BackgroundIndex;
									if tmp > 255 then
										dst[h] := 255
									else
										dst[h] := tmp;
								end
						else
							for h := 0 to nCols - 1 do begin
									tmp := Dst[h] + Src[h];
									if tmp > 255 then
										dst[h] := 255
									else
										dst[h] := tmp;
								end;
					SubtractOp: 
						if ScaleArithmetic then
							for h := 0 to nCols - 1 do begin
									tmp := Dst[h] - Src[h] - min;
									if range <> 0 then
										tmp := tmp * 256 div range
									else
										tmp := BackgroundIndex;
									if tmp > 255 then
										dst[h] := 255
									else
										dst[h] := tmp;
								end
						else
							for h := 0 to nCols - 1 do begin
									tmp := Dst[h] - Src[h];
									if tmp < 0 then
										dst[h] := 0
									else
										dst[h] := tmp;
								end;
					MultiplyOp: 
						if ScaleArithmetic then
							for h := 0 to nCols - 1 do begin
									tmp := Dst[h];
									tmp := tmp * Src[h] - min;
									if range <> 0 then
										tmp := tmp * 256 div range
									else
										tmp := BackgroundIndex;
									if tmp > 255 then
										dst[h] := 255
									else
										dst[h] := tmp;
								end
						else
							for h := 0 to nCols - 1 do begin
									tmp := Dst[h];
									tmp := tmp * Src[h];
									if tmp > 255 then
										dst[h] := 255
									else
										dst[h] := tmp;
								end;
					DivideOp: 
						if ScaleArithmetic then
							for h := 0 to nCols - 1 do begin
									tmp := Src[h];
									if tmp = 0 then
										tmp := 1;
									x := Dst[h] / tmp - xmin;
									if xrange <> 0.0 then
										tmp := trunc(x * xxscale)
									else
										tmp := BackgroundIndex;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									dst[h] := tmp;
								end
						else
							for h := 0 to nCols - 1 do begin
									tmp := Src[h];
									if tmp = 0 then
										tmp := 1;
									dst[h] := Dst[h] div tmp;
								end;
				end;
				PutLine(hDstStart, vDst, nCols, Dst);
				vDst := vDst + 1;
				PixelCount := PixelCount + ncols;
				if PixelCount > PixelsPerUpdate then begin
						UpdateScreen(info^.RoiRect);
						if CommandPeriod then begin
								UpdateScreen(info^.RoiRect);
								beep;
								exit(DoPasteMath)
							end;
						PixelCount := 0;
					end;
			end;
		with info^ do begin
				ShowTime(StartTicks, RoiRect, '');
				UpdateScreen(RoiRect);
			end;
	end;


	procedure SetPasteMode (item: integer);
		var
			SavePort: GrafPtr;
			BlendColor: rgbColor;
	begin
		if not macro then begin
				SetForegroundColor(BlackIndex);
				SetBackGroundColor(WhiteIndex);
			end;
		case Item of
			CopyModeItem: 
				PasteTransferMode := SrcCopy;
			AndItem: 
				PasteTransferMode := NotSrcBic; {And}
			OrItem: 
				PasteTransferMode := SrcOr;
			XorItem: 
				PasteTransferMode := SrcXor;
			ReplaceItem: 
				PasteTransferMode := Transparent;
			BlendItem:  begin
					GetPort(SavePort);
					with BlendColor do begin
							red := 32767;
							blue := 32767;
							green := 32767;
						end;
					SetPort(GrafPtr(info^.osPort));
					OpColor(BlendColor);
					SetPort(SavePort);
					PasteTransferMode := Blend;
				end;
			otherwise
		end; {case}
	end;


	function GetTransferModeItem: integer;
	begin
		case PasteTransferMode of
			SrcCopy: 
				GetTransferModeItem := CopyModeItem;
			NotSrcBic: 
				GetTransferModeItem := AndItem;
			SrcOr: 
				GetTransferModeItem := OrItem;
			SrcXor: 
				GetTransferModeItem := XorItem;
			Transparent: 
				GetTransferModeItem := ReplaceItem;
			Blend: 
				GetTransferModeItem := BlendItem;
		end;
	end;


	procedure DrawPasteControl;
		const
			bWidth = 64;
			bHeight = 14;
			vinc = 18;
			bhloc = 114;
			bvloc = 6;
		var
			tPort: GrafPtr;
			i, hloc, vloc, item: integer;
			tType: pcItemType;
			tRect, TriangleRect: rect;
			ItemStr: str255;
	begin
		GetPort(tPort);
		SetPort(PasteControl);
		with PcItem[1] do begin
				SetRect(r, 15, 22, 95, 40);
				itype := pcPopupMenu;
				str := 'Transfer Mode';
			end;
		with pcItem[2] do begin
				SetRect(r, 88, 50, 100, 62);
				itype := pcCheckBox;
				str := 'Scale Math';
			end;
		with pcItem[3] do begin
				SetRect(r, 88, 65, 100, 77);
				itype := pcCheckBox;
				str := 'Live Paste';
			end;
		hloc := bhloc;
		vloc := bvloc;
		tType := pcButton;
		with pcItem[4] do begin
				SetRect(r, hloc, vloc, hloc + bWidth, vloc + bHeight);
				itype := tType;
				str := 'Add';
			end;
		vloc := vloc + vinc;
		with pcItem[5] do begin
				SetRect(r, hloc, vloc, hloc + bWidth, vloc + bHeight);
				itype := tType;
				str := 'Subtract';
			end;
		vloc := vloc + vinc;
		with pcItem[6] do begin
				SetRect(r, hloc, vloc, hloc + bWidth, vloc + bHeight);
				itype := tType;
				str := 'Multiply';
			end;
		vloc := vloc + vinc;
		with pcItem[7] do begin
				SetRect(r, hloc, vloc, hloc + bWidth, vloc + bHeight);
				itype := tType;
				str := 'Divide';
			end;
		TextFont(SystemFont);
		TextSize(12);
		for i := 1 to npcItems do
			with pcItem[i] do
				case iType of
					pcPopupMenu: 
						with r do begin
								MoveTo(r.left - 10, r.top - 4);
								DrawString(str);
								DrawDropBox(r);
								item := GetTransferModeItem;
								GetMenuItemText(TransferModeMenuH, item, ItemStr);
								MoveTo(left + 13, bottom - 5);
								DrawString(ItemStr);
							end;
					pcCheckBox: 
						with r do begin
								MoveTo(left - StringWidth(str) - 4, bottom - 2);
								DrawString(str);
								EraseRect(r);
								FrameRect(r);
								if ((i = 2) and ScaleArithmetic) or ((i = 3) and LivePasteMode) then begin
										MoveTo(left, top);
										LineTo(right - 1, bottom - 1);
										MoveTo(left, bottom - 1);
										LineTo(right - 1, top);
									end;
							end;
					pcButton:  begin
							FrameRoundRect(r, 6, 6);
							with r do
								MoveTo(left + ((right - left) - StringWidth(str)) div 2, bottom - 3);
							DrawString(str);
						end;
				end; {case}
		SetPort(tPort);
	end;


	procedure DoMouseDownInPasteControl; {(loc:point)}
		var
			nItem, i, MenuItem: integer;
			tr: rect;
	begin
		if not (OpPending and (CurrentOp = PasteOp)) then begin
				PutError('Paste Control is only available during paste operations.');
				exit(DoMouseDownInPasteControl);
			end;
		SetPort(PasteControl);
		GlobalToLocal(loc);
		nItem := 0;
		for i := 1 to npcItems do
			if PtInRect(loc, pcItem[i].r) then
				nitem := i;
		if nItem > 0 then begin
				case pcItem[nItem].itype of
					pcPopUpMenu: 
						with pcItem[1].r do begin
								MenuItem := PopUpMenu(TransferModeMenuH, left, top, GetTransferModeItem);
								SetPasteMode(MenuItem);
							end;
					pcCheckBox:  begin
							tr := pcItem[nItem].r;
							InsetRect(tr, 1, 1);
							FrameRect(tr);
							if nitem = 2 then
								ScaleArithmetic := not ScaleArithmetic;
							if nitem = 3 then begin
									LivePasteMode := not LivePasteMode;
									if LivePasteMode then begin
											ExternalTrigger := false;
											UpdateVideoControl
										end;
								end;
						end;
					pcButton:  begin
							InvertRoundRect(pcItem[nitem].r, 6, 6);
							while Button and (nitem > 0) do begin
									GetMouse(loc);
									if not PtInRect(loc, pcItem[nitem].r) then begin
											InvertRoundRect(pcItem[nitem].r, 6, 6);
											nItem := 0;
										end;
								end;
						end;
				end; {case}
				repeat
				until not button;
				if nItem > 0 then
					with pcItem[nitem] do begin
							case itype of
								pcPopupMenu: 
									;
								pcCheckBox:  begin
									end;
								pcButton:  begin
										InvertRoundRect(pcItem[nitem].r, 6, 6);
										if info^.RoiType = RectRoi then begin
												case nitem of
													4: 
														CurrentOp := AddOp;
													5: 
														CurrentOp := SubtractOp;
													6: 
														CurrentOp := MultiplyOp;
													7: 
														CurrentOp := DivideOp;
												end;
												DoPasteMath;
											end; {if}
									end; {pcButton}
							end; {case}
						end; {with}
			end; {if nitem>0}
		if LivePasteMode and (((WhatsOnClip <> CameraPic) and (WhatsOnClip <> LivePic)) or (FrameGrabber =NoFrameGrabber)) then begin
				PutError('"Live Paste" requires that a rectangular selection be first copied from the Camera window to the Clipboard.');
				LivePasteMode := false;
			end;
		if LivePasteMode and (info^.PictureType = FrameGrabberType) then begin
				PutError('Live pasting into the Camera window is not supported.');
				LivePasteMode := false;
			end;
		DrawPasteControl;
	end;


	procedure ShowPasteControl;
		var
			tPort: GrafPtr;
			trect: rect;
			wp: ^WindowPtr;
	begin
		SetRect(trect, PasteControlLeft, PasteControlTop, PasteControlLeft + pcwidth, PasteControlTop + pcheight);
		PasteControl := NewWindow(nil, trect, 'Paste Control', true, rDocProc, nil, true, 0);
		WindowPeek(PasteControl)^.WindowKind := PasteControlKind;
		wp := pointer(GhostWindow);
		wp^ := PasteControl;
		PasteTransferMode := SrcCopy;
		LivePasteMode := false;
	end;


	function GetRealC (message: str255; default: extended; precision: integer; var Canceled: boolean): extended;
		const
			NumberID = 3;
			CalibrateID = 5;
		var
			mylog: DialogPtr;
			item: integer;
	begin
		InitCursor;
		ParamText(message, '', '', '');
		mylog := GetNewDialog(290, nil, pointer(-1));
		SetDReal(MyLog, NumberID, default, precision);
		SetDlogItem(MyLog, CalibrateID, ord(RealArithmetic));
		SelectdialogItemText(MyLog, NumberID, 0, 32767);
		OutlineButton(MyLog, ok, 16);
		repeat
			ModalDialog(nil, item);
			if item = CalibrateID then begin
				RealArithmetic := not RealArithmetic;
				SetDlogItem(MyLog, CalibrateID, ord(RealArithmetic));
			end;
		until (item = ok) or (item = cancel);
		if item = ok then begin
				GetRealC := GetDReal(MyLog, NumberID);
				Canceled := false;
			end
		else begin
				GetRealC := default;
				Canceled := true;
			end;
		DisposeDialog(mylog);
	end;


	procedure DoArithmetic (MenuItem: integer; constant: extended);
		var
			table: LookupTable;
			i, iConst, pNum, digits: integer;
			tmp: LongInt;
			LogScale: extended;
			Canceled: boolean;
			result: str255;
	begin
		canceled := false;
		if info^.fit <> Uncalibrated then
			RealArithmetic := true;
		if not macro then
			case menuItem of
				AddItem: 
					constant := GetRealC('Constant to add:', 25, 0, Canceled);
				SubtractItem: 
					constant := GetRealC('Constant to subtract:', 25, 0, Canceled);
				MultiplyItem:  begin
						constant := GetRealC('Constant to multiply by:', 1.25, 2, Canceled);
						if constant < 0.0 then begin
								PutError('Constant must be positive.');
								exit(DoArithmetic);
							end;
					end;
				DivideItem:  begin
						constant := GetRealC('Constant to divide by:', 1.25, 2, Canceled);
						if constant <= 0.0 then begin
								PutError('Constant must be nonzero and positive.');
								exit(DoArithmetic);
							end;
					end;
				AndItem2: 
					constant := GetInt('AND image with:', 240, Canceled);
				OrItem2: 
					constant := GetInt('OR image with:', 31, Canceled);
				XorItem2: 
					constant := GetInt('XOR image with:', 31, Canceled);
				LogItem:  begin
						if not CheckCalibration then
							exit(DoArithmetic);
						constant := 0.0;
						LogScale := 255.0 / ln(255.0);
					end;
			end; {case}
		if Canceled then
			exit(DoArithmetic);
		if RealArithmetic and (MenuItem >= AddItem) and (MenuItem <= DivideItem)  and not macro then begin
			if constant < 1.0 then
				digits := 3
			else if constant <10.0 then
				digits := 2
			else digits := 1;
			if trunc(constant) = constant then
				digits := 0;
			with info^ do case MenuItem of
				AddItem: begin
						MathGain := 1.0;
						MathOffset := constant;
						result := StringOf(title, '+', constant:1:digits) 
					end;
				SubtractItem: begin
						MathGain := 1.0;
						MathOffset := -constant;
						result := StringOf(title, '-', constant:1:digits) 
					end;
				MultiplyItem: begin
						MathGain := constant;
						MathOffset := 0.0;
						result := StringOf(title, '*', constant:1:digits) 
					end;
				DivideItem: begin
						MathGain := 1.0 / constant;
						MathOffset := 0.0;
						result := StringOf(title, '/', constant:1:digits) 
					end;
			end;
			pNum := info^.picNum;
			with info^.PicRect do
				if not NewRealWindow(result, right-left, bottom-top) then
				exit(DoArithmetic);
			CurrentMathOp := CopyMath;
			RealImageMath := true;
			DoMath(pNum, pNum, Info, info^.picRect);
			exit(DoArithmetic);
		end;
		iconst := trunc(constant);
		for i := 0 to 255 do begin
				case MenuItem of
					AddItem: 
						tmp := round(i + constant);
					SubtractItem: 
						tmp := round(i - constant);
					MultiplyItem: 
						tmp := round(i * constant);
					DivideItem: 
						tmp := round(i / constant);
					AndItem2: 
						tmp := band(i, iconst);
					OrItem2: 
						tmp := bor(i, iconst);
					XorItem2:
						tmp := bxor(i, iconst);
					LogItem:
						if i = 0 then
							tmp := 0
						else
							tmp := round(ln(i) * LogScale);
				end;
				if tmp < 0 then
					tmp := 0;
				if tmp > 255 then
					tmp := 255;
				table[i] := tmp;
			end;
		ApplyTable(table);
		if (MenuItem >= AddItem) and (MenuItem <= DivideItem) then
			RemoveDensityCalibration;
	end;


	function GetInfoPtr (PicN: integer): InfoPtr;
  {Converts a pic number or pid number to an Info ptr.}
		var
			i: integer;
	begin
		i := 0;
		while (PicN < 0) and (i < nPics) do begin
				i := i + 1;
				if InfoPtr(WindowPeek(PicWindow[i])^.RefCon)^.pidNum = PicN then
					PicN := i;
			end;
		if (PicN >= 1) and (PicN <= nPics) then
			GetInfoPtr := pointer(WindowPeek(PicWindow[PicN])^.RefCon)
		else
			GetInfoPtr := nil;
	end;


	procedure ShowRoi(roi:rect);
	begin
		with roi do ShowMessage(StringOf(left:4, top:4, right-left:4, bottom-top:4));
		wait(200);
	end;


	function GetMathRoi(Src1PicNum, Src2PicNum:integer; var roi:rect):boolean;
		var
			Src1Info, Src2Info: InfoPtr;
			ignore:boolean;
	begin
		KillRoi;
		GetMathRoi:=false;
		Src1Info := GetInfoPtr(Src1PicNum);
		Src2Info := GetInfoPtr(Src2PicNum);
 		if (Src1Info = nil) or (Src2Info = nil) then begin
				PutError('Bad pic num or pid num.');
				AbortMacro;
				exit(GetMathRoi);
			end;
		roi := Src1Info^.PicRect;
		ignore := SectRect(roi, Src2Info^.PicRect, roi);
		GetMathRoi:=true;
	end;
	
	

	procedure DoFFTMath (Src1Info, Src2Info, DstInfo:InfoPtr);
	{Does image math when both source images are FFTs.}
	var
		i, StartTicks: LongInt;
		r, c, modAdd, rowMod, colMod, maxN: LongInt;
		H2e, H2o, mag: real;
		h1, h2, Out: rImagePtr;
		roi:rect;
	begin
		maxN := Src1Info^.pixelsPerLine;
		h1 :=  rImagePtr(Src1Info^.DataH^);
		h2 :=  rImagePtr(Src2Info^.DataH^);
		out := rImagePtr(DstInfo^.DataH^);
		StartTicks := TickCount;
		case CurrentMathOp of
			MulMath, cMulMath: {Point by point Hartley multiplication or conjugate multiplication}
				for r := 0 to maxN - 1 do begin
					rowMod := (maxN - r) mod maxN;
					for c := 0 to maxN - 1 do begin
						colMod := (maxN - c) mod maxN;
						H2e := (h2^[r * maxN + c] + h2^[rowMod * maxN + colMod]) / 2;
						H2o := (h2^[r * maxN + c] - h2^[rowMod * maxN + colMod]) / 2;
						if CurrentMathOp = cMulMath then
							Out^[r * maxN + c] := h1^[r * maxN + c] * H2e - h1^[rowMod * maxN + colMod] * H2o
						else
							Out^[r * maxN + c] := h1^[r * maxN + c] * H2e + h1^[rowMod * maxN + colMod] * H2o;
					end;
				end; {for}
			DivMath:  {Point by point Hartley division}
				for r := 0 to maxN - 1 do begin
					rowMod := (maxN - r) mod maxN;
					for c := 0 to maxN - 1 do begin
						colMod := (maxN - c) mod maxN;
						mag := H2^[r * maxN + c] * H2^[r * maxN + c] + H2^[rowMod * maxN + colMod] * H2^[rowMod * maxN + colMod];
						if mag < 1e-20 then
							mag := 1e-20;
						H2e := (H2^[r * maxN + c] + H2^[rowMod * maxN + colMod]);
						H2o := (H2^[r * maxN + c] - H2^[rowMod * maxN + colMod]);
						Out^[r * maxN + c] := (H1^[r * maxN + c] * H2e - H1^[rowMod * maxN + colMod] * H2o) / mag;
					end;
				end; {for}
		end; {case}
		info := dstInfo;
		DisplayPowerSpectrum(out);
		SetFFTWindowName(false);
		hunlock(src1Info^.dataH);
		hunlock(src2Info^.dataH);		
		hunlock(dstInfo^.dataH);
		SetRect(roi, 0, 0, maxN, maxN);
		ShowTime(StartTicks, roi, '');
	end; {DoFFTMath}

	
	
	procedure DoRealRealMath (Src1Info, Src2Info, DstInfo:InfoPtr; roi:rect);
	{Does image math when both source images are real.}
		var
			nrows, ncols, hStart, vStart: LongInt;
			h, v, i: LongInt;
			tmp, tmp1, tmp2, StartTicks: LongInt;
			SaveRow:LongInt;
			NextUpdate: LongInt;
			MaskRect:rect;
			min, max, x, scale, x1, x2: extended;
			s1Data, s2Data, rData: rImagePtr;
			base: LongInt;
			isFFT: boolean;
			saveInfo: InfoPtr;

	begin
		with src1Info^ do begin
			if DataH = nil then begin
				saveInfo := info;
				info := src1Info;
				if not ConvertToReal then begin
					info := saveInfo;
					exit(DoRealRealMath)
				end;
			end;
			hlock(DataH);
			s1Data := rImagePtr(DataH^);
			isFFT := pos('FFT', title) = 1;
		end;
		if (CurrentMathOp <> CopyMath) then
			with src2Info^ do begin
				if DataH = nil then begin
					saveInfo := info;
					info := src2Info;
					if not ConvertToReal then begin
						info := saveInfo;
						hunlock(DataH);
						exit(DoRealRealMath)
					end;
				end;
				hlock(DataH);
				s2Data := rImagePtr(DataH^);
				isFFT := isFFT and (pos('FFT', title) = 1);
			end;
		with DstInfo^ do begin
			hlock(DataH);
			rData := rImagePtr(DataH^);
			isFFT := isFFT and (pixelsPerLine = nlines) and isPowerOf2(pixelsPerLine);
		end;
		if isFFT and (src1Info^.PixelsPerLine = src2Info^.PixelsPerLine)
			and (src1Info^.PixelsPerLine = src1Info^.nLines)
			and (CurrentMathOp in [MulMath, cMulMath, DivMath]) then begin
				doFFTMath(Src1Info, Src2Info, DstInfo);
				exit(DoRealRealMath);
			end;
		with roi do begin
				ncols := right - left;
				nrows := bottom - top;
				hStart := left;
				vStart := top;
			end;
		StartTicks := TickCount;
		NextUpdate:=TickCount+10;
		min:=10e99;
		max:=-10e99;
		
		for v := vStart to vStart + nRows - 1 do begin
			base := v * ncols;
			case CurrentMathOp of
				AddMath: 
					for h := 0 to nCols - 1 do begin
							x := s1Data^[base + h] + s2Data^[base + h];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				SubMath: 
					for h := 0 to nCols - 1 do begin
							x := s1Data^[base + h] - s2Data^[base + h];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				MulMath, cMulMath: 
					for h := 0 to nCols - 1 do begin
							x := s1Data^[base + h] * s2Data^[base + h];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				DivMath: 
					for h := 0 to nCols - 1 do begin
							x := s1Data^[base + h] / s2Data^[base + h];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				MaxMath: 
					for h := 0 to nCols - 1 do begin
							x1 := s1Data^[base + h];
							x2 := s2Data^[base + h];
							if x1 >= x2 then
								x := x1
							else
								x := x2;
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				MinMath: 
					for h := 0 to nCols - 1 do begin
							x1 := s1Data^[base + h];
							x2 := s2Data^[base + h];
							if x1 <= x2 then
								x := x1
							else
								x := x2;
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				CopyMath: 
					for h := 0 to nCols - 1 do begin
							x := s1Data^[base + h];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
			end; {case}
			if TickCount >= NextUpdate then begin
				ShowAnimatedWatch;
				NextUpdate:=TickCount+10;
				if CommandPeriod then begin
						beep;
						AbortMacro;
						exit(DoRealRealMath)
					end;
			end;
		end; {for}
		info := dstInfo;
		if isFFT then begin
			DisplayPowerSpectrum(rData);
			SetFFTWindowName(false);
		end else
			DisplayRealImage(rData, min, max, false);
		hunlock(src1Info^.dataH);
		if CurrentMathOp <> copyMath then	
			hunlock(src2Info^.dataH);		
		hunlock(dstInfo^.dataH);		
		ShowTime(StartTicks, roi, '');
	end; {DoRealRealMath}



	procedure DoRealMath (Src1PicNum, Src2PicNum: integer; DstInfo:InfoPtr; roi:rect);
		var
			nrows, ncols, hStart, vStart: LongInt;
			Src1Info, Src2Info: InfoPtr;
			h, v, i: LongInt;
			src1, src2, dst: LineType;
			tmp, tmp1, tmp2, StartTicks: LongInt;
			SaveRow:LongInt;
			NextUpdate: LongInt;
			MaskRect:rect;
			min, max, x, scale, x1, x2: extended;
			BlackIsZero:boolean;
			cvalue1, cvalue2: array[0..255] of extended;
			rData: rImagePtr;
			base: LongInt;

	begin
		Src1Info := GetInfoPtr(Src1PicNum);
		Src2Info := GetInfoPtr(Src2PicNum);
		ShowWatch;
		if DstInfo^.DataH = nil then begin
			PutError('Output image must be real.');
			exit(DoRealMath);
		end;
		if (src1Info^.dataH <> nil) or (src2Info^.dataH <> nil) then begin
			doRealRealMath(Src1Info, Src2Info, DstInfo, roi);
			exit(DoRealMath);
		end;
		with DstInfo^ do begin
			if DataH = nil then
				exit(DoRealMath);
			hlock(DataH);
			rData := rImagePtr(DataH^);
		end;
		BlackIsZero:=((Src1Info^.fit=StraightLine) and (Src1Info^.coefficient[2]<0))
					or ((Src2Info^.fit=StraightLine) and (Src2Info^.coefficient[2]<0));
		with roi do begin
				ncols := right - left;
				nrows := bottom - top;
				hStart := left;
				vStart := top;
			end;
		info := Src2Info;
		GenerateValues;
		for i:=0 to 255 do cvalue2[i]:=cvalue[i];
		info := Src1Info;
		GenerateValues;
		for i:=0 to 255 do cvalue1[i]:=cvalue[i];
		StartTicks := TickCount;
		NextUpdate:=TickCount+10;
		min:=10e99;
		max:=-10e99;
		
		for v := vStart to vStart + nRows - 1 do begin
			info := Src1Info;
			GetLine(hStart, v, nCols, src1);
			Info := Src2Info;
			GetLine(hStart, v, nCols, src2);
			base := v * ncols;
			case CurrentMathOp of
				AddMath: 
					for h := 0 to nCols - 1 do begin
							x := cvalue1[src1[h]] + cvalue2[src2[h]];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				SubMath: 
					for h := 0 to nCols - 1 do begin
							x := cvalue1[src1[h]] - cvalue2[src2[h]];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				MulMath, cMulMath: 
					for h := 0 to nCols - 1 do begin
							x := cvalue1[src1[h]] * cvalue2[src2[h]];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				DivMath: 
					for h := 0 to nCols - 1 do begin
							x := cvalue1[src1[h]] / cvalue2[src2[h]];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				MaxMath: 
					for h := 0 to nCols - 1 do begin
							x1 := cvalue1[src1[h]];
							x2 := cvalue2[src2[h]];
							if x1 >= x2 then
								x := x1
							else
								x := x2;
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				MinMath: 
					for h := 0 to nCols - 1 do begin
							x1 := cvalue1[src1[h]];
							x2 := cvalue2[src2[h]];
							if x1 <= x2 then
								x := x1
							else
								x := x2;
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
				CopyMath: 
					for h := 0 to nCols - 1 do begin
							x := cvalue1[src1[h]];
							x := x * MathGain + MathOffset;
							rData^[base + h] := x;
							if x < min then
								min := x;
							if x > max then
								max := x;
						end;
			end; {case}
			if TickCount >= NextUpdate then begin
				ShowAnimatedWatch;
				NextUpdate:=TickCount+10;
				if CommandPeriod then begin
						beep;
						AbortMacro;
						exit(DoRealMath)
					end;
			end;
		end; {for}
		info := dstInfo;
		DisplayRealImage(rData, min, max, BlackisZero);
		hunlock(dstInfo^.dataH);		
		ShowTime(StartTicks, roi, '');
	end; {DoRealMath}



	procedure DoMath (Src1PicNum, Src2PicNum: integer; DstInfo:InfoPtr; roi:rect);
		var
			nrows, ncols, hStart, vStart: integer;
			Src1Info, Src2Info: InfoPtr;
			h, v: integer;
			src1, src2, dst: LineType;
			tmp, tmp1, tmp2, StartTicks, scale, ScaledGain: LongInt;
			rtmp,rtmp2: extended;
			DoScaling: boolean;
			SaveRow:integer;
			NextUpdate: LongInt;
			MaskRect:rect;
			IntegerOffset: LongInt;
	begin
		if TooWide then
			exit(DoMath);
		if RealImageMath and not (CurrentMathOp in [AndMath, OrMath, XorMath]) then begin
			DoRealMath(Src1PicNum, Src2PicNum, DstInfo, roi);
			exit(DoMath);
		end;
		Src1Info := GetInfoPtr(Src1PicNum);
		Src2Info := GetInfoPtr(Src2PicNum);
		ShowWatch;
		with roi do begin
				ncols := right - left;
				nrows := bottom - top;
				hStart := left;
				vStart := top;
			end;
		StartTicks := TickCount;
		scale := 10000;
		ScaledGain := round(MathGain * scale);
		IntegerOffset := round(MathOffset);
		DoScaling := (MathGain <> 1.0) or (MathOffset <> 0.0);
		SaveRow:=vStart;
		NextUpdate:=TickCount+6; {Update screen 10 times per second}
		for v := vStart to vStart + nRows - 1 do begin
				info := Src1Info;
				GetLine(hStart, v, nCols, src1);
				Info := Src2Info;
				GetLine(hStart, v, nCols, src2);
				case CurrentMathOp of
					AddMath: 
						for h := 0 to nCols - 1 do begin
								tmp := src1[h] + src2[h];
								tmp := (tmp * ScaledGain) div scale + IntegerOffset;
								if tmp > 255 then
									tmp := 255;
								if tmp < 0 then
									tmp := 0;
								dst[h] := tmp;
							end;
					SubMath: 
						for h := 0 to nCols - 1 do begin
								tmp := src1[h] - src2[h];
								tmp := (tmp * ScaledGain) div scale + IntegerOffset;
								if tmp > 255 then
									tmp := 255;
								if tmp < 0 then
									tmp := 0;
								dst[h] := tmp;
							end;
					MulMath, cMulMath: 
						for h := 0 to nCols - 1 do begin
								tmp := src1[h];
								tmp := tmp * src2[h];
								tmp := (tmp * ScaledGain) div scale + IntegerOffset;
								if tmp > 255 then
									tmp := 255;
								if tmp < 0 then
									tmp := 0;
								dst[h] := tmp;
							end;
					DivMath: 
						for h := 0 to nCols - 1 do begin
								rtmp2:=src2[h];
								rtmp := src1[h] / rtmp2;
								tmp := round(rtmp * MathGain) + IntegerOffset;
								if tmp > 255 then
									tmp := 255;
								if tmp < 0 then
									tmp := 0;
								dst[h] := tmp;
							end;
					AndMath: 
						for h := 0 to nCols - 1 do begin
								tmp := band(src1[h], src2[h]);
								if DoScaling then begin
										tmp := (tmp * ScaledGain) div scale + IntegerOffset;
										if tmp > 255 then
											tmp := 255;
										if tmp < 0 then
											tmp := 0;
									end;
								dst[h] := tmp;
							end;
					OrMath: 
						for h := 0 to nCols - 1 do begin
								tmp := bor(src1[h], src2[h]);
								if DoScaling then begin
										tmp := (tmp * ScaledGain) div scale + IntegerOffset;
										if tmp > 255 then
											tmp := 255;
										if tmp < 0 then
											tmp := 0;
									end;
								dst[h] := tmp;
							end;
					XorMath: 
						for h := 0 to nCols - 1 do begin
								tmp := bxor(src1[h], src2[h]);
								if DoScaling then begin
										tmp := (tmp * ScaledGain) div scale + IntegerOffset;
										if tmp > 255 then
											tmp := 255;
										if tmp < 0 then
											tmp := 0;
									end;
								dst[h] := tmp;
							end;
					MaxMath: 
						for h := 0 to nCols - 1 do begin
								tmp1 := src1[h];
								tmp2 := src2[h];
								if tmp1 >= tmp2 then
									tmp := tmp1
								else
									tmp := tmp2;
								if DoScaling then begin
										tmp := (tmp * ScaledGain) div scale + IntegerOffset;
										if tmp > 255 then
											tmp := 255;
										if tmp < 0 then
											tmp := 0;
									end;
								dst[h] := tmp;
							end;
					MinMath: 
						for h := 0 to nCols - 1 do begin
								tmp1 := src1[h];
								tmp2 := src2[h];
								if tmp1 <= tmp2 then
									tmp := tmp1
								else
									tmp := tmp2;
								if DoScaling then begin
										tmp := (tmp * ScaledGain) div scale + IntegerOffset;
										if tmp > 255 then
											tmp := 255;
										if tmp < 0 then
											tmp := 0;
									end;
								dst[h] := tmp;
							end;
					CopyMath: 
						for h := 0 to nCols - 1 do begin
								tmp := src1[h];
								if DoScaling then begin
										tmp := (tmp * ScaledGain) div scale + IntegerOffset;
										if tmp > 255 then
											tmp := 255;
										if tmp < 0 then
											tmp := 0;
									end;
								dst[h] := tmp;
							end;
				end;
				Info := DstInfo;
				PutLine(0, v - vstart, nCols, Dst);
				if TickCount>=NextUpdate then begin
					SetRect(MaskRect, hStart, SaveRow, hStart+ncols, v + 1);
					UpdateScreen(MaskRect);
					SaveRow:=v+1;
					NextUpdate:=TickCount+6;
					ShowAnimatedWatch;
					if CommandPeriod then begin
							UpdateScreen(info^.RoiRect);
							beep;
							AbortMacro;
							exit(DoMath)
						end;
				end;
			end;
		with info^ do begin
				ShowTime(StartTicks, RoiRect, '');
				SetRect(MaskRect, hStart, SaveRow, hStart+ncols, vStart+nRows);
				UpdateScreen(MaskRect);
				Changes := true;
				RemoveDensityCalibration;
			end;
	end; {DoMath}


	function ImageTitle (var PicNumber: integer): str255;
		var
			TempInfo: InfoPtr;
	begin
		if (PicNumber < 1) or (PicNumber > nPics) then
			PicNumber := 1;
		TempInfo := pointer(WindowPeek(PicWindow[PicNumber])^.RefCon);
		ImageTitle := TempInfo^.title;
	end;


	procedure ImageMathUProc (d: DialogPtr; item: integer);
     {User proc for Image Math dialog box}
		const
			ops = '';
		var
			str: str255;
			VersInfo: str255;
			r: rect;
	begin
		SetPort(d);
		GetDItemRect(d, item, r);
		DrawDropBox(r);
		case item of
			Src1Item: 
				DrawPopUpText(ImageTitle(MathSrc1), r);
			Src2Item: 
				DrawPopUpText(ImageTitle(MathSrc2), r);
			OpItem:  begin
					GetMenuItemText(ImageMathOpsMenuH, ord(CurrentMathOp) + 1, str);
					DrawPopUpText(str, r);
				end;
		end;
	end;


	function PopUpImageList (r: rect; CurrentImage: integer): integer;
		var
			i: integer;
	begin
		for i := 1 to nPics do begin
				AppendMenu(ImageListMenuH, ' ');
				SetMenuItemText(ImageListMenuH, i, ImageTitle(i));
			end;
		PopUpImageList := PopUpMenu(ImageListMenuH, r.left, r.top, CurrentImage);
		for i := 1 to nPics do
			DeleteMenuItem(ImageListMenuH, 1);
	end;


	procedure DoImageMath;
		const
			ScaleItem = 10;
			OffSetMenuItemText = 11;
			ResultItem = 12;
			RealMathItem=14;
		var
			d: DialogPtr;
			item, i, MenuItem: integer;
			r: rect;
			str: str255;
			ScaleOffEdited: boolean;
			roi:rect;
			DstInfo:InfoPtr;

		procedure ShowScaleAndOffset;
		begin
			SetDReal(d, ScaleItem, MathGain, 4);
			if RealImageMath then
				SetDReal(d, OffSetMenuItemText, MathOffset, 2)
			else
				SetDNum(d, OffSetMenuItemText, round(MathOffset));
		end;

		procedure ResetScaleOff;
		begin
			if not ScaleOffEdited then begin
					MathGain := 1.0;
					MathOffset := 0.0;
					ShowScaleAndOffset;
				end;
		end;

		procedure CheckForRealImage(picNum: integer);
		var
			srcInfo:InfoPtr;
		begin
			srcInfo := GetInfoPtr(picNum);
			if srcInfo^.dataH <> nil then begin
				RealImageMath := true;
				SetDlogItem(d, RealMathItem, ord(RealImageMath));
			end;
		end;
	
	begin
		if ImageMathUserProc=nil
			then ImageMathUserProc:=NewRoutineDescriptor(@ImageMathUProc, uppUserItemProcInfo, GetCurrentISA);
		InitCursor;
		ScaleOffEdited := false;
		d := GetNewDialog(200, nil, pointer(-1));
		SetUProc(d, Src1Item, handle(ImageMathUserProc));
		SetUProc(d, Src2Item, handle(ImageMathUserProc));
		SetUProc(d, OpItem, handle(ImageMathUserProc));
		ShowScaleAndOffset;
		SetDString(d, ResultItem, MathResult);
		SetDlogItem(d, RealMathItem, ord(RealImageMath));
		if (MathSrc1 = 1) and (MathSrc2 = 1) then
			MathSrc1 := info^.PicNum;
		if MathSrc1 = MathSrc2 then begin
				if MathSrc1 = info^.PicNum then begin
						MathSrc2 := MathSrc2 + 1;
						if MathSrc2 > nPics then
							MathSrc2 := 1;
					end
				else
					MathSrc2 := info^.PicNum;
			end;
		CheckForRealImage(MathSrc1);
		CheckForRealImage(MathSrc2);
		repeat
			ModalDialog(nil, item);
			if item = Src1Item then begin
					setport(d);
					GetDItemRect(d, item, r);
					MenuItem := PopUpImageList(r, MathSrc1);
					if MenuItem <> 0 then
						MathSrc1 := MenuItem;
					CheckForRealImage(MathSrc1);
					InvalRect(r);
				end;
			if item = Src2Item then begin
					setport(d);
					GetDItemRect(d, item, r);
					MenuItem := PopUpImageList(r, MathSrc2);
					if MenuItem <> 0 then
						MathSrc2 := MenuItem;
					CheckForRealImage(MathSrc2);
					InvalRect(r);
				end;
			if item = OpItem then begin
					setport(d);
					GetDItemRect(d, item, r);
					MenuItem := PopUpMenu(ImageMathOpsMenuH, r.left, r.top, ord(CurrentMathOp) + 1);
					case MenuItem of
						1:  begin
								CurrentMathOp := AddMath;
								if not ScaleOffEdited then begin
										if RealImageMath then
											ResetScaleOff
										else begin
											MathGain := 0.5;
											MathOffset := 0.0;
											ShowScaleAndOffset;
										end;
									end;
							end;
						2:  begin
								CurrentMathOp := SubMath;
								if not ScaleOffEdited then begin
										if RealImageMath then
											ResetScaleOff
										else begin
											MathGain := MathSubGain;
											MathOffset := MathSubOffset;
											ShowScaleAndOffset;
										end;
									end;
							end;
						3:  begin
								CurrentMathOp := MulMath;
								if not ScaleOffEdited then begin
										if RealImageMath then
											ResetScaleOff
										else begin
											MathGain := 1.0 / 255.0;
											MathOffset := 0.0;
											ShowScaleAndOffset;
										end;
									end;
							end;
						4:  begin
								CurrentMathOp := DivMath;
								if not ScaleOffEdited then begin
										if RealImageMath then
											ResetScaleOff
										else begin
											MathGain := 255.0;
											MathOffset := 0.0;
											ShowScaleAndOffset;
										end;
									end;
							end;
						5:  begin
								CurrentMathOp := AndMath;
								ResetScaleOff;
							end;
						6:  begin
								CurrentMathOp := OrMath;
								ResetScaleOff;
							end;
						7:  begin
								CurrentMathOp := XorMath;
								ResetScaleOff;
							end;
						8:  begin
								CurrentMathOp := MaxMath;
								ResetScaleOff;
							end;
						9:  begin
								CurrentMathOp := MinMath;
								ResetScaleOff;
							end;
						10:	begin
								CurrentMathOp := cMulMath;
								ResetScaleOff;
							end;
						11:	begin
								CurrentMathOp := CopyMath;
								ResetScaleOff;
							end;
						otherwise
					end;
					InvalRect(r);
				end;
			if item = ScaleItem then begin
					MathGain := GetDReal(d, ScaleItem);
					ScaleOffEdited := true;
				end;
			if item = OffSetMenuItemText then begin
					MathOffset := GetDReal(d, OffSetMenuItemText);
					ScaleOffEdited := true;
				end;
			if item = RealMathItem then begin
					RealImageMath := not RealImageMath;
					SetDlogItem(d, RealMathItem, ord(RealImageMath));
					ResetScaleOff;
				end;
		until (item = ok) or (item = cancel);
		MathResult := GetDString(d, ResultItem);
		DisposeDialog(d);
		if item = cancel then
			exit(DoImageMath);
		if not GetMathRoi(MathSrc1, MathSrc2, roi) then
			exit(DoImageMath);
		with roi do
			if RealImageMath then begin
				if not NewRealWindow(MathResult, right-left, bottom-top) then
					exit(DoImageMath)
			end else begin
				if not NewPicWindow(MathResult, right-left, bottom-top) then
					exit(DoImageMath)
			end;
		DstInfo := Info;
		DoMath(MathSrc1, MathSrc2, DstInfo, roi);
		if CurrentMathOp=SubMath then begin
			MathSubGain:=MathGain;
			MathSubOffset:=MathOffset;
		end;
	end;



end.