Improve performances for ring extraction

[kif]

x4 has been noticed
close #2078

[kif]

Line-profile before the optimisation:

Timer unit: 1e-09 s

Total time: 0.475535 s
File: /users/kieffer/.venv/py311/lib/python3.11/site-packages/pyFAI/goniometer.py
Function: extract_cp at line 651

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   651                                               def extract_cp(self, max_rings=None, pts_per_deg=1.0, Imin=0):
   652                                                   """Performs an automatic keypoint extraction and update the geometry refinement part
   653                                           
   654                                                   :param max_ring: extract at most N rings from the image
   655                                                   :param pts_per_deg: number of control points per azimuthal degree (increase for better precision)
   656                                                   """
   657         1        993.0    993.0      0.0          if self.massif is None:
   658                                                       if self.detector:
   659                                                           mask = self.detector.dynamic_mask(self.image)
   660                                                       else:
   661                                                           mask = None
   662                                                       self.massif = Massif(self.image, mask)
   663                                           
   664         1      24545.0  24545.0      0.0          tth = numpy.array([i for i in self.calibrant.get_2th() if i is not None])
   665         1     105481.0 105481.0      0.0          tth = numpy.unique(tth)
   666         1      17914.0  17914.0      0.0          tth_min = numpy.zeros_like(tth)
   667         1      14076.0  14076.0      0.0          tth_max = numpy.zeros_like(tth)
   668         1      23244.0  23244.0      0.0          delta = (tth[1:] - tth[:-1]) / 4.0
   669         1       2439.0   2439.0      0.0          tth_max[:-1] = delta
   670         1       1074.0   1074.0      0.0          tth_max[-1] = delta[-1]
   671         1       4540.0   4540.0      0.0          tth_min[1:] = -delta
   672         1       1015.0   1015.0      0.0          tth_min[0] = -delta[0]
   673         1       3444.0   3444.0      0.0          tth_max += tth
   674         1        819.0    819.0      0.0          tth_min += tth
   675         1        858.0    858.0      0.0          shape = self.image.shape
   676         1      11601.0  11601.0      0.0          ttha = self.geometry_refinement.twoThetaArray(shape)
   677         1       6083.0   6083.0      0.0          chia = self.geometry_refinement.chiArray(shape)
   678         1        202.0    202.0      0.0          rings = 0
   679         1      54696.0  54696.0      0.0          cp = ControlPoints(calibrant=self.calibrant)
   680         1        234.0    234.0      0.0          if max_rings is None:
   681         1        878.0    878.0      0.0              max_rings = tth.size
   682                                                   
   683         1       6693.0   6693.0      0.0          if self.detector.mask is not None:
   684         1     272706.0 272706.0      0.1              detector_mask = numpy.logical_not(self.geometry_refinement.detector.mask)
   685                                                   else:
   686                                                       detector_mask = numpy.zeros(self.detector.shape, dtype=bool)
   687         1       3206.0   3206.0      0.0          mask = numpy.empty_like(detector_mask)
   688         1        972.0    972.0      0.0          mask2 = numpy.empty_like(detector_mask)
   689         2    1152219.0 576109.5      0.2          ms = marchingsquares.MarchingSquaresMergeImpl(ttha,
   690         1       1594.0   1594.0      0.0                                                        mask=self.geometry_refinement.detector.mask,
   691         1        189.0    189.0      0.0                                                        use_minmax_cache=True)
   692        50      22151.0    443.0      0.0          for i in range(tth.size):
   693        49      16112.0    328.8      0.0              if rings >= max_rings:
   694                                                           break
   695        49   40882052.0 834327.6      8.6              numpy.logical_and(ttha >= tth_min[i], ttha < tth_max[i], out=mask)
   696        49     125229.0   2555.7      0.0              if self.detector.mask is not None:
   697        49    4665875.0  95221.9      1.0                  numpy.logical_and(mask, detector_mask, out=mask)
   698        49   28991788.0 591669.1      6.1              size = mask.sum(dtype=int)
   699        49      40560.0    827.8      0.0              if (size > 0):
   700         7       4764.0    680.6      0.0                  rings += 1
   701         7    2878095.0 411156.4      0.6                  sub_data = self.image.ravel()[numpy.where(mask.ravel())]
   702         7     651790.0  93112.9      0.1                  mean = sub_data.mean(dtype=numpy.float64)
   703         7    1963965.0 280566.4      0.4                  std = sub_data.std(dtype=numpy.float64)
   704         7       6529.0    932.7      0.0                  upper_limit = mean + std
   705         7    5634584.0 804940.6      1.2                  numpy.logical_and(self.image > upper_limit, mask, out=mask2)
   706         7    4298486.0 614069.4      0.9                  size2 = mask2.sum(dtype=int)
   707         7       7199.0   1028.4      0.0                  if size2 < 1000:
   708                                                               upper_limit = mean
   709                                                               numpy.logical_and(self.image > upper_limit, mask, out=mask2)
   710                                                               size2 = mask2.sum()
   711                                                           # length of the arc:
   712         7   13022253.0    2e+06      2.7                  points = ms.find_pixels(tth[i])
   713         7    3561853.0 508836.1      0.7                  seeds = set((i[0], i[1]) for i in points if mask2[i[0], i[1]])
   714                                                           # max number of points: 360 points for a full circle
   715         7     404648.0  57806.9      0.1                  azimuthal = chia[points[:, 0].clip(0, shape[0]), points[:, 1].clip(0, shape[1])]
   716         7     478693.0  68384.7      0.1                  nb_deg_azim = numpy.unique(numpy.rad2deg(azimuthal).round()).size
   717         7      13485.0   1926.4      0.0                  keep = int(nb_deg_azim * pts_per_deg)
   718         7       2159.0    308.4      0.0                  if keep == 0:
   719                                                               continue
   720         7      10485.0   1497.9      0.0                  dist_min = len(seeds) / 2.0 / keep
   721                                                           # why 3.0, why not ?
   722                                           
   723        14      54166.0   3869.0      0.0                  logger.info("Extracting datapoint for ring %s (2theta = %.2f deg); " +
   724                                                                       "searching for %i pts out of %i with I>%.1f, dmin=%.1f",
   725         7      24804.0   3543.4      0.0                              i, numpy.degrees(tth[i]), keep, size2, upper_limit, dist_min)
   726         7  365222962.0    5e+07     76.8                  res = self.massif.peaks_from_area(mask2, Imin=Imin, keep=keep, dmin=dist_min, seed=seeds, ring=i)
   727         7     275832.0  39404.6      0.1                  cp.append(res, i)
   728         1       1130.0   1130.0      0.0          self.control_points = cp
   729         1     561283.0 561283.0      0.1          self.geometry_refinement.data = numpy.asarray(cp.getList(), dtype=numpy.float64)
   730         1        281.0    281.0      0.0          return cp

