在使用MATLAB进行图像滤波时如何处理边缘问题
我们在写滤波程序时一般会用矩阵模板与原图像做卷积,这时候在做图像边界的处理是一般都选择忽略边缘,不过要是模板比较大,那么处理的效果就不好了,图像四周就会是原图像,中间才是滤波后的结果,虽然用Matlab的imfilter就能解决,不过还是自己通过滤波的原理实践一下比较好。
模板和图像一共有如下16种关系,我粗略的画了一下,前三张小矩形的是模板、大的矩形是图像,最后一张大的是模板,小的是图像。这就是图像和模板卷积时的所有关系。
看似好像要写16个if判断,其实是不用的,我们只要判断卷积时模板的四个边界和图像的四个边界的关系就行了。这里有两对相对坐标,一个是表示图像的卷积范围,一个表示模板的卷积范围。先说怎么表示图像的卷积范围吧,如果当前处理的点是(i,j),模板大小都是2*r+1(我这里都用了对称的奇数模板,模板边界像素要是偶数会很难处理,这里我干脆把奇数也化成了偶数,原理差不多的)。八个边界可以这样表示,图像就是1表示图像上边缘,m表示图像下边缘,1表示图像左边缘,n表示图像右边缘,i-r表示模板上边缘,i+r表示模板下边缘,j-r表示模板左边缘,j+r表示模板右边缘。通过这四对的组合就能16个关系,具体还是见下面的代码吧,看注释结合代码比较清楚。
main.m
1 clear all;
2 close all;
3 clc;
4 r=20;
5 w=fspecial('average',[2*r+1 2*r+1]);
6
7 img=imread('lena.jpg');
8 img=mat2gray(img);
9 [m n]=size(img);
10 imshow(img);
11
12
13 imgn=filterim(img,w);
14 figure;
15 imshow(mat2gray(imgn));
16
17 imgn=img;
18
19 for i=r+1:m-r
20 for j=r+1:n-r
21
22 imgn(i,j)=sum(sum(img(i-r:i+r,j-r:j+r).*w));
23
24 end
25 end
26 figure;
27 imshow(mat2gray(imgn));
28
29
30 figure;
31 img=imfilter(img,w);
32 imshow(mat2gray(img))
filterim.m(实现主要功能):
1 function imgn=filterim(img,w)
2
3 [r r]=size(w);
4 [m n]=size(img);
5 if mod(r,2)==0
6 r=r+1;
7 w=imresize(w,[r r]);
8 end
9 imgn=zeros(m,n);
10 r=floor(r/2);
11
12 for i=1:m
13 for j=1:n
14 %图像需要获得四个边界的卷积范围,模板只需要获得最上面和最左面就可以了,因为图像和模板两个卷积范围是一样的。
15 if i-r<1 %判断模板上边缘和图像上边缘的关系
16 img_up=1; %如果当前像素的高小于模板的一半,那么选择图像的上边缘作为卷积图像的上边缘
17 mark_up=r-i+1; %模板的上边缘使用和图像相交的上边缘
18 else
19 img_up=i-r; %使用当前像素的高减去模板的一半作为卷积图像的上边缘
20 mark_up=1; %使用模板的最上边缘作为卷积模板的上边缘
21 end
22
23 if i+r>m %判断模板下边缘和图像下边缘的关系
24 img_down=m; %如果当前像素的高加上模板的一半超过整个图像的高,那么卷积图像的下边缘使用整个图像的下边缘
25 else
26 img_down=i+r; %否则卷积图像的下边缘使用当前像素的高加上模板的一半
27 end
28
29 if j-r<1 %判断模板左边缘和图像左边缘的关系
30 img_left=1; %如果当前像素的宽小于模板的一半,那么选择图像的左边缘作为卷积图像的左边缘
31 mark_left=r-j+1; %模板的左边缘使用和图像相交的左边缘
32 else
33 img_left=j-r; %使用当前像素的宽减去模板的一半作为卷积图像的左边缘
34 mark_left=1; %使用模板的最左边缘作为卷积模板的左边缘
35 end
36
37 if j+r>n %判断模板右边缘和图像右边缘的关系
38 img_right=n; %如果当前像素的宽加上模板的一般超过整个图像的宽,那么卷积图像的右边缘使用整个图像的右边缘
39 else
40 img_right=j+r; %否则卷积图像的右边缘使用当前像素的宽加上模板的一半
41 end
42
43 imgn(i,j)=sum(sum(img(img_up:img_down,img_left:img_right).*w(mark_up:mark_up+(img_down-img_up),mark_left:mark_left+(img_right-img_left))));%/((img_down-img_up+1)*(img_right-img_left+1));
44 %卷积图像上边缘:下边缘,左边缘:右边缘 %卷积模板上边缘:上边缘+(竖直卷积范围),卷积模板左边缘:左边缘+(水平卷积范围)
45 end
46 end
47 end
说实话,这个注释写的我很纠结,注释我已经尽力写清楚了,虽然我感觉还是没解释清楚。我果然需要锻炼一下写作与表达能力。
下面是效果图:
原图
自己写的滤波时边界处理的效果
通常的不做边界处理的效果
Matlab函数处理的结果
我这个效果基本和matlab自带的效果很接近了,不过速度好像慢很多,matlab自带的函数也许使用汇编处理的,总之,算法就是这样啦。
推荐阅读
-
epoll简介及触发模式(accept、read、send)-epoll的简单介绍 epoll在LT和ET模式下的读写方式 一、epoll的接口非常简单,一共就三个函数:1. int epoll_create(int size);创建一个epoll的句柄,size用来告诉内核这个监听的数目一共有多大。这个参数不同于select中的第一个参数,给出最大监听的fd+1的值。需要注意的是,当创建好epoll句柄后,它就是会占用一个fd值,在linux下如果查看/proc/进程id/fd/,是能够看到这个fd的,所以在使用完epoll后,必须调用close关闭,否则可能导致fd被耗尽。2. int epoll_ctl(int epfd, int op, int fd, struct epoll_event *event);epoll的事件注册函数,它不同与select是在监听事件时告诉内核要监听什么类型的事件,而是在这里先注册要监听的事件类型。第一个参数是epoll_create的返回值,第二个参数表示动作,用三个宏来表示:EPOLL_CTL_ADD:注册新的fd到epfd中;EPOLL_CTL_MOD:修改已经注册的fd的监听事件;EPOLL_CTL_DEL:从epfd中删除一个fd;第三个参数是需要监听的fd,第四个参数是告诉内核需要监听什么事,struct epoll_event结构如下:struct epoll_event { __uint32_t events; /* Epoll events */ epoll_data_t data; /* User data variable */};events可以是以下几个宏的集合:EPOLLIN :表示对应的文件描述符可以读(包括对端SOCKET正常关闭); EPOLLIN事件:EPOLLIN事件则只有当对端有数据写入时才会触发,所以触发一次后需要不断读取所有数据直到读完EAGAIN为止。否则剩下的数据只有在下次对端有写入时才能一起取出来了。现在明白为什么说epoll必须要求异步socket了吧?如果同步socket,而且要求读完所有数据,那么最终就会在堵死在阻塞里。 EPOLLOUT:表示对应的文件描述符可以写; EPOLLOUT事件:EPOLLOUT事件只有在连接时触发一次,表示可写,其他时候想要触发,那要先准备好下面条件:1.某次write,写满了发送缓冲区,返回错误码为EAGAIN。2.对端读取了一些数据,又重新可写了,此时会触发EPOLLOUT。简单地说:EPOLLOUT事件只有在不可写到可写的转变时刻,才会触发一次,所以叫边缘触发,这叫法没错的!其实,如果真的想强制触发一次,也是有办法的,直接调用epoll_ctl重新设置一下event就可以了,event跟原来的设置一模一样都行(但必须包含EPOLLOUT),关键是重新设置,就会马上触发一次EPOLLOUT事件。1. 缓冲区由满变空.2.同时注册EPOLLIN | EPOLLOUT事件,也会触发一次EPOLLOUT事件这个两个也会触发EPOLLOUT事件 EPOLLPRI:表示对应的文件描述符有紧急的数据可读(这里应该表示有带外数据到来);EPOLLERR:表示对应的文件描述符发生错误;EPOLLHUP:表示对应的文件描述符被挂断;EPOLLET: 将EPOLL设为边缘触发(Edge Triggered)模式,这是相对于水平触发(Level Triggered)来说的。EPOLLONESHOT:只监听一次事件,当监听完这次事件之后,如果还需要继续监听这个socket的话,需要再次把这个socket加入到EPOLL队列里3. int epoll_wait(int epfd, struct epoll_event * events, int maxevents, int timeout);等待事件的产生,类似于select调用。参数events用来从内核得到事件的集合,maxevents告之内核这个events有多大,这个maxevents的值不能大于创建epoll_create时的size,参数timeout是超时时间(毫秒,0会立即返回,-1将不确定,也有说法说是永久阻塞)。该函数返回需要处理的事件数目,如返回0表示已超时。-------------------------------------------------------------------------------------------- 从man手册中,得到ET和LT的具体描述如下EPOLL事件有两种模型:Edge Triggered (ET)Level Triggered (LT)假如有这样一个例子:1. 我们已经把一个用来从管道中读取数据的文件句柄(RFD)添加到epoll描述符2. 这个时候从管道的另一端被写入了2KB的数据3. 调用epoll_wait(2),并且它会返回RFD,说明它已经准备好读取操作4. 然后我们读取了1KB的数据5. 调用epoll_wait(2)......Edge Triggered 工作模式:如果我们在第1步将RFD添加到epoll描述符的时候使用了EPOLLET标志,那么在第5步调用epoll_wait(2)之后将有可能会挂起,因为剩余的数据还存在于文件的输入缓冲区内,而且数据发出端还在等待一个针对已经发出数据的反馈信息。只有在监视的文件句柄上发生了某个事件的时候 ET 工作模式才会汇报事件。因此在第5步的时候,调用者可能会放弃等待仍在存在于文件输入缓冲区内的剩余数据。在上面的例子中,会有一个事件产生在RFD句柄上,因为在第2步执行了一个写操作,然后,事件将会在第3步被销毁。因为第4步的读取操作没有读空文件输入缓冲区内的数据,因此我们在第5步调用 epoll_wait(2)完成后,是否挂起是不确定的。epoll工作在ET模式的时候,必须使用非阻塞套接口,以避免由于一个文件句柄的阻塞读/阻塞写操作把处理多个文件描述符的任务饿死。最好以下面的方式调用ET模式的epoll接口,在后面会介绍避免可能的缺陷。 i 基于非阻塞文件句柄 ii 只有当read(2)或者write(2)返回EAGAIN时才需要挂起,等待。但这并不是说每次read时都需要循环读,直到读到产生一个EAGAIN才认为此次事件处理完成,当read返回的读到的数据长度小于请求的数据长度时,就可以确定此时缓冲中已没有数据了,也就可以认为此事读事件已处理完成。Level Triggered 工作模式相反的,以LT方式调用epoll接口的时候,它就相当于一个速度比较快的poll(2),并且无论后面的数据是否被使用,因此他们具有同样的职能。因为即使使用ET模式的epoll,在收到多个chunk的数据的时候仍然会产生多个事件。调用者可以设定EPOLLONESHOT标志,在 epoll_wait(2)收到事件后epoll会与事件关联的文件句柄从epoll描述符中禁止掉。因此当EPOLLONESHOT设定后,使用带有 EPOLL_CTL_MOD标志的epoll_ctl(2)处理文件句柄就成为调用者必须作的事情。然后详细解释ET, LT:LT(level triggered)是缺省的工作方式,并且同时支持block和no-block socket.在这种做法中,内核告诉你一个文件描述符是否就绪了,然后你可以对这个就绪的fd进行IO操作。如果你不作任何操作,内核还是会继续通知你的,所以,这种模式编程出错误可能性要小一点。传统的select/poll都是这种模型的代表.ET(edge-triggered)是高速工作方式,只支持no-block socket。在这种模式下,当描述符从未就绪变为就绪时,内核通过epoll告诉你。然后它会假设你知道文件描述符已经就绪,并且不会再为那个文件描述符发送更多的就绪通知,直到你做了某些操作导致那个文件描述符不再为就绪状态了(比如,你在发送,接收或者接收请求,或者发送接收的数据少于一定量时导致了一个EWOULDBLOCK 错误)。但是请注意,如果一直不对这个fd作IO操作(从而导致它再次变成未就绪),内核不会发送更多的通知(only once),不过在TCP协议中,ET模式的加速效用仍需要更多的benchmark确认(这句话不理解)。在许多测试中我们会看到如果没有大量的idle -connection或者dead-connection,epoll的效率并不会比select/poll高很多,但是当我们遇到大量的idle- connection(例如WAN环境中存在大量的慢速连接),就会发现epoll的效率大大高于select/poll。(未测试)另外,当使用epoll的ET模型来工作时,当产生了一个EPOLLIN事件后,读数据的时候需要考虑的是当recv返回的大小如果等于请求的大小,那么很有可能是缓冲区还有数据未读完,也意味着该次事件还没有处理完,所以还需要再次读取: 这里只是说明思路(参考《UNIX网络编程》) while(rs) {buflen = recv(activeevents[i].data.fd, buf, sizeof(buf), 0);if(buflen < 0){// 由于是非阻塞的模式,所以当errno为EAGAIN时,表示当前缓冲区已无数据可读// 在这里就当作是该次事件已处理处.if(errno == EAGAIN)break; else return; }else if(buflen == 0) { // 这里表示对端的socket已正常关闭. } if(buflen == sizeof(buf) rs = 1; // 需要再次读取 else rs = 0; } 还有,假如发送端流量大于接收端的流量(意思是epoll所在的程序读比转发的socket要快),由于是非阻塞的socket,那么send函数虽然返回,但实际缓冲区的数据并未真正发给接收端,这样不断的读和发,当缓冲区满后会产生EAGAIN错误(参考man send),同时,不理会这次请求发送的数据.所以,需要封装socket_send的函数用来处理这种情况,该函数会尽量将数据写完再返回,返回-1表示出错。在socket_send内部,当写缓冲已满(send返回-1,且errno为EAGAIN),那么会等待后再重试.这种方式并不很完美,在理论上可能会长时间的阻塞在socket_send内部,但暂没有更好的办法. ssize_t socket_send(int sockfd, const char* buffer, size_t buflen) { ssize_t tmp; size_t total = buflen; const char *p = buffer; while(1) { tmp = send(sockfd, p, total, 0); if(tmp < 0) { // 当send收到信号时,可以继续写,但这里返回-1. if(errno == EINTR) return -1; // 当socket是非阻塞时,如返回此错误,表示写缓冲队列已满, // 在这里做延时后再重试. if(errno == EAGAIN) { usleep(1000); continue; } return -1; } if((size_t)tmp == total) return buflen; total -= tmp; p += tmp; } return tmp; } 二、epoll在LT和ET模式下的读写方式 在一个非阻塞的socket上调用read/write函数, 返回EAGAIN或者EWOULDBLOCK(注: EAGAIN就是EWOULDBLOCK) 从字面上看, 意思是: * EAGAIN: 再试一次 * EWOULDBLOCK: 如果这是一个阻塞socket, 操作将被block * perror输出: Resource temporarily unavailable 总结: 这个错误表示资源暂时不够, 可能read时, 读缓冲区没有数据, 或者, write时,写缓冲区满了 。 遇到这种情况, 如果是阻塞socket, read/write就要阻塞掉。 而如果是非阻塞socket, read/write立即返回-1, 同 时errno设置为EAGAIN. 所以, 对于阻塞socket, read/write返回-1代表网络出错了. 但对于非阻塞socket, read/write返回-1不一定网络真的出错了. 可能是Resource temporarily unavailable. 这时你应该再试, 直到Resource available. 综上, 对于non-blocking的socket, 正确的读写操作为: 读: 忽略掉errno = EAGAIN的错误, 下次继续读 写: 忽略掉errno = EAGAIN的错误, 下次继续写 对于select和epoll的LT模式, 这种读写方式是没有问题的. 但对于epoll的ET模式, 这种方式还有漏洞. epoll的两种模式 LT 和 ET
-
卷积的意义--我见过最生动易懂的解释--就是在图像处理中,将两组分辨率不同的图像进行卷积处理,从而形成易于处理的平滑图像。卷积甚至可以用在考试作弊中,为了让照片中的两个人同时像,只要对两个人的图像进行卷积处理就可以了,这是一种平滑处理,但我们如何才能真正把这个公式与实际建立一种联系,也就是说我们能不能从生活中找到一个很方便具体的例子来表达这个公式的物理意义呢? 有一个七品县令,喜欢打骂无赖,并有一个惯例:只要不犯大罪,只打一顿就放他回家,以示爱民如子。 有一种无赖,想扬名立万却又不抱多大希望,心想:既然扬不了好名,出了臭名也成啊。怎样才能出恶名呢?炒作!怎么炒作?找名人!他自然而然地想到了自己的长官--县令。 无赖于是在光天化日之下,站在县衙门口撒了泡尿,后果可想而知,自然是被请进堂上挨了板子,然后昂首挺胸地回家,躺了一天,哎!身体并无大碍!第二天照样如此,全然不顾行政长管的仁慈和衙门的尊严,第三天、第四天 ......每天去县衙领板子回来,还兴高采烈,坚持了一个月之久!这个无赖的名声像衙门口的臭气一样传遍了八方! 县太爷噤了噤鼻子,愣愣地望着惊堂木案,皱了皱眉头,思考着一个问题:这三十块大木板怎么会不好用呢?......想想也是,当年这位大人金榜题名的时候,我数学考了满分,所以这道题至少今天得解出来: --人(系统!)会怎么样(系统!)之后会怎么样(输出!)人(系统!)被打之后会怎么样? --有什么用,很疼! --我问的是:会发生什么? --取决于有多疼。就像这个无赖的体质,每天挨一板什么事都不会发生,连哼哼两声都不行,你看他那得意洋洋的样子(输出 0);如果一次连打他十板,他可能会皱着眉头,咬着牙,硬是不哼一声(输出 1);打到二十板,他会疼得脸都变形了,像猪一样哼哼唧唧(输出 3);打到三十板,他可能会像驴一样嚎叫,一把鼻涕一把泪,求你饶他一命(输出 5);打到四十板,他会大小便失禁,勉强哼哼(输出 1);打到五十板,他连哼哼都不能哼一下(输出 0)--死! 县官摊开坐标纸,绘制了一条以挨打次数为 X 轴、哼唱程度(输出)为 Y 轴的曲线: --"呜呼!这条曲线就像一座山,想不通,想不通。为什么那个无赖被打了三十天也不喊救命? --哦,你打的时间间隔(Δτ=24小时)太长了,这样无赖一天承受的痛苦程度,没有叠加,始终是个常数;如果缩短时间间隔(建议Δτ=0。5 秒),那么他的疼痛程度就可以迅速叠加;等到无赖挨了三十下(t=30)时,疼痛程度已经达到他叫喊能力的极限,就会收到最好的惩戒效果,再多挨几下也不会手下留情。 --还是不太明白,为什么疼痛程度会在小时间间隔内叠加? --这跟人(线性时变系统)对木板(脉冲、输入、激发)的反应有关。什么是响应?人收到板子后,疼痛的感觉会在一天内(假设,因人而异)慢慢消失(衰减),而不是突然消失。这样,只要中风的时间间隔较小,每次中风造成的疼痛就没有时间完全衰减,都会对最终的疼痛程度产生不同的影响: t 块大板造成的疼痛程度 = Σ(第 τ 块大板造成的疼痛程度 * 衰减系数)[衰减系数是 (t - τ) 的函数,请仔细品味] 数学表达式为:y(t) = ∫T(τ)H(t-τ)
-
在使用MATLAB进行图像滤波时如何处理边缘问题
-
[Halcon&拟合] 拟合直线边缘并计算距离-图像预处理: 一般是去噪或抠图(blob分析抠图或手绘ROI区域抠图)两方面 轮廓提取: 1)boundary:区域轮廓提取 2)edges_sub_pix:图像轮廓提取 3)threshold_sub_pix:图像轮廓提取 使用算子edges_sub_pix进行亚像素的边缘提取最为普遍。其用到的滤波器有Deriche, Lanser, Shen, or Canny filters。 关于这几个滤波器的对比,帮助文档有如下介绍: Deriche, Lanser, Shen为递归滤波器,Canny 为掩膜滤波器; 递归滤波器的执行时间不依赖滤波器的大小,Canny的执行时间与滤波器大小成正相关。 参数alpha数值越大,Deriche, Lanser, Shen滤波器宽度越小,平滑越差,细节越突出,而Canny效果相反。 分割、联合(根据情况而定) 分割算子: segment_contours_xld:可分割’lines’,‘lines_circles’,‘lines_ellipses’,原理是多边形逼近,逼近程度通过算子中后两个阀值参数控制。 联合算子: 临近:union_adjacent_contours_xld (Operator) 共线:union_collinear_contours_xld (Operator) 共圆:union_cocircular_contours_xld (Operator) 拟合 fit_line_contour_xld:拟合直线 fit_line_contour_xld:拟合圆 fit_ellipse_contour_xld:拟合椭圆 fit_rectangle2_contour_xld:拟合矩形 注:有时候在拟合轮廓之前需要判断一下轮廓属性,以确定应拟合成直线还是还是圆,可通过算子:get_contour_global_attrib_xld (SingleSegment,‘cont_approx’,)名字:获取轮廓属性描述:用于确定应拟合成直线还是还是圆参数:SingleSegment:输入轮廓(input_object)cont_approx:属性名称,即采用什么方式去计算 ,一般用这个参数就可以了(input_control)Attrib:属性值: Attrib>0:拟合圆,否则拟合直线(output_control) ) 求距离 二、示例: