ZJT's Blog

Keep your determination.

自动脸部替换工具的构想与实现

2018-5-18 22:312018-5-20 9:13
python图像处理

最近比较有空,这几天又一直在P图,感觉换脸的过程其实就是一个枯燥的重复过程,于是就产生了做一个自动换脸工具的想法。


先来看一下最终的效果:

效果图

やりますねぇ!(赞赏)


准备工作

要实现自动P脸,肯定会用到人脸的识别。这里我直接使用了github上的一个轮子:face_recognition,这里面用到了dlib,在安装它之前要先装好cmake,接着在命令行里用pip安装即可(当然还需要用到PIL和numpy,就假装已经装好了吧)

这个库可以支持在一张图片中识别人脸,并返回每张人脸的所有特征点。后面做特征匹配的时候用到的特征点就是用它找出来的。

人脸定位

现在我们有两张人脸,每张人脸都有许多不同的特征点...

feature

我们要对其中一张人脸进行变换,使得变换后的所有特征点与另一张人脸的特征点尽可能接近。这里的变换只涉及到旋转,平移,缩放,平行扭曲(仿射变换),可以用$6$个参数$a,b,c,d,e,f$来表示:

$$x_1=ax_0+by_0+c,y_1=dx_0+ey_0+f$$

我们定义代价函数为$S=\frac 12\sum_i{({x_1}_i-{x_2}_i)^2+({y_1}_i-{y_2}_i)^2}$,其中${x_1}_i,{y_1}_i$是人脸$1$经过变换后的第$i$个特征点的坐标,${x_2}_i,{y_2}_i$是人脸$2$中第$i$个特征点的坐标。这里乘上$\frac 12$是为了后面求导方便。

我们代入变换的公式:

$$ S=\frac 12\sum(ax_0+by_0+c-x_2)^2+(dx_0+ey_0+f-y_2)^2 $$

这里省略了下标$i$。我们对六个参数求偏导:

$$ \frac{\partial S}{\partial a}=\sum(ax_0+by_0+c-x_2)x_0\\ \frac{\partial S}{\partial b}=\sum(ax_0+by_0+c-x_2)y_0\\ \frac{\partial S}{\partial c}=\sum ax_0+by_0+c-x_2\\ \frac{\partial S}{\partial d}=\sum(dx_0+ey_0+f-y_2)x_0\\ \frac{\partial S}{\partial e}=\sum(dx_0+ey_0+f-y_2)y_0\\ \frac{\partial S}{\partial f}=\sum dx_0+ey_0+f-y_2\\ $$

本来想直接梯度下降最优化这个东西的,后来发现令所有偏导数都为$0$,可以直接得到关于参数的线性方程组,直接解方程组就能得到最优的仿射变换参数了。具体方程组为:

$$ \begin{cases} (\sum x_0^2)a&+(\sum y_0x_0)b&+(\sum x_0)c&=\sum x_2x_0\\ (\sum x_0y_0)a&+(\sum y_0^2)b&+(\sum y_0)c&=\sum x_2y_0\\ (\sum x_0)a&+(\sum y_0)b&+nc&=\sum x_2\\ \end{cases} $$

($d,e,f$类似)

这样我们就能在$O(n)$的时间内求出对应的仿射变换了(虽然这里对性能并没有什么要求。。因为特征点并不多)。找到仿射变换后,可以用PIL.Image.transform把变换作用到图片上。

为了解决面部中心特征点过于密集的问题,我们还可以让特征点带上权值,降低较密集的特征点的权值。具体可以看代码。

我们把所有经过对应变换的人脸叠加在一张新图中,后面需要把这张图和原图进行合成。

图片合成

现在我们有两张图片:底图(记为图1)和一张由新的人脸组成的图(记为图2)。

一个比较直接的方法是直接把图2按alpha通道混合在图1中(这意味着什么都没有处理)。这样的话新图的颜色可能会跟原图对不上,人脸的边缘也会比较粗糙。

