You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

378 lines
14 KiB

10 years ago
10 years ago
  1. # -*- coding: utf-8 -*-
  2. import array
  3. from datetime import datetime
  4. import numpy as np
  5. import re
  6. import struct
  7. import six
  8. class Nd2Parser(object):
  9. """
  10. Reads .nd2 files, provides an interface to the metadata, and generates numpy arrays from the image data.
  11. You should not ever need to instantiate this class manually unless you're a developer.
  12. """
  13. CHUNK_HEADER = 0xabeceda
  14. CHUNK_MAP_START = six.b("ND2 FILEMAP SIGNATURE NAME 0001!")
  15. CHUNK_MAP_END = six.b("ND2 CHUNK MAP SIGNATURE 0000001!")
  16. def __init__(self, filename):
  17. self._filename = filename
  18. self._fh = None
  19. self._chunk_map_start_location = None
  20. self._cursor_position = 0
  21. self._dimension_text = None
  22. self._label_map = {}
  23. self.metadata = {}
  24. self._read_map()
  25. self._parse_metadata()
  26. @property
  27. def _file_handle(self):
  28. if self._fh is None:
  29. self._fh = open(self._filename, "rb")
  30. return self._fh
  31. def _get_raw_image_data(self, image_group_number, channel_offset):
  32. """
  33. Reads the raw bytes and the timestamp of an image.
  34. :param image_group_number: groups are made of images with the same time index, field of view and z-level.
  35. :type image_group_number: int
  36. :param channel_offset: the offset in the array where the bytes for this image are found.
  37. :type channel_offset: int
  38. :return: (int, array.array()) or None
  39. """
  40. chunk = self._label_map[six.b("ImageDataSeq|%d!" % image_group_number)]
  41. data = self._read_chunk(chunk)
  42. # All images in the same image group share the same timestamp! So if you have complicated image data,
  43. # your timestamps may not be entirely accurate. Practically speaking though, they'll only be off by a few
  44. # seconds unless you're doing something super weird.
  45. timestamp = struct.unpack("d", data[:8])[0]
  46. image_group_data = array.array("H", data)
  47. image_data_start = 4 + channel_offset
  48. # The images for the various channels are interleaved within the same array. For example, the second image
  49. # of a four image group will be composed of bytes 2, 6, 10, etc. If you understand why someone would design
  50. # a data structure that way, please send the author of this library a message.
  51. image_data = image_group_data[image_data_start::self._channel_count]
  52. # Skip images that are all zeros! This is important, since NIS Elements creates blank "gap" images if you
  53. # don't have the same number of images each cycle. We discovered this because we only took GFP images every
  54. # other cycle to reduce phototoxicity, but NIS Elements still allocated memory as if we were going to take
  55. # them every cycle.
  56. if np.any(image_data):
  57. return timestamp, image_data
  58. return None
  59. @property
  60. def _dimensions(self):
  61. """
  62. While there are metadata values that represent a lot of what we want to capture, they seem to be unreliable.
  63. Sometimes certain elements don't exist, or change their data type randomly. However, the human-readable text
  64. is always there and in the same exact format, so we just parse that instead.
  65. :rtype: str
  66. """
  67. if self._dimension_text is None:
  68. for line in self.metadata[six.b('ImageTextInfo')][six.b('SLxImageTextInfo')].values():
  69. if six.b("Dimensions:") in line:
  70. metadata = line
  71. break
  72. else:
  73. raise ValueError("Could not parse metadata dimensions!")
  74. for line in metadata.split(six.b("\r\n")):
  75. if line.startswith(six.b("Dimensions:")):
  76. self._dimension_text = line
  77. break
  78. else:
  79. raise ValueError("Could not parse metadata dimensions!")
  80. if six.PY3:
  81. return self._dimension_text.decode("utf8")
  82. return self._dimension_text
  83. @property
  84. def _channels(self):
  85. """
  86. These are labels created by the NIS Elements user. Typically they may a short description of the filter cube
  87. used (e.g. "bright field", "GFP", etc.)
  88. :rtype: str
  89. """
  90. metadata = self.metadata[six.b('ImageMetadataSeq')][six.b('SLxPictureMetadata')][six.b('sPicturePlanes')]
  91. try:
  92. validity = self.metadata[six.b('ImageMetadata')][six.b('SLxExperiment')][six.b('ppNextLevelEx')][six.b('')][0][six.b('ppNextLevelEx')][six.b('')][0][six.b('pItemValid')]
  93. except KeyError:
  94. # If none of the channels have been deleted, there is no validity list, so we just make one
  95. validity = [True for _ in metadata]
  96. # Channel information is contained in dictionaries with the keys a0, a1...an where the number
  97. # indicates the order in which the channel is stored. So by sorting the dicts alphabetically
  98. # we get the correct order.
  99. for (label, chan), valid in zip(sorted(metadata[six.b('sPlaneNew')].items()), validity):
  100. if not valid:
  101. continue
  102. yield chan[six.b('sDescription')].decode("utf8")
  103. def _calculate_image_group_number(self, time_index, fov, z_level):
  104. """
  105. Images are grouped together if they share the same time index, field of view, and z-level.
  106. :type time_index: int
  107. :type fov: int
  108. :type z_level: int
  109. :rtype: int
  110. """
  111. return time_index * self._field_of_view_count * self._z_level_count + (fov * self._z_level_count + z_level)
  112. @property
  113. def _channel_offset(self):
  114. """
  115. Image data is interleaved for each image set. That is, if there are four images in a set, the first image
  116. will consist of pixels 1, 5, 9, etc, the second will be pixels 2, 6, 10, and so forth.
  117. :rtype: dict
  118. """
  119. channel_offset = {}
  120. for n, channel in enumerate(self._channels):
  121. channel_offset[channel] = n
  122. return channel_offset
  123. @property
  124. def _absolute_start(self):
  125. """
  126. The date and time when acquisition began.
  127. :rtype: datetime.datetime()
  128. """
  129. for line in self.metadata[six.b('ImageTextInfo')][six.b('SLxImageTextInfo')].values():
  130. line = line.decode("utf8")
  131. absolute_start_12 = None
  132. absolute_start_24 = None
  133. # ND2s seem to randomly switch between 12- and 24-hour representations.
  134. try:
  135. absolute_start_24 = datetime.strptime(line, "%m/%d/%Y %H:%M:%S")
  136. except (TypeError, ValueError):
  137. pass
  138. try:
  139. absolute_start_12 = datetime.strptime(line, "%m/%d/%Y %I:%M:%S %p")
  140. except (TypeError, ValueError):
  141. pass
  142. if not absolute_start_12 and not absolute_start_24:
  143. continue
  144. return absolute_start_12 if absolute_start_12 else absolute_start_24
  145. raise ValueError("This ND2 has no recorded start time. This is probably a bug.")
  146. @property
  147. def _channel_count(self):
  148. """
  149. The number of different channels used, including bright field.
  150. :rtype: int
  151. """
  152. pattern = r""".*?λ\((\d+)\).*?"""
  153. try:
  154. count = int(re.match(pattern, self._dimensions).group(1))
  155. except AttributeError as e:
  156. return 1
  157. else:
  158. return count
  159. @property
  160. def _field_of_view_count(self):
  161. """
  162. The metadata contains information about fields of view, but it contains it even if some fields
  163. of view were cropped. We can't find anything that states which fields of view are actually
  164. in the image data, so we have to calculate it. There probably is something somewhere, since
  165. NIS Elements can figure it out, but we haven't found it yet.
  166. :rtype: int
  167. """
  168. pattern = r""".*?XY\((\d+)\).*?"""
  169. try:
  170. count = int(re.match(pattern, self._dimensions).group(1))
  171. except AttributeError:
  172. return 1
  173. else:
  174. return count
  175. @property
  176. def _time_index_count(self):
  177. """
  178. The number of cycles.
  179. :rtype: int
  180. """
  181. pattern = r""".*?T'\((\d+)\).*?"""
  182. try:
  183. count = int(re.match(pattern, self._dimensions).group(1))
  184. except AttributeError:
  185. return 1
  186. else:
  187. return count
  188. @property
  189. def _z_level_count(self):
  190. """
  191. The number of different levels in the Z-plane.
  192. :rtype: int
  193. """
  194. pattern = r""".*?Z\((\d+)\).*?"""
  195. try:
  196. count = int(re.match(pattern, self._dimensions).group(1))
  197. except AttributeError:
  198. return 1
  199. else:
  200. return count
  201. @property
  202. def _image_count(self):
  203. """
  204. The total number of images in the ND2. Warning: this may be inaccurate as it includes "gap" images.
  205. :rtype: int
  206. """
  207. return self.metadata[six.b('ImageAttributes')][six.b('SLxImageAttributes')][six.b('uiSequenceCount')]
  208. def _parse_metadata(self):
  209. """
  210. Reads all metadata.
  211. """
  212. for label in self._label_map.keys():
  213. if label.endswith(six.b("LV!")) or six.b("LV|") in label:
  214. data = self._read_chunk(self._label_map[label])
  215. stop = label.index(six.b("LV"))
  216. self.metadata[label[:stop]] = self._read_metadata(data, 1)
  217. def _read_map(self):
  218. """
  219. Every label ends with an exclamation point, however, we can't directly search for those to find all the labels
  220. as some of the bytes contain the value 33, which is the ASCII code for "!". So we iteratively find each label,
  221. grab the subsequent data (always 16 bytes long), advance to the next label and repeat.
  222. """
  223. self._file_handle.seek(-8, 2)
  224. chunk_map_start_location = struct.unpack("Q", self._file_handle.read(8))[0]
  225. self._file_handle.seek(chunk_map_start_location)
  226. raw_text = self._file_handle.read(-1)
  227. label_start = raw_text.index(Nd2Parser.CHUNK_MAP_START) + 32
  228. while True:
  229. data_start = raw_text.index(six.b("!"), label_start) + 1
  230. key = raw_text[label_start: data_start]
  231. location, length = struct.unpack("QQ", raw_text[data_start: data_start + 16])
  232. if key == Nd2Parser.CHUNK_MAP_END:
  233. # We've reached the end of the chunk map
  234. break
  235. self._label_map[key] = location
  236. label_start = data_start + 16
  237. def _read_chunk(self, chunk_location):
  238. """
  239. Gets the data for a given chunk pointer
  240. """
  241. self._file_handle.seek(chunk_location)
  242. # The chunk metadata is always 16 bytes long
  243. chunk_metadata = self._file_handle.read(16)
  244. header, relative_offset, data_length = struct.unpack("IIQ", chunk_metadata)
  245. if header != Nd2Parser.CHUNK_HEADER:
  246. raise ValueError("The ND2 file seems to be corrupted.")
  247. # We start at the location of the chunk metadata, skip over the metadata, and then proceed to the
  248. # start of the actual data field, which is at some arbitrary place after the metadata.
  249. self._file_handle.seek(chunk_location + 16 + relative_offset)
  250. return self._file_handle.read(data_length)
  251. def _parse_unsigned_char(self, data):
  252. return struct.unpack("B", data.read(1))[0]
  253. def _parse_unsigned_int(self, data):
  254. return struct.unpack("I", data.read(4))[0]
  255. def _parse_unsigned_long(self, data):
  256. return struct.unpack("Q", data.read(8))[0]
  257. def _parse_double(self, data):
  258. return struct.unpack("d", data.read(8))[0]
  259. def _parse_string(self, data):
  260. value = data.read(2)
  261. while not value.endswith(six.b("\x00\x00")):
  262. # the string ends at the first instance of \x00\x00
  263. value += data.read(2)
  264. return value.decode("utf16")[:-1].encode("utf8")
  265. def _parse_char_array(self, data):
  266. array_length = struct.unpack("Q", data.read(8))[0]
  267. return array.array("B", data.read(array_length))
  268. def _parse_metadata_item(self, data):
  269. """
  270. Reads hierarchical data, analogous to a Python dict.
  271. """
  272. new_count, length = struct.unpack("<IQ", data.read(12))
  273. length -= data.tell() - self._cursor_position
  274. next_data_length = data.read(length)
  275. value = self._read_metadata(next_data_length, new_count)
  276. # Skip some offsets
  277. data.read(new_count * 8)
  278. return value
  279. def _get_value(self, data, data_type):
  280. """
  281. ND2s use various codes to indicate different data types, which we translate here.
  282. """
  283. parser = {1: self._parse_unsigned_char,
  284. 2: self._parse_unsigned_int,
  285. 3: self._parse_unsigned_int,
  286. 5: self._parse_unsigned_long,
  287. 6: self._parse_double,
  288. 8: self._parse_string,
  289. 9: self._parse_char_array,
  290. 11: self._parse_metadata_item}
  291. return parser[data_type](data)
  292. def _read_metadata(self, data, count):
  293. """
  294. Iterates over each element some section of the metadata and parses it.
  295. """
  296. data = six.BytesIO(data)
  297. metadata = {}
  298. for _ in range(count):
  299. self._cursor_position = data.tell()
  300. header = data.read(2)
  301. if not header:
  302. # We've reached the end of some hierarchy of data
  303. break
  304. if six.PY3:
  305. header = header.decode("utf8")
  306. data_type, name_length = map(ord, header)
  307. name = data.read(name_length * 2).decode("utf16")[:-1].encode("utf8")
  308. value = self._get_value(data, data_type)
  309. if name not in metadata.keys():
  310. metadata[name] = value
  311. else:
  312. if not isinstance(metadata[name], list):
  313. # We have encountered this key exactly once before. Since we're seeing it again, we know we
  314. # need to convert it to a list before proceeding.
  315. metadata[name] = [metadata[name]]
  316. # We've encountered this key before so we're guaranteed to be dealing with a list. Thus we append
  317. # the value to the already-existing list.
  318. metadata[name].append(value)
  319. return metadata