-- DynamicBitsXXXImpl.meta
-- Copyright (C) 1984, 1985, Xerox Corporation.  All rights reserved.
-- Michael Plass, November 6, 1985 9:52:38 am PST
-- PARAMETERS sKernelRadius, fKernelRadius, sModelRadius, fModelRadius, swathSize

DIRECTORY DynamicBits, Basics, ImagerPixelMap, ImagerSample, PixelMapOps;

DynamicBits001Impl: CEDAR PROGRAM
IMPORTS Basics, ImagerPixelMap, ImagerSample, PixelMapOps, DynamicBits
~ BEGIN OPEN DynamicBits;
	
	
	KernelArray: TYPE ~ ARRAY [0..1) OF ARRAY [0..1) OF INTEGER;
	
	SampleBuffer: TYPE ~ ImagerSample.SampleBuffer;
	
	PixelRow: TYPE ~ RECORD[sOrigin, fOrigin: INTEGER, sample: SampleBuffer];
	
	CreatePixelRow: PROC [s, f: INTEGER, size: NAT] RETURNS [PixelRow] ~ {
		buffer: SampleBuffer ~ ImagerSample.NewBuffer[iSize: 1, jSize: size];
		RETURN [[sOrigin: s, fOrigin: f, sample: buffer]]
		};
	
	CopyPixelRow: PROC [old: PixelRow] RETURNS [PixelRow] ~ {
		new: PixelRow ~ CreatePixelRow[old.sOrigin, old.fOrigin, old.sample.jSize];
		PixelMapOps.CopySamples[buffer: new.sample, bi: 0, bj: 0, count: old.sample.jSize, source: old.sample, si: 0, sj: 0];
		RETURN [new]
		};
	
	ModelData: TYPE ~ REF ModelDataRep;
	ModelDataRep: TYPE ~ ARRAY [0..2) OF REF ModelDataItem;
	ModelDataItem: TYPE ~ RECORD [
		intensity: Intensity,
		noisePenalty: NoisePenalty,
		a: KernelArray
		];
	
	ModelParameterFault: PUBLIC ERROR [data: REF] ~ CODE;
	
	neighborhoodBounds: DeviceRectangle ~ [0, 0, 1, 1];
	CheckNeighborhood: PROC [neighborhood: DeviceRectangle] ~ {
		IF neighborhood # neighborhoodBounds THEN ModelParameterFault[$neighborhood];
		};
	
	kernelBounds: DeviceRectangle ~ [0, 0, 1, 1];
	CheckKernel: PROC [bounds: DeviceRectangle] ~ {
		IF bounds # kernelBounds THEN ModelParameterFault[$kernel];
		};
	
	ArrayFromKernel: PROC [kernel: Kernel] RETURNS [kernelArray: KernelArray]  ~ {
		CheckKernel[kernel.Window];
		FOR s: INTEGER IN [0..1) DO
			FOR f: INTEGER IN [0..1) DO
				kernelArray[s][f] _ LOOPHOLE[kernel.GetPixel[s, f], INTEGER];
				ENDLOOP;
			ENDLOOP;
		};
	
	CreatePrinterModel: PROC [neighborhood: DeviceRectangle, printerModel: PrinterModel, kernel: PixelMap] RETURNS [model: Model] ~ {
		bitmap: PixelMap _ ImagerPixelMap.Create[0, neighborhood];
		data: ModelData _ NEW[ModelDataRep];
		kernelArray: KernelArray _ ArrayFromKernel[kernel];
		CheckNeighborhood[neighborhood];
		model _ NEW[ModelRep];
		FOR encoding: CARDINAL IN [0..2) DO
			intensity: Intensity;
			noisePenalty: NoisePenalty;
			scaledArray: KernelArray;
			TRUSTED {
				t: CARDINAL _ Basics.BITSHIFT[encoding, Basics.bitsPerWord-1];
				bitmap _ MakePixelMapFromBits[@t, 1, neighborhoodBounds, bitmap.refRep];
				};
			[intensity, noisePenalty] _ printerModel[bitmap, encoding];
			FOR s: INTEGER IN [0..1) DO
				FOR f: INTEGER IN [0..1) DO
					k: INTEGER _ kernelArray[s][f];
					scaledArray[s][f] _
					IF k>=0 THEN (Basics.LongMult[k, intensity]+127)/255
					ELSE -((Basics.LongMult[-k, intensity]+127)/255);
					ENDLOOP;
				ENDLOOP;
			data[encoding] _ NEW[ModelDataItem _ [intensity, noisePenalty, scaledArray]];
			ENDLOOP;
		model.neighborhood _ neighborhood;
		model.kernel _ kernel.Copy;
		model.data _ data;
		};
	
	DoWithPrinterModel: PROC [model: Model, action: PROC[PrinterModel]] ~ {
		data: ModelData ~ NARROW[model.data];
		printerModel: PROC [bitmap: PixelMap, encoding: CARDINAL] RETURNS [Intensity, NoisePenalty] ~ {
			item: REF ModelDataItem ~ data^[encoding];
			RETURN [item.intensity, item.noisePenalty];
			};
		action[printerModel];
		};
	
	StorePixelRow: PROC [p: PixelRow, pixelMap: PixelMap] ~ {
		fMin: INTEGER _ MAX[p.fOrigin, pixelMap.fOrigin+pixelMap.fMin];
		fMax: INTEGER _ MIN[p.fOrigin+p.sample.jSize, pixelMap.fOrigin+pixelMap.fMin+pixelMap.fSize];
		IF fMax<=fMin THEN RETURN;
		PixelMapOps.PutF[pixelMap: pixelMap, s: p.sOrigin, f: fMin, buffer: p.sample, bi: 0, bj: fMin-p.fOrigin, count: fMax-fMin];
		};
	
	Load: PROC [p: PixelRow, pixelMap: PixelMap] ~ {
		-- Takes care of boundary conditions
		fMin: INTEGER _ MAX[p.fOrigin, pixelMap.fOrigin+pixelMap.fMin];
		fMax: INTEGER _ MIN[p.fOrigin+p.sample.jSize, pixelMap.fOrigin+pixelMap.fMin+pixelMap.fSize];
		sMin: INTEGER _ MIN[MAX[p.sOrigin, pixelMap.sOrigin+pixelMap.sMin], pixelMap.sOrigin+pixelMap.sMin+pixelMap.sSize-1];
		IF fMax<=fMin THEN {PixelMapOps.ClearSamples[buffer: p.sample, count: p.sample.jSize]; RETURN};
		PixelMapOps.GetF[pixelMap: pixelMap, s: sMin, f: fMin, buffer: p.sample, bi: 0, bj: fMin-p.fOrigin, count: fMax-fMin];
		FOR i: INTEGER IN [0..fMin-p.fOrigin) DO
			p.sample[i] _ p.sample[fMin-p.fOrigin];
			ENDLOOP;
		FOR i: INTEGER IN [fMax-p.fOrigin..p.sample.jSize) DO
			p.sample[i] _ p.sample[fMax-1-p.fOrigin];
			ENDLOOP;
		};
	
	Convolve: PROC [pixelMap: PixelMap, kernel: Kernel, unitIntensity: CARDINAL] ~ {
		bounds: DeviceRectangle _ pixelMap.BoundedWindow;
		kernelArray: KernelArray _ ArrayFromKernel[kernel];
		p: ARRAY [0..1) OF PixelRow _ CreatePixelRows[];
		CreatePixelRows: PROC RETURNS[p: ARRAY[0..1) OF PixelRow]~{
			FOR i: INTEGER IN [0..1) DO
				p[i] _ CreatePixelRow[bounds.sMin+i, bounds.fMin+0, bounds.fSize+1-1];
				ENDLOOP;
			};
		q: PixelRow _ CopyPixelRow[p[0]];
		kernelMag: INTEGER _ KernelWeight[kernel];
		halfUnit: CARDINAL _ unitIntensity/2;
		maxMag: INT _ Basics.LongMult[ABS[kernelMag], INT[unitIntensity]]+halfUnit;
		Cd: PROC [i: INTEGER] RETURNS [CARDINAL] ~ INLINE {RETURN [LOOPHOLE[i]]};
		short: BOOLEAN _ maxMag < CARDINAL.LAST;
		FOR ss: INT IN [0..1) DO
			FOR ff: INT IN [0..1) DO
				short _ short AND kernelArray[ss][ff] >= 0;
				ENDLOOP;
			ENDLOOP;
		FOR i: INTEGER IN [0..1) DO Load[p[i], pixelMap] ENDLOOP;
		THROUGH [0..bounds.sSize) DO
			IF short THEN {
				-- Use 16-bit arithemetic for greater speed.
				FOR f: NAT IN [-0..bounds.fSize-0) DO
					pix: CARDINAL _ p[-0].sample[f-0]*Cd[kernelArray[0][0]] +  0;
					q.sample[f] _ (pix+halfUnit)/unitIntensity;
					ENDLOOP;
				}
			ELSE {
				FOR f: NAT IN [-0..bounds.fSize-0) DO
					pix: INT _ 0;
					FOR ss: INT IN [0..1) DO
						FOR ff: INT IN [0..1) DO
							pix _ pix + p[-ss].sample[f-ff]*INT[kernelArray[ss][ff]];
							ENDLOOP;
						ENDLOOP;
					q.sample[f] _ (pix+halfUnit)/unitIntensity;
					ENDLOOP;
				};
			StorePixelRow[q, pixelMap];
			{t: PixelRow _ p[0]; p[1-1] _ t};
			p[1-1].sOrigin _ p[1-1].sOrigin + (1-0);
			Load[p[1-1], pixelMap];
			q.sOrigin _ p[0].sOrigin;
			ENDLOOP;
		};
	
	ApplyModel: PROC [bitmap: PixelMap, model: Model] RETURNS [PixelMap] ~ {
		w: DeviceRectangle _ bitmap.Window;
		gray: PixelMap _ ImagerPixelMap.Create[3, w];
		bitRow: PixelRow _ CreatePixelRow[w.sMin+0, w.fMin+0, w.fSize+1];
		indexRow: PixelRow _ CreatePixelRow[w.sMin, w.fMin, w.fSize];
		destRow: PixelRow _ CopyPixelRow[indexRow];
		data: ModelData _ NARROW[model.data];
		MergeBits: PROC ~ {
			FOR f: NAT IN [0..indexRow.sample.jSize) DO
				k: CARDINAL _ 0;
				FOR ff: NAT IN [f..f+1) DO
					k _ k*2+bitRow.sample[ff];
					ENDLOOP;
				indexRow.sample[f] _ (indexRow.sample[f] * 2 + k) MOD 2;
				ENDLOOP;
			};
		PixelMapOps.ClearSamples[buffer: bitRow.sample, count: bitRow.sample.jSize];
		PixelMapOps.ClearSamples[buffer: indexRow.sample, count: indexRow.sample.jSize];
		PixelMapOps.ClearSamples[buffer: destRow.sample, count: destRow.sample.jSize];
		FOR i: NAT IN [0..1) DO
			Load[bitRow, bitmap];
			MergeBits[];
			bitRow.sOrigin _ bitRow.sOrigin + 1;
			ENDLOOP;
		FOR s: NAT IN [0..w.sSize) DO
			FOR f: NAT IN [0..destRow.sample.jSize) DO
				destRow.sample[f] _ data[indexRow.sample[f]].intensity;
				ENDLOOP;
			StorePixelRow[destRow, gray];
			destRow.sOrigin _ destRow.sOrigin + 1;
			Load[bitRow, bitmap];
			MergeBits[];
			bitRow.sOrigin _ bitRow.sOrigin + 1;
			ENDLOOP;
		RETURN [gray]
		};
	
	TuneSwath: PROC [gray: PixelMap, bitmap: PixelMap, fixedBits: PixelMap, swath: DeviceRectangle, model: Model, scratch: REF] RETURNS [newScratch: REF] ~ {
		
		
		-- Make sure the swath has the right fSize
		xxx: [0..0] _ swath.fSize-1;
		-- sOrigin and fOrigin give the position of the origin of the local swatch in the bitmap coordinate system.
		sOrigin: INTEGER _ swath.sMin-1-0;
		fOrigin: INTEGER _ swath.fMin;
		pixelRow: PixelRow _ CreatePixelRow[sOrigin, fOrigin+0, 1-0];
		swathRow: PixelRow _ CreatePixelRow[sOrigin, fOrigin, 1];
		bit: ARRAY [0..1) OF ARRAY [0..1) OF [0..1];
		-- The current bits from the bitmap within the swatch
		index: ARRAY [0..1) OF ARRAY [0..1) OF [0..2);
		-- The current bits, packed up into indices into the model array
		fixed: ARRAY [0..1) OF ARRAY [0..1) OF [0..1];
		-- A local copy of a piece of fixedBits
		target: ARRAY [0..1) OF INTEGER;
		-- A local copy of a piece of gray
		trial: ARRAY [0..1) OF INTEGER;
		-- The partially computed approximation to the gray
		modelData: ModelData _ NARROW[model.data];
		
		nonFixedBitCount: NAT _ 0;
		SF: TYPE ~ RECORD [s: [0..1), f: [0..1)];
		nonFixedBitIndex: ARRAY [0..(1-0)*1) OF SF;
		
		questionBits: [0..2);
		-- The bits that need to be determined for a given context.
		contextBits: [0..1);
		-- The context bits in the swath.
		prevContextBits: [0..1);
		-- The context bits for the previous subproblem.
		
		SubproblemSolution: TYPE ~ RECORD [
			setting: [0..2),
			-- The setting of questionBits that achieves the lowest total penalty for the swath so far.
			minErrHighHalf: [0..CARDINAL.LAST/2],
			minErrLowHalf: CARDINAL
			-- The lowest total penalty for the swath so far.
			-- Packed so the whole thing fits in two words, since there are a lot of them.
			];
		
		MakeSubproblemSolution: PROC [setting: [0..2), minErr: INT] RETURNS [SubproblemSolution] ~ INLINE BEGIN
			long: Basics.LongNumber _ [li[minErr]];
			RETURN [[setting: setting, minErrHighHalf: long.highbits, minErrLowHalf: long.lowbits]]
			END;
		
		MinErr: PROC [s: SubproblemSolution] RETURNS [INT] ~ INLINE BEGIN
			long: Basics.LongNumber _ [num[lowbits: s.minErrLowHalf, highbits: s.minErrHighHalf]];
			RETURN [long.li]
			END;
		
		PerSwatchTable: TYPE ~ ARRAY [0..1) OF SubproblemSolution;
		
		avail: LIST OF REF PerSwatchTable _ NARROW[scratch];
		
		subproblemTable: LIST OF REF PerSwatchTable _ NIL;
		
		Alloc: PROC RETURNS [new: LIST OF REF PerSwatchTable] ~ BEGIN
			IF avail = NIL THEN new _ LIST[NEW[PerSwatchTable]]
			ELSE {
				new _ avail;
				avail _ avail.rest;
				new.rest _ NIL;
				};
			new.first^ _ ALL[[0, CARDINAL.LAST/2, CARDINAL.LAST]];
			END;
		
		ComputeFixedIndex: PROC ~ BEGIN
			nonFixedBitCount _ 0;
			FOR s: [0..1) IN [0..1) DO
				FOR f: [0..1) IN [0..1) DO
					IF fixed[s][f] = 0 THEN {
						nonFixedBitIndex[nonFixedBitCount] _ [s,f];
						nonFixedBitCount _ nonFixedBitCount + 1;
						};
					ENDLOOP;
				ENDLOOP;
			END;
		
		Initialize: PROC ~ BEGIN
			FOR s: [0..1) IN [0..1) DO
				pixelRow.sOrigin _ sOrigin+s;
				Load[pixelRow, bitmap];
				FOR f: [0..1) IN [0..1) DO
					bit[s][f] _ pixelRow.sample[f+fOrigin-pixelRow.fOrigin];
					ENDLOOP;
				IF s+sOrigin NOT IN [swath.sMin..swath.sMin+swath.sSize)
				THEN fixed[s] _ ALL[1]
				ELSE {
					swathRow.sOrigin _ sOrigin+s;
					Load[swathRow, fixedBits];
					FOR f: [0..1) IN [0..1) DO
						fixed[s][f] _ swathRow.sample[f+fOrigin-swathRow.fOrigin];
						ENDLOOP;
					};
				ENDLOOP;
			pixelRow.sOrigin _ sOrigin;
			Load[pixelRow, gray];
			FOR f: INTEGER IN [0..1) DO
				target[f] _ pixelRow.sample[f+fOrigin-pixelRow.fOrigin];
				ENDLOOP;
			ComputeOuterModelIndices[];
			ComputePartialTrial[];
			ComputeFixedIndex[];
			END;
		
		Advance: PROC ~ BEGIN
			sOrigin _ sOrigin + 1;
			FOR s: [0..1) IN [0..1-1) DO
				bit[s] _ bit[s+1];
				fixed[s] _ fixed[s+1];
				ENDLOOP;
			pixelRow.sOrigin _ sOrigin+1-1;
			Load[pixelRow, bitmap];
			FOR f: [0..1) IN [0..1) DO
				bit[1-1][f] _ pixelRow.sample[f+fOrigin-pixelRow.fOrigin];
				ENDLOOP;
			IF 1-1+sOrigin NOT IN [swath.sMin..swath.sMin+swath.sSize)
			THEN fixed[1-1] _ ALL[1]
			ELSE {
				swathRow.sOrigin _ sOrigin+1-1;
				Load[swathRow, fixedBits];
				FOR f: [0..1) IN [0..1) DO
					fixed[1-1][f] _ swathRow.sample[f+fOrigin-swathRow.fOrigin];
					ENDLOOP;
				};
			pixelRow.sOrigin _ sOrigin;
			Load[pixelRow, gray];
			FOR f: INTEGER IN [0..1) DO
				target[f] _ pixelRow.sample[f+fOrigin-pixelRow.fOrigin];
				ENDLOOP;
			ComputeOuterModelIndices[];
			ComputePartialTrial[];
			ComputeFixedIndex[];
			END;
		
		ComputeOuterModelIndices: PROC ~ BEGIN
			
			END;
		
		ComputeInnerModelIndices: PROC ~ BEGIN
			index[0][0] _ bit[0][0];
			
			END;
		
		ComputePartialTrial: PROC ~ BEGIN
			-- From the index values
			t: REF ModelDataItem _ NIL;
			trial _ ALL[0];
			 
			END;
		
		ComputeRowPenalty: PROC RETURNS [penalty: CARDINAL] ~ BEGIN
			-- From the index values
			t: REF ModelDataItem _ NIL;
			try: ARRAY [0..1) OF INTEGER _ trial;
			t _ modelData[index[0][0]];
			try[0] _ try[0] + t.a[0][0];
			penalty _ 0 + ABS[try[0]-target[0]] + modelData[index[0][0]].noisePenalty;
			END;
		
		SetCombination: PROC [k: CARDINAL] ~ BEGIN
			-- k is a packed version of the non-fixed bits
			-- Sets questionBits and contextBits
			FOR i: NAT IN [0..nonFixedBitCount) DO
				sf: SF _ nonFixedBitIndex[i];
				bit[sf.s][sf.f] _ k MOD 2;
				k _ k/2;
				ENDLOOP;
			questionBits _ bit[0][0];
			contextBits _ 0;
			
			IF k # 0 THEN ERROR;
			prevContextBits _ contextBits/2 + 
			questionBits*(1/2);
			END;
		
		UnpackBits: PROC [questionBits, contextBits: CARDINAL] ~ BEGIN
			-- Unpacks the args into the swath portion of the bit array
			k: CARDINAL _ contextBits+questionBits*1;
			bit[0][0] _ k MOD 2;
			k _ k / 2;
			
			END;
		
		DoDynamicProgramming: PROC ~ {
			-- If subproblemTable#NIL, subproblemTable.first gives the best ways to color the pixels up to sOrigin
			UNTIL sOrigin = swath.sMin+swath.sSize DO
				Initialize[];
				BEGIN
					nComb: NAT _ Basics.BITSHIFT[1, nonFixedBitCount];
					new: LIST OF REF PerSwatchTable _ Alloc[];
					pst: REF PerSwatchTable _ new.first;
					prevPST: REF PerSwatchTable _ NIL;
					IF subproblemTable # NIL THEN prevPST _ subproblemTable.first;
					FOR k: NAT IN [0..nComb) DO
						totBad: INT _ 0;
						SetCombination[k];
						ComputeInnerModelIndices[];
						IF prevPST # NIL THEN totBad _ MinErr[prevPST[prevContextBits]];
						totBad _ totBad+ComputeRowPenalty[];
						IF totBad < MinErr[pst[contextBits]] THEN {
							pst[contextBits] _ MakeSubproblemSolution[questionBits, totBad];
							};
						ENDLOOP;
					new.rest _ subproblemTable;
					subproblemTable _ new;
					END;
				sOrigin _ sOrigin + 1;
				ENDLOOP;
			};
		
		Traceback: PROC ~ {
			lowestPenalty: INT _ MinErr[subproblemTable.first[0]];
			contextBits _ 0;
			-- First look through the most recent subproblem to find the best entry.
			FOR i: NAT IN (0..1) DO
				err: INT _ MinErr[subproblemTable.first[i]];
				IF err < lowestPenalty THEN {
					lowestPenalty _ err;
					contextBits _ i;
					};
				ENDLOOP;
			sOrigin _ swath.sMin+swath.sSize;
			FOR t: LIST OF REF PerSwatchTable _ subproblemTable, t.rest WHILE sOrigin-0 > swath.sMin DO
				questionBits _ t.first[contextBits].setting;
				sOrigin _ sOrigin - 1;
				UnpackBits[questionBits, contextBits];
				swathRow.sOrigin _ sOrigin;
				FOR j: NAT IN [0..swathRow.sample.jSize) DO
					swathRow.sample[j] _ bit[0][j];
					ENDLOOP;
				StorePixelRow[swathRow, bitmap];
				contextBits _ contextBits/2 + 
				questionBits*(1/2);
				ENDLOOP;
			};
		
		DoDynamicProgramming[];
		Traceback[];
		newScratch _ subproblemTable;
		 
		};
	
	
	
	DynamicBits.RegisterImpl[[
			sKernelRadius: 0,
			fKernelRadius: 0,
			sModelRadius: 0,
			fModelRadius: 0,
			swathSize: 1,
			CreatePrinterModel: CreatePrinterModel,
			DoWithPrinterModel: DoWithPrinterModel,
			Convolve: Convolve,
			ApplyModel: ApplyModel,
			TuneSwath: TuneSwath
			]];
	
	END.