为了解决这两个问题,我们做两个额外的步骤:颜色调整和边缘羽化。

为了让颜色更加自然,我们需要调整图2中的颜色,让其接近图1的同时尽可能地保留细节。解决方案是把图1对图2遮罩得到图3(即削除了图1中不含人脸的部分),把图2和图3分别高斯模糊得到图4和图5,把图2对应像素除以图4再乘以图5,就能得到这个像素调整后的颜色。这么做是因为高斯模糊消除了图片细节,所以这样调整后即不会表现出覆盖掉的内容中的细节,也能体现其原有的大体颜色。

处理完颜色后,我们还需要对边缘进行羽化,让过渡更加自然。我没想到什么特别好的方法,直接对图2的alpha通道做高斯模糊后,再令新的alpha值$\alpha'$为$\min(0,\max(1,(\alpha-0.5)\times 2))$。这样做的原理是把让边缘的$alpha$值取$0$,且往内部扩展时$alpha$值变大,在边缘平滑的时候效果比较好,因为平滑边缘处的高斯模糊值可近似看作$0.5$,但在边缘鬼畜的时候可能会有一些奇妙问题。。不过一般情况下效果还是挺不错的。


具体实现

https://github.com/Alif-01/AutoHeadReplace

求一波star QwQ

神秘链接

查看详细内容

FFT相关的一些优化技巧

2018-4-17 17:132018-4-17 17:13
FFT卡常数

大概就是myy论文里的那些东西吧,闲着没事整理一下。

记FFT变换的矩阵为$F$(即单位根的范德蒙德矩阵),将FFT视作$F$乘上待变换的行向量:$FFT(a)=Fa$。我们记$\overline F$为$F$的复共轭矩阵(矩阵的每个元素都共轭),$a'$为$a$除了第$0$项外都翻转后的结果(本文的所有矩阵、向量的下标均从$0$开始),即第$a_i$与$a_{n-i}$交换($n$为$a$的长度)。显然$a'$可以写成$a$乘上一个翻转矩阵$T$,即$a'=aT$。

  • 性质 1:$\overline Fa=(Fa)'$

这个根据单位根的性质$\omega^{-k}=\omega^{n-k}$可以轻易得到。

  • 性质 2:$F^{-1}=\overline F/n$

这其实就是IFFT最后一步干的事情。。

技巧 1:快速求两个实多项式的点值(以及插值)

现在有两个实向量$a,b$,我们想计算$Fa,Fb$。按照一般的做法,直接跑两遍FFT就行了。其实这个过程可以做到只用一次FFT。

我们构造复向量$a+bi$,考虑把这个向量FFT: $$F(a+bi)=Fa+F(bi)=Fa+iFb$$

我们把这整个东西取复共轭,看看会发生什么: $$ \begin{aligned} \overline{F(a+bi)}&=\overline {Fa}+\overline {iFb}\\ &=\overline Fa-i\overline Fb\\ &=(Fa)'-i(Fb)'\\ \end{aligned} $$

两边同时乘上$T$,再联立之前FFT的结果: $$ Fa-iFb=(\overline{F(a+bi)})'\\ Fa+iFb=F(a+bi) $$

知道这两个东西之后,可以直接解出$Fa,Fb$。也就是说,我们直接算出$a+bi$的FFT,把它复共轭之后翻转系数(从第$1$项开始翻转),之后就可以解出$Fa$和$Fb$了。

同理,我们也可以同时做两个实多项式点值的插值。把这整个过程反过来做一遍即可。

技巧 2:4次FFT的任意模数多项式乘法

这其实是上面那个东西的一个应用。现在有两个整系数多项式$A(x),B(x)$,我们要算它们循环卷积的结果(如果要做普通卷积,直接把次数扩大一倍即可),所有系数都对一个$10^9$级别的整数$P$取模(这个整数可以什么性质都没有)。我们考虑选取整数$M\approx\sqrt P$,我们把$A,B$的系数对$M$作带余除法,得到$A(x)=A_0(x)+MA_1(x),B(x)=B_0(x)+MB_1(x)$。这里$A_0,A_1,B_0,B_1$的系数都是$M$级别的,他们之间乘起来之后系数是$O(M^2n)=O(Pn)$级别的,$n$不太大的话(大概$5*10^5$以内)可以认为是在double的精度范围内,可以直接用double做乘法后取模。

