[题解]-CF1626E Black and White Tree

Posted on 2022-01-19 Edited on 2022-02-15 In 题解

题意

给出一个大小为$n$ $(n\leq 3 \cdot 10^5)$的树, 树上每个节点被染成白色或者黑色. 假设当前在点$k$有一个棋子, 你每次可以选择一个黑点$b$作为目标点, 将这个棋子向黑点按从$k$到$b$的最短路径移动一格(一条边). 连续两次操作选择的黑点不能相同. 对于所有的$k = 1...n$, 判断该棋子是否能在有限步操作后移动到任意一个黑点. 能移动到则该点答案为$1$否则为$0$.

题解

首先显然黑点和与黑点相邻的白点答案均为$1$. 对于其它白点, 当且仅当它所在的某条链上有至少两个黑点时它才能移动到黑点上. 这种情况有一个特例:

如图, 尽管1不在任意一条有两个黑点的链上, 但它仍然能到达黑点5, 只要交替选择5, 7两个点即可.

为了解决这种特殊情况, 我们把黑点及与黑点相邻的白点缩成一个点(如果有两个相邻的黑点那么显然答案全为1. 如果不进行特判的话就需要把在同一联通分量里的黑点缩成一个点)

缩点之后如图. 这样1便在有两个黑点的链上了.

判断每个点是否在有不少于两个黑点的链上可以用换根dp求. 首先指定一个根$r$, $siz[u]$表示以$u$为根的子树里拥有黑点数最大的链的黑点数目. 换根dp的具体细节见代码中dfs2(). 需要注意的是缩点之后可能不存在$1$这个点.

时间复杂度$O(n)$.

代码

//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
mt19937 Rand(0);
const int maxn = 3e5+10;//CHANGE IT!
const int maxm = maxn*2;
const int p = 1e9+7;//CHANGE IT!
ll Pow(ll x,ll d,ll p){ll ta=1;if(d==0)return 1%p;x %= p;ll a=Pow(x,d/2,p);ta=a*a%p;if(d%2)ta=ta*x%p;return ta%p;}
vector<int>G[maxn],G2[maxn];
int cli[maxn],w[maxn],ans[maxn],siz[maxn];
void pre(int u,int fa,int r){
	cli[u] = r;
	ans[u] = 1;
	if(r != u)w[r] += w[u];
	for(int v:G2[u]){
		if(v == fa)continue;
		if(w[u] and !ans[v])pre(v,u,r);
	}
}
void dfs(int u,int fa){
	for(int v:G[u]){
		if(v == fa)continue;
		dfs(v,u);
		siz[u] = max(siz[u],siz[v]);
	}
	siz[u] += w[u];
}
void dfs2(int u,int fa,int res){
	int tot = 0,mx1 = 0,mx2 = 0;
	for(int v:G[u]){
		if(v == fa)continue;
		tot = max(tot,siz[v]);
		if(siz[v] >= mx1){
			mx2 = mx1;
			mx1 = siz[v];
		}
		else if(siz[v] > mx2)mx2 = siz[v];
	}
	tot = max(tot,res) + w[u];
	if(res >= mx1){
		mx2 = mx1;
		mx1 = res;
	}
	else if(res > mx2)mx2 = res;
	if(tot > 1)ans[u] = 1;
	for(int v:G[u]){
		if(v == fa)continue;
		int rw = w[u] + (mx1 == siz[v]?mx2:mx1);
		dfs2(v,u,rw);
	}
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	int n;
	cin >> n;
	for(int i = 1;i <= n;i++)cin >> w[i],cli[i] = i;
	for(int i = 1;i < n;i++){
		int f,t;
		cin >> f >> t;
		G2[f].push_back(t);
		G2[t].push_back(f);
	}
	for(int i = 1;i <= n;i++)if(w[i] and cli[i] == i)pre(i,-1,i);
	for(int u = 1;u <= n;u++){
		for(int v:G2[u]){
			if(cli[u] == cli[v])continue;
			G[cli[u]].push_back(cli[v]);
		}
	}
	for(int i = 1;i <= n;i++){
		if(cli[i] == i){
			dfs(i,-1);
			dfs2(i,-1,0);
			break;
		}
	}
	for(int i = 1;i <= n;i++)cout << ans[i] << " ";
	return 0;
}

[题解]CF1625D-Binary Spiders

Posted on 2022-01-13 In 题解

题意

给出$n$个数$a_i$, 和一个整数$k$. 从这$n$个数中选出尽可能多的数, 使得它们之间任意两者按位异或后大于等于$k$. 输出数字的个数和下标, 如果不存在则输出$-1$.

$2 \le n \le 3 \cdot 10^5, 0 \leq a_i,k \leq 2^{30} - 1$

题解

设$k$的最高位是第$h$位. 首先会很自然的想到——如果把每个数位数小于等于$h$的部分变为0, 按此分类, 从每种里选出一个数字加入答案里一定合法. 因为任意两个数字异或后都至少有一个高于$h$的位为$1$. 以样例$k = 8(1000)$为例:

原始	2(10)	8(1000)	4(100)	16(10000)	10(1010)	14(1110)
处理后组别	0	0	0	16(10000)	0	0

此时我们选择$16$, 再从$2,8,4,10,14$里任选一个都能组成一个解，但它不是最优的: 例如在$0$这一组里选择$2$和$10$加入答案也是合法解.

进一步思考可以发现, 在每一组里有可能可以选择两个数字且最多只能选择两个. 如果选择了超过两个数字, 那么至少有两个数字在第$h$位上是一样的, 这会导致异或后的结果小于$k$. 使用Trie判断是否有解即可.

注意需要特判$k = 0$的情况.

代码

//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
const int maxn = 3e5+10;//CHANGE IT!
const int maxm = maxn*2;
const int maxc = 1e7+10;
const int p = 1e9+7;//CHANGE IT!
ll Pow(ll x,ll d,ll p){ll ta=1;if(d==0)return 1%p;x %= p;ll a=Pow(x,d/2,p);ta=a*a%p;if(d%2)ta=ta*x%p;return ta%p;}
ll raw[maxn],u[maxn],siz;
int nxt[maxc][2];
int w[maxc*2];
map<int,vector<pair<int,int>> >hsh;
void insert(int s,int idx){
	int u = 0;
	for(int i = 30;i >= 0;i--){
		int c = ((s>>i)&1);
		if(nxt[u][c] == -1)nxt[u][c] = ++siz;
		u = nxt[u][c];
	}
	w[u] = idx;
}
int query(int x,int k){
	int u = 0,res = 0;
	for(int i = 30;i >= 0;i--){
		int c = ((x>>i)&1);
		if(nxt[u][c^1] != -1){
			res += (1<<i);
			u = nxt[u][c^1];
		}
		else if(nxt[u][c] != -1){
			u = nxt[u][c];
		}
		else return -1;
	}
	if(res < k)return -1;
	else return w[u];
}
void clear(){
	for(int i = 0;i <= siz;i++){
		nxt[i][0] = nxt[i][1] = -1;
		w[i] = 0;
	}
	siz = 0;
}
signed main(){
	memset(nxt,-1,sizeof(nxt));
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	int n,k;
	cin >> n >> k;
	if(!k){
		cout << n << endl;
		for(int i = 1;i <= n;i++)cout << i << " ";
		return 0;
	}
	int h = 0;
	for(int i = 30;i >= 0;i--)if((1<<i)&k){
		h = i;
		break;
	}
	set<int>great,res;
	vector<int>gidx;
	for(int i = 1;i <= n;i++){
		cin >> raw[i];
		u[i] = raw[i];
		for(int j = h;j >= 0;j--)if((1<<j)&u[i])u[i] -= (1<<j);
		cerr << u[i] << " ";
		hsh[u[i]].push_back({raw[i],i});
	}
	vector<int>ans;
	for(auto it:hsh){
		vector<pair<int,int> > & r = it.second;
		clear();
		for(auto it2:r)insert(it2.first,it2.second);
		bool ok = 0;
		for(auto it2:r){
			int x = query(it2.first,k);
			if(x != -1){
				ans.push_back(it2.second);
				ans.push_back(x);
				ok = 1;
				break;
			}
		}
		if(!ok and r.size())ans.push_back(r[0].second);
	}
	if(ans.size() < 2){
		cout << -1;
	}
	else{
		cout << ans.size() << endl;
		for(int x:ans)cout << x << " ";
	}
	return 0;
}

[题解]CF1610D-Not Quite Lee

Posted on 2021-11-24 Edited on 2022-01-13 In -题解

经典我是傻逼环节.

题意

给出长度为$n$的数组$A = \{a_i\}$, $1 \leq n \leq 2 \cdot 10^5,1 \leq a_i \leq 10^9$.

对于$A$的任意长度为$m$的非空子序列$B = \{b_1,b_2,..,b_m\}$, 定义$f(b_i,k_i) = \sum_{j = 0}^{b_i} (j + k_i)$, 其中$k_i$为任意整数. 一个子序列$B$合法当且仅当 $ {k_i} |_{i = 1}^mf(b_i,k_i) = 0$.

求$A$中全部合法子序列数量模$10^9+7$.

题解

首先我们考虑如何判断某个序列是否合法. 有$f(b_i,k_i) = \frac{b_i\cdot(b_i + 1)}{2} + k_i\cdot b_i$, 那么$\sum_{i = 1}^mf(b_i,k_i) = 0$可以转化为$\sum_{i = 0}^m{\frac{b_i\cdot (b_i + 1)}{2}} = \sum_{i = 0}^m k_i \cdot b_i$. 这是一个一次不定方程. 设$g = gcd(b_1,b_2,...,b_m)$, 当且仅当$g | \sum_{i = 0}^m \frac{b_i \cdot (b_i+1)}{2}$时有整数解.

手画一下几组样例可以看出, 当序列中存在奇数时必定合法. 当序列中存在至少一个奇数时, 显然$g$也是奇数. 那么对于每个$\frac{b_i \cdot (b_i + 1)}{2}$, 如果$b_i$是偶数有$g | \frac{b_i}{2}$, 否则有$g|(b_i + 1)$ 且$\frac{b_i}{2}$为整数. 因而序列中存在奇数时必有$g | \sum_{i = 0}^m \frac{b_i \cdot (b_i+1)}{2}$.

设$A$中奇数个数为$odd$, 这一部分对答案的贡献为$(2^{odd} - 1) \cdot 2^{n-odd}$

现在考虑序列中全为偶数的情况. 设$h$是最大的使得$2^h | g$的非负整数. 我们可以把$g$拆分为$g' \cdot 2^h$. 显然$g'$为奇数. 根据上一段的讨论, 有$g' | \sum_{i = 0}^m \frac{b_i \cdot (b_i + 1)}{2}$. 因而全为偶数的序列是否合法取决于$2^h$能否整除$\sum_{i = 0}^m \frac{b_i \cdot (b_i + 1)}{2}$. 对于每个$b_i$, 设$l$是最大的使得$2^l|b_i$的非负整数. $l$可以分三种情况来讨论:

$l > h$.

此时有$2^h | \frac{b_i}{2}$, 所以$2^h | \frac{b_i \cdot(b_i + 1)}{2}$.
$l = h$.

