-- 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;

DynamicBits202Impl: CEDAR PROGRAM
IMPORTS Basics, ImagerPixelMap, ImagerSample, PixelMapOps, DynamicBits
~ BEGIN OPEN DynamicBits;
	
	
	KernelArray: TYPE ~ ARRAY [(-2)..3) OF ARRAY [(-2)..3) 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 ~ [(-2), (-2), 5, 5];
	CheckKernel: PROC [bounds: DeviceRectangle] ~ {
		IF bounds # kernelBounds THEN ModelParameterFault[$kernel];
		};
	
	ArrayFromKernel: PROC [kernel: Kernel] RETURNS [kernelArray: KernelArray]  ~ {
		CheckKernel[kernel.Window];
		FOR s: INTEGER IN [(-2)..3) DO
			FOR f: INTEGER IN [(-2)..3) 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 [(-2)..3) DO
				FOR f: INTEGER IN [(-2)..3) 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 [(-2)..3) OF PixelRow _ CreatePixelRows[];
		CreatePixelRows: PROC RETURNS[p: ARRAY[(-2)..3) OF PixelRow]~{
			FOR i: INTEGER IN [(-2)..3) DO
				p[i] _ CreatePixelRow[bounds.sMin+i, bounds.fMin+(-2), bounds.fSize+5-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 [(-2)..3) DO
			FOR ff: INT IN [(-2)..3) DO
				short _ short AND kernelArray[ss][ff] >= 0;
				ENDLOOP;
			ENDLOOP;
		FOR i: INTEGER IN [(-2)..3) 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 [-(-2)..bounds.fSize-(-2)) DO
					pix: CARDINAL _ p[-(-2)].sample[f-(-2)]*Cd[kernelArray[(-2)][(-2)]] + p[-(-2)].sample[f-(-1)]*Cd[kernelArray[(-2)][(-1)]] + p[-(-2)].sample[f-0]*Cd[kernelArray[(-2)][0]] + p[-(-2)].sample[f-1]*Cd[kernelArray[(-2)][1]] + p[-(-2)].sample[f-2]*Cd[kernelArray[(-2)][2]] + p[-(-1)].sample[f-(-2)]*Cd[kernelArray[(-1)][(-2)]] + p[-(-1)].sample[f-(-1)]*Cd[kernelArray[(-1)][(-1)]] + p[-(-1)].sample[f-0]*Cd[kernelArray[(-1)][0]] + p[-(-1)].sample[f-1]*Cd[kernelArray[(-1)][1]] + p[-(-1)].sample[f-2]*Cd[kernelArray[(-1)][2]] + p[-0].sample[f-(-2)]*Cd[kernelArray[0][(-2)]] + p[-0].sample[f-(-1)]*Cd[kernelArray[0][(-1)]] + p[-0].sample[f-0]*Cd[kernelArray[0][0]] + p[-0].sample[f-1]*Cd[kernelArray[0][1]] + p[-0].sample[f-2]*Cd[kernelArray[0][2]] + p[-1].sample[f-(-2)]*Cd[kernelArray[1][(-2)]] + p[-1].sample[f-(-1)]*Cd[kernelArray[1][(-1)]] + p[-1].sample[f-0]*Cd[kernelArray[1][0]] + p[-1].sample[f-1]*Cd[kernelArray[1][1]] + p[-1].sample[f-2]*Cd[kernelArray[1][2]] + p[-2].sample[f-(-2)]*Cd[kernelArray[2][(-2)]] + p[-2].sample[f-(-1)]*Cd[kernelArray[2][(-1)]] + p[-2].sample[f-0]*Cd[kernelArray[2][0]] + p[-2].sample[f-1]*Cd[kernelArray[2][1]] + p[-2].sample[f-2]*Cd[kernelArray[2][2]] +  0;
					q.sample[f] _ (pix+halfUnit)/unitIntensity;
					ENDLOOP;
				}
			ELSE {
				FOR f: NAT IN [-(-2)..bounds.fSize-(-2)) DO
					pix: INT _ 0;
					FOR ss: INT IN [(-2)..3) DO
						FOR ff: INT IN [(-2)..3) 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[(-2)]; p[(-2)] _ p[(-1)]; p[(-1)] _ p[0]; p[0] _ p[1]; p[1] _ p[2]; p[3-1] _ t};
			p[3-1].sOrigin _ p[3-1].sOrigin + (3-(-2));
			Load[p[3-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-2;
		-- sOrigin and fOrigin give the position of the origin of the local swatch in the bitmap coordinate system.
		sOrigin: INTEGER _ swath.sMin-1-2;
		fOrigin: INTEGER _ swath.fMin;
		pixelRow: PixelRow _ CreatePixelRow[sOrigin, fOrigin+(-4), 6-(-4)];
		swathRow: PixelRow _ CreatePixelRow[sOrigin, fOrigin, 2];
		bit: ARRAY [(-2)..3) OF ARRAY [(-4)..6) OF [0..1];
		-- The current bits from the bitmap within the swatch
		index: ARRAY [(-2)..3) OF ARRAY [(-4)..6) OF [0..2);
		-- The current bits, packed up into indices into the model array
		fixed: ARRAY [(-2)..3) OF ARRAY [0..2) OF [0..1];
		-- A local copy of a piece of fixedBits
		target: ARRAY [(-2)..4) OF INTEGER;
		-- A local copy of a piece of gray
		trial: ARRAY [(-2)..4) OF INTEGER;
		-- The partially computed approximation to the gray
		modelData: ModelData _ NARROW[model.data];
		
		nonFixedBitCount: NAT _ 0;
		SF: TYPE ~ RECORD [s: [(-2)..3), f: [0..2)];
		nonFixedBitIndex: ARRAY [0..(3-(-2))*2) OF SF;
		
		questionBits: [0..4);
		-- The bits that need to be determined for a given context.
		contextBits: [0..256);
		-- The context bits in the swath.
		prevContextBits: [0..256);
		-- The context bits for the previous subproblem.
		
		SubproblemSolution: TYPE ~ RECORD [
			setting: [0..4),
			-- The setting of questionBits that achieves the lowest total penalty for the swath so far.
			minErrHighHalf: [0..CARDINAL.LAST/4],
			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..4), 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..256) 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/4, CARDINAL.LAST]];
			END;
		
		ComputeFixedIndex: PROC ~ BEGIN
			nonFixedBitCount _ 0;
			FOR s: [(-2)..3) IN [(-2)..3) DO
				FOR f: [0..2) IN [0..2) DO
					IF fixed[s][f] = 0 THEN {
						nonFixedBitIndex[nonFixedBitCount] _ [s,f];
						nonFixedBitCount _ nonFixedBitCount + 1;
						};
					ENDLOOP;
				ENDLOOP;
			END;
		
		Initialize: PROC ~ BEGIN
			FOR s: [(-2)..3) IN [(-2)..3) DO
				pixelRow.sOrigin _ sOrigin+s;
				Load[pixelRow, bitmap];
				FOR f: [(-4)..6) IN [(-4)..6) 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..2) IN [0..2) DO
						fixed[s][f] _ swathRow.sample[f+fOrigin-swathRow.fOrigin];
						ENDLOOP;
					};
				ENDLOOP;
			pixelRow.sOrigin _ sOrigin;
			Load[pixelRow, gray];
			FOR f: INTEGER IN [(-2)..4) DO
				target[f] _ pixelRow.sample[f+fOrigin-pixelRow.fOrigin];
				ENDLOOP;
			ComputeOuterModelIndices[];
			ComputePartialTrial[];
			ComputeFixedIndex[];
			END;
		
		Advance: PROC ~ BEGIN
			sOrigin _ sOrigin + 1;
			FOR s: [(-2)..3) IN [(-2)..3-1) DO
				bit[s] _ bit[s+1];
				fixed[s] _ fixed[s+1];
				ENDLOOP;
			pixelRow.sOrigin _ sOrigin+3-1;
			Load[pixelRow, bitmap];
			FOR f: [(-4)..6) IN [(-4)..6) DO
				bit[3-1][f] _ pixelRow.sample[f+fOrigin-pixelRow.fOrigin];
				ENDLOOP;
			IF 3-1+sOrigin NOT IN [swath.sMin..swath.sMin+swath.sSize)
			THEN fixed[3-1] _ ALL[1]
			ELSE {
				swathRow.sOrigin _ sOrigin+3-1;
				Load[swathRow, fixedBits];
				FOR f: [0..2) IN [0..2) DO
					fixed[3-1][f] _ swathRow.sample[f+fOrigin-swathRow.fOrigin];
					ENDLOOP;
				};
			pixelRow.sOrigin _ sOrigin;
			Load[pixelRow, gray];
			FOR f: INTEGER IN [(-2)..4) DO
				target[f] _ pixelRow.sample[f+fOrigin-pixelRow.fOrigin];
				ENDLOOP;
			ComputeOuterModelIndices[];
			ComputePartialTrial[];
			ComputeFixedIndex[];
			END;
		
		ComputeOuterModelIndices: PROC ~ BEGIN
			index[(-2)][(-4)] _ bit[(-2)][(-4)];
			index[(-2)][(-3)] _ bit[(-2)][(-3)];
			index[(-2)][(-2)] _ bit[(-2)][(-2)];
			index[(-2)][(-1)] _ bit[(-2)][(-1)];
			index[(-2)][2] _ bit[(-2)][2];
			index[(-2)][3] _ bit[(-2)][3];
			index[(-2)][4] _ bit[(-2)][4];
			index[(-2)][5] _ bit[(-2)][5];
			index[(-1)][(-4)] _ bit[(-1)][(-4)];
			index[(-1)][(-3)] _ bit[(-1)][(-3)];
			index[(-1)][(-2)] _ bit[(-1)][(-2)];
			index[(-1)][(-1)] _ bit[(-1)][(-1)];
			index[(-1)][2] _ bit[(-1)][2];
			index[(-1)][3] _ bit[(-1)][3];
			index[(-1)][4] _ bit[(-1)][4];
			index[(-1)][5] _ bit[(-1)][5];
			index[0][(-4)] _ bit[0][(-4)];
			index[0][(-3)] _ bit[0][(-3)];
			index[0][(-2)] _ bit[0][(-2)];
			index[0][(-1)] _ bit[0][(-1)];
			index[0][2] _ bit[0][2];
			index[0][3] _ bit[0][3];
			index[0][4] _ bit[0][4];
			index[0][5] _ bit[0][5];
			index[1][(-4)] _ bit[1][(-4)];
			index[1][(-3)] _ bit[1][(-3)];
			index[1][(-2)] _ bit[1][(-2)];
			index[1][(-1)] _ bit[1][(-1)];
			index[1][2] _ bit[1][2];
			index[1][3] _ bit[1][3];
			index[1][4] _ bit[1][4];
			index[1][5] _ bit[1][5];
			index[2][(-4)] _ bit[2][(-4)];
			index[2][(-3)] _ bit[2][(-3)];
			index[2][(-2)] _ bit[2][(-2)];
			index[2][(-1)] _ bit[2][(-1)];
			index[2][2] _ bit[2][2];
			index[2][3] _ bit[2][3];
			index[2][4] _ bit[2][4];
			index[2][5] _ bit[2][5];
			
			END;
		
		ComputeInnerModelIndices: PROC ~ BEGIN
			index[(-2)][0] _ bit[(-2)][0];
			index[(-2)][1] _ bit[(-2)][1];
			index[(-1)][0] _ bit[(-1)][0];
			index[(-1)][1] _ bit[(-1)][1];
			index[0][0] _ bit[0][0];
			index[0][1] _ bit[0][1];
			index[1][0] _ bit[1][0];
			index[1][1] _ bit[1][1];
			index[2][0] _ bit[2][0];
			index[2][1] _ bit[2][1];
			
			END;
		
		ComputePartialTrial: PROC ~ BEGIN
			-- From the index values
			t: REF ModelDataItem _ NIL;
			trial _ ALL[0];
			t _ modelData[index[(-2)][(-4)]];
			trial[(-2)] _ trial[(-2)] + t.a[2][2];
			t _ modelData[index[(-2)][(-3)]];
			trial[(-2)] _ trial[(-2)] + t.a[2][1];
			trial[(-1)] _ trial[(-1)] + t.a[2][2];
			t _ modelData[index[(-2)][(-2)]];
			trial[(-2)] _ trial[(-2)] + t.a[2][0];
			trial[(-1)] _ trial[(-1)] + t.a[2][1];
			trial[0] _ trial[0] + t.a[2][2];
			t _ modelData[index[(-2)][(-1)]];
			trial[(-2)] _ trial[(-2)] + t.a[2][(-1)];
			trial[(-1)] _ trial[(-1)] + t.a[2][0];
			trial[0] _ trial[0] + t.a[2][1];
			trial[1] _ trial[1] + t.a[2][2];
			t _ modelData[index[(-2)][2]];
			trial[0] _ trial[0] + t.a[2][(-2)];
			trial[1] _ trial[1] + t.a[2][(-1)];
			trial[2] _ trial[2] + t.a[2][0];
			trial[3] _ trial[3] + t.a[2][1];
			t _ modelData[index[(-2)][3]];
			trial[1] _ trial[1] + t.a[2][(-2)];
			trial[2] _ trial[2] + t.a[2][(-1)];
			trial[3] _ trial[3] + t.a[2][0];
			t _ modelData[index[(-2)][4]];
			trial[2] _ trial[2] + t.a[2][(-2)];
			trial[3] _ trial[3] + t.a[2][(-1)];
			t _ modelData[index[(-2)][5]];
			trial[3] _ trial[3] + t.a[2][(-2)];
			t _ modelData[index[(-1)][(-4)]];
			trial[(-2)] _ trial[(-2)] + t.a[1][2];
			t _ modelData[index[(-1)][(-3)]];
			trial[(-2)] _ trial[(-2)] + t.a[1][1];
			trial[(-1)] _ trial[(-1)] + t.a[1][2];
			t _ modelData[index[(-1)][(-2)]];
			trial[(-2)] _ trial[(-2)] + t.a[1][0];
			trial[(-1)] _ trial[(-1)] + t.a[1][1];
			trial[0] _ trial[0] + t.a[1][2];
			t _ modelData[index[(-1)][(-1)]];
			trial[(-2)] _ trial[(-2)] + t.a[1][(-1)];
			trial[(-1)] _ trial[(-1)] + t.a[1][0];
			trial[0] _ trial[0] + t.a[1][1];
			trial[1] _ trial[1] + t.a[1][2];
			t _ modelData[index[(-1)][2]];
			trial[0] _ trial[0] + t.a[1][(-2)];
			trial[1] _ trial[1] + t.a[1][(-1)];
			trial[2] _ trial[2] + t.a[1][0];
			trial[3] _ trial[3] + t.a[1][1];
			t _ modelData[index[(-1)][3]];
			trial[1] _ trial[1] + t.a[1][(-2)];
			trial[2] _ trial[2] + t.a[1][(-1)];
			trial[3] _ trial[3] + t.a[1][0];
			t _ modelData[index[(-1)][4]];
			trial[2] _ trial[2] + t.a[1][(-2)];
			trial[3] _ trial[3] + t.a[1][(-1)];
			t _ modelData[index[(-1)][5]];
			trial[3] _ trial[3] + t.a[1][(-2)];
			t _ modelData[index[0][(-4)]];
			trial[(-2)] _ trial[(-2)] + t.a[0][2];
			t _ modelData[index[0][(-3)]];
			trial[(-2)] _ trial[(-2)] + t.a[0][1];
			trial[(-1)] _ trial[(-1)] + t.a[0][2];
			t _ modelData[index[0][(-2)]];
			trial[(-2)] _ trial[(-2)] + t.a[0][0];
			trial[(-1)] _ trial[(-1)] + t.a[0][1];
			trial[0] _ trial[0] + t.a[0][2];
			t _ modelData[index[0][(-1)]];
			trial[(-2)] _ trial[(-2)] + t.a[0][(-1)];
			trial[(-1)] _ trial[(-1)] + t.a[0][0];
			trial[0] _ trial[0] + t.a[0][1];
			trial[1] _ trial[1] + t.a[0][2];
			t _ modelData[index[0][2]];
			trial[0] _ trial[0] + t.a[0][(-2)];
			trial[1] _ trial[1] + t.a[0][(-1)];
			trial[2] _ trial[2] + t.a[0][0];
			trial[3] _ trial[3] + t.a[0][1];
			t _ modelData[index[0][3]];
			trial[1] _ trial[1] + t.a[0][(-2)];
			trial[2] _ trial[2] + t.a[0][(-1)];
			trial[3] _ trial[3] + t.a[0][0];
			t _ modelData[index[0][4]];
			trial[2] _ trial[2] + t.a[0][(-2)];
			trial[3] _ trial[3] + t.a[0][(-1)];
			t _ modelData[index[0][5]];
			trial[3] _ trial[3] + t.a[0][(-2)];
			t _ modelData[index[1][(-4)]];
			trial[(-2)] _ trial[(-2)] + t.a[(-1)][2];
			t _ modelData[index[1][(-3)]];
			trial[(-2)] _ trial[(-2)] + t.a[(-1)][1];
			trial[(-1)] _ trial[(-1)] + t.a[(-1)][2];
			t _ modelData[index[1][(-2)]];
			trial[(-2)] _ trial[(-2)] + t.a[(-1)][0];
			trial[(-1)] _ trial[(-1)] + t.a[(-1)][1];
			trial[0] _ trial[0] + t.a[(-1)][2];
			t _ modelData[index[1][(-1)]];
			trial[(-2)] _ trial[(-2)] + t.a[(-1)][(-1)];
			trial[(-1)] _ trial[(-1)] + t.a[(-1)][0];
			trial[0] _ trial[0] + t.a[(-1)][1];
			trial[1] _ trial[1] + t.a[(-1)][2];
			t _ modelData[index[1][2]];
			trial[0] _ trial[0] + t.a[(-1)][(-2)];
			trial[1] _ trial[1] + t.a[(-1)][(-1)];
			trial[2] _ trial[2] + t.a[(-1)][0];
			trial[3] _ trial[3] + t.a[(-1)][1];
			t _ modelData[index[1][3]];
			trial[1] _ trial[1] + t.a[(-1)][(-2)];
			trial[2] _ trial[2] + t.a[(-1)][(-1)];
			trial[3] _ trial[3] + t.a[(-1)][0];
			t _ modelData[index[1][4]];
			trial[2] _ trial[2] + t.a[(-1)][(-2)];
			trial[3] _ trial[3] + t.a[(-1)][(-1)];
			t _ modelData[index[1][5]];
			trial[3] _ trial[3] + t.a[(-1)][(-2)];
			t _ modelData[index[2][(-4)]];
			trial[(-2)] _ trial[(-2)] + t.a[(-2)][2];
			t _ modelData[index[2][(-3)]];
			trial[(-2)] _ trial[(-2)] + t.a[(-2)][1];
			trial[(-1)] _ trial[(-1)] + t.a[(-2)][2];
			t _ modelData[index[2][(-2)]];
			trial[(-2)] _ trial[(-2)] + t.a[(-2)][0];
			trial[(-1)] _ trial[(-1)] + t.a[(-2)][1];
			trial[0] _ trial[0] + t.a[(-2)][2];
			t _ modelData[index[2][(-1)]];
			trial[(-2)] _ trial[(-2)] + t.a[(-2)][(-1)];
			trial[(-1)] _ trial[(-1)] + t.a[(-2)][0];
			trial[0] _ trial[0] + t.a[(-2)][1];
			trial[1] _ trial[1] + t.a[(-2)][2];
			t _ modelData[index[2][2]];
			trial[0] _ trial[0] + t.a[(-2)][(-2)];
			trial[1] _ trial[1] + t.a[(-2)][(-1)];
			trial[2] _ trial[2] + t.a[(-2)][0];
			trial[3] _ trial[3] + t.a[(-2)][1];
			t _ modelData[index[2][3]];
			trial[1] _ trial[1] + t.a[(-2)][(-2)];
			trial[2] _ trial[2] + t.a[(-2)][(-1)];
			trial[3] _ trial[3] + t.a[(-2)][0];
			t _ modelData[index[2][4]];
			trial[2] _ trial[2] + t.a[(-2)][(-2)];
			trial[3] _ trial[3] + t.a[(-2)][(-1)];
			t _ modelData[index[2][5]];
			trial[3] _ trial[3] + t.a[(-2)][(-2)];
			 
			END;
		
		ComputeRowPenalty: PROC RETURNS [penalty: CARDINAL] ~ BEGIN
			-- From the index values
			t: REF ModelDataItem _ NIL;
			try: ARRAY [(-2)..4) OF INTEGER _ trial;
			t _ modelData[index[(-2)][0]];
			try[(-2)] _ try[(-2)] + t.a[2][(-2)];
			try[(-1)] _ try[(-1)] + t.a[2][(-1)];
			try[0] _ try[0] + t.a[2][0];
			try[1] _ try[1] + t.a[2][1];
			try[2] _ try[2] + t.a[2][2];
			t _ modelData[index[(-2)][1]];
			try[(-1)] _ try[(-1)] + t.a[2][(-2)];
			try[0] _ try[0] + t.a[2][(-1)];
			try[1] _ try[1] + t.a[2][0];
			try[2] _ try[2] + t.a[2][1];
			try[3] _ try[3] + t.a[2][2];
			t _ modelData[index[(-1)][0]];
			try[(-2)] _ try[(-2)] + t.a[1][(-2)];
			try[(-1)] _ try[(-1)] + t.a[1][(-1)];
			try[0] _ try[0] + t.a[1][0];
			try[1] _ try[1] + t.a[1][1];
			try[2] _ try[2] + t.a[1][2];
			t _ modelData[index[(-1)][1]];
			try[(-1)] _ try[(-1)] + t.a[1][(-2)];
			try[0] _ try[0] + t.a[1][(-1)];
			try[1] _ try[1] + t.a[1][0];
			try[2] _ try[2] + t.a[1][1];
			try[3] _ try[3] + t.a[1][2];
			t _ modelData[index[0][0]];
			try[(-2)] _ try[(-2)] + t.a[0][(-2)];
			try[(-1)] _ try[(-1)] + t.a[0][(-1)];
			try[0] _ try[0] + t.a[0][0];
			try[1] _ try[1] + t.a[0][1];
			try[2] _ try[2] + t.a[0][2];
			t _ modelData[index[0][1]];
			try[(-1)] _ try[(-1)] + t.a[0][(-2)];
			try[0] _ try[0] + t.a[0][(-1)];
			try[1] _ try[1] + t.a[0][0];
			try[2] _ try[2] + t.a[0][1];
			try[3] _ try[3] + t.a[0][2];
			t _ modelData[index[1][0]];
			try[(-2)] _ try[(-2)] + t.a[(-1)][(-2)];
			try[(-1)] _ try[(-1)] + t.a[(-1)][(-1)];
			try[0] _ try[0] + t.a[(-1)][0];
			try[1] _ try[1] + t.a[(-1)][1];
			try[2] _ try[2] + t.a[(-1)][2];
			t _ modelData[index[1][1]];
			try[(-1)] _ try[(-1)] + t.a[(-1)][(-2)];
			try[0] _ try[0] + t.a[(-1)][(-1)];
			try[1] _ try[1] + t.a[(-1)][0];
			try[2] _ try[2] + t.a[(-1)][1];
			try[3] _ try[3] + t.a[(-1)][2];
			t _ modelData[index[2][0]];
			try[(-2)] _ try[(-2)] + t.a[(-2)][(-2)];
			try[(-1)] _ try[(-1)] + t.a[(-2)][(-1)];
			try[0] _ try[0] + t.a[(-2)][0];
			try[1] _ try[1] + t.a[(-2)][1];
			try[2] _ try[2] + t.a[(-2)][2];
			t _ modelData[index[2][1]];
			try[(-1)] _ try[(-1)] + t.a[(-2)][(-2)];
			try[0] _ try[0] + t.a[(-2)][(-1)];
			try[1] _ try[1] + t.a[(-2)][0];
			try[2] _ try[2] + t.a[(-2)][1];
			try[3] _ try[3] + t.a[(-2)][2];
			penalty _ 0 + ABS[try[(-2)]-target[(-2)]] + modelData[index[0][(-2)]].noisePenalty + ABS[try[(-1)]-target[(-1)]] + modelData[index[0][(-1)]].noisePenalty + ABS[try[0]-target[0]] + modelData[index[0][0]].noisePenalty + ABS[try[1]-target[1]] + modelData[index[0][1]].noisePenalty + ABS[try[2]-target[2]] + modelData[index[0][2]].noisePenalty + ABS[try[3]-target[3]] + modelData[index[0][3]].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[(-2)][0]*2 + bit[(-2)][1]);
			contextBits _ (((((((bit[(-1)][0]*2 + bit[(-1)][1])*2 + bit[0][0])*2 + bit[0][1])*2 + bit[1][0])*2 + bit[1][1])*2 + bit[2][0])*2 + bit[2][1]);
			
			IF k # 0 THEN ERROR;
			prevContextBits _ contextBits/4 + 
			questionBits*(256/4);
			END;
		
		UnpackBits: PROC [questionBits, contextBits: CARDINAL] ~ BEGIN
			-- Unpacks the args into the swath portion of the bit array
			k: CARDINAL _ contextBits+questionBits*256;
			bit[2][1] _ k MOD 2;
			k _ k / 2;
			bit[2][0] _ k MOD 2;
			k _ k / 2;
			bit[1][1] _ k MOD 2;
			k _ k / 2;
			bit[1][0] _ k MOD 2;
			k _ k / 2;
			bit[0][1] _ k MOD 2;
			k _ k / 2;
			bit[0][0] _ k MOD 2;
			k _ k / 2;
			bit[(-1)][1] _ k MOD 2;
			k _ k / 2;
			bit[(-1)][0] _ k MOD 2;
			k _ k / 2;
			bit[(-2)][1] _ k MOD 2;
			k _ k / 2;
			bit[(-2)][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..256) 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-2 > 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/4 + 
				questionBits*(256/4);
				ENDLOOP;
			};
		
		DoDynamicProgramming[];
		Traceback[];
		newScratch _ subproblemTable;
		 
		};
	
	
	
	DynamicBits.RegisterImpl[[
			sKernelRadius: 2,
			fKernelRadius: 2,
			sModelRadius: 0,
			fModelRadius: 0,
			swathSize: 2,
			CreatePrinterModel: CreatePrinterModel,
			DoWithPrinterModel: DoWithPrinterModel,
			Convolve: Convolve,
			ApplyModel: ApplyModel,
			TuneSwath: TuneSwath
			]];
	
	END.