Finally got blackman-harris to work correctly.

src/gui/CMakeLists.txt

@@ -6,6 +6,7 @@ add_executable(${PROJECT_NAME}
 	block.cpp
 	icon_fetch.cpp
 	grid.cpp
+	resize_harris.cpp
 )
 
 target_link_libraries(${PROJECT_NAME}

src/gui/grid.cpp

@@ -4,6 +4,10 @@
 #include "icon_fetch.hpp"
 #include <cassert>
 
+namespace {
+	const constexpr int g_icon_width = 16 * 5;
+} //unnamed namespace
+
 namespace implem {
 	void draw_icons (
 		int icon_size,
@@ -52,8 +56,8 @@ void grid_drawer_trigger::attached (
 	tabstop(wd);
 
 	for (int z = 0; z < static_cast<int>(m_icons.size()); ++z) {
-		const int x = (16 + 2) * (z % 3);
-		const int y = (16 + 2) * (z / 3);
+		const int x = (g_icon_width + 2) * (z % 3);
+		const int y = (g_icon_width + 2) * (z / 3);
 		m_icons[z].output(*m_widget, nana::point(x, y));
 	}
 }
@@ -89,6 +93,6 @@ void grid_drawer_trigger::refresh (nana::paint::graphics& graph) {
 
 	if (m_memorycard and not m_icons_added) {
 		m_icons_added = true;
-		implem::draw_icons(16, *m_memorycard, m_icons, *m_widget);
+		implem::draw_icons(g_icon_width, *m_memorycard, m_icons, *m_widget);
 	}
 }

src/gui/icon_fetch.cpp

@@ -1,4 +1,5 @@
 #include "icon_fetch.hpp"
+#include "resize_harris.hpp"
 #include <cstdint>
 #include <cassert>
 #include <climits>
@@ -63,10 +64,6 @@ namespace {
 		std::transform(data.begin(), data.end(), data.begin(), &to_srgb);
 	}
 
-	std::vector<float> bilinear_resize (int in_w, int in_h, const std::vector<float>& data, int out_w, int out_h) {
-		return data;
-	}
-
 	std::vector<float> rgb4_to_float (const std::vector<char>& data, const std::vector<uint8_t>& palette, int w, int h) {
 		assert(w * h);
 		std::vector<float> retval(w * h * 3);
@@ -105,11 +102,12 @@ namespace {
 		const int scanl_sz = scanline_size(w, 24);
 		std::vector<char> retval(scanl_sz * h);
 		assert(retval.size() >= data.size());
+		assert(data.size() == w * h * 3);
 
 		for (int y = 0; y < h; ++y) {
 			for (int x = 0; x < w; ++x) {
 				for (int c = 0; c < 3; ++c) {
-					retval[y * scanl_sz + x * 3 + c] = static_cast<char>(static_cast<uint8_t>(data[y * scanl_sz + x * 3 + c] * 255.0f));
+					retval[y * scanl_sz + x * 3 + c] = static_cast<char>(static_cast<uint8_t>(data[y * w * 3 + x * 3 + c] * 255.0f));
 				}
 			}
 		}
@@ -172,11 +170,10 @@ std::vector<std::vector<char>> icon_fetch (const ConstBlock& block, int width, i
 		std::vector<char> resized_rgb;
 		if (scale) {
 			std::vector<float> float_data = rgb4_to_float(orig_rgb, palette, in_width, in_height);
-			//srgb_to_linear(float_data, in_width, in_height);
-			//auto float_resized_rgb = bilinear_resize(in_width, in_height, float_data, width, height);
-			//linear_to_srgb(float_resized_rgb, width, height);
-			//resized_rgb = float_to_rgb24(float_resized_rgb, width, height);
-			resized_rgb = float_to_rgb24(float_data, width, height);
+			srgb_to_linear(float_data, in_width, in_height);
+			auto float_resized_rgb = blackmanharris_resize(in_width, in_height, float_data, width, height);
+			linear_to_srgb(float_resized_rgb, width, height);
+			resized_rgb = float_to_rgb24(float_resized_rgb, width, height);
 		}
 		else {
 			std::swap(resized_rgb, orig_rgb);
@@ -185,9 +182,9 @@ std::vector<std::vector<char>> icon_fetch (const ConstBlock& block, int width, i
 
 		assert(frame.size() == data_offset);
 		frame.resize(data_offset + scan_sz * height);
-		for (int32_t y = 0; y < in_height; ++y) {
+		for (int32_t y = 0; y < height; ++y) {
 			std::copy_n(
-				resized_rgb.begin() + scan_sz * (in_height - 1 - y),
+				resized_rgb.begin() + scan_sz * (height - 1 - y),
 				scan_sz,
 				frame.begin() + data_offset + y * scan_sz
 			);

src/gui/resize_harris.cpp

@@ -0,0 +1,309 @@
+/*
+  Copyright 2008 Larry Gritz and the other authors and contributors.
+  All Rights Reserved.
+  Based on BSD-licensed software Copyright 2004 NVIDIA Corp.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+  * Redistributions of source code must retain the above copyright
+	notice, this list of conditions and the following disclaimer.
+  * Redistributions in binary form must reproduce the above copyright
+	notice, this list of conditions and the following disclaimer in the
+	documentation and/or other materials provided with the distribution.
+  * Neither the name of the software's owners nor the names of its
+	contributors may be used to endorse or promote products derived from
+	this software without specific prior written permission.
+  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+  OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+  LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+  DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+  THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+  OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+  (This is the Modified BSD License)
+*/
+
+#include "resize_harris.hpp"
+#include <array>
+#include <cmath>
+#include <cassert>
+#include <algorithm>
+#include <numeric>
+
+namespace {
+	template <typename T>
+	[[gnu::pure]]
+	T clamp (T v, T m, T M) {
+		if (v < m)
+			return m;
+		else if (v > M)
+			return M;
+		else
+			return v;
+	}
+
+	class ClampedCoords {
+	public:
+		ClampedCoords (int tw, int th) :
+			m_xpos(0), m_ypos(0),
+			m_xbegin(0),
+			m_xend(0),
+			m_total_width(tw), m_total_height(th)
+		{
+		}
+
+		void reset (int xbegin, int xend, int ybegin) {
+			m_xpos = m_xbegin = xbegin;
+			m_ypos = ybegin;
+			m_xend = xend;
+		}
+
+		ClampedCoords& operator++() {
+			if (++m_xpos < m_xend) {
+				// Special case: we only incremented x, didn't change y
+				// or z, and the previous position was within the data
+				// window.  Call a shortcut version of pos.
+			}
+			else {
+				// Wrap to the next scanline
+				m_xpos = m_xbegin;
+				++m_ypos;
+			}
+		}
+
+		void advance_row() { ++m_ypos; }
+		int index (int nchannels) const { return xpos() * nchannels + ypos() * m_total_width * nchannels; }
+
+	private:
+		int xpos() const { return clamp(m_xpos, 0, m_total_width - 1); }
+		int ypos() const { return clamp(m_ypos, 0, m_total_height - 1); }
+
+		int m_xpos, m_ypos;
+		int m_xbegin;
+		int m_xend;
+		int m_total_width, m_total_height;
+	};
+
+	[[gnu::pure]]
+	/// Return (x-floor(x)) and put (int)floor(x) in *xint.  This is similar
+	/// to the built-in modf, but returns a true int, always rounds down
+	/// (compared to modf which rounds toward 0), and always returns
+	/// frac >= 0 (comapred to modf which can return <0 if x<0).
+	float floorfrac (float x, int *xint)
+	{
+#if 1
+		float f = std::floor(x);
+		*xint = int(f);
+		return x - f;
+#else /* vvv This idiom is slower */
+		int i = ifloor(x);
+		*xint = i;
+		return x - static_cast<float>(i);   // Return the fraction left over
+#endif
+	}
+
+	[[gnu::pure]]
+	int ifloor (float x) {
+		return static_cast<int>(std::floor(x));
+	}
+
+	[[gnu::pure]]
+	float blackman_harris_bh1d (float x) {
+		if (x < -1.0f || x > 1.0f)  // Early out if outside filter range
+			return 0.0f;
+		// Compute BH.  Straight from classic BH paper, but the usual
+		// formula assumes that the filter is centered at 0.5, so scale:
+		x = (x + 1.0f) * 0.5f;
+		const float A0 =  0.35875f;
+		const float A1 = -0.48829f;
+		const float A2 =  0.14128f;
+		const float A3 = -0.01168f;
+		const float m_pi(M_PI);
+#if 0
+		// original -- three cos calls!
+		return A0 + A1 * cosf(2.f * m_pi * x)
+			 + A2 * cosf(4.f * m_pi * x) + A3 * cosf(6.f * m_pi * x);
+#else
+		// Use trig identintities to reduce to just one cos.
+		// https://en.wikipedia.org/wiki/List_of_trigonometric_identities
+		// cos(2x) = 2 cos^2(x) - 1
+		// cos(3x) = 4 cos^3(x) − 3 cos(x)
+		float cos2pix = std::cos(2.f * m_pi * x);
+		float cos4pix = 2.0f * cos2pix * cos2pix - 1.0f;
+		float cos6pix = cos2pix * (2.0f * cos4pix - 1.0f);
+		return A0 + A1 * cos2pix + A2 * cos4pix + A3 * cos6pix;
+#endif
+	}
+
+	[[gnu::pure]]
+	float blackman_harris_xfilt (float x, float wrad_inv) {
+		return blackman_harris_bh1d(x * wrad_inv);
+	}
+
+	template <typename IT, typename T>
+	void normalize (IT beg, IT end, T div) {
+		while (beg != end) {
+			*beg /= div;
+			++beg;
+		}
+	}
+} //unnamed namespace
+
+std::vector<float> blackmanharris_resize (
+	int in_w,
+	int in_h,
+	const std::vector<float>& data,
+	int out_w,
+	int out_h
+) {
+	assert(in_w > 0 and in_h > 0);
+	assert(out_w > 0 and out_h > 0);
+	//if (in_w == out_w and in_h == out_h)
+	//	return data;
+
+	const constexpr int nchannels = 3;
+
+	// Local copies of the source image window, converted to float
+	const constexpr float srcfx = 0.0f;
+	const constexpr float srcfy = 0.0f;
+	const float srcfw = static_cast<float>(in_w);
+	const float srcfh = static_cast<float>(in_h);
+
+	// Ratios of dst/src size.  Values larger than 1 indicate that we
+	// are maximizing (enlarging the image), and thus want to smoothly
+	// interpolate.  Values less than 1 indicate that we are minimizing
+	// (shrinking the image), and thus want to properly filter out the
+	// high frequencies.
+	const float xratio = static_cast<float>(out_w) / srcfw; // 2 upsize, 0.5 downsize
+	const float yratio = static_cast<float>(out_h) / srcfh;
+
+	const float filter_width = 3.0f * std::max(1.0f, xratio);
+	const float& filter_height = filter_width;
+	const constexpr float dstfx = 0.0f;
+	const constexpr float dstfy = 0.0f;
+	const float dstfw = static_cast<float>(out_w);
+	const float dstfh = static_cast<float>(out_h);
+	const float dstpixelwidth = 1.0f / dstfw;
+	const float dstpixelheight = 1.0f / dstfh;
+	std::array<float, nchannels> pel;
+	const float filterrad = filter_width / 2.0f;
+
+	// radi,radj is the filter radius, as an integer, in source pixels.  We
+	// will filter the source over [x-radi, x+radi] X [y-radj,y+radj].
+	const int radi = static_cast<int>(std::ceil (filterrad/xratio));
+	const int radj = static_cast<int>(std::ceil (filterrad/yratio));
+	const int xtaps = 2*radi + 1;
+	const int ytaps = 2*radj + 1;
+	std::vector<float> yfiltval(ytaps);
+
+	// For separable filters, horizontal tap weights will be the same
+	// for every column. So we precompute all the tap weights for every
+	// x position we'll need. We do the same thing in y, but row by row
+	// inside the loop (since we never revisit a y row). This
+	// substantially speeds up resize.
+	std::vector<float> xfiltval_all(xtaps * out_w);
+	for (int x = 0;  x < out_w;  ++x) {
+		float * const xfiltval = xfiltval_all.data() + x * xtaps;
+		const float s = (x-dstfx+0.5f)*dstpixelwidth;
+		const float src_xf = srcfx + s * srcfw;
+		int src_x;
+		float src_xf_frac = floorfrac (src_xf, &src_x);
+		std::fill(pel.begin(), pel.end(), 0.0f);
+		for (int i = 0;  i < xtaps;  ++i) {
+			const float w = blackman_harris_xfilt (xratio * (i-radi-(src_xf_frac-0.5f)), 2.0f / filter_width);
+			xfiltval[i] = w;
+		}
+		const float totalweight_x = std::accumulate(xfiltval, xfiltval + xtaps, 0.0f);
+		if (totalweight_x != 0.0f) {
+			normalize(xfiltval, xfiltval + xtaps, totalweight_x);
+		}
+	}
+
+	// We're going to loop over all output pixels we're interested in.
+	//
+	// (s,t) = NDC space coordinates of the output sample we are computing.
+	//	 This is the "sample point".
+	// (src_xf, src_xf) = source pixel space float coordinates of the
+	//	 sample we're computing. We want to compute the weighted sum
+	//	 of all the source image pixels that fall under the filter when
+	//	 centered at that location.
+	// (src_x, src_y) = image space integer coordinates of the floor,
+	//	 i.e., the closest pixel in the source image.
+	// src_xf_frac and src_yf_frac are the position within that pixel
+	//	 of our sample.
+	//
+	// Separate cases for separable and non-separable filters.
+
+#if !defined(NDEBUG)
+	std::vector<float> retval(nchannels * out_w * out_h, 1.0f);
+#else
+	std::vector<float> retval(nchannels * out_w * out_h);
+#endif
+
+	auto out = retval.begin();
+	ClampedCoords srcpel(in_w, in_h);
+
+	for (int y = 0;  y < out_h;  ++y) {
+		const float t = (y-dstfy+0.5f)*dstpixelheight;
+		const float src_yf = srcfy + t * srcfh;
+		int src_y;
+		const float src_yf_frac = floorfrac (src_yf, &src_y);
+		// If using separable filters, our vertical set of filter tap
+		// weights will be the same for the whole scanline we're on.  Just
+		// compute and normalize them once.
+		for (int j = 0;  j < ytaps;  ++j) {
+			const float w = blackman_harris_xfilt(yratio * (j-radj-(src_yf_frac-0.5f)), 2.0f / filter_height);
+			yfiltval[j] = w;
+		}
+		const float totalweight_y = std::accumulate(yfiltval.begin(), yfiltval.end(), 0.0f);
+		if (totalweight_y != 0.0f) {
+			normalize(yfiltval.begin(), yfiltval.end(), totalweight_y);
+		}
+
+		for (int x = 0; x < out_w; ++x, out += nchannels) {
+			const float s = (x-dstfx+0.5f)*dstpixelwidth;
+			const float src_xf = srcfx + s * srcfw;
+			const int src_x = ifloor (src_xf);
+			std::fill_n(pel.begin(), pel.size(), 0.0f);
+			const float * const xfiltval = xfiltval_all.data() + x * xtaps;
+			const float totalweight_x = std::accumulate(xfiltval, xfiltval + xtaps, 0.0f);
+			if (totalweight_x != 0.0f) {
+				srcpel.reset(src_x-radi, src_x+radi+1, src_y-radj/*, src_y+radj+1*/);
+				for (int j = -radj;  j <= radj;  ++j) {
+					float wy = yfiltval[j+radj];
+					if (wy == 0.0f) {
+						// 0 weight for this y tap -- move to next line
+						//srcpel.pos (srcpel.x(), srcpel.y()+1, srcpel.z());
+						srcpel.advance_row();
+						continue;
+					}
+					for (int i = 0; i < xtaps; ++i, ++srcpel) {
+						const float w = wy * xfiltval[i];
+						if (w) {
+							const auto index = srcpel.index(nchannels);
+							for (int c = 0;  c < nchannels;  ++c) {
+								pel[c] += w * data[index + c];
+							}
+						}
+					}
+				}
+			}
+			// Copy the pixel value (already normalized) to the output.
+			if (totalweight_y == 0.0f) {
+				std::fill_n(out, nchannels, 0.0f);
+			}
+			else {
+				assert(pel.size() == nchannels);
+				std::copy_n(pel.begin(), pel.size(), out);
+			}
+		}
+	}
+
+	return retval;
+}

src/gui/resize_harris.hpp

@@ -0,0 +1,11 @@
+#pragma once
+
+#include <vector>
+
+std::vector<float> blackmanharris_resize (
+	int in_w,
+	int in_h,
+	const std::vector<float>& data,
+	int out_w,
+	int out_h
+);