此时$2^h \nmid \frac{b_i}{2}$且$b_i + 1$是奇数, 因此$2^h \nmid \frac{b_i \cdot(b_i + 1)}{2}$. 由于$\frac{b_i \cdot(b_i + 1)}{2 \cdot g'} \mod 2^h = 2^{h-1}$, 如果这种情况下的$b_i$个数有偶数个, 那它们的和仍然可以被$2^h$整除.
$l < h$. 由于$g \geq 2^h$, 不可能存在这种情况.

综上可知, 当序列中全为偶数时, 按能被$2^l$整除的最大$l$来对$b_i$分类, 当且仅当$l$最小那组的个数为偶数时序列合法. 在统计时可以枚举最小的$l$.

总时间复杂度$n \log (\max\{a_i\})$

代码

//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define int long long
#define endl '\n'
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
const int maxn = 2e5+10;//CHANGE IT!
const int maxm = maxn*2;
const int p = 1e9+7;//CHANGE IT!
ll fac[maxn],cnt[maxn];
ll Pow(ll x,ll d,ll p = p){ll ta=1;if(d==0)return 1%p;x %= p;ll a=Pow(x,d/2,p);ta=a*a%p;if(d%2)ta=ta*x%p;return ta%p;}
ll C(ll n,ll m){
	if(m > n)return 0;
	return (fac[n]*Pow(fac[m],p-2)%p*Pow(fac[n-m],p-2)%p)%p;
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	fac[0] = 1;
	for(int i = 1;i < maxn;i++)fac[i] = fac[i-1]*i%p;
	int n;
	cin >> n;
	int odd = 0;
	for(int i = 1;i <= n;i++){
		int x;
		cin >> x;
		if(x&1)odd++;
		else{
			int qwq = 0;
			while(x%2 == 0)qwq++,x /= 2;
			cnt[qwq+1]++;
		}
	}
	ll ans = (Pow(2,odd) - 1 + p)%p*(Pow(2,n - odd)) % p;
	ll res = n - odd;
	for(int k = 1;k <= 32;k++)if(cnt[k]){
		ll tmp = 0;
		res -= cnt[k];
		ans = (ans + (Pow(2,cnt[k] - 1) - 1 + p)%p * Pow(2,res) % p) % p;
	}
	cout << ans;
	return 0;
}

solution-ABC213E

Posted on 2021-10-17

[题解]HDU7059-Counting Stars

Posted on 2021-09-16 In 题解

题意

本题是2021东北四省赛D-Lowbit 的强化版.

给出一个长度为$n$的数列$raw[]$, 维护下列三种操作:

1 l r查询$raw[l] + raw[l+1] + ... + raw[r]$对$998244353$取模的值

2 l r对区间$[l,r]$里的每个值分别减去其lowbit

3 l r对区间$[l,r]$里的每个值分别加上其最高位的值

$1 \leq n,q \leq 10^5$

题解

容易发现对于操作2, 每个数$a_i$在减去最多$\log a_i$次后就会变成0. 因而我们暴力模拟这个过程, 同时用$full[x] = 1$表示该节点下的所有值全部为0, 在之后的操作2中将这种节点忽略即可.

对于操作3, 只需把最高位拆出来单独维护即可. 注意操作2会影响到拆出来的最高位.

代码

//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
#define int long long
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
const int maxn = 2e5+10;//CHANGE IT!
const int maxm = maxn*2;
const int p = 998244353;//CHANGE IT!
int p2[maxn];
ll Pow(ll x,ll d){ll ta=1;if(d==0)return 1%p;x %= p;ll a=Pow(x,d/2);ta=a*a%p;if(d%2)ta=ta*x%p;return ta%p;}
struct Segtree{
	int siz = 1;
	vector<int>hbit,addv,subv,sumv;
	vector<bool>full;
	void init(int n){
		while(siz < n)siz *= 2;
		hbit.assign(siz*2,0);
		addv.assign(siz*2,0);
		sumv.assign(siz*2,0);
		full.assign(siz*2,0);
	}
	inline int ls(int x){return (x<<1)|1;}
	inline int rs(int x){return (x<<1)+2;}
	int lowbit(int x){return x & -x;}
	int highbit(int x){
		for(int i = 30;i >= 0;i--)if((1<<i)&x)return (1<<i);
		return 0;
	}
	void pushup(int x){
		if(full[ls(x)] and full[rs(x)])full[x] = 1;
		sumv[x] = (sumv[ls(x)] + sumv[rs(x)])%p;
		hbit[x] = (hbit[ls(x)] + hbit[rs(x)])%p;
	}
	void pushdown(int x){
		if(addv[x]){
			int & v = addv[x];
			addv[ls(x)] += v,addv[rs(x)] += v;
			hbit[ls(x)] = (hbit[ls(x)] * p2[v])%p;
			hbit[rs(x)] = (hbit[rs(x)] * p2[v])%p;
			v = 0;
		}
	}
	void build(vector<int> & in,int x,int lx,int rx){
		if(rx - lx == 1){
			if(lx < (int)in.size()){
				int h = highbit(in[lx]);
				hbit[x] = h;
				in[lx] -= h;
				sumv[x] = in[lx];
				if(!sumv[x] and !hbit[x])full[x] = 1;
			}
			return;
		}
		int m = (lx+rx)/2;
		build(in,ls(x),lx,m);
		build(in,rs(x),m,rx);
		pushup(x);
	}
	void build(vector<int> & in){
		build(in,0,0,siz);
	}
	void add(int l,int r,int x,int lx,int rx){
		if(l >= rx or lx >= r)return;
		if(rx - lx > 1)pushdown(x);
		if(l <= lx and rx <= r){
			hbit[x] = (hbit[x]*2)%p;
			addv[x]++;
			return;
		}
		int m = (lx+rx)/2;
		add(l,r,ls(x),lx,m);
		add(l,r,rs(x),m,rx);
		pushup(x);
	}
	void add(int l,int r){
		r++;
		add(l,r,0,0,siz);
	}
	void sub(int l,int r,int x,int lx,int rx){
		if(l >= rx or lx >= r)return;
		if(rx - lx > 1)pushdown(x);
		if(l <= lx and rx <= r){
			if(full[x]){
				return;
			}
			if(rx - lx == 1){
				if(!sumv[x]){
					full[x] = 1;
					hbit[x] = 0;
				}
				sumv[x] -= lowbit(sumv[x]);
				return;
			}
		}
		int m = (lx+rx)/2;
		sub(l,r,ls(x),lx,m);
		sub(l,r,rs(x),m,rx);
		pushup(x);
	}
	void sub(int l,int r){
		r++;
		sub(l,r,0,0,siz);
	}
	ll sum(int l,int r,int x,int lx,int rx){
		if(l >= rx or lx >= r)return 0;
		if(rx - lx > 1)pushdown(x);
		if(l <= lx and rx <= r){
			return (sumv[x] + hbit[x])%p;
		}
		int m = (lx+rx)/2;
		ll s1 = sum(l,r,ls(x),lx,m);
		ll s2 = sum(l,r,rs(x),m,rx);
		return (s1+s2)%p;
	}
	ll sum(int l,int r){
		r++;
		return sum(l,r,0,0,siz);
	}
};
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);
	p2[0] = 1;
	for(int i = 1;i < maxn;i++)p2[i] = p2[i-1]*2%p;	 
	int __;
	cin >> __;
	while(__--){
		int n;
		cin >> n;
		vector<int>in(n+1);
		for(int i = 1;i <= n;i++)cin >> in[i];
		Segtree seg;
		seg.init(n+10);
		seg.build(in);
		int q;
		cin >> q;
		while(q--){
			int opt,l,r;
			cin >> opt >> l >> r;
			if(opt == 1)cout << seg.sum(l,r) << endl;
			else if(opt == 2)seg.sub(l,r);
			else seg.add(l,r);
		}
	}
	return 0;
}

[题解]-CF1567E Non-Decreasing Dilemma

Posted on 2021-09-14 Edited on 2021-09-16 In 线段树

题意

给出$raw[]$和两种操作:

$x\space v$, 将$raw[x]$修改为$v$
$l \space r$, 询问$raw[l],raw[l+1],...,raw[r]$中有多少合法的数对$(i,j)$满足$l \leq i \leq j \leq r$且$raw[i]...raw[j]$单调不降

$1 \leq n,q \leq 2\cdot10^5$

题解

线段树维护区间合并. 以下$x$代表线段树节点标号, $ls(x)$和$rs(x)$分别代表左右儿子节点标号

容易发现对于每段长度为$m$的单调不降子序列, 其对答案的贡献是$\frac{m\cdot (m+1)}2$

建一颗线段树, 维护如下信息:

$prev[x]$: 自节点$x$的最左端开始的最长单调不降子序列的长度

$sufv[x]$: 自节点$x$的最右端结束的最长单调不降子序列的长度

$totv[x]$: 节点$x$有多少合法数对$(i,j)$满足题意

有几点细节需要注意:

对于$prev[x]$, 只有其左儿子的$prev[ls(x)]$等于区间长度, 并且左儿子和右儿子的连接处不降(即$raw[m-1] <= raw[m])$时才有$prev[x] = prev[ls(x)] + prev[rs(x)]$. $sufv[]$的维护同理
对于$totv[x]$, 显然有$totv[x] += totv[ls(x)] + totv[rs(x)]$. 而当其左儿子和右儿子的连接处不降时, 还需加上$sufv[ls(x)] \cdot prev[rs(x)]$, 即新增合法区间的左端点全部从左儿子中选, 右端点全部从右儿子中选.
在统计答案时, 需要处理当前节点区间比询问区间大的情况

代码

//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
#define int long long
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
const int maxn = 8e5+10;//CHANGE IT!
const int maxm = maxn*2;
const int p = 1e9+7;//CHANGE IT!
int raw[maxn];
ll Pow(ll x,ll d,ll p){ll ta=1;if(d==0)return 1%p;x %= p;ll a=Pow(x,d/2,p);ta=a*a%p;if(d%2)ta=ta*x%p;return ta%p;}
struct Segtree{
	int siz = 1;
	vector<int>prev,sufv,totv;
	void init(int n){
		while(siz < n)siz *= 2;
		prev.assign(siz*2,0);
		sufv.assign(siz*2,0);
		totv.assign(siz*2,0);
	}
	inline int ls(int x){return (x<<1)|1;}
	inline int rs(int x){return (x<<1)+2;}
	inline int f(int x){return x*(x+1)/2;}
	void pushup(int x,int lx,int rx){
		int m = (lx+rx)/2;
		int len = (rx - lx)/2;
		totv[x] = totv[ls(x)] + totv[rs(x)];
		totv[x] += (raw[m-1] <= raw[m]?sufv[ls(x)]*prev[rs(x)]:0);
		prev[x] = max(prev[ls(x)],(prev[ls(x)] == len and raw[m-1] <= raw[m])?prev[ls(x)] + prev[rs(x)]:0);
		sufv[x] = max(sufv[rs(x)],(sufv[rs(x)] == len and raw[m-1] <= raw[m])?sufv[rs(x)] + sufv[ls(x)]:0);
	}
	void Set(int idx,int v,int x,int lx,int rx){
		if(rx - lx == 1){
			raw[idx] = v;
			prev[x] = sufv[x] = totv[x] = 1;
			return;
		}
		int m = (lx+rx)/2;
		if(lx <= idx and idx < m)Set(idx,v,ls(x),lx,m);
		else Set(idx,v,rs(x),m,rx);
		pushup(x,lx,rx);
	}
	void Set(int idx,int v){
		Set(idx,v,0,0,siz);
	}
	ll cal(int l,int r,int x,int lx,int rx){
		if(l >= rx or lx >= r)return 0;
		if(l <= lx and rx <= r)return totv[x];
		int m = (lx+rx)/2;
		ll s1 = cal(l,r,ls(x),lx,m);
		ll s2 = cal(l,r,rs(x),m,rx);
		ll ans = s1+s2;
		if(raw[m-1] <= raw[m]){
			ans += max(min(m - l,sufv[ls(x)]),0LL)*max(min(r - m,prev[rs(x)]),0LL);
		}
		return ans;
	}
	ll cal(int l,int r){
		r++;
		return cal(l,r,0,0,siz);
	}
};
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	int n,q;
	Segtree seg;
	cin >> n >> q;
	for(int i = 1;i <= n;i++)raw[i] = 0;
	seg.init(n+1);
	for(int i = 1;i <= n;i++){
		int v;
		cin >> v;
		seg.Set(i,v);
	}
	while(q--){
		int opt,l,r;
		cin >> opt >> l >> r;
		if(opt == 1)seg.Set(l,r);
		else cout << seg.cal(l,r) << endl;
	}
	return 0;
}