我们用上面提到的那个技巧,用$2$次FFT算出$A_0,A_1,B_0,B_1$的FFT点值$a_0,a_1,b_0,b_1$,之后根据$$\begin{aligned}A(x)B(x)&=A_0(x)B_0(x)\\&+M(A_0(x)B_1(x)+A_1(x)B_0(x))\\&+M^2A_1(x)B_1(x)\end{aligned}$$我们可以再次使用上面那个技巧,用$2$次FFT算出$a_0\otimes b_0,a_0\otimes b_1,a_1\otimes b_0,a_1\otimes b_1$的IFFT插值,就可以直接进行计算了。

总共$4$次FFT。

技巧 3:2次FFT的实多项式乘法

这个用技巧$1$也可以做到,但较技巧$1$更容易实现。
有两个实多项式$A(x),B(x)$,我们要计算它们的循环卷积。

我们可以用两次FFT计算$(A(x)+iB(x))^2=A^2(x)-B^2(x)+2iA(x)B(x)$。

直接取虚部并除以2即可。

查看详细内容

动态图连通性问题的在线解法

2018-4-16 22:252018-4-16 22:25
动态图Euler-tour Tree毒瘤

最近居然在省选模拟赛里面做到了这个毒瘤问题。。

动态图连通性问题是指维护一个不断加边删边的无向图,要求查询两个点之间的连通性,以及支持对连通块的各种信息的维护。离线的话非常好做,直接处理出每条边的删除时间后,维护关于删除时间的最大生成树即可。强制在线的话,也有一个$O(\log^2n)$的做法,是通过将图分层并用ETT维护每一层的生成树来解决的。不过这个算法实现起来比较复杂,第一次写完+调试完大概花了6h,估计场上想写是不太可能的了。。

为了维护连通性的信息,我们可以考虑维护生成森林。我们给每条边钦点一个非负整数的等级(level),定义$G_i$为所有等级为$i$的边构成的子图,$T_i$为所有等级大于等于$i$的边构成的图中,关于等级的最大生成森林。最终询问的时候,直接在$T_0$上做询问就行了。

当一条新边加入时,我们先令它的等级为$0$。这时,我们会在$G_0$中加入这条边,如果$T_0$中加入这条边不会形成环,我们也把它加入$T_0$。

删除的情况比较复杂。如果删掉的是一条不在$T_0$中的边,那它也不可能在任何$T_i$中出现,我们直接在对应等级的$G_i$中把它删掉就行了。然而如果我们要删掉一条树边,则可能存在一条非树边,使得删掉这条树边形成两个连通块可以被这条非树边连接。这时候这条非树边就会取代删去的边,成为新的树边。

考虑如何在图中找到这条非树边。我们记删去的边为$e_0$,等级为$l_0$,替换它的非树边为$e_1$,等级为$l_1$。显然$l_1\leq l_0$,否则最大生成森林的性质就不满足了($e_1$可以替代$e_0$得到一个更大的生成森林)。

我们考虑从$l_0$到$0$枚举$l_1$,看看是否存在等级为$l_1$的非树边$e_1$能替换删去的$e_0$。一个比较显然的结论是$e_1$连接的两个端点一定在$T_{l_1}$上与$e_0$的端点连通,这个可以通过MST的性质得到。记$e_0$在$T_{l_1}$上连接的两个连通块为$X,Y$,且$|X|\leq|Y|$,考虑一个暴力做法,我们直接在$G_{l_1}$中枚举每条从$X$连出去的边,如果连出去的边位于$Y$中,则找到了$e_1$。如果枚举完每条边之后还是没找到,则说明不存在等级为$l_1$的满足条件的$e_1$。

