From 1f27c925a64180ac16f2153a6c22da9b3cf20dc2 Mon Sep 17 00:00:00 2001 From: ishmal Date: Sun, 23 Apr 2006 07:35:38 +0000 Subject: [PATCH] svd hacks --- src/extension/internal/odf.cpp | 147 ++++++++++++++++++++++++++------- 1 file changed, 117 insertions(+), 30 deletions(-) diff --git a/src/extension/internal/odf.cpp b/src/extension/internal/odf.cpp index edbc26c98..da3a4068b 100644 --- a/src/extension/internal/odf.cpp +++ b/src/extension/internal/odf.cpp @@ -147,9 +147,25 @@ public: NR::Matrix getV(); /** - * Return the three singular values along the diagonal + * Return the right singular vectors + * @return U x Vtransposed + */ + NR::Matrix getUVt(); + + /** + * Return the s[0] value */ - void getSingularValues(double &s0, double &s1, double &s2); + double getS0(); + + /** + * Return the s[1] value + */ + double getS1(); + + /** + * Return the s[2] value + */ + double getS2(); /** * Two norm @@ -182,7 +198,7 @@ private: }; -static double hypot(double a, double b) +static double svd_hypot(double a, double b) { double r; @@ -239,7 +255,7 @@ void SingularValueDecomposition::calculate() // Compute 2-norm of k-th column without under/overflow. s[k] = 0; for (int i = k; i < m; i++) { - s[k] = hypot(s[k],A[i][k]); + s[k] = svd_hypot(s[k],A[i][k]); } if (s[k] != 0.0) { if (A[k][k] < 0.0) { @@ -288,7 +304,7 @@ void SingularValueDecomposition::calculate() // Compute 2-norm without under/overflow. e[k] = 0; for (int i = k+1; i < n; i++) { - e[k] = hypot(e[k],e[i]); + e[k] = svd_hypot(e[k],e[i]); } if (e[k] != 0.0) { if (e[k+1] < 0.0) { @@ -409,8 +425,12 @@ void SingularValueDecomposition::calculate() int pp = p-1; int iter = 0; - double eps = pow(2.0,-52.0); - double tiny = pow(2.0,-966.0); + //double eps = pow(2.0,-52.0); + //double tiny = pow(2.0,-966.0); + //let's just calculate these now + //a double can be e ± 308.25, so this is safe + double eps = 2.22e-16; + double tiny = 1.6e-291; while (p > 0) { int k,kase; @@ -472,7 +492,7 @@ void SingularValueDecomposition::calculate() double f = e[p-2]; e[p-2] = 0.0; for (int j = p-2; j >= k; j--) { - double t = hypot(s[j],f); + double t = svd_hypot(s[j],f); double cs = s[j]/t; double sn = f/t; s[j] = t; @@ -497,7 +517,7 @@ void SingularValueDecomposition::calculate() double f = e[k-1]; e[k-1] = 0.0; for (int j = k; j < p; j++) { - double t = hypot(s[j],f); + double t = svd_hypot(s[j],f); double cs = s[j]/t; double sn = f/t; s[j] = t; @@ -520,10 +540,8 @@ void SingularValueDecomposition::calculate() // Calculate the shift. - double scale = 0.0; - double d = fabs(s[p-1]); - if (d>scale) scale=d; - d = fabs(s[p-2]); + double scale = fabs(s[p-1]); + double d = fabs(s[p-2]); if (d>scale) scale=d; d = fabs(e[p-2]); if (d>scale) scale=d; @@ -552,7 +570,7 @@ void SingularValueDecomposition::calculate() // Chase zeros. for (int j = k; j < p-1; j++) { - double t = hypot(f,g); + double t = svd_hypot(f,g); double cs = f/t; double sn = g/t; if (j != k) { @@ -569,7 +587,7 @@ void SingularValueDecomposition::calculate() V[i][j] = t; } } - t = hypot(f,g); + t = svd_hypot(f,g); cs = f/t; sn = g/t; s[j] = t; @@ -645,7 +663,7 @@ void SingularValueDecomposition::calculate() NR::Matrix SingularValueDecomposition::getU() { NR::Matrix mat(U[0][0], U[1][0], U[0][1], - U[1][1], U[2][0], U[2][1]); + U[1][1], U[0][2], U[1][2]); return mat; } @@ -657,19 +675,51 @@ NR::Matrix SingularValueDecomposition::getU() NR::Matrix SingularValueDecomposition::getV() { NR::Matrix mat(V[0][0], V[1][0], V[0][1], - V[1][1], V[2][0], V[2][1]); + V[1][1], V[0][2], V[1][2]); return mat; } /** - * Return the three singular values along the diagonal + * Return the right singular vectors + * @return U x Vtransposed + */ + +NR::Matrix SingularValueDecomposition::getUVt() +{ + //instead of sum(row*column), sum(column, column) + double a = U[0][0] * V[0][0] + U[1][0] * V[1][0]; + double b = U[0][0] * V[0][1] + U[1][0] * V[1][1]; + double c = U[0][1] * V[0][0] + U[1][1] * V[1][0]; + double d = U[0][1] * V[0][1] + U[1][1] * V[1][1]; + double e = U[0][2] * V[0][0] + U[1][2] * V[1][0]; + double f = U[0][2] * V[0][1] + U[1][2] * V[1][1]; + + NR::Matrix mat(a, b, c, d, e, f); + return mat; +} + +/** + * Return the s[0] value + */ +double SingularValueDecomposition::getS0() +{ + return s[0]; +} + +/** + * Return the s[1] value */ -void SingularValueDecomposition::getSingularValues( - double &s0, double &s1, double &s2) +double SingularValueDecomposition::getS1() { - s0 = s[0]; - s1 = s[1]; - s2 = s[2]; + return s[1]; +} + +/** + * Return the s[2] value + */ +double SingularValueDecomposition::getS2() +{ + return s[2]; } /** @@ -766,6 +816,45 @@ static std::string formatTransform(NR::Matrix &tf) } +/** + * An affine transformation Q may be decomposed via + * singular value decomposition into + * + * T = UDVt + * = (UDUt)UVt + * ('t' means transposed) + * where U and V are orthonormal matrices and D is a diagonal + * matrix. The decomposition may be interpreted as such: + * the image is firstly rotated by UVt. The image is then + * rotated by Ut, stretched in the coordinate directions + * by D then rotated back by U. The net effect is a slant operation in + * some tilt direction followed by an isotropic scale. If rot(x) + * is a matrix that rotates by x we can rewrite this as + * T = rot(-tau)( k 0, 0 1) rot(tau) rot(theta) S + * where S is a scaling matrix, k is a multiplying (contraction) + * factor related to slant, tau is the tilt direction and theta + * is the initial rotation angle. + */ +static void analyzeTransform(NR::Matrix &tf) +{ + SingularValueDecomposition svd(tf); + double scale1 = svd.getS0(); + double rotate = svd.getS1(); + double scale2 = svd.getS2(); + NR::Matrix u = svd.getU(); + NR::Matrix v = svd.getV(); + NR::Matrix uvt = svd.getUVt(); + g_message("s1:%f rot:%f s2:%f", scale1, rotate, scale2); + std::string us = formatTransform(u); + g_message("u:%s", us.c_str()); + std::string vs = formatTransform(v); + g_message("v:%s", vs.c_str()); + std::string uvts = formatTransform(uvt); + g_message("uvt:%s", uvts.c_str()); +} + + + /** * Method descends into the repr tree, converting image and style info @@ -1209,11 +1298,8 @@ bool OdfOutput::writeTree(Writer &outs, Inkscape::XML::Node *node) NR::Matrix itemTransform = item->transform; std::string itemTransformString = formatTransform(itemTransform); - - SingularValueDecomposition svd(itemTransform); - double scale1, rotate, scale2; - svd.getSingularValues(scale1, rotate, scale2); - g_message("s1:%f rot:%f s2:%f", scale1, rotate, scale2); + g_message("trans:%s", itemTransformString.c_str()); + analyzeTransform(itemTransform); std::string href = getAttribute(node, "xlink:href"); std::map::iterator iter = imageTable.find(href); @@ -1440,13 +1526,14 @@ bool OdfOutput::writeContent(ZipFile &zf, Inkscape::XML::Node *node) void OdfOutput::save(Inkscape::Extension::Output *mod, SPDocument *doc, gchar const *uri) { + //g_message("native file:%s\n", uri); + documentUri = URI(uri); + ZipFile zf; styleTable.clear(); styleLookupTable.clear(); imageTable.clear(); preprocess(zf, doc->rroot); - g_message("native file:%s\n", uri); - documentUri = URI(uri); if (!writeManifest(zf)) { -- 2.30.2