MATLAB中处理任意大图像的函数blockproc

rocwoods

MATLAB分块处理矩阵的函数很早就有了——blkproc，但缺点是blkproc必须一次性把要处理的矩阵全部导入内存中，这样大大限制了其应用范围，对一些超大型的图像就无能为力了。幸运的是随着MATLAB使用范围越来越广，Mathworks也紧跟用户需求，新版本的MATLAB推出了可以处理任意大图像的函数blockproc，其用法如下：

1. 
   
2. B = blockproc(A,[M N],fun)
   
3. B = blockproc(src_filename,[M N],fun)
   
4. B = blockproc(adapter,[M N],fun)
   
5. blockproc(...,param,val,...)
   

复制代码

A是要处理的图像矩阵，如果图像太大不能完全导入内存，也可以用图像文件名src_filename来表示。[M,N]是希望每次分块处理的矩阵大小，fun是函数句柄，即对每块矩阵的处理函数。
需要说明的是blockproc默认支持tiff/tif和jpeg2000格式的任意大图像处理，如果要读取其他格式的大图像需要针对该图像格式再写一个继承自MATLAB中ImageAdapter这个抽象类的子类adapter，来满足blockproc的输入要求。MATLAB帮助文档中有一个读取lan格式的LanAdapter示例类，大家可以参照那个格式来构造任意图像格式的Adapter类来实现blockproc函数对任意大图像文件的支持。
最后一种调用格式可以实现读取大图像文件，分块处理后再在指定路径写入处理后的图像文件，这个非常有用。下面给一个简单的例子,更多的用法希望大家讨论完善。

1. 
   
2. 
   
3. fun = @(block_struct) repmat(block_struct.data,5,5) ;%注意处理函数的输入必须是结构体，其data域是我们的矩阵数据，这是由blockproc分块后的机制决定
   
4. blockproc('moon.tif',[16 16],fun,'Destination','C:\moonmoon.tif');
   
5. imshow('C:\moonmoon.tif')
   

复制代码

可以看到分块处理后的效果，当然这是简单的把原图像分块（每一子块大小16*16）复制了25倍后的效果。

rocwoods

大家有没有试验一下，随便找个2G以上的tiff文件（手头没有可以利用这个函数把MATLAB自带的tiff文件repmat成2G左右）操作一下。我现在的疑问是，对于“操作本身”和“图像分块”彼此是独立的情形可以很方便利用blockproc，但是反之则不大方便，譬如对图像整体转置旋转等操作，如果简单用blockproc处理后得到的是每个块的转置和旋转。大家有好的办法没有？

bainhome

这个命令没用过，2G也怪吓人的...
感觉能不能先旋转或者转置，然后再用blockproc来分块？

rocwoods

问题难就难在一些大的图像（遥感、卫星影像可能达几G几十G）无法全部导入内存，因此先转置就比较困难了。在Mathworks新闻组我也问了这个问题，如下：
http://www.mathworks.com/matlabc ... /view_thread/298278
一个叫brendan的老外（不知道是不是Mathworks内部工程师）给出了他的办法，结果非常给力。简单说一下他的思路：
先构造一个“辅助图像矩阵”充当块计数器，利用blockproc遍历这个辅助图像的像素(其实对应的就是实际图像的每一个块)。充分利用用户自定义的那个函数，那个函数可以随心所欲，实现最后大图像相应的块处显示哪个图像。他给的那个自定义函数blockTranspose思路如下：
blockproc的分块机制会产生一个结构体，记录当前块的数据以及位置等信息，因此先构造一个100*100的辅助矩阵，分成100*100块，实际就是一个像素分一块，对应10000*10000大图像的100*100小块。blockTranspose接收辅助矩阵的每一像素，这个像素包含了其代表的原图像相应块的位置信息。然后调用imread函数实现读取大图像中相应位置的图像块实现转置。
当然Brendan也承认这种用法不同于blockproc的常规用法，甚至是剑走偏锋了，仔细看了他的回复后，不得不感叹思路之巧妙，自定义函数的用法之灵活，让人十分佩服！更加体会了思想有多远，MATLAB就能玩得有多神！下面给出我对moon.tif图像转置的示例（当然这个图像不大，可以全部读进内存再转置，这里是拿它来做例子，几个G的图像也可以按这种方法处理）

1. 
   
2. function GlobalBlocDemo
   
3. function output_block  = myfun(block_struct)
   
4. start_index = (block_struct.location - 1) * p + 1;
   
5. end_index = start_index + p-1;
   
6. output_rows = [start_index(1) end_index(1)];
   
7. output_cols = [start_index(2) end_index(2)];
   