这样做的复杂度显然是错的,我们考虑对它进行优化。我们在暴力枚举每条边的时候 ,如果这条边不能满足条件,就把它的等级$+1$(在$G_{l_1}$中删去它,在$G_{l_1+1}$中加入它,如果它是树边还要在$T_{l_1+1}$中也加入它)。这样每次访问到一条边时,要么找到了答案,要么这条边的等级增加,所以我们暴力枚举的复杂度就只跟每条边的最大等级有关了(等级增加次数肯定不会超过最大等级)。同时,我们还可以证明,$T_i$中的最大连通块大小是$O(n/2^i)$级别的。因为我们每次都是只处理从较小的连通块中连出的边($|X|\leq|Y|$),且在我们把所有边升级完之后,下一等级的生成树上形成的连通块大小也不会超过$|X|$(因为当前层的生成树一定包含下一层的生成树),所以每一层的连通块大小的上界都会减半。这也说明了边的最大等级是$O(\log n)$级别的,即每条边都最多只会被升级$O(\log n)$次。

我们用ETT维护$T_i$中有哪些点有$G_i$的连边,这样每次搜索的复杂度就是$O($搜到的边数$\times\log n)$的,在ETT中的link和cut也是$O(\log n)$的,所以总复杂度为$O(m\log^2 n)$。这样就成功解决了这个问题。

