Experimental mesh smoothing code.

(bzr r13890)

9 files changed, 910 insertions, 27 deletions

diff --git a/src/attributes.cpp b/src/attributes.cpp
index 81b12c658..0b18b80cf 100644
--- a/src/attributes.cpp
+++ b/src/attributes.cpp
@@ -291,6 +291,7 @@ static SPStyleProp const props[] = {
     {SP_ATTR_FY, "fy"},
     /* SPMeshPatch */
     {SP_ATTR_TENSOR, "tensor"},
+    {SP_ATTR_SMOOTH, "smooth"},
     /* SPPattern */
     {SP_ATTR_PATTERNUNITS, "patternUnits"},
     {SP_ATTR_PATTERNCONTENTUNITS, "patternContentUnits"},
diff --git a/src/attributes.h b/src/attributes.h
index 660145a54..645aad00b 100644
--- a/src/attributes.h
+++ b/src/attributes.h
@@ -293,6 +293,7 @@ enum SPAttributeEnum {
     SP_ATTR_FY,
     /* SPMeshPatch */
     SP_ATTR_TENSOR,
+    SP_ATTR_SMOOTH,
     /* SPPattern */
     SP_ATTR_PATTERNUNITS,
     SP_ATTR_PATTERNCONTENTUNITS,
diff --git a/src/sp-gradient.h b/src/sp-gradient.h
index fbb48df62..f80d806fc 100644
--- a/src/sp-gradient.h
+++ b/src/sp-gradient.h
@@ -141,7 +141,8 @@ public:
 
     /** Composed array (for mesh gradients) */
     SPMeshNodeArray array;
-
+    SPMeshNodeArray array_smoothed; // Smoothed version of array
+    
     bool hasPatches() const;
 
 
diff --git a/src/sp-mesh-array.cpp b/src/sp-mesh-array.cpp
index 8bfe23656..c1ead565f 100644
--- a/src/sp-mesh-array.cpp
+++ b/src/sp-mesh-array.cpp
@@ -32,7 +32,7 @@
  * Authors:
  *   Tavmjong Bah <tavmjong@free.fr>
  *
- * Copyrigt  (C) 2012 Tavmjong Bah
+ * Copyright  (C) 2012, 2015 Tavmjong Bah
  *
  * Released under GNU GPL, read the file 'COPYING' for more information
  */
@@ -580,6 +580,47 @@ SPMeshNodeArray::SPMeshNodeArray( SPMeshGradient *mg ) {
 
 };
 
+
+// Copy constructor
+SPMeshNodeArray::SPMeshNodeArray( const SPMeshNodeArray& rhs ) {
+
+    built = false;
+    mg = NULL;
+    drag_valid = false;
+
+    nodes = rhs.nodes; // This only copies the pointers but it does size the vector of vectors.
+
+    for( unsigned i=0; i < nodes.size(); ++i ) {
+        for( unsigned j=0; j < nodes[i].size(); ++j ) {
+            nodes[i][j] = new SPMeshNode( *rhs.nodes[i][j] ); // Copy data.
+        }
+    }
+};
+
+
+// Copy assignment operator
+SPMeshNodeArray& SPMeshNodeArray::operator=( const SPMeshNodeArray& rhs ) {
+
+    if( this == &rhs ) return *this;
+
+    clear(); // Clear any existing array.
+
+    built = false;
+    mg = NULL;
+    drag_valid = false;
+
+    nodes = rhs.nodes; // This only copies the pointers but it does size the vector of vectors.
+
+    for( unsigned i=0; i < nodes.size(); ++i ) {
+        for( unsigned j=0; j < nodes[i].size(); ++j ) {
+            nodes[i][j] = new SPMeshNode( *rhs.nodes[i][j] ); // Copy data.
+        }
+    }
+    
+    return *this;
+};
+
+
 void SPMeshNodeArray::read( SPMeshGradient *mg_in ) {
 
     mg = mg_in;
@@ -1382,6 +1423,543 @@ void SPMeshNodeArray::print() {
 };
 
 
+
+// Find the slopes at start and end for Hermite interpolation.
+// Smooth using Hermite interpolation.
+// Inputs are:
+//  pb: color value before patch
+//  p0: color value start of patch
+//  p1: color value end of patch
+//  pa: color value after patch
+//  is_first:  If first patch in row/column
+//  is_last:   If last patch in row/column
+//  type: Type smoothing
+// Output:
+//  m0: slope of Hermite function at start.
+//  m1: slope of Hermite function at end.
+void find_slopes( const double &pb, const double &p0, const double &p1, const double &pa,
+                  const bool &is_first, const bool &is_last, const SPMeshSmooth type,
+                  double &m0, double &m1 ) {
+
+    // We use Hermite interpolation. We have end points, we need tangents.
+
+    // Try various ways of finding tangents m0, m1
+
+    // Default to Catmul-Rom (assumes pb and pa already calculatedd)
+    m0 = (p1 - pb)/2.0;
+    m1 = (pa - p0)/2.0;
+
+    bool parabolic = false;  // Require end patches to be parabolic
+    switch (type) {
+        case SP_MESH_SMOOTH_SMOOTH1:
+            // Flat
+            m0 = 0.0;
+            m1 = 0.0;
+            break;
+        case SP_MESH_SMOOTH_SMOOTH2:
+            // Catmul-Rom, standard end treatment. Double first/last point.
+            if( is_first ) {
+                m0 = (p1-p0)/2.0;
+            }
+            if( is_last ) {
+                m1 = (p1-p0)/2.0;
+            }
+            break;
+        case SP_MESH_SMOOTH_SMOOTH3:
+            // Catmul-Rom, standard end treatment. Reflect first/last point.
+            if( is_first ) {
+                m0 = (p1-p0);
+            }
+            if( is_last ) {
+                m1 = (p1-p0);
+            }
+            break;
+        case SP_MESH_SMOOTH_SMOOTH4:
+            // Catmul-Rom, Parabolic ends
+            parabolic = true;
+            break;
+        case SP_MESH_SMOOTH_SMOOTH:
+        case SP_MESH_SMOOTH_SMOOTH5:
+            // Catmul-Rom, Parabolic ends, no color min/max in middle of patch.
+            parabolic = true;
+
+            if( (pb > p0 && p1 > p0) ||
+                (pb < p0 && p1 < p0) ) {
+                // tangents flat at min/max
+                m0 = 0;
+            } else {
+                // ensure we don't overshoot
+                if( fabs(m0) > fabs(3*(p1-p0)) ) {
+                    m0 = 3*(p1-p0);
+                }
+                if( fabs(m0) > fabs(3*(p0-pb)) ) {
+                    m0 = 3*(p0-pb);
+                }
+            }
+            if( (p0 > p1 && pa > p1) ||
+                (p0 < p1 && pa < p1) ) {
+                // tangents flat at min/max
+                m1 = 0;
+            } else {
+                // ensure we don't overshoot
+                if( fabs(m1) > fabs(3*(pa-p1)) ) {
+                    m1 = 3*(pa-p1);
+                }
+                if( fabs(m1) > fabs(3*(p1-p0)) ) {
+                    m1 = 3*(p1-p0);
+                }
+            }
+            break;
+        case SP_MESH_SMOOTH_NONE:
+        default:
+            std::cerr << "find_slopes() Invalid smoothing type." << std::endl;
+            break;
+    }
+
+    // Force end patches to be parabolic
+    if( parabolic ) {
+        if( is_first ) {
+            // Constraint for parabola
+            m0 = 2.0*(p1-p0) - m1;
+            if ( ((p1-p0) < 0 && m0 > 0) || ((p1-p0) > 0 && m0 < 0 ) ) {
+                m0 = 0;  // Prevent overshooting start value;
+            }
+        } else if( is_last ) {
+            // Constraint for parabola
+            m1 = 2.0*(p1-p0) - m0;
+            if ( ((p1-p0) < 0 && m1 > 0) || ((p1-p0) > 0 && m1 < 0 ) ) {
+                m1 = 0;  // Prevent overshooting end value;
+            }
+        }
+    }
+
+    // std::cout << "  pb: " << pb
+    //           << "  p0: " << p0
+    //           << "  p1: " << p1
+    //           << "  pa: " << pa
+    //           << "  m0: " << m0
+    //           << "  m1: " << m1 << std::endl;
+}
+
+double hermite( const double p0, const double p1, const double m0, const double m1, const double t ) {
+    double t2 = t*t;
+    double t3 = t2*t;
+
+    double result = (2.0*t3 - 3.0*t2 +1.0) * p0
+                  + (t3 - 2.0*t2 + t)      * m0
+                  + (-2.0*t3 + 3.0*t2)     * p1
+                  + (t3 -t2)               * m1;
+
+    return result;
+}
+
+
+/**
+   Fill 'smooth' with a smoothed version of the array by subdividing each patch into smaller patches.
+*/
+void SPMeshNodeArray::smooth( SPMeshNodeArray* smooth, SPMeshSmooth type ) {
+
+    *smooth = *this;  // Deep copy via copy assignment constructor, smooth cleared before copy
+    // std::cout << "SPMeshNodeArray::smooth(): " << this->patch_rows() << " " << smooth->patch_rows() << std::endl;
+    // std::cout << "   " << smooth << " " << this << std::endl;
+    // Next split each patch into 8x8 smaller patches.
+    
+    // Do rows first.
+
+    // Split each row into eight rows.
+    // Must do it from end so inserted rows don't mess up indexing 
+    for( int i = smooth->patch_rows() - 1; i >= 0; --i ) {
+        smooth->split_row( i, unsigned(8) );
+    }
+
+    // Update color values (every third node is a corner)
+    for( unsigned i = 0; i < this->patch_rows(); ++i ) { // i is orignal patch index
+
+        bool is_first_row = (i == 0);
+        bool is_last_row = (i == this->patch_rows() - 1 );
+        //std::cout << "  last row: " << smooth->patch_rows()/8 - 1 << " " << is_last_row << std::endl;
+        for( unsigned j = 0; j < smooth->patch_columns()+1; ++j ) { // j is smooth patch index
+
+            // Can't use guint32 since delta can be negative
+            float pb[3]; // Point before patch
+            float p0[3]; // Point at start of patch
+            float p1[3]; // Point at end of patch
+            float pa[3]; // Point after patch
+            float result[3][8];
+            sp_color_get_rgb_floatv( &this->nodes[  i   *3 ][ j*3 ]->color, p0 ); 
+            sp_color_get_rgb_floatv( &this->nodes[ (i+1)*3 ][ j*3 ]->color, p1 );
+            if( !is_first_row ) {
+                sp_color_get_rgb_floatv( &this->nodes[ (i-1)*3 ][ j*3 ]->color, pb );
+            } else {
+                pb[0] = 2.0*p0[0] - p1[0];
+                pb[1] = 2.0*p0[1] - p1[1];
+                pb[2] = 2.0*p0[2] - p1[2];
+            }
+            if( !is_last_row ) {
+                sp_color_get_rgb_floatv( &this->nodes[ (i+2)*3 ][ j*3 ]->color, pa );
+            } else {
+                pa[0] = 2.0*p1[0] - p0[0];
+                pa[1] = 2.0*p1[1] - p0[1];
+                pa[2] = 2.0*p1[2] - p0[2];
+            }
+
+            for( unsigned n = 0; n < 3; ++n ) { // Loop over colors
+
+                // We use Hermite interpolation. We have end points, we need tangents.
+                double m0 = 0;
+                double m1 = 0;
+                find_slopes( pb[n], p0[n], p1[n], pa[n], is_first_row, is_last_row, type, m0, m1 );
+        
+                for( unsigned k = 1; k < 8; ++k ) {
+                    double t = k/8.0;
+                    // Cubic Hermite (four constraints)
+                    result[n][k] = hermite( p0[n], p1[n], m0, m1, t );
+                    // Clamp to allowed values
+                    if( result[n][k] > 1.0 )
+                        result[n][k] = 1.0;
+                    if( result[n][k] < 0.0 )
+                        result[n][k] = 0.0;
+                }
+            }
+
+            for( unsigned k = 1; k < 8; ++k ) {
+                smooth->nodes[ (i*8+k)*3 ][ j*3 ]->color.set( result[0][k], result[1][k], result[2][k] ); 
+            }
+        }
+    }
+
+    // Split each column into eight columns.
+    // Must do it from end so inserted columns don't mess up indexing 
+    for( int i = smooth->patch_columns() - 1; i >= 0; --i ) {
+        smooth->split_column( i, (unsigned)8 );
+    }
+
+    // Update color values (every third node is a corner)
+    for( unsigned i = 0; i < this->patch_columns(); ++i ) { // i is orignal patch index
+
+        bool is_first_column = (i == 0);
+        bool is_last_column = (i == this->patch_columns() - 1 );
+        //std::cout << "  last column: " << smooth->patch_columns()/8 - 1 << " " << is_last_column << std::endl;
+        for( unsigned j = 0; j < smooth->patch_rows()+1; ++j ) { // j is smooth patch index
+
+            // Can't use guint32 since delta can be negative
+            float pb[3]; // Point before patch
+            float p0[3]; // Point at start of patch
+            float p1[3]; // Point at end of patch
+            float pa[3]; // Point after patch
+            float result[3][8];
+            sp_color_get_rgb_floatv( &smooth->nodes[ j*3 ][  i   *3*8 ]->color, p0 ); 
+            sp_color_get_rgb_floatv( &smooth->nodes[ j*3 ][ (i+1)*3*8 ]->color, p1   );
+            if( !is_first_column ) {
+                sp_color_get_rgb_floatv( &smooth->nodes[ j*3 ][ (i-1)*3*8 ]->color, pb );
+            } else {
+                pb[0] = 2.0*p0[0] - p1[0];
+                pb[1] = 2.0*p0[1] - p1[1];
+                pb[2] = 2.0*p0[2] - p1[2];
+            }
+            if( !is_last_column ) {
+                sp_color_get_rgb_floatv( &smooth->nodes[ j*3 ][ (i+2)*3*8 ]->color, pa );
+            } else {
+                pa[0] = 2.0*p1[0] - p0[0];
+                pa[1] = 2.0*p1[1] - p0[1];
+                pa[2] = 2.0*p1[2] - p0[2];
+            }
+
+            for( unsigned n = 0; n < 3; ++n ) { // Loop over colors
+
+                // We use Hermite interpolation. We have end points, we need tangents.
+                double m0 = 0;
+                double m1 = 0;
+                find_slopes( pb[n], p0[n], p1[n], pa[n], is_first_column, is_last_column, type, m0, m1 );
+        
+                for( unsigned k = 1; k < 8; ++k ) {
+                    double t = k/8.0;
+                    // Cubic Hermite (four constraints)
+                    result[n][k] = hermite( p0[n], p1[n], m0, m1, t );
+                    // Clamp to allowed values
+                    if( result[n][k] > 1.0 )
+                        result[n][k] = 1.0;
+                    if( result[n][k] < 0.0 )
+                        result[n][k] = 0.0;
+                }
+            }
+
+            for( unsigned k = 1; k < 8; ++k ) {
+                smooth->nodes[ j*3 ][ (i*8+k)*3 ]->color.set( result[0][k], result[1][k], result[2][k] ); 
+            }
+        }
+    }
+}
+
+class SPMeshSmoothCorner {
+
+public:
+    SPMeshSmoothCorner() {
+        for( unsigned i = 0; i < 3; ++i ) {
+            for( unsigned j = 0; j < 4; ++j ) {
+                g[i][j] = 0;
+            }
+        }
+    }
+    
+    double g[3][4]; // 3 colors, 4 parameters: f, f_x, f_y, f_xy
+    // double x[4]; // Lengths between neigboring points  x, y, xy, yx
+};
+
+// Find slope at point 1 given values at previous and next points
+double find_slope1( double p0, double p1, double p2 ) {
+
+    double slope = (p2 - p0)/2.0;  // Simple Catmul-Rom condition
+    if( ( p0 > p1 && p1 < p2 ) ||
+        ( p0 < p1 && p1 > p2 ) ) {
+        // At minimum or maximum, use slope of zero
+        slope = 0;
+    } else {
+        // Ensure we don't overshoot
+        if( fabs(slope) > fabs(3*(p1-p0)) ) {
+            slope = 3*(p1-p0);
+        }
+        if( fabs(slope) > fabs(3*(p2-p1)) ) {
+            slope = 3*(p2-p1);
+        }
+    }
+    return slope;
+};
+
+
+// Find slope at point 0 given values at previous and next points
+double find_slope2( double pmm, double ppm, double pmp, double ppp, double p0 ) {
+
+    // pmm == d[i-1][j-1], ...  'm' is minus, 'p' is plus
+    double slope = (ppp - ppm - pmp + pmm)/2.0;
+    if( (ppp > p0 && ppm > p0 && pmp > p0 && pmm > 0) ||
+        (ppp < p0 && ppm < p0 && pmp < p0 && pmm < 0) ) {
+        // At minimum or maximum, use slope of zero
+        slope = 0;
+    } else {
+        // Don't really know what to do here
+        if( fabs(slope) > fabs(3*(ppp-p0)) ) {
+            slope = 3*(ppp-p0);
+        }
+        if( fabs(slope) > fabs(3*(pmp-p0)) ) {
+            slope = 3*(pmp-p0);
+        }
+        if( fabs(slope) > fabs(3*(ppm-p0)) ) {
+            slope = 3*(ppm-p0);
+        }
+        if( fabs(slope) > fabs(3*(pmm-p0)) ) {
+            slope = 3*(pmm-p0);
+        }
+    }
+    return slope;
+}
+
+// https://en.wikipedia.org/wiki/Bicubic_interpolation
+void invert( const double v[16], double alpha[16] ) {
+
+    const double  A[16][16] = {
+
+        { 1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0 },
+        { 0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0 },
+        {-3, 3, 0, 0, -2,-1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0 },
+        { 2,-2, 0, 0,  1, 1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0 },
+        { 0, 0, 0, 0,  0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0 },
+        { 0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  1, 0, 0, 0 },
+        { 0, 0, 0, 0,  0, 0, 0, 0, -3, 3, 0, 0, -2,-1, 0, 0 },
+        { 0, 0, 0, 0,  0, 0, 0, 0,  2,-2, 0, 0,  1, 1, 0, 0 },
+        {-3, 0, 3, 0,  0, 0, 0, 0, -2, 0,-1, 0,  0, 0, 0, 0 },
+        { 0, 0, 0, 0, -3, 0, 3, 0,  0, 0, 0, 0, -2, 0,-1, 0 },
+        { 9,-9,-9, 9,  6, 3,-6,-3,  6,-6, 3,-3,  4, 2, 2, 1 },
+        {-6, 6, 6,-6, -3,-3, 3, 3, -4, 4,-2, 2, -2,-2,-1,-1 },
+        { 2, 0,-2, 0,  0, 0, 0, 0,  1, 0, 1, 0,  0, 0, 0, 0 },
+        { 0, 0, 0, 0,  2, 0,-2, 0,  0, 0, 0, 0,  1, 0, 1, 0 },
+        {-6, 6, 6,-6, -4,-2, 4, 2, -3, 3,-3, 3, -2,-1,-2,-1 },
+        { 4,-4,-4, 4,  2, 2,-2,-2,  2,-2, 2,-2,  1, 1, 1, 1 }
+    };
+
+    for( unsigned i = 0; i < 16; ++i ) {
+        alpha[i] = 0;
+        for( unsigned j = 0; j < 16; ++j ) {
+            alpha[i] += A[i][j]*v[j];
+        }
+    }
+}
+
+double sum( const double alpha[16], const double& x, const double& y ) {
+
+    double result = 0;
+    
+    double xx = x*x;
+    double xxx = xx * x;
+    double yy = y*y;
+    double yyy = yy * y;
+
+    result += alpha[  0 ];
+    result += alpha[  1 ] * x;
+    result += alpha[  2 ] * xx;
+    result += alpha[  3 ] * xxx;
+    result += alpha[  4 ] * y;
+    result += alpha[  5 ] * y * x;
+    result += alpha[  6 ] * y * xx;
+    result += alpha[  7 ] * y * xxx;
+    result += alpha[  8 ] * yy;
+    result += alpha[  9 ] * yy * x;
+    result += alpha[ 10 ] * yy * xx;
+    result += alpha[ 11 ] * yy * xxx;
+    result += alpha[ 12 ] * yyy;
+    result += alpha[ 13 ] * yyy * x;
+    result += alpha[ 14 ] * yyy * xx;
+    result += alpha[ 15 ] * yyy * xxx;
+    
+    return result;
+}
+
+/**
+   Fill 'smooth' with a smoothed version of the array by subdividing each patch into smaller patches.
+*/
+void SPMeshNodeArray::smooth2( SPMeshNodeArray* smooth, SPMeshSmooth type ) {
+
+    
+    *smooth = *this;  // Deep copy via copy assignment constructor, smooth cleared before copy
+    // std::cout << "SPMeshNodeArray::smooth2(): " << this->patch_rows() << " " << smooth->patch_columns() << std::endl;
+    // std::cout << "  " << smooth << " " << this << std::endl;
+
+    // Find derivatives at corners
+
+    // Create array of corner points
+    std::vector< std::vector <SPMeshSmoothCorner> > d; 
+    d.resize( smooth->patch_rows() + 1 );
+    for( unsigned i = 0; i < d.size(); ++i ) {
+        d[i].resize( smooth->patch_columns() + 1 );
+        for( unsigned j = 0; j < d[i].size(); ++j ) {
+            float rgb_color[3];
+            sp_color_get_rgb_floatv( &this->nodes[ i*3 ][ j*3 ]->color, rgb_color );
+            d[i][j].g[0][0] =  rgb_color[ 0 ];
+            d[i][j].g[1][0] =  rgb_color[ 1 ];
+            d[i][j].g[2][0] =  rgb_color[ 2 ];
+        }
+    }
+
+    // Calculate interior derivatives
+    for( unsigned i = 1; i < d.size()-1; ++i ) {
+        for( unsigned j = 1; j < d[i].size()-1; ++j ) {
+            for( unsigned k = 0; k < 3; ++k ) { // Loop over colors
+                if( type == SP_MESH_SMOOTH_SMOOTH7 || type == SP_MESH_SMOOTH_SMOOTH ) {
+                    d[i][j].g[k][1] = find_slope1( d[i-1][j].g[k][0], d[i][j].g[k][0], d[i+1][j].g[k][0] );
+                    d[i][j].g[k][2] = find_slope1( d[i][j-1].g[k][0], d[i][j].g[k][0], d[i][j+1].g[k][0] );
+                    d[i][j].g[k][3] = find_slope2( d[i-1][j-1].g[k][0], d[i+1][j-1].g[k][0],
+                                                   d[i-1][j+1].g[k][0], d[i-1][j-1].g[k][0],
+                                                   d[i][j].g[k][0] );
+                } else {
+                    // Catmul-Rom
+                    d[i][j].g[k][1] = (d[i+1][j].g[k][0] - d[i-1][j].g[k][0])/2.0;
+                    d[i][j].g[k][2] = (d[i][j+1].g[k][0] - d[i][j-1].g[k][0])/2.0;
+                }
+            }
+        }
+    }
+    
+    // Calculate exterior derivatives
+    for( unsigned j = 1; j< d[0].size()-1; ++j ) {
+        for( unsigned k = 0; k < 3; ++k ) { // Loop over colors
+            unsigned z = d.size()-1;
+            if( type == SP_MESH_SMOOTH_SMOOTH7 || type == SP_MESH_SMOOTH_SMOOTH ) {
+                // Parabolic
+                d[0][j].g[k][1] = 2.0*(d[1][j].g[k][0] - d[0  ][j].g[k][0]) - d[1][j].g[k][1];
+                d[z][j].g[k][1] = 2.0*(d[z][j].g[k][0] - d[z-1][j].g[k][0]) - d[z][j].g[k][1];
+            } else {
+                // Catmul-Rom
+                d[0][j].g[k][1] = (d[1][j].g[k][0] - d[0  ][j].g[k][0])/2.0;
+                d[z][j].g[k][1] = (d[z][j].g[k][0] - d[z-1][j].g[k][0])/2.0;
+            }
+        }
+    }
+
+    for( unsigned i = 1; i< d.size()-1; ++i ) {
+        for( unsigned k = 0; k < 3; ++k ) { // Loop over colors
+            unsigned z = d[0].size()-1;
+            if( type == SP_MESH_SMOOTH_SMOOTH7 || type == SP_MESH_SMOOTH_SMOOTH ) {
+                // Parabolic
+                d[i][0].g[k][2] = 2.0*(d[i][1].g[k][0] - d[i][0  ].g[k][0]) - d[i][0].g[k][1];
+                d[i][z].g[k][2] = 2.0*(d[i][z].g[k][0] - d[i][z-1].g[k][0]) - d[i][z].g[k][1];
+            } else {
+                // Catmul-Rom
+                d[i][0].g[k][2] = (d[i][1].g[k][0] - d[i][0  ].g[k][0])/2.0;
+                d[i][z].g[k][2] = (d[i][z].g[k][0] - d[i][z-1].g[k][0])/2.0;
+            }
+        }
+    }
+
+    // Leave outside corner derivatives at zero.
+    
+    // Next split each patch into 8x8 smaller patches.
+    
+    // Split each row into eight rows.
+    // Must do it from end so inserted rows don't mess up indexing 
+    for( int i = smooth->patch_rows() - 1; i >= 0; --i ) {
+        smooth->split_row( i, unsigned(8) );
+    }
+
+    // Split each column into eight columns.
+    // Must do it from end so inserted columns don't mess up indexing 
+    for( int i = smooth->patch_columns() - 1; i >= 0; --i ) {
+        smooth->split_column( i, (unsigned)8 );
+    }
+
+    // Fill new patches
+    for( unsigned i = 0; i < this->patch_rows(); ++i ) {
+        for( unsigned j = 0; j < this->patch_columns(); ++j ) {
+
+            // Temp loop over 0..8 to get last column/row edges
+            float r[3][9][9]; // result
+            for( unsigned m = 0; m < 3; ++m ) {
+
+                double v[16];
+                v[ 0] = d[i  ][j  ].g[m][0];
+                v[ 1] = d[i+1][j  ].g[m][0];
+                v[ 2] = d[i  ][j+1].g[m][0];
+                v[ 3] = d[i+1][j+1].g[m][0];
+                v[ 4] = d[i  ][j  ].g[m][1];
+                v[ 5] = d[i+1][j  ].g[m][1];
+                v[ 6] = d[i  ][j+1].g[m][1];
+                v[ 7] = d[i+1][j+1].g[m][1];
+                v[ 8] = d[i  ][j  ].g[m][2];
+                v[ 9] = d[i+1][j  ].g[m][2];
+                v[10] = d[i  ][j+1].g[m][2];
+                v[11] = d[i+1][j+1].g[m][2];
+                v[12] = d[i  ][j  ].g[m][3];
+                v[13] = d[i+1][j  ].g[m][3];
+                v[14] = d[i  ][j+1].g[m][3];
+                v[15] = d[i+1][j+1].g[m][3];
+
+                double alpha[16];
+                invert( v, alpha );
+                
+                for( unsigned k = 0; k < 9; ++k ) {
+                    for( unsigned l = 0; l < 9; ++l ) {
+                        double x = k/8.0;
+                        double y = l/8.0;
+                        r[m][k][l] = sum( alpha, x, y );
+                        // Clamp to allowed values
+                        if( r[m][k][l] > 1.0 )
+                            r[m][k][l] = 1.0;
+                        if( r[m][k][l] < 0.0 )
+                            r[m][k][l] = 0.0;
+                    }
+                }
+
+            } // Loop over colors
+
+            for( unsigned k = 0; k < 9; ++k ) {
+                for( unsigned l = 0; l < 9; ++l ) {
+                    // Every third node is a corner node
+                    smooth->nodes[ (i*8+k)*3 ][(j*8+l)*3 ]->color.set( r[0][k][l], r[1][k][l], r[2][k][l] );
+                }
+            }
+        }
+    }
+}
+
 /**
    Number of patch rows.
 */
diff --git a/src/sp-mesh-array.h b/src/sp-mesh-array.h
index 330f6d87b..19f8634f8 100644
--- a/src/sp-mesh-array.h
+++ b/src/sp-mesh-array.h
@@ -45,6 +45,18 @@
 // For color picking
 #include "sp-item.h"
 
+enum SPMeshSmooth {
+  SP_MESH_SMOOTH_NONE,
+  SP_MESH_SMOOTH_SMOOTH,
+  SP_MESH_SMOOTH_SMOOTH1,
+  SP_MESH_SMOOTH_SMOOTH2,
+  SP_MESH_SMOOTH_SMOOTH3,
+  SP_MESH_SMOOTH_SMOOTH4,
+  SP_MESH_SMOOTH_SMOOTH5,
+  SP_MESH_SMOOTH_SMOOTH6,
+  SP_MESH_SMOOTH_SMOOTH7
+};
+
 enum NodeType {
   MG_NODE_TYPE_UNKNOWN,
   MG_NODE_TYPE_CORNER,
@@ -148,6 +160,9 @@ public:
 
   SPMeshNodeArray() { built = false; mg = NULL; drag_valid = false; };
   SPMeshNodeArray( SPMeshGradient *mg );
+  SPMeshNodeArray( const SPMeshNodeArray& rhs );
+  SPMeshNodeArray& operator=(const SPMeshNodeArray& rhs);
+
   ~SPMeshNodeArray() { clear(); };
   bool built;
 
@@ -157,6 +172,10 @@ public:
   void clear();
   void print();
 
+  // Fill 'smooth' with a smoothed version by subdividing each patch.
+  void smooth( SPMeshNodeArray* smooth, SPMeshSmooth type);
+  void smooth2( SPMeshNodeArray* smooth, SPMeshSmooth type);
+
   // Get size of patch
   unsigned int patch_rows();
   unsigned int patch_columns();
diff --git a/src/sp-mesh-gradient.cpp b/src/sp-mesh-gradient.cpp
index bf28164a9..65fbe2e46 100644
--- a/src/sp-mesh-gradient.cpp
+++ b/src/sp-mesh-gradient.cpp
@@ -21,7 +21,7 @@ namespace {
  * Mesh Gradient
  */
 //#define MESH_DEBUG
-SPMeshGradient::SPMeshGradient() : SPGradient() {
+SPMeshGradient::SPMeshGradient() : SPGradient(), smooth(SP_MESH_SMOOTH_NONE), smooth_set(false) {
     // Start coordinate of mesh
     this->x.unset(SVGLength::NONE, 0.0, 0.0);
     this->y.unset(SVGLength::NONE, 0.0, 0.0);
@@ -36,6 +36,8 @@ void SPMeshGradient::build(SPDocument *document, Inkscape::XML::Node *repr) {
     // Start coordinate of mesh
     this->readAttr( "x" );
     this->readAttr( "y" );
+
+    this->readAttr( "smooth" );
 }
 
 
@@ -57,6 +59,39 @@ void SPMeshGradient::set(unsigned key, gchar const *value) {
             this->requestModified(SP_OBJECT_MODIFIED_FLAG);
             break;
 
+        case SP_ATTR_SMOOTH:
+	    if (value) {
+	      if (!strcmp(value, "none")) {
+		this->smooth = SP_MESH_SMOOTH_NONE;
+	      } else if (!strcmp(value, "smooth")) {
+		this->smooth = SP_MESH_SMOOTH_SMOOTH;
+	      } else if (!strcmp(value, "smooth1")) {
+		this->smooth = SP_MESH_SMOOTH_SMOOTH1;
+	      } else if (!strcmp(value, "smooth2")) {
+		this->smooth = SP_MESH_SMOOTH_SMOOTH2;
+	      } else if (!strcmp(value, "smooth3")) {
+		this->smooth = SP_MESH_SMOOTH_SMOOTH3;
+	      } else if (!strcmp(value, "smooth4")) {
+		this->smooth = SP_MESH_SMOOTH_SMOOTH4;
+	      } else if (!strcmp(value, "smooth5")) {
+		this->smooth = SP_MESH_SMOOTH_SMOOTH5;
+	      } else if (!strcmp(value, "smooth6")) {
+		this->smooth = SP_MESH_SMOOTH_SMOOTH6;
+	      } else if (!strcmp(value, "smooth7")) {
+		this->smooth = SP_MESH_SMOOTH_SMOOTH7;
+	      } else {
+		std::cout << "SPMeshGradient::set(): invalid value " << value << std::endl;
+	      }
+	      this->smooth_set = TRUE;
+	    } else {
+	      std::cout << "SPMeshGradient::set() No value " << std::endl;
+	      this->smooth = SP_MESH_SMOOTH_NONE;
+	      this->smooth_set = FALSE;
+	    }
+
+            this->requestModified(SP_OBJECT_MODIFIED_FLAG);
+            break;
+
         default:
             SPGradient::set(key, value);
             break;
@@ -83,6 +118,41 @@ Inkscape::XML::Node* SPMeshGradient::write(Inkscape::XML::Document *xml_doc, Ink
     	sp_repr_set_svg_double(repr, "y", this->y.computed);
     }
 
+    if ((flags & SP_OBJECT_WRITE_ALL) || this->smooth_set) {
+        switch (this->smooth) {
+	    case SP_MESH_SMOOTH_SMOOTH:
+	      repr->setAttribute("smooth", "smooth");
+	      break;
+	    case SP_MESH_SMOOTH_SMOOTH1:
+	      repr->setAttribute("smooth", "smooth1");
+	      break;
+	    case SP_MESH_SMOOTH_SMOOTH2:
+	      repr->setAttribute("smooth", "smooth2");
+	      break;
+	    case SP_MESH_SMOOTH_SMOOTH3:
+	      repr->setAttribute("smooth", "smooth3");
+	      break;
+	    case SP_MESH_SMOOTH_SMOOTH4:
+	      repr->setAttribute("smooth", "smooth4");
+	      break;
+	    case SP_MESH_SMOOTH_SMOOTH5:
+	      repr->setAttribute("smooth", "smooth5");
+	      break;
+	    case SP_MESH_SMOOTH_SMOOTH6:
+	      repr->setAttribute("smooth", "smooth6");
+	      break;
+	    case SP_MESH_SMOOTH_SMOOTH7:
+	      repr->setAttribute("smooth", "smooth7");
+	      break;
+	    case SP_MESH_SMOOTH_NONE:
+	      repr->setAttribute("smooth", "none");
+	      break;
+	    default:
+	      // Do nothing
+	      break;
+	}
+    }
+
     SPGradient::write(xml_doc, repr, flags);
 
     return repr;
@@ -117,14 +187,38 @@ cairo_pattern_t* SPMeshGradient::pattern_new(cairo_t * /*ct*/,
     cairo_pattern_t *cp = NULL;
 
 #if CAIRO_VERSION >= CAIRO_VERSION_ENCODE(1, 11, 4)
-    SPMeshNodeArray* array = &(this->array);
+    SPMeshNodeArray* my_array = &array;
+
+    if( smooth_set ) {
+      switch (smooth) {
+        case SP_MESH_SMOOTH_NONE:
+	  // std::cout << "SPMeshGradient::pattern_new: no smoothing" << std::endl;
+	  break;
+	case SP_MESH_SMOOTH_SMOOTH1:
+	case SP_MESH_SMOOTH_SMOOTH2:
+	case SP_MESH_SMOOTH_SMOOTH3:
+	case SP_MESH_SMOOTH_SMOOTH4:
+	case SP_MESH_SMOOTH_SMOOTH5:
+	  // std::cout << "SPMeshGradient::pattern_new: calling array.smooth" << std::endl;
+	  array.smooth( &array_smoothed, smooth );
+	  my_array = &array_smoothed;
+	  break;
+	case SP_MESH_SMOOTH_SMOOTH:
+	case SP_MESH_SMOOTH_SMOOTH6:
+	case SP_MESH_SMOOTH_SMOOTH7:
+	  // std::cout << "SPMeshGradient::pattern_new: calling array.smooth2" << std::endl;
+	  array.smooth2( &array_smoothed, smooth );
+	  my_array = &array_smoothed;
+	  break;
+	}
+    }
 
     cp = cairo_pattern_create_mesh();
 
-    for( unsigned int i = 0; i < array->patch_rows(); ++i ) {
-        for( unsigned int j = 0; j < array->patch_columns(); ++j ) {
+    for( unsigned int i = 0; i < my_array->patch_rows(); ++i ) {
+        for( unsigned int j = 0; j < my_array->patch_columns(); ++j ) {
 
-            SPMeshPatchI patch( &(array->nodes), i, j );
+            SPMeshPatchI patch( &(my_array->nodes), i, j );
 
             cairo_mesh_pattern_begin_patch( cp );
             cairo_mesh_pattern_move_to( cp, patch.getPoint( 0, 0 )[X], patch.getPoint( 0, 0 )[Y] );
diff --git a/src/sp-mesh-gradient.h b/src/sp-mesh-gradient.h
index 4df753f62..343cae840 100644
--- a/src/sp-mesh-gradient.h
+++ b/src/sp-mesh-gradient.h
@@ -19,7 +19,8 @@ public:
 
     SVGLength x;  // Upper left corner of mesh
     SVGLength y;  // Upper right corner of mesh
-
+    SPMeshSmooth smooth;
+    bool smooth_set;
     virtual cairo_pattern_t* pattern_new(cairo_t *ct, Geom::OptRect const &bbox, double opacity);
 
 protected:
diff --git a/src/widgets/mesh-toolbar.cpp b/src/widgets/mesh-toolbar.cpp
index e9e9fa344..674d0f350 100644
--- a/src/widgets/mesh-toolbar.cpp
+++ b/src/widgets/mesh-toolbar.cpp
@@ -74,34 +74,119 @@ static bool blocked = false;
 //########################
 
 /*
+ * Get the current selection and dragger status from the desktop
+ */
+void ms_read_selection( Inkscape::Selection *selection,
+                        SPMeshGradient *&ms_selected,
+                        bool &ms_selected_multi,
+                        SPMeshSmooth &ms_smooth,
+                        bool &ms_smooth_multi )
+{
+
+    // Read desktop selection
+    bool first = true;
+    ms_smooth = SP_MESH_SMOOTH_NONE;
+    
+    for (GSList const* i = selection->itemList(); i; i = i->next) {
+        SPItem *item = SP_ITEM(i->data);
+        SPStyle *style = item->style;
+
+        if (style && (style->fill.isPaintserver())) {
+            SPPaintServer *server = item->style->getFillPaintServer();
+            if ( SP_IS_MESHGRADIENT(server) ) {
+
+                SPMeshGradient *gradient = SP_MESHGRADIENT(server); // ->getVector();
+                SPMeshSmooth smooth = gradient->smooth;
+                bool smooth_set = gradient->smooth_set;
+
+                if (gradient != ms_selected) {
+                    if (ms_selected) {
+                        ms_selected_multi = true;
+                    } else {
+                        ms_selected = gradient;
+                    }
+                }
+                if( smooth != ms_smooth ) {
+                    if (ms_smooth != SP_MESH_SMOOTH_NONE && !first) {
+                        ms_smooth_multi = true;
+                    } else {
+                        ms_smooth = smooth;
+                    }
+                }
+                first = false;
+            }
+        }
+
+        if (style && (style->stroke.isPaintserver())) {
+            SPPaintServer *server = item->style->getStrokePaintServer();
+            if ( SP_IS_MESHGRADIENT(server) ) {
+
+                SPMeshGradient *gradient = SP_MESHGRADIENT(server); // ->getVector();
+                SPMeshSmooth smooth = gradient->smooth;
+                bool smooth_set = gradient->smooth_set;
+
+                if (gradient != ms_selected) {
+                    if (ms_selected) {
+                        ms_selected_multi = true;
+                    } else {
+                        ms_selected = gradient;
+                    }
+                }
+                if( smooth != ms_smooth ) {
+                    if (ms_smooth != SP_MESH_SMOOTH_NONE && !first) {
+                        ms_smooth_multi = true;
+                    } else {
+                        ms_smooth = smooth;
+                    }
+                }
+                first = false;
+            }
+        }
+    }
+ }
+
+/*
  * Core function, setup all the widgets whenever something changes on the desktop
  */
-static void ms_tb_selection_changed(Inkscape::Selection * /*selection*/, gpointer /*data*/)
+static void ms_tb_selection_changed(Inkscape::Selection * /*selection*/, gpointer data)
 {
-    // DOES NOTHING AT MOMENT
 
     // std::cout << "ms_tb_selection_changed" << std::endl;
 
-    // if (blocked)
-    //     return;
-
-    // GtkWidget *widget = GTK_WIDGET(data);
+    if (blocked)
+        return;
 
-    // SPDesktop *desktop = static_cast<SPDesktop *>(g_object_get_data(G_OBJECT(widget), "desktop"));
-    // if (!desktop) {
-    //     return;
-    // }
+    GtkWidget *widget = GTK_WIDGET(data);
 
-    // Inkscape::Selection *selection = desktop->getSelection(); // take from desktop, not from args
-    // if (selection) {
-    //     ToolBase *ev = sp_desktop_event_context(desktop);
-    //     GrDrag *drag = NULL;
-    //     if (ev) {
-    //         drag = ev->get_drag();
-    //         // Hide/show handles?
-    //     }
+    SPDesktop *desktop = static_cast<SPDesktop *>(g_object_get_data(G_OBJECT(widget), "desktop"));
+    if (!desktop) {
+        return;
+    }
 
-    // }
+    Inkscape::Selection *selection = desktop->getSelection(); // take from desktop, not from args
+    if (selection) {
+        // ToolBase *ev = sp_desktop_event_context(desktop);
+        // GrDrag *drag = NULL;
+        // if (ev) {
+        //     drag = ev->get_drag();
+        //     // Hide/show handles?
+        // }
+
+        SPMeshGradient *ms_selected = 0;
+        SPMeshSmooth ms_smooth = SP_MESH_SMOOTH_NONE;
+        bool ms_selected_multi = false;
+        bool ms_smooth_multi = false; 
+        ms_read_selection( selection, ms_selected, ms_selected_multi, ms_smooth, ms_smooth_multi );
+        // std::cout << "   smooth: " << ms_smooth << std::endl;
+        
+        EgeSelectOneAction* smooth = (EgeSelectOneAction *) g_object_get_data(G_OBJECT(widget), "mesh_select_smooth_action");
+        gtk_action_set_sensitive( GTK_ACTION(smooth), (ms_selected && !ms_selected_multi) );
+        if (ms_selected) {
+            blocked = TRUE;
+            ege_select_one_action_set_active( smooth, ms_smooth );
+            blocked = FALSE;
+        }
+    }
 }
 
 
@@ -126,6 +211,33 @@ static void ms_defs_modified(SPObject * /*defs*/, guint /*flags*/, GObject *widg
     ms_tb_selection_changed(NULL, widget);
 }
 
+void ms_get_dt_selected_gradient(Inkscape::Selection *selection, SPMeshGradient *&ms_selected)
+{
+    SPMeshGradient *gradient = 0;
+
+    for (GSList const* i = selection->itemList(); i; i = i->next) {
+         SPItem *item = SP_ITEM(i->data); // get the items gradient, not the getVector() version
+         SPStyle *style = item->style;
+         SPPaintServer *server = 0;
+
+         if (style && (style->fill.isPaintserver())) {
+             server = item->style->getFillPaintServer();
+         }
+         if (style && (style->stroke.isPaintserver())) {
+             server = item->style->getStrokePaintServer();
+         }
+
+         if ( SP_IS_MESHGRADIENT(server) ) {
+             gradient = SP_MESHGRADIENT(server);
+         }
+    }
+
+    if (gradient) {
+        ms_selected = gradient;
+    }
+}
+
+
 /*
  * Callback functions for user actions
  */
@@ -178,6 +290,30 @@ static void ms_col_changed(GtkAdjustment *adj, GObject * /*tbl*/ )
     blocked = FALSE;
 }
 
+static void ms_smooth_changed(EgeSelectOneAction *act, GtkWidget *widget)
+{
+    // std::cout << "ms_smooth_changed" << std::endl;
+    if (blocked) {
+        return;
+    }
+
+    SPDesktop *desktop = static_cast<SPDesktop *>(g_object_get_data(G_OBJECT(widget), "desktop"));
+    Inkscape::Selection *selection = desktop->getSelection();
+    SPMeshGradient *gradient = 0;
+    ms_get_dt_selected_gradient(selection, gradient);
+
+    if (gradient) {
+        SPMeshSmooth smooth = (SPMeshSmooth) ege_select_one_action_get_active(act);
+        // std::cout << "   smooth: " << smooth << std::endl;
+        gradient->smooth = smooth;
+        gradient->smooth_set = true;
+        gradient->updateRepr();
+
+        DocumentUndo::done(desktop->getDocument(), SP_VERB_CONTEXT_GRADIENT,
+                   _("Set mesh smoothing"));
+    }
+}
+
 static void mesh_toolbox_watch_ec(SPDesktop* dt, Inkscape::UI::Tools::ToolBase* ec, GObject* holder);
 
 /**
@@ -322,6 +458,56 @@ void sp_mesh_toolbox_prep(SPDesktop * desktop, GtkActionGroup* mainActions, GObj
     g_object_set_data(holder, "desktop", desktop);
 
     desktop->connectEventContextChanged(sigc::bind(sigc::ptr_fun(mesh_toolbox_watch_ec), holder));
+
+    /* Warning */
+    {
+        GtkAction* act = gtk_action_new( "MeshWarningAction",
+          _("WARNING: Mesh SVG Syntax Subject to Change, Smoothing Experimental"), NULL, NULL );
+        gtk_action_group_add_action( mainActions, GTK_ACTION(act) );
+    }
+
+    /* Smoothing method */
+    {
+        GtkListStore* model = gtk_list_store_new( 2, G_TYPE_STRING, G_TYPE_INT );
+
+        GtkTreeIter iter;
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, C_("Smoothing", "None"), 1, SP_MESH_SMOOTH_NONE, -1 );
+
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, _("Default"), 1, SP_MESH_SMOOTH_SMOOTH, -1 );
+
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, _("Smooth1"), 1, SP_MESH_SMOOTH_SMOOTH1, -1 );
+
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, _("Smooth2"), 1, SP_MESH_SMOOTH_SMOOTH2, -1 );
+
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, _("Smooth3"), 1, SP_MESH_SMOOTH_SMOOTH3, -1 );
+
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, _("Smooth4"), 1, SP_MESH_SMOOTH_SMOOTH4, -1 );
+
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, _("Smooth5"), 1, SP_MESH_SMOOTH_SMOOTH5, -1 );
+
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, _("Smooth6"), 1, SP_MESH_SMOOTH_SMOOTH6, -1 );
+
+        gtk_list_store_append( model, &iter );
+        gtk_list_store_set( model, &iter, 0, _("Smooth7"), 1, SP_MESH_SMOOTH_SMOOTH7, -1 );
+
+        EgeSelectOneAction* act = ege_select_one_action_new( "MeshSmoothAction", _("None"),
+               _("If the mesh should be smoothed across patch boundaries."),
+                NULL, GTK_TREE_MODEL(model) );
+        g_object_set( act, "short_label", _("Smoothing:"), NULL );
+        ege_select_one_action_set_appearance( act, "compact" );
+        gtk_action_set_sensitive( GTK_ACTION(act), FALSE );
+        g_signal_connect( G_OBJECT(act), "changed", G_CALLBACK(ms_smooth_changed), holder );
+        gtk_action_group_add_action( mainActions, GTK_ACTION(act) );
+        g_object_set_data( holder, "mesh_select_smooth_action", act );
+    }
 }
 
 static void mesh_toolbox_watch_ec(SPDesktop* desktop, Inkscape::UI::Tools::ToolBase* ec, GObject* holder)
diff --git a/src/widgets/toolbox.cpp b/src/widgets/toolbox.cpp
index 4f5682f7a..5d52db6f2 100644
--- a/src/widgets/toolbox.cpp
+++ b/src/widgets/toolbox.cpp
@@ -521,6 +521,8 @@ static gchar const * ui_descr =
 //        "    <toolitem action='MeshEditFillAction' />"
 //        "    <toolitem action='MeshEditStrokeAction' />"
 //        "    <toolitem action='MeshShowHandlesAction' />"
+        "    <toolitem action='MeshWarningAction' />"
+        "    <toolitem action='MeshSmoothAction' />"
         "    <separator />"
         "  </toolbar>"