[题解]CF1561D/CF1560B-Up the Strip

Posted on 2021-08-25 In 题解

题解

设$f[i]$为从$n$走到$i$的方案数. 分减法和整除两种情况分别讨论: \[ \begin{align} 减法:&\\ f[i] &= \sum_{j = i}^nf[j]\\ 整除:&\\ f[i] &= \sum_{j = i}^{i*j<= n}\sum_{\lfloor \frac{k}{j} \rfloor = i}f[k]\\ &= \sum_{j = i}^{i*j<= n}\sum_{k = i*j}^{(i+1)*j - 1}f[k]\\ &= \sum_{j = i}^{i*j<= n}suf[j] - suf[(i+1)*j]\\ 其中&suf[i] = \sum_{j = i}^n f[i]\\ \end{align} \]

综上, 递推式为 \[ \begin{align} f[n] &= suf[n] = 1\\ suf[i] &=suf[i+1] + f[i]\\ f[i] &= suf[i+1]+\sum_{j = 1}^{i*j <= n}suf[j] - suf[(i+1)*j] \end{align} \] 注意处理$(i+1)*j$越界的问题. 复杂度$O(n\ln n)$

代码

//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define int long long
#define endl '\n'
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
const int maxn = 4e6+10;//CHANGE IT!
const int maxm = maxn*2;
int p;
ll Pow(ll x,ll d){ll ta=1;if(d==0)return 1%p;ll a=Pow(x,d/2);ta=a*a%p;if(d%2)ta=ta*x%p;return ta%p;}
ll f[maxn],pre[maxn];
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	int n;
	cin >> n >> p;
	f[n] = pre[n] = 1;
	for(int i = n-1;i >= 1;i--){
		f[i] = pre[i+1];
		for(int j = 2;i*j <= n;j++){
			f[i] = (f[i] + (pre[i*j] - ((i+1)*(j) > maxn?0:pre[(i+1)*(j)])))%p;
		}
		pre[i] = (f[i] + pre[i+1])%p;
	}
	cout << f[1];
	return 0;
}

[题解]HDU7047-Link with Balls

Posted on 2021-08-11 Edited on 2021-08-16 In 题解

题意

给你两类盒子, 每类$n$个, 编号分别从$1$到$n$.

第一类: 可以从第$i$个盒子里拿出$i$的任意非负整数倍数个小球

第二类: 可以从第$i$个盒子里拿出最少$0$个最多$i$个小球.

多组数据$(1 \leq t \leq 10^5)$, 每组数据给出$1 \leq n,m \leq 10^6$, 求从这$2n$个盒子里拿出$m$个小球的方案数模$1e9+7$.

题解

懒得从组合意义的角度考虑, 可以用生成函数去做.

第一类第$i$个盒子对应的生成函数是 \[ \sum_{j = 0}^{∞}x^{i\cdot j} \]

第二类对应的是 \[ \sum_{j = 0}^ix^i \]

那么答案即为: \[ [m](\prod_{i = 1}^{n}\sum_{j = 0}^{∞}x^{i\cdot j})\cdot(\prod_{i = 1}^{n}\sum_{j = 0}^ix^i) \\ =[m](\prod_{i = 1}^{n}\frac{1}{1 - x^i})\cdot(\prod_{i = 1}^{n}\frac{1-x^{i+1}}{1-x})\\ =[m]\frac{1-x^{n+1}}{(1-x)^{n+1}}\\ =[m](1-x^{n+1})\cdot \sum_{i = 0}^{∞}C_{n+i}^ix^i\\ = C_{n+m}^m - C_{m-1}^{n} \]

[题解]CF958F3 Lightsabers (hard)

Posted on 2021-04-04 Edited on 2021-04-05 In 题解

题意

给出$m$个集合, 第$i$个集合有$c_i$种颜色为$i$的小球, 保证$\sum_{i = 1}^m c_i = n$. 现从这些集合中取出$k$个小球, 问这些小球构成的集合有多少种. 答案模$1009$.

$1 \leq m,k \leq n \leq 2\cdot 10^5$

题解

普通型生成函数裸题.

设$G(x) = \prod_{i = 1}^m(\sum_{j = 0}^{c_i}x^j)$

那么$G(x)[k]$即为答案. 问题转换为$m$个多项式相乘. 使用分治法可在$O(n\log^2n)$时间内用FFT求解.

代码

//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
const int maxn = 4e6+10;
const int p = 1009;
int raw[maxn],cnt[maxn];
typedef double ld;
typedef vector<int> Poly;
const ld PI = acos(-1.0);
int rev[maxn];
struct Comp{
    double x, y;
    Comp (double x = 0, double y = 0):x(x),y(y){}
    Comp operator * (Comp B){
   		return Comp(x*B.x - y*B.y,x*B.y + y*B.x);
	}
	Comp operator + (Comp B){
	    return Comp(x + B.x,y + B.y);
	}
	Comp operator - (Comp B){
	    return Comp(x - B.x,y - B.y);
	}
};
Comp F1[maxn],F2[maxn];
void FFT(Comp * A,int siz,int type){
	int n = siz;
	int S = log2(n);
	for(int i = 0;i < n;i++)rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (S - 1));
	for(int i = 0;i < n;i++)if(i < rev[i])swap(A[i],A[rev[i]]);
	for(int i = 1;i < n;i *= 2){
		Comp Wn(cos(PI/i),type*sin(PI/i));
		for(int j = 0;j < n;j += i*2){
			Comp W(1,0);
			for(int k = 0;k < i;k++){
				Comp facx = A[j+k],facy = W*A[j+k+i];
				A[j+k] = facx + facy;
				A[j+k+i] = facx - facy;
				W = W*Wn;
			}
		}
	}
    if(type == -1)for(int i = 0;i < n;i++)A[i].x = ((A[i].x/n + 0.5));
}
Poly mul(Poly A,Poly B){
	int n = A.size(),m = B.size();
	int siz = n + m - 1;
	Poly C(siz);
	if(siz < 64){
		for(int i = 0;i < n;i++){
			for(int j = 0;j < m;j++)C[i+j] = (C[i+j] + 1LL*A[i]*B[j])%p;
		}
		return C;
	}
	int fsiz = 1;
	while(fsiz <= siz)fsiz *= 2;
	for(int i = 0;i < fsiz;i++)F1[i] = F2[i] = 0;
	for(int i = 0;i < n;i++)F1[i] = A[i];
	for(int i = 0;i < m;i++)F2[i] = B[i];
	FFT(F1,fsiz,1);
	FFT(F2,fsiz,1);
	for(int i = 0;i < fsiz;i++)F1[i] = F1[i]*F2[i];
	FFT(F1,fsiz,-1);
	for(int i = 0;i < siz;i++){
		C[i] = ((ll)F1[i].x)%p;
	}
	return C;
}
Poly solve(int l,int r){
	if(l == r){
		return Poly(cnt[l]+1,1);
	}
	int m = (l+r)/2;
	return mul(solve(l,m),solve(m+1,r));
}
Poly ans;
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	int n,m,k;
	cin >> n >> m >> k;
	for(int i = 1;i <= n;i++){
		int x;
		cin >> x;
		cnt[x]++;
	}
	Poly ans = solve(1,m);
	cout << ans[k] << endl;
	return 0;
}

GAMES101笔记

Posted on 2021-04-04 Edited on 2023-02-26 In 笔记

部分图片来自GAMES101

向量

基础

可以直接跳过这部分.

定义

符号表示: $\vec{a}$ 或$\boldsymbol{a}$, 或者写成$\vec{AB}$, 表示从$A$指向$B$的向量($\vec{AB} = B - A$).

本文中向量均表示列向量, 即$\vec{a} = \begin{pmatrix}x\\y\end{pmatrix}$. 行向量则$\vec{a}^T = (x,y)$

向量有方向和长度, 一般不指定起点.

长度

$|\vec{a}|$ 或$\Vert\vec{a}\Vert$(二范数)表示向量的长度. 对于二维向量$\vec{a} = \begin{pmatrix}x\\y\end{pmatrix}$, $\vert a\vert = \sqrt{x^2 + y^2}$.

单位向量

$\hat{a}$表示$\vec{a}$的单位向量, 读作a hat. \[ \hat{a} = \frac{\vec{a}}{\vert a \vert} \]

一般用单位向量来表示方向

基本运算

加法

\[ \vec{a} = \begin{pmatrix}x_1\\y_1\end{pmatrix}, \vec{b} = \begin{pmatrix}x_2\\y_2\end{pmatrix}. \space \vec{a} + \vec{b} = \begin{pmatrix}x_1+x_2\\y_1+y_2\end{pmatrix} \]

向量的加法遵循平行四边形法则或三角形法则(如下图)

对于多个向量相加, 按三角形法则将它们依序首尾拼接即可.

点积

\[ \vec{a} \cdot \vec{b} = \Vert\vec{a}\Vert\Vert\vec{b}\Vert\cos\theta \]

\[ \vec{a} \cdot \vec{b} = \begin{pmatrix}x_1\\y_1\end{pmatrix} \cdot \begin{pmatrix}x_2\\y_2\end{pmatrix} = x_1\cdot x_2 + y_1 \cdot y_2 \]

可以写成矩阵形式: \[ \vec a \cdot \vec b = \vec{a}^T\vec{b} \\ = \left( \begin{array}{cccc} x_1 & y_1 & z_1\\ \end{array} \right) \left( \begin{array}{cccc} x_2\\ y_2\\ z_2 \end{array} \right) \] 点积运算符合交换、结合、分配率: \[ \begin{align} \vec{a} \cdot \vec{b} &= \vec{b}\cdot\vec{a}\\ \vec{a} \cdot(\vec{b} + \vec{c}) &= \vec{a} \cdot \vec{b} + \vec{a} \cdot \vec{c}\\ (k \vec{a}) \cdot \vec{b} &= \vec{a} \cdot(k\vec{b}) = k(\vec{a}\cdot\vec{b}) \end{align} \]

两向量夹角

\[ \cos\theta = \frac{\vec{a}\cdot\vec{b}}{\Vert\vec{a}\Vert\Vert\vec{b}\Vert} = \hat{a}\cdot\hat{b} \]

投影

$\vec{b}$ 在$\vec a$上的投影记作$\vec{b}_{\perp}$, 读作b perp:

(垂直: perpendicular)

\[ \vec{b}_{\perp} = k\hat{a} \] 其中$k = \Vert \vec{b}_{\perp}\Vert = \Vert \vec{b}\Vert \cos \theta$

投影可以将向量$\vec b$分解成平行于$\vec a$($\vec{b}_{\perp}$)和垂直于$\vec a$($\vec{b} - \vec{b}_{\perp}$)两部分, 如下图:

判断两向量方向

如图, 两向量的点积大于$0$则相向, 小于$0$则反向, 等于$0$则互相垂直.

叉积

\[ \vec{a} = \begin{pmatrix}x_1\\y_1\\z_1\end{pmatrix},\vec{b} = \begin{pmatrix}x_2\\y_2\\z_2\end{pmatrix} \]

