% a class which provides an efficient method spline_integral()
% for evaluating all the edge responses within a strip

% Written by Yi-Qing Wang, 2016

classdef FullStrip<handle

	properties
		max_scale;
		max_slope;
		interp_scale;

		% lstrip;
		% rstrip;
		midcol;

		max_num_lps;

		chi_test;

		lhalf;
		rhalf;
		fullstrip;
	end

	methods

		% constructor
		function obj = FullStrip(fullstrip, interp_scale, max_scale, max_slope, max_num_lps, chi_test)

			if nargin == 6

				assert(2^max_scale+1 == size(fullstrip, 2));

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

				obj.max_num_lps = max_num_lps;

				% divide the full strip into left and right half
				midcol = 2^(max_scale-1) + 1;
				lhalf = fliplr(fullstrip(:, 1:midcol));
				rhalf = fullstrip(:, midcol:end);

				obj.midcol = midcol;

				obj.chi_test = chi_test;

				obj.lhalf = lhalf;
				obj.rhalf = rhalf;
				obj.fullstrip = fullstrip;

			else
				error('6 input arguments expected for FullStrip');
			end
		end

		% calculate spline edge response at any scale
		function maxret = spline_integral(obj, scale, controls, start_row, num_rowls, threshold)

			% divide controls into two
			middle = (size(controls, 2) + 1)/2;

			lconts = fliplr(controls(:, 1:middle));
			rconts = controls(:, middle:end);

			% calculate the responses
			% arranged as contid * num_rowls + rowl_id

			lintegrals = csweep(min(scale, obj.interp_scale), obj.interp_scale, obj.max_slope, ...
						start_row-1, num_rowls, 0, lconts, obj.lhalf, 1, obj.max_num_lps);

			% lintegralss = sweep(obj.lstrip, min(scale, obj.interp_scale), lconts, start_row, num_rowls);
			% assert(sum(abs(lintegrals-lintegralss)) < 1e-2);
			% clear lintegralss;

			rintegrals = csweep(min(scale, obj.interp_scale), obj.interp_scale, obj.max_slope, ...
						start_row-1, num_rowls, 0, rconts, obj.rhalf, 1, obj.max_num_lps);

			% rintegralss = sweep(obj.rstrip, min(scale, obj.interp_scale), rconts, start_row, num_rowls);
			% assert(sum(abs(rintegrals-rintegralss)) < 1e-2);
			% clear rintegralss;

			integrals = lintegrals + rintegrals;

			% complete with the endpoints
			num_controls = size(controls, 1);
			integrals = integrals + half_endpoints(obj, scale, lconts, rconts, num_controls, start_row, num_rowls);

			maxret.maxval = max(abs(integrals));
			if maxret.maxval > threshold

				ids = find(abs(integrals) > threshold);

				rowl_ids = mod(ids -1, num_rowls);
				cont_ids = (ids - 1 - rowl_ids)/num_rowls + 1;
				rows = rowl_ids + start_row;
				pixel_vals = get_intensity(obj, rows', controls(cont_ids, :), scale);

				f_vals = sum(pixel_vals, 2);
				abs_vs = abs(f_vals);

				if sum(abs_vs > threshold) == 0
					maxret.maxval = -1;
					return;
				end

				% filter once again with accurate pixel values
				cont_ids = cont_ids(abs_vs > threshold);
				rows = rows(abs_vs > threshold);
				f_vals = f_vals(abs_vs > threshold);
				abs_vs = abs(f_vals);

				fired = [f_vals'; abs_vs'; rows; cont_ids]';
				vals = sortrows(fired, [3, 2]);

				% powerline debugging
				% vals = vals(vals(:, 3) > 1630, :);

				% soccer field debugging
				% vals = vals(vals(:, 3) > 840, :);

				% a robustness test: a valid curve ought to be detected at least min_duration times
				% beneficial
				min_duration = 2;
				torf = true(1, size(vals, 1));
				multiplicity = diff(vals(:, 3));
				% new curves
				pos = find(multiplicity > 2);
				if ~isempty(pos)
					if pos(1) == 1
					    torf(1) = false;
					end
					for k = 1:numel(pos)-1
					    if pos(k+1) - pos(k) <= min_duration
					        torf(pos(k)+1:pos(k+1)) = false;
					    end
					end
					vals = vals(torf, :);
				end

				num_fires = size(vals, 1);

				% choose the most significant response among those test splines having the same rows
				change = find(diff(vals(:, 3)) > 0);
				if isempty(change)
					% all the test splines go through the same point
					maxret.vals = vals(end, :);
				else
					change = [change; num_fires];
					vals = vals(change, :);
					num_fires = size(vals, 1);

					% group edge responses based on their signs and row positions
					% http://blogs.mathworks.com/steve/2010/09/07/almost-connected-component-labeling/

					minrow = vals(1, 3);
					maxrow = vals(num_fires, 3);
					binary = zeros(maxrow-minrow+1, 2);
					reals = zeros(size(binary));

					signs = (sign(vals(:, 1)) + 1)/2 + 1;
					row_pos = vals(:, 3) - minrow + 1;

					binary(sub2ind(size(binary), row_pos, signs)) = 1;
					reals(sub2ind(size(binary), row_pos, signs)) = vals(:, 2);

					neighborhood_size = 0;
					dilated = bwdist(binary) <= neighborhood_size;
					labels = labelmatrix(bwconncomp(dilated,4));
					labels(~binary) = 0;

					stats = regionprops(labels, reals, 'MaxIntensity');

					idx = arrayfun(@(x)find(vals(:, 2)==x,1), [stats.MaxIntensity]);
					maxret.vals = vals(idx, :);
				end

				% go on reducing false alarms
				% use the max (pixel responses) - min (pixel responses)
				% along the detected curves

				% a remedy to edge contrast inconsistency with high SNR
				if obj.chi_test && scale > 1
					contid = maxret.vals(:, 4);
					maxret = gap_test(obj, maxret, controls(contid, :), scale);
				end

			end
		end

		% get the pixel intensities along a specified curved
		function vals = get_intensity(obj, rows, controls, scale)

			num_controls = size(controls, 1);

			half_size = 2^(scale-1);
			colids = ones(num_controls, 1) * ((-half_size:half_size) + obj.midcol);

			rowdif = diff(controls, 1, 2);
			interp_length = 2^(min(scale, obj.interp_scale) - 1);
			increment = ones(1, interp_length)/interp_length;

			rowids = floor(cumsum(horzcat(controls(:, 1), kron(rowdif, increment)), 2));
			rowids = bsxfun(@plus, rows, rowids);

			vals = reshape(obj.fullstrip(sub2ind(size(obj.fullstrip), rowids(:), colids(:))), size(rowids));

		end

		% chi square test for rejecting broken edges
		function maxret = gap_test(obj, maxret, controls, scale)

			num_controls = size(controls, 1);

			vals = get_intensity(obj, maxret.vals(:, 3), controls, scale);

			stats = var(vals, 0, 2) - 1;

			% stat2 = (sum((vals-repmat(mean(vals, 2), 1, size(vals, 2))).^2, 2)/(2^scale)) - 1;
			% assert(isequal(stats, stat2));

			falsePositive = 1e-2;
			threshold = log(num_controls/falsePositive)/(2^scale);
			threshold = (threshold + sqrt(threshold)) * 2;

			maxret.vals = maxret.vals(stats < threshold, :);
			if isempty(maxret.vals)
				maxret.maxval = -1;
			end
		end


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

			colid = 2^(scale-1) + 1;
			colids = colid * ones(1, num_controls*num_rowls);

			pos = 2^(max(scale - obj.interp_scale, 0))+1;
			shifts = repmat(start_row:start_row+num_rowls-1, 1, num_controls);
			lcont = lconts(:, pos);  % not necessary (enable assert)
			% assert(pos == size(lconts, 2));
			lrowids = kron(lcont', ones(1, num_rowls)) + shifts;
			rcont = rconts(:, pos);  % not necessary (enable assert)
			% assert(pos == size(rconts, 2));
			rrowids = kron(rcont', ones(1, num_rowls)) + shifts;

			vals = 0.5*(obj.lhalf(sub2ind(size(obj.lhalf), lrowids, colids)) + ...
					obj.rhalf(sub2ind(size(obj.rhalf), rrowids, colids)));

		end

	end


end
