几何建模与处理-CVT应用于stippling

点画

接上篇《几何建模与处理-实现平面点集CVT的Lloyd算法》，这篇基于CVT的Lloyd算法，实现stippling(点画)效果。参考一篇British Columbia大学的paper，发表于2002年。Weighted Voronoi Stippling。 代码参考余畅的实现 实现平面点击CVT的Lloyd算法-余畅

代码地址github:
github.com/summer-go/G…

下面是一些stippling效果

看起来很简单，把一张灰色图离散的采样不就可以了么？

实际上，这些点要满足一些要求，比如

1. 具备随机性，否则可能会产生非预期的视觉效果
2. 随机的点要保留原图的轮廓、明暗特征

基于CVT-Lloyd算法能很好的满足这两点要求

变密度CVT

回顾下上篇中的CVT图

将vonoroi图中的每个点移到网格的中心，迭代至收敛，就能得到CVT。

stippling的实现原理和CVT过程本质是一样的。区别在于：

• CVT中点的移动是取每个网格顶点的几何中心
• stippling 以图片的灰度值作为权重来计算质心

可以理解为，将点朝着颜色更深的方向移动

原理及代码说明

核心逻辑还是CVT，这里主要讲"变密度"及图片生成的逻辑

• 基于C++ Clion实现
• 依赖库 JCash/Voronoi

主流程

int main(int argc, char** argv)
{
    jcv_rect bounding_box = {{0.0f, 0.0f}, {1.0f, 1.0f}};
    std::string density_image = argv[1];
    // 1. 加载灰度图
    CentroidVoronoiDiagram cvd( bounding_box, 100, density_image);
    
    // 2. 迭代1000次变密度CVT
    for(int i = 0; i < 1000; i++)
    {
        double max_diff = cvd.lloyd();
        if((i+1)%4 == 0)
            // 导出svg矢量图
            cvd.export_svg(argv[2]+std::to_string(int(i/4))+".svg");
        // 判断收敛，退出
        if(max_diff < 1e-7) break;
        std::cout << "\r" << "lloyd: " << i << ' '<< max_diff << std::endl ;
    }
}

变密度

看代码和注释吧，按照调用时序来说明

// 变密度lloyd
double lloyd()
{
    // 每次操作前，先释放上次的上下文内存占用
    if(diagram_.internal)
        jcv_diagram_free(&diagram_);

    // 养成好习惯，初始化内存
    memset(&diagram_, 0, sizeof(jcv_diagram));

    // 调用Jcash/Voronoi库，生成voronoi图
    jcv_diagram_generate(points_.size(), points_.data(), &bbox_, 0, &diagram_);
    // 进行CVT操作，这个过程也叫relax(放松~~形象)
    return relax_points();
}

迭代进行变密度CVT

// relax操作
double relax_points()
{
    // 获取vonoroi图，每个site代表一个网格
    const jcv_site* sites = jcv_diagram_get_sites(&diagram_);
    // 判断收敛的阈值
    double max_diff = 1e-9;
    for (int i = 0; i < diagram_.numsites; ++i)
    {
        const jcv_site* site = &sites[i];
        const jcv_point pre_p = points_[site->index];
        //  以灰度图为权重来计算，即变密度，这篇文章对应的逻辑走这里
        if (!density_.is_constant())
            // 记录每个点的新坐标
            points_[site->index] = varying_polygon_centroid(site->edges);
        else
            // is_constant表示均匀CVT，即不加权重
            points_[site->index] = polygon_centroid(site->edges);

        max_diff = std::max(max_diff, jcv_point_dist_sq(&pre_p, &points_[site->index]) );
    }
    return max_diff;
}

最重要的逻辑，单点变密度CVT计算

// 先计算每个网格的几何中心pc，即不加权重
// 以pc为原点，剖分该网格，生成多个三角形，对每个三角形进行带权重的计算质心
// 权重从输入的灰度图采样，每个三角形采样4 * 4 个点(注意判断该采样点必须在三角形内)，逻辑很像三角形光栅化时的抗锯齿操作
jcv_point CentroidVoronoiDiagram::varying_polygon_centroid(const jcv_graphedge* graph_edge,
                                                                  size_t n_sample)
{
    // 先计算每个网格的几何中心pc，即不加权重
    jcv_point pc = polygon_centroid(graph_edge);
    jcv_point center{0, 0};
    double W = 0;
    for(;graph_edge;graph_edge = graph_edge->next)
    {
        jcv_point p1 = graph_edge->pos[0], p2 = graph_edge->pos[1];
        // 以pc为原点，剖分该网格，生成多个三角形，对每个三角形进行带权重的计算质心
        Eigen::RowVector3d centroid = density_.centric(p1.x, p1.y, p2.x, p2.y, pc.x, pc.y);
        center.x += centroid.x();
        center.y += centroid.y();
        W += centroid.z();
    }
    center.x /= W;
    center.y /= W;
    return center;
}

对灰度图采样计算权重