Total time: 0.3563 s
File: /users/kieffer/.venv/py311/lib/python3.11/site-packages/pyFAI/massif.py
Function: peaks_from_area at line 211

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   211                                               def peaks_from_area(self, mask, Imin=numpy.finfo(numpy.float64).min,
   212                                                                   keep=1000, dmin=0.0, seed=None, **kwarg):
   213                                                   """
   214                                                   Return the list of peaks within an area
   215                                           
   216                                                   :param mask: 2d array with mask.
   217                                                   :param Imin: minimum of intensity above the background to keep the point
   218                                                   :param keep: maximum number of points to keep
   219                                                   :param kwarg: ignored parameters
   220                                                   :param dmin: minimum distance to another peak
   221                                                   :param seed: list of good guesses to start with
   222                                                   :return: list of peaks [y,x], [y,x], ...]
   223                                                   """
   224                                                   
   225         7       5709.0    815.6      0.0          res = []
   226         7       2382.0    340.3      0.0          cnt = 0
   227         7       2590.0    370.0      0.0          dmin2 = dmin * dmin
   228                                           
   229         7   16941895.0    2e+06      4.8          all_points = numpy.vstack(numpy.where(mask)).T
   230         7       7962.0   1137.4      0.0          if len(all_points) > 0:
   231         7   32494826.0    5e+06      9.1              numpy.random.shuffle(all_points)
   232                                                   
   233         7       3838.0    548.3      0.0          if seed:
   234         6     148350.0  24725.0      0.0              seeds = numpy.array(list(seed))
   235         6       3844.0    640.7      0.0              if len(seeds) > 0:
   236         6     409239.0  68206.5      0.1                  numpy.random.shuffle(seeds)
   237         6     187146.0  31191.0      0.1              all_points = numpy.concatenate((seeds, all_points))
   238                                           
   239      3012    1060518.0    352.1      0.3          for idx in all_points:
   240      3012   16487380.0   5473.9      4.6              out = self.nearest_peak(idx)
   241      3012     561993.0    186.6      0.2              if out is not None:
   242      3012     534516.0    177.5      0.2                  msg = "[ %3i, %3i ] -> [ %.1f, %.1f ]"
   243      3012    4651766.0   1544.4      1.3                  logger.debug(msg, idx[1], idx[0], out[1], out[0])
   244      3012    2415009.0    801.8      0.7                  p0, p1 = int(round(out[0])), int(round(out[1]))
   245      3012    1078631.0    358.1      0.3                  if mask[p0, p1]:
   246      2823  276943158.0  98102.4     77.7                      if (self.data[p0, p1] > Imin) and is_far_from_group(out, res, dmin2):
   247       565     247846.0    438.7      0.1                          res.append(out)
   248       565     104835.0    185.5      0.0                          cnt = 0
   249      3012    1234471.0    409.9      0.3              if len(res) >= keep or cnt > keep:
   250         7       2198.0    314.0      0.0                  break
   251                                                       else:
   252      3005     768127.0    255.6      0.2                  cnt += 1
   253         7       1362.0    194.6      0.0          return res