实现上可能有一些细节问题要注意,首先在搜索$X$的出边时,要先搜所有的树边。如果先搜了非树边,可能会因为非树边等级的提高而出现生成树心态发生改变的情况,这种情况是不好处理的。其次对边的搜索要一边搜索一边检查是否可行,如果先搜出了所有边再一个个判断,复杂度就是错的,在loj上会被卡(

这是模板题

附上代码实现一份(loj垫底qwq):

#include <bits/stdc++.h>
#define MAXN 500010
#define MAXL 20
using namespace std;

int read(){
    char c;
    while((c=getchar())<'0' || c>'9');
    int res=c-'0';
    while((c=getchar())>='0' && c<='9') res=res*10+c-'0';
    return res;
}

int n,m,online;
int v[MAXN];

namespace Duliu{
    int ls[MAXN],numls;

    struct ETT{
        struct node{
            int lc,rc,p,key,size,pt;
            bool tag,stag;
        }a[MAXN<<1];

        int blist[MAXN],numb;
        int p[MAXN];
        map<int,int> e[MAXN];
        int tagp[MAXN];

        ETT(int n){
            for(int i=2*n;i>=1;i--) a[blist[++numb]=i].key=rand();
        }

        int newnode(int pt){
            int x=blist[numb--];
            a[x].p=a[x].lc=a[x].rc=a[x].tag=a[x].stag=0;
            a[x].size=1;
            a[x].pt=pt;
            return x;
        }

        void delnode(int x){
            blist[++numb]=x;
        }

        inline void pushup(int x){
            a[x].stag=a[a[x].lc].stag||a[a[x].rc].stag||a[x].tag;
            a[x].size=a[a[x].lc].size+a[a[x].rc].size+1;
        }

        int merge(int x,int y){
            if(!x || !y) return x+y;
            if(a[x].key<a[y].key){
                a[x].rc=merge(a[x].rc,y);
                a[a[x].rc].p=x;
                pushup(x);
                return x;
            }else{
                a[y].lc=merge(x,a[y].lc);
                a[a[y].lc].p=y;
                pushup(y);
                return y;
            }
        }

        pair<int,int> split(int x,int d){
            int s1=a[x].lc,s2=a[x].rc;
            if(d==-1){
                a[x].rc=0;
                pushup(x);
                s1=x;
                a[s2].p=0;
            }else if(d==1){
                a[x].lc=0;
                pushup(x);
                s2=x;
                a[s1].p=0;
            }else a[s1].p=a[s2].p=0;
            while(a[x].p){
                int p=a[x].p;
                if(a[p].lc==x){
                    a[p].lc=s2;
                    a[s2].p=p;
                    a[s1].p=0;
                    s2=p;
                }else{
                    a[p].rc=s1;
                    a[s1].p=p;
                    a[s2].p=0;
                    s1=p;
                }
                pushup(x=p);
            }
            return make_pair(s1,s2);
        }

        int getRoot(int x){
            while(a[x].p) x=a[x].p;
            return x;
        }

        int get_p(int x){
            if(e[x].empty()) return 0;
            return e[e[x].begin()->first][x];
        }

        void link(int x,int y){
            int p1=newnode(y),p2=newnode(x);
            if(tagp[x]==-1){
                tagp[x]=p2;
                a[p2].tag=a[p2].stag=1;
            }
            if(tagp[y]==-1){
                tagp[y]=p1;
                a[p1].tag=a[p1].stag=1;
            }
            int t1,t2,t3;
            if(e[x].empty()) t1=t3=0;
            else{
                pair<int,int> temp=split(e[x].begin()->second,1);
                t1=temp.first; t3=temp.second;
            }
            if(e[y].empty()) t2=0;
            else{
                pair<int,int> temp=split(e[y].begin()->second,1);
                t2=merge(temp.second,temp.first);
            }
            e[x][y]=p1; e[y][x]=p2;
            merge(merge(merge(t1,p1),t2),merge(p2,t3));
        }

        void cut(int x,int y){
            int p1=e[x][y],p2=e[y][x];
            e[x].erase(e[x].find(y));
            e[y].erase(e[y].find(x));
            pair<int,int> t1=split(p1,0);
            if(getRoot(p2)==t1.first){
                pair<int,int> t2=split(p2,0);
                merge(t2.first,t1.second);
            }else{
                pair<int,int> t2=split(p2,0);
                merge(t1.first,t2.second);
            }
            if(a[p1].tag){
                if(e[y].empty()) tagp[y]=-1;
                else{
                    tagp[y]=get_p(y);
                    a[tagp[y]].tag=1;
                    for(int i=tagp[y];i;i=a[i].p) pushup(i);
                }
            }
            if(a[p2].tag){
                if(e[x].empty()) tagp[x]=-1;
                else{
                    tagp[x]=get_p(x);
                    a[tagp[x]].tag=1;
                    for(int i=tagp[x];i;i=a[i].p) pushup(i);
                }
            }
            delnode(p1); delnode(p2);
        }

        void setTag(int x,int v){
            if(v==1){
                int p=get_p(x);
                if(!p){
                    tagp[x]=-1;
                    return;
                }
                x=tagp[x]=p;
            }else{
                int t=tagp[x];
                tagp[x]=0;
                x=t;
                if(x==-1){
                    tagp[x]=0;
                    return;
                }
            }
            a[x].tag=v;
            while(x){
                pushup(x);
                x=a[x].p;
            }
        }

        int get_size(int x){
            if(!x) return 1;
            return a[getRoot(x)].size/2+1;
        }

        void dfs1(int x,int i);
        bool dfs2(int x,int i,int py);
    }*s[MAXL];

    int l;
    map<int,int> gl[MAXN];
    set<int> g[MAXL][MAXN];
    pair<int,int> res_e;

    void init(int n){
        for(l=0;(1<<l)<n;l++);
        for(int i=0;i<=l;i++) s[i]=new ETT(n);
    }

    bool link(int x,int y){
        gl[x][y]=0;
        gl[y][x]=0;
        bool flagx=g[0][x].empty(),flagy=g[0][y].empty();
        g[0][x].insert(y);
        g[0][y].insert(x);
        int px=s[0]->get_p(x),py=s[0]->get_p(y);
        if(!px || !py || s[0]->getRoot(px)!=s[0]->getRoot(py)){
            s[0]->link(x,y);
            if(flagx) s[0]->setTag(x,1);
            if(flagy) s[0]->setTag(y,1);
            return 1;
        }
        if(flagx) s[0]->setTag(x,1);
        if(flagy) s[0]->setTag(y,1);
        return 0;
    }

    void level_up(int x,int y,int curl){
        if(g[curl][x].empty()) s[curl]->setTag(x,0);
        if(g[curl][y].empty()) s[curl]->setTag(y,0);
        gl[x][y]++;
        gl[y][x]++;
        bool flagx=g[curl+1][x].empty(),flagy=g[curl+1][y].empty();
        g[curl+1][x].insert(y);
        g[curl+1][y].insert(x);
        int px=s[curl+1]->get_p(x),py=s[curl+1]->get_p(y);
        if(!px || !py || s[curl+1]->getRoot(px)!=s[curl+1]->getRoot(py)) s[curl+1]->link(x,y);
        if(flagx) s[curl+1]->setTag(x,1);
        if(flagy) s[curl+1]->setTag(y,1);
    }

    void gao_1(int x,int i){
        auto it=g[i][x].begin();
        while(it!=g[i][x].end()){
            if(s[i]->e[x].count(*it)){
                int y=*it;
                g[i][y].erase(g[i][y].find(x));
                it=g[i][x].erase(it);
                level_up(x,y,i);
            }else it++;
        }
    }

    bool gao_2(int x,int i,int py){
        auto it=g[i][x].begin();
        while(it!=g[i][x].end()){
            if(s[i]->getRoot(s[i]->get_p(*it))==py){
                for(int k=0;k<=i;k++)
                    s[k]->link(x,*it);
                res_e=make_pair(x,*it);
                return 1;
            }else{
                int y=*it;
                g[i][y].erase(g[i][y].find(x));
                it=g[i][x].erase(it);
                level_up(x,y,i);
            }
        }
        return 0;
    }

    pair<int,int> cut(int x,int y){
        int curl=gl[x][y];
        gl[x].erase(gl[x].find(y));
        gl[y].erase(gl[y].find(x));
        g[curl][x].erase(g[curl][x].find(y));
        g[curl][y].erase(g[curl][y].find(x));
        if(s[0]->e[x].count(y)){
            for(int i=0;i<=curl;i++) s[i]->cut(x,y);
            for(int i=curl;i>=0;i--){
                int px=s[i]->get_p(x),py=s[i]->get_p(y);
                int sx=s[i]->get_size(px),sy=s[i]->get_size(py);
                if(sy<sx) swap(x,y),swap(sx,sy),swap(px,py);
                py=s[i]->getRoot(py);
                px=s[i]->getRoot(px);
                if(px){
                    s[i]->dfs1(px,i);
                    if(s[i]->dfs2(px,i,py)) return res_e;
                }else{
                    gao_1(x,i);
                    if(gao_2(x,i,py)) return res_e;
                }
            }
            return make_pair(0,0);
        }
        return make_pair(-1,0);
    }

    bool check(){
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++){
                bool flag=0;
                for(int k=l;k>=0;k--)
                    if(s[k]->e[i].count(j)) flag=1;
                    else if(flag) return 0;
            }
        return 1;
    }

    void ETT::dfs1(int x,int i){
        if(!a[x].stag) return;
        if(a[x].tag) gao_1(a[x].pt,i);
        if(a[a[x].lc].stag) dfs1(a[x].lc,i);
        if(a[a[x].rc].stag) dfs1(a[x].rc,i);
    }

    bool ETT::dfs2(int x,int i,int py){
        if(!a[x].stag) return 0;
        if(a[x].tag) if(gao_2(a[x].pt,i,py)) return 1;
        return dfs2(a[x].lc,i,py) || dfs2(a[x].rc,i,py);
    }
}