$\vec{a} \times \vec{b}$可用如下矩阵的行列式来计算: \[ \vec{a} \times \vec{b} = \left| \begin{array}{cccc} \vec{i} & \vec{j} & \vec{k} \\ x_1 & y_1 & z_1\\ x_2 & y_2 & z_2 \end{array} \right| \\= (y_1z_2 - y_2z_1) \vec{i} + (x_2z_1 - x_1z_2)\vec{j} + (x_1y_2 - x_2y_1)\vec{k} \]

写成矩阵形式: \[ \vec a \times \vec b = A^* \vec b = \left( \begin{array}{cccc} 0 & -z_1 & y_1\\ z_1 & 0 & -x_1\\ -y_1 & x_1 & 0 \end{array} \right) \left( \begin{array}{cccc} x_2\\ y_2\\ z_2 \end{array} \right) \] $A^*$称为dual martix.

两向量$\vec a$ $\vec b$ 叉积得到一个垂直于二者的向量.

该向量的方向可用右手螺旋定则判断: 伸出右手, 将四指按从$\vec a$扫向$\vec b$的方向弯曲, 拇指的指向就是该向量的方向. 遵从该定则的坐标系称为右手系.

叉积的运算不满足交换律, 基本关系如下:

\[ \vec a \times \vec b = - \vec b \times \vec a\\ \vec a \times \vec a = \vec 0\\ \vec a \times (\vec b + \vec c) = \vec a \times \vec b + \vec a \times \vec c\\ \vec a \times (k \vec b) = k (\vec a \times \vec b) \]

判断两向量左右关系

如图, 计算$\vec{c} = \vec{a} \times \vec{b}$, 若$\vec{c}$与$z$轴同向($z$的系数为正)则$\vec{b}$在$\vec{a}$的左边, 否则$\vec{b}$在$\vec{a}$的右边.

判断点是否在凸多边形内部

将凸多边形的点逆时针/顺时针排列, 将待判断的点与每条边进行叉乘, 若点均在边的左侧/右侧则该点在凸多边形内部. (注意这是$O(n)$的,$O(\log n)$的做法右转计算几何板子).

正交坐标系

假设向量$\vec{x},\vec{y},\vec{z}$满足 \[ \Vert \vec{x} \Vert = \Vert \vec{y} \Vert = \Vert \vec{z} \Vert = 1\\ \vec{x} \cdot \vec{y} = \vec{x} \cdot \vec{z} = \vec{y} \cdot \vec{z} = 0\\ \vec{x} \times \vec{y} = \vec{z} \] 那么$\vec{x},\vec{y},\vec{z}$组成一个正交坐标系, 且是右手系.

对于任意向量$\vec p$, 可通过投影将其分解为三个基底的和: \[ \vec p = (\vec p \cdot \vec x)\vec x + (\vec p \cdot \vec y )\vec y + (\vec p \cdot \vec z)\vec z \]

矩阵

矩阵乘法

设$A$是一个$n \times m$的矩阵, $B$是一个$m \times k$的矩阵, 那么$A \times B$得到一个$n \times k$的矩阵$C$, 其中 \[ C_{i,j} = \sum_{k = 1}^m A_{i,k} \cdot B_{k,j} \] 矩阵乘法不存在交换律 \[ AB(C) = A(BC)\\ A(B+C) = AB + AC\\ (A+B)C = AC + BC \]

转置

将行和列互换. 例如 \[ \left( \begin{array}{cccc} 1 & 2 & 3\\ 4 & 5 & 6 \end{array} \right)^T = \left( \begin{array}{cccc} 1 & 4\\ 2 & 5\\ 3 & 6 \end{array} \right) \]

\[ (AB)^T = B^TA^T \]

单位矩阵

