Instrument Neutral Distributed Interface INDI 2.2.5
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
27dsp_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 }
156fail:
157 if(status) {
158 fits_get_errstatus(status, error_status);
159 perr("FITS Error: %s\n", error_status);
160 }
161 return NULL;
162}
163
164void 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
254fail_fptr:
255 if(status) {
256 fits_get_errstatus(status, error_status);
257 perr("FITS Error: %s\n", error_status);
258 }
259fail:
260 dsp_stream_free_buffer(tmp);
261 dsp_stream_free(tmp);
262 free(naxes);
263 free (buf);
264}
265
266void 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 }
335 dsp_stream_free_buffer(tmp);
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
364fail_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
374void 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]);
393 dsp_stream_free_buffer(tmp);
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
506fail_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}
514dsp_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
552void 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
603void 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
648dsp_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
744dsp_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
850dsp_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
946dsp_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
1040dsp_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
1064dsp_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
1139void 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 }
1168 dsp_stream_free_buffer(in);
1169 dsp_stream_free(in);
1170 }
1171}
max
double max(void)
perr
#define perr(...)
Definition dsp.h:161
dsp_t
double dsp_t
Definition dsp.h:76
dsp_t_max
#define dsp_t_max
Definition dsp.h:78
Min
#define Min(a, b)
if min() is not present you can use this one
Definition dsp.h:180
dsp_stream_t::sizes
int * sizes
Sizes of each dimension.
Definition dsp.h:389
dsp_stream_t::red
int red
Red pixel (Bayer)
Definition dsp.h:417
dsp_stream_t::dims
int dims
Number of dimensions of the buffers.
Definition dsp.h:387
dsp_stream_t::buf
dsp_t * buf
buffer
Definition dsp.h:391
dsp_stream_t::len
int len
The buffers length.
Definition dsp.h:385
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:1054
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:1075
dsp_buffer_stretch
#define dsp_buffer_stretch(buf, len, _mn, _mx)
Stretch minimum and maximum values of the input stream.
Definition dsp.h:808
dsp_buffer_pow1
void dsp_buffer_pow1(dsp_stream_p stream, double val)
Expose elements of the input stream to the given power.
Definition buffer.c:205
dsp_stream_free
DLL_EXPORT void dsp_stream_free(dsp_stream_p stream)
Free the DSP stream passed as argument.
Definition stream.c:167
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:200
dsp_stream_new
DLL_EXPORT dsp_stream_p dsp_stream_new(void)
Allocate a new DSP stream type.
Definition stream.c:128
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:235
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:78
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:112
dsp_stream_from_components
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_buffer_components_to_rgb
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
dsp_file_write_jpeg_composite
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_buffer_rgb_to_components
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_file_read_jpeg
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_file_write_fits_bayer
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
dsp_file_read_fits
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
dsp_file_write_fits_composite
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_file_bayer_2_composite
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 contiguous component array.
Definition file.c:850
dsp_file_write_jpeg
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_file_bayer_2_rgb
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
dsp_file_bayer_2_gray
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
dsp_file_composite_2_bayer
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_file_write_fits
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_stream_t
Contains a set of information and data relative to a buffer and how to use it.
Definition dsp.h:379
TRUE
#define TRUE
Definition stvdriver.c:54
Generated on Thu Oct 22 2026 09:22:00 for Instrument Neutral Distributed Interface INDI by doxygen 1.9.8