void link(int x,int y){ Duliu::link(x,y); }
void cut(int x,int y){ Duliu::cut(x,y); }

bool query(int x,int y){
    using namespace Duliu;
    if(x==y) return 1;
    int px=s[0]->get_p(x),py=s[0]->get_p(y);
    return px && py && s[0]->getRoot(px)==s[0]->getRoot(py);
}

int main(){
#ifdef DEBUG
    freopen("122.in","r",stdin);
#endif
    n=read(); m=read(); online=1;
    Duliu::init(n);
    int lastans=0;
    while(m--){
        int t=read(),x=read(),y=read();
        if(online){
            x^=lastans;
            y^=lastans;
        }
        if(t==0) link(x,y);
        if(t==1) cut(x,y);
        if(t==2){
            if(query(x,y)){
                lastans=x;
                puts("Y");
            }else{
                lastans=y;
                puts("N");
            }
        }
    }
    return 0;
}

查看详细内容

一般图最大匹配的一些细节

2018-3-15 13:02018-3-15 15:53
一般图匹配带花树

两年前就听hgr说了这个东西,一直打算找时间学习一个,结果一直鸽到了现在。。

我们都知道一般图最大匹配可以用带花树算法来做。关于带花树算法,网上有很多非常不错的学习资源(比如这个),所以算法的大体思路这里就不再赘述了。