主对角线为$1$, 其余元素为$0$的方阵, 记为$I_n$ 例如 \[ I_3 = \left( \begin{array}{cccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{array} \right) \]

\[ AI_m = I_nA = A \]

可以用来定义矩阵的逆: \[ AA^{-1} = I \]

\[ (AB)^{-1} = B^{-1}A^{-1} \]

平面

三维平面的一般式方程: \[ Ax + By + Cz + D = 0 \] 其中平面法线为$\vec n = (A,B,C)^T$

对于任意一点$P(x_0,y_0,z_0)$, 其到平面的距离为$d = Ax_0 + By_0 + Cz_0 + D$. 可见$D$为原点到平面的距离. 当$d > 0$时说明点$P$在平面法线所指一侧.

若两点$P_1,P_2$分列于平面两侧, 且它们到平面的距离为$d_1,d_2$, 则两点连线与平面的交点可以插值得到: \[ \lambda = \frac{|d_1|}{|d_1 + d_2|} = \frac{d_1}{d_1 - d_2} \] \[ P_{inter} = \lambda P_2 + (1 - \lambda)P_1 \]

变换基础知识

线性变换

可以表示为$\vec{x'} = A\vec{x}$的变换均为线性变换.

缩放(Scale)

将横纵坐标按比例缩放: \[ x' = sx\\ y' = sy \] 写为矩阵形式: \[ \left[ \begin{array}{cccc} x'\\ y' \end{array} \right] = \left[ \begin{array}{cccc} s & 0\\ 0 & s \end{array} \right]\left[ \begin{array}{cccc} x\\ y \end{array} \right] \]

横纵坐标按不同比例缩放: \[ \left[ \begin{array}{cccc} x'\\ y' \end{array} \right] = \left[ \begin{array}{cccc} s_x & 0\\ 0 & s_y \end{array} \right]\left[ \begin{array}{cccc} x\\ y \end{array} \right] \]

反射(Reflection)

\[ \left[ \begin{array}{cccc} x'\\ y' \end{array} \right] = \left[ \begin{array}{cccc} -1 & 0\\ 0 & 1 \end{array} \right]\left[ \begin{array}{cccc} x\\ y \end{array} \right] \]

切变(Shear)

\[ x' = x + ay\\ y' = y \]

\[ \left[ \begin{array}{cccc} x'\\ y' \end{array} \right] = \left[ \begin{array}{cccc} 1 & a\\ 0 & 1 \end{array} \right]\left[ \begin{array}{cccc} x\\ y \end{array} \right] \]

旋转(Rotate)

默认为以$(0,0)$为中心逆时针旋转.

\[ R_\theta = \left[ \begin{array}{cccc} \cos\theta & -\sin\theta\\ \sin\theta & \cos\theta \end{array} \right] \]

\[ R_{-\theta} = \left[ \begin{array}{cccc} \cos\theta & \sin\theta\\ -\sin\theta & \cos\theta \end{array} \right] = R_\theta^T = R_\theta^{-1} \]

说明旋转矩阵是一个正交矩阵.

如果旋转的中心不在原点$(0,0)$, 需要先将中心点平移至原点后旋转, 旋转后再将中心点平移至原位置.

平移(Translation)

\[ x' = x + t_x\\ y' = y + t_y \] 我们无法把它写成矩阵相乘的形式, 因此平移不是线性变换.

使用齐次坐标表示后(见下文), 平移也可表示为矩阵形式: \[ \left( \begin{array}{cccc} x'\\ y'\\ w' \end{array} \right) = \left( \begin{array}{cccc} 1 & 0 & t_x\\ 0 & 1 & t_y\\ 0 & 0 & 1 \end{array} \right)\cdot \left( \begin{array}{cccc} x\\ y\\ 1 \end{array} \right) = \left( \begin{array}{cccc} x + t_x\\ y + t_y\\ 1 \end{array} \right) \]

齐次坐标(Homogeneous Coordinates)

我们把二维的点和向量添加一维, 变成齐次坐标$(x,y,w)^T$. \[ \tt{point} = \left( \begin{array}{cccc} x\\ y\\ 1 \end{array} \right)\tt{,vector = }\left( \begin{array}{cccc} x\\ y\\ 0 \end{array} \right) \] 对于$w \neq 0$, \[ \left( \begin{array}{cccc} x\\ y\\ w \end{array} \right) 表示点\left( \begin{array}{cccc} \frac{x}{w} \\ \frac{y}{w} \\ 1 \end{array} \right) \] 根据$w$的值可以轻易地判断两个齐次坐标相加后的结果:

点 + 向量 = 点
向量 + 向量 = 向量
点 - 点 = 向量
点 + 点在齐次坐标下表示两点的中点

仿射变换(Affine Transformation)

对一个向量空间进行一次线性变换后再进行平移, 称为仿射变换. \[ \left( \begin{array}{cccc} x'\\ y'\\ \end{array} \right) = \left( \begin{array}{cccc} a & b\\ c & d\\ \end{array} \right) \cdot \left( \begin{array}{cccc} x\\ y\\ \end{array} \right) + \left( \begin{array}{cccc} t_x\\ t_y\\ \end{array} \right) \] 写成齐次坐标形式: \[ \left( \begin{array}{cccc} x'\\ y'\\ 1 \end{array} \right) = \left( \begin{array}{cccc} a & b & t_x\\ c & d & t_y\\ 0 & 0 & 1 \end{array} \right) \cdot \left( \begin{array}{cccc} x\\ y\\ 1 \end{array} \right) \]

其中, 左上的2x2矩阵为原线性变换矩阵, $t_x$, $t_y$表示平移. 三维与其类似. 注意: 这个矩阵表示先进行线性变换, 再进行平移.

逆变换(Inverse Transform)

求逆即可

变换的组合

假设对$\vec x$依次应用变换$M_1,M_2,...,M_n$, 那么有 \[ \vec{x'} = M_n \cdot M_{n-1} \cdots M_2 \cdot M_1 \cdot \vec{x} \]

变换的分解

例:按某点旋转

假设我们要按点$C$旋转$\alpha$度, 可以先将点$C$平移回原点$T(-c)$, 进行旋转$R(\alpha)$, 再将点$C$移动回原位. \[ \vec{x'} = T(c) \cdot R(\alpha) \cdot T(-c) \cdot \vec{x} \]

三维变换

此处仅列出和二维有较大不同的部分.

旋转

按坐标轴旋转

按$x$轴旋转$\alpha$度 \[ R_x(\alpha) = \left( \begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & \cos \alpha & -\sin\alpha & 0\\ 0 & \sin\alpha & \cos\alpha & 0\\ 0 & 0 & 0 & 1 \end{array} \right) \]
按$y$轴(注意和其它两个的区别, $ x z = -y$) \[ R_y(\alpha) = \left( \begin{array}{cccc} \cos\alpha & 0 & \sin\alpha & 0\\ 0 & 1 & 0 & 1\\ -\sin\alpha & 0 & \cos\alpha & 0\\ 0 & 0 & 0 & 1 \end{array} \right) \]
按$z$轴 \[ R_z(\alpha) = \left( \begin{array}{cccc} \cos\alpha & -\sin\alpha & 0 & 0 \\ \sin\alpha & \cos\alpha & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{array} \right) \]

一般性的旋转

我们用$R_{xyz}(\alpha, \beta, \gamma)$表示绕$x$轴旋转$\alpha$, 绕$y$轴旋转$\beta$, 绕$z$轴旋转$\gamma$. $\alpha\beta\gamma$也被叫做欧拉角. \[ R_{xyz}(\alpha, \beta, \gamma) = R_x(\alpha)R_y(\beta)R_z(\gamma) \]

Rodrigues旋转公式

$R(\vec n,\alpha)$表示绕轴$\vec n$(起点在原点)按右手方 m 向旋转角度为$\alpha$的变换. \[ R(\vec n,\alpha) = \cos(\alpha)\vec I + (1 - \cos(\alpha))\vec n \vec n^T + \sin(\alpha)\left( \begin{array}{cccc} 0 & -\vec{n}_{z} & \vec{n}_y\\ \vec{n}_z & 0 & -\vec{n}_x\\ -\vec{n}_y & \vec{n}_x & 0\\ \end{array} \right) \]

欧拉角

(待补充)四元数(Quaternion)

四元数是一种有三个虚部的复数: \[ q = a + b i + c j + d k \] 其中有 \[ i^2 = j^2 = k^2 = ijk = -1 \] 四元数的加法只需把对应系数相加. 乘法每个虚部遵从如下乘法表:

和虚数的共轭类似, 有四元数$q$的共轭$q^*$: \[ q^* = a - bi - cj - dk \] 注意: $(pq)^* = q^*p^*$

一个四元数的绝对值$|q|$定义为: \[ |q| = \sqrt{q\cdot q^*} = \sqrt{a^2 + b^2 + c^2 + d^2} \] 其中绝对值为1的四元数称为单位四元数.

四元数的逆(乘逆)$q^{-1}$: \[ q^{-1} = \frac{q^*}{|q|^2} \]

实部为0($a = 0$)的四元数称为纯四元数. 一个三维空间的坐标$(x,y,z)$可以用纯四元数$xi + yj +zk$来表示.

渲染过程中的各种变换

推荐阅读:http://www.codinglabs.net/article_world_view_projection_matrix.aspx

假设我们要拍一张照片，需要以下三个步骤:

安排好要拍的人或物. 对应模型变换(model transformation), 即将物体坐标转换为世界坐标.
寻找一个放置相机的位置和角度. 对应视图变换(view transformation).
拍照. 对应投影变换(projection transformation).

模型变换(Model Transformation)

这一步是将模型坐标转换为世界坐标的过程, 转换矩阵记为$M_{model}$

视图变换(View Transformation)

要进行视图变换, 我们首先需要定义相机的摆放:

位置(Position)$\vec e$
视角(Look-at direction/gaze direction)$\hat g$
向上方向(Up direction)$\hat t$

因为我们写的是一个渲染引擎而不是在现实世界，所以完全可以搞一个以相机为中心的世界. 约定相机永远位于原点，向上方向是$\text Y$轴, 指向$\text Z$轴的反方向. 再次提醒坐标系是右手系, 且$\text Y$是纵坐标轴:

那么对于一个给定的相机位置，需要进行一系列变换将其变成标准的相机位置. 我们把这个变换矩阵记作$M_{view}$, 具体需要进行如下变换:

将$\vec e$平移到原点
将$\vec g$旋转到$-\vec Z$
将$\hat t$旋转到$\vec Y$
将$\hat g \times \hat t$旋转到$\vec X$

根据前面的内容, 我们需要先平移($T_{view}$)再旋转$(R_{view})$, 即$M_{view} = R_{view}T_{view}$ \[ T_{view} = \left( \begin{array}{cccc} 1 & 0 & 0 & -x_e\\ 0 & 1 & 0 & -y_e\\ 0 & 0 & 1 & -z_e\\ 0 & 0 & 0 & 1 \end{array} \right) \] 同样根据前面的内容, 我们有$R_{-\theta} = R_\theta^{-1} = R_\theta^T$. 因此与其正向进行旋转, 不如先把$X\space Y\space-Z$轴旋转到$etg$, 再求它的逆矩阵. 因为我们要旋转的向量都很简单(比如$(1,0,0)^T$), 所以这个矩阵可以直接构造出来: \[ R_{view}^{-1} = \left( \begin{array}{cccc} x_{\hat g \times \hat t} & x_\hat t & x_{-\hat g} & 0\\ y_{\hat g \times \hat t} & y_\hat t & y_{-\hat g} & 0\\ z_{\hat g \times \hat t} & z_\hat t & z_{-\hat g} & 0\\ 0 & 0 & 0 & 1 \end{array} \right) \] 于是有 \[ R_{view} = (R_{view}^{-1})^T \left( \begin{array}{cccc} x_{\hat g \times \hat t} & y_{\hat g \times \hat t} & z_{\hat g \times \hat t} & 0\\ x_\hat t & y_\hat t & z_\hat t & 0\\ x_{-\hat g} & y_{-\hat g} & z_{-\hat g} & 0\\ 0 & 0 & 0 & 1 \end{array} \right) \]

在实现时, 我们通常不指定视角$\vec g$, 而是通过视角所指向的点(设为$C$)和相机位置($E$)来确定$\vec g$. 这时有: \[ \vec g = \vec{EC}\\ \] 有了$\vec g,\vec t$, 便可以计算出相机坐标系: \[ \vec w = -\frac{\vec g}{|\vec g|}\\ \vec u = \frac{\vec t \times \vec w}{|\vec t \times \vec w|}\\ \vec v = \vec w \times \vec u \] 这时$R_{view}$为 \[ R_{view} = (R_{view}^{-1})^T \left( \begin{array}{cccc} x_{\vec u} & y_{\vec u} & z_{\vec u} & 0\\ x_{\vec v} & y_{\vec v} & z_{\vec v} & 0\\ x_{\vec w} & y_{\vec w} & z_{\vec w} & 0\\ 0 & 0 & 0 & 1 \end{array} \right) \]

LookAt

投影变换(Projection Transformation)

投影变换将3D的图形转换到2D上，分为正交投影(Orthographic Projection)和透视投影(Perspective Projection)

正交投影(Orthographic Projection)

其实直接将所有的$z$值忽略掉就可以获得正交投影. 但一般采用如下方法:

我们先定义一个$[l,r] \times [b,t] \times [f,n]$的立方体(left,right,bottom,top,far,near). (注意, $z$值越小离相机越远, 因为相机的视线是-Z) 然后我们将这个立方体先平移再缩放, 得到中心位于原点, $[-1,1] \times[-1,1] \times [-1,1]$的标准立方体(canonical cube).

可以很简单地写出变换矩阵: \[ M_{ortho} = \left( \begin{array}{cccc} \frac{2}{r - l} & 0 & 0 & 0\\ 0 & \frac{2}{t - b} & 0 & 0\\ 0 & 0 & \frac{2}{n - f} & 0\\ 0 & 0 & 0 & 1 \end{array} \right)\left( \begin{array}{cccc} 1 & 0 & 0 & -\frac{l+r}{2} \\ 0 & 1 & 0 & -\frac{b+t}{2} \\ 0 & 0 & 1 & -\frac{f+n}{2} \\ 0 & 0 & 0 & 1 \end{array} \right) \]

相乘之后的结果:

透视投影(Perspective Projection)

进行透视投影时, 我们假设相机是一个点, 而所有光线都汇聚到这个点上(或者从这个点发出):

在三维空间上我们的视野便形成了一个四棱锥(因为你的电脑屏幕是长方形而不是圆形). 如果以待渲染的画面所在平面为顶将这个四棱锥的头截去, 便得到了一个截头四棱锥(frustum). 将这个四棱锥挤压成一个立方体, 再进行正交投影, 便能得到想要的画面:

考虑在挤压时Frustum中每个点的$x,y$会如何变化.

上图可以看作Frustum的一个侧视图. 一个点的$y$在挤压后变成了$y'$. 根据相似三角形很容易得出 \[ y' = \frac{n}{z}y \] 同样 \[ x' = \frac{n}{z}x \] 假设挤压过程对应的矩阵是$M_{persp\to ortho}$. 此时我们还不知道$z$坐标会如何变化(注意它不是不变).

注意：可以利用矩阵的最后一行对$x,y,z$进行除法运算.

举例:

根据 \[ M_{persp \to ortho} \cdot \left( \begin{array}{cccc} x \\ y \\ z \\ 1 \end{array} \right) = \left( \begin{array}{cccc} \frac{n}{z}x \\ \frac{n}{z}y \\ ? \\ 1 \end{array} \right) = \left( \begin{array}{cccc} nx \\ ny \\ ? \\ z \end{array} \right) \] 可以先填出矩阵的部分内容: \[ M_{persp \to ortho} = \left( \begin{array}{cccc} n & 0 & 0 & 0\\ 0 & n & 0 & 0\\ ? & ? & ? & ?\\ 0 & 0 & 1 & 0 \end{array} \right) \] 这个矩阵的第三行需要进行一点推导. 尽管我们不知道每个点的$z$坐标如何变换, 但我们知道, 在最前面平面($z = n$)上的点, 其$z$坐标永远等于$n$, 即: \[ \left( \begin{array}{cccc} n & 0 & 0 & 0\\ 0 & n & 0 & 0\\ ? & ? & ? & ?\\ 0 & 0 & 1 & 0 \end{array} \right) \cdot \left( \begin{array}{cccc} x \\ y \\ n \\ 1 \end{array} \right) = \left( \begin{array}{cccc} x \\ y \\ n \\ 1 \end{array} \right) = \left( \begin{array}{cccc} nx \\ ny \\ n^2 \\ n \end{array} \right) \] 可知$M[2][0] = M[2][1] = 0$. 设$M[2][2]$为$A$, $M[2][3]$为$B$, 则 \[ (0,0,A,B) \cdot \left( \begin{array}{cccc} x \\ y \\ n \\ 1 \end{array} \right) = n^2 \]

\[ An + B = n^2 \]

同理有 \[ Af + B = f^2 \] 解得 \[ A = n + f\\ B = -nf \] 于是我们求出了整个矩阵: \[ M_{persp \to ortho} = \left( \begin{array}{cccc} n & 0 & 0 & 0\\ 0 & n & 0 & 0\\ 0 & 0 & n+f & -nf\\ 0 & 0 & 1 & 0 \end{array} \right) \] 整个透视投影的过程便是先挤压再正交投影的过程, 即 \[ M_{persp} = M_{ortho}M_{persp\to ortho} \]

做一下矩阵乘法, 这个矩阵的结果是: \[ M_{persp } = \left( \begin{array}{cccc} \frac{2n}{r - l} & 0 & \frac{l + r}{l - r} & 0 \\ 0 & \frac{2n}{t - b} & \frac{b + t}{b - t} & 0 \\ 0 & 0 & \frac{f+n}{n - f} & -\frac{2 f n}{n - f} \\ 0 & 0 & 1 & 0 \end{array} \right) \] 可以看到, 这个矩阵的最后一行是$[0,0,1,0]$. 也就是说, 一个点在进行透视投影后, 它的深度信息$z$就是齐次坐标的$w$.

上述矩阵还可以进行优化. 一般情况下, 我们的视锥是对称的(即$l = -r, t = -b$). 这种情况下有: \[ M_{persp } = \left( \begin{array}{cccc} \frac{n}{r} & 0 & 0 & 0 \\ 0 & \frac{n}{t} & 0 & 0 \\ 0 & 0 & \frac{f+n}{n - f} & -\frac{2 f n}{n - f} \\ 0 & 0 & 1 & 0 \end{array} \right) \] 注意: 如果在OpenGL下使用上述矩阵, 需要进行转置.

构造投影矩阵

来自作业1. 参考资料: https://zhuanlan.zhihu.com/p/361156478

https://en.wikipedia.org/wiki/Field_of_view_in_video_games

有些时候, 我们没有$l,r,b,t,n,f$六个参数, 而是需要用以下四个参数，构造投影矩阵:

float eyefov: (垂直)可视角度

float aspect_ratio: 宽高比

float zNear: 近平面的$z$轴坐标

float zFar: 远平面的$z$轴坐标

如图, 约定$V$代表垂直视角, $H$代表水平视角.$w$和$h$代表显示器的宽度和高度. $l,r,b,t,n,f$见投影变换部分. 在此假设所有角已经转换成弧度制. 易得: \[ t = \tan(\frac{V}{2}) \cdot \text{zNear}\\ b = -t \] 因为$\text{aspect\_ratio} = \frac{w}{h} = \frac{r}{t}$, 所以有 \[ r = \text{aspect\_ratio} \cdot t\\ l = -r \] 这样便可构建出投影变换矩阵.

裁剪坐标系(clip coordinate system)

在经过$M_{persp\to ortho}$的变换后, 每个顶点均位于裁剪坐标系中. 此时每个点仍然是以齐次坐标系$(x_c,y_c,z_c,w_c)$的形式表示的. 对每个点进行透视除法(即将$xyz$分量都除以$w$), 便可将其转换为标准化设备坐标系(normalized device coordinates). 但在进行透视除法前, 我们需要对顶点进行裁剪. 对于任意$xyz$坐标, 如果其绝对值大于$|w|$, 则这个顶点需要被裁剪.

在裁剪后并执行透视除法后, 每个点便被变换到了$[-1,1]^3$的空间中. 且这里的每个点都曾位于透视投影定义的截头锥体内.

视口变换(Viewport Transformation)

屏幕像素有一些约定, 见下图

视口变换是这一系列变换的最后一步, 它的作用是将上一步正交投影得到的$[-1,1] \times [-1,1]$的画面转换到显示器$[0,width] \times [0,height]$的画面. 由于像素$(x,y)$在屏幕上的实际坐标是$(x+0.5,y+0.5)$, 最终我们要将正交投影的矩阵转化为$[-0.5,width - 0.5] \times [-0.5,height - 0.5]$ 易得转移矩阵如下: \[ M_{viewport} = \left( \begin{array}{cccc} \frac{width}{2} & 0 & 0 & \frac{width-1}{2} \\ 0 & \frac{height}{2} & 0 & \frac{height-1}{2}\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{array} \right) \] 注意在上面的变换中, 我们保留了$z$的值供z-buffer使用.

综上所述, 一个完整的变换矩阵为: \[ M = M_{viewport}M_{persp}M_{view}M_{model} \]

光栅化(Rasterizing)

光栅化时处理的一般是三角形, 因为它有如下优秀的性质:

将一个三角形变为屏幕上的像素便为光栅化:

在光栅化时, 先求出三角形的边界, 再遍历边界里的每个像素, 判断其是否在三角形的内部(可以用叉积或者重心坐标), 绘制所有在三角形内部的像素即可.

抗锯齿(Antialiasing)理论知识

本节相当多的内容引用自Wikipedia.

混叠（英语：Aliasing），在信号频谱上可称作叠频；在影像上可称作叠影，主要来自于对连续时间信号作取样以数字化时，取样频率低于两倍奈奎斯特频率。

在统计、信号处理和相关领域中，混叠是指取样信号被还原成连续信号时产生彼此交叠而失真的现象。当混叠发生时，原始信号无法从取样信号还原。而混叠可能发生在时域上，称做时间混叠，或是发生在频域上，被称作空间混叠。

在视觉影像的模拟数字转换或音乐信号领域，混叠都是相当重要的议题。因为在做模拟-数字转换时若取样频率选取不当将造成高频信号和低频信号混叠在一起，因此无法完美地重建出原始的信号。为了避免此情形发生，取样前必须先做滤波的操作。

先对原图像进行模糊处理(滤波), 再进行采样, 可以得到较好的抗锯齿效果:

但是如果先采样再进行模糊, 只会得到模糊的锯齿:

下面解释原因.(日后再补)

频域(Frequency Domain)

频率（frequency）又称周率，是物理学上描述某具规律周期性的现象或事件，在每单位时间内（即每秒）重复发生的次数.

设$\tau$时间内某事件重复发生$n$次, 则该事件发生的频率$f$为 \[ f = \frac{n}{\tau} \text{Hz} \] 又因为周期定义为重复事件发生的最小时间间隔，故频率也可以表示为周期（$T$）的倒数： \[ f = \frac{1}{T} \text{Hz} \]

MSAA

MSAA先进行超采样, 即原来采样一次的单个像素, 现在采样$n \times n$次得到$n^2$个小像素. 超采样后, 将每个像素对应的在三角形内部的小像素进行平均(这一步相当于模糊操作), 即可得到最终的采样结果.

Z-buffer

着色(Shading)

Bling-Phong模型

如图, 现实世界中的光线可以分为三个部分:

反射高光(Specular highlights)
漫反射(Diffuse reflection)
环境光(Ambient lighting)

分别对这三种光进行模拟, 便可得到较为真实的效果.

在计算光线反射时, 我们是对每个点(Shading point)进行单独计算的. 尽管待计算的点可能位于曲面上, 我们仍然认为极小的局部是一个平面. 待处理的输入参数有以下几个:

观测方向$\hat v$(Viewer direction)
平面法线$\hat n$(Surface normal)
光线方向$\hat l$(Light direction), 它从表面指向光源, 与光线的方向相反.
表面参数(Surface parameters), 包括颜色, 反射率等

注意: 以上的向量均只表示方向, 为单位向量. 另外着色只考虑自身, 不会考虑其他物体的存在, 要将着色(Shading)和阴影(Shadow)进行区分.

漫反射(Diffuse reflection)

当一束光射向粗糙的物体表面时, 粗糙表面会把入射光向各个方向进行反射, 称为漫反射. 需要注意无论我们的视角角度如何, 漫反射的效果都是一样的. 因此漫反射的光照强度只和入射角、平面法线、平面距光源的距离有关.

首先来考虑入射角:

设$\theta$是$\hat l$和$\hat n$的夹角, 那么平面接收到的能量是与$\cos \theta = \hat l \cdot \hat n$成正比的. 这被称为兰伯特余弦定律(Lambert's cosine law).

再考虑距离:

能量与光源的距离成平方反比关系, 即距离光源$r$的点接收到的光强是光源的$\frac{1}{r^2}$. (半径为$r$的球壳表面积为$4\pi r^2$, 结合能量守恒即可推导出).

因此我们可以得出漫反射光强的计算公式, 其中$k_d$为漫反射参数, 一般代表颜色, 可以定义成一个三维的RGB向量. \[ L_d = k_d \cdot \frac{I}{r^2} \cdot \max(0,\hat n \cdot \hat l) \]

反射高光(Specular highlights)

反射高光的光强与我们的观察角度有很大关系. 对于一个比较光滑的物体, 它的反射高光只集中在反射向量($\hat R$)附近的一小块区域(图中黄色部分). 要进行计算, 我们就需要先计算出$\hat R$. 而Bling-Phong模型使用半程向量(half vector)来进行近似计算, 避免了对$\hat R$的计算:

半程向量$\hat h$是位于$\hat l$和$\hat v$夹角的角平分线上的单位向量. \[ \hat h = \frac{\hat v + \hat l}{|\hat v + \hat l|} \] 假设$\hat h$与$\hat n$的夹角$\alpha$近似等于$\hat R$与$\hat v$的夹角, 便可以得出反射光的光强计算公式: \[ L_s = k_s \cdot \frac{I}{r^2} \cdot \max(0,\hat n \cdot \hat h)^p \] 其中$k_s$一般为白色, 也可取光源颜色. 指数$p$可以加速高光的衰减, 让模型更加真实. 这个参数一般取$100 \sim 200$

环境光(Ambient lighting)

环境光非常的简单粗暴: \[ L_a = k_a \cdot I_a \] 把三种光全部加起来, 就能得到Bling-Phong模型的光照: \[ L = L_d + L_s + L_a\\ = k_d \cdot \frac{I}{r^2} \cdot \max(0,\hat n \cdot \hat l) + k_s \cdot \frac{I}{r^2} \cdot \max(0,\hat n \cdot \hat h)^p + k_a \cdot I_a \]

着色频率(待补全)

Bling-Phong模型只给出了着色方法, 但Shading Point是什么并没有给出定义. 本节便讨论不同着色频率的实现和效果.

着色频率可以分为平面着色(Flat shading), Gouraud shading和Phong shading.

在这里先看一下各种着色频率的区别:

平面着色(Flat Shading)

平面着色对每个三角形面求出一个法线, 并据此求出整个三角形的shading结果.

Gouraud Shading

Gouraud着色先求出三角形每个顶点的法线, 计算出颜色, 再对三角形内部的每个点进行插值(interpolation).

在一般的obj模型中我们一般只能获得三角形的法线, 因而每个顶点的法线由与其相邻的三角形法线(加权)平均得来: \[ \vec N_v = \frac{\sum_{i}\vec N_i}{\sum_{i}|\vec N_i|} \]

Phong Shading

Phong shading先插值计算出每个像素的法线, 再对每个像素进行shading.

渲染管线(Graphics Pipeline)

纹理映射(Texture Mapping)

一般默认$u,v \in [0,1]$

重心坐标(Barycentric coordinates)

一般定义

设$V_1,...,V_n$是$n$维向量空间$V$中单纯形的顶点, 对于$V$中任意一点$P$, 有 \[ (\sum_{i = 1}^n \lambda_i)\cdot P = \sum_{i = 1}^n \lambda_i \cdot V_i \] 则系数$\lambda_1,...,\lambda_n$称为$P$关于$V_1,...V_n$的重心坐标. 一般规定$\sum_{i = 1}^n \lambda_i = 1$, 此时称为正规化的重心坐标. 以下讨论的重心坐标均为正规化的.

二维平面中的重心坐标系

定义

二维平面上, 单纯形为三角形. 假设给出一三角形$\Delta ABC$和三角形所在平面上的一点$P$, 则$P$可被如下重心坐标表示: \[ P = \alpha A + \beta B + \gamma C \]

考虑下图(来源见链接):

将$\vec{AB},\vec{AC}$作为坐标系, 可以得到 \[ \begin{align} P =& A + \beta(\vec{AB}) + \gamma(\vec{AC})\\ =& A + \beta(B - A) + \gamma(C - A)\\ =& (1 - \beta - \gamma)A + \beta B + \gamma C \end{align} \] 这样便可以从坐标系的角度来解释重心坐标, 同时能得到 \[ \alpha = 1 - \beta - \gamma \]

求解:解方程组

将点带入重心坐标的定义, 可以得到: \[ \left \{ \begin{array}{cccc} x_p = (1 - \beta - \gamma) x_a + \beta x_b + \gamma x_c \\ y_p = (1 - \beta - \gamma) y_a + \beta y_b + \gamma y_c \end{array} \right. \] 提取出$\beta \space \gamma$: \[ \left \{ \begin{array}{cccc} x_p - x_a = \beta(x_b - x_a) + \gamma (x_c - x_a)\\ y_p - y_a = \beta(y_b - y_a) + \gamma (y_c - y_a) \end{array} \right. \] 转化为矩阵形式: \[ \left[ \begin{array}{} x_b - x_a & x_c - x_a\\ y_b - y_a & y_c - y_a \end{array} \right] \left[ \begin{array}{} \beta \\ \gamma \end{array} \right] = \left[ \begin{array}{} x_p - x_a\\ y_p - y_a \end{array} \right] \] 这样便可求解出$\beta \space \gamma$进而求出$\alpha$

求解:向量叉积

把坐标系的表达式变换一下: \[ \beta \vec{AB} + \gamma \vec{AC} + \vec{PA} = 0 \] 拆开后: \[ \left \{ \begin{array}{cccc} \beta\vec{AB}_x + \gamma \vec{AC}_x + \vec{PA}_x = 0\\ \beta\vec{AB}_y + \gamma \vec{AC}_y + \vec{PA}_y = 0 \end{array} \right. \] 转化为矩阵形式: \[ \left[ \begin{array}{} \beta & \gamma & 1 \end{array} \right] \left[ \begin{array}{} \vec{AB}_x\\ \vec{AC}_x\\ \vec{PA}_x \end{array} \right] = 0\\ \left[ \begin{array}{} \beta & \gamma & 1 \end{array} \right] \left[ \begin{array}{} \vec{AB}_y\\ \vec{AC}_y\\ \vec{PA}_y \end{array} \right] = 0\\ \] 从几何意义考虑, 这相当于$[\beta ,\gamma,1]$这个向量和后两个向量分别垂直. 因而我们求出后两个向量的叉积(假设结果为$[x,y,z]$), 则$[\frac{x}{z}, \frac{y}{z},1]$等于$[\beta,\gamma,1]$.

需要注意$z = 0$的特殊情况, 根据叉乘的定义, 此时有$\vec{AB}_x \cdot \vec{AC}_y - \vec{AB}_y \cdot \vec{AC}_x = 0$, 即$\vec{AB} \cdot \vec{AC} = 0$, 说明$ABC$三点共线. 这时返回一个任意的负重心坐标即可.

求解:三角形面积

如图, 假设三角形面积分别为$A_A,A_B,A_C$. 则 \[ \alpha = \frac{A_A}{A_A+A_B+A_C}\\ \beta = \frac{A_B}{A_A+A_B+A_C}\\ \gamma = \frac{A_C}{A_A+A_B+A_C} \]

应用

判断点与三角形的关系

重心坐标有很多应用, 在这里我们能用到的结论是: 若$\alpha,\beta,\gamma$均大于$0$, 则点$P$在三角形$\Delta ABC$的内部; 若$\alpha,\beta,\gamma$有一个等于$0$, 则点$P$在三角形的边上; 若$\alpha,\beta,\gamma$有两个等于$0$, 则点$P$在三角形的顶点上.

下面的代码采用了向量叉积法:

static bool insideTriangle(int x, int y, const Vector3f* _v)
{
    bool inside = true;
    Vector3f result =  //cal cross
            Vector3f(_v[1].x() - _v[0].x(),_v[2].x() - _v[0].x(),_v[0].x() - x)
            .cross(Vector3f(_v[1].y() - _v[0].y(),_v[2].y() - _v[0].y(),_v[0].y() - y));
    if(fabs(result.z()) < eps)inside = false;//z == 0
    result.x() /= result.z(),result.y() /= result.z();
    float alpha = 1 - result.x() - result.y(),beta = result.x(),gamma = result.y();
    if(alpha < eps or beta < eps or gamma < eps)inside = false;
    return inside;
}

重心

根据顶点进行线性插值

求出重心坐标后, 只需计算 \[ V = \alpha V_A + \beta V_B + \gamma V_C \] 其中$V$可以是任何想要插值的值. (注意投影会改变重心坐标, 在插值时需要先在三维计算重心坐标再投影到二维)

透视矫正插值

阅读资料: https://zhuanlan.zhihu.com/p/403259571

令 \[ Z_n = \frac{1}{\alpha + \beta + \gamma} \] 顶点的$z$坐标和纹理坐标在经过重心插值后, 还需乘上这个系数进行矫正.

纹理放大(Texture Magnification)

通常情况下我们希望像素(pixel)和纹理元素(texel)是一一对应的关系, 但当pixel数量多于/少于texel时, 便需要一种方法来将纹理进行变换.

双线性插值(Bilinear Interpolation)

如下图, 假设我们需要采样红点处的texel值$f(x,y)$

按如下方法选取红点临近的四个texel并将其标号($u$即该texel的值, 如颜色值, 深度等等)

并记红点到$u_{00}$的中心点的距离为$s,t$. 这里假设$s,t \in [0,1]$

我们按如下方法定义一维的线性插值: \[ \text{lerp}(x,v_0,v_1) = v_0 + x (v_1 - v_0) \]

首先从一个方向(这里是水平)进行插值: \[ u_0 = \text{lerp}(s,u_{00},u_{10})\\ u_1 = \text{lerp}(s,u_{01},u_{11}) \] 再从另一个方向(竖直)插值即可得到$f(x,y)$: \[ f(x,y) = lerp(t,u_0,u_1) \] 这样一个pixel便综合考虑了它临近的四个texel的值.

当纹理过大时

如图，当一个pixel覆盖了很多个texel时(地平线处的格子), 仍然直接进行插值便会出现很严重的问题. 在之前的反走样部分, 通过超采样(例如MSAA)可以解决这一问题. 但超采样对性能的消耗很高. 为了解决性能问题, 我们直接求每个pixel覆盖的texel的平均值.

Mipmap

如上图, mipmap从第$0$层(原图)开始, 每次将图像的边长缩小一半, 形成新的一层(记为$D$), 直到不能缩小为止. 这样每一层都存储了一个方形区域texel平均值的近似值. 查询时, 只要根据pixel的坐标和覆盖大小, 在mipmap中找到对应位置, 便可快速求出平均值, 容易计算这样只会占用$\frac{1}{3}$的额外存储空间.

在计算$D$时, 我们取当前待求点所在的一个正方形(上图红色部分), 将其投影到纹理空间中. 令$L$为当前点到相邻两个点在纹理空间中距离的最大值, 则$D = \log_2L$

用这样的方法求出来的$D$将其取整到整数后, 有很大的问题: 不连续(图中的色块).

三线性插值(Trilinear Interpolation)可以解决这个问题. 我们先对$D$和$D+1$层分别做一次双线性插值, 再根据取整之前$D$的值再做一次线性插值, 便能得到不错的结果:

各项异性过滤(Anisotropic Filtering)

如上图.屏幕空间到纹理空间的映射过程中, 并不是每个像素都对应一个理想的正方形, 有的会被拉伸成细长的矩形. Mipmap在处理这样的像素时, 因为平均范围过大变得过度模糊而失去了细节.

各向异性过滤在Mipmap的基础上, 计算了不同长度宽度的组合(Mipmap只计算了上图中对角线部分). 这样在查询时, 对于细长的、近似水平或竖直的矩形可以得到更好的效果, 但仍然难处理细长的斜向矩形.

立方体贴图(Cube Map)

在记录环境光时, 可以用一个光滑的金属球, 记录其表面的反射光, 再将其展开成一个平面. 但这样有一个问题: 展开图的上部和下部会发生比较严重的变形. 立方体贴图便是解决这种问题的方法.

如图, 我们用一个包围盒把球围起来. 从球心出发向外的每一条光线最终会打在这个包围盒上. 记录这些打在包围盒上的光线, 再将包围盒展开, 便得到了立方体贴图. 立方体贴图的变形明显更少, 但是需要计算一下每个方向对应的面.

凹凸/法线贴图(Bump/Normal Mapping)

推荐阅读: https://zhuanlan.zhihu.com/p/412555049

在之前的Bling Phong模型中, 我们利用了顶点的法线数据来插值得到每个像素的法线. 这种方式计算量大, 却不能在模型面数较少时取得较好的效果. 凹凸贴图/法线贴图也是材质的一种, 它让我们可以用贴图的信息来计算每个顶点的法线. 右图便是一个简单的球形应用了法线贴图后的效果. 可以看出,凹凸/法线贴图可以在不增加三角形面的前提下提升表现细节.

凹凸贴图

凹凸贴图是一张二维的黑白图片, 一般越白的部分表示高度越高. 它通过重新定义(或者说扰动)高度的方法, 来让我们能计算新的法线. 在上图中, 物体本身的法线为$\vec{p}$, 经过凹凸贴图扰动后(橙色曲线)计算出的新法线为$\vec{n}$.

先考虑如何在一维求扰动后的法线. 假设原物体是一条和$x$轴平行的直线, 这样它的法线$\vec{n} = (0,1)$. 蓝色曲线是被凹凸贴图扰动后的结果. 我们先近似求出蓝色点切线的斜率$dp$. 因为凹凸贴图也是一种材质, 能细分的最小单位是一个像素, 所以$dp$可以简单的用相邻两点间高度$h$的差值得到(其中系数$c$定义了凹凸贴图的影响有多大): \[ dp = c(h(p+1) - h(p)) \] 这样切线便是$(1,dp)$. 将其旋转90°可得法线$\vec{n} = (-dp,1)\text{.normalized()}$. 时刻注意表示方向的向量必须归一化.

三维情况下, 我们假设最初的法线是$(0,0,1)$. 类比(这里有较为复杂的推导步骤, 在这里省略)一维情况, 分别求出$u,v$两个方向上的斜率, 可得法线: \[ \vec{n} = (-du,-dv,1)\text{.normalized()} \] 注意: 上述坐标是在假设法线为$(0,0,1)$, 且构成坐标系的一个轴的前提下计算的. 要真正使用凹凸贴图定义的法线, 还需将其变换回世界坐标.

法线贴图

凹凸贴图存储了每个像素扰动后的高度, 在使用时仍需进行计算来得到每个像素的法线. 而法线贴图则帮我们完成了这一步, 它直接将法线信息存储在材质中. 材质的$R,G,B$三维便对应了法线的$x,y,z$. 需要注意的是, RGB分量是无符号值, 而法线是有符号的. 因此我们在写入($\text{normal} \to \text{color}$)和读取($\text{color} \to \text{normal}$)时要进行一定的转换: \[ \text{color} = \frac{\text{normal}}{2} + 0.5\\ \text{normal} = 2 \cdot (\text{color} - 0.5) \] 另外值得注意的一点是, 大多数法线贴图看上去都是蓝色的, 这是因为默认法线是$(0,0,1)$, 编码成颜色信息后是$\text{RGB}(0.5,0.5,1)$. 蓝色部分说明对法线信息没有太大改动.

TBN矩阵(待补充)

在使用凹凸贴图/法线贴图时, 我们均假设了法线指向$(0,0,1)$

(等学到光线追踪再补)

位移贴图(Displacement Mapping)

阴影贴图(Shadow mapping)

阴影贴图的工作流程如下:

从光源出发(将光源当作相机), 生成一张物体的深度图
从相机出发，对于看到的每一个点, 将其投影回上一步的深度图对应的位置. 如果该点的实际深度与深度图记录的相同, 那么这个点就是可见的, 否则不可见.

Shadow Mapping也会存在一些问题:

硬阴影，且只支持点光源
质量与深度图的分辨率有关
存在数值精度问题 Z-fighting

几何(Geometry)

基本表示方法

隐式表述

普通方程

普通方程直接给出$x,y,z$之间的关系, 例如$f(x,y,z) = 0$. 观察上图的公式, 很难直接看出普通方程对应的图形是什么, 即普通方程的采样十分困难. 但判断某个点在普通方程定义的几何体的内部/外部非常简单: 将该点的坐标带入$f$, 与等式右边的值比较即可. 同时, 隐式表述的几何体很容易做光线和几何的交.

构造实体几何(Constructive Solid Geometry)

构造实体几何(CSG)使用逻辑运算来将简单的几何体组合成复杂的形体.

距离函数(Distance Functions)

水平集方法(Level Set Methods)

分形(Fractals)

显式表述

参数方程

显式表述的参数方程直接定义一个从一个坐标向另一个坐标转换的方式. 它采样起来十分简单: 只需将原坐标的值代入$f$就能得到新坐标的值. 但它判断某个点在几何体的内部/外部十分困难.

点云(Point Cloud)

多边形面(Polygon Mesh)

贝塞尔曲线

见曲线章节

曲线

贝塞尔曲线(Bézier Curves)

贝塞尔曲线是一种由控制点来确定形状的曲线. $n$个控制点描述了一个$n-1$阶的贝塞尔曲线. 下面是正在绘制的不同阶贝塞尔曲线的一些例子:

注意图中的$t = 0.49$. 遍历$t \in[0,1]$, 计算并绘制某一点的运动轨迹(图中红点), 这条轨迹便是贝塞尔曲线. 从上图可以看出, $k$阶贝塞尔曲线可以从$k-1$阶递推得到. 这一递推算法称作德卡斯特里奥算法(De Casteljau's algorithm).

在解释这一算法之前, 我们先来看一下一阶贝塞尔曲线的表达式: \[ B_1(t) = P_0 + (P_1 - P_0)t \] 很显然, 这是一个线性插值. 观察一阶贝塞尔曲线的图像, 可以看出它就是一条直线.

在时刻为$t$时, 图中红色部分与整个线段的长度比为$t:1$. 现在我们扩展到二阶贝塞尔曲线:

可以看到, 我们对$\vec{b_0b_1}$和$\vec{b_1b_2}$分别进行线性插值(简称lerp)得到了$b_0^1$和$b_0^2$两个点, 这两个点可以作为一个新的一阶贝塞尔曲线的控制点. 再对其进行线性插值, 便得到了要绘制的点$b_0^2$.

三阶曲线也类似. 这样, 对于$n$阶的贝塞尔曲线, 每插值一轮, 其阶数便减$1$, 直到阶数等于$1$. 这便是德卡斯特里奥算法的核心思想.

贝塞尔曲线也有非递推的公式($P_i$为控制点): \[ b^n(t) = \sum_{i = 0}^nP_i \cdot B_i^n(t), t\in[0,1] \] 这个方法更快, 但数值稳定性不如递推算法高. 其中$B_i^n(t)$称为伯恩施坦多项式 (Bernstein polynomial) \[ B_i^n(t) = \binom{n}{i}t^i(1-t)^{n-i} \]

伯恩施坦多项式有一个概率意义上的解释: 某独立事件发生一次的概率为$p$, 那么进行$n$次事件, 恰好发生$i$次的概率为$B_i^n(p)$. 既然我们能将其解释为概率, 那么同一阶伯恩施坦多项式在任何时间下的和必为$1$, 即$\sum_{i = 0}^nB_i^n(t) = 1$.

贝塞尔曲线有如下性质:

一条贝塞尔曲线一定从第一个控制点开始, 在最后一个控制点结束. 需要注意的是贝塞尔曲线不必经过全部控制点.
对于三阶贝塞尔曲线, 其插值得到的两个二阶曲线为: $b'(0) = 3(P_1 - P_0),b'(1) = 3(P_3 - P_2)$.
对控制点做仿射变换等价于对贝塞尔曲线做仿射变换, 但对投影变换不满足.
贝塞尔曲线一定位于控制点构成的凸包内.

分段贝塞尔曲线(Piecewise Bézier Curves)

高阶的贝塞尔曲线看上去很不直观, 也很难再用控制点来控制它的形状.

将控制点每$k$个点(比如每四个)分段, 便是分段贝塞尔曲线. 但直接进行分段可能会出现尖锐的部分(如上图). 对贝塞尔曲线的连通性有如下定义:

$C^0$连续性: 第一段曲线的终点和第二段曲线的起点相同.

$C^1$连续性: 两曲线的交点及其相邻的左右两个点共线, 且交点是这一线段的中点.

样条(splines)

B样条(待补充)

曲面

贝塞尔曲面(Bézier Surfaces)

我们用类似双线性插值的思想, 将$4\times 4$个控制点分成四组, 每组四个点, 求出每一组在时间$t_1$时所绘制的点. 再将这四个点作为一条新的贝塞尔曲线的控制点, 遍历$t_2$. 这样在$t_1 \cdot t_2$时间内, 绘制的点便能形成一个贝塞尔曲面.

对网格面的几何操作

曲面细分(Mesh subdivision)

表面细分做了两件事: 1. 增加三角形数量 2. 调整三角形的位置

Loop细分(Loop Subdivision)

注意这里的Loop和循环没关系, Loop是发明者的姓.

Loop细分只能处理三角形网格. 它将一个三角形分成四个, 在划分的过程中会产生新的点. 我们需要将新生成的点和原先存在的点区分进行处理:

对于新生成的每个点, 记这个点所在的三角形边上两个顶点为$A,B$. 这条边相对的顶点为$C,D$. 则新生成的顶点$P_{new}$位置调整到: \[ P_{new} = \frac{3}{8} \cdot(A + B) + \frac{1}{8} \cdot (C + D) \]

对于旧顶点, 按下式更新: \[ P_{old} = (1 - n\cdot u) \cdot P_{origin} + u\cdot \sum P_{neighbor} \] 其中:

$P_{origin}$: 该点原先的位置
$\sum P_{neighbor}$: 该点相邻点的和
$n$: 该点的度数(相邻点的数量)
$u = \frac{3}{8n}$. 若$n = 3, u = \frac{3}{16}$

Catmull-Clark Subdivision (General Mesh)

Catmull-Clark Subdivision可以处理一般化的网格. 它将每个面分为四边形面和非四边形面两部分. 将点按度数分类: 所有度数不是$4$的点为奇异点(Extraordinary vertex).

在细分时, 我们在每个面内选取一个点, 再将其和每条边的中点相连. 可以看出, 做了一次细分之后, 新增的奇异点个数等于原来非四边形面的个数, 且这些非四边形面都会变成四边形面. 也就是说, 奇异点个数仅可能在第一次细分时增加.

计算方式见上图(待详细补充).

曲面简化(Mesh Simplification)

二次误差度量(Quadric Error Metrics)

(待补充)

我们在删掉一个顶点之后, 需要调整其他顶点的位置来让简化后的模型与原模型误差尽量的小. 这个误差可以用二次误差度量来表示.

光线追踪(Ray Tracing)

我们对光线做如下假设:

光线沿直线传播
光线不会发生碰撞
光线从光源出发, 经过一系列折射反射等后, 进入观察者的眼睛, 这个过程可逆.(reciprocity)

光线投射算法(Ray Casting)

Recursive (Whitted-Style) Ray Tracing

如上图. 对于每个像素, 从观测点发出一条光线. 这一条光线可以经过多次反射、折射、吸收, 在多个地方与多个物体有交点. 我们对这每一个交点将其与光源连线, 判断可见性并着色, 再将着色结果全部加回到像素中. 这便是Recursive (Whitted-Style) Ray Tracing的基本思想.

光线与表面求交(Ray-Surface Intersection)

光线定义

与隐式表面求交

先看如何与球求交

球的定义: \[ \bold p : (\bold p - \bold c)^2 - R^2 = 0 \]

上图$\bold o$为光源位置, $\bold d$为光线方向. $\bold p$为光线与球的交点, $\bold c$为球心, $R$为球的半径.

若光线与球相交, 则交点必须满足光线和球的表达式(即$\bold p = \bold o + t\bold d$). 代入得: \[ (\bold o + t \bold d - \bold c)^2 - R^2 = 0 \]

直接展开求解即可. 注意$t$必须为非负实数才有意义.

隐式表面: \[ \bold p: f(\bold p) = 0 \] 则光线与隐式表面求交: \[ f(\bold o + t\bold d) = 0 \] 求根即可.

与三角面(显式表面)求交

可以将光线与三角形求交拆成两步:

求光线与三角形所在平面的交点
判断交点是否在三角形内部

平面可以由一条法线$\vec n$和屏幕上一点$\bold{p'}$来定义. \[ \bold p : (\bold p - \bold p') \cdot \vec n = 0 \] 展开为$ax + by + cz + d = 0$的形式.

将光线定义代入:

求出$t$后判断是否在三角形内部即可.

Möller Trumbore Algorithm

回忆一下我们用重心坐标判断点与三角形关系的过程: 若$\alpha,\beta,\gamma$均大于$0$, 则点$P$在三角形$\Delta ABC$的内部. 那么可以写出如下式子: \[ \bold o + t\bold d = (1 - b_1 - b_2)\bold{P_0} + b_1 \bold{P_1} + b_2 \bold{P_2} \]

这是一个有三个变量$(t,b_1,b_2)$三个方程的线性方程组. 很容易求解.

作业5(待补充)

参考: https://zhuanlan.zhihu.com/p/431092843

生成每个像素对应的光线

求交过程中的加速

包围盒(Bounding Volumes)

在之前的算法中, 我们对每个像素发出的光线都要和全部三角形面求交. 使用包围盒后, 便可以先和包围盒求交. 如果二者不相交, 包围盒内的所有三角形都可以忽略不计.

我们可以把包围盒理解为三对平面的交集. 这三对平面一般和坐标轴构成的平面平行. 下面主要讨论如何与轴对齐包围盒求交:

先考虑二维的情况. 先将光线与$x$面求交, 可以得到光线进入和射出的时间(这两个时间可以求出一条线段), $y$面同理. 将这两个线段求交(即$t_{min}$取$\text{max}$,$t_{max}$取$\text{min}$), 即可求出光线在包围盒内的部分. 可以看出, 光线在包围盒内当且仅当光线进入矩形的所有对面($x$和$y$). 光线离开包围盒只需满足离开任意对面.

现在推广到三维情况: 我们对三组对面分别求出它们的$t_{min}$和$t_{max}$.

可得进入时间: $t_{enter} = \max(t_{min})$

离开时间: $t_{exit} = \min(t_{max})$

当$t_{enter} \leq t_{exit}$且$t_{exit} \ge 0$时, 光线与包围盒相交.

对于AABB, 求交十分简单

空间分割(Uniform Spatial Partitions (Grids))

在之前的优化中, 我们通过包围盒来避免了每条光线与每个三角面进行求交. 空间分割是对每个包围盒内部的判断进行进一步优化.

在一个包围盒内(图中最外层的边界即为包围盒), 将这个包围盒分成$n \times n$个小格子. 如果一个格子包含了物体的表面, 就将其进行标记. 对于一条打进包围盒内部的光线, 我们首先计算它的行进路径上经过了哪些格子(Bresenham直线算法), 再对经过的格子中被标记的格子进行光线与物体的求交. 这样便避免了将光线与包围盒内的每个物体进行求交运算. 在三维空间中, 格子的个数一般取$27 \times 物体个数$

这种方法适合物体较均匀分布的场景.

空间划分(Spatial Partitions)

K-D Tree

https://oi-wiki.org/ds/kdt/

K-D Tree有一些问题. 其一是物体与包围盒之间的关系很难进行判断(例如包围盒被三角形完全包围的情况). 其二是一个物体可能会同时存在于多个叶子节点中. 目前实践上已经很少使用K-D Tree, 而是改用BVH.

物体划分与包围盒层次结构(Object Partitions & Bounding Volume Hierarchy (BVH))

PPT的信息已经足够充分, 这里不做过多解释. 详细实现留在作业6中.

辐射度量学(Radiometry)

基础

辐射能和通量(Radiant Energy and Flux (Power))

辐射能即电磁辐射的能量, 用字母$Q$表示. 单位是焦耳($J$,Joule). \[ Q [\text J = \text{Joule}] \]

辐射通量(Radiant flux), 也称作辐射功率(Radiant power), 是单位时间内发射/反射/传送/接受等等的能量. 用符号$\Phi$表示. 单位是瓦或者流明. \[ \Phi \equiv \frac{\text dQ}{\text d t}[\text{W = Watt}][\text{lm = lumen}] \]

上面是衡量光线的几种概念.

Radiant Intensity: 辐射强度

irradiance: 辐照度

radiance: 辐射

(中文译名不重要)

辐射强度(Radiant Intensity,I)

\[ I(\omega) = \frac{\text d \Phi}{\text d \omega} [\text{candela}] \]

立体角(Solid Angles)

在二维平面上, 角度$\theta$被定义为扇形对应弧长$l$与半径$r$的比值. 在三维空间上, 令锥体所对的球面面积为$A$, 则立体角: \[ \Omega = \frac{A}{r^2} \] 单位是球面度(steradians). 注意立体角和角度一样, 也是无量纲量. 球的弧面度是$4\pi$