Commit: 4986e10d08106c1adcd49e5dc918735fa3ca7736 - color_quantizer (git)

減色プログラム

Commit MetaInfo

• Revisión: 4986e10d08106c1adcd49e5dc918735fa3ca7736 (tree)
• Tiempo: 2011-05-22 07:22:26
• Autor: berupon <berupon@gmai...>
• Commiter: berupon

Log Message

optimized

Cambiar Resumen

• modified: Color4d.h
• modified: Color4f.h
• modified: quantize.cpp

Diferencia incremental

--- a/Color4d.h
+++ b/Color4d.h
@@ -38,11 +38,18 @@ struct Color4d
 
 	double dot_product(const Color4d& rhs) {
 // http://www.icnet.ne.jp/~nsystem/simd_tobira/dpps.html
+#if 1
+		Color4d result = (*this) * rhs;
+		__m128d v = _mm_add_pd(result.v[0], result.v[1]);
+		v = _mm_hadd_pd(v, v);
+		return v.m128d_f64[0];
+#else
 		double result = 0;
 		for (int i=0; i<3; i++) {
 			result += (*this)[i] * rhs[i];
 		}
 		return result;
+#endif
 	}
 
 	Color4d& operator += (const Color4d& rhs) {
@@ -51,16 +58,16 @@ struct Color4d
 		return *this;
 	}
 
+/*
 	Color4d& operator += (const Color4f& rhs) {
 		v[0] = _mm_add_pd(v[0], _mm_cvtps_pd(rhs.v));
 		v[1] = _mm_add_pd(v[1], _mm_cvtps_pd(_mm_movehl_ps(rhs.v, rhs.v)));
 		return *this;
 	}
+*/
 
 	Color4d operator + (const Color4d& rhs) {
-		Color4d result(*this);
-		result += rhs;
-		return result;
+		return Color4d(*this) += rhs;
 	}
 
 	Color4d& operator -= (const Color4d& rhs) {
@@ -70,9 +77,7 @@ struct Color4d
 	}
 
 	Color4d operator - (const Color4d& rhs) {
-		Color4d result(*this);
-		result -= rhs;
-		return result;
+		return Color4d(*this) -= rhs;
 	}
 	
 	Color4d& operator *= (const Color4d& rhs) {
@@ -82,9 +87,7 @@ struct Color4d
 	}
 
 	Color4d operator * (const Color4d& rhs) {
-		Color4d result(*this);
-		result *= rhs;
-		return result;
+		return Color4d(*this) *= rhs;
 	}
 
 	Color4d& operator *= (double scalar) {
@@ -95,11 +98,9 @@ struct Color4d
 	}
 
 	Color4d operator * (double scalar) {
-		Color4d result(*this);
-		result *= scalar;
-		return result;
+		return Color4d(*this) *= scalar;
 	}
-
+	
 	double& operator[] (int idx) {
 		return ((double*)&v)[idx];
 	}
@@ -108,13 +109,22 @@ struct Color4d
 	}
 	
 	double norm_squared() {
+#if 1
+		__m128d t = _mm_add_pd(
+			_mm_mul_pd(v[0], v[0]),
+			_mm_mul_pd(v[1], v[1])
+			);
+		t = _mm_hadd_pd(t, t);
+		return t.m128d_f64[0];
+#else
 		double result = 0;
 		for (int i=0; i<3; i++) {
 			result += (*this)[i] * (*this)[i];
 		}
 		return result;
+#endif
 	}
-
+	
 	void zero() {
 		v[0] = _mm_setzero_pd();
 		v[1] = _mm_setzero_pd();
@@ -122,7 +132,10 @@ struct Color4d
 };
 
 inline Color4d operator * (double scalar, const Color4d& c) {
-	Color4d tmp = c;
-	return tmp * scalar;
+	return Color4d(c) *= scalar;
+}
+
+inline Color4d operator * (const Color4d& c, double scalar) {
+	return Color4d(c) *= scalar;
 }
 

--- a/Color4f.h
+++ b/Color4f.h
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <emmintrin.h>
+#include <intrin.h>
 
 // http://www2.kobe-u.ac.jp/~lerl2/l_cc_p_10.1.008/doc/main_cls/mergedProjects/intref_cls/
 