这个算法的细节比较多,没有一个比较深刻的认识的话,感觉好像很容易写错的样子。为了明确细节,我们先来强调一下算法的基本概念。在这个算法当中,我们有两个概念(这个概念是我自己瞎脑补的,所以听起来可能比较扯淡),一个是“原交错树”,一个是“缩花后的交错树”。原交错树就是未经过缩花的交错树,我们bfs扩展的就是这棵交错树。在这棵交错树上,是可能会出现两个偶点相邻的情况的。为了避免这种情况,我们还需要另一棵“缩花后的交错树”,原交错树中的一个花,在这棵交错树上就会被缩成一个点,缩花后的点可以用这朵花的花托来表示。这样的话,就不会出现两个偶点相邻,就可以跟二分图匹配一样去找增广路了。具体实现上,缩花后的交错树中的每个节点,都对应了原交错树中的一个点集,可以用并查集来维护。

我们知道,同一朵花的每个点,不管是奇点还是偶点,都可以通过绕花的方式找到一条到花托的交错路。所以,在找交错路径的时候,可以把这些点都看作同一个偶点,这也就是“缩花”的本质。对于一个点,它在并查集中的根,表示了这个点通过绕花能绕到的最高点(最高指的是在交错树中的深度最小)。同时,对于一个被缩起来的奇点,我们可以在原交错树上记录这个点应该如何绕花,也就是记下把对应花的横叉边。这样,我们只要利用记下的信息,就能找到原交错树中的某个点到根节点的交错路径。

现在,这些细节就很明确了:所有跟缩花以及寻找增广路有关的操作,都应该在使用并查集的缩花后的交错树上进行。而发现增广路后,我们需要在原交错树上沿着增广路翻转匹配边,这里可以利用之前记录的花的信息来处理。注意后一步操作是跟缩花后的交错树无关的,所以不需要用到并查集。

另外,在广搜的过程中,队列中应该只记录待扩展的偶点。对于那些在缩花过程中被缩起来的奇点,他们也可以被看成偶点,所以缩花时也要把缩起来的所有奇点加入队列。关于如何找交错树上的LCA,我的做法是类似树链剖分那样,记录每个点在原交错树中的深度,每次选取对应花托较深的那个点往上跳,直到两个点跳到了同一个花内。这样写感觉细节不多,比沿路打标记的做法要好写一点。

附上代码实现一份。由于写法奇特,在uoj上垫底了。。

#include <bits/stdc++.h>
#define MAXN 3010
using namespace std;

struct edge{
    int to,next;
    edge(int _to=0,int _next=0):to(_to),next(_next){}
}e[1000000];

int n,m;
int g[MAXN],nume;
//访问状态,0:未访问,1:奇点,2:偶点
int visit[MAXN];
int ans;
int nump;
//p表示每个点匹配点(0表示未匹配),pre表示每个点在原交错树上的父亲
int p[MAXN],pre[MAXN];
//c1,c2记录的是每个奇点对应的花的横叉边的两个端点
int c1[MAXN],c2[MAXN];
//每个点在原交错树上的深度
int dep[MAXN];
static queue<int> Q;

