pyplotsuite/imageanalyzer/integral.py

   1 #
   2 #   Copyright (C) 2006-2007 Antonino Ingargiola <tritemio@gmail.com>
   3 #
   4 #   This file is part of PyPlotSuite <http://pyplotsuite.sourceforge.net>.
   5 #
   6 #   PyPlotSuite is free software; you can redistribute it and/or modify
   7 #   it under the terms of the GNU General Public License as published by
   8 #   the Free Software Foundation; either version 2 of the License, or
   9 #   (at your option) any later version.
  10
  11 from numpy import logical_and, logical_or, logical_not, arange, sqrt
  12
  13 def peak_integral1(a, t, xc, yc, r):
  14
  15     # Creating distance function using numpy broadcasting rules
  16     x = arange(a.shape[1])
  17     y = arange(a.shape[0]).reshape(a.shape[0],1)
  18     distance = sqrt((x-xc)**2 + (y-yc)**2)
  19
  20     # Creating domain and residual_domain for the integral
  21     domain_inside_circle = (distance < r)
  22     if t > 0:
  23         domain_above_threshold = (a > t)
  24         domain_below_threshold_abs = logical_and((a < t), (a > -t))
  25     else:
  26         # for t <= 0 select all the array
  27         domain_above_threshold = True
  28         domain_below_threshold_abs = True
  29     domain_index = logical_and(domain_above_threshold,domain_inside_circle)
  30     domain_size = domain_index.sum()
  31     domain = a[domain_index]
  32     residual_domain = \
  33             a[logical_and(logical_not(domain_inside_circle),
  34                     domain_below_threshold_abs)]
  35
  36     # Calculating integrals
  37     raw_integral = float(domain.sum())
  38     raw_residual = float(residual_domain.sum())
  39
  40     # Offset compensation
  41     offset = raw_residual / residual_domain.size
  42     offset_integral = raw_integral - offset*domain.size
  43
  44     return raw_integral, domain_size, raw_residual, offset_integral, offset
  45
  46 def peak_integral2t(a, t, rt, xc, yc, r):
  47
  48     # Creating distance function using numpy broadcasting rules
  49     x = arange(a.shape[1])
  50     y = arange(a.shape[0]).reshape(a.shape[0],1)
  51     distance = sqrt((x-xc)**2 + (y-yc)**2)
  52
  53     # Creating domain and residual_domain for the integral
  54     domain_inside_circle = (distance < r)
  55     if t > 0:
  56         domain_above_threshold = (a > t)
  57     else:
  58         # for t <= 0 select all the array
  59         domain_above_threshold = True
  60     if rt > 0:
  61         domain_below_threshold_abs = logical_and((a < rt), (a > -rt))
  62     else:
  63         # for rt <= 0 select all the array
  64         domain_below_threshold_abs = True
  65     domain_index = logical_and(domain_above_threshold,domain_inside_circle)
  66     domain_size = domain_index.sum()
  67     domain = a[domain_index]
  68     residual_domain = \
  69             a[logical_and(logical_not(domain_inside_circle),
  70                     domain_below_threshold_abs)]
  71
  72     # Calculating integrals
  73     raw_integral = float(domain.sum())
  74     raw_residual = float(residual_domain.sum())
  75
  76     # Offset compensation
  77     offset = raw_residual / residual_domain.size
  78     offset_integral = raw_integral - offset*domain.size
  79
  80     return raw_integral, domain_size, raw_residual, offset_integral, offset
  81
  82 if __name__ == "__main__":
  83     from PIL import Image
  84     from pylab import *
  85     import profile
  86     import time
  87     import scipy.ndimage.measurements as M
  88     from numpy import log
  89
  90     print 'Begin'
  91     fname='samples/test-img.tif'
  92     a = image2array(Image.open(fname)) + 0.1
  93     i = Image.open(fname)
  94     m = a.max()
  95     m2 = int(round(m+1))
  96     t1 = time.time()
  97     raw_integral, domain_size, raw_residual, offset_integral, offset = \
  98             peak_integral2(a, 335, 264, 227, 100)
  99     t2 = time.time()
 100     print 'tempo: ', t2-t1
 101     print a.dtype
 102