Instrument Neutral Distributed Interface INDI  1.9.5
file.c
Go to the documentation of this file.
1 /*
2  * libDSPAU - a digital signal processing library for astronoms usage
3  * Copyright (C) 2017 Ilia Platone <info@iliaplatone.com>
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  */
18 
19 #include "dsp.h"
20 #include <fitsio.h>
21 #include <time.h>
22 #include <locale.h>
23 #include <unistd.h>
24 #include <jpeglib.h>
25 
26 dsp_stream_p* dsp_file_read_fits(char *filename, int *channels, int stretch)
27 {
28  fitsfile *fptr = (fitsfile*)malloc(sizeof(fitsfile));
29  memset(fptr, 0, sizeof(fitsfile));
30  int bpp = 16;
31  int status = 0;
32  char value[150];
33  char comment[150];
34  char error_status[64];
35 
36  fits_open_file(&fptr, filename, READONLY_FILE, &status);
37 
38  if (status)
39  {
40  goto fail;
41  }
42 
43  fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status);
44  if(status)
45  {
46  goto fail;
47  }
48 
49  int dims;
50  long naxes[3] = { 1, 1, 1 };
51  fits_get_img_param(fptr, 3, &bpp, &dims, naxes, &status);
52  if(status)
53  {
54  goto fail;
55  }
56 
57  int dim, nelements = 1;
58  for(dim = 0; dim < dims; dim++) {
59  nelements *= naxes[dim];
60  }
61  void *array = malloc((unsigned int)(abs(bpp) * nelements / 8));
62  int anynul = 0;
63  dsp_t* buf = (dsp_t*)malloc(sizeof(dsp_t)*(unsigned int)nelements);
64  switch(bpp) {
65  case BYTE_IMG:
66  fits_read_img(fptr, TBYTE, 1, (long)nelements, NULL, array, &anynul, &status);
67  dsp_buffer_stretch(((unsigned char*)(array)), (long)nelements, 0, dsp_t_max);
68  dsp_buffer_copy(((unsigned char*)array), buf, nelements);
69  break;
70  case SHORT_IMG:
71  fits_read_img(fptr, TUSHORT, 1, (long)nelements, NULL, array, &anynul, &status);
72  if(abs(bpp) > 8*(int)sizeof(dsp_t))
73  dsp_buffer_stretch(((short*)(array)), (long)nelements, 0, dsp_t_max);
74  dsp_buffer_copy(((short*)array), buf, nelements);
75  break;
76  case USHORT_IMG:
77  fits_read_img(fptr, TUSHORT, 1, (long)nelements, NULL, array, &anynul, &status);
78  if(abs(bpp) > 8*(int)sizeof(dsp_t))
79  dsp_buffer_stretch(((unsigned short*)(array)), (long)nelements, 0, dsp_t_max);
80  dsp_buffer_copy(((unsigned short*)array), buf, nelements);
81  break;
82  case LONG_IMG:
83  fits_read_img(fptr, TULONG, 1, (long)nelements, NULL, array, &anynul, &status);
84  if(abs(bpp) > 8*(int)sizeof(dsp_t))
85  dsp_buffer_stretch(((int*)(array)), (long)nelements, 0, dsp_t_max);
86  dsp_buffer_copy(((int*)array), buf, nelements);
87  break;
88  case ULONG_IMG:
89  fits_read_img(fptr, TULONG, 1, (long)nelements, NULL, array, &anynul, &status);
90  if(abs(bpp) > 8*(int)sizeof(dsp_t))
91  dsp_buffer_stretch(((unsigned int*)(array)), (long)nelements, 0, dsp_t_max);
92  dsp_buffer_copy(((unsigned int*)array), buf, nelements);
93  break;
94  case LONGLONG_IMG:
95  fits_read_img(fptr, TLONGLONG, 1, (long)nelements, NULL, array, &anynul, &status);
96  if(abs(bpp) > 8*(int)sizeof(dsp_t))
97  dsp_buffer_stretch(((long*)(array)), (long)nelements, 0, dsp_t_max);
98  dsp_buffer_copy(((long*)array), buf, nelements);
99  break;
100  case FLOAT_IMG:
101  fits_read_img(fptr, TFLOAT, 1, (long)nelements, NULL, array, &anynul, &status);
102  if(abs(bpp) > 8*(int)sizeof(dsp_t))
103  dsp_buffer_stretch(((float*)(array)), (long)nelements, 0, dsp_t_max);
104  dsp_buffer_copy(((float*)array), buf, nelements);
105  break;
106  case DOUBLE_IMG:
107  fits_read_img(fptr, TDOUBLE, 1, (long)nelements, NULL, array, &anynul, &status);
108  if(abs(bpp) > 8*(int)sizeof(dsp_t))
109  dsp_buffer_stretch(((double*)(array)), (long)nelements, 0, dsp_t_max);
110  dsp_buffer_copy(((double*)array), buf, nelements);
111  break;
112  }
113  free(array);
114  if(status||anynul)
115  goto fail;
116  if(stretch)
117  dsp_buffer_stretch(buf, nelements, 0, dsp_t_max);
118 
119  int red = -1;
120  ffgkey(fptr, "XBAYROFF", value, comment, &status);
121  if (!status)
122  {
123  red = atoi(value);
124  }
125 
126  ffgkey(fptr, "YBAYROFF", value, comment, &status);
127 
128  if (!status)
129  {
130  red |= atoi(value)<<1;
131  }
132  fits_close_file(fptr, &status);
133  if(status)
134  goto fail;
135  int sizes[3] = { (int)naxes[0], (int)naxes[1], (int)naxes[2] };
136  *channels = (int)naxes[2];
137  return dsp_stream_from_components(buf, dims, sizes, sizes[2]);
138 fail:
139  if(status) {
140  fits_report_error(stderr, status); /* print out any error messages */
141  fits_get_errstatus(status, error_status);
142  fprintf(stderr, "FITS Error: %s\n", error_status);
143  }
144  return NULL;
145 }
146 
147 void* dsp_file_write_fits(int bpp, size_t* memsize, dsp_stream_p stream)
148 {
149  dsp_stream_p tmp = dsp_stream_copy(stream);
150  int img_type = USHORT_IMG;
151  int byte_type = TUSHORT;
152  char bit_depth[64] = "16 bits per sample";
153  void* buf = malloc((unsigned int)(tmp->len * abs(bpp) / 8 + 512));
154  void *memfile = NULL;
155  int status = 0;
156  int naxis = tmp->dims;
157  long *naxes = (long*)malloc(sizeof(long) * (unsigned int)tmp->dims);
158  long nelements = tmp->len;
159  char error_status[64];
160  unsigned int i;
161 
162  for (i = 0; i < (unsigned int)tmp->dims; i++)
163  naxes[i] = tmp->sizes[i];
164  switch (bpp)
165  {
166  case 8:
167  byte_type = TBYTE;
168  img_type = BYTE_IMG;
169  dsp_buffer_copy(tmp->buf, ((unsigned char*)buf), tmp->len);
170  strcpy(bit_depth, "8 bits unsigned integer per sample");
171  break;
172 
173  case 16:
174  byte_type = TUSHORT;
175  img_type = USHORT_IMG;
176  dsp_buffer_copy(tmp->buf, ((unsigned short*)buf), tmp->len);
177  strcpy(bit_depth, "16 bits unsigned integer per sample");
178  break;
179 
180  case 32:
181  byte_type = TULONG;
182  img_type = ULONG_IMG;
183  dsp_buffer_copy(tmp->buf, ((unsigned int*)buf), tmp->len);
184  strcpy(bit_depth, "32 bits unsigned integer per sample");
185  break;
186 
187  case 64:
188  byte_type = TLONGLONG;
189  img_type = LONGLONG_IMG;
190  dsp_buffer_copy(tmp->buf, ((unsigned long*)buf), tmp->len);
191  strcpy(bit_depth, "64 bits unsigned integer per sample");
192  break;
193 
194  case -32:
195  byte_type = TFLOAT;
196  img_type = FLOAT_IMG;
197  dsp_buffer_copy(tmp->buf, ((float*)buf), tmp->len);
198  strcpy(bit_depth, "32 bits floating point per sample");
199  break;
200 
201  case -64:
202  byte_type = TDOUBLE;
203  img_type = DOUBLE_IMG;
204  dsp_buffer_copy(tmp->buf, ((double*)buf), tmp->len);
205  strcpy(bit_depth, "64 bits floating point per sample");
206  break;
207 
208  default:
209  fprintf(stderr, "Unsupported bits per sample value %d", bpp);
210  goto fail;
211  }
212 
213  fitsfile *fptr = NULL;
214  size_t len = 5760;
215  memfile = malloc(len);
216  fits_create_memfile(&fptr, &memfile, &len, 2880, realloc, &status);
217 
218  if (status)
219  {
220  goto fail_fptr;
221  }
222 
223  fits_create_img(fptr, img_type, naxis, naxes, &status);
224 
225  if (status)
226  {
227  goto fail_fptr;
228  }
229 
230  fits_write_img(fptr, byte_type, 1, nelements, buf, &status);
231 
232  if (status)
233  {
234  goto fail_fptr;
235  }
236 
237  fits_close_file(fptr, &status);
238 
239  *memsize = len;
240 
241 fail_fptr:
242  if(status) {
243  fits_report_error(stderr, status);
244  fits_get_errstatus(status, error_status);
245  fprintf(stderr, "FITS Error: %s\n", error_status);
246  free(memfile);
247  }
248 fail:
250  dsp_stream_free(tmp);
251  free(naxes);
252  free(buf);
253  return memfile;
254 }
255 
256 dsp_stream_p* dsp_file_read_jpeg(char *filename, int *channels, int stretch)
257 {
258  int width, height;
259  unsigned int components;
260  unsigned int bpp = 8;
261  unsigned char * buf;
262  struct jpeg_decompress_struct info;
263  struct jpeg_error_mgr err;
264 
265  info.err = jpeg_std_error(& err);
266  jpeg_create_decompress(& info);
267 
268  FILE *jpeg = fopen (filename, "r");
269  if(jpeg == NULL)
270  return NULL;
271 
272  jpeg_stdio_src(&info, jpeg);
273  jpeg_read_header(&info, TRUE);
274  info.dct_method = JDCT_FLOAT;
275  jpeg_start_decompress(&info);
276  width = (int)info.output_width;
277  height = (int)info.output_height;
278  components = (unsigned int)info.num_components;
279 
280  int row_stride = (int)components * (int)width;
281  buf = (unsigned char *)malloc((unsigned int)(width * height) * components);
282  unsigned char *image = buf;
283  unsigned int row;
284  for (row = 0; row < (unsigned int)height; row++)
285  {
286  jpeg_read_scanlines(&info, &image, 1);
287  image += row_stride;
288  }
289  jpeg_finish_decompress(&info);
290  *channels = (int)components;
291  return dsp_buffer_rgb_to_components(buf, 2, (int[]){width, height}, (int)components, (int)bpp, stretch);
292 }
293 
294 void dsp_file_write_jpeg_composite(char *filename, int components, int quality, dsp_stream_p* stream)
295 {
296  int bpp = 8;
297  unsigned int row_stride;
298  unsigned int width = (unsigned int)stream[0]->sizes[0];
299  unsigned int height = (unsigned int)stream[0]->sizes[1];
300  void *buf = malloc((unsigned int)(stream[0]->len*components*bpp/8));
301  unsigned char *image = (unsigned char *)buf;
302  dsp_buffer_components_to_rgb(stream, buf, components, bpp);
303 
304  struct jpeg_compress_struct cinfo;
305  struct jpeg_error_mgr jerr;
306  FILE * outfile;
307  cinfo.err = jpeg_std_error(&jerr);
308  jpeg_create_compress(&cinfo);
309  if ((outfile = fopen(filename, "wb")) == NULL) {
310  fprintf(stderr, "can't open %s\n", filename);
311  return;
312  }
313  jpeg_stdio_dest(&cinfo, outfile);
314  cinfo.image_width = width;
315  cinfo.image_height = height;
316  cinfo.in_color_space = (components == 1 ? JCS_GRAYSCALE : JCS_RGB);
317  cinfo.input_components = components;
318  jpeg_set_defaults(&cinfo);
319  cinfo.dct_method = JDCT_FLOAT;
320  cinfo.optimize_coding = TRUE;
321  cinfo.restart_in_rows = 1;
322 
323  jpeg_set_quality(&cinfo, quality, TRUE);
324  jpeg_start_compress(&cinfo, TRUE);
325  row_stride = width * (unsigned int)components;
326  int row;
327  for (row = 0; row < (int)height; row++)
328  {
329  jpeg_write_scanlines(&cinfo, &image, 1);
330  image += row_stride;
331  }
332 
333  jpeg_finish_compress(&cinfo);
334  jpeg_destroy_compress(&cinfo);
335  fclose(outfile);
336  free(buf);
337 }
338 
339 dsp_stream_p *dsp_stream_from_components(dsp_t* buf, int dims, int *sizes, int components)
340 {
341  int d, y, x;
342  dsp_stream_p* picture = (dsp_stream_p*) malloc(sizeof(dsp_stream_p)*(unsigned int)(components+1));
343  for(y = 0; y <= components; y++) {
344  picture[y] = dsp_stream_new();
345  for(d = 0; d < dims; d++)
346  dsp_stream_add_dim(picture[y], sizes[d]);
347  dsp_stream_alloc_buffer(picture[y], picture[y]->len);
348  if(y < components) {
349  dsp_buffer_copy(((dsp_t*)&buf[y*picture[y]->len]), picture[y]->buf, picture[y]->len);
350  } else {
351  for(x = 0; x < picture[y]->len; x++) {
352  double val = 0;
353  for(d = 0; d < components; d++) {
354  val += buf[x+d*picture[y]->len];
355  }
356  picture[y]->buf[x] = (dsp_t)val/components;
357  }
358  }
359  }
360  free (buf);
361  return picture;
362 }
363 
364 dsp_stream_p *dsp_buffer_rgb_to_components(void* buf, int dims, int *sizes, int components, int bpp, int stretch)
365 {
366  dsp_stream_p* picture = (dsp_stream_p*) malloc(sizeof(dsp_stream_p)*(unsigned int)(components+1));
367  int x, y, z, d;
368  for(y = 0; y < components; y++) {
369  picture[y] = dsp_stream_new();
370  for(d = 0; d < dims; d++)
371  dsp_stream_add_dim(picture[y], sizes[d]);
372  dsp_stream_alloc_buffer(picture[y], picture[y]->len);
373  switch(bpp) {
374  case 8:
375  dsp_buffer_copy_stepping(((unsigned char*)&((unsigned char*)buf)[y]), picture[y]->buf, picture[y]->len*components, picture[y]->len, components, 1);
376  break;
377  case 16:
378  dsp_buffer_copy_stepping(((unsigned short*)&((unsigned short*)buf)[y]), picture[y]->buf, picture[y]->len*components, picture[y]->len, components, 1);
379  break;
380  case 32:
381  dsp_buffer_copy_stepping(((unsigned int*)&((unsigned int*)buf)[y]), picture[y]->buf, picture[y]->len*components, picture[y]->len, components, 1);
382  break;
383  case 64:
384  dsp_buffer_copy_stepping(((unsigned long*)&((unsigned long*)buf)[y]), picture[y]->buf, picture[y]->len*components, picture[y]->len, components, 1);
385  break;
386  case -32:
387  dsp_buffer_copy_stepping(((float*)&((float*)buf)[y]), picture[y]->buf, picture[y]->len*components, picture[y]->len, components, 1);
388  break;
389  case -64:
390  dsp_buffer_copy_stepping(((double*)&((double*)buf)[y]), picture[y]->buf, picture[y]->len*components, picture[y]->len, components, 1);
391  break;
392  default:
393  break;
394  }
395  if(stretch) {
396  dsp_buffer_stretch(picture[y]->buf, picture[y]->len, 0, dsp_t_max);
397  }
398  }
399  picture[y] = dsp_stream_new();
400  for(d = 0; d < dims; d++)
401  dsp_stream_add_dim(picture[y], sizes[d]);
402  dsp_stream_alloc_buffer(picture[y], picture[y]->len);
403  for(x = 0; x < picture[y]->len; x++) {
404  double v = 0;
405  switch(bpp) {
406  case 8:
407  for(z = 0; z < y; z++) v += ((unsigned char*)buf)[x*y+z];
408  break;
409  case 16:
410  for(z = 0; z < y; z++) v += ((unsigned short*)buf)[x*y+z];
411  break;
412  case 32:
413  for(z = 0; z < y; z++) v += ((unsigned int*)buf)[x*y+z];
414  break;
415  case 64:
416  for(z = 0; z < y; z++) v += ((unsigned long*)buf)[x*y+z];
417  break;
418  case -32:
419  for(z = 0; z < y; z++) v += (double)((float*)buf)[x*y+z];
420  break;
421  case -64:
422  for(z = 0; z < y; z++) v += ((double*)buf)[x*y+z];
423  break;
424  default:
425  break;
426  }
427  picture[y]->buf[x] = (dsp_t)(v/y);
428  }
429  if(stretch) {
430  dsp_buffer_stretch(picture[y]->buf, picture[y]->len, 0, dsp_t_max);
431  }
432  free (buf);
433  return picture;
434 }
435 
436 void dsp_buffer_components_to_rgb(dsp_stream_p *stream, void* rgb, int components, int bpp)
437 {
438  size_t y;
439  int len = stream[0]->len * components;
440  dsp_t max = (dsp_t)((double)((1<<abs(bpp))-1));
441  max = Min(max, (dsp_t)~0);
442  for(y = 0; y < (size_t)components; y++) {
443  dsp_stream_p in = dsp_stream_copy(stream[y]);
444  dsp_buffer_stretch(in->buf, in->len, 0, max);
445  switch(bpp) {
446  case 8:
447  dsp_buffer_copy_stepping(in->buf, ((unsigned char*)&((unsigned char*)rgb)[y]), in->len, len, 1, components);
448  break;
449  case 16:
450  dsp_buffer_copy_stepping(in->buf, ((unsigned short*)&((unsigned short*)rgb)[y]), in->len, len, 1, components);
451  break;
452  case 32:
453  dsp_buffer_copy_stepping(in->buf, ((unsigned int*)&((unsigned int*)rgb)[y]), in->len, len, 1, components);
454  break;
455  case 64:
456  dsp_buffer_copy_stepping(in->buf, ((unsigned long*)&((unsigned long*)rgb)[y]), in->len, len, 1, components);
457  break;
458  case -32:
459  dsp_buffer_copy_stepping(in->buf, ((float*)&((float*)rgb)[y]), in->len, len, 1, components);
460  break;
461  case -64:
462  dsp_buffer_copy_stepping(in->buf, ((double*)&((double*)rgb)[y]), in->len, len, 1, components);
463  break;
464  }
466  dsp_stream_free(in);
467  }
468 }
dsp_file_write_jpeg_composite
void dsp_file_write_jpeg_composite(char *filename, int components, int quality, dsp_stream_p *stream)
Write the components dsp_stream_p array into a JPEG file,.
Definition: file.c:294
dsp.h
dsp_stream_new
DLL_EXPORT dsp_stream_p dsp_stream_new()
Allocate a new DSP stream type.
Definition: stream.c:48
dsp_file_read_jpeg
dsp_stream_p * dsp_file_read_jpeg(char *filename, int *channels, int stretch)
Read a JPEG file and fill a array of dsp_stream_p with its content, each color channel has its own st...
Definition: file.c:256
dsp_t_max
#define dsp_t_max
Definition: dsp.h:58
dsp_stream_from_components
dsp_stream_p * dsp_stream_from_components(dsp_t *buf, int dims, int *sizes, int components)
Definition: file.c:339
dsp_stream_free
DLL_EXPORT void dsp_stream_free(dsp_stream_p stream)
Free the DSP stream passed as argument.
Definition: stream.c:76
max
double max(void)
dsp_stream_t
Contains a set of informations and data relative to a buffer and how to use it.
Definition: dsp.h:175
dsp_stream_copy
DLL_EXPORT dsp_stream_p dsp_stream_copy(dsp_stream_p stream)
Create a copy of the DSP stream passed as argument.
Definition: stream.c:89
TRUE
#define TRUE
Definition: stvdriver.c:54
Min
#define Min(a, b)
if min() is not present you can use this one
Definition: dsp.h:65
dsp_file_read_fits
dsp_stream_p * dsp_file_read_fits(char *filename, int *channels, int stretch)
Read a FITS file and fill a dsp_stream_p with its content.
Definition: file.c:26
dsp_buffer_rgb_to_components
dsp_stream_p * dsp_buffer_rgb_to_components(void *buf, int dims, int *sizes, int components, int bpp, int stretch)
Definition: file.c:364
dsp_stream_t::len
int len
The buffers length.
Definition: dsp.h:178
dsp_buffer_copy
#define dsp_buffer_copy(in, out, len)
Fill the output buffer with the values of the elements of the input stream by casting them to the out...
Definition: dsp.h:791
dsp_stream_free_buffer
DLL_EXPORT void dsp_stream_free_buffer(dsp_stream_p stream)
Free the buffer of the DSP Stream passed as argument.
Definition: stream.c:41
dsp_stream_t::sizes
int * sizes
Sizes of each dimension.
Definition: dsp.h:182
dsp_buffer_copy_stepping
#define dsp_buffer_copy_stepping(in, out, inlen, outlen, instep, outstep)
Fill the output buffer with the values of the elements of the input stream by casting them to the out...
Definition: dsp.h:811
dsp_file_write_fits
void * dsp_file_write_fits(int bpp, size_t *memsize, dsp_stream_p stream)
Write the components dsp_stream_p array into a FITS file,.
Definition: file.c:147
dsp_stream_alloc_buffer
DLL_EXPORT void dsp_stream_alloc_buffer(dsp_stream_p stream, int len)
Allocate a buffer with length len on the stream passed as argument.
Definition: stream.c:21
dsp_buffer_components_to_rgb
void dsp_buffer_components_to_rgb(dsp_stream_p *stream, void *rgb, int components, int bpp)
Definition: file.c:436
dsp_stream_t::dims
int dims
Number of dimensions of the buffers.
Definition: dsp.h:180
dsp_stream_add_dim
DLL_EXPORT void dsp_stream_add_dim(dsp_stream_p stream, int len)
Add a dimension with length len to a DSP stream.
Definition: stream.c:111
dsp_t
#define dsp_t
Definition: dsp.h:57
dsp_buffer_stretch
#define dsp_buffer_stretch(buf, len, _mn, _mx)
Stretch minimum and maximum values of the input stream.
Definition: dsp.h:557
dsp_stream_t::buf
dsp_t * buf
buffer
Definition: dsp.h:184