//并查集
int f[MAXN];
int getf(int x){ if(f[x]==x) return x; return f[x]=getf(f[x]); } 
void merge(int x,int y){ f[getf(x)]=getf(y); }

void addEdge(int u,int v){
    e[nume]=edge(v,g[u]);
    g[u]=nume++;
    e[nume]=edge(u,g[v]);
    g[v]=nume++;
}

int getLCA(int x,int y){
    while(getf(x)^getf(y)){
        if(dep[f[x]]<dep[f[y]]) x^=y^=x^=y;
        x=pre[f[x]];
    }
    return getf(x);
}

//把x到y的交错路径翻转(为了好写用了递归,所以比较慢...)
void go(int x,int y,int d){
    //GO is GOD
    if(x==y) return;
    if(!d){
        p[x]=pre[x];
        p[pre[x]]=x;
        go(pre[x],y,1);
    }else{
        if(visit[x]==2) go(pre[x],y,0);
        else{
            //绕花
            int t1=c1[x],t2=c2[x];
            p[t1]=t2;
            p[t2]=t1;
            go(t1,x,1);
            go(t2,y,1);
        }
    }
}

//缩花,并把奇点加入队列
void gao_blossom(int x,int y,int lca){
    int _x=x,_y=y;
    x=getf(x);
    while(x^lca){
        int t=pre[x];
        Q.push(t);
        c1[t]=_x; c2[t]=_y;
        int t2=getf(pre[t]);
        merge(x,t);
        merge(t,t2);
        x=t2;
    }
}

bool gao(int rt){
    for(int i=1;i<=n;i++) visit[i]=0,f[i]=i;
    while(!Q.empty()) Q.pop();
    Q.push(rt);
    visit[rt]=2;
    dep[rt]=0;
    while(!Q.empty()){
        int x=Q.front();
        Q.pop();
        for(int i=g[x];~i;i=e[i].next){
            int y=e[i].to;
            if(getf(x)==getf(y)) continue;
            if(!visit[y]){
                visit[y]=1;
                if(!p[y]){
                    //找到增广路
                    pre[y]=x;
                    go(y,rt,0);
                    return 1;
                }else{
                    //扩展交错树
                    int z=p[y];
                    visit[z]=2;
                    pre[y]=x;
                    pre[z]=y;
                    dep[y]=dep[x]+1;
                    dep[z]=dep[y]+1;
                    Q.push(z);
                }
            }else if(visit[getf(y)]==2){
                //发现花
                int lca=getLCA(getf(x),getf(y));
                gao_blossom(x,y,lca);
                gao_blossom(y,x,lca);
            }
        }
    }
    return 0;
}

int main(){
    memset(g,-1,sizeof g);
    ans=0;
    scanf("%d%d",&n,&m);
    for(int i=1;i<=m;i++){
        int u,v;
        scanf("%d%d",&u,&v);
        addEdge(u,v);
    }
    for(int i=1;i<=n;i++)
        if(!p[i] && gao(i))
            ans++;
    printf("%d\n",ans);
    for(int i=1;i<=n;i++) printf("%d ",p[i]);
    return 0;
}

查看详细内容

扩展拉格朗日反演

2018-3-4 9:432018-4-6 16:8
拉格朗日反演FFT生成函数

突然发现拉格朗日反演其实还有扩展的形式。。。

定义$F^{-1}$为$F$的复合逆:

$$ F^{-1}(F(x))=x $$

$$[x^n]G(F^{-1}(x))=\frac1n[x^{-1}] ( \frac{dG(x)}{dx}\frac1{F^n(x)})$$

(这个好像叫Lagrange–Bürmann formula

注意一般$G$的表达式都同时含有$x$和$F^{-1}(x)$,我们把$F(x)$代入$x$后就能得到$G(x)$的表达式,进而求出$\frac {dG(x)}{dx}$。

感觉还是蛮有用的(大雾)

查看详细内容