全景图拼接——opencv之stitching

时间:2024-02-19 21:13:07
  1 """
  2 Stitching sample (advanced)
  3 ===========================
  4 
  5 Show how to use Stitcher API from python.
  6 """
  7 
  8 # Python 2/3 compatibility
  9 from __future__ import print_function
 10 
 11 import numpy as np
 12 import cv2 as cv
 13 
 14 import sys
 15 import argparse
 16 #创建parser实例
 17 parser = argparse.ArgumentParser(prog=\'stitching_detailed.py\', description=\'Rotation model images stitcher\')
 18 #继续增加参数
 19 parser.add_argument(\'img_names\', nargs=\'+\',help=\'files to stitch\',type=str)
 20 #两个必须的入参,执行文件名,图片名
 21 parser.add_argument(\'--preview\',help=\'Run stitching in the preview mode. Works faster than usual mode but output image will have lower resolution.\',type=bool,dest = \'preview\' )
 22 parser.add_argument(\'--try_cuda\',action = \'store\', default = False,help=\'Try to use CUDA. The default value is no. All default values are for CPU mode.\',type=bool,dest = \'try_cuda\' )
 23 #默认不用cuda
 24 #存储值到变量,action = store
 25 parser.add_argument(\'--work_megapix\',action = \'store\', default = 0.6,help=\' Resolution for image registration step. The default is 0.6 Mpx\',type=float,dest = \'work_megapix\' )
 26 parser.add_argument(\'--features\',action = \'store\', default = \'orb\',help=\'Type of features used for images matching. The default is orb.\',type=str,dest = \'features\' )
 27 parser.add_argument(\'--matcher\',action = \'store\', default = \'homography\',help=\'Matcher used for pairwise image matching.\',type=str,dest = \'matcher\' )
 28 parser.add_argument(\'--estimator\',action = \'store\', default = \'homography\',help=\'Type of estimator used for transformation estimation.\',type=str,dest = \'estimator\' )
 29 parser.add_argument(\'--match_conf\',action = \'store\', default = 0.3,help=\'Confidence for feature matching step. The default is 0.65 for surf and 0.3 for orb.\',type=float,dest = \'match_conf\' )
 30 parser.add_argument(\'--conf_thresh\',action = \'store\', default = 1.0,help=\'Threshold for two images are from the same panorama confidence.The default is 1.0.\',type=float,dest = \'conf_thresh\' )
 31 parser.add_argument(\'--ba\',action = \'store\', default = \'ray\',help=\'Bundle adjustment cost function. The default is ray.\',type=str,dest = \'ba\' )
 32 parser.add_argument(\'--ba_refine_mask\',action = \'store\', default = \'xxxxx\',help=\'Set refinement mask for bundle adjustment.  mask is "xxxxx"\',type=str,dest = \'ba_refine_mask\' )
 33 parser.add_argument(\'--wave_correct\',action = \'store\', default = \'horiz\',help=\'Perform wave effect correction. The default is "horiz"\',type=str,dest = \'wave_correct\' )
 34 parser.add_argument(\'--save_graph\',action = \'store\', default = None,help=\'Save matches graph represented in DOT language to <file_name> file.\',type=str,dest = \'save_graph\' )
 35 parser.add_argument(\'--warp\',action = \'store\', default = \'plane\',help=\'Warp surface type. The default is "spherical".\',type=str,dest = \'warp\' )
 36 parser.add_argument(\'--seam_megapix\',action = \'store\', default = 0.1,help=\' Resolution for seam estimation step. The default is 0.1 Mpx.\',type=float,dest = \'seam_megapix\' )
 37 #缝合线尺度0.1M大概10万像素,320*320
 38 parser.add_argument(\'--seam\',action = \'store\', default = \'no\',help=\'Seam estimation method. The default is "gc_color".\',type=str,dest = \'seam\' )
 39 parser.add_argument(\'--compose_megapix\',action = \'store\', default = -1,help=\'Resolution for compositing step. Use -1 for original resolution.\',type=float,dest = \'compose_megapix\' )
 40 parser.add_argument(\'--expos_comp\',action = \'store\', default = \'no\',help=\'Exposure compensation method. The default is "gain_blocks".\',type=str,dest = \'expos_comp\' )
 41 parser.add_argument(\'--expos_comp_nr_feeds\',action = \'store\', default = 1,help=\'Number of exposure compensation feed.\',type=np.int32,dest = \'expos_comp_nr_feeds\' )
 42 parser.add_argument(\'--expos_comp_nr_filtering\',action = \'store\', default = 2,help=\'Number of filtering iterations of the exposure compensation gains\',type=float,dest = \'expos_comp_nr_filtering\' )
 43 parser.add_argument(\'--expos_comp_block_size\',action = \'store\', default = 32,help=\'BLock size in pixels used by the exposure compensator.\',type=np.int32,dest = \'expos_comp_block_size\' )
 44 parser.add_argument(\'--blend\',action = \'store\', default = \'multiband\',help=\'Blending method. The default is "multiband".\',type=str,dest = \'blend\' )
 45 parser.add_argument(\'--blend_strength\',action = \'store\', default = 5,help=\'Blending strength from [0,100] range.\',type=np.int32,dest = \'blend_strength\' )
 46 parser.add_argument(\'--output\',action = \'store\', default = \'result.jpg\',help=\'The default is "result.jpg"\',type=str,dest = \'output\' )
 47 parser.add_argument(\'--timelapse\',action = \'store\', default = None,help=\'Output warped images separately as frames of a time lapse movie, with "fixed_" prepended to input file names.\',type=str,dest = \'timelapse\' )
 48 parser.add_argument(\'--rangewidth\',action = \'store\', default = -1,help=\'uses range_width to limit number of images to match with.\',type=int,dest = \'rangewidth\' )
 49 
 50 __doc__ += \'\n\' + parser.format_help()
 51 # __doc__可以用函数名调用,会输出函数下的"""xxxxxx """三引号注释中的内容
 52 
 53 def main():
 54     args = parser.parse_args()
 55     #参数解析
 56     #参数应用
 57     img_names=args.img_names  
 58 #这种传参方法才是合理的,因为服务器传来的params是一个整体的json文件。因此参数解析要用类字典的形式。来的是整体,即便只有一个参数也不能直接类似params = 1如此传参
 59     print(img_names)
 60     preview = args.preview
 61     try_cuda = args.try_cuda
 62     work_megapix = args.work_megapix
 63     seam_megapix = args.seam_megapix
 64     compose_megapix = args.compose_megapix
 65     conf_thresh = args.conf_thresh
 66     features_type = args.features
 67     matcher_type = args.matcher
 68     estimator_type = args.estimator
 69     ba_cost_func = args.ba
 70     ba_refine_mask = args.ba_refine_mask
 71     wave_correct = args.wave_correct
 72     if wave_correct==\'no\':
 73         do_wave_correct= False
 74     else:
 75         do_wave_correct=True
 76     if args.save_graph is None:
 77         save_graph = False
 78     else:
 79         save_graph =True
 80         save_graph_to = args.save_graph
 81     warp_type = args.warp
 82     #曝光补偿供4种方式:增量补偿,块增量补偿,通道补偿,块通道补偿
 83     if args.expos_comp==\'no\':
 84         expos_comp_type = cv.detail.ExposureCompensator_NO
 85     elif  args.expos_comp==\'gain\':
 86         expos_comp_type = cv.detail.ExposureCompensator_GAIN
 87     elif  args.expos_comp==\'gain_blocks\':
 88         expos_comp_type = cv.detail.ExposureCompensator_GAIN_BLOCKS
 89     elif  args.expos_comp==\'channel\':
 90         expos_comp_type = cv.detail.ExposureCompensator_CHANNELS
 91     elif  args.expos_comp==\'channel_blocks\':
 92         expos_comp_type = cv.detail.ExposureCompensator_CHANNELS_BLOCKS
 93     else:
 94         print("Bad exposure compensation method")
 95         exit()
 96     expos_comp_nr_feeds = args.expos_comp_nr_feeds
 97     expos_comp_nr_filtering = args.expos_comp_nr_filtering
 98     expos_comp_block_size = args.expos_comp_block_size
 99     match_conf = args.match_conf
100     seam_find_type = args.seam
101     blend_type = args.blend
102     blend_strength = args.blend_strength
103     result_name = args.output
104     if args.timelapse is not None:
105         timelapse = True
106         if args.timelapse=="as_is":
107             timelapse_type = cv.detail.Timelapser_AS_IS
108         elif args.timelapse=="crop":
109             timelapse_type = cv.detail.Timelapser_CROP
110         else:
111             print("Bad timelapse method")
112             exit()
113     else:
114         timelapse= False
115     range_width = args.rangewidth
116     #特征描述方法,三种方式 orb                surf      sift, 参数解析默认值是orb
117     if features_type==\'orb\':
118         finder= cv.ORB.create()
119     elif features_type==\'surf\':
120         finder= cv.xfeatures2d_SURF.create()
121     elif features_type==\'sift\':
122         finder= cv.xfeatures2d_SIFT.create()
123     else:
124         print ("Unknown descriptor type")
125         exit()
126     seam_work_aspect = 1
127     full_img_sizes=[]
128     features=[]
129     images=[]
130     is_work_scale_set = False
131     is_seam_scale_set = False
132     is_compose_scale_set = False;
133     for name in img_names:#循环 读图放到list中
134         full_img = cv.imread(cv.samples.findFile(name))
135         
136         if full_img is None:
137             print("Cannot read image ", name)
138             exit()
139         full_img_sizes.append((full_img.shape[1],full_img.shape[0]))
140         if work_megapix < 0:
141             img = full_img
142             work_scale = 1
143             is_work_scale_set = True
144         else:
145             if is_work_scale_set is False:
146                 work_scale = min(1.0, np.sqrt(work_megapix * 1e6 / (full_img.shape[0]*full_img.shape[1])))
147                 is_work_scale_set = True
148             img = cv.resize(src=full_img, dsize=None, fx=work_scale, fy=work_scale, interpolation=cv.INTER_LINEAR_EXACT)
149         if is_seam_scale_set is False:
150             seam_scale = min(1.0, np.sqrt(seam_megapix * 1e6 / (full_img.shape[0]*full_img.shape[1])))
151             seam_work_aspect = seam_scale / work_scale
152             #缝合区域大小/配准
153             is_seam_scale_set = True
154         imgFea= cv.detail.computeImageFeatures2(finder,img)
155         #单张图片的特征提取,finder是特征提取器对象
156         features.append(imgFea)
157         #根据以上来调整图片大小,保存到images图片列表中
158         img = cv.resize(src=full_img, dsize=None, fx=seam_scale, fy=seam_scale, interpolation=cv.INTER_LINEAR_EXACT)
159         images.append(img)
160     if matcher_type== "affine":
161     #匹配类型  仿射
162         matcher = cv.detail_AffineBestOf2NearestMatcher(False, try_cuda, match_conf)
163     elif range_width==-1:
164         matcher = cv.detail.BestOf2NearestMatcher_create(try_cuda, match_conf)
165     else:
166         matcher = cv.detail.BestOf2NearestRangeMatcher_create(range_width, try_cuda, match_conf)
167     #此时的features是所有入参图片的特征集合    
168     p=matcher.apply2(features)
169     matcher.collectGarbage()
170     if save_graph:
171         f = open(save_graph_to,"w")
172         f.write(cv.detail.matchesGraphAsString(img_names, p, conf_thresh))
173         f.close()
174     indices=cv.detail.leaveBiggestComponent(features,p,0.3)
175     img_subset =[]
176     img_names_subset=[]
177     full_img_sizes_subset=[]
178     num_images=len(indices)
179     for i in range(len(indices)):
180         img_names_subset.append(img_names[indices[i,0]])
181         img_subset.append(images[indices[i,0]])
182         full_img_sizes_subset.append(full_img_sizes[indices[i,0]])
183     images = img_subset;
184     img_names = img_names_subset;
185     full_img_sizes = full_img_sizes_subset;
186     num_images = len(img_names)
187     if num_images < 2:
188         print("Need more images")
189         exit()
190 
191     if estimator_type == "affine":
192         estimator = cv.detail_AffineBasedEstimator()#基于仿射变换特征
193     else:
194         estimator = cv.detail_HomographyBasedEstimator()#基于透视变换的单应阵
195     b, cameras =estimator.apply(features,p,None)#变换估计出变换矩阵和相机参数
196     if not b:
197         print("Homography estimation failed.")
198         exit()
199     for cam in cameras:
200     #相机的旋转矩阵
201         cam.R=cam.R.astype(np.float32)
202 
203     if ba_cost_func == "reproj":
204         adjuster = cv.detail_BundleAdjusterReproj()
205     elif ba_cost_func == "ray":#默认,光束平差的损失函数类型
206         adjuster = cv.detail_BundleAdjusterRay()
207     elif ba_cost_func == "affine":
208         adjuster = cv.detail_BundleAdjusterAffinePartial()
209     elif ba_cost_func == "no":
210         adjuster = cv.detail_NoBundleAdjuster()
211     else:
212         print( "Unknown bundle adjustment cost function: ", ba_cost_func )
213         exit()
214     adjuster.setConfThresh(1)#确定两幅图片属于同一张全景图的阈值
215     refine_mask=np.zeros((3,3),np.uint8)
216     #掩膜,3*3的0阵,如果要优化则将对应位置赋成1
217 """数字图像处理中,掩模为二维矩阵数组,有时也用多值图像。数字图像处理中,图像掩模主要用于:
218 
219 ①提取感兴趣区,用预先制作的感兴趣区掩模与待处理图像相乘,得到感兴趣区图像,感兴趣区内图像值保持不变,而区外图像值都为0。 
220 ②屏蔽作用,用掩模对图像上某些区域作屏蔽,使其不参加处理或不参加处理参数的计算,或仅对屏蔽区作处理或统计。 
221 ③结构特征提取,用相似性变量或图像匹配方法检测和提取图像中与掩模相似的结构特征。 
222 """
223 
224     if ba_refine_mask[0] == \'x\':
225         refine_mask[0,0] = 1
226     if ba_refine_mask[1] == \'x\':
227         refine_mask[0,1] = 1
228     if ba_refine_mask[2] == \'x\':
229         refine_mask[0,2] = 1
230     if ba_refine_mask[3] == \'x\':
231         refine_mask[1,1] = 1
232     if ba_refine_mask[4] == \'x\':
233         refine_mask[1,2] = 1
234     adjuster.setRefinementMask(refine_mask)
235     b,cameras = adjuster.apply(features,p,cameras)
236     #光束平差,精修相机姿态等
237 #光束是指相机光心发出的射线。损失函数一般有重映射误差和射线距离误差两种
238     if not b:
239         print("Camera parameters adjusting failed.")
240         exit()
241     focals=[]
242     for cam in cameras:
243         focals.append(cam.focal) #焦距
244     sorted(focals)
245     #可以通过两张图的匹配点或单应阵求出焦距,用所有焦距的平均作为全景图焦距的近似。用这个焦距作为柱面的半径,先将图像投射到柱面,在柱面上区拼接
246     if len(focals)%2==1:
247     #根据焦距确定 图片投影的尺度 warped_image_scale
248         warped_image_scale = focals[len(focals) // 2]
249     else:
250         warped_image_scale = (focals[len(focals) // 2]+focals[len(focals) // 2-1])/2
251     if do_wave_correct:
252     #波形修正得到修正后的旋转矩阵
253         rmats=[]
254         for cam in cameras:
255             rmats.append(np.copy(cam.R))
256         rmats    =    cv.detail.waveCorrect(    rmats,  cv.detail.WAVE_CORRECT_HORIZ)
257         for idx,cam in enumerate(cameras):
258             cam.R = rmats[idx]
259     corners=[]
260     mask=[]
261     masks_warped=[]
262     images_warped=[]
263     sizes=[]
264     masks=[]
265     for i in range(0,num_images):
266         um=cv.UMat(255*np.ones((images[i].shape[0],images[i].shape[1]),np.uint8))
267         masks.append(um)
268     #处理相机旋转造成的扭曲
269     warper = cv.PyRotationWarper(warp_type,warped_image_scale*seam_work_aspect) # warper peut etre nullptr?
270     for idx in range(0,num_images):
271         K = cameras[idx].K().astype(np.float32)#相机内参
272         swa = seam_work_aspect
273         K[0,0] *= swa
274         K[0,2] *= swa
275         K[1,1] *= swa
276         K[1,2] *= swa
277         corner,image_wp =warper.warp(images[idx],K,cameras[idx].R,cv.INTER_LINEAR, cv.BORDER_REFLECT)
278         corners.append(corner)
279         sizes.append((image_wp.shape[1],image_wp.shape[0]))
280         images_warped.append(image_wp)
281 
282         p,mask_wp =warper.warp(masks[idx],K,cameras[idx].R,cv.INTER_NEAREST, cv.BORDER_CONSTANT)#线性差值,连续边界处理
283         masks_warped.append(mask_wp.get())
284     images_warped_f=[]
285     for img in images_warped:#对投影后的图片,单精度存储
286         imgf=img.astype(np.float32)
287         images_warped_f.append(imgf)
288     if cv.detail.ExposureCompensator_CHANNELS == expos_comp_type:
289         compensator = cv.detail_ChannelsCompensator(expos_comp_nr_feeds)
290     #    compensator.setNrGainsFilteringIterations(expos_comp_nr_filtering)
291     elif cv.detail.ExposureCompensator_CHANNELS_BLOCKS == expos_comp_type:
292         compensator=cv.detail_BlocksChannelsCompensator(expos_comp_block_size, expos_comp_block_size,expos_comp_nr_feeds)
293     #    compensator.setNrGainsFilteringIterations(expos_comp_nr_filtering)
294     else:
295         compensator=cv.detail.ExposureCompensator_createDefault(expos_comp_type)
296         #曝光补偿需要        分块 投影图片 和 掩膜
297 #对于分块补偿,得到分块曝光补偿系数后,此时直接补偿会出现明显的视觉分"块"效果还要对系数进行分段线性滤波处理,得到全局最终的补偿系数
298     compensator.feed(corners=corners, images=images_warped, masks=masks_warped)
299     #缝合线(两图最佳的一条融合曲线)查找描述子
300 #一般有三种方法:1基于距离的逐点法  2.动态规划查找法(快,省内存)   3.最大流图割法
301     if seam_find_type == "no":
302         seam_finder = cv.detail.SeamFinder_createDefault(cv.detail.SeamFinder_NO)
303     elif seam_find_type == "voronoi":
304         seam_finder = cv.detail.SeamFinder_createDefault(cv.detail.SeamFinder_VORONOI_SEAM);
305     elif seam_find_type == "gc_color":
306         seam_finder = cv.detail_GraphCutSeamFinder("COST_COLOR")
307     elif seam_find_type == "gc_colorgrad":
308         seam_finder = cv.detail_GraphCutSeamFinder("COST_COLOR_GRAD")
309     elif seam_find_type == "dp_color":
310         seam_finder = cv.detail_DpSeamFinder("COLOR")
311     elif seam_find_type == "dp_colorgrad":
312         seam_finder = cv.detail_DpSeamFinder("COLOR_GRAD")
313     if seam_finder is None:
314         print("Can\'t create the following seam finder ",seam_find_type)
315         exit()
316     seam_finder.find(images_warped_f, corners,masks_warped )
317     imgListe=[]
318     compose_scale=1
319     corners=[]
320     sizes=[]
321     images_warped=[]
322     images_warped_f=[]
323     masks=[]
324     blender= None
325     timelapser=None
326     compose_work_aspect=1
327     for idx,name in enumerate(img_names): # https://github.com/opencv/opencv/blob/master/samples/cpp/stitching_detailed.cpp#L725 ?
328         full_img  = cv.imread(name)
329         if not is_compose_scale_set:
330             if compose_megapix > 0:
331             #exp是2.718; e是10^
332                 compose_scale = min(1.0, np.sqrt(compose_megapix * 1e6 / (full_img.shape[0]*full_img.shape[1])))
333             is_compose_scale_set = True;
334             compose_work_aspect = compose_scale / work_scale;
335             warped_image_scale *= compose_work_aspect
336             warper =  cv.PyRotationWarper(warp_type,warped_image_scale)
337             for i in range(0,len(img_names)):
338                 cameras[i].focal *= compose_work_aspect
339                 cameras[i].ppx *= compose_work_aspect
340                 cameras[i].ppy *= compose_work_aspect
341                 sz = (full_img_sizes[i][0] * compose_scale,full_img_sizes[i][1]* compose_scale)
342                 K = cameras[i].K().astype(np.float32)
343                 roi = warper.warpRoi(sz, K, cameras[i].R);
344                 #缝合区获得
345                 corners.append(roi[0:2])
346                 sizes.append(roi[2:4])
347         if abs(compose_scale - 1) > 1e-1:
348             img =cv.resize(src=full_img, dsize=None, fx=compose_scale, fy=compose_scale, interpolation=cv.INTER_LINEAR_EXACT)
349         else:
350             img = full_img;
351         img_size = (img.shape[1],img.shape[0]);
352         K=cameras[idx].K().astype(np.float32)
353         corner,image_warped =warper.warp(img,K,cameras[idx].R,cv.INTER_LINEAR, cv.BORDER_REFLECT)
354         mask =255*np.ones((img.shape[0],img.shape[1]),np.uint8)
355         p,mask_warped =warper.warp(mask,K,cameras[idx].R,cv.INTER_NEAREST, cv.BORDER_CONSTANT)
356         #块增益补偿
357         compensator.apply(idx,corners[idx],image_warped,mask_warped)
358         image_warped_s = image_warped.astype(np.int16)
359         image_warped=[]
360         #膨胀
361         dilated_mask = cv.dilate(masks_warped[idx],None)
362         seam_mask = cv.resize(dilated_mask,(mask_warped.shape[1],mask_warped.shape[0]),0,0,cv.INTER_LINEAR_EXACT)
363         mask_warped = cv.bitwise_and(seam_mask,mask_warped)#二进制的与运算,1&1 == 1    ,1&0 == 0
364         if blender==None and not timelapse:
365             blender = cv.detail.Blender_createDefault(cv.detail.Blender_NO)
366             dst_sz = cv.detail.resultRoi(corners=corners,sizes=sizes)
367             blend_width = np.sqrt(dst_sz[2]*dst_sz[3]) * blend_strength / 100
368             if blend_width < 1:
369                 blender = cv.detail.Blender_createDefault(cv.detail.Blender_NO)
370             elif blend_type == "multiband":
371                 blender = cv.detail_MultiBandBlender()
372                 blender.setNumBands((np.log(blend_width)/np.log(2.) - 1.).astype(np.int))
373             elif blend_type == "feather":#羽化融合
374                 blender = cv.detail_FeatherBlender()
375                 blender.setSharpness(1./blend_width)
376             blender.prepare(dst_sz)
377         elif timelapser==None  and timelapse:
378             timelapser = cv.detail.Timelapser_createDefault(timelapse_type)
379             timelapser.initialize(corners, sizes)
380         if timelapse:
381         #延时处理
382             matones=np.ones((image_warped_s.shape[0],image_warped_s.shape[1]), np.uint8)
383             timelapser.process(image_warped_s, matones, corners[idx])
384             pos_s = img_names[idx].rfind("/");
385             if pos_s == -1:
386                 fixedFileName = "fixed_" + img_names[idx];
387             else:
388                 fixedFileName = img_names[idx][:pos_s + 1 ]+"fixed_" + img_names[idx][pos_s + 1: ]
389             cv.imwrite(fixedFileName, timelapser.getDst())
390         else:
391             blender.feed(cv.UMat(image_warped_s), mask_warped, corners[idx])
392     if not timelapse:
393         result=None
394         result_mask=None
395         result,result_mask = blender.blend(result,result_mask)
396         cv.imwrite(result_name,result)
397         zoomx = 600.0 / result.shape[1]
398         dst=cv.normalize(src=result,dst=None,alpha=255.,norm_type=cv.NORM_MINMAX,dtype=cv.CV_8U)
399         dst=cv.resize(dst,dsize=None,fx=zoomx,fy=zoomx)
400         cv.imshow(result_name,dst)
401         cv.waitKey()
402 
403     print(\'Done\')
404 
405 
406 if __name__ == \'__main__\':
407     print(__doc__)
408     main()
409     cv.destroyAllWindows()