Total time: 0.187934 s
File: /users/kieffer/.venv/py311/lib/python3.11/site-packages/pyFAI/utils/mathutil.py
Function: is_far_from_group at line 719

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   719                                           def is_far_from_group(pt, lst_pts, d2):
   720                                               """
   721                                               Tells if a point is far from a group of points, distance greater than d2 (distance squared)
   722                                           
   723                                               :param pt: point of interest
   724                                               :param lst_pts: list of points
   725                                               :param d2: minimum distance squarred
   726                                               :return: True If the point is far from all others.
   727                                           
   728                                               """
   729     99233   16214394.0    163.4      8.6      for apt in lst_pts:
   730     98668  152082194.0   1541.4     80.9          dsq = sum((i - j) * (i - j) for i, j in zip(apt, pt))
   731     98668   19107071.0    193.7     10.2          if dsq <= d2:
   732      2258     432428.0    191.5      0.2              return False
   733       565      98260.0    173.9      0.1      return True

[kif]

And after this PR has been merged:

Timer unit: 1e-09 s

Total time: 0.116329 s
File: /users/kieffer/.venv/py311/lib/python3.11/site-packages/pyFAI/goniometer.py
Function: extract_cp at line 651

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   651                                               def extract_cp(self, max_rings=None, pts_per_deg=1.0, Imin=0):
   652                                                   """Performs an automatic keypoint extraction and update the geometry refinement part
   653                                           
   654                                                   :param max_ring: extract at most N rings from the image
   655                                                   :param pts_per_deg: number of control points per azimuthal degree (increase for better precision)
   656                                                   """
   657         1       2421.0   2421.0      0.0          if self.massif is None:
   658                                                       if self.detector:
   659                                                           mask = self.detector.dynamic_mask(self.image)
   660                                                       else:
   661                                                           mask = None
   662                                                       self.massif = Massif(self.image, mask)
   663                                           
   664         1      59689.0  59689.0      0.1          tth = numpy.array([i for i in self.calibrant.get_2th() if i is not None])
   665         1     163671.0 163671.0      0.1          tth = numpy.unique(tth)
   666         1      33360.0  33360.0      0.0          tth_min = numpy.zeros_like(tth)
   667         1      12448.0  12448.0      0.0          tth_max = numpy.zeros_like(tth)
   668         1      20553.0  20553.0      0.0          delta = (tth[1:] - tth[:-1]) / 4.0
   669         1       4864.0   4864.0      0.0          tth_max[:-1] = delta
   670         1       2316.0   2316.0      0.0          tth_max[-1] = delta[-1]
   671         1       8700.0   8700.0      0.0          tth_min[1:] = -delta
   672         1       2375.0   2375.0      0.0          tth_min[0] = -delta[0]
   673         1       6233.0   6233.0      0.0          tth_max += tth
   674         1       1918.0   1918.0      0.0          tth_min += tth
   675         1       1684.0   1684.0      0.0          shape = self.image.shape
   676         1      25829.0  25829.0      0.0          ttha = self.geometry_refinement.twoThetaArray(shape)
   677         1      12858.0  12858.0      0.0          chia = self.geometry_refinement.chiArray(shape)
   678         1        450.0    450.0      0.0          rings = 0
   679         1     120331.0 120331.0      0.1          cp = ControlPoints(calibrant=self.calibrant)
   680         1        525.0    525.0      0.0          if max_rings is None:
   681         1       1342.0   1342.0      0.0              max_rings = tth.size
   682                                                   
   683         1    5027806.0    5e+06      4.3          qmask, count = build_qmask(ttha, tth_min, tth_max, self.geometry_refinement.detector.mask) 
   684                                                   
   685                                                   # if self.detector.mask is not None:
   686                                                   #     detector_mask = numpy.logical_not(self.geometry_refinement.detector.mask)
   687                                                   # else:
   688                                                   #     detector_mask = numpy.zeros(self.detector.shape, dtype=bool)
   689                                                   # mask = numpy.empty_like(detector_mask)
   690         1       7728.0   7728.0      0.0          mask2 = numpy.empty(qmask.shape, dtype=bool)
   691         2    1703180.0 851590.0      1.5          ms = marchingsquares.MarchingSquaresMergeImpl(ttha,
   692         1       6900.0   6900.0      0.0                                                        mask=self.geometry_refinement.detector.mask,
   693         1        145.0    145.0      0.0                                                        use_minmax_cache=True)
   694        50      11634.0    232.7      0.0          for i in range(tth.size):
   695        49       8536.0    174.2      0.0              if rings >= max_rings:
   696                                                           break
   697                                                       # numpy.logical_and(ttha >= tth_min[i], ttha < tth_max[i], out=mask)
   698                                                       # if self.detector.mask is not None:
   699                                                       #     numpy.logical_and(mask, detector_mask, out=mask)
   700                                                       # size = count[i] #mask.sum(dtype=int)
   701        49      18388.0    375.3      0.0              if count[i]:
   702         7       2045.0    292.1      0.0                  rings += 1
   703         7    2137193.0 305313.3      1.8                  mask = qmask==i
   704         7    1645640.0 235091.4      1.4                  sub_data = self.image[mask]
   705         7     524770.0  74967.1      0.5                  mean = sub_data.mean(dtype=numpy.float64)
   706         7    1620892.0 231556.0      1.4                  std = sub_data.std(dtype=numpy.float64)
   707         7       5276.0    753.7      0.0                  upper_limit = mean + std
   708         7    5016760.0 716680.0      4.3                  numpy.logical_and(self.image > upper_limit, mask, out=mask2)
   709         7    4514142.0 644877.4      3.9                  size2 = mask2.sum(dtype=int)
   710         7       6459.0    922.7      0.0                  if size2 < 1000:
   711                                                               upper_limit = mean
   712                                                               numpy.logical_and(self.image > upper_limit, mask, out=mask2)
   713                                                               size2 = mask2.sum()
   714                                                           # length of the arc:
   715         7   11818323.0    2e+06     10.2                  points = ms.find_pixels(tth[i])
   716         7    3126439.0 446634.1      2.7                  seeds = set((i[0], i[1]) for i in points if mask2[i[0], i[1]])
   717                                                           # max number of points: 360 points for a full circle
   718         7     315693.0  45099.0      0.3                  azimuthal = chia[points[:, 0].clip(0, shape[0]), points[:, 1].clip(0, shape[1])]
   719         7     381053.0  54436.1      0.3                  nb_deg_azim = numpy.unique(numpy.rad2deg(azimuthal).round()).size
   720         7       9274.0   1324.9      0.0                  keep = int(nb_deg_azim * pts_per_deg)
   721         7       2470.0    352.9      0.0                  if keep == 0:
   722                                                               continue
   723         7       9059.0   1294.1      0.0                  dist_min = len(seeds) / 2.0 / keep
   724                                                           # why 3.0, why not ?
   725                                           
   726        14      34912.0   2493.7      0.0                  logger.info("Extracting datapoint for ring %s (2theta = %.2f deg); " +
   727                                                                       "searching for %i pts out of %i with I>%.1f, dmin=%.1f",
   728         7      19011.0   2715.9      0.0                              i, numpy.degrees(tth[i]), keep, size2, upper_limit, dist_min)
   729         7   77208628.0    1e+07     66.4                  res = self.massif.peaks_from_area(mask2, Imin=Imin, keep=keep, dmin=dist_min, seed=seeds, ring=i)
   730         7     212305.0  30329.3      0.2                  cp.append(res, i)
   731         1        669.0    669.0      0.0          self.control_points = cp
   732         1     451700.0 451700.0      0.4          self.geometry_refinement.data = numpy.asarray(cp.getList(), dtype=numpy.float64)
   733         1        191.0    191.0      0.0          return cp