8. input_block = imread('moon.tif','PixelRegion',{output_cols,output_rows});
   
9. output_block = input_block';
   
10. end
    
11. info = imfinfo('moon.tif');
    
12. p = gcd(info.Width,info.Height);
    
13. blockproc(ones(info.Width/p,info.Height/p),[1 1],@myfun,'Destination','C:\moon2.tif')
    
14. end
    

复制代码

运行完毕后可以看下结果：

1. 
   
2. imshow('moon.tif');
   
3. figure;imshow('C:\moon2.tif')
   

复制代码

qibbxxt

4# rocwoods
由于分块矩阵的转置，在处理子矩阵的转置的同时，还是考虑到块本身的转置，
所以很难实现，的确brendan的程序完美的解答了这一问题，
而且构造辅助矩阵的想法确实体现了对matlab函数的深刻的理解，
将图像的处理放在自定义的函数中，确实给人一种高屋建瓴之感。
感谢roc对这个问题的深入讨论，使得我们在分块矩阵的转置有了新的认识。
-----------------------------------------------------------------------------------------------------------
我看了下blockproc，对没有小块的处理顺序，刚开始的时候感觉有点乱
好在blockproc的一些代码是可见的，才发现

1. 
   
2. [r,c] = meshgrid(2:mblocks,2:nblocks);
   
3. rr = r(:);
   
4. cc = c(:);
   
5. [r,c] = meshgrid(1,2:end_col);
   
6. rr = [rr;r(:)];
   
7. cc = [cc;c(:)];
   
8. [r,c] = meshgrid(2:end_row,1);
   
9. rr = [rr;r(:)];
   
10. cc = [cc;c(:)];
    

复制代码

这些代码产生了执行小块的先后顺序，之所以这样，可以有通用性的考虑吧

qibbxxt

为了能够理解4#roc给出的程序，我写了一个矩阵转置的小例子，供大家参考
主程序：

1. clear;clc;close all
   
2. m=2;% 子矩阵的行数
   
3. n=4;% 子矩阵的列数
   
4. M=3;% 大矩阵的块的行数
   
5. N=2;% 大矩阵的块的列数
   
6. maxM=3; % 子矩阵的最大数
   
7. A=randi(maxM,[m,n]);% 生成子矩阵
   
8. %  为了索引子矩阵处理方便，我们建立一个辅助的单元数组
   
9. Bcell=cell(M,N);
   
10. Bcell(:)=arrayfun(@(x)x*A,1:6,'UniformOutput',false);
    
11. Bnum=cell2mat(Bcell);
    
12. C=largeMatrixTransposeDemo(Bcell,M,N);
    
13. disp('原矩阵');
    
14. disp(Bnum);
    
15. disp('用blockproc函数转置后的矩阵');
    
16. disp(C);
    
17. disp('用matlab的transpose转置的结果');
    
18. disp(Bnum');

复制代码

调用的函数

1. function C=largeMatrixTransposeDemo(B,M,N)
   
2.     function output_a  = myfun(input_a)
   
3.        pos=input_a.location;
   
4.        output_a=B{pos(2),pos(1)}';
   
5.     end
   
6. C=blockproc(zeros(N,M)+1,[1 1],@myfun);
   
7. end

复制代码

运行的结果：

1. 原矩阵
   
2.      3     3     3     2    12    12    12     8
   
3.      1     2     1     1     4     8     4     4
   
4.      6     6     6     4    15    15    15    10
   
5.      2     4     2     2     5    10     5     5
   
6.      9     9     9     6    18    18    18    12
   
7.      3     6     3     3     6    12     6     6
   
8. 
   
9. 用blockproc函数转置后的矩阵
   
10.      3     1     6     2     9     3
    
11.      3     2     6     4     9     6
    
12.      3     1     6     2     9     3
    
13.      2     1     4     2     6     3
    
14.     12     4    15     5    18     6
    
15.     12     8    15    10    18    12
    
16.     12     4    15     5    18     6
    
17.      8     4    10     5    12     6
    
18. 
    
19. 用matlab的transpose转置的结果
    
20.      3     1     6     2     9     3
    
21.      3     2     6     4     9     6
    
22.      3     1     6     2     9     3
    
23.      2     1     4     2     6     3
    
24.     12     4    15     5    18     6
    
25.     12     8    15    10    18    12
    
26.     12     4    15     5    18     6
    
27.      8     4    10     5    12     6

复制代码

mxlzhenzhu

qibbxxt 发表于 2010-12-9 17:42
为了能够理解4#roc给出的程序，我写了一个矩阵转置的小例子，供大家参考
主程序：调用的函数运行的结果： ...

大矩阵转置？
太好了，我正要。