// 注意，灰度图单个像素点灰度值越高越白，最后生成的SVG图是以black来填充，所以采样时用1 - color/255，即取反，和svg保持一致
Eigen::Vector3d ImageSampler::centric(double x1, double y1, double x2, double y2, double x3, double y3) {
    // 生成采样矩形 (max - min)/3 为采样间距
    auto [max_x, min_x] = std::minmax({x1, x2, x3});
    auto [max_y, min_y] = std::minmax({y1, y2, y3});
    double area = area_2(x1, y1, x2, y2, x3, y3);
    double dx = (max_x - min_x)/3;
    double dy = (max_y - min_y)/3;
    double integral_x = 0, integral_y = 0;
    double total_d = 0;
    for(size_t i = 0; i <= 3; i++ ) {
        for(size_t j = 0; j <= 3; j++) {
            double x = min_x + i * dx;
            double y = min_y + j * dy;
            // 计算叉乘，以判断采样点是否在三角形内部
            double alpha = area_2(x, y, x2, y2, x3, y3)/area;
            double beta = area_2(x1, y1, x, y, x3, y3)/area;
            double gamma = (1-alpha -beta);
            if( 0 <= alpha && alpha <= 1 &&  beta >= 0 && beta <= 1.0 && gamma >= 0 && gamma <= 1.0 ){
                // 对采样点的值平方，生成每个采样点的权重
                // 这里的平方我没有太理解，猜测可能是为了增加权重的差异度，使点图的轮廓更明显，类似锐化吧
                double d = std::max(0.001, std::pow(1 - this->operator()(y, x)/255, 2));
                integral_x += d * x;
                integral_y += d * y;
                total_d += d;
            }
        }
    }
    return {integral_x, integral_y, total_d};
}

图片采样算法，不是这篇文章的重点，但是也很有技巧，值得了解

// 对灰度图采样
// 实际上是采样了四个点，很经典的采样插值计算，不好解释，看代码吧
double ImageSampler::operator()(double y, double x) const
{
    // 判断边界
    x = x > 0 ? x : 0;
    x = x > 1.0 ? 1.0 : x;
    y = y > 0 ? y : 0;
    y = y > 1.0 ? 1.0 : y;
    x = (image_.cols() - 1) * x;
    y = (image_.rows() - 1) * y;
    
    size_t floor_x = std::floor(x), floor_y = std::floor(y);
    size_t ceil_x = std::ceil(x), ceil_y = std::ceil(y);

    double _x = x - floor_x, _y = y - floor_y;

    // 这儿 _x  _y 弄反了？？
    double floor_pixel = pixel(floor_y, floor_x) * (1 - _x) + pixel(ceil_y, floor_x) * _x;
    double ceil_pixel = pixel(floor_y, ceil_x) * (1 - _x) + pixel(ceil_y, ceil_x) * _x;
    return (floor_pixel * (1 - _y) + ceil_pixel * _y);
}

导出svg

svg图的格式，有兴趣可以详细了解下。stippling用黑色圆表示，坐标(cx cy)，半径(r)，填充色(black).

// 导出svg图，导出512的分辨率图
std::string  CentroidVoronoiDiagram::export_svg(std::string output_file) {
    std::stringstream svg_str;
    double w = 512;
    
    // 均值输出512 * 512， 这里有输入的灰度图，即变密度
    double h = density_.is_constant() ? 512: density_.height()/double(density_.width())*512;
    // svg 头格式
    svg_str << "<svg width=\"" << w << "\" height=\"" << h << "\" viewBox=\"0 0 " << w << ' ' << h
            << R"(" xmlns="http://www.w3.org/2000/svg" >)"
            << "<rect width=\"100%\" height=\"100%\" fill=\"white\"/>\n";
    
    // 获取voronoi图
    const jcv_site* sites = jcv_diagram_get_sites(&diagram_);
    // 如果是均匀的cvt，输出网格图边框，这次又输入的灰度图，不走这里的逻辑
    for (size_t i = 0; i < diagram_.numsites && density_.is_constant(); i++)
    {
        svg_str << "<polygon points=\"";
        jcv_graphedge* graph_edge = sites[i].edges;
        while (graph_edge)
        {
            jcv_point p1 = graph_edge->pos[0], p2 = graph_edge->pos[1];
            svg_str << p1.x*w << ',' << p1.y*h << ' ' << p2.x*w << ',' << p2.y*h << ' ';
            graph_edge = graph_edge->next;
        }
        svg_str << R"(" fill="#1d1d9b" stroke="white" stroke-width="2" />)" << "\n";
    }
    
    // 网格值心点转成svg的point值
    for( int i = 0; i < diagram_.numsites; ++i )
    {
        const jcv_site* site = &sites[i];
        jcv_point p = site->p, pc = polygon_centroid(site->edges);
        // 均匀cvt逻辑
        if(density_.is_constant()){
            svg_str << "<circle cx=\"" << pc.x * w  << "\" cy=\"" << pc.y*h << R"(" r="1" fill="red"/>)" << '\n';
            svg_str << "<circle cx=\"" << p.x * w  << "\" cy=\"" << p.y*h << R"(" r="1" fill="yellow"/>)" << '\n';
        }else
            // 变密度cvt逻辑
            svg_str << "<circle cx=\"" << p.x * w  << "\" cy=\"" << p.y*h << R"(" r="1" fill="black"/>)" << '\n';
    }

    svg_str << "</svg>";

    // 写文件
    if(output_file.size() != 0)
    {
        std::ofstream svg_file(output_file);
        svg_file << svg_str.str();
        svg_file.close();
    }
    
    return svg_str.str();
}

来个彩蛋 用老婆孩子的照片生成了一张stippling

欢迎关注公众号:sumsmile /专注图形渲染、Android开发

参考资料：

[1]Weighted Voronoi Stippling: www.cs.ubc.ca/labs/imager…

[2]实现平面点击CVT的Lloyd算法-余畅: github.com/g1n0st/GAME…

[3]JCash-voronoi库: github.com/JCash/voron…

[4]markdown 数学公式大全: blog.csdn.net/weixin_4278…

[5]点刻画(非真实渲染NPR-Weighted Voronoi Stippling): zhuanlan.zhihu.com/p/338395642

[6]markdown 数学公式大全 markdown 数学公式大全

[7]点刻画(非真实渲染NPR-Weighted Voronoi Stippling) 点刻画(非真实渲染NPR-Weighted Voronoi Stippling)