Total time: 0.0712687 s
File: /users/kieffer/.venv/py311/lib/python3.11/site-packages/pyFAI/massif.py
Function: peaks_from_area at line 211

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   211                                               def peaks_from_area(self, mask, Imin=numpy.finfo(numpy.float64).min,
   212                                                                   keep=1000, dmin=0.0, seed=None, **kwarg):
   213                                                   """
   214                                                   Return the list of peaks within an area
   215                                           
   216                                                   :param mask: 2d array with mask.
   217                                                   :param Imin: minimum of intensity above the background to keep the point
   218                                                   :param keep: maximum number of points to keep
   219                                                   :param kwarg: ignored parameters
   220                                                   :param dmin: minimum distance to another peak
   221                                                   :param seed: list of good guesses to start with
   222                                                   :return: list of peaks [y,x], [y,x], ...]
   223                                                   """
   224                                                   
   225         7       4143.0    591.9      0.0          res = []
   226         7       1513.0    216.1      0.0          cnt = 0
   227         7       2300.0    328.6      0.0          dmin2 = dmin * dmin
   228                                           
   229         7   16040919.0    2e+06     22.5          all_points = numpy.vstack(numpy.where(mask)).T
   230         7       5838.0    834.0      0.0          if len(all_points) > 0:
   231         7   29550031.0    4e+06     41.5              numpy.random.shuffle(all_points)
   232                                                   
   233         7       3448.0    492.6      0.0          if seed:
   234         6     127058.0  21176.3      0.2              seeds = numpy.array(list(seed))
   235         6       3588.0    598.0      0.0              if len(seeds) > 0:
   236         6     362113.0  60352.2      0.5                  numpy.random.shuffle(seeds)
   237         6     163995.0  27332.5      0.2              all_points = numpy.concatenate((seeds, all_points))
   238                                           
   239         7       1566.0    223.7      0.0          msg = "[ %3i, %3i ] -> [ %.1f, %.1f ]"
   240      2491     741518.0    297.7      1.0          for idx in all_points:
   241      2491   10426753.0   4185.8     14.6              out = self.nearest_peak(idx)
   242      2491     387967.0    155.7      0.5              if out is not None:
   243                                                           
   244      2491    3142218.0   1261.4      4.4                  logger.debug(msg, idx[1], idx[0], out[1], out[0])
   245      2491    1653504.0    663.8      2.3                  p0, p1 = int(round(out[0])), int(round(out[1]))
   246      2491     696156.0    279.5      1.0                  if mask[p0, p1]:
   247      2335    6280320.0   2689.6      8.8                      if (self.data[p0, p1] > Imin) and is_far_from_group(out, res, dmin2):
   248       564     168572.0    298.9      0.2                          res.append(out)
   249       564      86185.0    152.8      0.1                          cnt = 0
   250      2491     895589.0    359.5      1.3              if len(res) >= keep or cnt > keep:
   251         7       2354.0    336.3      0.0                  break
   252                                                       else:
   253      2484     519871.0    209.3      0.7                  cnt += 1
   254         7       1193.0    170.4      0.0          return res