Source

pyrankfilter / pyrankfilter / rankfilter.py

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
# -*- coding: utf-8 -*-
'''
*Test the scipy weave module to the filt function (rankfilter) 8 and 10 bit per pixel*

.. note: use Weave for C compilation

the following example read an image from file, processed it applying a median filter (radius = 30):

.. code-block:: python

   from scipy import misc
   from pyrankfilter import rankfilter

   im = misc.imread('../test/data/cameraman.tif')
   f = rankfilter(im,filtName = 'median',radius = 30,verbose = False)
   print f
'''
__author__ = 'Copyright (C) 2012, Olivier Debeir <odebeir@ulb.ac.be>'
__license__ ="""
pyrankfilter is a python module that implements 2D numpy arrays rank filters, the filter core is C-code
compiled on the fly (with an ad-hoc kernel).

Copyright (C) 2012  Olivier Debeir

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""
import numpy as npy
from scipy.weave import inline
import os.path as path


def _dtype2ctype(array):
    """convert numpy type in C equivalent type
    """
    types = {
        npy.dtype(npy.float64): 'double',
        npy.dtype(npy.float32): 'float',
        npy.dtype(npy.int32): 'int',
        npy.dtype(npy.int16): 'short',
        npy.dtype(npy.uint8): 'unsigned char',
        npy.dtype(npy.uint16): 'unsigned short',
    }
    return types.get(array.dtype)

kernel = {}
kernel['volume'] = '*ptemp=pop;'
kernel['meansubstraction'] = 'sum=0.0f; for(m=0;m<HISTOSIZE;m++) sum+=(m*histo[m]); *ptemp=(((float)*p-sum/pop)/2.0f+MEDIF);'
kernel['egalise'] = 'sum=0.0f; for(m=0;m<=*p;m++) sum+=histo[m]; *ptemp=sum*MAXIF/pop;'
kernel['minimum'] = 'for(m=0;(m<HISTOSIZE)&&(histo[m]==0);m++); *ptemp=m;'
kernel['maximum'] = 'for(m=HISTOSIZE;(--m)&&(histo[m]==0);); *ptemp=m;'
kernel['tophat'] = 'for(m=HISTOSIZE;(--m)&&(histo[m]==0);); *ptemp=m-*p;'
kernel['bottomhat'] = 'for(m=0;(m<HISTOSIZE)&&(histo[m]==0);m++); *ptemp=*p-m;'
kernel['gradient'] = 'for(m=HISTOSIZE;(--m)&&(histo[m]==0);); max=m; for(m=0;(m<HISTOSIZE)&&(histo[m]==0);m++); min=m; *ptemp=max-min;'
kernel['morph_contr_enh'] = 'for(m=HISTOSIZE;(--m)&&(histo[m]==0);); max=m; for(m=0;(m<HISTOSIZE)&&(histo[m]==0);m++); min=m; if(abs(*p-min)<abs(*p-max)) *ptemp=min; else *ptemp=max;'
kernel['autolevel'] = 'for(m=HISTOSIZE;(--m)&&(histo[m]==0);); max=m; for(m=0;(m<HISTOSIZE)&&(histo[m]==0);m++); min=m; if(max-min){*ptemp=((float)(*p-min)/(float)(max-min)*MAXIF);}else{*ptemp=(float)*p;}'
kernel['lowest'] = 'for(m=0;(m<HISTOSIZE)&&(histo[m]==0);m++); min=m; if(abs(*p-min)<=delta)*ptemp=MAXI;'
kernel['highest'] = 'for(m=HISTOSIZE;(--m)&&(histo[m]==0);); max=m; if(abs(*p-max)<=delta)*ptemp=MAXI;'
kernel['extrema'] = 'for(m=HISTOSIZE;(--m)&&(histo[m]==0);); max=m; for(m=0;(m<HISTOSIZE)&&(histo[m]==0);m++); min=m; if((abs(*p-min)<=delta)||(abs(*p-max)<=delta))*ptemp=MAXI;'
kernel['threshold'] = 'sum=0.0f;for(m=0;m<HISTOSIZE;m++)sum+=(m*histo[m]); if(*p+delta>sum/pop)*ptemp=MAXI;'
kernel['median_threshold'] = 'sum=0; for(m=0;m<HISTOSIZE;m++){sum+=(float)histo[m]; if(sum>pop*0.5){if(*p+delta>m)*ptemp=MAXI;m=HISTOSIZE; }};'
kernel['rank']   = 'sum=0; for(m=0;m<HISTOSIZE;m++){sum+=(float)histo[m]; if(sum>pop*rank){*ptemp=m;m=HISTOSIZE; }};'
kernel['median'] = 'sum=0; for(m=0;m<HISTOSIZE;m++){sum+=(float)histo[m]; if(sum>=pop*.5){*ptemp=m;m=HISTOSIZE; }};'
kernel['modal'] = 'max=0; pmax=0; for(m=0;m<HISTOSIZE;m++){if(max<histo[m]){max=histo[m]; pmax=m; }}; *ptemp=pmax;'
kernel['soft_morph_contr_enh'] = 'sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*inf){min=m; m=HISTOSIZE; }}; sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*sup){max=m; m=HISTOSIZE; }} if(abs(*p-min)<abs(*p-max)) *ptemp=min; else *ptemp=max;'
kernel['soft_autolevel'] = 'sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*inf){min=m; m=HISTOSIZE;}}; sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*sup){max=m; m=HISTOSIZE; }}; if(max-min){f=((float)(*p-min)*MAXIF)/(float)(max-min); f=MAX(0,MIN(MAXIF,f)); *ptemp=f;}else{*ptemp=(float)*p;}'
kernel['soft_gradient'] = 'sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*inf){min=m; m=HISTOSIZE;}}; sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*sup){max=m; m=HISTOSIZE; }}; f=(float)(max-min); *ptemp=f;'
kernel['soft_mean'] = 'sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*inf){min=m; m=HISTOSIZE;}}; sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*sup){max=m; m=HISTOSIZE; }}; sum=0; count=0; for(m=min;m<max;m++){sum+=(m*histo[m]); count+=histo[m]; }; if(pop){f=(float)(sum)/(float)count;}else{f=0; }; *ptemp=f;'
kernel['soft_meansubstraction'] = 'sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*inf){min=m; m=HISTOSIZE; }}; sum=0; for(m=0;m<HISTOSIZE;m++){sum+=histo[m]; if(sum>pop*sup){max=m; m=HISTOSIZE; }}; sum=0; count=0; for(m=min;m<max;m++){sum+=(m*histo[m]); count+=histo[m]; }; if(pop){f=((float)(*p)-(float)(sum)/(float)count)/2.0f+MEDIF; }else{f=0; }; *ptemp=f;'
#spectral kernels
kernel['spectral_mean'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); stot=0; sum=0.0f; for(m=sga;m<=sgb;m++){sum+=(m*histo[m]); stot+=histo[m];};*ptemp=(float)sum/(float)stot;'
kernel['spectral_median'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); stot=0; sum=0; for(m=sga;m<=sgb;m++){stot+=histo[m];}; for(m=sga;m<=sgb;m++){sum+=(float)histo[m];if(sum>=(float)stot*.5){*ptemp=m;m=HISTOSIZE; }};'
kernel['spectral_gradient'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); for(m=sgb+1;(--m>=sga)&&(histo[m]==0);); max=m; for(m=sga;(m<=sgb)&&(histo[m]==0);m++); min=m; *ptemp=max-min;'
kernel['spectral_morph_contr_enh'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); for(m=sgb;(--m>=sga)&&(histo[m]==0);); max=m; for(m=sga;(m<=sgb)&&(histo[m]==0);m++); min=m; if(abs(*p-min)<abs(*p-max)) *ptemp=min; else *ptemp=max;'
kernel['spectral_minimum'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); for(m=sga;(m<=sgb)&&(histo[m]==0);m++); *ptemp=m;'
kernel['spectral_maximum'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); for(m=sgb+1;(--m>=sga)&&(histo[m]==0);); *ptemp=m;'
kernel['spectral_volume'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); stot=0; for(m=sga;m<=sgb;m++){stot+=histo[m];}; *ptemp=255*(float)stot/(float)pop;'
kernel['spectral_autolevel'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); for(m=sgb+1;(--m>=sga)&&(histo[m]==0);); max=m; for(m=sga;(m<=sgb)&&(histo[m]==0);m++); min=m; if(max-min){*ptemp=((float)(*p-min)/(float)(max-min)*MAXIF);}else{*ptemp=(float)*p;}'
kernel['spectral_lowest'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); for(m=sga;(m<=sgb)&&(histo[m]==0);m++); min=m; if(abs(*p-min)<=delta)*ptemp=MAXI;'
kernel['spectral_highest'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); for(m=sgb+1;(--m>=sga)&&(histo[m]==0);); max=m; if(abs(*p-max)<=delta)*ptemp=MAXI;'
#spectral soft kernels
kernel['spectral_soft_autolevel'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); for(m=sgb+1;--m&&(histo[m]==0);); max=m; for(m=sga;(m<=sgb)&&(histo[m]==0);m++); min=m; if(max-min){f=((float)(*p-min)*MAXIF)/(float)(max-min); f=MAX(0,MIN(MAXIF,f)); *ptemp=f; }else{*ptemp=(float)*p; }'
kernel['spectral_soft_morph_contr_enh'] = 'sga=MAX(0,*p-sa); sgb=MIN(MAXIF,*p+sb); stot=0; for(m=sga;m<=sgb;m++){stot+=histo[m]; }; sum=0; for(m=sga;m<=sgb;m++){sum+=histo[m]; if(sum>(float)stot*inf){min=m; m=HISTOSIZE; }}; sum=0; for(m=sga;m<=sgb;m++){sum+=histo[m]; if(sum>(float)stot*sup){max=m; m=HISTOSIZE; }} if(abs(*p-min)<abs(*p-max)) *ptemp=min; else *ptemp=max;'
#statistical kernels
kernel['mean'] = 'sum=0.0f; for(m=0;m<HISTOSIZE;m++) sum+=(m*histo[m]); if (pop){*ptemp=sum/pop;}'
kernel['mean2'] = 'sum=0.0f; for(m=0;m<HISTOSIZE;m++) sum+=(m*m*histo[m]); if (pop){*ptemp=sum/pop;}'
kernel['entropy'] = 'sum=0.0f; for(m=0;m<HISTOSIZE;m++){if(histo[m]){sum-=((float)histo[m]*log((float)histo[m]));}}  *ptemp=sum;'

def kernel_list():
    """
    :param: no param   
    :returns: list of the available kernel for the rank filter function
    :rtype: list of stings

    returns the available kernel for the rank filter function

    example:

    >>> from pyrankfilter import kernel_list
    """
    list = [fname for fname in kernel]
    list.sort()
    return list

def rankfilter(ima, filtName, radius=None, struct_elem=None, struct_elem_center=None, infSup = (.1,.9),
               mask = None,bitDepth = 8,verbose = False, spectral_interval = (5,5),force = False, float_output = False):

    if not  kernel.has_key(filtName):
        raise ValueError('unknown filter name %s. Use kernelList() for a complete list.'%(filtName))

    c_kernel = kernel[filtName]

    return rankfilter_raw(ima, c_kernel, radius=radius, struct_elem=struct_elem, struct_elem_center=struct_elem_center,
        infSup = infSup,mask = mask,bitDepth = bitDepth,verbose = verbose, spectral_interval = spectral_interval,
        force = force, float_output = float_output)

def rankfilter_raw(ima, c_kernel, radius=None, struct_elem=None, struct_elem_center=None, infSup = (.1,.9),
               mask = None,bitDepth = 8,verbose = False, spectral_interval = (5,5),force = False, float_output = False):
    """compute the rank filtered 2D image on a radius wide circular window
    could eventually be limited on a given mask
    border of the image are filtered too
    the list of the available filters is given using kernelList()

    :param ima: image array
    :type ima: uint8 or uint16 numpy.ndarray
    :param c_kernel: string containing the kernel C-code (single line without newlines)
    :type filtName: string
    :param bitDepth: number of bits used max 16
    :type bitDepth: int[1-16]
    :param radius: radius of the neighborhood used
    :type radius: float [1-100]
    :param struct_elem: 2D image (if None : disk is used)
    :type struct_elem: numpy.ndarray (max size is 100x100)
    :param struct_elem_center: 2D image (non zero values are the element)
    :type struct_elem_center: int tuple(,). If None, the center is center of the structuring element
    :param infSup: inf and sup value used for 'soft' filters e.g. (.01,.99)
    :type infSup: float tuple (,)
    :param mask: mask image (if None : complete image area is filtered)
    :type mask: bool numpy.ndarray
    :param verbose: True : displays details during C code executes (default is False)
    :type verbose: Bool
    :param spectral_interval: if (a,b) the interval [g-a,g+b] is used to limit the neighborhood to similar points (w.r.t. gray level g) this option is only used with 'spectral' kernels (e.g. spectral_mean)
    :type infSup: float tuple (,)
    :param force: force the rankfilter re-compilation. If force is False, cached version, if present) will be used.
    :type force: bool
    :param float_output: False : force the result to be DOUBLE precision float
    :type float_output: Bool

    :returns: filtered image (same size and type as the original image)
    :raises: ValueError


    example 1: with a circular structuring element

    >>> from scipy import misc
    >>> import numpy as npy
    >>> from pyrankfilter import rankfilter
    >>> im = npy.zeros((100,100),dtype='uint8')
    >>> im[0:2,0:2] = 255
    >>> f = rankfilter(im,filtName = 'mean',radius = 3)
    >>> print f[:10,:10]
    [[92 72 60 42 14  0  0  0  0  0]
     [72 56 46 33 11  0  0  0  0  0]
     [60 46 37 18  0  0  0  0  0  0]
     [42 33 18  8  0  0  0  0  0  0]
     [14 11  0  0  0  0  0  0  0  0]
     [ 0  0  0  0  0  0  0  0  0  0]
     [ 0  0  0  0  0  0  0  0  0  0]
     [ 0  0  0  0  0  0  0  0  0  0]
     [ 0  0  0  0  0  0  0  0  0  0]
     [ 0  0  0  0  0  0  0  0  0  0]]


    example 2: with a custom structuring element

    >>> from scipy import misc
    >>> import numpy as npy
    >>> from pyrankfilter import rankfilter
    >>> im = npy.zeros((100,100),dtype='uint8')
    >>> im[1,1] = 255
    >>> elem = npy.asarray([[0,1,0],[1,1,1],[0,1,0]])
    >>> f = rankfilter(im,filtName = 'maximum',struct_elem = elem)
    >>> print f[:10,:10]
    [[  0 255   0   0   0   0   0   0   0   0]
     [255 255 255   0   0   0   0   0   0   0]
     [  0 255   0   0   0   0   0   0   0   0]
     [  0   0   0   0   0   0   0   0   0   0]
     [  0   0   0   0   0   0   0   0   0   0]
     [  0   0   0   0   0   0   0   0   0   0]
     [  0   0   0   0   0   0   0   0   0   0]
     [  0   0   0   0   0   0   0   0   0   0]
     [  0   0   0   0   0   0   0   0   0   0]
     [  0   0   0   0   0   0   0   0   0   0]]

    example 3: returning a float value (using custom structuring element)

    >>> from scipy import misc
    >>> import numpy as npy
    >>> from pyrankfilter import rankfilter
    >>> im = npy.zeros((100,100),dtype='uint8')
    >>> im[1,1] = 255
    >>> elem = npy.asarray([[0,1,0],[1,1,1],[0,1,0]])
    >>> f = rankfilter(im,filtName = 'mean',struct_elem = elem,float_output=True)
    >>> #f = rankfilter(im,filtName = 'mean',radius=2,float_output=True)
    >>> print f[:5,:5]
    [[   0.   255.     0.     0.     0. ]
     [ 255.   127.5  127.5    0.     0. ]
     [   0.   127.5    0.     0.     0. ]
     [   0.     0.     0.     0.     0. ]
     [   0.     0.     0.     0.     0. ]]

    example 4: returning a float value (using a disk)

    >>> f = rankfilter(im,filtName = 'mean',radius=2,float_output=True)
    >>> print f[:5,:5]
    [[ 42.5         31.875       28.33333397   0.           0.        ]
     [ 31.875       23.18181801  21.25        21.25         0.        ]
     [ 28.33333397  21.25        19.61538506   0.           0.        ]
     [  0.          21.25         0.           0.           0.        ]
     [  0.           0.           0.           0.           0.        ]]
     """

    if struct_elem is None:
        return _rankfilter_disk_raw(ima, c_kernel, radius, infSup=infSup, mask = mask,bitDepth = bitDepth,verbose = verbose,
            spectral_interval = spectral_interval,force = force, float_output = float_output)
    else:
        return _rankfilter_general_raw(ima, c_kernel, struct_elem=struct_elem, struct_elem_center=struct_elem_center,
            infSup = infSup, mask = mask, bitDepth = bitDepth, verbose = verbose, spectral_interval = spectral_interval,
            force = force, float_output = float_output)

def _rankfilter_disk(ima,filtName,radius,infSup = (.1,.9), mask = None,bitDepth = 8,verbose = False,
                     spectral_interval = (5,5),force = False, float_output = False):

    if not  kernel.has_key(filtName):
        raise ValueError('unknown filter name %s. Use kernelList() for a complete list.'%(filtName))

    c_kernel = kernel[filtName]
    return _rankfilter_disk_raw(ima,c_kernel,radius,infSup = infSup, mask = mask,bitDepth = bitDepth,verbose = verbose,
        spectral_interval = spectral_interval,force = force, float_output = float_output)

def _rankfilter_disk_raw(ima,c_kernel,radius,infSup = (.1,.9), mask = None,bitDepth = 8,verbose = False,
               spectral_interval = (5,5),force = False, float_output = False):
    """compute the rank filtered 2D image on a radius wide circular window (see rankfilter)
    c_kernel is the C-code snippet to be used
     """

    if not isinstance(ima,npy.ndarray):
        raise TypeError('2D numpy.array expected')
    if not ( (ima.dtype == npy.uint8) | (ima.dtype == npy.uint16)):
        raise TypeError('uint8 or uint16 numpy.array expected')
    if not(len(ima.shape) == 2):
        raise TypeError('2D numpy.array expected')

    DATATYPE =  _dtype2ctype(ima)
    if float_output:
        DATATYPE_OUT = 'double'
    else:
        DATATYPE_OUT = DATATYPE
    MAXRADIUS = 100

    #parameters
    if  ((radius > MAXRADIUS) | (radius<1)) :
        raise ValueError('radius must be [0,%d]'%MAXRADIUS)

    (infRange,supRange) = infSup
    (infSpectralRange,supSpectralRange) = spectral_interval

    if ima.dtype == npy.uint8:
        HISTOSIZE = 256
        MAXI = 255
        MEDI = 127
        MAXIF = '255.0f'
        MEDIF = '127.0f'
    else:
        if bitDepth not in range(0,17):
            raise ValueError('bitDepth must be in [0,16] range')
        HISTOSIZE = 2**bitDepth
        MAXI = HISTOSIZE-1
        MEDI = HISTOSIZE/2 - 1
        MAXIF = '%d.0f'%MAXI
        MEDIF = '%d.0f'%MEDI
        ima = npy.minimum(ima,HISTOSIZE-1) #avoid value >10bit

    #extend image border and create mask
    (m,n) = ima.shape
    extIn   = npy.zeros((m+2*radius, n+2*radius),dtype = ima.dtype, order='C')
    extMask = npy.zeros((m+2*radius, n+2*radius),dtype = ima.dtype, order='C')
    if float_output:
        extOut  = npy.zeros((m+2*radius, n+2*radius),dtype = float, order='C')
    else:
        extOut  = npy.zeros((m+2*radius, n+2*radius),dtype = ima.dtype, order='C')
    extIn[radius:m+radius,radius:n+radius] = ima

    if mask is not None:
        if not isinstance(mask,npy.ndarray):
            raise TypeError('2D numpy.array expected for mask')
        if not( ima.shape == mask.shape):
            raise ValueError('ima(%dx%d) and mask(%dx%d) must have same sizes.'%(ima.shape[0],ima.shape[1],mask.shape[0],mask.shape[1]))
        extMask[radius:m+radius,radius:n+radius] = mask
    else: #default mask covers the entire image
        extMask[radius:m+radius,radius:n+radius] = 1

    #activate or de-activate verbose in C code
    if verbose:
        VERBOSE = 'True'
    else:
        VERBOSE = ''

    #code contains the main circular window displacement loop
    #specific rank filter function is injected via the _PROCESSING_ preprocessing directive

    params = {'verbose':VERBOSE,
              'datatype':DATATYPE,
              'datatype_out':DATATYPE_OUT,
              'histosize':str(HISTOSIZE),
              'histosize_1':str(HISTOSIZE),
              'maxi':str(MAXI),
              'medi':str(MEDI),
              'maxif':MAXIF,
              'medif':MEDIF,
              'maxradius':str(MAXRADIUS),
              'processing':c_kernel}

    #read C code from .c file

    try:
        import pkgutil
        c_code = pkgutil.get_data(__name__, 'c-code/core_default.c')
    except ImportError:
        import pkg_resources
        c_code = pkg_resources.resource_string(__name__, 'c-code/core_default.c')

    #inject parameters(#define, etc) into the C source code
    code = c_code%params

    if verbose:
        print 'C code [%s] ="%s"\nprocess = "%s"'%(filtName,code,kernel[filtName])

    # compile C code
    inline(code,
        ['extIn','extMask','extOut','radius','infRange','supRange','infSpectralRange','supSpectralRange'],
        force=force)

    return extOut[radius:m+radius,radius:n+radius]

def _rankfilter_general(ima, filtName, struct_elem, struct_elem_center=None,infSup = (.1,.9), mask = None, bitDepth = 8,
                        verbose = False, spectral_interval = (5,5), force = False, float_output = False):
    if not  kernel.has_key(filtName):
        raise ValueError('unknown filter name %s. Use kernelList() for a complete list.'%(filtName))

    c_kernel = kernel[filtName]
    return _rankfilter_general_raw(ima, c_kernel, struct_elem=struct_elem, struct_elem_center=struct_elem_center,
        infSup = infSup, mask = mask, bitDepth = bitDepth,
        verbose = verbose, spectral_interval = spectral_interval, force = force, float_output = float_output)

def _rankfilter_general_raw(ima, c_kernel, struct_elem, struct_elem_center=None,infSup = (.1,.9), mask = None, bitDepth = 8,
                      verbose = False, spectral_interval = (5,5), force = False, float_output = False):
    """compute the rank filtered 2D image on a custom structuring element (see rankfilter)
    """

    if not isinstance(ima,npy.ndarray):
        raise TypeError('2D numpy.array expected')
    if not ( (ima.dtype == npy.uint8) | (ima.dtype == npy.uint16)):
        raise TypeError('uint8 or uint16 numpy.array expected')
    if not(len(ima.shape) == 2):
        raise TypeError('2D numpy.array expected')

    M_elem,N_elem = struct_elem.shape

    if struct_elem_center is None:
        center_m = M_elem/2
        center_n = N_elem/2
    else:
        center_m,center_n = struct_elem_center

    DATATYPE =  _dtype2ctype(ima)

    if float_output:
        DATATYPE_OUT = 'double'
    else:
        DATATYPE_OUT = DATATYPE

    MAXEDGE = 200

    #parameters
    if  (struct_elem.shape[0] > MAXEDGE) | (struct_elem.shape[1] > MAXEDGE) :
        raise ValueError('struct_elem.shape must be < [%d,%d]'%(MAXEDGE,MAXEDGE))

    (infRange,supRange) = infSup
    (infSpectralRange,supSpectralRange) = spectral_interval

    if ima.dtype == npy.uint8:
        HISTOSIZE = 256
        MAXI = 255
        MEDI = 127
        MAXIF = '255.0f'
        MEDIF = '127.0f'
    else:
        if bitDepth not in range(0,17):
            raise ValueError('bitDepth must be in [0,16] range')
        HISTOSIZE = 2**bitDepth
        MAXI = HISTOSIZE-1
        MEDI = HISTOSIZE/2 - 1
        MAXIF = '%d.0f'%MAXI
        MEDIF = '%d.0f'%MEDI
        ima = npy.minimum(ima,HISTOSIZE-1) #avoid value >10bit

    #extend image border and create mask
    (m,n) = ima.shape
    extIn   = npy.zeros((m+M_elem-1, n+N_elem-1),dtype = ima.dtype, order='C')
    extMask = npy.zeros((m+M_elem-1, n+N_elem-1),dtype = ima.dtype, order='C')
    if float_output:
        extOut  = npy.zeros((m+M_elem-1, n+N_elem-1),dtype = float, order='C')
    else:
        extOut  = npy.zeros((m+M_elem-1, n+N_elem-1),dtype = ima.dtype, order='C')
    extIn[center_m:center_m+ima.shape[0],center_n:center_n+ima.shape[1]] = ima

    struct_elem = struct_elem.astype('uint8')

    if mask is not None:
        if not isinstance(mask,npy.ndarray):
            raise TypeError('2D numpy.array expected for mask')
        if not( ima.shape == mask.shape):
            raise ValueError('ima(%dx%d) and mask(%dx%d) must have same sizes.'%(ima.shape[0],ima.shape[1],mask.shape[0],mask.shape[1]))
        extMask[center_m:center_m+ima.shape[0],center_n:center_n+ima.shape[1]] = mask
    else: #default mask covers the entire image
        extMask[center_m:center_m+ima.shape[0],center_n:center_n+ima.shape[1]] = 1

    #activate or de-activate verbose in C code
    if verbose:
        VERBOSE = 'True'
    else:
        VERBOSE = ''

    #code contains the main circular window displacement loop
    #specific rank filter function is injected via the _PROCESSING_ preprocessing directive

    params = {'verbose':VERBOSE,
              'datatype':DATATYPE,
              'datatype_out':DATATYPE_OUT,
              'histosize':str(HISTOSIZE),
              'histosize_1':str(HISTOSIZE),
              'maxi':str(MAXI),
              'medi':str(MEDI),
              'maxif':MAXIF,
              'medif':MEDIF,
              'maxedge':str(MAXEDGE**2),
              'processing':c_kernel}

    #read C code from .c file
    try:
        import pkgutil
        c_code = pkgutil.get_data(__name__, 'c-code/core_general.c')
    except ImportError:
        import pkg_resources
        c_code = pkg_resources.resource_string(__name__, 'c-code/core_general.c')

    #inject parameters(#define, etc) into the C source code
    code = c_code%params

    if verbose:
        print 'C code [%s] ="%s"\nprocess = "%s"'%(filtName,code,kernel[filtName])

    # compile C code
    inline(code,
        ['extIn','extMask','extOut','struct_elem','center_m','center_n','infRange','supRange','infSpectralRange','supSpectralRange'],
        force=force)

    return extOut[center_m:center_m+ima.shape[0],center_n:center_n+ima.shape[1]]




if __name__ == "__main__":
    #automatic testing with doctest module see >>> lines in the function doc
    print __license__
    import doctest
    doctest.testmod()
#    print help(rankfilter)