ransac算法思想
Ransac: Random Sample Consensus, 隨機抽樣一致性。RANSAC算法在1981年由Fischler和Bolles首次提出。
Ransac是一種通過使用觀測到的數(shù)據(jù)點來估計數(shù)學(xué)模型參數(shù)的迭代方法。其中數(shù)據(jù)點包括內(nèi)點(inlier),外點(outlier)。outlier對模型的估計沒有價值,因此該方法也可以叫做outlier檢測方法。其中inliners指樣本中正確的的樣本點,常稱為內(nèi)點,內(nèi)群點,正確樣本,outliners指樣本中的噪聲點(如測量時誤差很大的點,數(shù)據(jù)值不準(zhǔn)確的點等),常稱為外點,離群點,噪聲點。
所以通俗點將,就是我們有一組觀測數(shù)據(jù),但是數(shù)據(jù)點中有很多噪聲點,如果我們直接用模型來擬合數(shù)據(jù),模型肯定會受噪聲點影響,而且噪聲點比例越大,擬合出的模型也會越不準(zhǔn)確。Ransac就是用來解決樣本中的外點問題(噪聲),最多可處理50%的外點情況。
1. Ransac的假設(shè)和思想
再深入點討論,有一組觀測數(shù)據(jù),要建立一個模型來擬合數(shù)據(jù),我們第一步要做的,肯定是用一個標(biāo)準(zhǔn)對數(shù)據(jù)進行篩選,去除噪聲點,讓數(shù)據(jù)盡量干凈和準(zhǔn)確。但是如果我們沒有一個合適的篩選標(biāo)準(zhǔn),該怎么辦呢?我們可以假設(shè):觀測數(shù)據(jù)中除了外點(噪聲),肯定存在內(nèi)點(準(zhǔn)確點),而且內(nèi)點的個數(shù)夠我們用來擬合一個模型(比如擬合一條直線,至少要兩個數(shù)據(jù)點)。有了這個假設(shè),我們可以從觀測數(shù)據(jù)中隨機挑選出n個數(shù)據(jù)點,用這n個點來擬合模型,然后對這個模型進行評價。如果通過迭代,反復(fù)重復(fù)這個過程,只要迭代次數(shù)過大,我們隨機挑選出來的n個樣本點是有可能全部是內(nèi)點的,再配合上模型評價,就能找到最優(yōu)的擬合模型。這就是ransac的主要假設(shè)和思想,官方敘述如下:
RANSAC的基本假設(shè)是 “內(nèi)群”數(shù)據(jù)可以通過幾組模型參數(shù)來敘述其數(shù)據(jù)分布,而“離群”數(shù)據(jù)則 是不適合模型化的數(shù)據(jù)。 數(shù)據(jù)會受噪聲影響, 噪聲指的是離群:例如從極端的噪聲, 錯誤解釋有關(guān)數(shù)據(jù)的測量, 或不正確的假設(shè)。RANSAC假定,給定一組(通常很小的)內(nèi)群,存在一個 程序,這個程序可以估算最佳解釋或最適用于這一數(shù)據(jù)模型的參數(shù)。
宏觀點來看,Ransac是一種思想,一種利用迭代來增加概率的方法,即:
RANSAC是一種思想,一個求解已知模型的參數(shù)的框架。它不限定某一特定的問題,可以是計算機視覺的問題,同樣也可以是統(tǒng)計數(shù)學(xué),甚至可以是經(jīng)濟學(xué)領(lǐng)域的模型參數(shù)估計問題。
RANSAC是一種迭代的方法,用來在一組包含離群的被觀測數(shù)據(jù)中估算出數(shù)學(xué)模型的參數(shù)。 RANSAC 是一個非確定性算法,在某種意義上說,它會產(chǎn)生一個在一定概率下合理的結(jié)果,其允許使用更多次的迭代來使其概率增加。
2. Ransac基本流程和參數(shù)設(shè)置
2.1 Ransac步驟
RANSAC算法的輸入:
1.一組觀測數(shù)據(jù)(往往含有較大的噪聲或無效點
2.一個用于解釋觀測數(shù)據(jù)的參數(shù)化模型
3.一些可信的參數(shù)
迭代過程:
在數(shù)據(jù)中隨機選擇n個點設(shè)定為內(nèi)群
計算適合內(nèi)群的模型,如線性直線模型y=ax+b
把其它剛才沒選到的點帶入剛才建立的模型中,計算是否為內(nèi)群點
記下內(nèi)群數(shù)量
重復(fù)以上步驟, 迭代k次
比較哪次計算中內(nèi)群數(shù)量最多,內(nèi)群最多的那次所建的模型就是我們所要求的解
注意:不同問題對應(yīng)的數(shù)學(xué)模型不同,因此在計算模型參數(shù)時方法必定不同,RANSAC的作用不在于 計算模型參數(shù)。(這導(dǎo)致ransac的缺點在于要求數(shù)學(xué)模型已知)
2.2 Ransac參數(shù)設(shè)置
在上面的步驟中,有兩個參數(shù)比較重要,需要提前確定,即:
一開始的時候我們要隨機選擇多少點(n)
以及要重復(fù)迭代多少次(k)
假設(shè)每個點是真正內(nèi)群的概率為 w:
[w = 內(nèi)群的數(shù)目/(內(nèi)群數(shù)目+外群數(shù)目)
]
通常我們不知道 w 是多少, (w^n)是所選擇的n個點都是內(nèi)群的機率, (1-w^n)是所選擇的n個點至少有一個不是內(nèi)群的機率, ((1 ? w^n)^k)是表示重復(fù) k 次都沒有全部的n個點都是內(nèi)群的機率, 假設(shè)算法跑 k 次以后成功的機率是p,那么
[1 ? p = (1 ? w^n)^k\
p = 1 ? (1 ? w^n)^k
]
我們可以通過P反算得到抽取次數(shù)K:(K=log(1-P)/log(1-w^n)。)
所以如果希望成功機率高:
當(dāng)n不變時,k越大,則p越大;
當(dāng)w不變時,n越大,所需的k就越大。
通常w未知,所以n選小一點比較好。
實際設(shè)置經(jīng)驗
可以先設(shè)置好P和n,然后K設(shè)置為一個很大的數(shù),然后在迭代的過程中,每次迭代完都可以計算一次w,然后對K進行動態(tài)修改。 opencv的findHomography函數(shù)中ransac就是采用的這一策略:
H = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5) # 采用Ransac方法計算單應(yīng)性矩陣
參考:https://blog.csdn.net/laobai1015/article/details/51683076/
3. Ransac應(yīng)用
在計算機視覺方面,Ransac最主要的應(yīng)用主要包括直線擬合,平面擬合,單應(yīng)性矩陣擬合等。這里簡單介紹下直線擬合和透視矩陣擬合。
3.1 Ransac直線擬合
假設(shè)我們知道兩個變量X與Y之間呈線性關(guān)系,Y=aX+b,我們想確定參數(shù)a與b的具體值。通過實驗, 可以得到一組X與Y的測試值。雖然理論上兩個未知數(shù)的方程只需要兩組值即可確認(rèn),但由于系統(tǒng)誤差的原因,任意取兩點算出的a與b的值都不盡相同。通常情況下,我們可以采用最小二乘法擬合出直線方程(最小二乘法:通過計算最小均方差關(guān)于參數(shù)a、b的偏導(dǎo)數(shù)為零時的值)。但是最小二乘法只適合于誤差較小的情況,如果測量數(shù)據(jù)中外點較多,誤差很大,就需要采用Ransac算法。
在下圖中分別是Ransac和最小二乘法擬合的直線,可以看出兩者的差別。直接采用最小二乘法擬合直線,直線會受離群點影響,偏離中間的內(nèi)群數(shù)據(jù)點,而ransac只是挑選部分?jǐn)?shù)據(jù)點進行直線擬合,所以不會受離群點影響。
下面是一段示例代碼,比較了ransac和最小二乘法擬合直線的區(qū)別:
#include<iostream>
#include<opencv2/opencv.hpp>
#include<vector>
#include<set>
using namespace std;
bool ransacFitLine(vector<cv::Point2d>& points, cv::Vec4f& bestParam, vector<cv::Point2d>& inlinerPoints, int initInlinerNums, int iter, double thresh, int EndInlinerNums){
/*
@param points: 所有的樣本點.
@param bestParam: 擬合最好的直線方程
@param inlinerPoints: 擬合最好時的內(nèi)點集合
@param initInlinerNums: 用來擬合直線的初始化內(nèi)點個數(shù)
@param iter: 最大迭代次數(shù)
@param thresh: 閾值,用來判斷某個點是否在擬合的直線上
@param EndInlinerNums: 閾值,內(nèi)點個數(shù)超過閾值時停止迭代
*/
set<int> randomIndex;
vector<cv::Point2d> randomPoints;
vector<cv::Point2d> tempInliners;
int n = points.size();
double minErr=DBL_MAX;
bool flag = false;
while(iter>=0){
randomIndex.clear();
randomPoints.clear();
tempInliners.clear();
// 隨機選取initInlinerNums個點,擬合直線,其他點作為測試點
for(int i=0; i<initInlinerNums; i++){
int index = rand()%n;
randomIndex.insert(index);
randomPoints.push_back(points[index]);
}
// 采用隨機點擬合直線
cv::Vec4f lineParam;
cv::fitLine(randomPoints, lineParam, cv::DIST_L2, 0, 0.01, 0.01);
double linek = lineParam[1]/lineParam[0];
// 計算測試點是否也在擬合直線上, 在直線上的點當(dāng)作inliner
for(int i=0; i<n; i++){
if(randomIndex.find(i)==randomIndex.end()){
double err = linek*(points[i].x - lineParam[2]) + lineParam[3] - points[i].y;
err = err*err; // 差值的平方作為誤差
if (err < thresh){
tempInliners.push_back(points[i]);
}
}
}
// 當(dāng)前inliner內(nèi)點個數(shù)超過閾值時
if(tempInliners.size()>EndInlinerNums){
//采用所有的內(nèi)點,擬合直線并計算平均誤差
tempInliners.insert(tempInliners.end(), randomPoints.begin(), randomPoints.end());
cv::fitLine(randomPoints, lineParam, cv::DIST_L2, 0, 0.01, 0.01);
double linek = lineParam[1]/lineParam[0];
double sumErr = 0;
for(int i=0; i<tempInliners.size(); i++){
double err = linek*(points[i].x - lineParam[2]) + lineParam[3] - points[i].y;
sumErr += err*err; // 差值的平方作為誤差
}
sumErr = sumErr/tempInliners.size(); // 誤差平方和的平均值
if(minErr>sumErr){ // 記錄最小的平均誤差,
minErr = sumErr;
bestParam = lineParam;
inlinerPoints.assign(tempInliners.begin(), tempInliners.end());
flag = true;
}
}
iter--;
}
if(!flag) cout << "did't meet fit acceptance criteria" << endl;
return flag;
}
int main(int argc, char* argv[]){
cv::RNG rng(200);//seed為200
vector<cv::Point2d> points;
// 生成500個隨機樣本點
int nSamples = 500; // 生成樣本數(shù)據(jù)個數(shù)
double k = 1 + rng.gaussian(1); // 0-2之間的隨機斜率
cout << k << endl;
for(int i=0; i<nSamples; i++){
double x = rng.uniform(40, 600);
double y = k*x + rng.uniform(-30, 30); // 加入隨機噪聲
x += rng.uniform(-30, 30); // 加入隨機噪聲
points.push_back(cv::Point2d(x, y));
}
// 加入outlier
for(int j=0; j< 100; j++){
double x = rng.uniform(200, 600);
double y = rng.uniform(10, 300);
points.push_back(cv::Point2d(x, y));
}
cv::Mat img = cv::Mat(720, 720, CV_8UC3, cv::Scalar(255,255,255));
//繪制point
for(int i=0; i<points.size(); i++){
cv::circle(img, points[i], 3, cv::Scalar(0, 0, 0), 2);
}
// 最小二乘法
if(1){
cv::Vec4f lineParam;
fitLine(points, lineParam, cv::DIST_L2, 0.0, 0.01, 0.01);
//獲取點斜式的點和斜率, y=k(x-x0) + y0
cv::Point point0;
point0.x = lineParam[2];
point0.y = lineParam[3];
double linek = lineParam[1] / lineParam[0];
//繪制最小二乘法擬合的直線
cv::Point2d point1, point2;
point1.x=10;
point1.y = linek*(point1.x-point0.x) + point0.y;
point2.x = 600;
point2.y = linek*(point2.x-point0.x) + point0.y;
cv::line(img, point1, point2, cv::Scalar(0, 255, 0), 4);
}
// ransac擬合直線
if(1){ // 為了避免命名沖突,所以放在if局部范圍內(nèi)
cv::Vec4f lineParam;
vector<cv::Point2d> inlinerPoints;
bool ret = ransacFitLine(points, lineParam, inlinerPoints, 50, 1000, 7e3, 300);
if(ret){
cout << ret << endl;
cout << lineParam << endl;
cout << inlinerPoints.size() << endl;
//獲取點斜式的點和斜率, y=k(x-x0) + y0
cv::Point point0;
point0.x = lineParam[2];
point0.y = lineParam[3];
double linek = lineParam[1] / lineParam[0];
//繪制擬合的直線
cv::Point2d point1, point2;
point1.x=10;
point1.y = linek*(point1.x-point0.x) + point0.y;
point2.x = 600;
point2.y = linek*(point2.x-point0.x) + point0.y;
cv::line(img, point1, point2, cv::Scalar(0, 0, 255), 4);
}else{
cout << "Fitting failed, try to change parameters for ransac!!!" << endl;
}
}
cv::imshow("img", img);
cv::waitKey(0);
cv::destroyAllWindows();
return 0;
}
代碼結(jié)果如下:
3.2 Ransac單應(yīng)性矩陣擬合
關(guān)于單應(yīng)矩陣?yán)斫猓篽ttps://zhuanlan.zhihu.com/p/49435367
在進行圖像匹配,全景拼接等時,常會用到單應(yīng)性矩陣。例如在做全景拼接時,對于同一個場景,相機在不同角度拍攝了兩張照片,一般先尋找兩幅圖片的匹配特征點,然后通過匹配特征點的對應(yīng)關(guān)系計算出一個矩陣,這個矩陣就是單應(yīng)性矩陣,利用這個矩陣就能將兩張圖片組合在一起。所以單應(yīng)矩陣描述的是兩組坐標(biāo)點的關(guān)系,即:
[egin{bmatrix}
x_0 \
y_0 \
1
end{bmatrix} =Hegin{bmatrix}
X'_0 \
Y'_0 \
1
end{bmatrix} \
]
上面H即為單應(yīng)性矩陣,((x_0, y_0))為第一幅圖片上的坐標(biāo)點,((X'_0, Y'_0))為第二幅圖片上的坐標(biāo)點,兩個點是一組匹配點。H一般為3x3的矩陣,形式如下:
[H = egin{bmatrix}
a_{11} & a_{12} & a_{13} \
a_{21} & a_{22} & a_{23} \
a_{31} & a_{32} & a_{33} \
end{bmatrix}
]
這個矩陣有9個未知數(shù)(自由度為8),為了求解矩陣我們一般令(a_{33}=1), 則一組匹配點可以列如下兩個方程:
若有四組匹配點,就能列出8個方程,便能計算出單應(yīng)矩陣H, 如下圖所示:
opencv的getPerspectiveTransform函數(shù)就是通過四組匹配點來計算單應(yīng)性矩陣(投影矩陣/透視矩陣),如下:
H = cv2.getPerspectiveTransform(src, dst) # 返回3*3的矩陣
參數(shù):
src:原圖像中的四組坐標(biāo),如 np.float32([[56,65],[368,52],[28,387],[389,390]])
dst: 轉(zhuǎn)換后的對應(yīng)四組坐標(biāo),如np.float32([[0,0],[300,0],[0,300],[300,300]])
如果我們只有四組匹配點,而且能夠保證匹配點都是正確的,的確可以采用上述方法求解單應(yīng)矩陣,但是如果我們有100組匹配點,而且有些匹配點可能是錯誤的,該怎么求解一個擬合最好的矩陣呢?
仔細(xì)觀察上面八個方程的形式,發(fā)現(xiàn)是一個矩陣求解的問題:
[AH=B
]
其中,A是8x8的矩陣,H是8x1的矩陣,B也是8x1的矩陣。若有100組匹配點,則A是200x8的矩陣,B是200x1的矩陣。
聯(lián)想到上面最小二乘法擬合直線的問題:(y=wx),那么在這里x就是矩陣A的每一行(8x1的向量),(w)就是矩陣H,(y)
就是矩陣B的每一行(1x1), 如下圖所示:
所以100組匹配點就是對應(yīng)200條數(shù)據(jù),采用最小二乘法進行擬合(不算是直線了,算平面擬合了)。
有了求解方法,另一個問題就是:100組匹配點中可能存在錯誤的匹配,很明顯錯誤的匹配就是所謂的outlier了。很自然的能想到:采用ransac代替上面的最小二乘法即可。
上述就是采用Ransac擬合單應(yīng)性矩陣的思想,在opencv中findHomography函數(shù)實現(xiàn)了上述過程,如下所示:
#findHomography(srcPoints,dstPoints,method=0,ransacReprojThreshold=3,mask=noArray(),maxIters=2000,confidence=0.995)
H = cv2.findHomography(src_pts, dst_pts, method=cv2.RANSAC, ransacReprojThreshold=5) # 生成變換矩陣
參數(shù):
src_pts:原圖像中的坐標(biāo)(大于等于四組)
dst_pts: 轉(zhuǎn)換后的對應(yīng)坐標(biāo)(大于等于四組)
method:求解單應(yīng)矩陣的方法,取值如下:
0:最小二乘法
RANSAC:RANSAC擬合方法
LMEDS: Least-Median robust method
RHO:PROSAC-based robust method
ransacReprojThreshold:ransac中判斷是否是內(nèi)點的閾值
下面看一個全景拼接的示例,全景拼接的常規(guī)流程 如下:
1、針對某個場景拍攝多張/序列圖像
2、通過匹配特征(sift匹配)計算下一張圖像與上一張圖像之間的匹配特征點
3、對匹配特征點進行篩選,剔除掉部分匹配點。采用ransanc擬合匹配點之間的單應(yīng)性矩陣
3、利用單應(yīng)性矩陣進行圖像映射,將下一張圖像疊加到上一張圖像的坐標(biāo)系中
4、變換后的融合/合成
下面是一個代碼示例:
# coding: utf-8
import numpy as np
import cv2
left_img = cv2.imdecode(np.fromfile(r"D:筆記cppimage_processing圖像拼接left.jpg", dtype=np.uint8), 1)
left_img = cv2.resize(left_img, (600, 400))
right_img = cv2.imdecode(np.fromfile(r"D:筆記cppimage_processing圖像拼接ight.jpg", dtype=np.uint8), 1)
right_img = cv2.resize(right_img, (600, 400))
left_gray = cv2.cvtColor(left_img, cv2.COLOR_BGR2GRAY)
right_gray = cv2.cvtColor(right_img, cv2.COLOR_BGR2GRAY)
hessian = 300
surf = cv2.xfeatures2d.SIFT_create(hessian) # 將Hessian Threshold設(shè)置為400,閾值越大能檢測的特征就越少
kp1, des1 = surf.detectAndCompute(left_gray, None) # 查找關(guān)鍵點和描述符
kp2, des2 = surf.detectAndCompute(right_gray, None)
# kp1s = np.float32([kp.pt for kp in kp1])
# kp2s = np.float32([kp.pt for kp in kp2])
# 顯示特征點
img_with_drawKeyPoint_left = cv2.drawKeypoints(left_gray, kp1, np.array([]), flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow("img_with_drawKeyPoint_left", img_with_drawKeyPoint_left)
img_with_drawKeyPoint_right = cv2.drawKeypoints(right_gray, kp2, np.array([]), flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow("img_with_drawKeyPoint_right", img_with_drawKeyPoint_right)
'''
BFMatcher簡稱暴力匹配,意思就是嘗試所有可能匹配,實現(xiàn)最佳匹配。
FlannBasedMatcher簡稱最近鄰近似匹配。
是一種近似匹配方法,并不追求完美!,因此速度更快。
可以調(diào)整FlannBasedMatcher參數(shù)改變匹配精度或改變算法速度。
參考:https://blog.csdn.net/claroja/article/details/83411108
'''
FLANN_INDEX_KDTREE = 0 # 建立FLANN匹配器的參數(shù)
indexParams = dict(algorithm=FLANN_INDEX_KDTREE, trees=5) # 配置索引,密度樹的數(shù)量為5
searchParams = dict(checks=50) # 指定遞歸次數(shù)
# FlannBasedMatcher:是目前最快的特征匹配算法(最近鄰搜索)
flann = cv2.FlannBasedMatcher(indexParams, searchParams) # 建立匹配器
# 參考https://docs.opencv.org/master/db/d39/classcv_1_1DescriptorMatcher.html#a378f35c9b1a5dfa4022839a45cdf0e89
'''
int queryIdx –>是測試圖像的特征點描述符(descriptor)的下標(biāo),同時也是描述符對應(yīng)特征點(keypoint)的下標(biāo)。
int trainIdx –> 是樣本圖像的特征點描述符的下標(biāo),同樣也是相應(yīng)的特征點的下標(biāo)。
int imgIdx –>當(dāng)樣本是多張圖像的話有用。
float distance –>代表這一對匹配的特征點描述符(本質(zhì)是向量)的歐氏距離,數(shù)值越小也就說明兩個特征點越相像。
'''
matches = flann.knnMatch(des1, des2, k=2) # 得出匹配的關(guān)鍵點
good = []
# 提取優(yōu)秀的特征點
for m, n in matches:
if m.distance < 0.7 * n.distance: # 如果第一個鄰近距離比第二個鄰近距離的0.7倍小,則保留
good.append(m)
src_pts = np.array([kp1[m.queryIdx].pt for m in good]) # 查詢圖像的特征描述子索引
dst_pts = np.array([kp2[m.trainIdx].pt for m in good]) # 訓(xùn)練(模板)圖像的特征描述子索引
# findHomography參考https://docs.opencv.org/master/d9/d0c/group__calib3d.html#ga4abc2ece9fab9398f2e560d53c8c9780
# 單應(yīng)矩陣:https://www.cnblogs.com/wangguchangqing/p/8287585.html
H = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5) # 生成變換矩陣
h1, w1 = left_gray.shape[:2]
h2, w2 = right_gray.shape[:2]
shift = np.array([[1.0, 0, w1], [0, 1.0, 0], [0, 0, 1.0]])
# 點積
M = np.dot(shift, H[0]) # 獲取左邊圖像到右邊圖像的投影映射關(guān)系, 向右移動w1?
dst = cv2.warpPerspective(left_img, M, (w1+w2, max(h1, h2))) # 透視變換,新圖像可容納完整的兩幅圖
cv2.imshow('left_img', dst) # 顯示,第一幅圖已在標(biāo)準(zhǔn)位置
dst[0:h2, w1:w1+w2] = right_img # 將第二幅圖放在右側(cè)
# cv2.imwrite('tiled.jpg',dst_corners)
cv2.imshow('total_img', dst)
cv2.imshow('leftgray', left_img)
cv2.imshow('rightgray', right_img)
cv2.waitKey(0)
cv2.destroyAllWindows()
效果如下:
參考:
https://blog.csdn.net/hzwwpgmwy/article/details/83578694
https://zhuanlan.zhihu.com/p/82793501
4 Ransac優(yōu)缺點
優(yōu)點:
它能魯棒的估計模型參數(shù)。例如,它能從包含大量局外點的數(shù)據(jù)集中估計出高精度的參數(shù)。
缺點:
它計算參數(shù)的迭代次數(shù)沒有上限;如果設(shè)置迭代次數(shù)的上限,得到的結(jié)果可能不是最優(yōu)的結(jié)果,甚至可能得到錯誤的結(jié)果。
RANSAC只有一定的概率得到可信的模型,概率與迭代次數(shù)成正比。
它要求設(shè)置跟問題相關(guān)的閥值。
RANSAC只能從特定的數(shù)據(jù)集中估計出一個模型,如果存在兩個(或多個)模型,RANSAC不能找到別的模型。
要求數(shù)學(xué)模型已知
參考:https://zhuanlan.zhihu.com/p/36301702/
https://blog.csdn.net/leonardohaig/article/details/104570965
總結(jié)
以上是生活随笔為你收集整理的ransac算法思想的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 生活常识:如何正确使用洗衣皂?
- 下一篇: 仿制福特号还欠缺哪些技术?