% a class which performs one-sided line integral and thus gives edge responses in the tracking process
% whereas FullStrip calculates the two-sided line integral and is used in the detection stage
% the key method from this class is spline_integral()

% Written by Yi-Qing Wang, 2016


classdef SubStrip<handle

	properties
		input_strip;	% a full strip
		start_col;	% which substrip

		original;	% input_strip is conv2(input_strip)
		mask_width;

		calculated;	% storage for recursion (no longer in use)

		num_offsets;	% maximal number of offsets per scale

		interp_scale;
		interp_unitL;	% interpolation length
		increment;

		max_scale;

		max_slope;

		variance_test;	% reduce edge jumping

		max_num_lps;
		sigma_squared;
	end

	methods

		% constructor
		function obj = SubStrip(input_strip, start_col, interp_scale, max_scale, max_slope, ...
					variance_test, original, mask_width, max_num_lps, estimated_sigma)

			if nargin >= 5

				assert(start_col + 2^(max_scale-1) <= size(input_strip, 2));
				assert(max_scale >= interp_scale);

				obj.input_strip = input_strip;
				obj.start_col = start_col;
				obj.num_offsets = zeros(1, interp_scale);

				obj.interp_unitL = 2^(interp_scale-1);
				obj.increment = ones(1, obj.interp_unitL)/obj.interp_unitL;

				obj.interp_scale = interp_scale;
				obj.max_scale = max_scale;

				obj.max_slope = max_slope;

				if nargin == 10
					obj.original = original;
					obj.mask_width = mask_width;
					obj.variance_test = variance_test;
					obj.max_num_lps = max_num_lps;
					obj.sigma_squared = estimated_sigma * estimated_sigma;
				else
					obj.variance_test = false;
				end

			end
		end


		% calculate the integrals for a whole scale
		function trapezoidal_integrals = sweep(obj, scale, controls, start_row, num_rowls)



			trapezoidal_integrals = csweep(scale, obj.interp_scale, obj.max_slope, ...
						start_row-1, num_rowls, obj.start_col-1, ...
						controls, obj.input_strip, 1);

		end

		% calculate spline edge response
		function maxret = spline_integral(obj, edge_sign, controls, start_row, num_rowls, free_memory)

			%%% the slow matlab routine (replaced by csweep)
			%%% integrals = sweep(obj, obj.interp_scale, controls, start_row, num_rowls);

			maxret = [-1, 0, 0];
			integrals = csweep(obj.interp_scale, obj.interp_scale, obj.max_slope, ...
						start_row-1, num_rowls, obj.start_col-1, ...
						controls, obj.input_strip, free_memory, obj.max_num_lps);

			num_controls = size(controls, 1);
			if num_controls > 0

				% complete with the endpoints
				integrals = integrals + half_endpoints(obj, controls, num_controls, start_row, num_rowls);

				% sintegrals is the signed integrals
				sintegrals = integrals * edge_sign;
				% it might be possible that edge switching is required for an extension, deny that
				num = min(sum(sintegrals > 0), 10);

				if num > 0

					if obj.variance_test

						maxret = varianceTest(obj, sintegrals, num, start_row, num_rowls, controls);

					else
						% finding the only maximum might be totally intensity inconsistent

						[maxval, maxid] = max(sintegrals);
						k = maxid - 1; % conform to the set index convention
						rowl_id = mod(k, num_rowls);
						contid = (k - rowl_id)/num_rowls + 1;
						row = rowl_id + start_row;
						maxret = [maxval, row, contid];
					end
				end

			end

		end

		function [vals, pixel_vals] = get_intensity(obj, rows, controls)

			num_controls = size(controls, 1);
			num_contopts = size(controls, 2);

			% controls exclusively interpreted at scale obj.interp_scale
			num_colids = (num_contopts-1)*obj.interp_unitL;
			colids = ones(num_controls, 1) * (obj.start_col : (obj.start_col+num_colids));

			valdiff = diff(controls, 1, 2);
			rowids = floor(cumsum(horzcat(rows, kron(valdiff, obj.increment)), 2));
			vals = reshape(obj.input_strip(sub2ind(size(obj.input_strip), rowids(:), colids(:))), size(rowids));

			if nargout == 2
				rowids = rowids + obj.mask_width;
				pixel_vals = reshape(obj.original(sub2ind(size(obj.original), rowids(:), colids(:))), size(rowids));
			end

		end

		function maxret = varianceTest(obj, sintegrals, topk, start_row, num_rowls, controls)

			% first choose the topk signed integrals
			[~,sortIndex] = sort(sintegrals, 'descend');
			idx = 1:topk;
			sindx = sortIndex(idx);

			% translate into coordinates
			rowl_idx = mod(sindx-1, num_rowls);
			contidx = (sindx-1 - rowl_idx)/num_rowls + 1;
			rowx = rowl_idx' + start_row;
			[vals, pixel_vals] = get_intensity(obj, rowx, controls(contidx, :));

			% a remedy to the toy example edge contrast inconsistent at low noise
			if size(controls, 2) > 2

			% subject edge contrast to consistency test
			falsePositive = 7e-7;
			threshold = log(size(vals, 1)/falsePositive)/(size(vals, 2)-1);
			threshold = (threshold + sqrt(threshold)) * 2;
			threshold = threshold * obj.sigma_squared;
			stats = var(vals, 0, 2) - obj.sigma_squared;

			idx = idx(threshold > stats);

			end

			maxret = [-1, 0, 0];

			if ~isempty(idx)

				% intensity test
				falsePositive = 1e-5;
				threshold = log(1/falsePositive)/(size(pixel_vals, 2)-1);
				threshold = (threshold + sqrt(threshold)) * 2;
				threshold = threshold * obj.sigma_squared;

				pstat = var(pixel_vals(idx, :), 0, 2) - obj.sigma_squared;
				% return the first color consistent
				for i = 1:numel(idx)
					if pstat(i) < threshold
						minid = idx(i);
						abs_vals = abs(sum(vals(minid, :), 2));
						maxret = [abs_vals, rowx(minid), contidx(minid)];
						break;
					end
				end
			end
		end

		% add back half of the intensities left out by the recursion
		function vals = half_endpoints(obj, controls, num_controls, start_row, num_rowls)

			% vectorize for accessing matrix entries

			% fshifts = repmat(start_row:start_row+num_rowls-1, 1, num_controls);
			% rconts = controls(:, end);
			% frrowids = kron(rconts', ones(1, num_rowls)) + fshifts;

			shifts = (ones(num_controls, 1) * (start_row:start_row+num_rowls-1))';
			rrowids = ones(num_rowls, 1) * controls(:, end)' + shifts;

			% assert(isequal(fshifts', shifts(:)));
			% assert(isequal(frrowids', rrowids(:)));

			rcolid = 2^(obj.max_scale-1) + obj.start_col;
			vals = 0.5*(obj.input_strip(shifts(:), obj.start_col) + obj.input_strip(rrowids(:), rcolid))';

		end
	end
end
