title: Algorithm to locate local maxima
body: <p>I have data that always looks something like this:</p>
  
  <p><a href="http://michaelfogleman.com/static/images/chart.png">alt text http://michaelfogleman.com/static/images/chart.png</a></p>
  
  <p>I need an algorithm to locate the three peaks.</p>
  
  <p>The x-axis is actually a camera position and the y-axis is a measure of image focus/contrast at that position.  There are features at three different distances that can be in focus and I need to determine the x-values for these three points.</p>
  
  <p>The middle hump is always a little harder to pick out, even for a human.</p>
  
  <p>I have a homemade algorithm that mostly works, but I'm wondering if there's a standard way to grab local maxima from a function that can have a little noise in it.  The peaks overcome the noise easily though.</p>
  
  <p>Also, being camera data, an algorithm that doesn't require scanning the full range could be useful.</p>
  
  <p><strong>Edit</strong>: Posting the Python code that I ended up using.  It uses my original code that finds maxima given a search threshold and does a binary search to find a threshold that results in the desired number of maxima.</p>
  
  <p><strong>Edit</strong>: Sample data included in code below. New code is O(n) instead of O(n^2).</p>
  
  <pre><code>def find_n_maxima(data, count):
      low = 0
      high = max(data) - min(data)
      for iteration in xrange(100): # max iterations
          mid = low + (high - low) / 2.0
          maxima = find_maxima(data, mid)
          if len(maxima) == count:
              return maxima
          elif len(maxima) &lt; count: # threshold too high
              high = mid
          else: # threshold too low
              low = mid
      return None # failed
  
  def find_maxima(data, threshold):
      def search(data, threshold, index, forward):
          max_index = index
          max_value = data[index]
          if forward:
              path = xrange(index + 1, len(data))
          else:
              path = xrange(index - 1, -1, -1)
          for i in path:
              if data[i] &gt; max_value:
                  max_index = i
                  max_value = data[i]
              elif max_value - data[i] &gt; threshold:
                  break
          return max_index, i
      # forward pass
      forward = set()
      index = 0
      while index &lt; len(data) - 1:
          maximum, index = search(data, threshold, index, True)
          forward.add(maximum)
          index += 1
      # reverse pass
      reverse = set()
      index = len(data) - 1
      while index &gt; 0:
          maximum, index = search(data, threshold, index, False)
          reverse.add(maximum)
          index -= 1
      return sorted(forward &amp; reverse)
  
  data = [
      1263.900, 1271.968, 1276.151, 1282.254, 1287.156, 1296.513,
      1298.799, 1304.725, 1309.996, 1314.484, 1321.759, 1323.988,
      1331.923, 1336.100, 1340.007, 1340.548, 1343.124, 1353.717,
      1359.175, 1364.638, 1364.548, 1357.525, 1362.012, 1367.190,
      1367.852, 1376.275, 1374.726, 1374.260, 1392.284, 1382.035,
      1399.418, 1401.785, 1400.353, 1418.418, 1420.401, 1423.711,
      1425.214, 1436.231, 1431.356, 1435.665, 1445.239, 1438.701,
      1441.988, 1448.930, 1455.066, 1455.047, 1456.652, 1456.771,
      1459.191, 1473.207, 1465.788, 1488.785, 1491.422, 1492.827,
      1498.112, 1498.855, 1505.426, 1514.587, 1512.174, 1525.244,
      1532.235, 1543.360, 1543.985, 1548.323, 1552.478, 1576.477,
      1589.333, 1610.769, 1623.852, 1634.618, 1662.585, 1704.127,
      1758.718, 1807.490, 1852.097, 1969.540, 2243.820, 2354.224,
      2881.420, 2818.216, 2552.177, 2355.270, 2033.465, 1965.328,
      1824.853, 1831.997, 1779.384, 1764.789, 1704.507, 1683.615,
      1652.712, 1646.422, 1620.593, 1620.235, 1613.024, 1607.675,
      1604.015, 1574.567, 1587.718, 1584.822, 1588.432, 1593.377,
      1590.533, 1601.445, 1667.327, 1739.034, 1915.442, 2128.835,
      2147.193, 1970.836, 1755.509, 1653.258, 1613.284, 1558.576,
      1552.720, 1541.606, 1516.091, 1503.747, 1488.797, 1492.021,
      1466.720, 1457.120, 1462.485, 1451.347, 1453.224, 1440.477,
      1438.634, 1444.571, 1428.962, 1431.486, 1421.721, 1421.367,
      1403.461, 1415.482, 1405.318, 1399.041, 1399.306, 1390.486,
      1396.746, 1386.178, 1376.941, 1369.880, 1359.294, 1358.123,
      1353.398, 1345.121, 1338.808, 1330.982, 1324.264, 1322.147,
      1321.098, 1313.729, 1310.168, 1304.218, 1293.445, 1285.296,
      1281.882, 1280.444, 1274.795, 1271.765, 1266.857, 1260.161,
      1254.380, 1247.886, 1250.585, 1246.901, 1245.061, 1238.658,
      1235.497, 1231.393, 1226.241, 1223.136, 1218.232, 1219.658,
      1222.149, 1216.385, 1214.313, 1211.167, 1208.203, 1206.178,
      1206.139, 1202.020, 1205.854, 1206.720, 1204.005, 1205.308,
      1199.405, 1198.023, 1196.419, 1194.532, 1194.543, 1193.482,
      1197.279, 1196.998, 1194.489, 1189.537, 1188.338, 1184.860,
      1184.633, 1184.930, 1182.631, 1187.617, 1179.873, 1171.960,
      1170.831, 1167.442, 1177.138, 1166.485, 1164.465, 1161.374,
      1167.185, 1174.334, 1186.339, 1202.136, 1234.999, 1283.328,
      1347.111, 1679.050, 1927.083, 1860.902, 1602.791, 1350.454,
      1274.236, 1207.727, 1169.078, 1138.025, 1117.319, 1109.169,
      1080.018, 1073.837, 1059.876, 1050.209, 1050.859, 1035.003,
      1029.214, 1024.602, 1017.932, 1006.911, 1010.722, 1005.582,
      1000.332, 998.0721, 992.7311, 992.6507, 981.0430, 969.9936,
      972.8696, 967.9463, 970.1519, 957.1309, 959.6917, 958.0536,
      954.6357, 954.9951, 947.8299, 953.3991, 949.2725, 948.9012,
      939.8549, 940.1641, 942.9881, 938.4526, 937.9550, 929.6279,
      935.5402, 921.5773, 933.6365, 918.7065, 922.5849, 939.6088,
      911.3251, 923.7205, 924.8227, 911.3192, 936.7066, 915.2046,
      919.0274, 915.0533, 910.9783, 913.6773, 916.6287, 907.9267,
      908.0421, 908.7398, 911.8401, 914.5696, 912.0115, 919.4418,
      917.0436, 920.5495, 917.6138, 907.5037, 908.5145, 919.5846,
      917.6047, 926.8447, 910.6347, 912.8305, 907.7085, 911.6889,
  ]
  
  for n in xrange(1, 6):
      print 'Looking for %d maxima:' % n
      indexes = find_n_maxima(data, n)
      print indexes
      print ', '.join(str(data[i]) for i in indexes)
      print
  </code></pre>
  
  <p>Output:</p>
  
  <pre><code>Looking for 1 maxima:
  [78]
  2881.42
  
  Looking for 2 maxima:
  [78, 218]
  2881.42, 1927.083
  
  Looking for 3 maxima:
  [78, 108, 218]
  2881.42, 2147.193, 1927.083
  
  Looking for 4 maxima:
  [78, 108, 218, 274]
  2881.42, 2147.193, 1927.083, 936.7066
  
  Looking for 5 maxima:
  [78, 108, 218, 269, 274]
  2881.42, 2147.193, 1927.083, 939.6088, 936.7066
  </code></pre>