1
2
3
4
5
6 """
7
8 Authors: Henning O. Sorensen & Erik Knudsen
9 Center for Fundamental Research: Metal Structures in Four Dimensions
10 Risoe National Laboratory
11 Frederiksborgvej 399
12 DK-4000 Roskilde
13 email:erik.knudsen@risoe.dk
14
15 and Jon Wright, ESRF
16
17 """
18
19 import numpy as N, math, os, cStringIO, gzip, bz2
20 import Image
21 import fabioutils
22
23
24
25
26
27
29 """
30 A common object for images in fable
31 Contains a numpy array (.data) and dict of meta data (.header)
32 """
33
34 _need_a_seek_to_read = False
35 _need_a_real_file = False
36
37 - def __init__(self, data=None , header=None):
38 """
39 Set up initial values
40 """
41 if type(data) == type("string"):
42 raise Exception("fabioimage.__init__ bad argument - " + \
43 "data should be numpy array")
44 self.data = data
45 self.pilimage = None
46 if header is None:
47 self.header = {}
48 else:
49 self.header = header
50 self.header_keys = self.header.keys()
51 if data is not None:
52 self.dim1, self.dim2 = data.shape
53 else:
54 self.dim1 = self.dim2 = 0
55 self.bytecode = None
56 self.bpp = 2
57
58 self.mean = self.maxval = self.stddev = self.minval = None
59
60 self.area_sum = None
61 self.slice = None
62
63 self.nframes = 1
64 self.currentframe = 0
65 self.filename = None
66
68 """ returns the file numbered 'num' in the series as a fabioimage """
69 if self.nframes == 1:
70
71 import openimage
72 return openimage.openimage(
73 fabio.jump_filename(self.filename, num))
74 raise Exception("getframe out of range")
75
81
87
88
90 """
91 Convert to Python Imaging Library 16 bit greyscale image
92
93 FIXME - this should be handled by the libraries now
94 """
95 if filename:
96 self.read(filename)
97 if self.pilimage is not None:
98 return self.pilimage
99
100 size = self.data.shape[:2][::-1]
101 typmap = {
102 'float32' : "F" ,
103 'int32' : "F;32S" ,
104 'uint32' : "F;32" ,
105 'int16' : "F;16S" ,
106 'uint16' : "F;16" ,
107 'int8' : "F;8S" ,
108 'uint8' : "F;8" }
109 if typmap.has_key(self.data.dtype.name):
110 mode2 = typmap[ self.data.dtype.name ]
111 mode1 = mode2[0]
112 else:
113 raise Exception("Unknown numpy type " + str(self.data.dtype.type))
114
115
116 testval = N.array((1, 0), N.uint8).view(N.uint16)[0]
117 if testval == 1:
118 dats = self.data.tostring()
119 elif testval == 256:
120 dats = self.data.byteswap().tostring()
121 else:
122 raise Exception("Endian unknown in fabioimage.toPIL16")
123
124 self.pilimage = Image.frombuffer(mode1,
125 size,
126 dats,
127 "raw",
128 mode2,
129 0,
130 1)
131
132 return self.pilimage
133
135 """ returns self.header """
136 return self.header
137
139 """ Find max value in self.data, caching for the future """
140 if self.maxval is None:
141 self.maxval = N.max(self.data)
142 return self.maxval
143
145 """ Find min value in self.data, caching for the future """
146 if self.minval is None:
147 self.minval = N.min(self.data)
148 return self.minval
149
151 """
152 Convert a len(4) set of coords into a len(2)
153 tuple (pair) of slice objects
154 the latter are immutable, meaning the roi can be cached
155 """
156 assert len(coords) == 4
157 if len(coords) == 4:
158
159 if coords[0] > coords[2]:
160 coords[0:3:2] = [coords[2], coords[0]]
161 if coords[1] > coords[3]:
162 coords[1:4:2] = [coords[3], coords[1]]
163
164
165
166
167 fixme = (self.dim2 - coords[3] - 1,
168 coords[0] ,
169 self.dim2 - coords[1] - 1,
170 coords[2])
171 return (slice(int(fixme[0]), int(fixme[2]) + 1) ,
172 slice(int(fixme[1]), int(fixme[3]) + 1))
173
174
176 """
177 Sums up a region of interest
178 if len(coords) == 4 -> convert coords to slices
179 if len(coords) == 2 -> use as slices
180 floor -> ? removed as unused in the function.
181 """
182 if self.data == None:
183
184 return 0
185 if len(coords) == 4:
186 sli = self.make_slice(coords)
187 elif len(coords) == 2 and isinstance(coords[0], slice) and \
188 isinstance(coords[1], slice):
189 sli = coords
190 if sli == self.slice and self.area_sum is not None:
191 return self.area_sum
192 self.slice = sli
193 self.area_sum = N.sum(
194 N.ravel(
195 self.data[ self.slice ].astype(N.float)))
196 return self.area_sum
197
199 """ return the mean """
200 if self.mean is None:
201 self.mean = N.mean(self.data)
202 return float(self.mean)
203
205 """ return the standard deviation """
206 if self.stddev == None:
207 self.stddev = N.std(self.data)
208 return float(self.stddev)
209
210 - def add(self, other):
211 """
212 Add another Image - warnign, does not clip to 16 bit images by default
213 """
214 if not hasattr(other, 'data'):
215 print 'edfimage.add() called with something that ' + \
216 'does not have a data field'
217 assert self.data.shape == other.data.shape , \
218 'incompatible images - Do they have the same size?'
219 self.data = self.data + other.data
220 self.resetvals()
221
222
224 """ Reset cache - call on changing data """
225 self.mean = self.stddev = self.maxval = self.minval = None
226 self.area_sum = None
227
228 - def rebin(self, x_rebin_fact, y_rebin_fact):
229 """ Rebin the data and adjust dims """
230 if self.data == None:
231 raise Exception('Please read in the file you wish to rebin first')
232 (mantis_x, exp_x) = math.frexp(x_rebin_fact)
233 (mantis_y, exp_y) = math.frexp(y_rebin_fact)
234
235 if (mantis_x != 0.5 or mantis_y != 0.5):
236 raise Exception('Rebin factors not power of 2 not supported (yet)')
237 if int(self.dim1 / x_rebin_fact) * x_rebin_fact != self.dim1 or \
238 int(self.dim2 / x_rebin_fact) * x_rebin_fact != self.dim2 :
239 raise('image size is not divisible by rebin factor - ' + \
240 'skipping rebin')
241 pass
242 i = 1
243 while i < x_rebin_fact:
244
245 self.data = ((self.data[:, ::2] + self.data[:, 1::2]) / 2)
246 i = i * 2
247 i = 1
248 while i < y_rebin_fact:
249 self.data = ((self.data[::2, :] + self.data[1::2, :]) / 2)
250 i = i * 2
251 self.resetvals()
252 self.dim1 = self.dim1 / x_rebin_fact
253 self.dim2 = self.dim2 / y_rebin_fact
254
255 self.update_header()
256
258 """
259 To be overwritten - write the file
260 """
261 raise Exception("Class has not implemented readheader method yet")
262
274
276 """
277 Must be overridden in classes
278 """
279 raise Exception("Class has not implemented _readheader method yet")
280
282 """
283 update the header entries
284 by default pass in a dict of key, values.
285 """
286 self.header.update(kwds)
287
288 - def read(self, filename):
289 """
290 To be overridden - fill in self.header and self.data
291 """
292 raise Exception("Class has not implemented read method yet")
293
294
295 - def _open(self, fname, mode="rb"):
296 """
297 Try to handle compressed files, streams, shared memory etc
298 Return an object which can be used for "read" and "write"
299 ... FIXME - what about seek ?
300 """
301 self.filename = fname
302 if hasattr(fname, "read") and hasattr(fname, "write"):
303
304 return fname
305 if isinstance(fname, (str, unicode)):
306 self.header["filename"] = fname
307 if os.path.splitext(fname)[1] == ".gz":
308 return self._compressed_stream(fname,
309 fabioutils.COMPRESSORS['.gz'],
310 gzip.GzipFile,
311 mode)
312 if os.path.splitext(fname)[1] == '.bz2':
313 return self._compressed_stream(fname,
314 fabioutils.COMPRESSORS['.bz2'],
315 bz2.BZ2File,
316 mode)
317
318
319
320
321
322 return open(fname, mode)
323
324 - def _compressed_stream(self,
325 fname,
326 system_uncompress,
327 python_uncompress,
328 mode='rb'):
329 """
330 Try to transparently handle gzip / bzip without always getting python
331 performance
332 """
333
334
335 if self._need_a_real_file and mode[0] == "r":
336 fo = python_uncompress(fname, mode)
337 fobj = os.tmpfile()
338 fobj.write(fo.read())
339 fo.close()
340 fobj.seek(0)
341 return fobj
342 if self._need_a_seek_to_read and mode[0] == "r":
343 fo = python_uncompress(fname, mode)
344 return cStringIO.StringIO(fo.read())
345 return python_uncompress(fname, mode)
346
347
348
349
351 """
352 check some basic fabioimage functionality
353 """
354 import time
355 start = time.time()
356
357 dat = N.ones((1024, 1024), N.uint16)
358 dat = (dat * 50000).astype(N.uint16)
359 assert dat.dtype.char == N.ones((1), N.uint16).dtype.char
360 hed = {"Title":"50000 everywhere"}
361 obj = fabioimage(dat, hed)
362
363 assert obj.getmax() == 50000
364 assert obj.getmin() == 50000
365 assert obj.getmean() == 50000 , obj.getmean()
366 assert obj.getstddev() == 0.
367
368 dat2 = N.zeros((1024, 1024), N.uint16, savespace=1)
369 cord = [ 256, 256, 790, 768 ]
370 slic = obj.make_slice(cord)
371 dat2[slic] = dat2[slic] + 100
372
373 obj = fabioimage(dat2, hed)
374
375
376 assert obj.maxval is None
377 assert obj.minval is None
378
379 assert obj.getmax() == 100, obj.getmax()
380 assert obj.getmin() == 0 , obj.getmin()
381 npix = (slic[0].stop - slic[0].start) * (slic[1].stop - slic[1].start)
382 obj.resetvals()
383 area1 = obj.integrate_area(cord)
384 obj.resetvals()
385 area2 = obj.integrate_area(slic)
386 assert area1 == area2
387 assert obj.integrate_area(cord) == obj.integrate_area(slic)
388 assert obj.integrate_area(cord) == npix * 100, obj.integrate_area(cord)
389
390
391 def clean():
392 """ clean up the created testfiles"""
393 for name in ["testfile", "testfile.gz", "testfile.bz2"]:
394 try:
395 os.remove(name)
396 except:
397 continue
398
399
400 clean()
401
402 gzip.open("testfile.gz", "wb").write("{ hello }")
403 fout = obj._open("testfile.gz")
404 readin = fout.read()
405 assert readin == "{ hello }", readin + " gzipped file"
406
407
408 bz2.BZ2File("testfilebz", "wb").write("{ hello }")
409 fout = obj._open("testfile.bz2")
410 readin = fout.read()
411 assert readin == "{ hello }", readin + " bzipped file"
412
413 ftest = open("testfile", "wb")
414 ftest.write("{ hello }")
415 assert ftest == obj._open(ftest)
416 ftest.close()
417 fout = obj._open("testfile")
418 readin = fout.read()
419 assert readin == "{ hello }", readin + "plain file"
420 fout.close()
421 ftest.close()
422 clean()
423
424 print "Passed in", time.time() - start, "s"
425
426 if __name__ == '__main__':
427 test()
428