UNIT PCMD1;  {Edge-preserving median 5x5 filter}

INTERFACE

{$R-}   { range checking off}
{$SC+}  { short circuit AND & OR statements in IF's}
USES memtypes, OSIntf, ToolIntf;

CONST
	MaxPixelsPerLine = 2048;
	WhiteC = 0;
	BlackC = 255;
	{returnCodes}
	No_Change			  = 	 0;
	CMDperiod			  = -10000;
	Changed				  = -20000;
	CMDperiod_and_changed = -30000;
	
TYPE
	 PCmdBlock = RECORD
	 				Primary : ptr; {all changes are made to primary}
					Secondary : ptr; {NIL if it couldn't load, or didn't exist}
					OptionKeyDownCall : boolean;
					ShiftKeyDownCall : boolean;
					NLines : integer;
					PixelsPerLine : integer;
					BytesPerRow : integer;
					sPixelsPerLine : integer;
					sBytesPerRow : integer;
					sNLines : integer;
					RoiRect : rect;
					Rectangular : boolean; {T/F}
					mask : ptr;
					ReturnCode : OSErr; {0=okay, <0 = changed,CommandPeriod, or both}
					FgColor : integer;
					iScale:extended;   {extended - # of pixels per unit}
					
    				iunits:string[2];  {2 characters - hold the current 
										calibration units}
					MyOSPort : GrafPtr;
	 			 END;
	 PCmdBlockPtr = ^PcmdBlock;
				 
  	UnsignedByte=0..255;
  	LineType=PACKED ARRAY[0..MaxPixelsperLine] OF UnsignedByte;
	SArray=Array[1..25] of integer;
	
				 
{ MAIN entry point}


PROCEDURE  DoPCommand (pData1:PCmdBlockptr);

IMPLEMENTATION

PROCEDURE DoPCommand (pData1:PCmdBlockptr);
VAR
	pData : PCMDBlock;
	
	PROCEDURE dumbBeep;
	BEGIN
		sysbeep(1);
	END;
	
	FUNCTION blankLine : linetype;
	VAR i:integer;
		jl : linetype;
	BEGIN
		for i:=0 to MaxPixelsPerLine-1 do
			jl[i]:=whiteC;
		blankline:=jl;
	END; {blankline}
	
	PROCEDURE GetLine(h,v,count:integer; VAR line:LineType);
	VAR
		offset:LongInt;
		p:ptr;
	BEGIN
		IF pdata.primary=NIL then
			BEGIN
				line:=BlankLine;
				exit(GetLine);
			END;			
		IF (h<0) OR (v<0) OR
		   ((h+count)>pdata.PixelsPerLine)
		   OR (v>=pdata.nlines) THEN
			BEGIN
				line:=BlankLine;
				exit(GetLine);
			END;
		offset:=LongInt(v)*pdata.BytesPerRow+h;
		p:=ptr(ord4(pdata.Primary)+offset);
		BlocKMove(p,@line,count);
	END; {getline}

	PROCEDURE GetBGLine(h,v,count:integer; VAR line:LineType);
	VAR
		offset:LongInt;
		p:ptr;
	BEGIN
		IF pdata.secondary=NIL then
			BEGIN
				line:=BlankLine;
				exit(GetBGLine);
			END;			
		IF (h<0) OR (v<0) OR
		   ((h+count)>pdata.sPixelsPerLine)
		   OR (v>=pdata.snlines) THEN
			BEGIN
				line:=BlankLine;
				exit(GetBGLine);
			END;
		offset:=LongInt(v)*pdata.sBytesPerRow+h;
		p:=ptr(ord4(pdata.secondary)+offset);
		BlocKMove(p,@line,count);
	END; {getBGline}

	PROCEDURE GetMaskLine(h,v,count:integer; VAR line:LineType);
	VAR
		offset:LongInt;
		p:ptr;
	BEGIN
		IF pdata.mask=NIL then
			BEGIN
				line:=BlankLine;
				exit(GetMaskLine);
			END;			
		IF (h<0) OR (v<0) OR
		   ((h+count)>pdata.PixelsPerLine)
		   OR (v>=pdata.nlines) THEN
			BEGIN
				line:=BlankLine;
				exit(GetMaskLine);
			END;
		offset:=LongInt(v)*pdata.BytesPerRow+h;
		p:=ptr(ord4(pdata.mask)+offset);
		BlocKMove(p,@line,count);
	END; {getMaskline}
	
	PROCEDURE PutLine(h,v,count:integer; VAR line:LineType);
	VAR
		offset:LongInt;
		p:ptr;
	BEGIN
		if pdata.primary=NIL then exit(putLine);
		IF (h<0) OR (v<0) OR (v>=pdata.nlines) THEN exit(PutLine);
		IF (h+count)>pdata.PixelsPerLine THEN count:=pdata.PixelsPerLine-h;
		offset:=LongInt(v)*pdata.BytesPerRow+h;
		p:=ptr(ord4(pdata.primary)+offset);
		BlocKMove(@line,p,count);
	END;{PutLine}

	PROCEDURE GetLineUsingMask(h,v,count:integer; var line:linetype;  padcolor:integer);
	VAR {only needed when eroding/dilating & the boundaries are 
	     assumed to be a special color}
		line2 : linetype;
		i: integer;
	BEGIN
		getline (h,v,count,line);
		getline(h,v,count,line2);
		FOR i:=0 to count-1 DO
			if line2[i]<>pdata.fgcolor THEN line[i]:=padcolor;
	END;{GetLineUsingMask}
	
	PROCEDURE PutLineUsingMask(h,v,count:integer; VAR line:LineType);
	VAR
		aLine,MaskLine:LineType;
		i:integer;
	BEGIN
		GetLine(h,v,count,aline);
		GetMaskLine(h,v,count,MaskLine);
		FOR i:=0 TO count-1 DO
		IF MaskLine[i]=pdata.fgcolor THEN aLine[i]:=line[i];
		PutLine(h,v,count,aLine);
	END;

	FUNCTION CommandPeriod:boolean;
	TYPE
		KeyPtrType=^KeyMap;
	VAR
		KeyPtr:KeyPtrType;
		keys:ARRAY[0..3] OF LongInt;
		event : eventrecord;
	BEGIN
	  systemtask;
	  KeyPtr:=KeyPtrType(@keys);
	  GetKeys(KeyPtr^);
	  CommandPeriod:=(BAND(keys[1],$808000))=$808000;
	END;

	PROCEDURE QuickSort(VAR X:SArray; N:Integer);	
		PROCEDURE QSort(VAR X:Sarray; M,N:integer);
		VAR
			i,j,hold:Integer;			
			PROCEDURE Partit(VAR A:SArray; VAR i,j:integer; left,right:integer);
			VAR
				pivot:integer;				
				PROCEDURE Swap(Var P,Q:integer);
				VAR
					Hold:integer;
				BEGIN
					Hold:=P;
					P:=Q;
					Q:=Hold;
				END; {Swap}
			BEGIN {Partit}
				Pivot:=A[(Left+Right) Div 2];
				i:=Left;
				j:=Right;
				While i<=j DO
					BEGIN
						While A[i]<Pivot DO
							i:=i+1;
						While Pivot<A[j] DO
							j:=j-1;
						IF i<=j THEN
							BEGIN
								{Swap(A[i],A[j]);}
								Hold:=A[i];
								A[i]:=A[j];
								A[j]:=Hold;
								i:=i+1;
								j:=j-1;
							END;
					END;{while}
			END;{Partit}
		BEGIN {QSort}
			IF M<N THEN
				BEGIN
					Partit(X,I,J,M,N); {divide in two}
					QSort(X,M,J); {sort left part}
					QSort(X,i,N); {Sort right part}
				END;
		END; {QSort}
	BEGIN{QuickSort}
		QSort(X,1,N);
	END;{QuickSort}

	PROCEDURE ShellSort(VAR A:SArray; N:Integer);
	VAR
		Done:Boolean;
		Jump,i,j,Hold:integer;
	BEGIN
		Jump:=N;
		While Jump>1 DO
			BEGIN
				Jump:=Jump DIV 2;
				Repeat
					Done:=True;
					For j:=1 to n-Jump DO
						BEGIN
							i:=j+Jump;
							IF A[j]>A[i] THEN
								BEGIN
									Hold:=A[j];
									A[j]:=A[i];
									A[i]:=Hold;
									Done:=False;
								END;
						END;
				Until Done;
			END;
	END;

VAR
	i,j : integer; {for scanning entire area}
	irect : rect;
	L1,L2,L3,L4,L5 : linetype; {space for 5 incoming lines}
	Result : linetype;         {result from processing}
	ptr1,ptr2,ptr3,ptr4,ptr5, shuffle : ^linetype;  {pointing at all 5 lines}
	width,height : integer;
	sortarray,medarray : SArray {array[1..9] of unsignedbyte};
	k,m,swap : integer; {for bubblesort}
	{Squares: Array[0..255] of extended;
	SumB:Longint;
	mean,variance,temp,sumB2,factor,scale:extended;}
	
BEGIN {DoMedian5x5}
	pdata:=pdata1^;
	pdata.ReturnCode:=changed; {if the image is going to change}
	
	irect:=pdata.roirect; {bottom is not inclusive in process area}
	with irect,pdata DO
		BEGIN
			if top<2 then top:=2; {avoid top 2 lines}
			if left<2 then left:=2;
			if bottom>nlines-2 then bottom:=nlines-2;
			if right>PixelsperLine-2 then right:=PixelsPerLine-2;
			width:=right-left;
			height:=bottom-top;
		END;
	
	{set up top 4 lines before entering loop}
	
	ptr1:=@L1;ptr2:=@L2;ptr3:=@L3;ptr4:=@L4;ptr5:=@L5; shuffle:=NIL;
	
			GetLine(irect.left-2,irect.top-2,width+4,ptr1^);
			GetLine(irect.left-2,irect.top-1,width+4,ptr2^);
			GetLine(irect.left-2,irect.top  ,width+4,ptr3^);
			GetLine(irect.left-2,irect.top+1,width+4,ptr4^);
		
	{ top of the loop}
	
	FOR i:=irect.top to irect.bottom-1 DO
		BEGIN
			GetLine(irect.left-2,i+2,width+4,ptr5^);
			
			for j:=2 to width+1 DO
				BEGIN { *** process each point, w/ Result[] ***}
					MedArray[1]:=ptr3^[j  ]; {central pixel}
				
				{First, sort the X pattern to get its median value}
					SortArray[ 1]:=ptr1^[j-2];
					SortArray[ 2]:=ptr1^[j+2];
					SortArray[ 3]:=ptr2^[j-1];
					SortArray[ 4]:=ptr2^[j+1];
					SortArray[ 5]:=ptr3^[j  ];
					SortArray[ 6]:=ptr4^[j-1];
					SortArray[ 7]:=ptr4^[j+1];
					SortArray[ 8]:=ptr5^[j-2];
					SortArray[ 9]:=ptr5^[j-2];

					ShellSort(SortArray,9);
					MedArray[2]:=SortArray[5];
								
				{Next, sort the + pattern to get its median value}
					SortArray[ 1]:=ptr1^[j  ];
					SortArray[ 2]:=ptr2^[j  ];
					SortArray[ 3]:=ptr3^[j-2];
					SortArray[ 4]:=ptr3^[j-1];
					SortArray[ 5]:=ptr3^[j  ];
					SortArray[ 6]:=ptr3^[j+1];
					SortArray[ 7]:=ptr3^[j+2];
					SortArray[ 8]:=ptr4^[j  ];
					SortArray[ 9]:=ptr5^[j  ];

					ShellSort(SortArray,9);
					MedArray[3]:=SortArray[5];
					
					ShellSort(MedArray,3);

					Result[j-2]:=MedArray[2]; {Median value}
					
					{JPL Speckle}
					{SumB:=0;
					SumB2:=0;
					For  k:=j-2 to j+2 DO
						BEGIN
							SumB:=SumB+longint(ptr1^[k])+longint(ptr2^[k])
								+longint(ptr3^[k])+longint(ptr4^[k])
								+longint(ptr5^[k]);
							SumB2:=SumB2+squares[ptr1^[k]]+squares[ptr2^[k]]
								+squares[ptr3^[k]]+squares[ptr4^[k]]
								+squares[ptr5^[k]];
						END;
					mean:=SumB/25/255;
					variance:=ABS(SumB2/25-mean*mean);
					temp:=ptr3^[j]/255;
					scale:=variance/36*factor;
					temp:=temp-(scale)*(mean-temp);
					if temp>1 then Result[j-2]:=255
						ELSE if temp<0 THEN Result[j-2]:=0
							ELSE Result[j-2]:=trunc(255*temp);}
				END; {for j}
				
			if pdata.rectangular THEN
				PutLine(irect.left,i,width,result)
			ELSE
				PutLineUsingMask(irect.left,i,width,result);
			
			shuffle:=ptr1;
			ptr1:=ptr2;
			ptr2:=ptr3;
			ptr3:=ptr4;
			ptr4:=ptr5;
			ptr5:=shuffle;
			
			IF commandperiod then
				BEGIN
					sysbeep(1);
					pdata.ReturnCode:=pdata.ReturnCode+CMDperiod;
					pdata1^:=pdata;
					Exit(DoPCommand);
				END;
		END; {for i}
		
		
		pData1^:=pdata;
	
END; {DoPCommand}



END. {UNIT}