Instrument Neutral Distributed Interface INDI  2.0.2
file.c
Go to the documentation of this file.
1 /*
2 * DSP API - a digital signal processing library for astronomy usage
3 * Copyright © 2017-2022 Ilia Platone
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 3 of the License, or (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 GNU
13 * Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with this program; if not, write to the Free Software Foundation,
17 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19 
20 #include "dsp.h"
21 #include <fitsio.h>
22 #include <time.h>
23 #include <locale.h>
24 #include <unistd.h>
25 #include <jpeglib.h>
26 
27 dsp_stream_p* dsp_file_read_fits(const char* filename, int *channels, int stretch)
28 {
29  fitsfile *fptr;
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, &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  bpp = fmax(16, bpp);
57  int dim, nelements = 1;
58  for(dim = 0; dim < dims; dim++) {
59  nelements *= naxes[dim];
60  }
61  void *array = malloc(((size_t)abs(bpp) * nelements / 8));
62  int anynul = 0;
63  dsp_t* buf = (dsp_t*)malloc(sizeof(dsp_t)*(size_t)(nelements+1));
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((size_t)abs(bpp) > 8*sizeof(dsp_t))
73  dsp_buffer_stretch(((unsigned short*)(array)), (long)nelements, 0, dsp_t_max);
74  dsp_buffer_copy(((unsigned 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((size_t)abs(bpp) > 8*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((size_t)abs(bpp) > 8*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((size_t)abs(bpp) > 8*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((size_t)abs(bpp) > 8*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((size_t)abs(bpp) > 8*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((size_t)abs(bpp) > 8*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  int red = -1;
117  ffgkey(fptr, "XBAYROFF", value, comment, &status);
118  if (!status)
119  {
120  red = atoi(value);
121  }
122  status = 0;
123  ffgkey(fptr, "YBAYROFF", value, comment, &status);
124 
125  if (!status)
126  {
127  red |= atoi(value)<<1;
128  }
129  status = 0;
130  fits_close_file(fptr, &status);
131  if(status)
132  goto fail;
133  dsp_stream_p* stream;
134  int x;
135  void* rgb = NULL;
136  int sizes[2] = { naxes[0], naxes[1] };
137  if(red > -1) {
138  *channels = 3;
139  rgb = dsp_file_bayer_2_rgb(buf, red, naxes[0], naxes[1]);
140  stream = dsp_buffer_rgb_to_components(rgb, dims, sizes, *channels, 32, 0);
141  } else {
142  *channels = naxes[2];
143  dims = 2;
144  stream = dsp_stream_from_components(buf, dims, sizes, *channels);
145  }
146  free(buf);
147  if(stream) {
148  if(stretch) {
149  for (x = 0; x < *channels; x++) {
150  dsp_buffer_pow1(stream[x], 0.5);
151  dsp_buffer_stretch(stream[x]->buf, stream[x]->len, 0, dsp_t_max);
152  }
153  }
154  return stream;
155  }
156 fail:
157  if(status) {
158  fits_get_errstatus(status, error_status);
159  perr("FITS Error: %s\n", error_status);
160  }
161  return NULL;
162 }
163 
164 void dsp_file_write_fits(const char* filename, int bpp, dsp_stream_p stream)
165 {
166  dsp_stream_p tmp = dsp_stream_copy(stream);
167  int img_type = USHORT_IMG;
168  int byte_type = TUSHORT;
169  char bit_depth[64] = "16 bits per sample";
170  void* buf = malloc((size_t)(tmp->len * (size_t)abs(bpp) / 8 + 512));
171  int status = 0;
172  int naxis = tmp->dims;
173  long *naxes = (long*)malloc(sizeof(long) * (size_t)tmp->dims);
174  long nelements = tmp->len;
175  char error_status[64];
176  int i;
177  dsp_buffer_stretch(tmp->buf, tmp->len, 0, dsp_t_max);
178  for (i = 0; i < tmp->dims; i++)
179  naxes[i] = tmp->sizes[i];
180  switch (bpp)
181  {
182  case 8:
183  byte_type = TBYTE;
184  img_type = BYTE_IMG;
185  dsp_buffer_copy(tmp->buf, ((unsigned char*)buf), tmp->len);
186  strcpy(bit_depth, "8 bits unsigned integer per sample");
187  break;
188 
189  case 16:
190  byte_type = TUSHORT;
191  img_type = USHORT_IMG;
192  dsp_buffer_copy(tmp->buf, ((unsigned short*)buf), tmp->len);
193  strcpy(bit_depth, "16 bits unsigned integer per sample");
194  break;
195 
196  case 32:
197  byte_type = TULONG;
198  img_type = ULONG_IMG;
199  dsp_buffer_copy(tmp->buf, ((unsigned int*)buf), tmp->len);
200  strcpy(bit_depth, "32 bits unsigned integer per sample");
201  break;
202 
203  case 64:
204  byte_type = TLONGLONG;
205  img_type = LONGLONG_IMG;
206  dsp_buffer_copy(tmp->buf, ((unsigned long*)buf), tmp->len);
207  strcpy(bit_depth, "64 bits unsigned integer per sample");
208  break;
209 
210  case -32:
211  byte_type = TFLOAT;
212  img_type = FLOAT_IMG;
213  dsp_buffer_copy(tmp->buf, ((float*)buf), tmp->len);
214  strcpy(bit_depth, "32 bits floating point per sample");
215  break;
216 
217  case -64:
218  byte_type = TDOUBLE;
219  img_type = DOUBLE_IMG;
220  dsp_buffer_copy(tmp->buf, ((double*)buf), tmp->len);
221  strcpy(bit_depth, "64 bits floating point per sample");
222  break;
223 
224  default:
225  perr("Unsupported bits per sample value %d", bpp);
226  goto fail;
227  }
228 
229  unlink(filename);
230  fitsfile *fptr;
231  fits_create_file(&fptr, filename, &status);
232 
233  if (status)
234  {
235  goto fail_fptr;
236  }
237 
238  fits_create_img(fptr, img_type, naxis, naxes, &status);
239 
240  if (status)
241  {
242  goto fail_fptr;
243  }
244 
245  fits_write_img(fptr, byte_type, 1, nelements, buf, &status);
246 
247  if (status)
248  {
249  goto fail_fptr;
250  }
251 
252  fits_close_file(fptr, &status);
253 
254 fail_fptr:
255  if(status) {
256  fits_get_errstatus(status, error_status);
257  perr("FITS Error: %s\n", error_status);
258  }
259 fail:
261  dsp_stream_free(tmp);
262  free(naxes);
263  free (buf);
264 }
265 
266 void dsp_file_write_fits_composite(const char* filename, int components, int bpp, dsp_stream_p* stream)
267 {
268  int x;
269  dsp_stream_p tmp = stream[components];
270  int img_type = USHORT_IMG;
271  int byte_type = TUSHORT;
272  char bit_depth[64] = "16 bits per sample";
273  void* buf = malloc((size_t)(tmp->len * components * (size_t)abs(bpp) / 8 + 512));
274  int status = 0;
275  int naxis = tmp->dims + 1;
276  long *naxes = (long*)malloc(sizeof(long) * (size_t)(tmp->dims + 1));
277  long nelements = tmp->len * components;
278  char error_status[64];
279  int i;
280  for (i = 0; i < tmp->dims; i++)
281  naxes[i] = tmp->sizes[i];
282  naxes[i] = components;
283  dsp_t max = (1<<(size_t)abs(bpp))/2-1;
284  for(x = 0; x < components; x++) {
285  tmp = dsp_stream_copy(stream[x]);
286  dsp_buffer_stretch(tmp->buf, tmp->len, 0, (dsp_t)max);
287  switch (bpp)
288  {
289  case 8:
290  byte_type = TBYTE;
291  img_type = BYTE_IMG;
292  dsp_buffer_copy(tmp->buf, ((unsigned char*)&((unsigned char*)buf)[tmp->len*x]), tmp->len);
293  strcpy(bit_depth, "8 bits unsigned integer per sample");
294  break;
295 
296  case 16:
297  byte_type = TUSHORT;
298  img_type = USHORT_IMG;
299  dsp_buffer_copy(tmp->buf, ((unsigned short*)&((unsigned short*)buf)[tmp->len*x]), tmp->len);
300  strcpy(bit_depth, "16 bits unsigned integer per sample");
301  break;
302 
303  case 32:
304  byte_type = TULONG;
305  img_type = ULONG_IMG;
306  dsp_buffer_copy(tmp->buf, ((unsigned int*)&((unsigned int*)buf)[tmp->len*x]), tmp->len);
307  strcpy(bit_depth, "32 bits unsigned integer per sample");
308  break;
309 
310  case 64:
311  byte_type = TLONGLONG;
312  img_type = LONGLONG_IMG;
313  dsp_buffer_copy(tmp->buf, ((unsigned long*)&((unsigned long*)buf)[tmp->len*x]), tmp->len);
314  strcpy(bit_depth, "64 bits unsigned integer per sample");
315  break;
316 
317  case -32:
318  byte_type = TFLOAT;
319  img_type = FLOAT_IMG;
320  dsp_buffer_copy(tmp->buf, ((float*)&((float*)buf)[tmp->len*x]), tmp->len);
321  strcpy(bit_depth, "32 bits floating point per sample");
322  break;
323 
324  case -64:
325  byte_type = TDOUBLE;
326  img_type = DOUBLE_IMG;
327  dsp_buffer_copy(tmp->buf, ((double*)&((double*)buf)[tmp->len*x]), tmp->len);
328  strcpy(bit_depth, "64 bits floating point per sample");
329  break;
330 
331  default:
332  perr("Unsupported bits per sample value %d", bpp);
333  break;
334  }
336  dsp_stream_free(tmp);
337  }
338 
339  unlink(filename);
340  fitsfile *fptr;
341  fits_create_file(&fptr, filename, &status);
342 
343  if (status)
344  {
345  goto fail_fptr;
346  }
347 
348  fits_create_img(fptr, img_type, naxis, naxes, &status);
349 
350  if (status)
351  {
352  goto fail_fptr;
353  }
354 
355  fits_write_img(fptr, byte_type, 1, nelements, buf, &status);
356 
357  if (status)
358  {
359  goto fail_fptr;
360  }
361 
362  fits_close_file(fptr, &status);
363 
364 fail_fptr:
365  if(status) {
366  fits_get_errstatus(status, error_status);
367  perr("FITS Error: %s\n", error_status);
368  }
369  free(naxes);
370  free (buf);
371 }
372 
373 
374 void dsp_file_write_fits_bayer(const char* filename, int components, int bpp, dsp_stream_p* stream)
375 {
376  int red = 0;
377  int x;
378  dsp_stream_p tmp = dsp_stream_copy(stream[components]);
379  int img_type = USHORT_IMG;
380  int byte_type = TUSHORT;
381  char bit_depth[64] = "16 bits per sample";
382  int len = tmp->len;
383  void* data = malloc((size_t)(len * (size_t)abs(bpp) / 8 + 512));
384  int status = 0;
385  int naxis = tmp->dims;
386  long *naxes = (long*)malloc(sizeof(long) * (size_t)(tmp->dims));
387  long nelements = tmp->len;
388  char error_status[64];
389  int i;
390  for (i = 0; i < tmp->dims; i++)
391  naxes[i] = tmp->sizes[i];
392  dsp_t *buf = dsp_file_composite_2_bayer(stream, red, tmp->sizes[0], tmp->sizes[1]);
394  dsp_stream_free(tmp);
395  for(x = 0; x < components; x++) {
396  dsp_buffer_stretch(buf, stream[components]->len, 0, (dsp_t)(1<<(size_t)abs(bpp))-1);
397  switch (bpp)
398  {
399  case 8:
400  byte_type = TBYTE;
401  img_type = BYTE_IMG;
402  dsp_buffer_copy(buf, ((unsigned char*)((unsigned char*)data)), stream[components]->len);
403  strcpy(bit_depth, "8 bits unsigned integer per sample");
404  break;
405 
406  case 16:
407  byte_type = TUSHORT;
408  img_type = USHORT_IMG;
409  dsp_buffer_copy(buf, ((unsigned short*)((unsigned short*)data)), stream[components]->len);
410  strcpy(bit_depth, "16 bits unsigned integer per sample");
411  break;
412 
413  case 32:
414  byte_type = TULONG;
415  img_type = ULONG_IMG;
416  dsp_buffer_copy(buf, ((unsigned int*)((unsigned int*)data)), stream[components]->len);
417  strcpy(bit_depth, "32 bits unsigned integer per sample");
418  break;
419 
420  case 64:
421  byte_type = TLONGLONG;
422  img_type = LONGLONG_IMG;
423  dsp_buffer_copy(buf, ((unsigned long*)((unsigned long*)data)), stream[components]->len);
424  strcpy(bit_depth, "64 bits unsigned integer per sample");
425  break;
426 
427  case -32:
428  byte_type = TFLOAT;
429  img_type = FLOAT_IMG;
430  dsp_buffer_copy(buf, ((float*)((float*)data)), stream[components]->len);
431  strcpy(bit_depth, "32 bits floating point per sample");
432  break;
433 
434  case -64:
435  byte_type = TDOUBLE;
436  img_type = DOUBLE_IMG;
437  dsp_buffer_copy(buf, ((double*)((double*)data)), stream[components]->len);
438  strcpy(bit_depth, "64 bits floating point per sample");
439  break;
440 
441  default:
442  perr("Unsupported bits per sample value %d", bpp);
443  break;
444  }
445  }
446 
447  unlink(filename);
448  fitsfile *fptr;
449  status = 0;
450  fits_create_file(&fptr, filename, &status);
451 
452  if (status)
453  {
454  goto fail_fptr;
455  }
456 
457  fits_create_img(fptr, img_type, naxis, naxes, &status);
458 
459  if (status)
460  {
461  goto fail_fptr;
462  }
463 
464  int redx = red&1;
465  int redy = (red>>1)&1;
466 
467  fits_write_key(fptr, TINT, "XBAYROFF", &redx, "X Bayer Offset", &status);
468  if (status)
469  {
470  goto fail_fptr;
471  }
472 
473  fits_write_key(fptr, TINT, "YBAYROFF", &redy, "Y Bayer Offset", &status);
474  if (status)
475  {
476  goto fail_fptr;
477  }
478  switch (red) {
479  case 0:
480  fits_write_key(fptr, TSTRING, "BAYERPAT", "RGGB", "Y Bayer Offset", &status);
481  break;
482  case 1:
483  fits_write_key(fptr, TSTRING, "BAYERPAT", "GRGB", "Y Bayer Offset", &status);
484  break;
485  case 2:
486  fits_write_key(fptr, TSTRING, "BAYERPAT", "GBRG", "Y Bayer Offset", &status);
487  break;
488  case 3:
489  fits_write_key(fptr, TSTRING, "BAYERPAT", "BGGR", "Y Bayer Offset", &status);
490  break;
491  }
492  if (status)
493  {
494  goto fail_fptr;
495  }
496 
497  fits_write_img(fptr, byte_type, 1, nelements, data, &status);
498 
499  if (status)
500  {
501  goto fail_fptr;
502  }
503 
504  fits_close_file(fptr, &status);
505 
506 fail_fptr:
507  if(status) {
508  fits_get_errstatus(status, error_status);
509  perr("FITS Error: %s\n", error_status);
510  }
511  free(naxes);
512  free (data);
513 }
514 dsp_stream_p* dsp_file_read_jpeg(const char* filename, int *channels, int stretch)
515 {
516  int width, height;
517  int components;
518  int bpp = 8;
519  unsigned char * buf;
520  struct jpeg_decompress_struct info;
521  struct jpeg_error_mgr err;
522 
523  info.err = jpeg_std_error(& err);
524  jpeg_create_decompress(& info);
525 
526  FILE *jpeg = fopen (filename, "r");
527  if(jpeg == NULL)
528  return NULL;
529 
530  jpeg_stdio_src(&info, jpeg);
531  jpeg_read_header(&info, TRUE);
532  info.dct_method = JDCT_FLOAT;
533  jpeg_start_decompress(&info);
534  width = (int)info.output_width;
535  height = (int)info.output_height;
536  components = (int)info.num_components;
537 
538  int row_stride = (components * width);
539  buf = (unsigned char *)malloc((size_t)(width * height * components));
540  unsigned char *image = buf;
541  int row;
542  for (row = 0; row < height; row++)
543  {
544  jpeg_read_scanlines(&info, &image, 1);
545  image += row_stride;
546  }
547  jpeg_finish_decompress(&info);
548  *channels = components;
549  return dsp_buffer_rgb_to_components(buf, 2, (int[]){width, height}, components, bpp, stretch);
550 }
551 
552 void dsp_file_write_jpeg(const char* filename, int quality, dsp_stream_p stream)
553 {
554  int width = stream->sizes[0];
555  int height = stream->sizes[1];
556  int components = (stream->red>=0) ? 3 : 1;
557  void *buf = malloc((size_t)(stream->len*(stream->red>=0?3:1)));
558  unsigned char *image = (unsigned char *)buf;
559  dsp_t* data = stream->buf;
560 
561  if(components > 1)
562  data = dsp_file_bayer_2_rgb(stream->buf, stream->red, width, height);
563  dsp_buffer_stretch(data, stream->len*(stream->red>=0?3:1), 0, 255);
564  dsp_buffer_copy(data, image, stream->len*(stream->red>=0?3:1));
565 
566  struct jpeg_compress_struct cinfo;
567  struct jpeg_error_mgr jerr;
568  FILE * outfile;
569  int row_stride;
570  cinfo.err = jpeg_std_error(&jerr);
571  jpeg_create_compress(&cinfo);
572  if ((outfile = fopen(filename, "wb")) == NULL) {
573  perr("can't open %s\n", filename);
574  return;
575  }
576  jpeg_stdio_dest(&cinfo, outfile);
577  cinfo.image_width = (unsigned int)width;
578  cinfo.image_height = (unsigned int)height;
579  cinfo.in_color_space = (components == 1 ? JCS_GRAYSCALE : JCS_RGB);
580  cinfo.input_components = components;
581  jpeg_set_defaults(&cinfo);
582  cinfo.dct_method = JDCT_FLOAT;
583  cinfo.optimize_coding = TRUE;
584  cinfo.restart_in_rows = 1;
585 
586  jpeg_set_quality(&cinfo, quality, TRUE);
587  jpeg_start_compress(&cinfo, TRUE);
588  row_stride = width * (stream->red>=0?3:1);
589 
590  int row;
591  for (row = 0; row < height; row++)
592  {
593  jpeg_write_scanlines(&cinfo, &image, 1);
594  image += row_stride;
595  }
596 
597  free(buf);
598  jpeg_finish_compress(&cinfo);
599  fclose(outfile);
600  jpeg_destroy_compress(&cinfo);
601 }
602 
603 void dsp_file_write_jpeg_composite(const char* filename, int components, int quality, dsp_stream_p* stream)
604 {
605  int bpp = 8;
606  unsigned int row_stride;
607  int width = stream[components]->sizes[0];
608  int height = stream[components]->sizes[1];
609  void *buf = malloc((size_t)(stream[components]->len*components));
610  unsigned char *image = (unsigned char *)buf;
611  dsp_buffer_components_to_rgb(stream, buf, components, bpp);
612 
613  struct jpeg_compress_struct cinfo;
614  struct jpeg_error_mgr jerr;
615  FILE * outfile;
616  cinfo.err = jpeg_std_error(&jerr);
617  if ((outfile = fopen(filename, "wb")) == NULL) {
618  perr("can't open %s\n", filename);
619  return;
620  }
621  jpeg_create_compress(&cinfo);
622  jpeg_stdio_dest(&cinfo, outfile);
623  cinfo.image_width = (unsigned int)width;
624  cinfo.image_height = (unsigned int)height;
625  cinfo.in_color_space = (components == 1 ? JCS_GRAYSCALE : JCS_RGB);
626  cinfo.input_components = components;
627  cinfo.dct_method = JDCT_FLOAT;
628  cinfo.optimize_coding = TRUE;
629  cinfo.restart_in_rows = 1;
630  jpeg_set_defaults(&cinfo);
631 
632  jpeg_set_quality(&cinfo, quality, TRUE);
633  jpeg_start_compress(&cinfo, TRUE);
634  row_stride = (unsigned int)width * (unsigned int)components;
635  int row;
636  for (row = 0; row < height; row++)
637  {
638  jpeg_write_scanlines(&cinfo, &image, 1);
639  image += row_stride;
640  }
641 
642  jpeg_finish_compress(&cinfo);
643  jpeg_destroy_compress(&cinfo);
644  fclose(outfile);
645  free(buf);
646 }
647 
648 dsp_t* dsp_file_bayer_2_gray(dsp_t *src, int width, int height)
649 {
650  int i;
651  dsp_t *rawpt, *scanpt;
652  int size;
653 
654  dsp_t * dst = (dsp_t*)malloc(sizeof(dsp_t)*(size_t)(width*height));
655  rawpt = src;
656  scanpt = dst;
657  size = width * height;
658  dsp_t val = 0;
659  for (i = 0; i < size; i++)
660  {
661  if ((i / width) % 2 == 0)
662  {
663  if ((i % 2) == 0)
664  {
665  /* B */
666  if ((i > width) && ((i % width) > 0))
667  {
668  val =
669  (*(rawpt - width - 1) + *(rawpt - width + 1) + *(rawpt + width - 1) + *(rawpt + width + 1)) /
670  4; /* R */
671  val += (*(rawpt - 1) + *(rawpt + 1) + *(rawpt + width) + *(rawpt - width)) / 4; /* G */
672  val += *rawpt; /* B */
673  }
674  else
675  {
676  val = *(rawpt + width + 1); /* R */
677  val += (*(rawpt + 1) + *(rawpt + width)) / 2; /* G */
678  val += *rawpt; /* B */
679  }
680  }
681  else
682  {
683  /* (B)G */
684  if ((i > width) && ((i % width) < (width - 1)))
685  {
686  val = (*(rawpt + width) + *(rawpt - width)) / 2; /* R */
687  val += *rawpt; /* G */
688  val += (*(rawpt - 1) + *(rawpt + 1)) / 2; /* B */
689  }
690  else
691  {
692  val = *(rawpt + width); /* R */
693  val += *rawpt; /* G */
694  val += *(rawpt - 1); /* B */
695  }
696  }
697  }
698  else
699  {
700  if ((i % 2) == 0)
701  {
702  /* G(R) */
703  if ((i < (width * (height - 1))) && ((i % width) > 0))
704  {
705  val = (*(rawpt - 1) + *(rawpt + 1)) / 2; /* R */
706  val += *rawpt; /* G */
707  val += (*(rawpt + width) + *(rawpt - width)) / 2; /* B */
708  }
709  else
710  {
711  /* bottom line or left column */
712  val = *(rawpt + 1); /* R */
713  val += *rawpt; /* G */
714  val += *(rawpt - width); /* B */
715  }
716  }
717  else
718  {
719  /* R */
720  if (i < (width * (height - 1)) && ((i % width) < (width - 1)))
721  {
722  val = *rawpt; /* R */
723  val += (*(rawpt - 1) + *(rawpt + 1) + *(rawpt - width) + *(rawpt + width)) / 4; /* G */
724  val +=
725  (*(rawpt - width - 1) + *(rawpt - width + 1) + *(rawpt + width - 1) + *(rawpt + width + 1)) /
726  4; /* B */
727  }
728  else
729  {
730  /* bottom line or right column */
731  val = *rawpt; /* R */
732  val += (*(rawpt - 1) + *(rawpt - width)) / 2; /* G */
733  val += *(rawpt - width - 1); /* B */
734  }
735  }
736  }
737  *scanpt++ = val;
738  rawpt++;
739  }
740  return dst;
741 }
742 
743 
744 dsp_t* dsp_file_composite_2_bayer(dsp_stream_p *src, int r, int width, int height)
745 {
746  int i;
747  dsp_t *rawpt, *red, *green, *blue;
748  int size;
749  dsp_t *dst = (dsp_t *)malloc(sizeof(dsp_t)*(size_t)(width*height));
750  rawpt = dst;
751  blue = src[0]->buf;
752  green = src[1]->buf;
753  red = src[2]->buf;
754  size = width * height;
755  for (i = 0; i < size; i++)
756  {
757  if ((i / width) % 2 == ((r&2)>>1))
758  {
759  if ((i % 2) == (r&1))
760  {
761  /* B */
762  if ((i > width) && ((i % width) > 0))
763  {
764  rawpt[i - width - 1] += red[i];
765  rawpt[i - width + 1] += red[i];
766  rawpt[i + width - 1] += red[i];
767  rawpt[i + width + 1] += red[i];
768  rawpt[i - 1] += green[i];
769  rawpt[i + 1] += green[i];
770  rawpt[i + width] += green[i];
771  rawpt[i - width] += green[i];
772  rawpt[i] += blue[i];
773  }
774  else
775  {
776  rawpt[i + width + 1] += red[i];
777  rawpt[i + 1] += green[i];
778  rawpt[i + width] += green[i];
779  rawpt[i] += blue[i];
780  }
781  }
782  else
783  {
784  /* (B)G */
785  if ((i > width) && ((i % width) < (width - 1)))
786  {
787  rawpt[i + width] += red[i];
788  rawpt[i - width] += red[i];
789  rawpt[i] += green[i];
790  rawpt[i - 1] += blue[i];
791  rawpt[i + 1] += blue[i];
792  }
793  else
794  {
795  rawpt[i + width] += red[i];
796  rawpt[i] += green[i];
797  rawpt[i - 1] += blue[i];
798  }
799  }
800  }
801  else
802  {
803  if ((i % 2) == (r&1))
804  {
805  /* G(R) */
806  if ((i < (width * (height - 1))) && ((i % width) > 0))
807  {
808  rawpt[i - 1] += red[i];
809  rawpt[i + 1] += red[i];
810  rawpt[i] += green[i];
811  rawpt[i + width] += blue[i];
812  rawpt[i - width] += blue[i];
813  }
814  else
815  {
816  rawpt[i + 1] += red[i];
817  rawpt[i] += green[i];
818  rawpt[i - width] += blue[i];
819  }
820  }
821  else
822  {
823  if (i < (width * (height - 1)) && ((i % width) < (width - 1)))
824  {
825  rawpt[i] = red[i];
826  rawpt[i - 1] += green[i];
827  rawpt[i + 1] += green[i];
828  rawpt[i - width] += green[i];
829  rawpt[i + width] += green[i];
830  rawpt[i - width - 1] += blue[i];
831  rawpt[i - width + 1] += blue[i];
832  rawpt[i + width + 1] += blue[i];
833  rawpt[i + width + 1] += blue[i];
834  }
835  else
836  {
837  /* bottom line or right column */
838  rawpt[i] += red[i];
839  rawpt[i - 1] += green[i];
840  rawpt[i - width] += green[i];
841  rawpt[i - width - 1] += blue[i];
842  }
843  }
844  }
845  }
846  return dst;
847 }
848 
849 
850 dsp_t* dsp_file_bayer_2_composite(dsp_t *src, int r, int width, int height)
851 {
852  int i;
853  dsp_t *rawpt, *red, *green, *blue;
854  int size;
855 
856  dsp_t *dst = (dsp_t *)malloc(sizeof(dsp_t)*(size_t)(width*height*3));
857  rawpt = src;
858  blue = &dst[0];
859  green = &dst[width*height];
860  red = &dst[width*height*2];
861  size = width * height;
862  for (i = 0; i < size; i++)
863  {
864  if ((i / width) % 2 == ((r&2)>>1))
865  {
866  if ((i % 2) == (r&1))
867  {
868  /* B */
869  if ((i > width) && ((i % width) > 0))
870  {
871  *red++ =
872  (*(rawpt - width - 1) + *(rawpt - width + 1) + *(rawpt + width - 1) + *(rawpt + width + 1)) /
873  4; /* R */
874  *green++ = (*(rawpt - 1) + *(rawpt + 1) + *(rawpt + width) + *(rawpt - width)) / 4; /* G */
875  *blue++ = *rawpt; /* B */
876  }
877  else
878  {
879  *red++ = *(rawpt + width + 1); /* R */
880  *green++ = (*(rawpt + 1) + *(rawpt + width)) / 2; /* G */
881  *blue++ = *rawpt; /* B */
882  }
883  }
884  else
885  {
886  /* (B)G */
887  if ((i > width) && ((i % width) < (width - 1)))
888  {
889  *red++ = (*(rawpt + width) + *(rawpt - width)) / 2; /* R */
890  *green++ = *rawpt; /* G */
891  *blue++ = (*(rawpt - 1) + *(rawpt + 1)) / 2; /* B */
892  }
893  else
894  {
895  *red++ = *(rawpt + width); /* R */
896  *green++ = *rawpt; /* G */
897  *blue++ = *(rawpt - 1); /* B */
898  }
899  }
900  }
901  else
902  {
903  if ((i % 2) == (r&1))
904  {
905  /* G(R) */
906  if ((i < (width * (height - 1))) && ((i % width) > 0))
907  {
908  *red++ = (*(rawpt - 1) + *(rawpt + 1)) / 2; /* R */
909  *green++ = *rawpt; /* G */
910  *blue++ = (*(rawpt + width) + *(rawpt - width)) / 2; /* B */
911  }
912  else
913  {
914  /* bottom line or left column */
915  *red++ = *(rawpt + 1); /* R */
916  *green++ = *rawpt; /* G */
917  *blue++ = *(rawpt - width); /* B */
918  }
919  }
920  else
921  {
922  /* R */
923  if (i < (width * (height - 1)) && ((i % width) < (width - 1)))
924  {
925  *red++ = *rawpt; /* R */
926  *green++ = (*(rawpt - 1) + *(rawpt + 1) + *(rawpt - width) + *(rawpt + width)) / 4; /* G */
927  *blue++ =
928  (*(rawpt - width - 1) + *(rawpt - width + 1) + *(rawpt + width - 1) + *(rawpt + width + 1)) /
929  4; /* B */
930  }
931  else
932  {
933  /* bottom line or right column */
934  *red++ = *rawpt; /* R */
935  *green++ = (*(rawpt - 1) + *(rawpt - width)) / 2; /* G */
936  *blue++ = *(rawpt - width - 1); /* B */
937  }
938  }
939  }
940  rawpt++;
941  }
942  return dst;
943 }
944 
945 
946 dsp_t* dsp_file_bayer_2_rgb(dsp_t *src, int red, int width, int height)
947 {
948  int i;
949  dsp_t *rawpt, *scanpt;
950  int size;
951 
952  dsp_t * dst = (dsp_t*)malloc(sizeof(dsp_t)*(size_t)(width*height*3));
953  rawpt = src;
954  scanpt = dst;
955  size = width * height;
956  (void)*scanpt++;
957  for (i = 0; i < size; i++)
958  {
959  if ((i / width) % 2 == ((red&2)>>1))
960  {
961  if ((i % 2) == (red&1))
962  {
963  /* B */
964  if ((i > width) && ((i % width) > 0))
965  {
966  *scanpt++ =
967  (*(rawpt - width - 1) + *(rawpt - width + 1) + *(rawpt + width - 1) + *(rawpt + width + 1)) /
968  4; /* R */
969  *scanpt++ = (*(rawpt - 1) + *(rawpt + 1) + *(rawpt + width) + *(rawpt - width)) / 4; /* G */
970  *scanpt++ = *rawpt; /* B */
971  }
972  else
973  {
974  *scanpt++ = *(rawpt + width + 1); /* R */
975  *scanpt++ = (*(rawpt + 1) + *(rawpt + width)) / 2; /* G */
976  *scanpt++ = *rawpt; /* B */
977  }
978  }
979  else
980  {
981  /* (B)G */
982  if ((i > width) && ((i % width) < (width - 1)))
983  {
984  *scanpt++ = (*(rawpt + width) + *(rawpt - width)) / 2; /* R */
985  *scanpt++ = *rawpt; /* G */
986  *scanpt++ = (*(rawpt - 1) + *(rawpt + 1)) / 2; /* B */
987  }
988  else
989  {
990  *scanpt++ = *(rawpt + width); /* R */
991  *scanpt++ = *rawpt; /* G */
992  *scanpt++ = *(rawpt - 1); /* B */
993  }
994  }
995  }
996  else
997  {
998  if ((i % 2) == (red&1))
999  {
1000  /* G(R) */
1001  if ((i < (width * (height - 1))) && ((i % width) > 0))
1002  {
1003  *scanpt++ = (*(rawpt - 1) + *(rawpt + 1)) / 2; /* R */
1004  *scanpt++ = *rawpt; /* G */
1005  *scanpt++ = (*(rawpt + width) + *(rawpt - width)) / 2; /* B */
1006  }
1007  else
1008  {
1009  /* bottom line or left column */
1010  *scanpt++ = *(rawpt + 1); /* R */
1011  *scanpt++ = *rawpt; /* G */
1012  *scanpt++ = *(rawpt - width); /* B */
1013  }
1014  }
1015  else
1016  {
1017  /* R */
1018  if (i < (width * (height - 1)) && ((i % width) < (width - 1)))
1019  {
1020  *scanpt++ = *rawpt; /* R */
1021  *scanpt++ = (*(rawpt - 1) + *(rawpt + 1) + *(rawpt - width) + *(rawpt + width)) / 4; /* G */
1022  *scanpt++ =
1023  (*(rawpt - width - 1) + *(rawpt - width + 1) + *(rawpt + width - 1) + *(rawpt + width + 1)) /
1024  4; /* B */
1025  }
1026  else
1027  {
1028  /* bottom line or right column */
1029  *scanpt++ = *rawpt; /* R */
1030  *scanpt++ = (*(rawpt - 1) + *(rawpt - width)) / 2; /* G */
1031  *scanpt++ = *(rawpt - width - 1); /* B */
1032  }
1033  }
1034  }
1035  rawpt++;
1036  }
1037  return dst;
1038 }
1039 
1040 dsp_stream_p *dsp_stream_from_components(dsp_t* buf, int dims, int *sizes, int components)
1041 {
1042  int d, y, x;
1043  dsp_stream_p* picture = (dsp_stream_p*) malloc(sizeof(dsp_stream_p)*(size_t)(components+1));
1044  for(y = 0; y <= components; y++) {
1045  picture[y] = dsp_stream_new();
1046  for(d = 0; d < dims; d++)
1047  dsp_stream_add_dim(picture[y], sizes[d]);
1048  dsp_stream_alloc_buffer(picture[y], picture[y]->len);
1049  if(y < components) {
1050  dsp_buffer_copy(((dsp_t*)&buf[y*picture[y]->len]), picture[y]->buf, picture[y]->len);
1051  } else {
1052  for(x = 0; x < picture[y]->len; x++) {
1053  double val = 0;
1054  for(d = 0; d < components; d++) {
1055  val += buf[x+d*picture[y]->len];
1056  }
1057  picture[y]->buf[x] = (dsp_t)val/components;
1058  }
1059  }
1060  }
1061  return picture;
1062 }
1063 
1064 dsp_stream_p *dsp_buffer_rgb_to_components(void* buf, int dims, int *sizes, int components, int bpp, int stretch)
1065 {
1066  dsp_stream_p* picture = (dsp_stream_p*) malloc(sizeof(dsp_stream_p)*(size_t)(components+1));
1067  int x, y, z, d;
1068  dsp_stream_p channel;
1069  for(y = 0; y < components; y++) {
1070  channel = dsp_stream_new();
1071  for(d = 0; d < dims; d++)
1072  dsp_stream_add_dim(channel, sizes[d]);
1073  dsp_stream_alloc_buffer(channel, channel->len);
1074  switch(bpp) {
1075  case 8:
1076  dsp_buffer_copy_stepping(((unsigned char*)&((unsigned char*)buf)[y]), channel->buf, channel->len*components, channel->len, components, 1);
1077  break;
1078  case 16:
1079  dsp_buffer_copy_stepping(((unsigned short*)&((unsigned short*)buf)[y]), channel->buf, channel->len*components, channel->len, components, 1);
1080  break;
1081  case 32:
1082  dsp_buffer_copy_stepping(((unsigned int*)&((unsigned int*)buf)[y]), channel->buf, channel->len*components, channel->len, components, 1);
1083  break;
1084  case 64:
1085  dsp_buffer_copy_stepping(((unsigned long*)&((unsigned long*)buf)[y]), channel->buf, channel->len*components, channel->len, components, 1);
1086  break;
1087  case -32:
1088  dsp_buffer_copy_stepping(((float*)&((float*)buf)[y]), channel->buf, channel->len*components, channel->len, components, 1);
1089  break;
1090  case -64:
1091  dsp_buffer_copy_stepping(((double*)&((double*)buf)[y]), channel->buf, channel->len*components, channel->len, components, 1);
1092  break;
1093  default:
1094  break;
1095  }
1096  if(stretch) {
1097  dsp_buffer_stretch(channel->buf, channel->len, 0, dsp_t_max);
1098  }
1099  picture[y] = channel;
1100  }
1101  channel = dsp_stream_new();
1102  for(d = 0; d < dims; d++)
1103  dsp_stream_add_dim(channel, sizes[d]);
1104  dsp_stream_alloc_buffer(channel, channel->len);
1105  for(x = 0; x < channel->len; x++) {
1106  double v = 0;
1107  switch(bpp) {
1108  case 8:
1109  for(z = 0; z < y; z++) v += ((unsigned char*)buf)[x*y+z];
1110  break;
1111  case 16:
1112  for(z = 0; z < y; z++) v += ((unsigned short*)buf)[x*y+z];
1113  break;
1114  case 32:
1115  for(z = 0; z < y; z++) v += ((unsigned int*)buf)[x*y+z];
1116  break;
1117  case 64:
1118  for(z = 0; z < y; z++) v += ((unsigned long*)buf)[x*y+z];
1119  break;
1120  case -32:
1121  for(z = 0; z < y; z++) v += ((float*)buf)[x*y+z];
1122  break;
1123  case -64:
1124  for(z = 0; z < y; z++) v += ((double*)buf)[x*y+z];
1125  break;
1126  default:
1127  break;
1128  }
1129  channel->buf[x] = (dsp_t)(v/y);
1130  }
1131  if(stretch) {
1132  dsp_buffer_stretch(channel->buf, channel->len, 0, dsp_t_max);
1133  }
1134  picture[y] = channel;
1135  free (buf);
1136  return picture;
1137 }
1138 
1139 void dsp_buffer_components_to_rgb(dsp_stream_p *stream, void* rgb, int components, int bpp)
1140 {
1141  ssize_t y;
1142  int len = stream[0]->len * components;
1143  dsp_t max = (dsp_t)((double)((1<<(size_t)abs(bpp))-1));
1144  max = Min(max, dsp_t_max);
1145  for(y = 0; y < components; y++) {
1146  dsp_stream_p in = dsp_stream_copy(stream[y]);
1147  dsp_buffer_stretch(in->buf, in->len, 0, max);
1148  switch(bpp) {
1149  case 8:
1150  dsp_buffer_copy_stepping(in->buf, ((unsigned char*)&((unsigned char*)rgb)[y]), in->len, len, 1, components);
1151  break;
1152  case 16:
1153  dsp_buffer_copy_stepping(in->buf, ((unsigned short*)&((unsigned short*)rgb)[y]), in->len, len, 1, components);
1154  break;
1155  case 32:
1156  dsp_buffer_copy_stepping(in->buf, ((unsigned int*)&((unsigned int*)rgb)[y]), in->len, len, 1, components);
1157  break;
1158  case 64:
1159  dsp_buffer_copy_stepping(in->buf, ((unsigned long*)&((unsigned long*)rgb)[y]), in->len, len, 1, components);
1160  break;
1161  case -32:
1162  dsp_buffer_copy_stepping(in->buf, ((float*)&((float*)rgb)[y]), in->len, len, 1, components);
1163  break;
1164  case -64:
1165  dsp_buffer_copy_stepping(in->buf, ((double*)&((double*)rgb)[y]), in->len, len, 1, components);
1166  break;
1167  }
1169  dsp_stream_free(in);
1170  }
1171 }
double max(void)
#define perr(...)
Definition: dsp.h:153
double dsp_t
Definition: dsp.h:69
#define dsp_t_max
Definition: dsp.h:71
#define Min(a, b)
if min() is not present you can use this one
Definition: dsp.h:172
int * sizes
Sizes of each dimension.
Definition: dsp.h:373
int red
Red pixel (Bayer)
Definition: dsp.h:401
int dims
Number of dimensions of the buffers.
Definition: dsp.h:371
dsp_t * buf
buffer
Definition: dsp.h:375
int len
The buffers length.
Definition: dsp.h:369
#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:1038
#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:1059
#define dsp_buffer_stretch(buf, len, _mn, _mx)
Stretch minimum and maximum values of the input stream.
Definition: dsp.h:792
void dsp_buffer_pow1(dsp_stream_p stream, double val)
Expose elements of the input stream to the given power.
Definition: buffer.c:205
DLL_EXPORT void dsp_stream_free(dsp_stream_p stream)
Free the DSP stream passed as argument.
Definition: stream.c:163
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:192
DLL_EXPORT dsp_stream_p dsp_stream_new(void)
Allocate a new DSP stream type.
Definition: stream.c:124
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:227
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:78
DLL_EXPORT void dsp_stream_free_buffer(dsp_stream_p stream)
Free the buffer of the DSP Stream passed as argument.
Definition: stream.c:112
void dsp_buffer_components_to_rgb(dsp_stream_p *stream, void *rgb, int components, int bpp)
Convert a component dsp_stream_p array into an RGB dsp_t array.
Definition: file.c:1139
void dsp_file_write_jpeg_composite(const char *filename, int components, int quality, dsp_stream_p *stream)
Write the components dsp_stream_p array into a JPEG file,.
Definition: file.c:603
dsp_stream_p * dsp_stream_from_components(dsp_t *buf, int dims, int *sizes, int components)
Convert a color component dsp_t array into a dsp_stream_p array each element containing the single co...
Definition: file.c:1040
dsp_t * dsp_file_bayer_2_gray(dsp_t *src, int width, int height)
Convert a bayer pattern dsp_t array into a grayscale array.
Definition: file.c:648
void dsp_file_write_fits_bayer(const char *filename, int components, int bpp, dsp_stream_p *stream)
Write a FITS file from a dsp_stream_p array.
Definition: file.c:374
void dsp_file_write_fits_composite(const char *filename, int components, int bpp, dsp_stream_p *stream)
Write the components dsp_stream_p array into a JPEG file,.
Definition: file.c:266
dsp_t * dsp_file_bayer_2_rgb(dsp_t *src, int red, int width, int height)
Convert a bayer pattern dsp_t array into a ordered 3 RGB array.
Definition: file.c:946
void dsp_file_write_jpeg(const char *filename, int quality, dsp_stream_p stream)
Write the stream into a JPEG file,.
Definition: file.c:552
dsp_stream_p * dsp_buffer_rgb_to_components(void *buf, int dims, int *sizes, int components, int bpp, int stretch)
Convert an RGB color dsp_t array into a dsp_stream_p array each element containing the single compone...
Definition: file.c:1064
dsp_t * dsp_file_composite_2_bayer(dsp_stream_p *src, int r, int width, int height)
Convert a component dsp_stream_p array into a bayer dsp_t array.
Definition: file.c:744
dsp_stream_p * dsp_file_read_jpeg(const 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:514
dsp_stream_p * dsp_file_read_fits(const char *filename, int *channels, int stretch)
Read a FITS file and fill a dsp_stream_p with its content.
Definition: file.c:27
void dsp_file_write_fits(const char *filename, int bpp, dsp_stream_p stream)
Write the dsp_stream_p into a FITS file,.
Definition: file.c:164
dsp_t * dsp_file_bayer_2_composite(dsp_t *src, int r, int width, int height)
Convert a bayer pattern dsp_t array into a contiguos component array.
Definition: file.c:850
Contains a set of informations and data relative to a buffer and how to use it.
Definition: dsp.h:363
#define TRUE
Definition: stvdriver.c:54