--- a/quantize.cpp
+++ b/quantize.cpp
@@ -293,18 +293,17 @@ void compute_initial_s(
 	)
 {
 	size_t palette_size  = s.width_;
-	int coarse_width  = coarse_variables.width_;
-	int coarse_height = coarse_variables.height_;
-	int center_x = (b.width_-1)/2, center_y = (b.height_-1)/2;
-	
-	Color center_b = b_value(b,0,0,0,0);
 	Color zero_vector;
 	zero_vector.zero();
 	for (size_t v=0; v<palette_size; ++v) {
-		for (size_t alpha=v; alpha<palette_size; ++alpha) {
-			s[alpha][v] = zero_vector;
+		for (size_t v2=0; v2<v+1; ++v2) {
+			s[v][v2] = zero_vector;
 		}
 	}
+	const int coarse_width = coarse_variables.width_;
+	const int coarse_height = coarse_variables.height_;
+	const int center_x = (b.width_-1)/2, center_y = (b.height_-1)/2;
+	const Color center_b = b_value(b,0,0,0,0);
 	for (int i_y=0; i_y<coarse_height; ++i_y) {
 		for (int i_x=0; i_x<coarse_width; ++i_x) {
 			const double* p_icv = &coarse_variables(i_x, i_y, 0);
@@ -313,15 +312,20 @@ void compute_initial_s(
 			for (size_t j_y=max<int>(0, i_y - center_y); j_y<max_j_y; ++j_y) {
 				for (int j_x=max<int>(0, i_x - center_x); j_x<max_j_x; ++j_x) {
 					if (i_x == j_x && i_y == j_y) continue;
-					Color b_ij = b_value(b,i_x,i_y,j_x,j_y);
+					const Color b_ij = b_value(b,i_x,i_y,j_x,j_y);
+					const double* p_jcv = &coarse_variables(j_x, j_y, 0);
 					for (size_t v=0; v<palette_size; ++v) {
-						Color b_ij2 = b_ij * p_icv[v];
-						const double* p_jcv = &coarse_variables(j_x, j_y, v);
+						const Color b_ij2 = b_ij * p_icv[v];
+						const double* p_jcv2 = p_jcv++;
+						double jcv = *p_jcv2;
 						Color* ps = s.pBuff_ + v * palette_size + v;
 // TODO: 変更画像、縦方向ではなく横方向に操作する。後で転置。
 						for (size_t alpha=v; alpha<palette_size; ++alpha) {
-							*ps += (*p_jcv++) * b_ij2;
+							++p_jcv2;
+							double njcv = *p_jcv2;
+							*ps += jcv * b_ij2;
 							ps += palette_size;
+							jcv = njcv;
 						}
 					}
 				}
@@ -343,27 +347,34 @@ void update_s(
 	const double delta
 	)
 {
-	const size_t palette_size  = s.width_;
-	const int coarse_width  = coarse_variables.width_;
-	const int coarse_height = coarse_variables.height_;
+	const size_t palette_size = s.width_;
 	const int center_x = (b.width_-1) / 2;
 	const int center_y = (b.height_-1) / 2;
-	const int max_i_x = min(coarse_width, j_x + center_x + 1);
-	const int max_i_y = min<int>(coarse_height, j_y + center_y + 1);
+	const size_t max_i_x = min<int>(coarse_variables.width_, j_x + center_x + 1);
+	const size_t max_i_y = min<int>(coarse_variables.height_, j_y + center_y + 1);
 	for (size_t i_y=max(0, j_y - center_y); i_y<max_i_y; ++i_y) {
 		for (size_t i_x=max(0, j_x - center_x); i_x<max_i_x; ++i_x) {
-			const Color delta_b_ij = delta * b_value(b,i_x,i_y,j_x,j_y);
 			if (i_x == j_x && i_y == j_y) continue;
+			const Color delta_b_ij = delta * b_value(b,i_x,i_y,j_x,j_y);
 			Color* ps = s[alpha];
 			const double* p_cv = &coarse_variables(i_x, i_y, 0);
+			double cv = *p_cv;
 			for (size_t v=0; v<=alpha; ++v) {
-				ps[v] += (*p_cv++) * delta_b_ij;
+				++p_cv;
+				double ncv = *p_cv;
+				*ps += cv * delta_b_ij;
+				++ps;
+				cv = ncv;
 			}
 			--p_cv;
-			ps += alpha;
+			--ps;
+			cv = *p_cv;
 			for (size_t v=alpha; v<palette_size; ++v) {
-				*ps += (*p_cv++) * delta_b_ij;
+				++p_cv;
+				double ncv = *p_cv;
+				*ps += cv * delta_b_ij;
 				ps += palette_size;
+				cv = ncv;
 			}
 		}
 	}
@@ -539,6 +550,8 @@ void spatial_color_quant(
 
 			// Compute 2*sum(j in extended neighborhood of i, j != i) b_ij
 
+			const int radius_width = (b.width_ - 1)/2;
+			const int radius_height = (b.height_ - 1)/2;
 			while (!visit_queue.empty()) {
 				// If we get to 10% above initial size, just revisit them all
 				if (visit_queue.size() > coarse_variables.width_*coarse_variables.height_*11.0/10) {
@@ -552,16 +565,45 @@ void spatial_color_quant(
 				// Compute (25)
 				Color p_i;
 				p_i.zero();
-				for (int y=0; y<b.height_; ++y) {
-					int j_y = y - center_y + i_y;
-					if (j_y < 0 || j_y >= coarse_variables.height_) continue;
-					for (int x=0; x<b.width_; ++x) {
-						int j_x = x - center_x + i_x;
-						if (i_x == j_x && i_y == j_y) continue;
-						if (j_x < 0 || j_x >= coarse_variables.width_) continue;
-						Color b_ij = b_value(b, i_x, i_y, j_x, j_y);
-						Color j_pal = (*j_palette_sum)[j_y][j_x];
-						p_i += b_ij * j_pal;
+				int sy = i_y - center_y;
+				int ey = sy + b.height_;
+				sy = max(sy, 0);
+				ey = min<int>(ey, coarse_variables.height_);
+				int ky = sy - i_y + radius_height;
+				if (ky < 0) {
+					sy += -ky;
+				}
+				ky = ey - 1 - i_y + radius_height;
+				if (ky >= b.height_) {
+					ey -= (ky - b.height_) + 1;
+				}
+				int sx = i_x - center_x;
+				int ex = sx + b.width_;
+				sx = max(sx, 0);
+				ex = min<int>(ex, coarse_variables.width_);
+				int kx = sx - i_x + radius_width;
+				if (kx < 0) {
+					sx += -kx;
+				}
+				kx = ex - 1 - i_x + radius_width;
+				if (kx >= b.width_) {
+					ex -= (kx - b.width_) + 1;
+				}
+				for (int y=sy; y<ey; ++y) {
+					int k_y = y - i_y + radius_height;
+					Color* cb = b[k_y] + sx - i_x + radius_width;
+					Color* cp = (*j_palette_sum)[y]+sx;
+					if (i_y == y) {
+						for (int x=sx; x<ex; ++x) {
+							Color c = (*cb++) * (*cp++);
+							if (i_x != x) {
+								p_i += c;
+							}
+						}
+					}else {
+						for (int x=sx; x<ex; ++x) {
+							p_i += (*cb++) * (*cp++);
+						}
 					}
 				}
 				p_i *= 2.0;
@@ -569,13 +611,13 @@ void spatial_color_quant(
 
 				double max_meanfield_log = -numeric_limits<double>::infinity();
 				double meanfield_sum = 0.0;
+				double minus_inv_temperature = -1.0 / temperature;
 				for (size_t v=0; v<num_colors; ++v) {
 					// Update m_{pi(i)v}^I according to (23)
 					// We can subtract an arbitrary factor to prevent overflow,
 					// since only the weight relative to the sum matters, so we
 					// will choose a value that makes the maximum e^100.
-					Color p_i2; p_i2 = p_i;
-					double m = -(palette[v].dot_product(p_i2 + middle_b.direct_product(palette[v]))) / temperature;
+					double m = palette[v].dot_product(p_i + middle_b.direct_product(palette[v])) * minus_inv_temperature;
 					meanfield_logs[v] = m;
 					if (m > max_meanfield_log) {
 						max_meanfield_log = m;