diff --git a/kstars/fitsviewer/bayer.c b/kstars/fitsviewer/bayer.c index 8267e0650..5bb2228db 100644 --- a/kstars/fitsviewer/bayer.c +++ b/kstars/fitsviewer/bayer.c @@ -1,2584 +1,2562 @@ /* * 1394-Based Digital Camera Control Library * * Bayer pattern decoding functions * * Written by Damien Douxchamps and Frederic Devernay * The original VNG and AHD Bayer decoding are from Dave Coffin's DCRAW. * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA */ #include "bayer.h" #include #include #include #include #include #define CLIP(in, out) \ in = in < 0 ? 0 : in; \ in = in > 255 ? 255 : in; \ out = in; #define CLIP16(in, out, bits) \ in = in < 0 ? 0 : in; \ in = in > ((1 << bits) - 1) ? ((1 << bits) - 1) : in; \ out = in; void ClearBorders(uint8_t *rgb, int sx, int sy, int w) { int i, j; /* black edges are added with a width w: */ i = 3 * sx * w - 1; j = 3 * sx * sy - 1; while (i >= 0) { rgb[i--] = 0; rgb[j--] = 0; } int low = sx * (w - 1) * 3 - 1 + w * 3; i = low + sx * (sy - w * 2 + 1) * 3; while (i > low) { j = 6 * w; while (j > 0) { rgb[i--] = 0; j--; } i -= (sx - 2 * w) * 3; } } void ClearBorders_uint16(uint16_t *rgb, int sx, int sy, int w) { int i, j; /* black edges: */ i = 3 * sx * w - 1; j = 3 * sx * sy - 1; while (i >= 0) { rgb[i--] = 0; rgb[j--] = 0; } int low = sx * (w - 1) * 3 - 1 + w * 3; i = low + sx * (sy - w * 2 + 1) * 3; while (i > low) { j = 6 * w; while (j > 0) { rgb[i--] = 0; j--; } i -= (sx - 2 * w) * 3; } } /************************************************************** * Color conversion functions for cameras that can * * output raw-Bayer pattern images, such as some Basler and * * Point Grey camera. Most of the algos presented here come * * from http://www-ise.stanford.edu/~tingchen/ and have been * * converted from Matlab to C and extended to all elementary * * patterns. * **************************************************************/ /* 8-bits versions */ /* insprired by OpenCV's Bayer decoding */ dc1394error_t dc1394_bayer_NearestNeighbor(const uint8_t *bayer, uint8_t *rgb, int sx, int sy, int tile) { const int bayerStep = sx; const int rgbStep = 3 * sx; int width = sx; int height = sy; int blue = tile == DC1394_COLOR_FILTER_BGGR || tile == DC1394_COLOR_FILTER_GBRG ? -1 : 1; int start_with_green = tile == DC1394_COLOR_FILTER_GBRG || tile == DC1394_COLOR_FILTER_GRBG; int i, imax, iinc; if ((tile > DC1394_COLOR_FILTER_MAX) || (tile < DC1394_COLOR_FILTER_MIN)) return DC1394_INVALID_COLOR_FILTER; /* add black border */ imax = sx * sy * 3; for (i = sx * (sy - 1) * 3; i < imax; i++) { rgb[i] = 0; } iinc = (sx - 1) * 3; for (i = (sx - 1) * 3; i < imax; i += iinc) { rgb[i++] = 0; rgb[i++] = 0; rgb[i++] = 0; } rgb += 1; width -= 1; height -= 1; for (; height--; bayer += bayerStep, rgb += rgbStep) { const uint8_t *bayerEnd = bayer + width; if (start_with_green) { rgb[-blue] = bayer[1]; rgb[0] = bayer[bayerStep + 1]; rgb[blue] = bayer[bayerStep]; bayer++; rgb += 3; } if (blue > 0) { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { rgb[-1] = bayer[0]; rgb[0] = bayer[1]; rgb[1] = bayer[bayerStep + 1]; rgb[2] = bayer[2]; rgb[3] = bayer[bayerStep + 2]; rgb[4] = bayer[bayerStep + 1]; } } else { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { rgb[1] = bayer[0]; rgb[0] = bayer[1]; rgb[-1] = bayer[bayerStep + 1]; rgb[4] = bayer[2]; rgb[3] = bayer[bayerStep + 2]; rgb[2] = bayer[bayerStep + 1]; } } if (bayer < bayerEnd) { rgb[-blue] = bayer[0]; rgb[0] = bayer[1]; rgb[blue] = bayer[bayerStep + 1]; bayer++; rgb += 3; } bayer -= width; rgb -= width * 3; blue = -blue; start_with_green = !start_with_green; } return DC1394_SUCCESS; } /* OpenCV's Bayer decoding */ dc1394error_t dc1394_bayer_Bilinear(const uint8_t *bayer, uint8_t *rgb, int sx, int sy, int tile) { const int bayerStep = sx; const int rgbStep = 3 * sx; int width = sx; int height = sy; /* the two letters of the OpenCV name are respectively the 4th and 3rd letters from the blinky name, and we also have to switch R and B (OpenCV is BGR) CV_BayerBG2BGR <-> DC1394_COLOR_FILTER_BGGR CV_BayerGB2BGR <-> DC1394_COLOR_FILTER_GBRG CV_BayerGR2BGR <-> DC1394_COLOR_FILTER_GRBG int blue = tile == CV_BayerBG2BGR || tile == CV_BayerGB2BGR ? -1 : 1; int start_with_green = tile == CV_BayerGB2BGR || tile == CV_BayerGR2BGR; */ int blue = tile == DC1394_COLOR_FILTER_BGGR || tile == DC1394_COLOR_FILTER_GBRG ? -1 : 1; int start_with_green = tile == DC1394_COLOR_FILTER_GBRG || tile == DC1394_COLOR_FILTER_GRBG; if ((tile > DC1394_COLOR_FILTER_MAX) || (tile < DC1394_COLOR_FILTER_MIN)) return DC1394_INVALID_COLOR_FILTER; ClearBorders(rgb, sx, sy, 1); rgb += rgbStep + 3 + 1; height -= 2; width -= 2; for (; height--; bayer += bayerStep, rgb += rgbStep) { int t0, t1; const uint8_t *bayerEnd = bayer + width; if (start_with_green) { /* OpenCV has a bug in the next line, which was t0 = (bayer[0] + bayer[bayerStep * 2] + 1) >> 1; */ t0 = (bayer[1] + bayer[bayerStep * 2 + 1] + 1) >> 1; t1 = (bayer[bayerStep] + bayer[bayerStep + 2] + 1) >> 1; rgb[-blue] = (uint8_t)t0; rgb[0] = bayer[bayerStep + 1]; rgb[blue] = (uint8_t)t1; bayer++; rgb += 3; } if (blue > 0) { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2] + bayer[bayerStep * 2 + 2] + 2) >> 2; t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2] + bayer[bayerStep * 2 + 1] + 2) >> 2; rgb[-1] = (uint8_t)t0; rgb[0] = (uint8_t)t1; rgb[1] = bayer[bayerStep + 1]; t0 = (bayer[2] + bayer[bayerStep * 2 + 2] + 1) >> 1; t1 = (bayer[bayerStep + 1] + bayer[bayerStep + 3] + 1) >> 1; rgb[2] = (uint8_t)t0; rgb[3] = bayer[bayerStep + 2]; rgb[4] = (uint8_t)t1; } } else { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2] + bayer[bayerStep * 2 + 2] + 2) >> 2; t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2] + bayer[bayerStep * 2 + 1] + 2) >> 2; rgb[1] = (uint8_t)t0; rgb[0] = (uint8_t)t1; rgb[-1] = bayer[bayerStep + 1]; t0 = (bayer[2] + bayer[bayerStep * 2 + 2] + 1) >> 1; t1 = (bayer[bayerStep + 1] + bayer[bayerStep + 3] + 1) >> 1; rgb[4] = (uint8_t)t0; rgb[3] = bayer[bayerStep + 2]; rgb[2] = (uint8_t)t1; } } if (bayer < bayerEnd) { t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2] + bayer[bayerStep * 2 + 2] + 2) >> 2; t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2] + bayer[bayerStep * 2 + 1] + 2) >> 2; rgb[-blue] = (uint8_t)t0; rgb[0] = (uint8_t)t1; rgb[blue] = bayer[bayerStep + 1]; bayer++; rgb += 3; } bayer -= width; rgb -= width * 3; blue = -blue; start_with_green = !start_with_green; } return DC1394_SUCCESS; } /* High-Quality Linear Interpolation For Demosaicing Of Bayer-Patterned Color Images, by Henrique S. Malvar, Li-wei He, and Ross Cutler, in ICASSP'04 */ dc1394error_t dc1394_bayer_HQLinear(const uint8_t *bayer, uint8_t *rgb, int sx, int sy, int tile) { const int bayerStep = sx; const int rgbStep = 3 * sx; int width = sx; int height = sy; int blue = tile == DC1394_COLOR_FILTER_BGGR || tile == DC1394_COLOR_FILTER_GBRG ? -1 : 1; int start_with_green = tile == DC1394_COLOR_FILTER_GBRG || tile == DC1394_COLOR_FILTER_GRBG; if ((tile > DC1394_COLOR_FILTER_MAX) || (tile < DC1394_COLOR_FILTER_MIN)) return DC1394_INVALID_COLOR_FILTER; ClearBorders(rgb, sx, sy, 2); rgb += 2 * rgbStep + 6 + 1; height -= 4; width -= 4; /* We begin with a (+1 line,+1 column) offset with respect to bilinear decoding, so start_with_green is the same, but blue is opposite */ blue = -blue; for (; height--; bayer += bayerStep, rgb += rgbStep) { int t0, t1; const uint8_t *bayerEnd = bayer + width; const int bayerStep2 = bayerStep * 2; const int bayerStep3 = bayerStep * 3; const int bayerStep4 = bayerStep * 4; if (start_with_green) { /* at green pixel */ rgb[0] = bayer[bayerStep2 + 2]; t0 = rgb[0] * 5 + ((bayer[bayerStep + 2] + bayer[bayerStep3 + 2]) << 2) - bayer[2] - bayer[bayerStep + 1] - bayer[bayerStep + 3] - bayer[bayerStep3 + 1] - bayer[bayerStep3 + 3] - bayer[bayerStep4 + 2] + ((bayer[bayerStep2] + bayer[bayerStep2 + 4] + 1) >> 1); t1 = rgb[0] * 5 + ((bayer[bayerStep2 + 1] + bayer[bayerStep2 + 3]) << 2) - bayer[bayerStep2] - bayer[bayerStep + 1] - bayer[bayerStep + 3] - bayer[bayerStep3 + 1] - bayer[bayerStep3 + 3] - bayer[bayerStep2 + 4] + ((bayer[2] + bayer[bayerStep4 + 2] + 1) >> 1); t0 = (t0 + 4) >> 3; CLIP(t0, rgb[-blue]); t1 = (t1 + 4) >> 3; CLIP(t1, rgb[blue]); bayer++; rgb += 3; } if (blue > 0) { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { /* B at B */ rgb[1] = bayer[bayerStep2 + 2]; /* R at B */ t0 = ((bayer[bayerStep + 1] + bayer[bayerStep + 3] + bayer[bayerStep3 + 1] + bayer[bayerStep3 + 3]) << 1) - (((bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) * 3 + 1) >> 1) + rgb[1] * 6; /* G at B */ t1 = ((bayer[bayerStep + 2] + bayer[bayerStep2 + 1] + bayer[bayerStep2 + 3] + bayer[bayerStep3 + 2]) << 1) - (bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) + (rgb[1] << 2); t0 = (t0 + 4) >> 3; CLIP(t0, rgb[-1]); t1 = (t1 + 4) >> 3; CLIP(t1, rgb[0]); /* at green pixel */ rgb[3] = bayer[bayerStep2 + 3]; t0 = rgb[3] * 5 + ((bayer[bayerStep + 3] + bayer[bayerStep3 + 3]) << 2) - bayer[3] - bayer[bayerStep + 2] - bayer[bayerStep + 4] - bayer[bayerStep3 + 2] - bayer[bayerStep3 + 4] - bayer[bayerStep4 + 3] + ((bayer[bayerStep2 + 1] + bayer[bayerStep2 + 5] + 1) >> 1); t1 = rgb[3] * 5 + ((bayer[bayerStep2 + 2] + bayer[bayerStep2 + 4]) << 2) - bayer[bayerStep2 + 1] - bayer[bayerStep + 2] - bayer[bayerStep + 4] - bayer[bayerStep3 + 2] - bayer[bayerStep3 + 4] - bayer[bayerStep2 + 5] + ((bayer[3] + bayer[bayerStep4 + 3] + 1) >> 1); t0 = (t0 + 4) >> 3; CLIP(t0, rgb[2]); t1 = (t1 + 4) >> 3; CLIP(t1, rgb[4]); } } else { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { /* R at R */ rgb[-1] = bayer[bayerStep2 + 2]; /* B at R */ t0 = ((bayer[bayerStep + 1] + bayer[bayerStep + 3] + bayer[bayerStep3 + 1] + bayer[bayerStep3 + 3]) << 1) - (((bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) * 3 + 1) >> 1) + rgb[-1] * 6; /* G at R */ t1 = ((bayer[bayerStep + 2] + bayer[bayerStep2 + 1] + bayer[bayerStep2 + 3] + bayer[bayerStep * 3 + 2]) << 1) - (bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) + (rgb[-1] << 2); t0 = (t0 + 4) >> 3; CLIP(t0, rgb[1]); t1 = (t1 + 4) >> 3; CLIP(t1, rgb[0]); /* at green pixel */ rgb[3] = bayer[bayerStep2 + 3]; t0 = rgb[3] * 5 + ((bayer[bayerStep + 3] + bayer[bayerStep3 + 3]) << 2) - bayer[3] - bayer[bayerStep + 2] - bayer[bayerStep + 4] - bayer[bayerStep3 + 2] - bayer[bayerStep3 + 4] - bayer[bayerStep4 + 3] + ((bayer[bayerStep2 + 1] + bayer[bayerStep2 + 5] + 1) >> 1); t1 = rgb[3] * 5 + ((bayer[bayerStep2 + 2] + bayer[bayerStep2 + 4]) << 2) - bayer[bayerStep2 + 1] - bayer[bayerStep + 2] - bayer[bayerStep + 4] - bayer[bayerStep3 + 2] - bayer[bayerStep3 + 4] - bayer[bayerStep2 + 5] + ((bayer[3] + bayer[bayerStep4 + 3] + 1) >> 1); t0 = (t0 + 4) >> 3; CLIP(t0, rgb[4]); t1 = (t1 + 4) >> 3; CLIP(t1, rgb[2]); } } if (bayer < bayerEnd) { /* B at B */ rgb[blue] = bayer[bayerStep2 + 2]; /* R at B */ t0 = ((bayer[bayerStep + 1] + bayer[bayerStep + 3] + bayer[bayerStep3 + 1] + bayer[bayerStep3 + 3]) << 1) - (((bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) * 3 + 1) >> 1) + rgb[blue] * 6; /* G at B */ t1 = (((bayer[bayerStep + 2] + bayer[bayerStep2 + 1] + bayer[bayerStep2 + 3] + bayer[bayerStep3 + 2])) << 1) - (bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) + (rgb[blue] << 2); t0 = (t0 + 4) >> 3; CLIP(t0, rgb[-blue]); t1 = (t1 + 4) >> 3; CLIP(t1, rgb[0]); bayer++; rgb += 3; } bayer -= width; rgb -= width * 3; blue = -blue; start_with_green = !start_with_green; } return DC1394_SUCCESS; } /* coriander's Bayer decoding */ /* Edge Sensing Interpolation II from http://www-ise.stanford.edu/~tingchen/ */ /* (Laroche,Claude A. "Apparatus and method for adaptively interpolating a full color image utilizing chrominance gradients" U.S. Patent 5,373,322) */ dc1394error_t dc1394_bayer_EdgeSense(const uint8_t *bayer, uint8_t *rgb, int sx, int sy, int tile) { uint8_t *outR, *outG, *outB; register int i3, j3, base; int i, j; int dh, dv; int tmp; int sx3 = sx * 3; /* sx and sy should be even */ switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_BGGR: outR = &rgb[0]; outG = &rgb[1]; outB = &rgb[2]; break; case DC1394_COLOR_FILTER_GBRG: case DC1394_COLOR_FILTER_RGGB: outR = &rgb[2]; outG = &rgb[1]; outB = &rgb[0]; break; default: return DC1394_INVALID_COLOR_FILTER; } switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_GBRG: /* copy original RGB data to output images */ for (i = 0, i3 = 0; i < sy * sx; i += (sx << 1), i3 += (sx3 << 1)) { for (j = 0, j3 = 0; j < sx; j += 2, j3 += 6) { base = i3 + j3; outG[base] = bayer[i + j]; outG[base + sx3 + 3] = bayer[i + j + sx + 1]; outR[base + 3] = bayer[i + j + 1]; outB[base + sx3] = bayer[i + j + sx]; } } /* process GREEN channel */ for (i3 = 3 * sx3; i3 < (sy - 2) * sx3; i3 += (sx3 << 1)) { for (j3 = 6; j3 < sx3 - 9; j3 += 6) { base = i3 + j3; dh = abs(((outB[base - 6] + outB[base + 6]) >> 1) - outB[base]); dv = abs(((outB[base - (sx3 << 1)] + outB[base + (sx3 << 1)]) >> 1) - outB[base]); tmp = (((outG[base - 3] + outG[base + 3]) >> 1) * (dh <= dv) + ((outG[base - sx3] + outG[base + sx3]) >> 1) * (dh > dv)); CLIP(tmp, outG[base]); } } for (i3 = 2 * sx3; i3 < (sy - 3) * sx3; i3 += (sx3 << 1)) { for (j3 = 9; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; dh = abs(((outR[base - 6] + outR[base + 6]) >> 1) - outR[base]); dv = abs(((outR[base - (sx3 << 1)] + outR[base + (sx3 << 1)]) >> 1) - outR[base]); tmp = (((outG[base - 3] + outG[base + 3]) >> 1) * (dh <= dv) + ((outG[base - sx3] + outG[base + sx3]) >> 1) * (dh > dv)); CLIP(tmp, outG[base]); } } /* process RED channel */ for (i3 = 0; i3 < (sy - 1) * sx3; i3 += (sx3 << 1)) { for (j3 = 6; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outR[base - 3] - outG[base - 3] + outR[base + 3] - outG[base + 3]) >> 1); CLIP(tmp, outR[base]); } } for (i3 = sx3; i3 < (sy - 2) * sx3; i3 += (sx3 << 1)) { for (j3 = 3; j3 < sx3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outR[base - sx3] - outG[base - sx3] + outR[base + sx3] - outG[base + sx3]) >> 1); CLIP(tmp, outR[base]); } for (j3 = 6; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outR[base - sx3 - 3] - outG[base - sx3 - 3] + outR[base - sx3 + 3] - outG[base - sx3 + 3] + outR[base + sx3 - 3] - outG[base + sx3 - 3] + outR[base + sx3 + 3] - outG[base + sx3 + 3]) >> 2); CLIP(tmp, outR[base]); } } /* process BLUE channel */ for (i3 = sx3; i3 < sy * sx3; i3 += (sx3 << 1)) { for (j3 = 3; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - 3] - outG[base - 3] + outB[base + 3] - outG[base + 3]) >> 1); CLIP(tmp, outB[base]); } } for (i3 = 2 * sx3; i3 < (sy - 1) * sx3; i3 += (sx3 << 1)) { for (j3 = 0; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - sx3] - outG[base - sx3] + outB[base + sx3] - outG[base + sx3]) >> 1); CLIP(tmp, outB[base]); } for (j3 = 3; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - sx3 - 3] - outG[base - sx3 - 3] + outB[base - sx3 + 3] - outG[base - sx3 + 3] + outB[base + sx3 - 3] - outG[base + sx3 - 3] + outB[base + sx3 + 3] - outG[base + sx3 + 3]) >> 2); CLIP(tmp, outB[base]); } } break; case DC1394_COLOR_FILTER_BGGR: case DC1394_COLOR_FILTER_RGGB: /* copy original RGB data to output images */ for (i = 0, i3 = 0; i < sy * sx; i += (sx << 1), i3 += (sx3 << 1)) { for (j = 0, j3 = 0; j < sx; j += 2, j3 += 6) { base = i3 + j3; outB[base] = bayer[i + j]; outR[base + sx3 + 3] = bayer[i + sx + (j + 1)]; outG[base + 3] = bayer[i + j + 1]; outG[base + sx3] = bayer[i + sx + j]; } } /* process GREEN channel */ for (i3 = 2 * sx3; i3 < (sy - 2) * sx3; i3 += (sx3 << 1)) { for (j3 = 6; j3 < sx3 - 9; j3 += 6) { base = i3 + j3; dh = abs(((outB[base - 6] + outB[base + 6]) >> 1) - outB[base]); dv = abs(((outB[base - (sx3 << 1)] + outB[base + (sx3 << 1)]) >> 1) - outB[base]); tmp = (((outG[base - 3] + outG[base + 3]) >> 1) * (dh <= dv) + ((outG[base - sx3] + outG[base + sx3]) >> 1) * (dh > dv)); CLIP(tmp, outG[base]); } } for (i3 = 3 * sx3; i3 < (sy - 3) * sx3; i3 += (sx3 << 1)) { for (j3 = 9; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; dh = abs(((outR[base - 6] + outR[base + 6]) >> 1) - outR[base]); dv = abs(((outR[base - (sx3 << 1)] + outR[base + (sx3 << 1)]) >> 1) - outR[base]); tmp = (((outG[base - 3] + outG[base + 3]) >> 1) * (dh <= dv) + ((outG[base - sx3] + outG[base + sx3]) >> 1) * (dh > dv)); CLIP(tmp, outG[base]); } } /* process RED channel */ for (i3 = sx3; i3 < (sy - 1) * sx3; i3 += (sx3 << 1)) { /* G-points (1/2) */ for (j3 = 6; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outR[base - 3] - outG[base - 3] + outR[base + 3] - outG[base + 3]) >> 1); CLIP(tmp, outR[base]); } } for (i3 = 2 * sx3; i3 < (sy - 2) * sx3; i3 += (sx3 << 1)) { for (j3 = 3; j3 < sx3; j3 += 6) { /* G-points (2/2) */ base = i3 + j3; tmp = outG[base] + ((outR[base - sx3] - outG[base - sx3] + outR[base + sx3] - outG[base + sx3]) >> 1); CLIP(tmp, outR[base]); } for (j3 = 6; j3 < sx3 - 3; j3 += 6) { /* B-points */ base = i3 + j3; tmp = outG[base] + ((outR[base - sx3 - 3] - outG[base - sx3 - 3] + outR[base - sx3 + 3] - outG[base - sx3 + 3] + outR[base + sx3 - 3] - outG[base + sx3 - 3] + outR[base + sx3 + 3] - outG[base + sx3 + 3]) >> 2); CLIP(tmp, outR[base]); } } /* process BLUE channel */ for (i = 0, i3 = 0; i < sy * sx; i += (sx << 1), i3 += (sx3 << 1)) { for (j = 1, j3 = 3; j < sx - 2; j += 2, j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - 3] - outG[base - 3] + outB[base + 3] - outG[base + 3]) >> 1); CLIP(tmp, outB[base]); } } for (i3 = sx3; i3 < (sy - 1) * sx3; i3 += (sx3 << 1)) { for (j3 = 0; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - sx3] - outG[base - sx3] + outB[base + sx3] - outG[base + sx3]) >> 1); CLIP(tmp, outB[base]); } for (j3 = 3; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - sx3 - 3] - outG[base - sx3 - 3] + outB[base - sx3 + 3] - outG[base - sx3 + 3] + outB[base + sx3 - 3] - outG[base + sx3 - 3] + outB[base + sx3 + 3] - outG[base + sx3 + 3]) >> 2); CLIP(tmp, outB[base]); } } break; } ClearBorders(rgb, sx, sy, 3); return DC1394_SUCCESS; } /* coriander's Bayer decoding */ dc1394error_t dc1394_bayer_Downsample(const uint8_t *bayer, uint8_t *rgb, int sx, int sy, int tile) { uint8_t *outR, *outG, *outB; register int i, j; int tmp; switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_BGGR: outR = &rgb[0]; outG = &rgb[1]; outB = &rgb[2]; break; case DC1394_COLOR_FILTER_GBRG: case DC1394_COLOR_FILTER_RGGB: outR = &rgb[2]; outG = &rgb[1]; outB = &rgb[0]; break; default: return DC1394_INVALID_COLOR_FILTER; } switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_GBRG: for (i = 0; i < sy * sx; i += (sx << 1)) { for (j = 0; j < sx; j += 2) { tmp = ((bayer[i + j] + bayer[i + sx + j + 1]) >> 1); CLIP(tmp, outG[((i >> 2) + (j >> 1)) * 3]); tmp = bayer[i + j + 1]; CLIP(tmp, outR[((i >> 2) + (j >> 1)) * 3]); tmp = bayer[i + sx + j]; CLIP(tmp, outB[((i >> 2) + (j >> 1)) * 3]); } } break; case DC1394_COLOR_FILTER_BGGR: case DC1394_COLOR_FILTER_RGGB: for (i = 0; i < sy * sx; i += (sx << 1)) { for (j = 0; j < sx; j += 2) { tmp = ((bayer[i + sx + j] + bayer[i + j + 1]) >> 1); CLIP(tmp, outG[((i >> 2) + (j >> 1)) * 3]); tmp = bayer[i + sx + j + 1]; CLIP(tmp, outR[((i >> 2) + (j >> 1)) * 3]); tmp = bayer[i + j]; CLIP(tmp, outB[((i >> 2) + (j >> 1)) * 3]); } } break; } return DC1394_SUCCESS; } /* this is the method used inside AVT cameras. See AVT docs. */ dc1394error_t dc1394_bayer_Simple(const uint8_t *bayer, uint8_t *rgb, int sx, int sy, int tile) { const int bayerStep = sx; const int rgbStep = 3 * sx; int width = sx; int height = sy; int blue = tile == DC1394_COLOR_FILTER_BGGR || tile == DC1394_COLOR_FILTER_GBRG ? -1 : 1; int start_with_green = tile == DC1394_COLOR_FILTER_GBRG || tile == DC1394_COLOR_FILTER_GRBG; int i, imax, iinc; if ((tile > DC1394_COLOR_FILTER_MAX) || (tile < DC1394_COLOR_FILTER_MIN)) return DC1394_INVALID_COLOR_FILTER; /* add black border */ imax = sx * sy * 3; for (i = sx * (sy - 1) * 3; i < imax; i++) { rgb[i] = 0; } iinc = (sx - 1) * 3; for (i = (sx - 1) * 3; i < imax; i += iinc) { rgb[i++] = 0; rgb[i++] = 0; rgb[i++] = 0; } rgb += 1; width -= 1; height -= 1; for (; height--; bayer += bayerStep, rgb += rgbStep) { const uint8_t *bayerEnd = bayer + width; if (start_with_green) { rgb[-blue] = bayer[1]; rgb[0] = (bayer[0] + bayer[bayerStep + 1] + 1) >> 1; rgb[blue] = bayer[bayerStep]; bayer++; rgb += 3; } if (blue > 0) { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { rgb[-1] = bayer[0]; rgb[0] = (bayer[1] + bayer[bayerStep] + 1) >> 1; rgb[1] = bayer[bayerStep + 1]; rgb[2] = bayer[2]; rgb[3] = (bayer[1] + bayer[bayerStep + 2] + 1) >> 1; rgb[4] = bayer[bayerStep + 1]; } } else { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { rgb[1] = bayer[0]; rgb[0] = (bayer[1] + bayer[bayerStep] + 1) >> 1; rgb[-1] = bayer[bayerStep + 1]; rgb[4] = bayer[2]; rgb[3] = (bayer[1] + bayer[bayerStep + 2] + 1) >> 1; rgb[2] = bayer[bayerStep + 1]; } } if (bayer < bayerEnd) { rgb[-blue] = bayer[0]; rgb[0] = (bayer[1] + bayer[bayerStep] + 1) >> 1; rgb[blue] = bayer[bayerStep + 1]; bayer++; rgb += 3; } bayer -= width; rgb -= width * 3; blue = -blue; start_with_green = !start_with_green; } return DC1394_SUCCESS; } /* 16-bits versions */ dc1394error_t dc1394_bayer16_RGBX_NearestNeighbor(const uint16_t *bayer, uint16_t *rgbx, int sx, int sy, int tile) { const int bayerStep = sx; const int rgbStep = 4 * sx; int width = sx; int height = sy; int blue = (tile == DC1394_COLOR_FILTER_BGGR || tile == DC1394_COLOR_FILTER_GBRG) ? -1 : 1; int start_with_green = (tile == DC1394_COLOR_FILTER_GBRG || tile == DC1394_COLOR_FILTER_GRBG) ? 1 : 0; if ((tile > DC1394_COLOR_FILTER_MAX) || (tile < DC1394_COLOR_FILTER_MIN)) return DC1394_INVALID_COLOR_FILTER; - /* add black border */ -#if 0 - int i, iinc, imax; - imax = sx * sy * 4; - for (i = sx * (sy - 1) * 4; i < imax; i++) - { - rgbx[i] = 0; - } - iinc = (sx - 1) * 4; - for (i = (sx - 1) * 4; i < imax; i += iinc) - { - rgbx[i++] = 0; - rgbx[i++] = 0; - rgbx[i++] = 0; - rgbx[i++] = 0; - } - - rgbx += 1; - height -= 1; - width -= 1; -#endif - for (; height--; bayer += bayerStep, rgbx += rgbStep) { const uint16_t *bayerEnd = bayer + width; if (start_with_green) { rgbx[-blue] = bayer[1]; rgbx[0] = bayer[bayerStep + 1]; rgbx[blue] = bayer[bayerStep]; rgbx[2] = 0xFFFF; bayer++; rgbx += 4; } if (blue > 0) { for (; bayer <= bayerEnd - 2; bayer += 2, rgbx += 8) { rgbx[-1] = bayer[0]; rgbx[0] = bayer[1]; rgbx[1] = bayer[bayerStep + 1]; rgbx[2] = 0xFFFF; rgbx[3] = bayer[2]; rgbx[4] = bayer[bayerStep + 2]; rgbx[5] = bayer[bayerStep + 1]; rgbx[6] = 0xFFFF; } } else { for (; bayer <= bayerEnd - 2; bayer += 2, rgbx += 8) { rgbx[2] = 0xFFFF; rgbx[1] = bayer[0]; rgbx[0] = bayer[1]; rgbx[-1] = bayer[bayerStep + 1]; rgbx[6] = 0xFFFF; rgbx[5] = bayer[2]; rgbx[4] = bayer[bayerStep + 2]; rgbx[3] = bayer[bayerStep + 1]; } } if (bayer < bayerEnd) { rgbx[-blue] = bayer[0]; rgbx[0] = bayer[1]; rgbx[blue] = bayer[bayerStep + 1]; rgbx[2] = 0xFFFF; bayer++; rgbx += 4; } bayer -= width; rgbx -= width * 4; blue = -blue; start_with_green = !start_with_green; } return DC1394_SUCCESS; } /* insprired by OpenCV's Bayer decoding */ dc1394error_t dc1394_bayer_NearestNeighbor_uint16(const uint16_t *bayer, uint16_t *rgb, int sx, int sy, int tile, int bits) { (void)bits; const int bayerStep = sx; const int rgbStep = 3 * sx; int width = sx; int height = sy; int blue = tile == DC1394_COLOR_FILTER_BGGR || tile == DC1394_COLOR_FILTER_GBRG ? -1 : 1; int start_with_green = tile == DC1394_COLOR_FILTER_GBRG || tile == DC1394_COLOR_FILTER_GRBG; int i, iinc, imax; if ((tile > DC1394_COLOR_FILTER_MAX) || (tile < DC1394_COLOR_FILTER_MIN)) return DC1394_INVALID_COLOR_FILTER; /* add black border */ imax = sx * sy * 3; for (i = sx * (sy - 1) * 3; i < imax; i++) { rgb[i] = 0; } iinc = (sx - 1) * 3; for (i = (sx - 1) * 3; i < imax; i += iinc) { rgb[i++] = 0; rgb[i++] = 0; rgb[i++] = 0; } rgb += 1; height -= 1; width -= 1; for (; height--; bayer += bayerStep, rgb += rgbStep) { const uint16_t *bayerEnd = bayer + width; if (start_with_green) { rgb[-blue] = bayer[1]; rgb[0] = bayer[bayerStep + 1]; rgb[blue] = bayer[bayerStep]; bayer++; rgb += 3; } if (blue > 0) { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { rgb[-1] = bayer[0]; rgb[0] = bayer[1]; rgb[1] = bayer[bayerStep + 1]; rgb[2] = bayer[2]; rgb[3] = bayer[bayerStep + 2]; rgb[4] = bayer[bayerStep + 1]; } } else { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { rgb[1] = bayer[0]; rgb[0] = bayer[1]; rgb[-1] = bayer[bayerStep + 1]; rgb[4] = bayer[2]; rgb[3] = bayer[bayerStep + 2]; rgb[2] = bayer[bayerStep + 1]; } } if (bayer < bayerEnd) { rgb[-blue] = bayer[0]; rgb[0] = bayer[1]; rgb[blue] = bayer[bayerStep + 1]; bayer++; rgb += 3; } bayer -= width; rgb -= width * 3; blue = -blue; start_with_green = !start_with_green; } return DC1394_SUCCESS; } /* OpenCV's Bayer decoding */ dc1394error_t dc1394_bayer_Bilinear_uint16(const uint16_t *bayer, uint16_t *rgb, int sx, int sy, int tile, int bits) { (void)bits; const int bayerStep = sx; const int rgbStep = 3 * sx; int width = sx; int height = sy; int blue = tile == DC1394_COLOR_FILTER_BGGR || tile == DC1394_COLOR_FILTER_GBRG ? -1 : 1; int start_with_green = tile == DC1394_COLOR_FILTER_GBRG || tile == DC1394_COLOR_FILTER_GRBG; if ((tile > DC1394_COLOR_FILTER_MAX) || (tile < DC1394_COLOR_FILTER_MIN)) return DC1394_INVALID_COLOR_FILTER; rgb += rgbStep + 3 + 1; height -= 2; width -= 2; for (; height--; bayer += bayerStep, rgb += rgbStep) { int t0, t1; const uint16_t *bayerEnd = bayer + width; if (start_with_green) { /* OpenCV has a bug in the next line, which was t0 = (bayer[0] + bayer[bayerStep * 2] + 1) >> 1; */ t0 = (bayer[1] + bayer[bayerStep * 2 + 1] + 1) >> 1; t1 = (bayer[bayerStep] + bayer[bayerStep + 2] + 1) >> 1; rgb[-blue] = (uint16_t)t0; rgb[0] = bayer[bayerStep + 1]; rgb[blue] = (uint16_t)t1; bayer++; rgb += 3; } if (blue > 0) { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2] + bayer[bayerStep * 2 + 2] + 2) >> 2; t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2] + bayer[bayerStep * 2 + 1] + 2) >> 2; rgb[-1] = (uint16_t)t0; rgb[0] = (uint16_t)t1; rgb[1] = bayer[bayerStep + 1]; t0 = (bayer[2] + bayer[bayerStep * 2 + 2] + 1) >> 1; t1 = (bayer[bayerStep + 1] + bayer[bayerStep + 3] + 1) >> 1; rgb[2] = (uint16_t)t0; rgb[3] = bayer[bayerStep + 2]; rgb[4] = (uint16_t)t1; } } else { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2] + bayer[bayerStep * 2 + 2] + 2) >> 2; t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2] + bayer[bayerStep * 2 + 1] + 2) >> 2; rgb[1] = (uint16_t)t0; rgb[0] = (uint16_t)t1; rgb[-1] = bayer[bayerStep + 1]; t0 = (bayer[2] + bayer[bayerStep * 2 + 2] + 1) >> 1; t1 = (bayer[bayerStep + 1] + bayer[bayerStep + 3] + 1) >> 1; rgb[4] = (uint16_t)t0; rgb[3] = bayer[bayerStep + 2]; rgb[2] = (uint16_t)t1; } } if (bayer < bayerEnd) { t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2] + bayer[bayerStep * 2 + 2] + 2) >> 2; t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2] + bayer[bayerStep * 2 + 1] + 2) >> 2; rgb[-blue] = (uint16_t)t0; rgb[0] = (uint16_t)t1; rgb[blue] = bayer[bayerStep + 1]; bayer++; rgb += 3; } bayer -= width; rgb -= width * 3; blue = -blue; start_with_green = !start_with_green; } return DC1394_SUCCESS; } /* High-Quality Linear Interpolation For Demosaicing Of Bayer-Patterned Color Images, by Henrique S. Malvar, Li-wei He, and Ross Cutler, in ICASSP'04 */ dc1394error_t dc1394_bayer_HQLinear_uint16(const uint16_t *bayer, uint16_t *rgb, int sx, int sy, int tile, int bits) { const int bayerStep = sx; const int rgbStep = 3 * sx; int width = sx; int height = sy; /* the two letters of the OpenCV name are respectively the 4th and 3rd letters from the blinky name, and we also have to switch R and B (OpenCV is BGR) CV_BayerBG2BGR <-> DC1394_COLOR_FILTER_BGGR CV_BayerGB2BGR <-> DC1394_COLOR_FILTER_GBRG CV_BayerGR2BGR <-> DC1394_COLOR_FILTER_GRBG int blue = tile == CV_BayerBG2BGR || tile == CV_BayerGB2BGR ? -1 : 1; int start_with_green = tile == CV_BayerGB2BGR || tile == CV_BayerGR2BGR; */ int blue = tile == DC1394_COLOR_FILTER_BGGR || tile == DC1394_COLOR_FILTER_GBRG ? -1 : 1; int start_with_green = tile == DC1394_COLOR_FILTER_GBRG || tile == DC1394_COLOR_FILTER_GRBG; if ((tile > DC1394_COLOR_FILTER_MAX) || (tile < DC1394_COLOR_FILTER_MIN)) return DC1394_INVALID_COLOR_FILTER; ClearBorders_uint16(rgb, sx, sy, 2); rgb += 2 * rgbStep + 6 + 1; height -= 4; width -= 4; /* We begin with a (+1 line,+1 column) offset with respect to bilinear decoding, so start_with_green is the same, but blue is opposite */ blue = -blue; for (; height--; bayer += bayerStep, rgb += rgbStep) { int t0, t1; const uint16_t *bayerEnd = bayer + width; const int bayerStep2 = bayerStep * 2; const int bayerStep3 = bayerStep * 3; const int bayerStep4 = bayerStep * 4; if (start_with_green) { /* at green pixel */ rgb[0] = bayer[bayerStep2 + 2]; t0 = rgb[0] * 5 + ((bayer[bayerStep + 2] + bayer[bayerStep3 + 2]) << 2) - bayer[2] - bayer[bayerStep + 1] - bayer[bayerStep + 3] - bayer[bayerStep3 + 1] - bayer[bayerStep3 + 3] - bayer[bayerStep4 + 2] + ((bayer[bayerStep2] + bayer[bayerStep2 + 4] + 1) >> 1); t1 = rgb[0] * 5 + ((bayer[bayerStep2 + 1] + bayer[bayerStep2 + 3]) << 2) - bayer[bayerStep2] - bayer[bayerStep + 1] - bayer[bayerStep + 3] - bayer[bayerStep3 + 1] - bayer[bayerStep3 + 3] - bayer[bayerStep2 + 4] + ((bayer[2] + bayer[bayerStep4 + 2] + 1) >> 1); t0 = (t0 + 4) >> 3; CLIP16(t0, rgb[-blue], bits); t1 = (t1 + 4) >> 3; CLIP16(t1, rgb[blue], bits); bayer++; rgb += 3; } if (blue > 0) { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { /* B at B */ rgb[1] = bayer[bayerStep2 + 2]; /* R at B */ t0 = ((bayer[bayerStep + 1] + bayer[bayerStep + 3] + bayer[bayerStep3 + 1] + bayer[bayerStep3 + 3]) << 1) - (((bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) * 3 + 1) >> 1) + rgb[1] * 6; /* G at B */ t1 = ((bayer[bayerStep + 2] + bayer[bayerStep2 + 1] + bayer[bayerStep2 + 3] + bayer[bayerStep * 3 + 2]) << 1) - (bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) + (rgb[1] << 2); t0 = (t0 + 4) >> 3; CLIP16(t0, rgb[-1], bits); t1 = (t1 + 4) >> 3; CLIP16(t1, rgb[0], bits); /* at green pixel */ rgb[3] = bayer[bayerStep2 + 3]; t0 = rgb[3] * 5 + ((bayer[bayerStep + 3] + bayer[bayerStep3 + 3]) << 2) - bayer[3] - bayer[bayerStep + 2] - bayer[bayerStep + 4] - bayer[bayerStep3 + 2] - bayer[bayerStep3 + 4] - bayer[bayerStep4 + 3] + ((bayer[bayerStep2 + 1] + bayer[bayerStep2 + 5] + 1) >> 1); t1 = rgb[3] * 5 + ((bayer[bayerStep2 + 2] + bayer[bayerStep2 + 4]) << 2) - bayer[bayerStep2 + 1] - bayer[bayerStep + 2] - bayer[bayerStep + 4] - bayer[bayerStep3 + 2] - bayer[bayerStep3 + 4] - bayer[bayerStep2 + 5] + ((bayer[3] + bayer[bayerStep4 + 3] + 1) >> 1); t0 = (t0 + 4) >> 3; CLIP16(t0, rgb[2], bits); t1 = (t1 + 4) >> 3; CLIP16(t1, rgb[4], bits); } } else { for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) { /* R at R */ rgb[-1] = bayer[bayerStep2 + 2]; /* B at R */ t0 = ((bayer[bayerStep + 1] + bayer[bayerStep + 3] + bayer[bayerStep * 3 + 1] + bayer[bayerStep3 + 3]) << 1) - (((bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) * 3 + 1) >> 1) + rgb[-1] * 6; /* G at R */ t1 = ((bayer[bayerStep + 2] + bayer[bayerStep2 + 1] + bayer[bayerStep2 + 3] + bayer[bayerStep3 + 2]) << 1) - (bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) + (rgb[-1] << 2); t0 = (t0 + 4) >> 3; CLIP16(t0, rgb[1], bits); t1 = (t1 + 4) >> 3; CLIP16(t1, rgb[0], bits); /* at green pixel */ rgb[3] = bayer[bayerStep2 + 3]; t0 = rgb[3] * 5 + ((bayer[bayerStep + 3] + bayer[bayerStep3 + 3]) << 2) - bayer[3] - bayer[bayerStep + 2] - bayer[bayerStep + 4] - bayer[bayerStep3 + 2] - bayer[bayerStep3 + 4] - bayer[bayerStep4 + 3] + ((bayer[bayerStep2 + 1] + bayer[bayerStep2 + 5] + 1) >> 1); t1 = rgb[3] * 5 + ((bayer[bayerStep2 + 2] + bayer[bayerStep2 + 4]) << 2) - bayer[bayerStep2 + 1] - bayer[bayerStep + 2] - bayer[bayerStep + 4] - bayer[bayerStep3 + 2] - bayer[bayerStep3 + 4] - bayer[bayerStep2 + 5] + ((bayer[3] + bayer[bayerStep4 + 3] + 1) >> 1); t0 = (t0 + 4) >> 3; CLIP16(t0, rgb[4], bits); t1 = (t1 + 4) >> 3; CLIP16(t1, rgb[2], bits); } } if (bayer < bayerEnd) { /* B at B */ rgb[blue] = bayer[bayerStep2 + 2]; /* R at B */ t0 = ((bayer[bayerStep + 1] + bayer[bayerStep + 3] + bayer[bayerStep3 + 1] + bayer[bayerStep3 + 3]) << 1) - (((bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) * 3 + 1) >> 1) + rgb[blue] * 6; /* G at B */ t1 = (((bayer[bayerStep + 2] + bayer[bayerStep2 + 1] + bayer[bayerStep2 + 3] + bayer[bayerStep3 + 2])) << 1) - (bayer[2] + bayer[bayerStep2] + bayer[bayerStep2 + 4] + bayer[bayerStep4 + 2]) + (rgb[blue] << 2); t0 = (t0 + 4) >> 3; CLIP16(t0, rgb[-blue], bits); t1 = (t1 + 4) >> 3; CLIP16(t1, rgb[0], bits); bayer++; rgb += 3; } bayer -= width; rgb -= width * 3; blue = -blue; start_with_green = !start_with_green; } return DC1394_SUCCESS; } /* coriander's Bayer decoding */ dc1394error_t dc1394_bayer_EdgeSense_uint16(const uint16_t *bayer, uint16_t *rgb, int sx, int sy, int tile, int bits) { uint16_t *outR, *outG, *outB; register int i3, j3, base; int i, j; int dh, dv; int tmp; int sx3 = sx * 3; /* sx and sy should be even */ switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_BGGR: outR = &rgb[0]; outG = &rgb[1]; outB = &rgb[2]; break; case DC1394_COLOR_FILTER_GBRG: case DC1394_COLOR_FILTER_RGGB: outR = &rgb[2]; outG = &rgb[1]; outB = &rgb[0]; break; default: return DC1394_INVALID_COLOR_FILTER; } switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_GBRG: /* copy original RGB data to output images */ for (i = 0, i3 = 0; i < sy * sx; i += (sx << 1), i3 += (sx3 << 1)) { for (j = 0, j3 = 0; j < sx; j += 2, j3 += 6) { base = i3 + j3; outG[base] = bayer[i + j]; outG[base + sx3 + 3] = bayer[i + j + sx + 1]; outR[base + 3] = bayer[i + j + 1]; outB[base + sx3] = bayer[i + j + sx]; } } /* process GREEN channel */ for (i3 = 3 * sx3; i3 < (sy - 2) * sx3; i3 += (sx3 << 1)) { for (j3 = 6; j3 < sx3 - 9; j3 += 6) { base = i3 + j3; dh = abs(((outB[base - 6] + outB[base + 6]) >> 1) - outB[base]); dv = abs(((outB[base - (sx3 << 1)] + outB[base + (sx3 << 1)]) >> 1) - outB[base]); tmp = (((outG[base - 3] + outG[base + 3]) >> 1) * (dh <= dv) + ((outG[base - sx3] + outG[base + sx3]) >> 1) * (dh > dv)); CLIP16(tmp, outG[base], bits); } } for (i3 = 2 * sx3; i3 < (sy - 3) * sx3; i3 += (sx3 << 1)) { for (j3 = 9; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; dh = abs(((outR[base - 6] + outR[base + 6]) >> 1) - outR[base]); dv = abs(((outR[base - (sx3 << 1)] + outR[base + (sx3 << 1)]) >> 1) - outR[base]); tmp = (((outG[base - 3] + outG[base + 3]) >> 1) * (dh <= dv) + ((outG[base - sx3] + outG[base + sx3]) >> 1) * (dh > dv)); CLIP16(tmp, outG[base], bits); } } /* process RED channel */ for (i3 = 0; i3 < (sy - 1) * sx3; i3 += (sx3 << 1)) { for (j3 = 6; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outR[base - 3] - outG[base - 3] + outR[base + 3] - outG[base + 3]) >> 1); CLIP16(tmp, outR[base], bits); } } for (i3 = sx3; i3 < (sy - 2) * sx3; i3 += (sx3 << 1)) { for (j3 = 3; j3 < sx3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outR[base - sx3] - outG[base - sx3] + outR[base + sx3] - outG[base + sx3]) >> 1); CLIP16(tmp, outR[base], bits); } for (j3 = 6; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outR[base - sx3 - 3] - outG[base - sx3 - 3] + outR[base - sx3 + 3] - outG[base - sx3 + 3] + outR[base + sx3 - 3] - outG[base + sx3 - 3] + outR[base + sx3 + 3] - outG[base + sx3 + 3]) >> 2); CLIP16(tmp, outR[base], bits); } } /* process BLUE channel */ for (i3 = sx3; i3 < sy * sx3; i3 += (sx3 << 1)) { for (j3 = 3; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - 3] - outG[base - 3] + outB[base + 3] - outG[base + 3]) >> 1); CLIP16(tmp, outB[base], bits); } } for (i3 = 2 * sx3; i3 < (sy - 1) * sx3; i3 += (sx3 << 1)) { for (j3 = 0; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - sx3] - outG[base - sx3] + outB[base + sx3] - outG[base + sx3]) >> 1); CLIP16(tmp, outB[base], bits); } for (j3 = 3; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - sx3 - 3] - outG[base - sx3 - 3] + outB[base - sx3 + 3] - outG[base - sx3 + 3] + outB[base + sx3 - 3] - outG[base + sx3 - 3] + outB[base + sx3 + 3] - outG[base + sx3 + 3]) >> 2); CLIP16(tmp, outB[base], bits); } } break; case DC1394_COLOR_FILTER_BGGR: case DC1394_COLOR_FILTER_RGGB: /* copy original RGB data to output images */ for (i = 0, i3 = 0; i < sy * sx; i += (sx << 1), i3 += (sx3 << 1)) { for (j = 0, j3 = 0; j < sx; j += 2, j3 += 6) { base = i3 + j3; outB[base] = bayer[i + j]; outR[base + sx3 + 3] = bayer[i + sx + (j + 1)]; outG[base + 3] = bayer[i + j + 1]; outG[base + sx3] = bayer[i + sx + j]; } } /* process GREEN channel */ for (i3 = 2 * sx3; i3 < (sy - 2) * sx3; i3 += (sx3 << 1)) { for (j3 = 6; j3 < sx3 - 9; j3 += 6) { base = i3 + j3; dh = abs(((outB[base - 6] + outB[base + 6]) >> 1) - outB[base]); dv = abs(((outB[base - (sx3 << 1)] + outB[base + (sx3 << 1)]) >> 1) - outB[base]); tmp = (((outG[base - 3] + outG[base + 3]) >> 1) * (dh <= dv) + ((outG[base - sx3] + outG[base + sx3]) >> 1) * (dh > dv)); CLIP16(tmp, outG[base], bits); } } for (i3 = 3 * sx3; i3 < (sy - 3) * sx3; i3 += (sx3 << 1)) { for (j3 = 9; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; dh = abs(((outR[base - 6] + outR[base + 6]) >> 1) - outR[base]); dv = abs(((outR[base - (sx3 << 1)] + outR[base + (sx3 << 1)]) >> 1) - outR[base]); tmp = (((outG[base - 3] + outG[base + 3]) >> 1) * (dh <= dv) + ((outG[base - sx3] + outG[base + sx3]) >> 1) * (dh > dv)); CLIP16(tmp, outG[base], bits); } } /* process RED channel */ for (i3 = sx3; i3 < (sy - 1) * sx3; i3 += (sx3 << 1)) { /* G-points (1/2) */ for (j3 = 6; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outR[base - 3] - outG[base - 3] + outR[base + 3] - outG[base + 3]) >> 1); CLIP16(tmp, outR[base], bits); } } for (i3 = 2 * sx3; i3 < (sy - 2) * sx3; i3 += (sx3 << 1)) { for (j3 = 3; j3 < sx3; j3 += 6) { /* G-points (2/2) */ base = i3 + j3; tmp = outG[base] + ((outR[base - sx3] - outG[base - sx3] + outR[base + sx3] - outG[base + sx3]) >> 1); CLIP16(tmp, outR[base], bits); } for (j3 = 6; j3 < sx3 - 3; j3 += 6) { /* B-points */ base = i3 + j3; tmp = outG[base] + ((outR[base - sx3 - 3] - outG[base - sx3 - 3] + outR[base - sx3 + 3] - outG[base - sx3 + 3] + outR[base + sx3 - 3] - outG[base + sx3 - 3] + outR[base + sx3 + 3] - outG[base + sx3 + 3]) >> 2); CLIP16(tmp, outR[base], bits); } } /* process BLUE channel */ for (i = 0, i3 = 0; i < sy * sx; i += (sx << 1), i3 += (sx3 << 1)) { for (j = 1, j3 = 3; j < sx - 2; j += 2, j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - 3] - outG[base - 3] + outB[base + 3] - outG[base + 3]) >> 1); CLIP16(tmp, outB[base], bits); } } for (i3 = sx3; i3 < (sy - 1) * sx3; i3 += (sx3 << 1)) { for (j3 = 0; j3 < sx3 - 3; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - sx3] - outG[base - sx3] + outB[base + sx3] - outG[base + sx3]) >> 1); CLIP16(tmp, outB[base], bits); } for (j3 = 3; j3 < sx3 - 6; j3 += 6) { base = i3 + j3; tmp = outG[base] + ((outB[base - sx3 - 3] - outG[base - sx3 - 3] + outB[base - sx3 + 3] - outG[base - sx3 + 3] + outB[base + sx3 - 3] - outG[base + sx3 - 3] + outB[base + sx3 + 3] - outG[base + sx3 + 3]) >> 2); CLIP16(tmp, outB[base], bits); } } break; } ClearBorders_uint16(rgb, sx, sy, 3); return DC1394_SUCCESS; } /* coriander's Bayer decoding */ dc1394error_t dc1394_bayer_Downsample_uint16(const uint16_t *bayer, uint16_t *rgb, int sx, int sy, int tile, int bits) { uint16_t *outR, *outG, *outB; register int i, j; int tmp; switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_BGGR: outR = &rgb[0]; outG = &rgb[1]; outB = &rgb[2]; break; case DC1394_COLOR_FILTER_GBRG: case DC1394_COLOR_FILTER_RGGB: outR = &rgb[2]; outG = &rgb[1]; outB = &rgb[0]; break; default: return DC1394_INVALID_COLOR_FILTER; } switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_GBRG: for (i = 0; i < sy * sx; i += (sx << 1)) { for (j = 0; j < sx; j += 2) { tmp = ((bayer[i + j] + bayer[i + sx + j + 1]) >> 1); CLIP16(tmp, outG[((i >> 2) + (j >> 1)) * 3], bits); tmp = bayer[i + sx + j + 1]; CLIP16(tmp, outR[((i >> 2) + (j >> 1)) * 3], bits); tmp = bayer[i + sx + j]; CLIP16(tmp, outB[((i >> 2) + (j >> 1)) * 3], bits); } } break; case DC1394_COLOR_FILTER_BGGR: case DC1394_COLOR_FILTER_RGGB: for (i = 0; i < sy * sx; i += (sx << 1)) { for (j = 0; j < sx; j += 2) { tmp = ((bayer[i + sx + j] + bayer[i + j + 1]) >> 1); CLIP16(tmp, outG[((i >> 2) + (j >> 1)) * 3], bits); tmp = bayer[i + sx + j + 1]; CLIP16(tmp, outR[((i >> 2) + (j >> 1)) * 3], bits); tmp = bayer[i + j]; CLIP16(tmp, outB[((i >> 2) + (j >> 1)) * 3], bits); } } break; } return DC1394_SUCCESS; } /* coriander's Bayer decoding */ dc1394error_t dc1394_bayer_Simple_uint16(const uint16_t *bayer, uint16_t *rgb, int sx, int sy, int tile, int bits) { uint16_t *outR, *outG, *outB; register int i, j; int tmp, base; /* sx and sy should be even */ switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_BGGR: outR = &rgb[0]; outG = &rgb[1]; outB = &rgb[2]; break; case DC1394_COLOR_FILTER_GBRG: case DC1394_COLOR_FILTER_RGGB: outR = &rgb[2]; outG = &rgb[1]; outB = &rgb[0]; break; default: return DC1394_INVALID_COLOR_FILTER; } switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_BGGR: outR = &rgb[0]; outG = &rgb[1]; outB = &rgb[2]; break; case DC1394_COLOR_FILTER_GBRG: case DC1394_COLOR_FILTER_RGGB: outR = &rgb[2]; outG = &rgb[1]; outB = &rgb[0]; break; default: outR = NULL; outG = NULL; outB = NULL; break; } switch (tile) { case DC1394_COLOR_FILTER_GRBG: case DC1394_COLOR_FILTER_GBRG: for (i = 0; i < sy - 1; i += 2) { for (j = 0; j < sx - 1; j += 2) { base = i * sx + j; tmp = ((bayer[base] + bayer[base + sx + 1]) >> 1); CLIP16(tmp, outG[base * 3], bits); tmp = bayer[base + 1]; CLIP16(tmp, outR[base * 3], bits); tmp = bayer[base + sx]; CLIP16(tmp, outB[base * 3], bits); } } for (i = 0; i < sy - 1; i += 2) { for (j = 1; j < sx - 1; j += 2) { base = i * sx + j; tmp = ((bayer[base + 1] + bayer[base + sx]) >> 1); CLIP16(tmp, outG[(base)*3], bits); tmp = bayer[base]; CLIP16(tmp, outR[(base)*3], bits); tmp = bayer[base + 1 + sx]; CLIP16(tmp, outB[(base)*3], bits); } } for (i = 1; i < sy - 1; i += 2) { for (j = 0; j < sx - 1; j += 2) { base = i * sx + j; tmp = ((bayer[base + sx] + bayer[base + 1]) >> 1); CLIP16(tmp, outG[base * 3], bits); tmp = bayer[base + sx + 1]; CLIP16(tmp, outR[base * 3], bits); tmp = bayer[base]; CLIP16(tmp, outB[base * 3], bits); } } for (i = 1; i < sy - 1; i += 2) { for (j = 1; j < sx - 1; j += 2) { base = i * sx + j; tmp = ((bayer[base] + bayer[base + 1 + sx]) >> 1); CLIP16(tmp, outG[(base)*3], bits); tmp = bayer[base + sx]; CLIP16(tmp, outR[(base)*3], bits); tmp = bayer[base + 1]; CLIP16(tmp, outB[(base)*3], bits); } } break; case DC1394_COLOR_FILTER_BGGR: case DC1394_COLOR_FILTER_RGGB: for (i = 0; i < sy - 1; i += 2) { for (j = 0; j < sx - 1; j += 2) { base = i * sx + j; tmp = ((bayer[base + sx] + bayer[base + 1]) >> 1); CLIP16(tmp, outG[base * 3], bits); tmp = bayer[base + sx + 1]; CLIP16(tmp, outR[base * 3], bits); tmp = bayer[base]; CLIP16(tmp, outB[base * 3], bits); } } for (i = 1; i < sy - 1; i += 2) { for (j = 0; j < sx - 1; j += 2) { base = i * sx + j; tmp = ((bayer[base] + bayer[base + 1 + sx]) >> 1); CLIP16(tmp, outG[(base)*3], bits); tmp = bayer[base + 1]; CLIP16(tmp, outR[(base)*3], bits); tmp = bayer[base + sx]; CLIP16(tmp, outB[(base)*3], bits); } } for (i = 0; i < sy - 1; i += 2) { for (j = 1; j < sx - 1; j += 2) { base = i * sx + j; tmp = ((bayer[base] + bayer[base + sx + 1]) >> 1); CLIP16(tmp, outG[base * 3], bits); tmp = bayer[base + sx]; CLIP16(tmp, outR[base * 3], bits); tmp = bayer[base + 1]; CLIP16(tmp, outB[base * 3], bits); } } for (i = 1; i < sy - 1; i += 2) { for (j = 1; j < sx - 1; j += 2) { base = i * sx + j; tmp = ((bayer[base + 1] + bayer[base + sx]) >> 1); CLIP16(tmp, outG[(base)*3], bits); tmp = bayer[base]; CLIP16(tmp, outR[(base)*3], bits); tmp = bayer[base + 1 + sx]; CLIP16(tmp, outB[(base)*3], bits); } } break; } /* add black border */ for (i = sx * (sy - 1) * 3; i < sx * sy * 3; i++) { rgb[i] = 0; } for (i = (sx - 1) * 3; i < sx * sy * 3; i += (sx - 1) * 3) { rgb[i++] = 0; rgb[i++] = 0; rgb[i++] = 0; } return DC1394_SUCCESS; } /* Variable Number of Gradients, from dcraw */ /* Ported to libdc1394 by Frederic Devernay */ #define FORC3 for (c = 0; c < 3; c++) #define SQR(x) ((x) * (x)) #define ABS(x) (((int)(x) ^ ((int)(x) >> 31)) - ((int)(x) >> 31)) #ifndef MIN #define MIN(a, b) ((a) < (b) ? (a) : (b)) #endif #ifndef MAX #define MAX(a, b) ((a) > (b) ? (a) : (b)) #endif #define LIM(x, min, max) MAX(min, MIN(x, max)) #define ULIM(x, y, z) ((y) < (z) ? LIM(x, y, z) : LIM(x, z, y)) /* In order to inline this calculation, I make the risky assumption that all filter patterns can be described by a repeating pattern of eight rows and two columns Return values are either 0/1/2/3 = G/M/C/Y or 0/1/2/3 = R/G1/B/G2 */ #define FC(row, col) (filters >> ((((row) << 1 & 14) + ((col)&1)) << 1) & 3) /* This algorithm is officially called: "Interpolation using a Threshold-based variable number of gradients" described in http://www-ise.stanford.edu/~tingchen/algodep/vargra.html I've extended the basic idea to work with non-Bayer filter arrays. Gradients are numbered clockwise from NW=0 to W=7. */ static const signed char bayervng_terms[] = { -2, -2, +0, -1, 0, 0x01, -2, -2, +0, +0, 1, 0x01, -2, -1, -1, +0, 0, 0x01, -2, -1, +0, -1, 0, 0x02, -2, -1, +0, +0, 0, 0x03, -2, -1, +0, +1, 1, 0x01, -2, +0, +0, -1, 0, 0x06, -2, +0, +0, +0, 1, 0x02, -2, +0, +0, +1, 0, 0x03, -2, +1, -1, +0, 0, 0x04, -2, +1, +0, -1, 1, 0x04, -2, +1, +0, +0, 0, 0x06, -2, +1, +0, +1, 0, 0x02, -2, +2, +0, +0, 1, 0x04, -2, +2, +0, +1, 0, 0x04, -1, -2, -1, +0, 0, 0x80, -1, -2, +0, -1, 0, 0x01, -1, -2, +1, -1, 0, 0x01, -1, -2, +1, +0, 1, 0x01, -1, -1, -1, +1, 0, 0x88, -1, -1, +1, -2, 0, 0x40, -1, -1, +1, -1, 0, 0x22, -1, -1, +1, +0, 0, 0x33, -1, -1, +1, +1, 1, 0x11, -1, +0, -1, +2, 0, 0x08, -1, +0, +0, -1, 0, 0x44, -1, +0, +0, +1, 0, 0x11, -1, +0, +1, -2, 1, 0x40, -1, +0, +1, -1, 0, 0x66, -1, +0, +1, +0, 1, 0x22, -1, +0, +1, +1, 0, 0x33, -1, +0, +1, +2, 1, 0x10, -1, +1, +1, -1, 1, 0x44, -1, +1, +1, +0, 0, 0x66, -1, +1, +1, +1, 0, 0x22, -1, +1, +1, +2, 0, 0x10, -1, +2, +0, +1, 0, 0x04, -1, +2, +1, +0, 1, 0x04, -1, +2, +1, +1, 0, 0x04, +0, -2, +0, +0, 1, 0x80, +0, -1, +0, +1, 1, 0x88, +0, -1, +1, -2, 0, 0x40, +0, -1, +1, +0, 0, 0x11, +0, -1, +2, -2, 0, 0x40, +0, -1, +2, -1, 0, 0x20, +0, -1, +2, +0, 0, 0x30, +0, -1, +2, +1, 1, 0x10, +0, +0, +0, +2, 1, 0x08, +0, +0, +2, -2, 1, 0x40, +0, +0, +2, -1, 0, 0x60, +0, +0, +2, +0, 1, 0x20, +0, +0, +2, +1, 0, 0x30, +0, +0, +2, +2, 1, 0x10, +0, +1, +1, +0, 0, 0x44, +0, +1, +1, +2, 0, 0x10, +0, +1, +2, -1, 1, 0x40, +0, +1, +2, +0, 0, 0x60, +0, +1, +2, +1, 0, 0x20, +0, +1, +2, +2, 0, 0x10, +1, -2, +1, +0, 0, 0x80, +1, -1, +1, +1, 0, 0x88, +1, +0, +1, +2, 0, 0x08, +1, +0, +2, -1, 0, 0x40, +1, +0, +2, +1, 0, 0x10 }, bayervng_chood[] = { -1, -1, -1, 0, -1, +1, 0, +1, +1, +1, +1, 0, +1, -1, 0, -1 }; dc1394error_t dc1394_bayer_VNG(const uint8_t *bayer, uint8_t *dst, int sx, int sy, dc1394color_filter_t pattern) { const int height = sy, width = sx; static const signed char *cp; /* the following has the same type as the image */ uint8_t(*brow[5])[3], *pix; /* [FD] */ int code[8][2][320], *ip, gval[8], gmin, gmax, sum[4]; int row, col, x, y, x1, x2, y1, y2, t, weight, grads, color, diag; int g, diff, thold, num, c; uint32_t filters; /* [FD] */ /* first, use bilinear bayer decoding */ dc1394_bayer_Bilinear(bayer, dst, sx, sy, pattern); switch (pattern) { case DC1394_COLOR_FILTER_BGGR: filters = 0x16161616; break; case DC1394_COLOR_FILTER_GRBG: filters = 0x61616161; break; case DC1394_COLOR_FILTER_RGGB: filters = 0x94949494; break; case DC1394_COLOR_FILTER_GBRG: filters = 0x49494949; break; default: return DC1394_INVALID_COLOR_FILTER; } for (row = 0; row < 8; row++) { /* Precalculate for VNG */ for (col = 0; col < 2; col++) { ip = code[row][col]; for (cp = bayervng_terms, t = 0; t < 64; t++) { y1 = *cp++; x1 = *cp++; y2 = *cp++; x2 = *cp++; weight = *cp++; grads = *cp++; color = FC(row + y1, col + x1); if ((int)FC(row + y2, col + x2) != color) continue; diag = ((int)FC(row, col + 1) == color && (int)FC(row + 1, col) == color) ? 2 : 1; if (abs(y1 - y2) == diag && abs(x1 - x2) == diag) continue; *ip++ = (y1 * width + x1) * 3 + color; /* [FD] */ *ip++ = (y2 * width + x2) * 3 + color; /* [FD] */ *ip++ = weight; for (g = 0; g < 8; g++) if (grads & 1 << g) *ip++ = g; *ip++ = -1; } *ip++ = INT_MAX; for (cp = bayervng_chood, g = 0; g < 8; g++) { y = *cp++; x = *cp++; *ip++ = (y * width + x) * 3; /* [FD] */ color = FC(row, col); if ((int)FC(row + y, col + x) != color && (int)FC(row + y * 2, col + x * 2) == color) *ip++ = (y * width + x) * 6 + color; /* [FD] */ else *ip++ = 0; } } } brow[4] = calloc(width * 3, sizeof **brow); for (row = 0; row < 3; row++) brow[row] = brow[4] + row * width; for (row = 2; row < height - 2; row++) { /* Do VNG interpolation */ for (col = 2; col < width - 2; col++) { pix = dst + (row * width + col) * 3; /* [FD] */ ip = code[row & 7][col & 1]; memset(gval, 0, sizeof gval); while ((g = ip[0]) != INT_MAX) { /* Calculate gradients */ diff = ABS(pix[g] - pix[ip[1]]) << ip[2]; gval[ip[3]] += diff; ip += 5; if ((g = ip[-1]) == -1) continue; gval[g] += diff; while ((g = *ip++) != -1) gval[g] += diff; } ip++; gmin = gmax = gval[0]; /* Choose a threshold */ for (g = 1; g < 8; g++) { if (gmin > gval[g]) gmin = gval[g]; if (gmax < gval[g]) gmax = gval[g]; } if (gmax == 0) { memcpy(brow[2][col], pix, 3 * sizeof *dst); /* [FD] */ continue; } thold = gmin + (gmax >> 1); memset(sum, 0, sizeof sum); color = FC(row, col); for (num = g = 0; g < 8; g++, ip += 2) { /* Average the neighbors */ if (gval[g] <= thold) { for (c = 0; c < 3; c++) /* [FD] */ if (c == color && ip[1]) sum[c] += (pix[c] + pix[ip[1]]) >> 1; else sum[c] += pix[ip[0] + c]; num++; } } for (c = 0; c < 3; c++) { /* [FD] Save to buffer */ t = pix[color]; if (c != color) t += (sum[c] - sum[color]) / num; CLIP(t, brow[2][col][c]); /* [FD] */ } } if (row > 3) /* Write buffer to image */ memcpy(dst + 3 * ((row - 2) * width + 2), brow[0] + 2, (width - 4) * 3 * sizeof *dst); /* [FD] */ for (g = 0; g < 4; g++) brow[(g - 1) & 3] = brow[g]; } memcpy(dst + 3 * ((row - 2) * width + 2), brow[0] + 2, (width - 4) * 3 * sizeof *dst); memcpy(dst + 3 * ((row - 1) * width + 2), brow[1] + 2, (width - 4) * 3 * sizeof *dst); free(brow[4]); return DC1394_SUCCESS; } dc1394error_t dc1394_bayer_VNG_uint16(const uint16_t *bayer, uint16_t *dst, int sx, int sy, dc1394color_filter_t pattern, int bits) { const int height = sy, width = sx; static const signed char *cp; /* the following has the same type as the image */ uint16_t(*brow[5])[3], *pix; /* [FD] */ int code[8][2][320], *ip, gval[8], gmin, gmax, sum[4]; int row, col, x, y, x1, x2, y1, y2, t, weight, grads, color, diag; int g, diff, thold, num, c; uint32_t filters; /* [FD] */ /* first, use bilinear bayer decoding */ dc1394_bayer_Bilinear_uint16(bayer, dst, sx, sy, pattern, bits); switch (pattern) { case DC1394_COLOR_FILTER_BGGR: filters = 0x16161616; break; case DC1394_COLOR_FILTER_GRBG: filters = 0x61616161; break; case DC1394_COLOR_FILTER_RGGB: filters = 0x94949494; break; case DC1394_COLOR_FILTER_GBRG: filters = 0x49494949; break; default: return DC1394_INVALID_COLOR_FILTER; } for (row = 0; row < 8; row++) { /* Precalculate for VNG */ for (col = 0; col < 2; col++) { ip = code[row][col]; for (cp = bayervng_terms, t = 0; t < 64; t++) { y1 = *cp++; x1 = *cp++; y2 = *cp++; x2 = *cp++; weight = *cp++; grads = *cp++; color = FC(row + y1, col + x1); if ((int)FC(row + y2, col + x2) != color) continue; diag = ((int)FC(row, col + 1) == color && (int)FC(row + 1, col) == color) ? 2 : 1; if (abs(y1 - y2) == diag && abs(x1 - x2) == diag) continue; *ip++ = (y1 * width + x1) * 3 + color; /* [FD] */ *ip++ = (y2 * width + x2) * 3 + color; /* [FD] */ *ip++ = weight; for (g = 0; g < 8; g++) if (grads & 1 << g) *ip++ = g; *ip++ = -1; } *ip++ = INT_MAX; for (cp = bayervng_chood, g = 0; g < 8; g++) { y = *cp++; x = *cp++; *ip++ = (y * width + x) * 3; /* [FD] */ color = FC(row, col); if ((int)FC(row + y, col + x) != color && (int)FC(row + y * 2, col + x * 2) == color) *ip++ = (y * width + x) * 6 + color; /* [FD] */ else *ip++ = 0; } } } brow[4] = calloc(width * 3, sizeof **brow); for (row = 0; row < 3; row++) brow[row] = brow[4] + row * width; for (row = 2; row < height - 2; row++) { /* Do VNG interpolation */ for (col = 2; col < width - 2; col++) { pix = dst + (row * width + col) * 3; /* [FD] */ ip = code[row & 7][col & 1]; memset(gval, 0, sizeof gval); while ((g = ip[0]) != INT_MAX) { /* Calculate gradients */ diff = ABS(pix[g] - pix[ip[1]]) << ip[2]; gval[ip[3]] += diff; ip += 5; if ((g = ip[-1]) == -1) continue; gval[g] += diff; while ((g = *ip++) != -1) gval[g] += diff; } ip++; gmin = gmax = gval[0]; /* Choose a threshold */ for (g = 1; g < 8; g++) { if (gmin > gval[g]) gmin = gval[g]; if (gmax < gval[g]) gmax = gval[g]; } if (gmax == 0) { memcpy(brow[2][col], pix, 3 * sizeof *dst); /* [FD] */ continue; } thold = gmin + (gmax >> 1); memset(sum, 0, sizeof sum); color = FC(row, col); for (num = g = 0; g < 8; g++, ip += 2) { /* Average the neighbors */ if (gval[g] <= thold) { for (c = 0; c < 3; c++) /* [FD] */ if (c == color && ip[1]) sum[c] += (pix[c] + pix[ip[1]]) >> 1; else sum[c] += pix[ip[0] + c]; num++; } } for (c = 0; c < 3; c++) { /* [FD] Save to buffer */ t = pix[color]; if (c != color) t += (sum[c] - sum[color]) / num; CLIP16(t, brow[2][col][c], bits); /* [FD] */ } } if (row > 3) /* Write buffer to image */ memcpy(dst + 3 * ((row - 2) * width + 2), brow[0] + 2, (width - 4) * 3 * sizeof *dst); /* [FD] */ for (g = 0; g < 4; g++) brow[(g - 1) & 3] = brow[g]; } memcpy(dst + 3 * ((row - 2) * width + 2), brow[0] + 2, (width - 4) * 3 * sizeof *dst); memcpy(dst + 3 * ((row - 1) * width + 2), brow[1] + 2, (width - 4) * 3 * sizeof *dst); free(brow[4]); return DC1394_SUCCESS; } /* AHD interpolation ported from dcraw to libdc1394 by Samuel Audet */ static dc1394bool_t ahd_inited = DC1394_FALSE; /* WARNING: not multi-processor safe */ #define CLIPOUT(x) LIM(x, 0, 255) #define CLIPOUT16(x, bits) LIM(x, 0, ((1 << bits) - 1)) static const double xyz_rgb[3][3] = { /* XYZ from RGB */ { 0.412453, 0.357580, 0.180423 }, { 0.212671, 0.715160, 0.072169 }, { 0.019334, 0.119193, 0.950227 } }; static const float d65_white[3] = { 0.950456, 1, 1.088754 }; static void cam_to_cielab(uint16_t cam[3], float lab[3]) /* [SA] */ { int c, i, j; float r, xyz[3]; static float cbrt[0x10000], xyz_cam[3][4]; if (cam == NULL) { for (i = 0; i < 0x10000; i++) { r = i / 65535.0; cbrt[i] = r > 0.008856 ? pow(r, 1 / 3.0) : 7.787 * r + 16 / 116.0; } for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) /* [SA] */ xyz_cam[i][j] = xyz_rgb[i][j] / d65_white[i]; /* [SA] */ } else { xyz[0] = xyz[1] = xyz[2] = 0.5; FORC3 { /* [SA] */ xyz[0] += xyz_cam[0][c] * cam[c]; xyz[1] += xyz_cam[1][c] * cam[c]; xyz[2] += xyz_cam[2][c] * cam[c]; } xyz[0] = cbrt[CLIPOUT16((int)xyz[0], 16)]; /* [SA] */ xyz[1] = cbrt[CLIPOUT16((int)xyz[1], 16)]; /* [SA] */ xyz[2] = cbrt[CLIPOUT16((int)xyz[2], 16)]; /* [SA] */ lab[0] = 116 * xyz[1] - 16; lab[1] = 500 * (xyz[0] - xyz[1]); lab[2] = 200 * (xyz[1] - xyz[2]); } } /* Adaptive Homogeneity-Directed interpolation is based on the work of Keigo Hirakawa, Thomas Parks, and Paul Lee. */ #define TS 256 /* Tile Size */ dc1394error_t dc1394_bayer_AHD(const uint8_t *bayer, uint8_t *dst, int sx, int sy, dc1394color_filter_t pattern) { int i, j, top, left, row, col, tr, tc, fc, c, d, val, hm[2]; /* the following has the same type as the image */ uint8_t(*pix)[3], (*rix)[3]; /* [SA] */ uint16_t rix16[3]; /* [SA] */ static const int dir[4] = { -1, 1, -TS, TS }; unsigned ldiff[2][4], abdiff[2][4], leps, abeps; float flab[3]; /* [SA] */ uint8_t(*rgb)[TS][TS][3]; short(*lab)[TS][TS][3]; char(*homo)[TS][TS], *buffer; /* start - new code for libdc1394 */ uint32_t filters; const int height = sy, width = sx; int x, y; if (ahd_inited == DC1394_FALSE) { /* WARNING: this might not be multi-processor safe */ cam_to_cielab(NULL, NULL); ahd_inited = DC1394_TRUE; } switch (pattern) { case DC1394_COLOR_FILTER_BGGR: filters = 0x16161616; break; case DC1394_COLOR_FILTER_GRBG: filters = 0x61616161; break; case DC1394_COLOR_FILTER_RGGB: filters = 0x94949494; break; case DC1394_COLOR_FILTER_GBRG: filters = 0x49494949; break; default: return DC1394_INVALID_COLOR_FILTER; } /* fill-in destination with known exact values */ for (y = 0; y < height; y++) { for (x = 0; x < width; x++) { int channel = FC(y, x); dst[(y * width + x) * 3 + channel] = bayer[y * width + x]; } } /* end - new code for libdc1394 */ /* start - code from border_interpolate (int border) */ { int border = 3; int row, col, y, x; unsigned f, c, sum[8]; for (row = 0; row < height; row++) for (col = 0; col < width; col++) { if (col == border && row >= border && row < height - border) col = width - border; memset(sum, 0, sizeof sum); for (y = row - 1; y != row + 2; y++) for (x = col - 1; x != col + 2; x++) if (y < height && x < width) { f = FC(y, x); sum[f] += dst[(y * width + x) * 3 + f]; /* [SA] */ sum[f + 4]++; } f = FC(row, col); FORC3 if (c != f && sum[c + 4]) /* [SA] */ dst[(row * width + col) * 3 + c] = sum[c] / sum[c + 4]; /* [SA] */ } } /* end - code from border_interpolate (int border) */ buffer = (char *)malloc(26 * TS * TS); /* 1664 kB */ /* merror (buffer, "ahd_interpolate()"); */ rgb = (uint8_t(*)[TS][TS][3])buffer; /* [SA] */ lab = (short(*)[TS][TS][3])(buffer + 12 * TS * TS); homo = (char(*)[TS][TS])(buffer + 24 * TS * TS); for (top = 0; top < height; top += TS - 6) for (left = 0; left < width; left += TS - 6) { memset(rgb, 0, 12 * TS * TS); /* Interpolate green horizontally and vertically: */ for (row = top < 2 ? 2 : top; row < top + TS && row < height - 2; row++) { col = left + (FC(row, left) == 1); if (col < 2) col += 2; for (fc = FC(row, col); col < left + TS && col < width - 2; col += 2) { pix = (uint8_t(*)[3])dst + (row * width + col); /* [SA] */ val = ((pix[-1][1] + pix[0][fc] + pix[1][1]) * 2 - pix[-2][fc] - pix[2][fc]) >> 2; rgb[0][row - top][col - left][1] = ULIM(val, pix[-1][1], pix[1][1]); val = ((pix[-width][1] + pix[0][fc] + pix[width][1]) * 2 - pix[-2 * width][fc] - pix[2 * width][fc]) >> 2; rgb[1][row - top][col - left][1] = ULIM(val, pix[-width][1], pix[width][1]); } } /* Interpolate red and blue, and convert to CIELab: */ for (d = 0; d < 2; d++) for (row = top + 1; row < top + TS - 1 && row < height - 1; row++) for (col = left + 1; col < left + TS - 1 && col < width - 1; col++) { pix = (uint8_t(*)[3])dst + (row * width + col); /* [SA] */ rix = &rgb[d][row - top][col - left]; if ((c = 2 - FC(row, col)) == 1) { c = FC(row + 1, col); val = pix[0][1] + ((pix[-1][2 - c] + pix[1][2 - c] - rix[-1][1] - rix[1][1]) >> 1); rix[0][2 - c] = CLIPOUT(val); /* [SA] */ val = pix[0][1] + ((pix[-width][c] + pix[width][c] - rix[-TS][1] - rix[TS][1]) >> 1); } else val = rix[0][1] + ((pix[-width - 1][c] + pix[-width + 1][c] + pix[+width - 1][c] + pix[+width + 1][c] - rix[-TS - 1][1] - rix[-TS + 1][1] - rix[+TS - 1][1] - rix[+TS + 1][1] + 1) >> 2); rix[0][c] = CLIPOUT(val); /* [SA] */ c = FC(row, col); rix[0][c] = pix[0][c]; rix16[0] = rix[0][0]; /* [SA] */ rix16[1] = rix[0][1]; /* [SA] */ rix16[2] = rix[0][2]; /* [SA] */ cam_to_cielab(rix16, flab); /* [SA] */ FORC3 lab[d][row - top][col - left][c] = 64 * flab[c]; } /* Build homogeneity maps from the CIELab images: */ memset(homo, 0, 2 * TS * TS); for (row = top + 2; row < top + TS - 2 && row < height; row++) { tr = row - top; for (col = left + 2; col < left + TS - 2 && col < width; col++) { tc = col - left; for (d = 0; d < 2; d++) for (i = 0; i < 4; i++) ldiff[d][i] = ABS(lab[d][tr][tc][0] - lab[d][tr][tc + dir[i]][0]); leps = MIN(MAX(ldiff[0][0], ldiff[0][1]), MAX(ldiff[1][2], ldiff[1][3])); for (d = 0; d < 2; d++) for (i = 0; i < 4; i++) if (i >> 1 == d || ldiff[d][i] <= leps) abdiff[d][i] = SQR(lab[d][tr][tc][1] - lab[d][tr][tc + dir[i]][1]) + SQR(lab[d][tr][tc][2] - lab[d][tr][tc + dir[i]][2]); abeps = MIN(MAX(abdiff[0][0], abdiff[0][1]), MAX(abdiff[1][2], abdiff[1][3])); for (d = 0; d < 2; d++) for (i = 0; i < 4; i++) if (ldiff[d][i] <= leps && abdiff[d][i] <= abeps) homo[d][tr][tc]++; } } /* Combine the most homogenous pixels for the final result: */ for (row = top + 3; row < top + TS - 3 && row < height - 3; row++) { tr = row - top; for (col = left + 3; col < left + TS - 3 && col < width - 3; col++) { tc = col - left; for (d = 0; d < 2; d++) for (hm[d] = 0, i = tr - 1; i <= tr + 1; i++) for (j = tc - 1; j <= tc + 1; j++) hm[d] += homo[d][i][j]; if (hm[0] != hm[1]) FORC3 dst[(row * width + col) * 3 + c] = CLIPOUT(rgb[hm[1] > hm[0]][tr][tc][c]); /* [SA] */ else FORC3 dst[(row * width + col) * 3 + c] = CLIPOUT((rgb[0][tr][tc][c] + rgb[1][tr][tc][c]) >> 1); /* [SA] */ } } } free(buffer); return DC1394_SUCCESS; } dc1394error_t dc1394_bayer_AHD_uint16(const uint16_t *bayer, uint16_t *dst, int sx, int sy, dc1394color_filter_t pattern, int bits) { int i, j, top, left, row, col, tr, tc, fc, c, d, val, hm[2]; /* the following has the same type as the image */ uint16_t(*pix)[3], (*rix)[3]; /* [SA] */ static const int dir[4] = { -1, 1, -TS, TS }; unsigned ldiff[2][4], abdiff[2][4], leps, abeps; float flab[3]; uint16_t(*rgb)[TS][TS][3]; /* [SA] */ short(*lab)[TS][TS][3]; char(*homo)[TS][TS], *buffer; /* start - new code for libdc1394 */ uint32_t filters; const int height = sy, width = sx; int x, y; if (ahd_inited == DC1394_FALSE) { /* WARNING: this might not be multi-processor safe */ cam_to_cielab(NULL, NULL); ahd_inited = DC1394_TRUE; } switch (pattern) { case DC1394_COLOR_FILTER_BGGR: filters = 0x16161616; break; case DC1394_COLOR_FILTER_GRBG: filters = 0x61616161; break; case DC1394_COLOR_FILTER_RGGB: filters = 0x94949494; break; case DC1394_COLOR_FILTER_GBRG: filters = 0x49494949; break; default: return DC1394_INVALID_COLOR_FILTER; } /* fill-in destination with known exact values */ for (y = 0; y < height; y++) { for (x = 0; x < width; x++) { int channel = FC(y, x); dst[(y * width + x) * 3 + channel] = bayer[y * width + x]; } } /* end - new code for libdc1394 */ /* start - code from border_interpolate(int border) */ { int border = 3; int row, col, y, x; unsigned f, c, sum[8]; for (row = 0; row < height; row++) for (col = 0; col < width; col++) { if (col == border && row >= border && row < height - border) col = width - border; memset(sum, 0, sizeof sum); for (y = row - 1; y != row + 2; y++) for (x = col - 1; x != col + 2; x++) if (y < height && x < width) { f = FC(y, x); sum[f] += dst[(y * width + x) * 3 + f]; /* [SA] */ sum[f + 4]++; } f = FC(row, col); FORC3 if (c != f && sum[c + 4]) /* [SA] */ dst[(row * width + col) * 3 + c] = sum[c] / sum[c + 4]; /* [SA] */ } } /* end - code from border_interpolate(int border) */ buffer = (char *)malloc(26 * TS * TS); /* 1664 kB */ /* merror (buffer, "ahd_interpolate()"); */ rgb = (uint16_t(*)[TS][TS][3])buffer; /* [SA] */ lab = (short(*)[TS][TS][3])(buffer + 12 * TS * TS); homo = (char(*)[TS][TS])(buffer + 24 * TS * TS); for (top = 0; top < height; top += TS - 6) for (left = 0; left < width; left += TS - 6) { memset(rgb, 0, 12 * TS * TS); /* Interpolate green horizontally and vertically: */ for (row = top < 2 ? 2 : top; row < top + TS && row < height - 2; row++) { col = left + (FC(row, left) == 1); if (col < 2) col += 2; for (fc = FC(row, col); col < left + TS && col < width - 2; col += 2) { pix = (uint16_t(*)[3])dst + (row * width + col); /* [SA] */ val = ((pix[-1][1] + pix[0][fc] + pix[1][1]) * 2 - pix[-2][fc] - pix[2][fc]) >> 2; rgb[0][row - top][col - left][1] = ULIM(val, pix[-1][1], pix[1][1]); val = ((pix[-width][1] + pix[0][fc] + pix[width][1]) * 2 - pix[-2 * width][fc] - pix[2 * width][fc]) >> 2; rgb[1][row - top][col - left][1] = ULIM(val, pix[-width][1], pix[width][1]); } } /* Interpolate red and blue, and convert to CIELab: */ for (d = 0; d < 2; d++) for (row = top + 1; row < top + TS - 1 && row < height - 1; row++) for (col = left + 1; col < left + TS - 1 && col < width - 1; col++) { pix = (uint16_t(*)[3])dst + (row * width + col); /* [SA] */ rix = &rgb[d][row - top][col - left]; if ((c = 2 - FC(row, col)) == 1) { c = FC(row + 1, col); val = pix[0][1] + ((pix[-1][2 - c] + pix[1][2 - c] - rix[-1][1] - rix[1][1]) >> 1); rix[0][2 - c] = CLIPOUT16(val, bits); /* [SA] */ val = pix[0][1] + ((pix[-width][c] + pix[width][c] - rix[-TS][1] - rix[TS][1]) >> 1); } else val = rix[0][1] + ((pix[-width - 1][c] + pix[-width + 1][c] + pix[+width - 1][c] + pix[+width + 1][c] - rix[-TS - 1][1] - rix[-TS + 1][1] - rix[+TS - 1][1] - rix[+TS + 1][1] + 1) >> 2); rix[0][c] = CLIPOUT16(val, bits); /* [SA] */ c = FC(row, col); rix[0][c] = pix[0][c]; cam_to_cielab(rix[0], flab); FORC3 lab[d][row - top][col - left][c] = 64 * flab[c]; } /* Build homogeneity maps from the CIELab images: */ memset(homo, 0, 2 * TS * TS); for (row = top + 2; row < top + TS - 2 && row < height; row++) { tr = row - top; for (col = left + 2; col < left + TS - 2 && col < width; col++) { tc = col - left; for (d = 0; d < 2; d++) for (i = 0; i < 4; i++) ldiff[d][i] = ABS(lab[d][tr][tc][0] - lab[d][tr][tc + dir[i]][0]); leps = MIN(MAX(ldiff[0][0], ldiff[0][1]), MAX(ldiff[1][2], ldiff[1][3])); for (d = 0; d < 2; d++) for (i = 0; i < 4; i++) if (i >> 1 == d || ldiff[d][i] <= leps) abdiff[d][i] = SQR(lab[d][tr][tc][1] - lab[d][tr][tc + dir[i]][1]) + SQR(lab[d][tr][tc][2] - lab[d][tr][tc + dir[i]][2]); abeps = MIN(MAX(abdiff[0][0], abdiff[0][1]), MAX(abdiff[1][2], abdiff[1][3])); for (d = 0; d < 2; d++) for (i = 0; i < 4; i++) if (ldiff[d][i] <= leps && abdiff[d][i] <= abeps) homo[d][tr][tc]++; } } /* Combine the most homogenous pixels for the final result: */ for (row = top + 3; row < top + TS - 3 && row < height - 3; row++) { tr = row - top; for (col = left + 3; col < left + TS - 3 && col < width - 3; col++) { tc = col - left; for (d = 0; d < 2; d++) for (hm[d] = 0, i = tr - 1; i <= tr + 1; i++) for (j = tc - 1; j <= tc + 1; j++) hm[d] += homo[d][i][j]; if (hm[0] != hm[1]) FORC3 dst[(row * width + col) * 3 + c] = CLIPOUT16(rgb[hm[1] > hm[0]][tr][tc][c], bits); /* [SA] */ else FORC3 dst[(row * width + col) * 3 + c] = CLIPOUT16((rgb[0][tr][tc][c] + rgb[1][tr][tc][c]) >> 1, bits); /* [SA] */ } } } free(buffer); return DC1394_SUCCESS; } dc1394error_t dc1394_bayer_decoding_8bit(const uint8_t *bayer, uint8_t *rgb, uint32_t sx, uint32_t sy, dc1394color_filter_t tile, dc1394bayer_method_t method) { switch (method) { case DC1394_BAYER_METHOD_NEAREST: return dc1394_bayer_NearestNeighbor(bayer, rgb, sx, sy, tile); case DC1394_BAYER_METHOD_SIMPLE: return dc1394_bayer_Simple(bayer, rgb, sx, sy, tile); case DC1394_BAYER_METHOD_BILINEAR: return dc1394_bayer_Bilinear(bayer, rgb, sx, sy, tile); case DC1394_BAYER_METHOD_HQLINEAR: return dc1394_bayer_HQLinear(bayer, rgb, sx, sy, tile); case DC1394_BAYER_METHOD_DOWNSAMPLE: return dc1394_bayer_Downsample(bayer, rgb, sx, sy, tile); case DC1394_BAYER_METHOD_EDGESENSE: return dc1394_bayer_EdgeSense(bayer, rgb, sx, sy, tile); case DC1394_BAYER_METHOD_VNG: return dc1394_bayer_VNG(bayer, rgb, sx, sy, tile); case DC1394_BAYER_METHOD_AHD: return dc1394_bayer_AHD(bayer, rgb, sx, sy, tile); default: return DC1394_INVALID_BAYER_METHOD; } } dc1394error_t dc1394_bayer_decoding_16bit(const uint16_t *bayer, uint16_t *rgb, uint32_t sx, uint32_t sy, dc1394color_filter_t tile, dc1394bayer_method_t method, uint32_t bits) { switch (method) { case DC1394_BAYER_METHOD_NEAREST: return dc1394_bayer_NearestNeighbor_uint16(bayer, rgb, sx, sy, tile, bits); case DC1394_BAYER_METHOD_SIMPLE: return dc1394_bayer_Simple_uint16(bayer, rgb, sx, sy, tile, bits); case DC1394_BAYER_METHOD_BILINEAR: return dc1394_bayer_Bilinear_uint16(bayer, rgb, sx, sy, tile, bits); case DC1394_BAYER_METHOD_HQLINEAR: return dc1394_bayer_HQLinear_uint16(bayer, rgb, sx, sy, tile, bits); case DC1394_BAYER_METHOD_DOWNSAMPLE: return dc1394_bayer_Downsample_uint16(bayer, rgb, sx, sy, tile, bits); case DC1394_BAYER_METHOD_EDGESENSE: return dc1394_bayer_EdgeSense_uint16(bayer, rgb, sx, sy, tile, bits); case DC1394_BAYER_METHOD_VNG: return dc1394_bayer_VNG_uint16(bayer, rgb, sx, sy, tile, bits); case DC1394_BAYER_METHOD_AHD: return dc1394_bayer_AHD_uint16(bayer, rgb, sx, sy, tile, bits); default: return DC1394_INVALID_BAYER_METHOD; } } diff --git a/kstars/fitsviewer/fitsdata.cpp b/kstars/fitsviewer/fitsdata.cpp index fb0978445..9fb0ab341 100644 --- a/kstars/fitsviewer/fitsdata.cpp +++ b/kstars/fitsviewer/fitsdata.cpp @@ -1,4874 +1,4255 @@ /*************************************************************************** FITSImage.cpp - FITS Image ------------------- begin : Thu Jan 22 2004 copyright : (C) 2004 by Jasem Mutlaq email : mutlaqja@ikarustech.com ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* * See http://members.aol.com/pkirchg for more details. * ***************************************************************************/ #include "fitsdata.h" #include "sep/sep.h" #include "fpack.h" #include "kstarsdata.h" #include "ksutils.h" #include "kspaths.h" #include "Options.h" #include "skymapcomposite.h" #include "auxiliary/ksnotification.h" #include #include #include #include #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB) #include #include #endif #ifndef KSTARS_LITE #include "fitshistogram.h" #endif #include #include #include #define ZOOM_DEFAULT 100.0 #define ZOOM_MIN 10 #define ZOOM_MAX 400 #define ZOOM_LOW_INCR 10 #define ZOOM_HIGH_INCR 50 const int MINIMUM_ROWS_PER_CENTER = 3; const QString FITSData::m_TemporaryPath = QStandardPaths::writableLocation(QStandardPaths::TempLocation); #define DIFFUSE_THRESHOLD 0.15 #define MAX_EDGE_LIMIT 10000 #define LOW_EDGE_CUTOFF_1 50 #define LOW_EDGE_CUTOFF_2 10 #define MINIMUM_EDGE_LIMIT 2 bool greaterThan(Edge * s1, Edge * s2) { //return s1->width > s2->width; return s1->sum > s2->sum; } FITSData::FITSData(FITSMode fitsMode): m_Mode(fitsMode) { debayerParams.method = DC1394_BAYER_METHOD_NEAREST; debayerParams.filter = DC1394_COLOR_FILTER_RGGB; debayerParams.offsetX = debayerParams.offsetY = 0; } FITSData::FITSData(const FITSData * other) { debayerParams.method = DC1394_BAYER_METHOD_NEAREST; debayerParams.filter = DC1394_COLOR_FILTER_RGGB; debayerParams.offsetX = debayerParams.offsetY = 0; this->m_Mode = other->m_Mode; this->m_DataType = other->m_DataType; this->m_Channels = other->m_Channels; memcpy(&stats, &(other->stats), sizeof(stats)); m_ImageBuffer = new uint8_t[stats.samples_per_channel * m_Channels * stats.bytesPerPixel]; memcpy(m_ImageBuffer, other->m_ImageBuffer, stats.samples_per_channel * m_Channels * stats.bytesPerPixel); } FITSData::~FITSData() { int status = 0; clearImageBuffers(); if (starCenters.count() > 0) qDeleteAll(starCenters); delete[] wcs_coord; if (objList.count() > 0) qDeleteAll(objList); if (fptr != nullptr) { fits_close_file(fptr, &status); fptr = nullptr; if (m_isTemporary && autoRemoveTemporaryFITS) QFile::remove(m_Filename); } qDeleteAll(records); } QFuture FITSData::loadFITS(const QString &inFilename, bool silent) { int status = 0; qDeleteAll(starCenters); starCenters.clear(); if (fptr != nullptr) { fits_close_file(fptr, &status); fptr = nullptr; // If current file is temporary AND // Auto Remove Temporary File is Set AND // New filename is different from existing filename // THen remove it. We have to check for name since we cannot delete // the same filename and try to open it below! if (m_isTemporary && autoRemoveTemporaryFITS && inFilename != m_Filename) QFile::remove(m_Filename); } m_Filename = inFilename; qCInfo(KSTARS_FITS) << "Loading FITS file " << m_Filename; QFuture result = QtConcurrent::run(this, &FITSData::privateLoad, silent); return result; } bool FITSData::privateLoad(bool silent) { int status = 0, anynull = 0; long naxes[3]; char error_status[512]; QString errMessage; if (m_Filename.startsWith(m_TemporaryPath)) m_isTemporary = true; if (m_Filename.endsWith(".fz")) { // Store so we don't lose. m_compressedFilename = m_Filename; QString uncompressedFile = QDir::tempPath() + QString("/%1").arg(QUuid::createUuid().toString().remove(QRegularExpression("[-{}]"))); fpstate fpvar; fp_init (&fpvar); if (fp_unpack(m_Filename.toLatin1().data(), uncompressedFile.toLatin1().data(), fpvar) < 0) { errMessage = i18n("Failed to unpack compressed fits"); if (!silent) KSNotification::error(errMessage, i18n("FITS Open")); qCCritical(KSTARS_FITS) << errMessage; return false; } // Remove compressed .fz if it was temporary if (m_isTemporary && autoRemoveTemporaryFITS) QFile::remove(m_Filename); m_Filename = uncompressedFile; m_isTemporary = true; m_isCompressed = true; } // Use open diskfile as it does not use extended file names which has problems opening // files with [ ] or ( ) in their names. if (fits_open_diskfile(&fptr, m_Filename.toLatin1(), READONLY, &status)) { fits_report_error(stderr, status); fits_get_errstatus(status, error_status); errMessage = i18n("Could not open file %1. Error %2", m_Filename, QString::fromUtf8(error_status)); if (!silent) KSNotification::error(errMessage, i18n("FITS Open")); qCCritical(KSTARS_FITS) << errMessage; return false; } stats.size = QFile(m_Filename).size(); if (fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status)) { fits_report_error(stderr, status); fits_get_errstatus(status, error_status); errMessage = i18n("Could not locate image HDU. Error %1", QString::fromUtf8(error_status)); if (!silent) KSNotification::error(errMessage, i18n("FITS Open")); qCCritical(KSTARS_FITS) << errMessage; return false; } if (fits_get_img_param(fptr, 3, &(stats.bitpix), &(stats.ndim), naxes, &status)) { fits_report_error(stderr, status); fits_get_errstatus(status, error_status); errMessage = i18n("FITS file open error (fits_get_img_param): %1", QString::fromUtf8(error_status)); if (!silent) KSNotification::error(errMessage, i18n("FITS Open")); qCCritical(KSTARS_FITS) << errMessage; return false; } if (stats.ndim < 2) { errMessage = i18n("1D FITS images are not supported in KStars."); if (!silent) KSNotification::error(errMessage, i18n("FITS Open")); qCCritical(KSTARS_FITS) << errMessage; return false; } switch (stats.bitpix) { case BYTE_IMG: m_DataType = TBYTE; stats.bytesPerPixel = sizeof(uint8_t); break; case SHORT_IMG: // Read SHORT image as USHORT m_DataType = TUSHORT; stats.bytesPerPixel = sizeof(int16_t); break; case USHORT_IMG: m_DataType = TUSHORT; stats.bytesPerPixel = sizeof(uint16_t); break; case LONG_IMG: // Read LONG image as ULONG m_DataType = TULONG; stats.bytesPerPixel = sizeof(int32_t); break; case ULONG_IMG: m_DataType = TULONG; stats.bytesPerPixel = sizeof(uint32_t); break; case FLOAT_IMG: m_DataType = TFLOAT; stats.bytesPerPixel = sizeof(float); break; case LONGLONG_IMG: m_DataType = TLONGLONG; stats.bytesPerPixel = sizeof(int64_t); break; case DOUBLE_IMG: m_DataType = TDOUBLE; stats.bytesPerPixel = sizeof(double); break; default: errMessage = i18n("Bit depth %1 is not supported.", stats.bitpix); if (!silent) KSNotification::error(errMessage, i18n("FITS Open")); qCCritical(KSTARS_FITS) << errMessage; return false; } if (stats.ndim < 3) naxes[2] = 1; if (naxes[0] == 0 || naxes[1] == 0) { errMessage = i18n("Image has invalid dimensions %1x%2", naxes[0], naxes[1]); if (!silent) KSNotification::error(errMessage, i18n("FITS Open")); qCCritical(KSTARS_FITS) << errMessage; return false; } stats.width = naxes[0]; stats.height = naxes[1]; stats.samples_per_channel = stats.width * stats.height; clearImageBuffers(); m_Channels = naxes[2]; // Channels always set to #1 if we are not required to process 3D Cubes // Or if mode is not FITS_NORMAL (guide, focus..etc) if (m_Mode != FITS_NORMAL || !Options::auto3DCube()) m_Channels = 1; //image_buffer = new float[stats.samples_per_channel * channels]; m_ImageBuffer = new uint8_t[stats.samples_per_channel * m_Channels * stats.bytesPerPixel]; //if (image_buffer == nullptr) if (m_ImageBuffer == nullptr) { qCWarning(KSTARS_FITS) << "FITSData: Not enough memory for image_buffer channel. Requested: " << stats.samples_per_channel * m_Channels * stats.bytesPerPixel << " bytes."; clearImageBuffers(); return false; } rotCounter = 0; flipHCounter = 0; flipVCounter = 0; long nelements = stats.samples_per_channel * m_Channels; if (fits_read_img(fptr, m_DataType, 1, nelements, nullptr, m_ImageBuffer, &anynull, &status)) { char errmsg[512]; fits_get_errstatus(status, errmsg); errMessage = i18n("Error reading image: %1", QString(errmsg)); if (!silent) KSNotification::error(errMessage, i18n("FITS Open")); fits_report_error(stderr, status); qCCritical(KSTARS_FITS) << errMessage; return false; } parseHeader(); if (Options::autoDebayer() && checkDebayer()) { bayerBuffer = m_ImageBuffer; if (debayer()) calculateStats(); } else calculateStats(); WCSLoaded = false; if (m_Mode == FITS_NORMAL || m_Mode == FITS_ALIGN) checkForWCS(); starsSearched = false; return true; } int FITSData::saveFITS(const QString &newFilename) { if (newFilename == m_Filename) return 0; if (m_isCompressed) { KSNotification::error(i18n("Saving compressed files is not supported.")); return -1; } int status = 0, exttype = 0; long nelements; fitsfile * new_fptr; if (HasDebayer) { /* close current file */ if (fits_close_file(fptr, &status)) { fits_report_error(stderr, status); return status; } // Skip "!" in the beginning of the new file name QString finalFileName(newFilename); finalFileName.remove('!'); // Remove first otherwise copy will fail below if file exists QFile::remove(finalFileName); if (!QFile::copy(m_Filename, finalFileName)) { qCCritical(KSTARS_FITS()) << "FITS: Failed to copy " << m_Filename << " to " << finalFileName; fptr = nullptr; return -1; } if (m_isTemporary && autoRemoveTemporaryFITS) { QFile::remove(m_Filename); m_isTemporary = false; } m_Filename = finalFileName; - - //fits_open_image(&fptr, filename.toLatin1(), READONLY, &status); - // Use open diskfile as it does not use extended file names which has problems opening // files with [ ] or ( ) in their names. fits_open_diskfile(&fptr, m_Filename.toLatin1(), READONLY, &status); fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status); return 0; } nelements = stats.samples_per_channel * m_Channels; /* Create a new File, overwriting existing*/ if (fits_create_file(&new_fptr, newFilename.toLatin1(), &status)) { fits_report_error(stderr, status); return status; } if (fits_movabs_hdu(fptr, 1, &exttype, &status)) { fits_report_error(stderr, status); return status; } - /*if (fits_copy_file(fptr, new_fptr, 0, 1, 1, &status)) - { - fits_report_error(stderr, status); - return status; - }*/ - if (fits_copy_header(fptr, new_fptr, &status)) { fits_report_error(stderr, status); return status; } /* close current file */ if (fits_close_file(fptr, &status)) { fits_report_error(stderr, status); return status; } status = 0; fptr = new_fptr; if (fits_movabs_hdu(fptr, 1, &exttype, &status)) { fits_report_error(stderr, status); return status; } /* Write Data */ if (fits_write_img(fptr, m_DataType, 1, nelements, m_ImageBuffer, &status)) { fits_report_error(stderr, status); return status; } /* Write keywords */ // Minimum if (fits_update_key(fptr, TDOUBLE, "DATAMIN", &(stats.min), "Minimum value", &status)) { fits_report_error(stderr, status); return status; } // Maximum if (fits_update_key(fptr, TDOUBLE, "DATAMAX", &(stats.max), "Maximum value", &status)) { fits_report_error(stderr, status); return status; } // NAXIS1 if (fits_update_key(fptr, TUSHORT, "NAXIS1", &(stats.width), "length of data axis 1", &status)) { fits_report_error(stderr, status); return status; } // NAXIS2 if (fits_update_key(fptr, TUSHORT, "NAXIS2", &(stats.height), "length of data axis 2", &status)) { fits_report_error(stderr, status); return status; } // ISO Date if (fits_write_date(fptr, &status)) { fits_report_error(stderr, status); return status; } QString history = QString("Modified by KStars on %1").arg(QDateTime::currentDateTime().toString("yyyy-MM-ddThh:mm:ss")); // History if (fits_write_history(fptr, history.toLatin1(), &status)) { fits_report_error(stderr, status); return status; } int rot = 0, mirror = 0; if (rotCounter > 0) { rot = (90 * rotCounter) % 360; if (rot < 0) rot += 360; } if (flipHCounter % 2 != 0 || flipVCounter % 2 != 0) mirror = 1; if ((rot != 0) || (mirror != 0)) rotWCSFITS(rot, mirror); rotCounter = flipHCounter = flipVCounter = 0; if (m_isTemporary && autoRemoveTemporaryFITS) { QFile::remove(m_Filename); m_isTemporary = false; } m_Filename = newFilename; fits_flush_file(fptr, &status); qCInfo(KSTARS_FITS) << "Saved FITS file:" << m_Filename; return status; } void FITSData::clearImageBuffers() { delete[] m_ImageBuffer; m_ImageBuffer = nullptr; bayerBuffer = nullptr; } void FITSData::calculateStats(bool refresh) { // Calculate min max calculateMinMax(refresh); // Get standard deviation and mean in one run switch (m_DataType) { case TBYTE: runningAverageStdDev(); break; case TSHORT: runningAverageStdDev(); break; case TUSHORT: runningAverageStdDev(); break; case TLONG: runningAverageStdDev(); break; case TULONG: runningAverageStdDev(); break; case TFLOAT: runningAverageStdDev(); break; case TLONGLONG: runningAverageStdDev(); break; case TDOUBLE: runningAverageStdDev(); break; default: return; } // FIXME That's not really SNR, must implement a proper solution for this value stats.SNR = stats.mean[0] / stats.stddev[0]; if (refresh && markStars) // Let's try to find star positions again after transformation starsSearched = false; } int FITSData::calculateMinMax(bool refresh) { int status, nfound = 0; status = 0; if ((fptr != nullptr) && !refresh) { if (fits_read_key_dbl(fptr, "DATAMIN", &(stats.min[0]), nullptr, &status) == 0) nfound++; if (fits_read_key_dbl(fptr, "DATAMAX", &(stats.max[0]), nullptr, &status) == 0) nfound++; // If we found both keywords, no need to calculate them, unless they are both zeros if (nfound == 2 && !(stats.min[0] == 0 && stats.max[0] == 0)) return 0; } stats.min[0] = 1.0E30; stats.max[0] = -1.0E30; stats.min[1] = 1.0E30; stats.max[1] = -1.0E30; stats.min[2] = 1.0E30; stats.max[2] = -1.0E30; switch (m_DataType) { case TBYTE: calculateMinMax(); break; case TSHORT: calculateMinMax(); break; case TUSHORT: calculateMinMax(); break; case TLONG: calculateMinMax(); break; case TULONG: calculateMinMax(); break; case TFLOAT: calculateMinMax(); break; case TLONGLONG: calculateMinMax(); break; case TDOUBLE: calculateMinMax(); break; default: break; } //qDebug() << "DATAMIN: " << stats.min << " - DATAMAX: " << stats.max; return 0; } template QPair FITSData::getParitionMinMax(uint32_t start, uint32_t stride) { auto * buffer = reinterpret_cast(m_ImageBuffer); T min = std::numeric_limits::max(); T max = std::numeric_limits::min(); uint32_t end = start + stride; for (uint32_t i = start; i < end; i++) { if (buffer[i] < min) min = buffer[i]; else if (buffer[i] > max) max = buffer[i]; } return qMakePair(min, max); } template void FITSData::calculateMinMax() { - //QTime timer; - //timer.start(); - //if (filename.contains("thread")) - //{ T min = std::numeric_limits::max(); T max = std::numeric_limits::min(); // Create N threads const uint8_t nThreads = 16; for (int n = 0; n < m_Channels; n++) { uint32_t cStart = n * stats.samples_per_channel; // Calculate how many elements we process per thread uint32_t tStride = stats.samples_per_channel / nThreads; // Calculate the final stride since we can have some left over due to division above uint32_t fStride = tStride + (stats.samples_per_channel - (tStride * nThreads)); // Start location for inspecting elements uint32_t tStart = cStart; // List of futures QList>> futures; for (int i = 0; i < nThreads; i++) { // Run threads futures.append(QtConcurrent::run(this, &FITSData::getParitionMinMax, tStart, (i == (nThreads - 1)) ? fStride : tStride)); tStart += tStride; } // Now wait for results for (int i = 0; i < nThreads; i++) { QPair result = futures[i].result(); if (result.first < min) min = result.first; if (result.second > max) max = result.second; } stats.min[n] = min; stats.max[n] = max; } -#if 0 -} -else -{ - T * buffer = reinterpret_cast(imageBuffer); - if (channels == 1) - { - for (unsigned int i = 0; i < stats.samples_per_channel; i++) - { - if (buffer[i] < stats.min[0]) - stats.min[0] = buffer[i]; - else if (buffer[i] > stats.max[0]) - stats.max[0] = buffer[i]; - } - } - else - { - int g_offset = stats.samples_per_channel; - int b_offset = stats.samples_per_channel * 2; - - for (unsigned int i = 0; i < stats.samples_per_channel; i++) - { - if (buffer[i] < stats.min[0]) - stats.min[0] = buffer[i]; - else if (buffer[i] > stats.max[0]) - stats.max[0] = buffer[i]; - - if (buffer[i + g_offset] < stats.min[1]) - stats.min[1] = buffer[i + g_offset]; - else if (buffer[i + g_offset] > stats.max[1]) - stats.max[1] = buffer[i + g_offset]; - - if (buffer[i + b_offset] < stats.min[2]) - stats.min[2] = buffer[i + b_offset]; - else if (buffer[i + b_offset] > stats.max[2]) - stats.max[2] = buffer[i + b_offset]; - } - } - -} - -qCInfo(KSTARS_FITS) << filename << "MinMax calculation took" << timer.elapsed() << "ms"; -#endif } template QPair FITSData::getSquaredSumAndMean(uint32_t start, uint32_t stride) { uint32_t m_n = 2; double m_oldM = 0, m_newM = 0, m_oldS = 0, m_newS = 0; auto * buffer = reinterpret_cast(m_ImageBuffer); uint32_t end = start + stride; for (uint32_t i = start; i < end; i++) { m_newM = m_oldM + (buffer[i] - m_oldM) / m_n; m_newS = m_oldS + (buffer[i] - m_oldM) * (buffer[i] - m_newM); m_oldM = m_newM; m_oldS = m_newS; m_n++; } return qMakePair(m_newM, m_newS); } template void FITSData::runningAverageStdDev() { - //QTime timer; - //timer.start(); - //if (filename.contains("thread")) - //{ // Create N threads const uint8_t nThreads = 16; for (int n = 0; n < m_Channels; n++) { uint32_t cStart = n * stats.samples_per_channel; // Calculate how many elements we process per thread uint32_t tStride = stats.samples_per_channel / nThreads; // Calculate the final stride since we can have some left over due to division above uint32_t fStride = tStride + (stats.samples_per_channel - (tStride * nThreads)); // Start location for inspecting elements uint32_t tStart = cStart; // List of futures QList>> futures; for (int i = 0; i < nThreads; i++) { // Run threads futures.append(QtConcurrent::run(this, &FITSData::getSquaredSumAndMean, tStart, (i == (nThreads - 1)) ? fStride : tStride)); tStart += tStride; } double mean = 0, squared_sum = 0; // Now wait for results for (int i = 0; i < nThreads; i++) { QPair result = futures[i].result(); mean += result.first; squared_sum += result.second; } double variance = squared_sum / stats.samples_per_channel; stats.mean[n] = mean / nThreads; stats.stddev[n] = sqrt(variance); } -#if 0 -} -else -{ - T * buffer = reinterpret_cast(imageBuffer); - - if (channels == 1) - { - int m_n = 2; - double m_oldM = 0, m_newM = 0, m_oldS = 0, m_newS = 0; - - for (unsigned int i = 1; i < stats.samples_per_channel; i++) - { - m_newM = m_oldM + (buffer[i] - m_oldM) / m_n; - m_newS = m_oldS + (buffer[i] - m_oldM) * (buffer[i] - m_newM); - - m_oldM = m_newM; - m_oldS = m_newS; - m_n++; - } - - double variance = (m_n == 2 ? 0 : m_newS / (m_n - 2)); - - stats.mean[0] = m_newM; - stats.stddev[0] = sqrt(variance); - } - else - { - int m_n[3] = {2, 2, 2}; - double m_oldM[3] = {0}, m_newM[3] = {0}, m_oldS[3] = {0}, m_newS[3] = {0}; - - T * rBuffer = buffer; - T * gBuffer = buffer + stats.samples_per_channel; - T * bBuffer = buffer + stats.samples_per_channel * 2; - - for (unsigned int i = 1; i < stats.samples_per_channel; i++) - { - m_newM[0] = m_oldM[0] + (rBuffer[i] - m_oldM[0]) / m_n[0]; - m_newS[0] = m_oldS[0] + (rBuffer[i] - m_oldM[0]) * (rBuffer[i] - m_newM[0]); - - m_oldM[0] = m_newM[0]; - m_oldS[0] = m_newS[0]; - m_n[0]++; - - m_newM[1] = m_oldM[1] + (gBuffer[i] - m_oldM[1]) / m_n[1]; - m_newS[1] = m_oldS[1] + (gBuffer[i] - m_oldM[1]) * (gBuffer[i] - m_newM[1]); - - m_oldM[1] = m_newM[1]; - m_oldS[1] = m_newS[1]; - m_n[1]++; - - m_newM[2] = m_oldM[2] + (bBuffer[i] - m_oldM[2]) / m_n[2]; - m_newS[2] = m_oldS[2] + (bBuffer[i] - m_oldM[2]) * (bBuffer[i] - m_newM[2]); - - m_oldM[2] = m_newM[2]; - m_oldS[2] = m_newS[2]; - m_n[2]++; - - } - - double variance = (m_n[0] == 2 ? 0 : m_newS[0] / (m_n[0] - 2)); - stats.mean[0] = m_newM[0]; - stats.stddev[0] = sqrt(variance); - - variance = (m_n[1] == 2 ? 0 : m_newS[1] / (m_n[1] - 2)); - stats.mean[1] = m_newM[1]; - stats.stddev[1] = sqrt(variance); - - variance = (m_n[2] == 2 ? 0 : m_newS[2] / (m_n[2] - 2)); - stats.mean[2] = m_newM[2]; - stats.stddev[2] = sqrt(variance); - } -} - -qCInfo(KSTARS_FITS) << filename << "runningMeanStdDev calculation took" << timer.elapsed() << "ms"; -#endif } void FITSData::setMinMax(double newMin, double newMax, uint8_t channel) { stats.min[channel] = newMin; stats.max[channel] = newMax; } bool FITSData::parseHeader() { char * header = nullptr; int status = 0, nkeys = 0; if (fits_hdr2str(fptr, 0, nullptr, 0, &header, &nkeys, &status)) { fits_report_error(stderr, status); free(header); return false; } QString recordList = QString(header); for (int i = 0; i < nkeys; i++) { Record * oneRecord = new Record; // Quotes cause issues for simplified below so we're removing them. QString record = recordList.mid(i * 80, 80).remove("'"); QStringList properties = record.split(QRegExp("[=/]")); // If it is only a comment if (properties.size() == 1) { oneRecord->key = properties[0].mid(0, 7); oneRecord->comment = properties[0].mid(8).simplified(); } else { oneRecord->key = properties[0].simplified(); oneRecord->value = properties[1].simplified(); if (properties.size() > 2) oneRecord->comment = properties[2].simplified(); // Try to guess the value. // Test for integer & double. If neither, then leave it as "string". bool ok = false; // Is it Integer? oneRecord->value.toInt(&ok); if (ok) oneRecord->value.convert(QMetaType::Int); else { // Is it double? oneRecord->value.toDouble(&ok); if (ok) oneRecord->value.convert(QMetaType::Double); } } records.append(oneRecord); } free(header); return true; } bool FITSData::getRecordValue(const QString &key, QVariant &value) const { for (Record * oneRecord : records) { if (oneRecord->key == key) { value = oneRecord->value; return true; } } return false; } bool FITSData::checkCollision(Edge * s1, Edge * s2) { int dis; //distance int diff_x = s1->x - s2->x; int diff_y = s1->y - s2->y; dis = std::abs(sqrt(diff_x * diff_x + diff_y * diff_y)); dis -= s1->width / 2; dis -= s2->width / 2; if (dis <= 0) //collision return true; //no collision return false; } int FITSData::findCannyStar(FITSData * data, const QRect &boundary) { switch (data->property("dataType").toInt()) { case TBYTE: return FITSData::findCannyStar(data, boundary); case TSHORT: return FITSData::findCannyStar(data, boundary); case TUSHORT: return FITSData::findCannyStar(data, boundary); case TLONG: return FITSData::findCannyStar(data, boundary); case TULONG: return FITSData::findCannyStar(data, boundary); case TFLOAT: return FITSData::findCannyStar(data, boundary); case TLONGLONG: return FITSData::findCannyStar(data, boundary); case TDOUBLE: return FITSData::findCannyStar(data, boundary); default: break; } return 0; } int FITSData::findStars(StarAlgorithm algorithm, const QRect &trackingBox) { int count = 0; starAlgorithm = algorithm; qDeleteAll(starCenters); starCenters.clear(); switch (algorithm) { case ALGORITHM_SEP: count = findSEPStars(trackingBox); break; case ALGORITHM_GRADIENT: count = findCannyStar(this, trackingBox); break; case ALGORITHM_CENTROID: count = findCentroid(trackingBox); break; case ALGORITHM_THRESHOLD: count = findOneStar(trackingBox); break; } starsSearched = true; return count; } int FITSData::filterStars(const float innerRadius, const float outerRadius) { long const sqDiagonal = this->width() * this->width() / 4 + this->height() * this->height() / 4; long const sqInnerRadius = std::lround(sqDiagonal * innerRadius * innerRadius); long const sqOuterRadius = std::lround(sqDiagonal * outerRadius * outerRadius); starCenters.erase(std::remove_if(starCenters.begin(), starCenters.end(), [&](Edge * edge) { long const x = edge->x - this->width() / 2; long const y = edge->y - this->height() / 2; long const sqRadius = x * x + y * y; return sqRadius < sqInnerRadius || sqOuterRadius < sqRadius; }), starCenters.end()); return starCenters.count(); } template int FITSData::findCannyStar(FITSData * data, const QRect &boundary) { int subX = qMax(0, boundary.isNull() ? 0 : boundary.x()); int subY = qMax(0, boundary.isNull() ? 0 : boundary.y()); int subW = (boundary.isNull() ? data->width() : boundary.width()); int subH = (boundary.isNull() ? data->height() : boundary.height()); int BBP = data->getBytesPerPixel(); uint16_t dataWidth = data->width(); // #1 Find offsets uint32_t size = subW * subH; uint32_t offset = subX + subY * dataWidth; // #2 Create new buffer auto * buffer = new uint8_t[size * BBP]; // If there is no offset, copy whole buffer in one go if (offset == 0) memcpy(buffer, data->getImageBuffer(), size * BBP); else { uint8_t * dataPtr = buffer; uint8_t * origDataPtr = data->getImageBuffer(); uint32_t lineOffset = 0; // Copy data line by line for (int height = subY; height < (subY + subH); height++) { lineOffset = (subX + height * dataWidth) * BBP; memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP); dataPtr += (subW * BBP); } } // #3 Create new FITSData to hold it auto * boundedImage = new FITSData(); boundedImage->stats.width = subW; boundedImage->stats.height = subH; boundedImage->stats.bitpix = data->stats.bitpix; boundedImage->stats.bytesPerPixel = data->stats.bytesPerPixel; boundedImage->stats.samples_per_channel = size; boundedImage->stats.ndim = 2; boundedImage->setProperty("dataType", data->property("dataType")); // #4 Set image buffer and calculate stats. boundedImage->setImageBuffer(buffer); boundedImage->calculateStats(true); // #5 Apply Median + High Contrast filter to remove noise and move data to non-linear domain boundedImage->applyFilter(FITS_MEDIAN); boundedImage->applyFilter(FITS_HIGH_CONTRAST); // #6 Perform Sobel to find gradients and their directions QVector gradients; QVector directions; // TODO Must trace neighbours and assign IDs to each shape so that they can be centered massed // and discarded whenever necessary. It won't work on noisy images unless this is done. boundedImage->sobel(gradients, directions); QVector ids(gradients.size()); int maxID = boundedImage->partition(subW, subH, gradients, ids); - //QVector thresholded = boundedImage->threshold(boundedImage->stats.mean[0], boundedImage->stats.max[0], gradients); - // Not needed anymore delete boundedImage; if (maxID == 0) return 0; typedef struct { float massX = 0; float massY = 0; float totalMass = 0; } massInfo; QMap masses; // #7 Calculate center of mass for all detected regions for (int y = 0; y < subH; y++) { for (int x = 0; x < subW; x++) { int index = x + y * subW; int regionID = ids[index]; if (regionID > 0) { float pixel = gradients[index]; masses[regionID].totalMass += pixel; masses[regionID].massX += x * pixel; masses[regionID].massY += y * pixel; } } } // Compare multiple masses, and only select the highest total mass one as the desired star int maxRegionID = 1; int maxTotalMass = masses[1].totalMass; double totalMassRatio = 1e6; for (auto key : masses.keys()) { massInfo oneMass = masses.value(key); if (oneMass.totalMass > maxTotalMass) { totalMassRatio = oneMass.totalMass / maxTotalMass; maxTotalMass = oneMass.totalMass; maxRegionID = key; } } // If image has many regions and there is no significant relative center of mass then it's just noise and no stars // are probably there above a useful threshold. if (maxID > 10 && totalMassRatio < 1.5) return 0; auto * center = new Edge; center->width = -1; center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5; center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5; center->HFR = 1; // Maximum Radius int maxR = qMin(subW - 1, subH - 1) / 2; for (int r = maxR; r > 1; r--) { int pass = 0; for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0) { int testX = center->x + std::cos(theta) * r; int testY = center->y + std::sin(theta) * r; // if out of bound, break; if (testX < 0 || testX >= subW || testY < 0 || testY >= subH) break; if (gradients[testX + testY * subW] > 0) //if (thresholded[testX + testY * subW] > 0) { if (++pass >= 24) { center->width = r * 2; // Break of outer loop r = 0; break; } } } } qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << center->x << " Y: " << center->y << " Width: " << center->width; // If no stars were detected if (center->width == -1) { delete center; return 0; } // 30% fuzzy //center->width += center->width*0.3 * (running_threshold / threshold); double FSum = 0, HF = 0, TF = 0; const double resolution = 1.0 / 20.0; int cen_y = qRound(center->y); double rightEdge = center->x + center->width / 2.0; double leftEdge = center->x - center->width / 2.0; QVector subPixels; subPixels.reserve(center->width / resolution); const T * origBuffer = reinterpret_cast(data->getImageBuffer()) + offset; - /*if (Options::fITSLogging()) - { - QDebug deb = qDebug(); - - for (int i=0; i < subW; i++) - deb << origBuffer[i + cen_y * dataWidth] << ","; - }*/ - for (double x = leftEdge; x <= rightEdge; x += resolution) { double slice = resolution * (origBuffer[static_cast(floor(x)) + cen_y * dataWidth]); FSum += slice; subPixels.append(slice); } // Half flux HF = FSum / 2.0; int subPixelCenter = (center->width / resolution) / 2; // Start from center TF = subPixels[subPixelCenter]; double lastTF = TF; // Integrate flux along radius axis until we reach half flux //for (double k=resolution; k < (center->width/(2*resolution)); k += resolution) for (int k = 1; k < subPixelCenter; k++) { TF += subPixels[subPixelCenter + k]; TF += subPixels[subPixelCenter - k]; if (TF >= HF) { // We overpassed HF, let's calculate from last TF how much until we reach HF // #1 Accurate calculation, but very sensitive to small variations of flux center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution; // #2 Less accurate calculation, but stable against small variations of flux //center->HFR = (k - 1) * resolution; break; } lastTF = TF; } // Correct center for subX and subY center->x += subX; center->y += subY; data->appendStar(center); qCDebug(KSTARS_FITS) << "Flux: " << FSum << " Half-Flux: " << HF << " HFR: " << center->HFR; return 1; } int FITSData::findOneStar(const QRect &boundary) { switch (m_DataType) { case TBYTE: return findOneStar(boundary); case TSHORT: return findOneStar(boundary); case TUSHORT: return findOneStar(boundary); case TLONG: return findOneStar(boundary); case TULONG: return findOneStar(boundary); case TFLOAT: return findOneStar(boundary); case TLONGLONG: return findOneStar(boundary); case TDOUBLE: return findOneStar(boundary); default: break; } return 0; } template int FITSData::findOneStar(const QRect &boundary) { if (boundary.isEmpty()) return -1; int subX = boundary.x(); int subY = boundary.y(); int subW = subX + boundary.width(); int subH = subY + boundary.height(); float massX = 0, massY = 0, totalMass = 0; auto * buffer = reinterpret_cast(m_ImageBuffer); // TODO replace magic number with something more useful to understand double threshold = stats.mean[0] * Options::focusThreshold() / 100.0; for (int y = subY; y < subH; y++) { for (int x = subX; x < subW; x++) { T pixel = buffer[x + y * stats.width]; if (pixel > threshold) { totalMass += pixel; massX += x * pixel; massY += y * pixel; } } } qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << massX / totalMass << " Y: " << massY / totalMass; auto * center = new Edge; center->width = -1; center->x = massX / totalMass + 0.5; center->y = massY / totalMass + 0.5; center->HFR = 1; // Maximum Radius int maxR = qMin(subW - 1, subH - 1) / 2; // Critical threshold double critical_threshold = threshold * 0.7; double running_threshold = threshold; while (running_threshold >= critical_threshold) { for (int r = maxR; r > 1; r--) { int pass = 0; for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 10.0) { int testX = center->x + std::cos(theta) * r; int testY = center->y + std::sin(theta) * r; // if out of bound, break; if (testX < subX || testX > subW || testY < subY || testY > subH) break; if (buffer[testX + testY * stats.width] > running_threshold) pass++; } //qDebug() << "Testing for radius " << r << " passes # " << pass << " @ threshold " << running_threshold; //if (pass >= 6) if (pass >= 5) { center->width = r * 2; break; } } if (center->width > 0) break; // Increase threshold fuzziness by 10% running_threshold -= running_threshold * 0.1; } // If no stars were detected if (center->width == -1) { delete center; return 0; } // 30% fuzzy //center->width += center->width*0.3 * (running_threshold / threshold); starCenters.append(center); double FSum = 0, HF = 0, TF = 0, min = stats.min[0]; const double resolution = 1.0 / 20.0; int cen_y = qRound(center->y); double rightEdge = center->x + center->width / 2.0; double leftEdge = center->x - center->width / 2.0; QVector subPixels; subPixels.reserve(center->width / resolution); for (double x = leftEdge; x <= rightEdge; x += resolution) { //subPixels[x] = resolution * (image_buffer[static_cast(floor(x)) + cen_y * stats.width] - min); double slice = resolution * (buffer[static_cast(floor(x)) + cen_y * stats.width] - min); FSum += slice; subPixels.append(slice); } // Half flux HF = FSum / 2.0; //double subPixelCenter = center->x - fmod(center->x,resolution); int subPixelCenter = (center->width / resolution) / 2; // Start from center TF = subPixels[subPixelCenter]; double lastTF = TF; // Integrate flux along radius axis until we reach half flux //for (double k=resolution; k < (center->width/(2*resolution)); k += resolution) for (int k = 1; k < subPixelCenter; k++) { TF += subPixels[subPixelCenter + k]; TF += subPixels[subPixelCenter - k]; if (TF >= HF) { // We have two ways to calculate HFR. The first is the correct method but it can get quite variable within 10% due to random fluctuations of the measured star. // The second method is not truly HFR but is much more resistant to noise. // #1 Approximate HFR, accurate and reliable but quite variable to small changes in star flux center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution; // #2 Not exactly HFR, but much more stable //center->HFR = (k*resolution) * (HF/TF); break; } lastTF = TF; } return 1; } /*** Find center of stars and calculate Half Flux Radius */ int FITSData::findCentroid(const QRect &boundary, int initStdDev, int minEdgeWidth) { switch (m_DataType) { case TBYTE: return findCentroid(boundary, initStdDev, minEdgeWidth); case TSHORT: return findCentroid(boundary, initStdDev, minEdgeWidth); case TUSHORT: return findCentroid(boundary, initStdDev, minEdgeWidth); case TLONG: return findCentroid(boundary, initStdDev, minEdgeWidth); case TULONG: return findCentroid(boundary, initStdDev, minEdgeWidth); case TFLOAT: return findCentroid(boundary, initStdDev, minEdgeWidth); case TLONGLONG: return findCentroid(boundary, initStdDev, minEdgeWidth); case TDOUBLE: return findCentroid(boundary, initStdDev, minEdgeWidth); default: return -1; } } template int FITSData::findCentroid(const QRect &boundary, int initStdDev, int minEdgeWidth) { double threshold = 0, sum = 0, avg = 0, min = 0; int starDiameter = 0; int pixVal = 0; int minimumEdgeCount = MINIMUM_EDGE_LIMIT; auto * buffer = reinterpret_cast(m_ImageBuffer); double JMIndex = 100; #ifndef KSTARS_LITE if (histogram) { if (!histogram->isConstructed()) histogram->constructHistogram(); JMIndex = histogram->getJMIndex(); } #endif float dispersion_ratio = 1.5; QList edges; if (JMIndex < DIFFUSE_THRESHOLD) { minEdgeWidth = JMIndex * 35 + 1; minimumEdgeCount = minEdgeWidth - 1; } else { minEdgeWidth = 6; minimumEdgeCount = 4; } while (initStdDev >= 1) { minEdgeWidth--; minimumEdgeCount--; minEdgeWidth = qMax(3, minEdgeWidth); minimumEdgeCount = qMax(3, minimumEdgeCount); if (JMIndex < DIFFUSE_THRESHOLD) { // Taking the average out seems to have better result for noisy images threshold = stats.max[0] - stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1); min = stats.min[0]; if (threshold - min < 0) { threshold = stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1); min = 0; } dispersion_ratio = 1.4 - (MINIMUM_STDVAR - initStdDev) * 0.08; } else { threshold = stats.mean[0] + stats.stddev[0] * initStdDev * (0.3 - (MINIMUM_STDVAR - initStdDev) * 0.05); min = stats.min[0]; // Ratio between centeroid center and edge dispersion_ratio = 1.8 - (MINIMUM_STDVAR - initStdDev) * 0.2; } qCDebug(KSTARS_FITS) << "SNR: " << stats.SNR; qCDebug(KSTARS_FITS) << "The threshold level is " << threshold << "(actual " << threshold - min << ") minimum edge width" << minEdgeWidth << " minimum edge limit " << minimumEdgeCount; threshold -= min; int subX, subY, subW, subH; if (boundary.isNull()) { if (m_Mode == FITS_GUIDE || m_Mode == FITS_FOCUS) { // Only consider the central 70% subX = round(stats.width * 0.15); subY = round(stats.height * 0.15); subW = stats.width - subX; subH = stats.height - subY; } else { // Consider the complete area 100% subX = 0; subY = 0; subW = stats.width; subH = stats.height; } } else { subX = boundary.x(); subY = boundary.y(); subW = subX + boundary.width(); subH = subY + boundary.height(); } // Detect "edges" that are above threshold for (int i = subY; i < subH; i++) { starDiameter = 0; for (int j = subX; j < subW; j++) { pixVal = buffer[j + (i * stats.width)] - min; // If pixel value > threshold, let's get its weighted average if (pixVal >= threshold) { avg += j * pixVal; sum += pixVal; starDiameter++; } // Value < threshold but avg exists else if (sum > 0) { // We found a potential centroid edge if (starDiameter >= minEdgeWidth) { float center = avg / sum + 0.5; if (center > 0) { int i_center = std::floor(center); // Check if center is 10% or more brighter than edge, if not skip if (((buffer[i_center + (i * stats.width)] - min) / (buffer[i_center + (i * stats.width) - starDiameter / 2] - min) >= dispersion_ratio) && ((buffer[i_center + (i * stats.width)] - min) / (buffer[i_center + (i * stats.width) + starDiameter / 2] - min) >= dispersion_ratio)) { qCDebug(KSTARS_FITS) << "Edge center is " << buffer[i_center + (i * stats.width)] - min << " Edge is " << buffer[i_center + (i * stats.width) - starDiameter / 2] - min << " and ratio is " << ((buffer[i_center + (i * stats.width)] - min) / (buffer[i_center + (i * stats.width) - starDiameter / 2] - min)) << " located at X: " << center << " Y: " << i + 0.5; auto * newEdge = new Edge(); newEdge->x = center; newEdge->y = i + 0.5; newEdge->scanned = 0; newEdge->val = buffer[i_center + (i * stats.width)] - min; newEdge->width = starDiameter; newEdge->HFR = 0; newEdge->sum = sum; edges.append(newEdge); } } } // Reset avg = sum = starDiameter = 0; } } } qCDebug(KSTARS_FITS) << "Total number of edges found is: " << edges.count(); // In case of hot pixels if (edges.count() == 1 && initStdDev > 1) { initStdDev--; continue; } if (edges.count() >= MAX_EDGE_LIMIT) { qCWarning(KSTARS_FITS) << "Too many edges, aborting... " << edges.count(); qDeleteAll(edges); return -1; } if (edges.count() >= minimumEdgeCount) break; qDeleteAll(edges); edges.clear(); initStdDev--; } int cen_count = 0; int cen_x = 0; int cen_y = 0; int cen_v = 0; int cen_w = 0; int width_sum = 0; // Let's sort edges, starting with widest std::sort(edges.begin(), edges.end(), greaterThan); // Now, let's scan the edges and find the maximum centroid vertically for (int i = 0; i < edges.count(); i++) { qCDebug(KSTARS_FITS) << "# " << i << " Edge at (" << edges[i]->x << "," << edges[i]->y << ") With a value of " << edges[i]->val << " and width of " << edges[i]->width << " pixels. with sum " << edges[i]->sum; // If edge scanned already, skip if (edges[i]->scanned == 1) { qCDebug(KSTARS_FITS) << "Skipping check for center " << i << " because it was already counted"; continue; } qCDebug(KSTARS_FITS) << "Investigating edge # " << i << " now ..."; // Get X, Y, and Val of edge cen_x = edges[i]->x; cen_y = edges[i]->y; cen_v = edges[i]->sum; cen_w = edges[i]->width; float avg_x = 0; float avg_y = 0; sum = 0; cen_count = 0; // Now let's compare to other edges until we hit a maxima for (int j = 0; j < edges.count(); j++) { if (edges[j]->scanned) continue; if (checkCollision(edges[j], edges[i])) { if (edges[j]->sum >= cen_v) { cen_v = edges[j]->sum; cen_w = edges[j]->width; } edges[j]->scanned = 1; cen_count++; avg_x += edges[j]->x * edges[j]->val; avg_y += edges[j]->y * edges[j]->val; sum += edges[j]->val; continue; } } int cen_limit = (MINIMUM_ROWS_PER_CENTER - (MINIMUM_STDVAR - initStdDev)); if (edges.count() < LOW_EDGE_CUTOFF_1) { if (edges.count() < LOW_EDGE_CUTOFF_2) cen_limit = 1; else cen_limit = 2; } qCDebug(KSTARS_FITS) << "center_count: " << cen_count << " and initstdDev= " << initStdDev << " and limit is " << cen_limit; if (cen_limit < 1) continue; // If centroid count is within acceptable range //if (cen_limit >= 2 && cen_count >= cen_limit) if (cen_count >= cen_limit) { // We detected a centroid, let's init it auto * rCenter = new Edge(); rCenter->x = avg_x / sum; rCenter->y = avg_y / sum; width_sum += rCenter->width; rCenter->width = cen_w; qCDebug(KSTARS_FITS) << "Found a real center with number with (" << rCenter->x << "," << rCenter->y << ")"; // Calculate Total Flux From Center, Half Flux, Full Summation double TF = 0; double HF = 0; double FSum = 0; cen_x = (int)std::floor(rCenter->x); cen_y = (int)std::floor(rCenter->y); if (cen_x < 0 || cen_x > stats.width || cen_y < 0 || cen_y > stats.height) { delete rCenter; continue; } // Complete sum along the radius //for (int k=0; k < rCenter->width; k++) for (int k = rCenter->width / 2; k >= -(rCenter->width / 2); k--) { FSum += buffer[cen_x - k + (cen_y * stats.width)] - min; //qDebug() << image_buffer[cen_x-k+(cen_y*stats.width)] - min; } // Half flux HF = FSum / 2.0; // Total flux starting from center TF = buffer[cen_y * stats.width + cen_x] - min; int pixelCounter = 1; // Integrate flux along radius axis until we reach half flux for (int k = 1; k < rCenter->width / 2; k++) { if (TF >= HF) { qCDebug(KSTARS_FITS) << "Stopping at TF " << TF << " after #" << k << " pixels."; break; } TF += buffer[cen_y * stats.width + cen_x + k] - min; TF += buffer[cen_y * stats.width + cen_x - k] - min; pixelCounter++; } // Calculate weighted Half Flux Radius rCenter->HFR = pixelCounter * (HF / TF); // Store full flux rCenter->val = FSum; qCDebug(KSTARS_FITS) << "HFR for this center is " << rCenter->HFR << " pixels and the total flux is " << FSum; starCenters.append(rCenter); } } if (starCenters.count() > 1 && m_Mode != FITS_FOCUS) { float width_avg = (float)width_sum / starCenters.count(); float lsum = 0, sdev = 0; for (auto ¢er : starCenters) lsum += (center->width - width_avg) * (center->width - width_avg); sdev = (std::sqrt(lsum / (starCenters.count() - 1))) * 4; // Reject stars > 4 * stddev foreach (Edge * center, starCenters) if (center->width > sdev) starCenters.removeOne(center); //foreach(Edge *center, starCenters) //qDebug() << center->x << "," << center->y << "," << center->width << "," << center->val << endl; } // Release memory qDeleteAll(edges); return starCenters.count(); } double FITSData::getHFR(HFRType type) { // This method is less susceptible to noise // Get HFR for the brightest star only, instead of averaging all stars // It is more consistent. // TODO: Try to test this under using a real CCD. if (starCenters.empty()) return -1; if (type == HFR_MAX) { maxHFRStar = nullptr; int maxVal = 0; int maxIndex = 0; for (int i = 0; i < starCenters.count(); i++) { if (starCenters[i]->HFR > maxVal) { maxIndex = i; maxVal = starCenters[i]->HFR; } } maxHFRStar = starCenters[maxIndex]; return static_cast(starCenters[maxIndex]->HFR); } QVector HFRs; for (auto center : starCenters) HFRs << center->HFR; std::sort(HFRs.begin(), HFRs.end()); double sum = std::accumulate(HFRs.begin(), HFRs.end(), 0.0); double m = sum / HFRs.size(); if (HFRs.size() > 3) { double accum = 0.0; std::for_each (HFRs.begin(), HFRs.end(), [&](const double d) { accum += (d - m) * (d - m); }); double stddev = sqrt(accum / (HFRs.size() - 1)); // Remove stars over 2 standard deviations away. auto end1 = std::remove_if(HFRs.begin(), HFRs.end(), [m, stddev](const double hfr) { return hfr > (m + stddev * 2); }); auto end2 = std::remove_if(HFRs.begin(), end1, [m, stddev](const double hfr) { return hfr < (m - stddev * 2); }); // New mean sum = std::accumulate(HFRs.begin(), end2, 0.0); const int num_remaining = std::distance(HFRs.begin(), end2); if (num_remaining > 0) m = sum / num_remaining; } return m; } double FITSData::getHFR(int x, int y) { if (starCenters.empty()) return -1; for (int i = 0; i < starCenters.count(); i++) { if (std::fabs(starCenters[i]->x - x) <= starCenters[i]->width / 2 && std::fabs(starCenters[i]->y - y) <= starCenters[i]->width / 2) { return starCenters[i]->HFR; } } return -1; } void FITSData::applyFilter(FITSScale type, uint8_t * image, QVector * min, QVector * max) { if (type == FITS_NONE) return; QVector dataMin(3); QVector dataMax(3); if (min) dataMin = *min; if (max) dataMax = *max; switch (type) { case FITS_AUTO_STRETCH: { for (int i = 0; i < 3; i++) { dataMin[i] = stats.mean[i] - stats.stddev[i]; dataMax[i] = stats.mean[i] + stats.stddev[i] * 3; } } break; case FITS_HIGH_CONTRAST: { for (int i = 0; i < 3; i++) { dataMin[i] = stats.mean[i] + stats.stddev[i]; dataMax[i] = stats.mean[i] + stats.stddev[i] * 3; } } break; case FITS_HIGH_PASS: { for (int i = 0; i < 3; i++) { dataMin[i] = stats.mean[i]; } } break; default: break; } switch (m_DataType) { case TBYTE: { for (int i = 0; i < 3; i++) { dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i]; dataMax[i] = dataMax[i] > UINT8_MAX ? UINT8_MAX : dataMax[i]; } applyFilter(type, image, &dataMin, &dataMax); } break; case TSHORT: { for (int i = 0; i < 3; i++) { dataMin[i] = dataMin[i] < INT16_MIN ? INT16_MIN : dataMin[i]; dataMax[i] = dataMax[i] > INT16_MAX ? INT16_MAX : dataMax[i]; } applyFilter(type, image, &dataMin, &dataMax); } break; case TUSHORT: { for (int i = 0; i < 3; i++) { dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i]; dataMax[i] = dataMax[i] > UINT16_MAX ? UINT16_MAX : dataMax[i]; } applyFilter(type, image, &dataMin, &dataMax); } break; case TLONG: { for (int i = 0; i < 3; i++) { dataMin[i] = dataMin[i] < INT_MIN ? INT_MIN : dataMin[i]; dataMax[i] = dataMax[i] > INT_MAX ? INT_MAX : dataMax[i]; } applyFilter(type, image, &dataMin, &dataMax); } break; case TULONG: { for (int i = 0; i < 3; i++) { dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i]; dataMax[i] = dataMax[i] > UINT_MAX ? UINT_MAX : dataMax[i]; } applyFilter(type, image, &dataMin, &dataMax); } break; case TFLOAT: { for (int i = 0; i < 3; i++) { dataMin[i] = dataMin[i] < FLT_MIN ? FLT_MIN : dataMin[i]; dataMax[i] = dataMax[i] > FLT_MAX ? FLT_MAX : dataMax[i]; } applyFilter(type, image, &dataMin, &dataMax); } break; case TLONGLONG: { for (int i = 0; i < 3; i++) { dataMin[i] = dataMin[i] < LLONG_MIN ? LLONG_MIN : dataMin[i]; dataMax[i] = dataMax[i] > LLONG_MAX ? LLONG_MAX : dataMax[i]; } applyFilter(type, image, &dataMin, &dataMax); } break; case TDOUBLE: { for (int i = 0; i < 3; i++) { dataMin[i] = dataMin[i] < DBL_MIN ? DBL_MIN : dataMin[i]; dataMax[i] = dataMax[i] > DBL_MAX ? DBL_MAX : dataMax[i]; } applyFilter(type, image, &dataMin, &dataMax); } break; default: return; } if (min != nullptr) *min = dataMin; if (max != nullptr) *max = dataMax; } template void FITSData::applyFilter(FITSScale type, uint8_t * targetImage, QVector * targetMin, QVector * targetMax) { bool calcStats = false; T * image = nullptr; if (targetImage) image = reinterpret_cast(targetImage); else { image = reinterpret_cast(m_ImageBuffer); calcStats = true; } T min[3], max[3]; for (int i = 0; i < 3; i++) { min[i] = (*targetMin)[i] < std::numeric_limits::min() ? std::numeric_limits::min() : (*targetMin)[i]; max[i] = (*targetMax)[i] > std::numeric_limits::max() ? std::numeric_limits::max() : (*targetMax)[i]; } // Create N threads const uint8_t nThreads = 16; uint32_t width = stats.width; uint32_t height = stats.height; //QTime timer; //timer.start(); switch (type) { case FITS_AUTO: case FITS_LINEAR: case FITS_AUTO_STRETCH: case FITS_HIGH_CONTRAST: case FITS_LOG: case FITS_SQRT: case FITS_HIGH_PASS: { // List of futures QList> futures; QVector coeff(3); if (type == FITS_LOG) { for (int i = 0; i < 3; i++) coeff[i] = max[i] / std::log(1 + max[i]); } else if (type == FITS_SQRT) { for (int i = 0; i < 3; i++) coeff[i] = max[i] / sqrt(max[i]); } for (int n = 0; n < m_Channels; n++) { if (type == FITS_HIGH_PASS) min[n] = stats.mean[n]; uint32_t cStart = n * stats.samples_per_channel; // Calculate how many elements we process per thread uint32_t tStride = stats.samples_per_channel / nThreads; // Calculate the final stride since we can have some left over due to division above uint32_t fStride = tStride + (stats.samples_per_channel - (tStride * nThreads)); T * runningBuffer = image + cStart; if (type == FITS_LOG) { for (int i = 0; i < nThreads; i++) { // Run threads futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max, coeff, n](T & a) { a = qBound(min[n], static_cast(round(coeff[n] * std::log(1 + qBound(min[n], a, max[n])))), max[n]); })); runningBuffer += tStride; } } else if (type == FITS_SQRT) { for (int i = 0; i < nThreads; i++) { // Run threads futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max, coeff, n](T & a) { a = qBound(min[n], static_cast(round(coeff[n] * a)), max[n]); })); } runningBuffer += tStride; } else { for (int i = 0; i < nThreads; i++) { // Run threads futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max, n](T & a) { a = qBound(min[n], a, max[n]); })); runningBuffer += tStride; } } } for (int i = 0; i < nThreads * m_Channels; i++) futures[i].waitForFinished(); if (calcStats) { for (int i = 0; i < 3; i++) { stats.min[i] = min[i]; stats.max[i] = max[i]; } //if (type != FITS_AUTO && type != FITS_LINEAR) runningAverageStdDev(); //QtConcurrent::run(this, &FITSData::runningAverageStdDev); } } break; case FITS_EQUALIZE: { #ifndef KSTARS_LITE if (histogram == nullptr) return; if (!histogram->isConstructed()) histogram->constructHistogram(); T bufferVal = 0; QVector cumulativeFreq = histogram->getCumulativeFrequency(); double coeff = 255.0 / (height * width); uint32_t row = 0; uint32_t index = 0; for (int i = 0; i < m_Channels; i++) { uint32_t offset = i * stats.samples_per_channel; for (uint32_t j = 0; j < height; j++) { row = offset + j * width; for (uint32_t k = 0; k < width; k++) { index = k + row; bufferVal = (image[index] - min[i]) / histogram->getBinWidth(i); if (bufferVal >= cumulativeFreq.size()) bufferVal = cumulativeFreq.size() - 1; image[index] = qBound(min[i], static_cast(round(coeff * cumulativeFreq[bufferVal])), max[i]); } } } #endif } if (calcStats) calculateStats(true); break; // Based on http://www.librow.com/articles/article-1 case FITS_MEDIAN: { uint8_t BBP = stats.bytesPerPixel; auto * extension = new T[(width + 2) * (height + 2)]; // Check memory allocation if (!extension) return; // Create image extension for (uint32_t ch = 0; ch < m_Channels; ch++) { uint32_t offset = ch * stats.samples_per_channel; uint32_t N = width, M = height; for (uint32_t i = 0; i < M; ++i) { memcpy(extension + (N + 2) * (i + 1) + 1, image + (N * i) + offset, N * BBP); extension[(N + 2) * (i + 1)] = image[N * i + offset]; extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1 + offset]; } // Fill first line of image extension memcpy(extension, extension + N + 2, (N + 2) * BBP); // Fill last line of image extension memcpy(extension + (N + 2) * (M + 1), extension + (N + 2) * M, (N + 2) * BBP); // Call median filter implementation N = width + 2; M = height + 2; // Move window through all elements of the image for (uint32_t m = 1; m < M - 1; ++m) for (uint32_t n = 1; n < N - 1; ++n) { // Pick up window elements int k = 0; float window[9]; memset(&window[0], 0, 9 * sizeof(float)); for (uint32_t j = m - 1; j < m + 2; ++j) for (uint32_t i = n - 1; i < n + 2; ++i) window[k++] = extension[j * N + i]; // Order elements (only half of them) for (uint32_t j = 0; j < 5; ++j) { // Find position of minimum element int mine = j; for (uint32_t l = j + 1; l < 9; ++l) if (window[l] < window[mine]) mine = l; // Put found minimum element in its place const float temp = window[j]; window[j] = window[mine]; window[mine] = temp; } // Get result - the middle element image[(m - 1) * (N - 2) + n - 1 + offset] = window[4]; } } // Free memory delete[] extension; if (calcStats) runningAverageStdDev(); } break; case FITS_ROTATE_CW: rotFITS(90, 0); rotCounter++; break; case FITS_ROTATE_CCW: rotFITS(270, 0); rotCounter--; break; case FITS_FLIP_H: rotFITS(0, 1); flipHCounter++; break; case FITS_FLIP_V: rotFITS(0, 2); flipVCounter++; break; default: break; } -#if 0 -} -else -{ - uint32_t index = 0, row = 0, offset = 0; - - switch (type) - { - case FITS_AUTO: - case FITS_LINEAR: - { - for (uint8_t i = 0; i < channels; i++) - { - offset = i * stats.samples_per_channel; - for (uint32_t j = 0; j < height; j++) - { - row = offset + j * width; - for (uint32_t k = 0; k < width; k++) - { - index = k + row; - image[index] = qBound(min, image[index], max); - } - } - } - - if (calcStats) - { - stats.min[0] = min; - stats.max[0] = max; - } - } - break; - - case FITS_LOG: - { - double coeff = max / log(1 + max); - - for (int i = 0; i < channels; i++) - { - offset = i * stats.samples_per_channel; - for (uint32_t j = 0; j < height; j++) - { - row = offset + j * width; - for (uint32_t k = 0; k < width; k++) - { - index = k + row; - image[index] = qBound(min, static_cast(round(coeff * log(1 + qBound(min, image[index], max)))), max); - } - } - } - - if (calcStats) - { - stats.min[0] = min; - stats.max[0] = max; - runningAverageStdDev(); - } - } - break; - - case FITS_SQRT: - { - double coeff = max / sqrt(max); - - for (int i = 0; i < channels; i++) - { - offset = i * stats.samples_per_channel; - for (uint32_t j = 0; j < height; j++) - { - row = offset + j * width; - for (uint32_t k = 0; k < width; k++) - { - index = k + row; - image[index] = qBound(min, static_cast(round(coeff * image[index])), max); - } - } - } - - if (calcStats) - { - stats.min[0] = min; - stats.max[0] = max; - runningAverageStdDev(); - } - } - break; - - // Only difference is how min and max are set - case FITS_AUTO_STRETCH: - case FITS_HIGH_CONTRAST: - { - for (uint32_t i = 0; i < channels; i++) - { - offset = i * stats.samples_per_channel; - for (uint32_t j = 0; j < height; j++) - { - row = offset + j * width; - for (uint32_t k = 0; k < width; k++) - image[k + row] = qBound(min, image[k + row], max); - } - } - if (calcStats) - { - stats.min[0] = min; - stats.max[0] = max; - runningAverageStdDev(); - } - } - break; - - case FITS_EQUALIZE: - { -#ifndef KSTARS_LITE - if (histogram == nullptr) - return; - - T bufferVal = 0; - QVector cumulativeFreq = histogram->getCumulativeFrequency(); - - double coeff = 255.0 / (height * width); - - for (uint32_t i = 0; i < channels; i++) - { - offset = i * stats.samples_per_channel; - for (uint32_t j = 0; j < height; j++) - { - row = offset + j * width; - for (uint32_t k = 0; k < width; k++) - { - index = k + row; - bufferVal = (image[index] - min) / histogram->getBinWidth(); - - if (bufferVal >= cumulativeFreq.size()) - bufferVal = cumulativeFreq.size() - 1; - - image[index] = qBound(min, static_cast(round(coeff * cumulativeFreq[bufferVal])), max); - } - } - } -#endif - } - if (calcStats) - calculateStats(true); - break; - - case FITS_HIGH_PASS: - { - min = stats.mean[0]; - for (uint32_t i = 0; i < channels; i++) - { - offset = i * stats.samples_per_channel; - for (uint32_t j = 0; j < height; j++) - { - row = offset + j * width; - for (uint32_t k = 0; k < width; k++) - { - index = k + row; - image[index] = qBound(min, image[index], max); - } - } - } - - if (calcStats) - { - stats.min[0] = min; - stats.max[0] = max; - runningAverageStdDev(); - } - } - break; - - // Based on http://www.librow.com/articles/article-1 - case FITS_MEDIAN: - { - int BBP = stats.bytesPerPixel; - T * extension = new T[(width + 2) * (height + 2)]; - // Check memory allocation - if (!extension) - return; - // Create image extension - for (uint32_t ch = 0; ch < channels; ch++) - { - offset = ch * stats.samples_per_channel; - uint32_t N = width, M = height; - - for (uint32_t i = 0; i < M; ++i) - { - memcpy(extension + (N + 2) * (i + 1) + 1, image + (N * i) + offset, N * BBP); - extension[(N + 2) * (i + 1)] = image[N * i + offset]; - extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1 + offset]; - } - // Fill first line of image extension - memcpy(extension, extension + N + 2, (N + 2) * BBP); - // Fill last line of image extension - memcpy(extension + (N + 2) * (M + 1), extension + (N + 2) * M, (N + 2) * BBP); - // Call median filter implementation - - N = width + 2; - M = height + 2; - - // Move window through all elements of the image - for (uint32_t m = 1; m < M - 1; ++m) - for (uint32_t n = 1; n < N - 1; ++n) - { - // Pick up window elements - int k = 0; - float window[9]; - - memset(&window[0], 0, 9 * sizeof(float)); - for (uint32_t j = m - 1; j < m + 2; ++j) - for (uint32_t i = n - 1; i < n + 2; ++i) - window[k++] = extension[j * N + i]; - // Order elements (only half of them) - for (uint32_t j = 0; j < 5; ++j) - { - // Find position of minimum element - int mine = j; - for (uint32_t l = j + 1; l < 9; ++l) - if (window[l] < window[mine]) - mine = l; - // Put found minimum element in its place - const float temp = window[j]; - window[j] = window[mine]; - window[mine] = temp; - } - // Get result - the middle element - image[(m - 1) * (N - 2) + n - 1 + offset] = window[4]; - } - } - - // Free memory - delete[] extension; - - if (calcStats) - runningAverageStdDev(); - } - break; - - case FITS_ROTATE_CW: - rotFITS(90, 0); - rotCounter++; - break; - - case FITS_ROTATE_CCW: - rotFITS(270, 0); - rotCounter--; - break; - - case FITS_FLIP_H: - rotFITS(0, 1); - flipHCounter++; - break; - - case FITS_FLIP_V: - rotFITS(0, 2); - flipVCounter++; - break; - - case FITS_CUSTOM: - default: - return; - break; - } -} - -qCInfo(KSTARS_FITS) << filename << "Apply Filter calculation took" << timer.elapsed() << "ms"; -#endif } QList FITSData::getStarCentersInSubFrame(QRect subFrame) const { QList starCentersInSubFrame; for (int i = 0; i < starCenters.count(); i++) { int x = static_cast(starCenters[i]->x); int y = static_cast(starCenters[i]->y); if(subFrame.contains(x, y)) { starCentersInSubFrame.append(starCenters[i]); } } return starCentersInSubFrame; } void FITSData::getCenterSelection(int * x, int * y) { if (starCenters.count() == 0) return; auto * pEdge = new Edge(); pEdge->x = *x; pEdge->y = *y; pEdge->width = 1; foreach (Edge * center, starCenters) if (checkCollision(pEdge, center)) { *x = static_cast(center->x); *y = static_cast(center->y); break; } delete (pEdge); } bool FITSData::checkForWCS() { #ifndef KSTARS_LITE #ifdef HAVE_WCSLIB int status = 0; char * header; int nkeyrec, nreject, nwcs; if (fits_hdr2str(fptr, 1, nullptr, 0, &header, &nkeyrec, &status)) { char errmsg[512]; fits_get_errstatus(status, errmsg); lastError = errmsg; return false; } if ((status = wcspih(header, nkeyrec, WCSHDR_all, -3, &nreject, &nwcs, &wcs)) != 0) { free(header); lastError = QString("wcspih ERROR %1: %2.").arg(status).arg(wcshdr_errmsg[status]); return false; } free(header); if (wcs == nullptr) { //fprintf(stderr, "No world coordinate systems found.\n"); lastError = i18n("No world coordinate systems found."); return false; } // FIXME: Call above goes through EVEN if no WCS is present, so we're adding this to return for now. if (wcs->crpix[0] == 0) { lastError = i18n("No world coordinate systems found."); return false; } if ((status = wcsset(wcs)) != 0) { lastError = QString("wcsset error %1: %2.").arg(status).arg(wcs_errmsg[status]); return false; } HasWCS = true; #endif #endif return HasWCS; } bool FITSData::loadWCS() { #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB) if (WCSLoaded) { qWarning() << "WCS data already loaded"; return true; } qCDebug(KSTARS_FITS) << "Started WCS Data Processing..."; int status = 0; char * header; int nkeyrec, nreject, nwcs, stat[2]; double imgcrd[2], phi = 0, pixcrd[2], theta = 0, world[2]; int w = width(); int h = height(); if (fits_hdr2str(fptr, 1, nullptr, 0, &header, &nkeyrec, &status)) { char errmsg[512]; fits_get_errstatus(status, errmsg); lastError = errmsg; return false; } if ((status = wcspih(header, nkeyrec, WCSHDR_all, -3, &nreject, &nwcs, &wcs)) != 0) { free(header); lastError = QString("wcspih ERROR %1: %2.").arg(status).arg(wcshdr_errmsg[status]); return false; } free(header); if (wcs == nullptr) { //fprintf(stderr, "No world coordinate systems found.\n"); lastError = i18n("No world coordinate systems found."); return false; } // FIXME: Call above goes through EVEN if no WCS is present, so we're adding this to return for now. if (wcs->crpix[0] == 0) { lastError = i18n("No world coordinate systems found."); return false; } if ((status = wcsset(wcs)) != 0) { lastError = QString("wcsset error %1: %2.").arg(status).arg(wcs_errmsg[status]); return false; } delete[] wcs_coord; wcs_coord = new wcs_point[w * h]; if (wcs_coord == nullptr) { lastError = "Not enough memory for WCS data!"; return false; } wcs_point * p = wcs_coord; for (int i = 0; i < h; i++) { for (int j = 0; j < w; j++) { pixcrd[0] = j; pixcrd[1] = i; if ((status = wcsp2s(wcs, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0])) != 0) { lastError = QString("wcsp2s error %1: %2.").arg(status).arg(wcs_errmsg[status]); } else { p->ra = world[0]; p->dec = world[1]; p++; } } } findObjectsInImage(&world[0], phi, theta, &imgcrd[0], &pixcrd[0], &stat[0]); WCSLoaded = true; HasWCS = true; qCDebug(KSTARS_FITS) << "Finished WCS Data processing..."; return true; #else return false; #endif } bool FITSData::wcsToPixel(SkyPoint &wcsCoord, QPointF &wcsPixelPoint, QPointF &wcsImagePoint) { #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB) int status = 0; int stat[2]; double imgcrd[2], worldcrd[2], pixcrd[2], phi[2], theta[2]; if (wcs == nullptr) { lastError = i18n("No world coordinate systems found."); return false; } worldcrd[0] = wcsCoord.ra0().Degrees(); worldcrd[1] = wcsCoord.dec0().Degrees(); if ((status = wcss2p(wcs, 1, 2, &worldcrd[0], &phi[0], &theta[0], &imgcrd[0], &pixcrd[0], &stat[0])) != 0) { lastError = QString("wcss2p error %1: %2.").arg(status).arg(wcs_errmsg[status]); return false; } wcsImagePoint.setX(imgcrd[0]); wcsImagePoint.setY(imgcrd[1]); wcsPixelPoint.setX(pixcrd[0]); wcsPixelPoint.setY(pixcrd[1]); return true; #else Q_UNUSED(wcsCoord); Q_UNUSED(wcsPixelPoint); Q_UNUSED(wcsImagePoint); return false; #endif } bool FITSData::pixelToWCS(const QPointF &wcsPixelPoint, SkyPoint &wcsCoord) { #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB) int status = 0; int stat[2]; double imgcrd[2], phi, pixcrd[2], theta, world[2]; if (wcs == nullptr) { lastError = i18n("No world coordinate systems found."); return false; } pixcrd[0] = wcsPixelPoint.x(); pixcrd[1] = wcsPixelPoint.y(); if ((status = wcsp2s(wcs, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0])) != 0) { lastError = QString("wcsp2s error %1: %2.").arg(status).arg(wcs_errmsg[status]); return false; } else { wcsCoord.setRA0(world[0] / 15.0); wcsCoord.setDec0(world[1]); } return true; #else Q_UNUSED(wcsPixelPoint); Q_UNUSED(wcsCoord); return false; #endif } #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB) void FITSData::findObjectsInImage(double world[], double phi, double theta, double imgcrd[], double pixcrd[], int stat[]) { int w = width(); int h = height(); int status = 0; char date[64]; KSNumbers * num = nullptr; if (fits_read_keyword(fptr, "DATE-OBS", date, nullptr, &status) == 0) { QString tsString(date); tsString = tsString.remove('\'').trimmed(); // Add Zulu time to indicate UTC tsString += "Z"; QDateTime ts = QDateTime::fromString(tsString, Qt::ISODate); if (ts.isValid()) num = new KSNumbers(KStarsDateTime(ts).djd()); } if (num == nullptr) num = new KSNumbers(KStarsData::Instance()->ut().djd()); //Set to current time if the above does not work. SkyMapComposite * map = KStarsData::Instance()->skyComposite(); wcs_point * wcs_coord = getWCSCoord(); if (wcs_coord != nullptr) { int size = w * h; objList.clear(); SkyPoint p1; p1.setRA0(dms(wcs_coord[0].ra)); p1.setDec0(dms(wcs_coord[0].dec)); p1.updateCoordsNow(num); SkyPoint p2; p2.setRA0(dms(wcs_coord[size - 1].ra)); p2.setDec0(dms(wcs_coord[size - 1].dec)); p2.updateCoordsNow(num); QList list = map->findObjectsInArea(p1, p2); foreach (SkyObject * object, list) { int type = object->type(); if (object->name() == "star" || type == SkyObject::PLANET || type == SkyObject::ASTEROID || type == SkyObject::COMET || type == SkyObject::SUPERNOVA || type == SkyObject::MOON || type == SkyObject::SATELLITE) { //DO NOT DISPLAY, at least for now, because these things move and change. } int x = -100; int y = -100; world[0] = object->ra0().Degrees(); world[1] = object->dec0().Degrees(); if ((status = wcss2p(wcs, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0])) != 0) { fprintf(stderr, "wcsp2s ERROR %d: %s.\n", status, wcs_errmsg[status]); } else { x = pixcrd[0]; //The X and Y are set to the found position if it does work. y = pixcrd[1]; } if (x > 0 && y > 0 && x < w && y < h) objList.append(new FITSSkyObject(object, x, y)); } } delete (num); } #endif QList FITSData::getSkyObjects() { return objList; } FITSSkyObject::FITSSkyObject(SkyObject * object, int xPos, int yPos) : QObject() { skyObjectStored = object; xLoc = xPos; yLoc = yPos; } SkyObject * FITSSkyObject::skyObject() { return skyObjectStored; } int FITSSkyObject::x() { return xLoc; } int FITSSkyObject::y() { return yLoc; } void FITSSkyObject::setX(int xPos) { xLoc = xPos; } void FITSSkyObject::setY(int yPos) { yLoc = yPos; } int FITSData::getFlipVCounter() const { return flipVCounter; } void FITSData::setFlipVCounter(int value) { flipVCounter = value; } int FITSData::getFlipHCounter() const { return flipHCounter; } void FITSData::setFlipHCounter(int value) { flipHCounter = value; } int FITSData::getRotCounter() const { return rotCounter; } void FITSData::setRotCounter(int value) { rotCounter = value; } /* Rotate an image by 90, 180, or 270 degrees, with an optional * reflection across the vertical or horizontal axis. * verbose generates extra info on stdout. * return nullptr if successful or rotated image. */ template bool FITSData::rotFITS(int rotate, int mirror) { int ny, nx; int x1, y1, x2, y2; uint8_t * rotimage = nullptr; int offset = 0; if (rotate == 1) rotate = 90; else if (rotate == 2) rotate = 180; else if (rotate == 3) rotate = 270; else if (rotate < 0) rotate = rotate + 360; nx = stats.width; ny = stats.height; int BBP = stats.bytesPerPixel; /* Allocate buffer for rotated image */ rotimage = new uint8_t[stats.samples_per_channel * m_Channels * BBP]; if (rotimage == nullptr) { qWarning() << "Unable to allocate memory for rotated image buffer!"; return false; } auto * rotBuffer = reinterpret_cast(rotimage); auto * buffer = reinterpret_cast(m_ImageBuffer); /* Mirror image without rotation */ if (rotate < 45 && rotate > -45) { if (mirror == 1) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (x1 = 0; x1 < nx; x1++) { x2 = nx - x1 - 1; for (y1 = 0; y1 < ny; y1++) rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } else if (mirror == 2) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { y2 = ny - y1 - 1; for (x1 = 0; x1 < nx; x1++) rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } else { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { for (x1 = 0; x1 < nx; x1++) rotBuffer[(y1 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } } /* Rotate by 90 degrees */ else if (rotate >= 45 && rotate < 135) { if (mirror == 1) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { x2 = ny - y1 - 1; for (x1 = 0; x1 < nx; x1++) { y2 = nx - x1 - 1; rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } } else if (mirror == 2) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { for (x1 = 0; x1 < nx; x1++) rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } else { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { x2 = ny - y1 - 1; for (x1 = 0; x1 < nx; x1++) { y2 = x1; rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } } stats.width = ny; stats.height = nx; } /* Rotate by 180 degrees */ else if (rotate >= 135 && rotate < 225) { if (mirror == 1) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { y2 = ny - y1 - 1; for (x1 = 0; x1 < nx; x1++) rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } else if (mirror == 2) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (x1 = 0; x1 < nx; x1++) { x2 = nx - x1 - 1; for (y1 = 0; y1 < ny; y1++) rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } else { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { y2 = ny - y1 - 1; for (x1 = 0; x1 < nx; x1++) { x2 = nx - x1 - 1; rotBuffer[(y2 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } } } /* Rotate by 270 degrees */ else if (rotate >= 225 && rotate < 315) { if (mirror == 1) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { for (x1 = 0; x1 < nx; x1++) rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } else if (mirror == 2) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { x2 = ny - y1 - 1; for (x1 = 0; x1 < nx; x1++) { y2 = nx - x1 - 1; rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } } else { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { x2 = y1; for (x1 = 0; x1 < nx; x1++) { y2 = nx - x1 - 1; rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } } stats.width = ny; stats.height = nx; } /* If rotating by more than 315 degrees, assume top-bottom reflection */ else if (rotate >= 315 && mirror) { for (int i = 0; i < m_Channels; i++) { offset = stats.samples_per_channel * i; for (y1 = 0; y1 < ny; y1++) { for (x1 = 0; x1 < nx; x1++) { x2 = y1; y2 = x1; rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset]; } } } } delete[] m_ImageBuffer; m_ImageBuffer = rotimage; return true; } void FITSData::rotWCSFITS(int angle, int mirror) { int status = 0; char comment[100]; double ctemp1, ctemp2, ctemp3, ctemp4, naxis1, naxis2; int WCS_DECIMALS = 6; naxis1 = stats.width; naxis2 = stats.height; if (fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status)) { // No WCS keywords return; } /* Reset CROTAn and CD matrix if axes have been exchanged */ if (angle == 90) { if (!fits_read_key_dbl(fptr, "CROTA1", &ctemp1, comment, &status)) fits_update_key_dbl(fptr, "CROTA1", ctemp1 + 90.0, WCS_DECIMALS, comment, &status); if (!fits_read_key_dbl(fptr, "CROTA2", &ctemp1, comment, &status)) fits_update_key_dbl(fptr, "CROTA2", ctemp1 + 90.0, WCS_DECIMALS, comment, &status); } status = 0; /* Negate rotation angle if mirrored */ if (mirror != 0) { if (!fits_read_key_dbl(fptr, "CROTA1", &ctemp1, comment, &status)) fits_update_key_dbl(fptr, "CROTA1", -ctemp1, WCS_DECIMALS, comment, &status); if (!fits_read_key_dbl(fptr, "CROTA2", &ctemp1, comment, &status)) fits_update_key_dbl(fptr, "CROTA2", -ctemp1, WCS_DECIMALS, comment, &status); status = 0; if (!fits_read_key_dbl(fptr, "LTM1_1", &ctemp1, comment, &status)) fits_update_key_dbl(fptr, "LTM1_1", -ctemp1, WCS_DECIMALS, comment, &status); status = 0; if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status)) fits_update_key_dbl(fptr, "CD1_1", -ctemp1, WCS_DECIMALS, comment, &status); if (!fits_read_key_dbl(fptr, "CD1_2", &ctemp1, comment, &status)) fits_update_key_dbl(fptr, "CD1_2", -ctemp1, WCS_DECIMALS, comment, &status); if (!fits_read_key_dbl(fptr, "CD2_1", &ctemp1, comment, &status)) fits_update_key_dbl(fptr, "CD2_1", -ctemp1, WCS_DECIMALS, comment, &status); } status = 0; /* Unbin CRPIX and CD matrix */ if (!fits_read_key_dbl(fptr, "LTM1_1", &ctemp1, comment, &status)) { if (ctemp1 != 1.0) { if (!fits_read_key_dbl(fptr, "LTM2_2", &ctemp2, comment, &status)) if (ctemp1 == ctemp2) { double ltv1 = 0.0; double ltv2 = 0.0; status = 0; if (!fits_read_key_dbl(fptr, "LTV1", <v1, comment, &status)) fits_delete_key(fptr, "LTV1", &status); if (!fits_read_key_dbl(fptr, "LTV2", <v2, comment, &status)) fits_delete_key(fptr, "LTV2", &status); status = 0; if (!fits_read_key_dbl(fptr, "CRPIX1", &ctemp3, comment, &status)) fits_update_key_dbl(fptr, "CRPIX1", (ctemp3 - ltv1) / ctemp1, WCS_DECIMALS, comment, &status); if (!fits_read_key_dbl(fptr, "CRPIX2", &ctemp3, comment, &status)) fits_update_key_dbl(fptr, "CRPIX2", (ctemp3 - ltv2) / ctemp1, WCS_DECIMALS, comment, &status); status = 0; if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp3, comment, &status)) fits_update_key_dbl(fptr, "CD1_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status); if (!fits_read_key_dbl(fptr, "CD1_2", &ctemp3, comment, &status)) fits_update_key_dbl(fptr, "CD1_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status); if (!fits_read_key_dbl(fptr, "CD2_1", &ctemp3, comment, &status)) fits_update_key_dbl(fptr, "CD2_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status); if (!fits_read_key_dbl(fptr, "CD2_2", &ctemp3, comment, &status)) fits_update_key_dbl(fptr, "CD2_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status); status = 0; fits_delete_key(fptr, "LTM1_1", &status); fits_delete_key(fptr, "LTM1_2", &status); } } } status = 0; /* Reset CRPIXn */ if (!fits_read_key_dbl(fptr, "CRPIX1", &ctemp1, comment, &status) && !fits_read_key_dbl(fptr, "CRPIX2", &ctemp2, comment, &status)) { if (mirror != 0) { if (angle == 0) fits_update_key_dbl(fptr, "CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status); else if (angle == 90) { fits_update_key_dbl(fptr, "CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status); } else if (angle == 180) { fits_update_key_dbl(fptr, "CRPIX1", ctemp1, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status); } else if (angle == 270) { fits_update_key_dbl(fptr, "CRPIX1", ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CRPIX2", ctemp1, WCS_DECIMALS, comment, &status); } } else { if (angle == 90) { fits_update_key_dbl(fptr, "CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CRPIX2", ctemp1, WCS_DECIMALS, comment, &status); } else if (angle == 180) { fits_update_key_dbl(fptr, "CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status); } else if (angle == 270) { fits_update_key_dbl(fptr, "CRPIX1", ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status); } } } status = 0; /* Reset CDELTn (degrees per pixel) */ if (!fits_read_key_dbl(fptr, "CDELT1", &ctemp1, comment, &status) && !fits_read_key_dbl(fptr, "CDELT2", &ctemp2, comment, &status)) { if (mirror != 0) { if (angle == 0) fits_update_key_dbl(fptr, "CDELT1", -ctemp1, WCS_DECIMALS, comment, &status); else if (angle == 90) { fits_update_key_dbl(fptr, "CDELT1", -ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CDELT2", -ctemp1, WCS_DECIMALS, comment, &status); } else if (angle == 180) { fits_update_key_dbl(fptr, "CDELT1", ctemp1, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CDELT2", -ctemp2, WCS_DECIMALS, comment, &status); } else if (angle == 270) { fits_update_key_dbl(fptr, "CDELT1", ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CDELT2", ctemp1, WCS_DECIMALS, comment, &status); } } else { if (angle == 90) { fits_update_key_dbl(fptr, "CDELT1", -ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CDELT2", ctemp1, WCS_DECIMALS, comment, &status); } else if (angle == 180) { fits_update_key_dbl(fptr, "CDELT1", -ctemp1, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CDELT2", -ctemp2, WCS_DECIMALS, comment, &status); } else if (angle == 270) { fits_update_key_dbl(fptr, "CDELT1", ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CDELT2", -ctemp1, WCS_DECIMALS, comment, &status); } } } /* Reset CD matrix, if present */ ctemp1 = 0.0; ctemp2 = 0.0; ctemp3 = 0.0; ctemp4 = 0.0; status = 0; if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status)) { fits_read_key_dbl(fptr, "CD1_2", &ctemp2, comment, &status); fits_read_key_dbl(fptr, "CD2_1", &ctemp3, comment, &status); fits_read_key_dbl(fptr, "CD2_2", &ctemp4, comment, &status); status = 0; if (mirror != 0) { if (angle == 0) { fits_update_key_dbl(fptr, "CD1_2", -ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status); } else if (angle == 90) { fits_update_key_dbl(fptr, "CD1_1", -ctemp4, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD1_2", -ctemp3, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_1", -ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_2", -ctemp1, WCS_DECIMALS, comment, &status); } else if (angle == 180) { fits_update_key_dbl(fptr, "CD1_1", ctemp1, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD1_2", ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_2", -ctemp4, WCS_DECIMALS, comment, &status); } else if (angle == 270) { fits_update_key_dbl(fptr, "CD1_1", ctemp4, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD1_2", ctemp3, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_1", ctemp3, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_2", ctemp1, WCS_DECIMALS, comment, &status); } } else { if (angle == 90) { fits_update_key_dbl(fptr, "CD1_1", -ctemp4, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD1_2", -ctemp3, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_1", ctemp1, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_2", ctemp1, WCS_DECIMALS, comment, &status); } else if (angle == 180) { fits_update_key_dbl(fptr, "CD1_1", -ctemp1, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD1_2", -ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_2", -ctemp4, WCS_DECIMALS, comment, &status); } else if (angle == 270) { fits_update_key_dbl(fptr, "CD1_1", ctemp4, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD1_2", ctemp3, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_1", -ctemp2, WCS_DECIMALS, comment, &status); fits_update_key_dbl(fptr, "CD2_2", -ctemp1, WCS_DECIMALS, comment, &status); } } } /* Delete any polynomial solution */ /* (These could maybe be switched, but I don't want to work them out yet */ status = 0; if (!fits_read_key_dbl(fptr, "CO1_1", &ctemp1, comment, &status)) { int i; char keyword[16]; for (i = 1; i < 13; i++) { sprintf(keyword, "CO1_%d", i); fits_delete_key(fptr, keyword, &status); } for (i = 1; i < 13; i++) { sprintf(keyword, "CO2_%d", i); fits_delete_key(fptr, keyword, &status); } } } uint8_t * FITSData::getImageBuffer() { return m_ImageBuffer; } void FITSData::setImageBuffer(uint8_t * buffer) { delete[] m_ImageBuffer; m_ImageBuffer = buffer; } bool FITSData::checkDebayer() { int status = 0; char bayerPattern[64]; // Let's search for BAYERPAT keyword, if it's not found we return as there is no bayer pattern in this image if (fits_read_keyword(fptr, "BAYERPAT", bayerPattern, nullptr, &status)) return false; if (stats.bitpix != 16 && stats.bitpix != 8) { KSNotification::error(i18n("Only 8 and 16 bits bayered images supported."), i18n("Debayer error")); return false; } QString pattern(bayerPattern); pattern = pattern.remove('\'').trimmed(); if (pattern == "RGGB") debayerParams.filter = DC1394_COLOR_FILTER_RGGB; else if (pattern == "GBRG") debayerParams.filter = DC1394_COLOR_FILTER_GBRG; else if (pattern == "GRBG") debayerParams.filter = DC1394_COLOR_FILTER_GRBG; else if (pattern == "BGGR") debayerParams.filter = DC1394_COLOR_FILTER_BGGR; // We return unless we find a valid pattern else { KSNotification::error(i18n("Unsupported bayer pattern %1.", pattern), i18n("Debayer error")); return false; } fits_read_key(fptr, TINT, "XBAYROFF", &debayerParams.offsetX, nullptr, &status); fits_read_key(fptr, TINT, "YBAYROFF", &debayerParams.offsetY, nullptr, &status); HasDebayer = true; return true; } void FITSData::getBayerParams(BayerParams * param) { param->method = debayerParams.method; param->filter = debayerParams.filter; param->offsetX = debayerParams.offsetX; param->offsetY = debayerParams.offsetY; } void FITSData::setBayerParams(BayerParams * param) { debayerParams.method = param->method; debayerParams.filter = param->filter; debayerParams.offsetX = param->offsetX; debayerParams.offsetY = param->offsetY; } bool FITSData::debayer() { if (bayerBuffer == nullptr) { int anynull = 0, status = 0; bayerBuffer = m_ImageBuffer; if (fits_read_img(fptr, m_DataType, 1, stats.samples_per_channel, nullptr, bayerBuffer, &anynull, &status)) { char errmsg[512]; fits_get_errstatus(status, errmsg); KSNotification::error(i18n("Error reading image: %1", QString(errmsg)), i18n("Debayer error")); return false; } } switch (m_DataType) { case TBYTE: return debayer_8bit(); case TUSHORT: return debayer_16bit(); default: return false; } } bool FITSData::debayer_8bit() { dc1394error_t error_code; int rgb_size = stats.samples_per_channel * 3 * stats.bytesPerPixel; auto * destinationBuffer = new uint8_t[rgb_size]; if (destinationBuffer == nullptr) { KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer."), i18n("Debayer error")); return false; } int ds1394_height = stats.height; uint8_t * dc1394_source = bayerBuffer; if (debayerParams.offsetY == 1) { dc1394_source += stats.width; ds1394_height--; } if (debayerParams.offsetX == 1) { dc1394_source++; } error_code = dc1394_bayer_decoding_8bit(dc1394_source, destinationBuffer, stats.width, ds1394_height, debayerParams.filter, debayerParams.method); if (error_code != DC1394_SUCCESS) { KSNotification::error(i18n("Debayer failed (%1)", error_code), i18n("Debayer error")); m_Channels = 1; delete[] destinationBuffer; return false; } if (m_Channels == 1) { delete[] m_ImageBuffer; m_ImageBuffer = new uint8_t[rgb_size]; if (m_ImageBuffer == nullptr) { delete[] destinationBuffer; KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer."), i18n("Debayer error")); return false; } } // Data in R1G1B1, we need to copy them into 3 layers for FITS uint8_t * rBuff = m_ImageBuffer; uint8_t * gBuff = m_ImageBuffer + (stats.width * stats.height); uint8_t * bBuff = m_ImageBuffer + (stats.width * stats.height * 2); int imax = stats.samples_per_channel * 3 - 3; for (int i = 0; i <= imax; i += 3) { *rBuff++ = destinationBuffer[i]; *gBuff++ = destinationBuffer[i + 1]; *bBuff++ = destinationBuffer[i + 2]; } m_Channels = (m_Mode == FITS_NORMAL) ? 3 : 1; delete[] destinationBuffer; bayerBuffer = nullptr; return true; } bool FITSData::debayer_16bit() { dc1394error_t error_code; int rgb_size = stats.samples_per_channel * 3 * stats.bytesPerPixel; auto * destinationBuffer = new uint8_t[rgb_size]; auto * buffer = reinterpret_cast(bayerBuffer); auto * dstBuffer = reinterpret_cast(destinationBuffer); if (destinationBuffer == nullptr) { KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer."), i18n("Debayer error")); return false; } int ds1394_height = stats.height; uint16_t * dc1394_source = buffer; if (debayerParams.offsetY == 1) { dc1394_source += stats.width; ds1394_height--; } if (debayerParams.offsetX == 1) { dc1394_source++; } error_code = dc1394_bayer_decoding_16bit(dc1394_source, dstBuffer, stats.width, ds1394_height, debayerParams.filter, debayerParams.method, 16); if (error_code != DC1394_SUCCESS) { KSNotification::error(i18n("Debayer failed (%1)", error_code), i18n("Debayer error")); m_Channels = 1; delete[] destinationBuffer; return false; } if (m_Channels == 1) { delete[] m_ImageBuffer; m_ImageBuffer = new uint8_t[rgb_size]; if (m_ImageBuffer == nullptr) { delete[] destinationBuffer; KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer."), i18n("Debayer error")); return false; } } buffer = reinterpret_cast(m_ImageBuffer); // Data in R1G1B1, we need to copy them into 3 layers for FITS uint16_t * rBuff = buffer; uint16_t * gBuff = buffer + (stats.width * stats.height); uint16_t * bBuff = buffer + (stats.width * stats.height * 2); int imax = stats.samples_per_channel * 3 - 3; for (int i = 0; i <= imax; i += 3) { *rBuff++ = dstBuffer[i]; *gBuff++ = dstBuffer[i + 1]; *bBuff++ = dstBuffer[i + 2]; } m_Channels = (m_Mode == FITS_NORMAL) ? 3 : 1; delete[] destinationBuffer; bayerBuffer = nullptr; return true; } double FITSData::getADU() const { double adu = 0; for (int i = 0; i < m_Channels; i++) adu += stats.mean[i]; return (adu / static_cast(m_Channels)); } /* CannyDetector, Implementation of Canny edge detector in Qt/C++. * Copyright (C) 2015 Gonzalo Exequiel Pedone * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see . * * Email : hipersayan DOT x AT gmail DOT com * Web-Site: http://github.com/hipersayanX/CannyDetector */ -#if 0 -void FITSData::sobel(const QImage &image, QVector &gradient, QVector &direction) -{ - int size = image.width() * image.height(); - gradient.resize(size); - direction.resize(size); - - for (int y = 0; y < image.height(); y++) - { - size_t yOffset = y * image.width(); - const quint8 * grayLine = image.constBits() + yOffset; - - const quint8 * grayLine_m1 = y < 1 ? grayLine : grayLine - image.width(); - const quint8 * grayLine_p1 = y >= image.height() - 1 ? grayLine : grayLine + image.width(); - - int * gradientLine = gradient.data() + yOffset; - int * directionLine = direction.data() + yOffset; - - for (int x = 0; x < image.width(); x++) - { - int x_m1 = x < 1 ? x : x - 1; - int x_p1 = x >= image.width() - 1 ? x : x + 1; - - int gradX = grayLine_m1[x_p1] - + 2 * grayLine[x_p1] - + grayLine_p1[x_p1] - - grayLine_m1[x_m1] - - 2 * grayLine[x_m1] - - grayLine_p1[x_m1]; - - int gradY = grayLine_m1[x_m1] - + 2 * grayLine_m1[x] - + grayLine_m1[x_p1] - - grayLine_p1[x_m1] - - 2 * grayLine_p1[x] - - grayLine_p1[x_p1]; - - gradientLine[x] = qAbs(gradX) + qAbs(gradY); - - /* Gradient directions are classified in 4 possible cases - * - * dir 0 - * - * x x x - * - - - - * x x x - * - * dir 1 - * - * x x / - * x / x - * / x x - * - * dir 2 - * - * \ x x - * x \ x - * x x \ - * - * dir 3 - * - * x | x - * x | x - * x | x - */ - if (gradX == 0 && gradY == 0) - directionLine[x] = 0; - else if (gradX == 0) - directionLine[x] = 3; - else - { - qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI; - - if (a >= -22.5 && a < 22.5) - directionLine[x] = 0; - else if (a >= 22.5 && a < 67.5) - directionLine[x] = 1; - else if (a >= -67.5 && a < -22.5) - directionLine[x] = 2; - else - directionLine[x] = 3; - } - } - } -} -#endif - template void FITSData::sobel(QVector &gradient, QVector &direction) { //int size = image.width() * image.height(); gradient.resize(stats.samples_per_channel); direction.resize(stats.samples_per_channel); for (int y = 0; y < stats.height; y++) { size_t yOffset = y * stats.width; const T * grayLine = reinterpret_cast(m_ImageBuffer) + yOffset; const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.width; const T * grayLine_p1 = y >= stats.height - 1 ? grayLine : grayLine + stats.width; float * gradientLine = gradient.data() + yOffset; float * directionLine = direction.data() + yOffset; for (int x = 0; x < stats.width; x++) { int x_m1 = x < 1 ? x : x - 1; int x_p1 = x >= stats.width - 1 ? x : x + 1; int gradX = grayLine_m1[x_p1] + 2 * grayLine[x_p1] + grayLine_p1[x_p1] - grayLine_m1[x_m1] - 2 * grayLine[x_m1] - grayLine_p1[x_m1]; int gradY = grayLine_m1[x_m1] + 2 * grayLine_m1[x] + grayLine_m1[x_p1] - grayLine_p1[x_m1] - 2 * grayLine_p1[x] - grayLine_p1[x_p1]; gradientLine[x] = qAbs(gradX) + qAbs(gradY); /* Gradient directions are classified in 4 possible cases * * dir 0 * * x x x * - - - * x x x * * dir 1 * * x x / * x / x * / x x * * dir 2 * * \ x x * x \ x * x x \ * * dir 3 * * x | x * x | x * x | x */ if (gradX == 0 && gradY == 0) directionLine[x] = 0; else if (gradX == 0) directionLine[x] = 3; else { qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI; if (a >= -22.5 && a < 22.5) directionLine[x] = 0; else if (a >= 22.5 && a < 67.5) directionLine[x] = 2; else if (a >= -67.5 && a < -22.5) directionLine[x] = 1; else directionLine[x] = 3; } } } } int FITSData::partition(int width, int height, QVector &gradient, QVector &ids) { int id = 0; for (int y = 1; y < height - 1; y++) { for (int x = 1; x < width - 1; x++) { int index = x + y * width; float val = gradient[index]; if (val > 0 && ids[index] == 0) { trace(width, height, ++id, gradient, ids, x, y); } } } // Return max id return id; } void FITSData::trace(int width, int height, int id, QVector &image, QVector &ids, int x, int y) { int yOffset = y * width; float * cannyLine = image.data() + yOffset; int * idLine = ids.data() + yOffset; if (idLine[x] != 0) return; idLine[x] = id; for (int j = -1; j < 2; j++) { int nextY = y + j; if (nextY < 0 || nextY >= height) continue; float * cannyLineNext = cannyLine + j * width; for (int i = -1; i < 2; i++) { int nextX = x + i; if (i == j || nextX < 0 || nextX >= width) continue; if (cannyLineNext[nextX] > 0) { // Trace neighbors. trace(width, height, id, image, ids, nextX, nextY); } } } } QString FITSData::getLastError() const { return lastError; } bool FITSData::getAutoRemoveTemporaryFITS() const { return autoRemoveTemporaryFITS; } void FITSData::setAutoRemoveTemporaryFITS(bool value) { autoRemoveTemporaryFITS = value; } -#if 0 -QVector FITSData::thinning(int width, int height, const QVector &gradient, const QVector &direction) -{ - QVector thinned(gradient.size()); - - for (int y = 0; y < height; y++) - { - int yOffset = y * width; - const int * gradientLine = gradient.constData() + yOffset; - const int * gradientLine_m1 = y < 1 ? gradientLine : gradientLine - width; - const int * gradientLine_p1 = y >= height - 1 ? gradientLine : gradientLine + width; - const int * directionLine = direction.constData() + yOffset; - int * thinnedLine = thinned.data() + yOffset; - - for (int x = 0; x < width; x++) - { - int x_m1 = x < 1 ? 0 : x - 1; - int x_p1 = x >= width - 1 ? x : x + 1; - - int direction = directionLine[x]; - int pixel = 0; - - if (direction == 0) - { - /* x x x - * - - - - * x x x - */ - if (gradientLine[x] < gradientLine[x_m1] - || gradientLine[x] < gradientLine[x_p1]) - pixel = 0; - else - pixel = gradientLine[x]; - } - else if (direction == 1) - { - /* x x / - * x / x - * / x x - */ - if (gradientLine[x] < gradientLine_m1[x_p1] - || gradientLine[x] < gradientLine_p1[x_m1]) - pixel = 0; - else - pixel = gradientLine[x]; - } - else if (direction == 2) - { - /* \ x x - * x \ x - * x x \ - */ - if (gradientLine[x] < gradientLine_m1[x_m1] - || gradientLine[x] < gradientLine_p1[x_p1]) - pixel = 0; - else - pixel = gradientLine[x]; - } - else - { - /* x | x - * x | x - * x | x - */ - if (gradientLine[x] < gradientLine_m1[x] - || gradientLine[x] < gradientLine_p1[x]) - pixel = 0; - else - pixel = gradientLine[x]; - } - - thinnedLine[x] = pixel; - } - } - - return thinned; -} - -QVector FITSData::threshold(int thLow, int thHi, const QVector &image) -{ - QVector thresholded(image.size()); - - for (int i = 0; i < image.size(); i++) - thresholded[i] = image[i] <= thLow ? 0 : - image[i] >= thHi ? 255 : - 127; - - return thresholded; -} - -QVector FITSData::hysteresis(int width, int height, const QVector &image) -{ - QVector canny(image); - - for (int y = 0; y < height; y++) - for (int x = 0; x < width; x++) - trace(width, height, canny, x, y); - - // Remaining gray pixels becomes black. - for (int i = 0; i < canny.size(); i++) - if (canny[i] == 127) - canny[i] = 0; - - return canny; -} - -#endif template void FITSData::convertToQImage(double dataMin, double dataMax, double scale, double zero, QImage &image) { #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wcast-align" auto * buffer = (T *)getImageBuffer(); #pragma GCC diagnostic pop const T limit = std::numeric_limits::max(); T bMin = dataMin < 0 ? 0 : dataMin; T bMax = dataMax > limit ? limit : dataMax; uint16_t w = width(); uint16_t h = height(); uint32_t size = w * h; double val; if (channels() == 1) { /* Fill in pixel values using indexed map, linear scale */ for (int j = 0; j < h; j++) { unsigned char * scanLine = image.scanLine(j); for (int i = 0; i < w; i++) { val = qBound(bMin, buffer[j * w + i], bMax); val = val * scale + zero; scanLine[i] = qBound(0, (unsigned char)val, 255); } } } else { double rval = 0, gval = 0, bval = 0; QRgb value; /* Fill in pixel values using indexed map, linear scale */ for (int j = 0; j < h; j++) { auto * scanLine = reinterpret_cast((image.scanLine(j))); for (int i = 0; i < w; i++) { rval = qBound(bMin, buffer[j * w + i], bMax); gval = qBound(bMin, buffer[j * w + i + size], bMax); bval = qBound(bMin, buffer[j * w + i + size * 2], bMax); value = qRgb(rval * scale + zero, gval * scale + zero, bval * scale + zero); scanLine[i] = value; } } } } QImage FITSData::FITSToImage(const QString &filename) { QImage fitsImage; double min, max; FITSData data; QFuture future = data.loadFITS(filename); // Wait synchronously future.waitForFinished(); if (future.result() == false) return fitsImage; data.getMinMax(&min, &max); if (min == max) { fitsImage.fill(Qt::white); return fitsImage; } if (data.channels() == 1) { fitsImage = QImage(data.width(), data.height(), QImage::Format_Indexed8); fitsImage.setColorCount(256); for (int i = 0; i < 256; i++) fitsImage.setColor(i, qRgb(i, i, i)); } else { fitsImage = QImage(data.width(), data.height(), QImage::Format_RGB32); } double dataMin = data.stats.mean[0] - data.stats.stddev[0]; double dataMax = data.stats.mean[0] + data.stats.stddev[0] * 3; double bscale = 255. / (dataMax - dataMin); double bzero = (-dataMin) * (255. / (dataMax - dataMin)); // Long way to do this since we do not want to use templated functions here switch (data.property("dataType").toInt()) { case TBYTE: data.convertToQImage(dataMin, dataMax, bscale, bzero, fitsImage); break; case TSHORT: data.convertToQImage(dataMin, dataMax, bscale, bzero, fitsImage); break; case TUSHORT: data.convertToQImage(dataMin, dataMax, bscale, bzero, fitsImage); break; case TLONG: data.convertToQImage(dataMin, dataMax, bscale, bzero, fitsImage); break; case TULONG: data.convertToQImage(dataMin, dataMax, bscale, bzero, fitsImage); break; case TFLOAT: data.convertToQImage(dataMin, dataMax, bscale, bzero, fitsImage); break; case TLONGLONG: data.convertToQImage(dataMin, dataMax, bscale, bzero, fitsImage); break; case TDOUBLE: data.convertToQImage(dataMin, dataMax, bscale, bzero, fitsImage); break; default: break; } return fitsImage; } bool FITSData::ImageToFITS(const QString &filename, const QString &format, QString &output) { if (QImageReader::supportedImageFormats().contains(format.toLatin1()) == false) { qCCritical(KSTARS_FITS) << "Failed to convert" << filename << "to FITS since format" << format << "is not supported in Qt"; return false; } QImage input; if (input.load(filename, format.toLatin1()) == false) { qCCritical(KSTARS_FITS) << "Failed to open image" << filename; return false; } output = QString(KSPaths::writableLocation(QStandardPaths::TempLocation) + QFileInfo(filename).fileName() + ".fits"); //This section sets up the FITS File fitsfile *fptr = nullptr; int status = 0; long fpixel = 1, naxis = input.allGray() ? 2 : 3, nelements, exposure; long naxes[3] = { input.width(), input.height(), naxis == 3 ? 3 : 1 }; char error_status[512] = {0}; if (fits_create_file(&fptr, QString('!' + output).toLatin1().data(), &status)) { qCCritical(KSTARS_FITS) << "Failed to create FITS file. Error:" << status; return false; } if (fits_create_img(fptr, BYTE_IMG, naxis, naxes, &status)) { qCWarning(KSTARS_FITS) << "fits_create_img failed:" << error_status; status = 0; fits_close_file(fptr, &status); return false; } exposure = 1; fits_update_key(fptr, TLONG, "EXPOSURE", &exposure, "Total Exposure Time", &status); // Gray image if (naxis == 2) { nelements = naxes[0] * naxes[1]; if (fits_write_img(fptr, TBYTE, fpixel, nelements, input.bits(), &status)) { fits_get_errstatus(status, error_status); qCWarning(KSTARS_FITS) << "fits_write_img GRAY failed:" << error_status; status = 0; fits_close_file(fptr, &status); return false; } } // RGB image, we have to convert from ARGB format to R G B for each plane else { nelements = naxes[0] * naxes[1] * 3; uint8_t *srcBuffer = input.bits(); // ARGB uint32_t srcBytes = naxes[0] * naxes[1] * 4 - 4; uint8_t *rgbBuffer = new uint8_t[nelements]; if (rgbBuffer == nullptr) { qCWarning(KSTARS_FITS) << "Not enough memory for RGB buffer"; fits_close_file(fptr, &status); return false; } uint8_t *subR = rgbBuffer; uint8_t *subG = rgbBuffer + naxes[0] * naxes[1]; uint8_t *subB = rgbBuffer + naxes[0] * naxes[1] * 2; for (uint32_t i = 0; i < srcBytes; i += 4) { *subB++ = srcBuffer[i]; *subG++ = srcBuffer[i + 1]; *subR++ = srcBuffer[i + 2]; } if (fits_write_img(fptr, TBYTE, fpixel, nelements, rgbBuffer, &status)) { fits_get_errstatus(status, error_status); qCWarning(KSTARS_FITS) << "fits_write_img RGB failed:" << error_status; status = 0; fits_close_file(fptr, &status); delete [] rgbBuffer; return false; } delete [] rgbBuffer; } if (fits_flush_file(fptr, &status)) { fits_get_errstatus(status, error_status); qCWarning(KSTARS_FITS) << "fits_flush_file failed:" << error_status; status = 0; fits_close_file(fptr, &status); return false; } if (fits_close_file(fptr, &status)) { fits_get_errstatus(status, error_status); qCWarning(KSTARS_FITS) << "fits_close_file failed:" << error_status; return false; } return true; } bool FITSData::createWCSFile(const QString &newWCSFile, double orientation, double ra, double dec, double pixscale) { int status = 0, exttype = 0; long nelements; fitsfile * new_fptr; char errMsg[512]; qCInfo(KSTARS_FITS) << "Creating new WCS file:" << newWCSFile << "with parameters Orientation:" << orientation << "RA:" << ra << "DE:" << dec << "Pixel Scale:" << pixscale; nelements = stats.samples_per_channel * m_Channels; /* Create a new File, overwriting existing*/ if (fits_create_file(&new_fptr, QString('!' + newWCSFile).toLatin1(), &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } if (fits_movabs_hdu(fptr, 1, &exttype, &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } if (fits_copy_file(fptr, new_fptr, 1, 1, 1, &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } /* close current file */ if (fits_close_file(fptr, &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } status = 0; if (m_isTemporary && autoRemoveTemporaryFITS) { QFile::remove(m_Filename); m_isTemporary = false; qCDebug(KSTARS_FITS) << "Removing FITS File: " << m_Filename; } m_Filename = newWCSFile; m_isTemporary = true; fptr = new_fptr; if (fits_movabs_hdu(fptr, 1, &exttype, &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } /* Write Data */ if (fits_write_img(fptr, m_DataType, 1, nelements, m_ImageBuffer, &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } /* Write keywords */ // Minimum if (fits_update_key(fptr, TDOUBLE, "DATAMIN", &(stats.min), "Minimum value", &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } // Maximum if (fits_update_key(fptr, TDOUBLE, "DATAMAX", &(stats.max), "Maximum value", &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } // NAXIS1 if (fits_update_key(fptr, TUSHORT, "NAXIS1", &(stats.width), "length of data axis 1", &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } // NAXIS2 if (fits_update_key(fptr, TUSHORT, "NAXIS2", &(stats.height), "length of data axis 2", &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } fits_update_key(fptr, TDOUBLE, "OBJCTRA", &ra, "Object RA", &status); fits_update_key(fptr, TDOUBLE, "OBJCTDEC", &dec, "Object DEC", &status); int epoch = 2000; fits_update_key(fptr, TINT, "EQUINOX", &epoch, "Equinox", &status); fits_update_key(fptr, TDOUBLE, "CRVAL1", &ra, "CRVAL1", &status); fits_update_key(fptr, TDOUBLE, "CRVAL2", &dec, "CRVAL1", &status); char radecsys[8] = "FK5"; char ctype1[16] = "RA---TAN"; char ctype2[16] = "DEC--TAN"; fits_update_key(fptr, TSTRING, "RADECSYS", radecsys, "RADECSYS", &status); fits_update_key(fptr, TSTRING, "CTYPE1", ctype1, "CTYPE1", &status); fits_update_key(fptr, TSTRING, "CTYPE2", ctype2, "CTYPE2", &status); double crpix1 = width() / 2.0; double crpix2 = height() / 2.0; fits_update_key(fptr, TDOUBLE, "CRPIX1", &crpix1, "CRPIX1", &status); fits_update_key(fptr, TDOUBLE, "CRPIX2", &crpix2, "CRPIX2", &status); // Arcsecs per Pixel double secpix1 = pixscale; double secpix2 = pixscale; fits_update_key(fptr, TDOUBLE, "SECPIX1", &secpix1, "SECPIX1", &status); fits_update_key(fptr, TDOUBLE, "SECPIX2", &secpix2, "SECPIX2", &status); double degpix1 = secpix1 / 3600.0; double degpix2 = secpix2 / 3600.0; fits_update_key(fptr, TDOUBLE, "CDELT1", °pix1, "CDELT1", &status); fits_update_key(fptr, TDOUBLE, "CDELT2", °pix2, "CDELT2", &status); // Rotation is CW, we need to convert it to CCW per CROTA1 definition double rotation = 360 - orientation; if (rotation > 360) rotation -= 360; fits_update_key(fptr, TDOUBLE, "CROTA1", &rotation, "CROTA1", &status); fits_update_key(fptr, TDOUBLE, "CROTA2", &rotation, "CROTA2", &status); // ISO Date if (fits_write_date(fptr, &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } QString history = QString("Modified by KStars on %1").arg(QDateTime::currentDateTime().toString("yyyy-MM-ddThh:mm:ss")); // History if (fits_write_history(fptr, history.toLatin1(), &status)) { fits_get_errstatus(status, errMsg); lastError = QString(errMsg); fits_report_error(stderr, status); return false; } fits_flush_file(fptr, &status); WCSLoaded = false; qCDebug(KSTARS_FITS) << "Finished creating WCS file: " << newWCSFile; return true; } bool FITSData::contains(const QPointF &point) const { return (point.x() >= 0 && point.y() >= 0 && point.x() <= stats.width && point.y() <= stats.height); } int FITSData::findSEPStars(const QRect &boundary) { int x = 0, y = 0, w = stats.width, h = stats.height, maxRadius = 50; if (!boundary.isNull()) { x = boundary.x(); y = boundary.y(); w = boundary.width(); h = boundary.height(); maxRadius = w; } auto * data = new float[w * h]; switch (stats.bitpix) { case BYTE_IMG: getFloatBuffer(data, x, y, w, h); break; case SHORT_IMG: getFloatBuffer(data, x, y, w, h); break; case USHORT_IMG: getFloatBuffer(data, x, y, w, h); break; case LONG_IMG: getFloatBuffer(data, x, y, w, h); break; case ULONG_IMG: getFloatBuffer(data, x, y, w, h); break; case FLOAT_IMG: delete [] data; data = reinterpret_cast(m_ImageBuffer); break; case LONGLONG_IMG: getFloatBuffer(data, x, y, w, h); break; case DOUBLE_IMG: getFloatBuffer(data, x, y, w, h); break; default: delete [] data; return -1; } float * imback = nullptr; double * flux = nullptr, *fluxerr = nullptr, *area = nullptr; short * flag = nullptr; short flux_flag = 0; int status = 0; sep_bkg * bkg = nullptr; sep_catalog * catalog = nullptr; float conv[] = {1, 2, 1, 2, 4, 2, 1, 2, 1}; double flux_fractions[2] = {0}; double requested_frac[2] = { 0.5, 0.99 }; QList edges; // #0 Create SEP Image structure sep_image im = {data, nullptr, nullptr, SEP_TFLOAT, 0, 0, w, h, 0.0, SEP_NOISE_NONE, 1.0, 0.0}; // #1 Background estimate status = sep_background(&im, 64, 64, 3, 3, 0.0, &bkg); if (status != 0) goto exit; // #2 Background evaluation imback = (float *)malloc((w * h) * sizeof(float)); status = sep_bkg_array(bkg, imback, SEP_TFLOAT); if (status != 0) goto exit; // #3 Background subtraction status = sep_bkg_subarray(bkg, im.data, im.dtype); if (status != 0) goto exit; // #4 Source Extraction // Note that we set deblend_cont = 1.0 to turn off deblending. status = sep_extract(&im, 2 * bkg->globalrms, SEP_THRESH_ABS, 10, conv, 3, 3, SEP_FILTER_CONV, 32, 1.0, 1, 1.0, &catalog); if (status != 0) goto exit; -#if 0 - // #4 Aperture photometry - im.noise = &(bkg->globalrms); /* set image noise level */ - im.ndtype = SEP_TFLOAT; - fluxt = flux = (double *)malloc(catalog->nobj * sizeof(double)); - fluxerrt = fluxerr = (double *)malloc(catalog->nobj * sizeof(double)); - areat = area = (double *)malloc(catalog->nobj * sizeof(double)); - flagt = flag = (short *)malloc(catalog->nobj * sizeof(short)); - for (int i = 0; i < catalog->nobj; i++, fluxt++, fluxerrt++, flagt++, areat++) - sep_sum_circle(&im, catalog->x[i], catalog->y[i], 10.0, 5, 0, fluxt, fluxerrt, areat, flagt); -#endif - // TODO // Must detect edge detection // Must limit to brightest 100 (by flux) centers // Should probably use ellipse to draw instead of simple circle? // Useful for galaxies and also elenogated stars. for (int i = 0; i < catalog->nobj; i++) { double flux = catalog->flux[i]; // Get HFR sep_flux_radius(&im, catalog->x[i], catalog->y[i], maxRadius, 5, 0, &flux, requested_frac, 2, flux_fractions, &flux_flag); auto * center = new Edge(); center->x = catalog->x[i] + x + 0.5; center->y = catalog->y[i] + y + 0.5; center->val = catalog->peak[i]; center->sum = flux; center->HFR = center->width = flux_fractions[0]; if (flux_fractions[1] < maxRadius) center->width = flux_fractions[1] * 2; edges.append(center); } // Let's sort edges, starting with widest std::sort(edges.begin(), edges.end(), [](const Edge * edge1, const Edge * edge2) -> bool { return edge1->width > edge2->width;}); // Take only the first 100 stars { int starCount = qMin(100, edges.count()); for (int i = 0; i < starCount; i++) starCenters.append(edges[i]); } edges.clear(); qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << "#" << "#X" << "#Y" << "#Flux" << "#Width" << "#HFR"; for (int i = 0; i < starCenters.count(); i++) qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << i << starCenters[i]->x << starCenters[i]->y << starCenters[i]->sum << starCenters[i]->width << starCenters[i]->HFR; exit: if (stats.bitpix != FLOAT_IMG) delete [] data; sep_bkg_free(bkg); sep_catalog_free(catalog); free(imback); free(flux); free(fluxerr); free(area); free(flag); if (status != 0) { char errorMessage[512]; sep_get_errmsg(status, errorMessage); qCritical(KSTARS_FITS) << errorMessage; return -1; } return starCenters.count(); } template void FITSData::getFloatBuffer(float * buffer, int x, int y, int w, int h) { auto * rawBuffer = reinterpret_cast(m_ImageBuffer); float * floatPtr = buffer; int x2 = x + w; int y2 = y + h; for (int y1 = y; y1 < y2; y1++) { int offset = y1 * stats.width; for (int x1 = x; x1 < x2; x1++) { *floatPtr++ = rawBuffer[offset + x1]; } } } void FITSData::saveStatistics(Statistic &other) { other = stats; } void FITSData::restoreStatistics(Statistic &other) { stats = other; } diff --git a/kstars/fitsviewer/fitsdata.h b/kstars/fitsviewer/fitsdata.h index b9456da6c..5481e0e0b 100644 --- a/kstars/fitsviewer/fitsdata.h +++ b/kstars/fitsviewer/fitsdata.h @@ -1,562 +1,556 @@ /*************************************************************************** fitsimage.cpp - FITS Image ------------------- begin : Tue Feb 24 2004 copyright : (C) 2004 by Jasem Mutlaq email : mutlaqja@ikarustech.com ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* * See http://members.aol.com/pkirchg for more details. * ***************************************************************************/ #pragma once #include "config-kstars.h" #include "bayer.h" #include "fitscommon.h" #ifdef WIN32 // This header must be included before fitsio.h to avoid compiler errors with Visual Studio #include #endif #include #include #include #include #include #ifndef KSTARS_LITE #include #ifdef HAVE_WCSLIB #include #endif #endif #define MINIMUM_PIXEL_RANGE 5 #define MINIMUM_STDVAR 5 class QProgressDialog; class SkyObject; class SkyPoint; class FITSHistogram; typedef struct { float ra; float dec; } wcs_point; class Edge { public: float x; float y; int val; int scanned; float width; float HFR; float sum; }; class FITSSkyObject : public QObject { Q_OBJECT public: explicit FITSSkyObject(SkyObject *object, int xPos, int yPos); SkyObject *skyObject(); int x(); int y(); void setX(int xPos); void setY(int yPos); private: SkyObject *skyObjectStored; int xLoc; int yLoc; }; class FITSData : public QObject { Q_OBJECT // Name of FITS file Q_PROPERTY(QString filename READ filename) // Size of file in bytes Q_PROPERTY(qint64 size READ size) // Width in pixels Q_PROPERTY(quint16 width READ width) // Height in pixels Q_PROPERTY(quint16 height READ height) // FITS MODE --> Normal, Focus, Guide..etc Q_PROPERTY(FITSMode mode MEMBER m_Mode) // 1 channel (grayscale) or 3 channels (RGB) Q_PROPERTY(quint8 channels READ channels) // Data type (BYTE, SHORT, INT..etc) Q_PROPERTY(quint32 dataType MEMBER m_DataType) // Bits per pixel Q_PROPERTY(quint8 bpp READ bpp WRITE setBPP) // Does FITS have WSC header? Q_PROPERTY(bool hasWCS READ hasWCS) // Does FITS have bayer data? Q_PROPERTY(bool hasDebyaer READ hasDebayer) public: explicit FITSData(FITSMode fitsMode = FITS_NORMAL); explicit FITSData(const FITSData *other); ~FITSData(); /** Structure to hold FITS Header records */ typedef struct { QString key; /** FITS Header Key */ QVariant value; /** FITS Header Value */ QString comment; /** FITS Header Comment, if any */ } Record; /// Stats struct to hold statisical data about the FITS data typedef struct { double min[3] = {0}, max[3] = {0}; double mean[3] = {0}; double stddev[3] = {0}; double median[3] = {0}; double SNR { 0 }; int bitpix { 8 }; int bytesPerPixel { 1 }; int ndim { 2 }; int64_t size { 0 }; uint32_t samples_per_channel { 0 }; uint16_t width { 0 }; uint16_t height { 0 }; } Statistic; /** * @brief loadFITS Loading FITS file asynchronously. * @param inFilename Path to FITS file (or compressed fits.gz) * @param silent If set, error messages are ignored. If set to false, the error message will get displayed in a popup. * @return A QFuture that can be watched until the async operation is complete. */ QFuture loadFITS(const QString &inFilename, bool silent = true); /* Save FITS */ int saveFITS(const QString &newFilename); /* Rescale image lineary from image_buffer, fit to window if desired */ int rescale(FITSZoom type); /* Calculate stats */ void calculateStats(bool refresh = false); /* Check if a particular point exists within the image */ bool contains(const QPointF &point) const; // Access functions void clearImageBuffers(); void setImageBuffer(uint8_t *buffer); uint8_t *getImageBuffer(); // Statistics void saveStatistics(Statistic &other); void restoreStatistics(Statistic &other); uint16_t width() const { return stats.width; } uint16_t height() const { return stats.height; } int64_t size() const { return stats.size; } int channels() const { return m_Channels; } double getMin(uint8_t channel = 0) const { return stats.min[channel]; } double getMax(uint8_t channel = 0) const { return stats.max[channel]; } void setMinMax(double newMin, double newMax, uint8_t channel = 0); void getMinMax(double *min, double *max, uint8_t channel = 0) const { *min = stats.min[channel]; *max = stats.max[channel]; } void setStdDev(double value, uint8_t channel = 0) { stats.stddev[channel] = value; } double getStdDev(uint8_t channel = 0) const { return stats.stddev[channel]; } void setMean(double value, uint8_t channel = 0) { stats.mean[channel] = value; } double getMean(uint8_t channel = 0) const { return stats.mean[channel]; } void setMedian(double val, uint8_t channel = 0) { stats.median[channel] = val; } double getMedian(uint8_t channel = 0) const { return stats.median[channel]; } int getBytesPerPixel() const { return stats.bytesPerPixel; } void setSNR(double val) { stats.SNR = val; } double getSNR() const { return stats.SNR; } void setBPP(uint8_t value) { stats.bitpix = value; } uint32_t bpp() const { return stats.bitpix; } double getADU() const; // FITS Record bool getRecordValue(const QString &key, QVariant &value) const; const QList &getRecords() const { return records; } // Star Detection - Native KStars implementation void setStarAlgorithm(StarAlgorithm algorithm) { starAlgorithm = algorithm; } int getDetectedStars() const { return starCenters.count(); } bool areStarsSearched() const { return starsSearched; } void appendStar(Edge *newCenter) { starCenters.append(newCenter); } QList getStarCenters() const { return starCenters; } QList getStarCentersInSubFrame(QRect subFrame) const; int findStars(StarAlgorithm algorithm = ALGORITHM_CENTROID, const QRect &trackingBox = QRect()); void getCenterSelection(int *x, int *y); int findOneStar(const QRect &boundary); // Star Detection - Partially customized Canny edge detection algorithm static int findCannyStar(FITSData *data, const QRect &boundary = QRect()); template static int findCannyStar(FITSData *data, const QRect &boundary); // Use SEP (Sextractor Library) to find stars template void getFloatBuffer(float *buffer, int x, int y, int w, int h); int findSEPStars(const QRect &boundary = QRect()); // Apply ring filter to searched stars int filterStars(const float innerRadius, const float outerRadius); // Half Flux Radius Edge *getMaxHFRStar() const { return maxHFRStar; } double getHFR(HFRType type = HFR_AVERAGE); double getHFR(int x, int y); // WCS // Check if image has valid WCS header information and set HasWCS accordingly. Call in loadFITS() bool checkForWCS(); // Does image have valid WCS? bool hasWCS() { return HasWCS; } // Load WCS data bool loadWCS(); // Is WCS Image loaded? bool isWCSLoaded() { return WCSLoaded; } wcs_point *getWCSCoord() { return wcs_coord; } /** * @brief wcsToPixel Given J2000 (RA0,DE0) coordinates. Find in the image the corresponding pixel coordinates. * @param wcsCoord Coordinates of target * @param wcsPixelPoint Return XY FITS coordinates * @param wcsImagePoint Return XY Image coordinates * @return True if conversion is successful, false otherwise. */ bool wcsToPixel(SkyPoint &wcsCoord, QPointF &wcsPixelPoint, QPointF &wcsImagePoint); /** * @brief pixelToWCS Convert Pixel coordinates to J2000 world coordinates * @param wcsPixelPoint Pixel coordinates in XY Image space. * @param wcsCoord Store back WCS world coordinate in wcsCoord * @return True if successful, false otherwise. */ bool pixelToWCS(const QPointF &wcsPixelPoint, SkyPoint &wcsCoord); /** * @brief createWCSFile Create a new FITS file given the WCS information supplied. Construct the necessary WCS keywords and save the * new file as the current active file * @param newWCSFile New file name * @param orientation Solver orientation, degrees E of N. * @param ra J2000 Right Ascension * @param dec J2000 Declination * @param pixscale Pixel scale in arcsecs per pixel * @return True if file is successfully created and saved, false otherwise. */ bool createWCSFile(const QString &newWCSFile, double orientation, double ra, double dec, double pixscale); // Debayer bool hasDebayer() { return HasDebayer; } bool debayer(); bool debayer_8bit(); bool debayer_16bit(); void getBayerParams(BayerParams *param); void setBayerParams(BayerParams *param); // Histogram #ifndef KSTARS_LITE void setHistogram(FITSHistogram *inHistogram) { histogram = inHistogram; } #endif // Filter void applyFilter(FITSScale type, uint8_t *image = nullptr, QVector *targetMin = nullptr, QVector *targetMax = nullptr); // Rotation counter. We keep count to rotate WCS keywords on save int getRotCounter() const; void setRotCounter(int value); // Filename const QString &filename() const { return m_Filename; } const QString &compressedFilename() const { return m_compressedFilename; } bool isTempFile() const { return m_isTemporary; } bool isCompressed() const { return m_isCompressed; } // Horizontal flip counter. We keep count to rotate WCS keywords on save int getFlipHCounter() const; void setFlipHCounter(int value); // Horizontal flip counter. We keep count to rotate WCS keywords on save int getFlipVCounter() const; void setFlipVCounter(int value); #ifndef KSTARS_LITE #ifdef HAVE_WCSLIB void findObjectsInImage(double world[], double phi, double theta, double imgcrd[], double pixcrd[], int stat[]); #endif #endif QList getSkyObjects(); QList objList; //Does this need to be public?? // Create autostretch image from FITS File static QImage FITSToImage(const QString &filename); /** * @brief ImageToFITS Convert an image file with supported format to a FITS file. * @param filename full path to filename without extension * @param format file extension. Supported formats are whatever supported by Qt (e.g. PNG, JPG,..etc) * @param output Output temporary file path. The created file is generated by the function and store in output. * @return True if conversion is successful, false otherwise. */ static bool ImageToFITS(const QString &filename, const QString &format, QString &output); bool getAutoRemoveTemporaryFITS() const; void setAutoRemoveTemporaryFITS(bool value); QString getLastError() const; signals: void converted(QImage); private: bool privateLoad(bool silent); void rotWCSFITS(int angle, int mirror); bool checkCollision(Edge *s1, Edge *s2); int calculateMinMax(bool refresh = false); bool checkDebayer(); void readWCSKeys(); // FITS Record bool parseHeader(); //int getFITSRecord(QString &recordList, int &nkeys); // Templated functions template bool debayer(); template bool rotFITS(int rotate, int mirror); // Apply Filter template void applyFilter(FITSScale type, uint8_t *targetImage, QVector * min = nullptr, QVector * max = nullptr); // Star Detect - Centroid template int findCentroid(const QRect &boundary, int initStdDev, int minEdgeWidth); int findCentroid(const QRect &boundary = QRect(), int initStdDev = MINIMUM_STDVAR, int minEdgeWidth = MINIMUM_PIXEL_RANGE); // Star Detect - Threshold template int findOneStar(const QRect &boundary); template void calculateMinMax(); template QPair getParitionMinMax(uint32_t start, uint32_t stride); /* Calculate running average & standard deviation using Welford’s method for computing variance */ template void runningAverageStdDev(); template QPair getSquaredSumAndMean(uint32_t start, uint32_t stride); // Sobel detector by Gonzalo Exequiel Pedone template void sobel(QVector &gradient, QVector &direction); template void convertToQImage(double dataMin, double dataMax, double scale, double zero, QImage &image); // Give unique IDs to each contiguous region int partition(int width, int height, QVector &gradient, QVector &ids); void trace(int width, int height, int id, QVector &image, QVector &ids, int x, int y); -#if 0 - QVector thinning(int width, int height, const QVector &gradient, const QVector &direction); - QVector threshold(int thLow, int thHi, const QVector &image); - QVector hysteresis(int width, int height, const QVector &image); -#endif - #ifndef KSTARS_LITE FITSHistogram *histogram { nullptr }; // Pointer to the FITS data histogram #endif /// Pointer to CFITSIO FITS file struct fitsfile *fptr { nullptr }; /// FITS image data type (TBYTE, TUSHORT, TINT, TFLOAT, TLONG, TDOUBLE) uint32_t m_DataType { 0 }; /// Number of channels uint8_t m_Channels { 1 }; /// Generic data image buffer uint8_t *m_ImageBuffer { nullptr }; /// Is this a temporary file or one loaded from disk? bool m_isTemporary { false }; /// is this file compress (.fits.fz)? bool m_isCompressed { false }; /// Did we search for stars yet? bool starsSearched { false }; ///Star Selection Algorithm StarAlgorithm starAlgorithm { ALGORITHM_GRADIENT }; /// Do we have WCS keywords in this FITS data? bool HasWCS { false }; /// Is the image debayarable? bool HasDebayer { false }; /// Is WCS data loaded? bool WCSLoaded { false }; /// Do we need to mark stars for the user? bool markStars { false }; /// Our very own file name QString m_Filename, m_compressedFilename; /// FITS Mode (Normal, WCS, Guide, Focus..etc) FITSMode m_Mode; /// How many times the image was rotated? Useful for WCS keywords rotation on save. int rotCounter { 0 }; /// How many times the image was flipped horizontally? int flipHCounter { 0 }; /// How many times the image was flipped vertically? int flipVCounter { 0 }; /// Pointer to WCS coordinate data, if any. wcs_point *wcs_coord { nullptr }; /// WCS Struct struct wcsprm *wcs { nullptr }; /// All the stars we detected, if any. QList starCenters; QList localStarCenters; /// The biggest fattest star in the image. Edge *maxHFRStar { nullptr }; uint8_t *bayerBuffer { nullptr }; /// Bayer parameters BayerParams debayerParams; Statistic stats; // A list of header records QList records; /// Remove temporary files after closing bool autoRemoveTemporaryFITS { true }; QString lastError; static const QString m_TemporaryPath; }; diff --git a/kstars/fitsviewer/fitsview.cpp b/kstars/fitsviewer/fitsview.cpp index 3528f2fea..753bedfcb 100644 --- a/kstars/fitsviewer/fitsview.cpp +++ b/kstars/fitsviewer/fitsview.cpp @@ -1,1965 +1,1743 @@ /* FITS View Copyright (C) 2003-2017 Jasem Mutlaq Copyright (C) 2016-2017 Robert Lancaster This application is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. */ #include "config-kstars.h" #include "fitsview.h" #include "fitsdata.h" #include "fitslabel.h" #include "kspopupmenu.h" #include "kstarsdata.h" #include "ksutils.h" #include "Options.h" #include "skymap.h" #include "fits_debug.h" #ifdef HAVE_INDI #include "basedevice.h" #include "indi/indilistener.h" #endif #include #include #include #include #include #include #include #define BASE_OFFSET 50 #define ZOOM_DEFAULT 100.0 #define ZOOM_MIN 10 #define ZOOM_MAX 400 #define ZOOM_LOW_INCR 10 #define ZOOM_HIGH_INCR 50 FITSView::FITSView(QWidget * parent, FITSMode fitsMode, FITSScale filterType) : QScrollArea(parent), zoomFactor(1.2) { grabGesture(Qt::PinchGesture); image_frame.reset(new FITSLabel(this)); filter = filterType; mode = fitsMode; setBackgroundRole(QPalette::Dark); markerCrosshair.setX(0); markerCrosshair.setY(0); setBaseSize(740, 530); connect(image_frame.get(), SIGNAL(newStatus(QString, FITSBar)), this, SIGNAL(newStatus(QString, FITSBar))); connect(image_frame.get(), SIGNAL(pointSelected(int, int)), this, SLOT(processPointSelection(int, int))); connect(image_frame.get(), SIGNAL(markerSelected(int, int)), this, SLOT(processMarkerSelection(int, int))); connect(&wcsWatcher, SIGNAL(finished()), this, SLOT(syncWCSState())); connect(&fitsWatcher, &QFutureWatcher::finished, this, &FITSView::loadInFrame); image_frame->setMouseTracking(true); setCursorMode( selectCursor); //This is the default mode because the Focus and Align FitsViews should not be in dragMouse mode noImageLabel = new QLabel(); noImage.load(":/images/noimage.png"); noImageLabel->setPixmap(noImage); noImageLabel->setAlignment(Qt::AlignCenter); this->setWidget(noImageLabel); redScopePixmap = QPixmap(":/icons/center_telescope_red.svg").scaled(32, 32, Qt::KeepAspectRatio, Qt::FastTransformation); magentaScopePixmap = QPixmap(":/icons/center_telescope_magenta.svg").scaled(32, 32, Qt::KeepAspectRatio, Qt::FastTransformation); - - //if (fitsMode == FITS_GUIDE) - //connect(image_frame.get(), SIGNAL(pointSelected(int,int)), this, SLOT(processPointSelection(int,int))); - - // Default size - //resize(INITIAL_W, INITIAL_H); } FITSView::~FITSView() { fitsWatcher.waitForFinished(); wcsWatcher.waitForFinished(); delete (imageData); } /** This method looks at what mouse mode is currently selected and updates the cursor to match. */ void FITSView::updateMouseCursor() { if (cursorMode == dragCursor) { if (horizontalScrollBar()->maximum() > 0 || verticalScrollBar()->maximum() > 0) { if (!image_frame->getMouseButtonDown()) viewport()->setCursor(Qt::PointingHandCursor); else viewport()->setCursor(Qt::ClosedHandCursor); } else viewport()->setCursor(Qt::CrossCursor); } else if (cursorMode == selectCursor) { viewport()->setCursor(Qt::CrossCursor); } else if (cursorMode == scopeCursor) { viewport()->setCursor(QCursor(redScopePixmap, 10, 10)); } else if (cursorMode == crosshairCursor) { viewport()->setCursor(QCursor(magentaScopePixmap, 10, 10)); } } /** This is how the mouse mode gets set. The default for a FITSView in a FITSViewer should be the dragMouse The default for a FITSView in the Focus or Align module should be the selectMouse The different defaults are accomplished by putting making the actual default mouseMode the selectMouse, but when a FITSViewer loads an image, it immediately makes it the dragMouse. */ void FITSView::setCursorMode(CursorMode mode) { cursorMode = mode; updateMouseCursor(); if (mode == scopeCursor && imageHasWCS()) { if (!imageData->isWCSLoaded() && !wcsWatcher.isRunning()) { QFuture future = QtConcurrent::run(imageData, &FITSData::loadWCS); wcsWatcher.setFuture(future); } } } void FITSView::resizeEvent(QResizeEvent * event) { if ((imageData == nullptr) && noImageLabel != nullptr) { noImageLabel->setPixmap( noImage.scaled(width() - 20, height() - 20, Qt::KeepAspectRatio, Qt::FastTransformation)); noImageLabel->setFixedSize(width() - 5, height() - 5); } QScrollArea::resizeEvent(event); } -#if 0 - -bool FITSView::loadFITS(const QString &inFilename, bool silent) -{ - if (floatingToolBar != nullptr) - { - floatingToolBar->setVisible(true); - } - - QProgressDialog fitsProg(this); - - bool setBayerParams = false; - - BayerParams param; - if ((imageData != nullptr) && imageData->hasDebayer()) - { - setBayerParams = true; - imageData->getBayerParams(¶m); - } - - // In case loadWCS is still running for previous image data, let's wait until it's over - wcsWatcher.waitForFinished(); - - delete imageData; - imageData = nullptr; - - filterStack.clear(); - filterStack.push(FITS_NONE); - if (filter != FITS_NONE) - filterStack.push(filter); - - imageData = new FITSData(mode); - - if (setBayerParams) - imageData->setBayerParams(¶m); - - if (mode == FITS_NORMAL) - { - fitsProg.setWindowModality(Qt::WindowModal); - fitsProg.setLabelText(i18n("Please hold while loading FITS file...")); - fitsProg.setWindowTitle(i18n("Loading FITS")); - fitsProg.setValue(10); - qApp->processEvents(); - } - - if (!imageData->loadFITS(inFilename, silent)) - return false; - - if (mode == FITS_NORMAL) - { - if (fitsProg.wasCanceled()) - return false; - else - { - fitsProg.setValue(65); - qApp->processEvents(); - } - } - - emit debayerToggled(imageData->hasDebayer()); - - currentWidth = imageData->width(); - currentHeight = imageData->height(); - - image_width = currentWidth; - image_height = currentHeight; - - image_frame->setSize(image_width, image_height); - - initDisplayImage(); - - // Rescale to fits window - if (firstLoad) - { - currentZoom = 100; - - if (rescale(ZOOM_FIT_WINDOW) != 0) - return false; - - firstLoad = false; - } - else - { - if (rescale(ZOOM_KEEP_LEVEL) != 0) - return false; - } - - if (mode == FITS_NORMAL) - { - if (fitsProg.wasCanceled()) - return false; - else - { - fitsProg.setValue(100); - qApp->processEvents(); - } - } - - setAlignment(Qt::AlignCenter); - - // Load WCS data now if selected and image contains valid WCS header - if (imageData->hasWCS() && Options::autoWCS() && (mode == FITS_NORMAL || mode == FITS_ALIGN) && !wcsWatcher.isRunning()) - { - QFuture future = QtConcurrent::run(imageData, &FITSData::loadWCS); - wcsWatcher.setFuture(future); - } - else - syncWCSState(); - - if (isVisible()) - emit newStatus(QString("%1x%2").arg(image_width).arg(image_height), FITS_RESOLUTION); - - if (showStarProfile) - { - if(floatingToolBar != nullptr) - toggleProfileAction->setChecked(true); - QTimer::singleShot(100, this, SLOT(viewStarProfile())); //Need to wait till the Focus module finds stars, if its the Focus module. - } - - updateFrame(); - - emit imageLoaded(); - - return true; -} -#endif void FITSView::loadFITS(const QString &inFilename, bool silent) { if (floatingToolBar != nullptr) { floatingToolBar->setVisible(true); } bool setBayerParams = false; BayerParams param; if ((imageData != nullptr) && imageData->hasDebayer()) { setBayerParams = true; imageData->getBayerParams(¶m); } // In case image is still loading, wait until it is done. fitsWatcher.waitForFinished(); // In case loadWCS is still running for previous image data, let's wait until it's over wcsWatcher.waitForFinished(); delete imageData; imageData = nullptr; filterStack.clear(); filterStack.push(FITS_NONE); if (filter != FITS_NONE) filterStack.push(filter); imageData = new FITSData(mode); if (setBayerParams) imageData->setBayerParams(¶m); fitsWatcher.setFuture(imageData->loadFITS(inFilename, silent)); } void FITSView::loadInFrame() { // Check if the loading was OK if (fitsWatcher.result() == false) { m_LastError = imageData->getLastError(); emit failed(); return; } // Notify if there is debayer data. emit debayerToggled(imageData->hasDebayer()); // Set current width and height currentWidth = imageData->width(); currentHeight = imageData->height(); image_width = currentWidth; image_height = currentHeight; image_frame->setSize(image_width, image_height); // Init the display image initDisplayImage(); uint8_t * ASImageBuffer = nullptr; - // if (Options::autoStretch() && (filter == FITS_NONE || (filter >= FITS_ROTATE_CW && filter <= FITS_FLIP_V))) - // { - // // If we perform autostretch, we need to create a buffer to save the raw image data before - // // autostretch filter operation changes the data. - // // After rescaling is done, we - // uint32_t totalBytes = image_width * image_height *imageData->channels() * imageData->getBytesPerPixel(); - // ASImageBuffer = new uint8_t[totalBytes]; - // memcpy(ASImageBuffer, imageData->getImageBuffer(), totalBytes); - // imageData->applyFilter(FITS_AUTO_STRETCH); - // } - // else - // imageData->applyFilter(filter); - imageData->applyFilter(filter); - // if (Options::autoStretch()) - // imageData->applyFilter(FITS_AUTO_STRETCH); // Rescale to fits window on first load if (firstLoad) { currentZoom = 100; if (rescale(ZOOM_FIT_WINDOW) == false) { m_LastError = i18n("Rescaling image failed."); delete [] ASImageBuffer; emit failed(); return; } firstLoad = false; } else { if (rescale(ZOOM_KEEP_LEVEL) == false) { m_LastError = i18n("Rescaling image failed."); delete [] ASImageBuffer; emit failed(); return; } } - // Restore original raw buffer after Autostretch if applicable - // if (ASImageBuffer) - // { - // imageData->setImageBuffer(ASImageBuffer); - // } - setAlignment(Qt::AlignCenter); // Load WCS data now if selected and image contains valid WCS header if (imageData->hasWCS() && Options::autoWCS() && (mode == FITS_NORMAL || mode == FITS_ALIGN) && !wcsWatcher.isRunning()) { QFuture future = QtConcurrent::run(imageData, &FITSData::loadWCS); wcsWatcher.setFuture(future); } else syncWCSState(); if (isVisible()) emit newStatus(QString("%1x%2").arg(image_width).arg(image_height), FITS_RESOLUTION); if (showStarProfile) { if(floatingToolBar != nullptr) toggleProfileAction->setChecked(true); //Need to wait till the Focus module finds stars, if its the Focus module. QTimer::singleShot(100, this, SLOT(viewStarProfile())); } scaledImage = QImage(); updateFrame(); emit loaded(); } int FITSView::saveFITS(const QString &newFilename) { return imageData->saveFITS(newFilename); } bool FITSView::rescale(FITSZoom type) { switch (imageData->property("dataType").toInt()) { case TBYTE: return rescale(type); case TSHORT: return rescale(type); case TUSHORT: return rescale(type); case TLONG: return rescale(type); case TULONG: return rescale(type); case TFLOAT: return rescale(type); case TLONGLONG: return rescale(type); case TDOUBLE: return rescale(type); default: break; } return false; } FITSView::CursorMode FITSView::getCursorMode() { return cursorMode; } void FITSView::enterEvent(QEvent * event) { Q_UNUSED(event) if ((floatingToolBar != nullptr) && (imageData != nullptr)) { QPointer eff = new QGraphicsOpacityEffect(this); floatingToolBar->setGraphicsEffect(eff); QPointer a = new QPropertyAnimation(eff, "opacity"); a->setDuration(500); a->setStartValue(0.2); a->setEndValue(1); a->setEasingCurve(QEasingCurve::InBack); a->start(QPropertyAnimation::DeleteWhenStopped); } } void FITSView::leaveEvent(QEvent * event) { Q_UNUSED(event) if ((floatingToolBar != nullptr) && (imageData != nullptr)) { QPointer eff = new QGraphicsOpacityEffect(this); floatingToolBar->setGraphicsEffect(eff); QPointer a = new QPropertyAnimation(eff, "opacity"); a->setDuration(500); a->setStartValue(1); a->setEndValue(0.2); a->setEasingCurve(QEasingCurve::OutBack); a->start(QPropertyAnimation::DeleteWhenStopped); } } template bool FITSView::rescale(FITSZoom type) { if (rawImage.isNull()) return false; uint8_t * imageBuffer = imageData->getImageBuffer(); uint8_t * displayBuffer = nullptr; uint32_t size = imageData->width() * imageData->height(); QVector min(3), max(3); if (Options::autoStretch()) { displayBuffer = new uint8_t[size * imageData->channels() * imageData->getBytesPerPixel()]; memcpy(displayBuffer, imageBuffer, size * imageData->channels() * imageData->getBytesPerPixel()); imageData->applyFilter(FITS_AUTO_STRETCH, displayBuffer, &min, &max); } else { displayBuffer = imageBuffer; for (int i = 0; i < 3; i++) { min[i] = imageData->getMin(i); max[i] = imageData->getMax(i); } } -#if 0 - int BBP = imageData->getBytesPerPixel(); - filter = filterStack.last(); - - if (Options::autoStretch() && (filter == FITS_NONE || (filter >= FITS_ROTATE_CW && filter <= FITS_FLIP_V))) - { - image_buffer = new uint8_t[size * imageData->channels() * BBP]; - memcpy(image_buffer, imageData->getImageBuffer(), size * imageData->channels() * BBP); - - displayBuffer = true; - - double data_min = -1; - double data_max = -1; - - imageData->applyFilter(FITS_AUTO_STRETCH, image_buffer, &data_min, &data_max); - - min = data_min; - max = data_max; - } - else - { - imageData->applyFilter(filter); - imageData->getMinMax(&min, &max); - } -#endif - scaledImage = QImage(); auto * buffer = reinterpret_cast(displayBuffer); if (min[0] == max[0]) { rawImage.fill(Qt::white); emit newStatus(i18n("Image is saturated."), FITS_MESSAGE); } else { if (image_height != imageData->height() || image_width != imageData->width()) { image_width = imageData->width(); image_height = imageData->height(); initDisplayImage(); if (isVisible()) emit newStatus(QString("%1x%2").arg(image_width).arg(image_height), FITS_RESOLUTION); } image_frame->setScaledContents(true); currentWidth = rawImage.width(); currentHeight = rawImage.height(); if (imageData->channels() == 1) { double range = max[0] - min[0]; double bscale = 255. / range; double bzero = (-min[0]) * (255. / range); QVector> futures; /* Fill in pixel values using indexed map, linear scale */ for (uint32_t j = 0; j < image_height; j++) { futures.append(QtConcurrent::run([ = ]() { T * runningBuffer = buffer + j * image_width; uint8_t * scanLine = rawImage.scanLine(j); for (uint32_t i = 0; i < image_width; i++) { //scanLine[i] = qBound(0, static_cast(runningBuffer[i] * bscale + bzero), 255); scanLine[i] = qBound(0.0, runningBuffer[i] * bscale + bzero, 255.0); } })); } for(QFuture future : futures) future.waitForFinished(); } else { QVector> futures; double bscale[3], bzero[3]; for (int i = 0; i < 3; i++) { bscale[i] = 255. / (max[i] - min[i]); bzero[i] = (-min[i]) * (255. / (max[i] - min[i])); } /* Fill in pixel values using indexed map, linear scale */ for (uint32_t j = 0; j < image_height; j++) { futures.append(QtConcurrent::run([ = ]() { auto * scanLine = reinterpret_cast((rawImage.scanLine(j))); T * runningBufferR = buffer + j * image_width; T * runningBufferG = buffer + j * image_width + size; T * runningBufferB = buffer + j * image_width + size * 2; for (uint32_t i = 0; i < image_width; i++) { scanLine[i] = qRgb(runningBufferR[i] * bscale[0] + bzero[0], runningBufferG[i] * bscale[1] + bzero[1], runningBufferB[i] * bscale[2] + bzero[2]); } })); } for(QFuture future : futures) future.waitForFinished(); } -#if 0 - if (imageData->getNumOfChannels() == 1) - { - /* Fill in pixel values using indexed map, linear scale */ - for (int j = 0; j < image_height; j++) - { - unsigned char * scanLine = display_image->scanLine(j); - - for (int i = 0; i < image_width; i++) - { - val = buffer[j * image_width + i] * bscale + bzero; - scanLine[i] = qBound(0.0, val, 255.0); - } - } - } - else - { - double rval = 0, gval = 0, bval = 0; - QRgb value; - /* Fill in pixel values using indexed map, linear scale */ - for (int j = 0; j < image_height; j++) - { - QRgb * scanLine = reinterpret_cast((display_image->scanLine(j))); - - for (int i = 0; i < image_width; i++) - { - rval = buffer[j * image_width + i]; - gval = buffer[j * image_width + i + size]; - bval = buffer[j * image_width + i + size * 2]; - - value = qRgb(rval * bscale + bzero, gval * bscale + bzero, bval * bscale + bzero); - - scanLine[i] = value; - } - } - } -#endif } // Clear memory if it was allocated. if (displayBuffer != imageBuffer) delete [] displayBuffer; switch (type) { case ZOOM_FIT_WINDOW: if ((rawImage.width() > width() || rawImage.height() > height())) { double w = baseSize().width() - BASE_OFFSET; double h = baseSize().height() - BASE_OFFSET; if (!firstLoad) { w = viewport()->rect().width() - BASE_OFFSET; h = viewport()->rect().height() - BASE_OFFSET; } // Find the zoom level which will enclose the current FITS in the current window size double zoomX = floor((w / static_cast(currentWidth)) * 100.); double zoomY = floor((h / static_cast(currentHeight)) * 100.); (zoomX < zoomY) ? currentZoom = zoomX : currentZoom = zoomY; currentWidth = image_width * (currentZoom / ZOOM_DEFAULT); currentHeight = image_height * (currentZoom / ZOOM_DEFAULT); if (currentZoom <= ZOOM_MIN) emit actionUpdated("view_zoom_out", false); } else { currentZoom = 100; currentWidth = image_width; currentHeight = image_height; } break; case ZOOM_KEEP_LEVEL: { currentWidth = image_width * (currentZoom / ZOOM_DEFAULT); currentHeight = image_height * (currentZoom / ZOOM_DEFAULT); } break; default: currentZoom = 100; break; } setWidget(image_frame.get()); if (type != ZOOM_KEEP_LEVEL) emit newStatus(QString("%1%").arg(currentZoom), FITS_ZOOM); return true; } void FITSView::ZoomIn() { if (currentZoom >= ZOOM_DEFAULT && Options::limitedResourcesMode()) { emit newStatus(i18n("Cannot zoom in further due to active limited resources mode."), FITS_MESSAGE); return; } if (currentZoom < ZOOM_DEFAULT) currentZoom += ZOOM_LOW_INCR; else currentZoom += ZOOM_HIGH_INCR; emit actionUpdated("view_zoom_out", true); if (currentZoom >= ZOOM_MAX) { currentZoom = ZOOM_MAX; emit actionUpdated("view_zoom_in", false); } currentWidth = image_width * (currentZoom / ZOOM_DEFAULT); currentHeight = image_height * (currentZoom / ZOOM_DEFAULT); updateFrame(); emit newStatus(QString("%1%").arg(currentZoom), FITS_ZOOM); } void FITSView::ZoomOut() { if (currentZoom <= ZOOM_DEFAULT) currentZoom -= ZOOM_LOW_INCR; else currentZoom -= ZOOM_HIGH_INCR; if (currentZoom <= ZOOM_MIN) { currentZoom = ZOOM_MIN; emit actionUpdated("view_zoom_out", false); } emit actionUpdated("view_zoom_in", true); currentWidth = image_width * (currentZoom / ZOOM_DEFAULT); currentHeight = image_height * (currentZoom / ZOOM_DEFAULT); updateFrame(); emit newStatus(QString("%1%").arg(currentZoom), FITS_ZOOM); } void FITSView::ZoomToFit() { if (rawImage.isNull() == false) { rescale(ZOOM_FIT_WINDOW); updateFrame(); } } void FITSView::setStarFilterRange(float const innerRadius, float const outerRadius) { starFilter.innerRadius = innerRadius; starFilter.outerRadius = outerRadius; } int FITSView::filterStars() { return starFilter.used() ? imageData->filterStars(starFilter.innerRadius, starFilter.outerRadius) : imageData->getStarCenters().count(); } void FITSView::updateFrame() { bool ok = false; if (currentZoom != ZOOM_DEFAULT) { // Only scale when necessary if (scaledImage.isNull() || currentWidth != lastWidth || currentHeight != lastHeight) { scaledImage = rawImage.scaled(currentWidth, currentHeight, Qt::KeepAspectRatio, Qt::SmoothTransformation); lastWidth = currentWidth; lastHeight = currentHeight; } ok = displayPixmap.convertFromImage(scaledImage); } else ok = displayPixmap.convertFromImage(rawImage); if (!ok) return; QPainter painter(&displayPixmap); drawOverlay(&painter); if (starFilter.used()) { double const diagonal = std::sqrt(currentWidth * currentWidth + currentHeight * currentHeight) / 2; int const innerRadius = std::lround(diagonal * starFilter.innerRadius); int const outerRadius = std::lround(diagonal * starFilter.outerRadius); QPoint const center(currentWidth / 2, currentHeight / 2); painter.save(); painter.setPen(QPen(Qt::blue, 1, Qt::DashLine)); painter.setOpacity(0.7); painter.setBrush(QBrush(Qt::transparent)); painter.drawEllipse(center, outerRadius, outerRadius); painter.setBrush(QBrush(Qt::blue, Qt::FDiagPattern)); painter.drawEllipse(center, innerRadius, innerRadius); painter.restore(); } image_frame->setPixmap(displayPixmap); image_frame->resize(currentWidth, currentHeight); } void FITSView::ZoomDefault() { if (image_frame != nullptr) { emit actionUpdated("view_zoom_out", true); emit actionUpdated("view_zoom_in", true); currentZoom = ZOOM_DEFAULT; currentWidth = image_width; currentHeight = image_height; updateFrame(); emit newStatus(QString("%1%").arg(currentZoom), FITS_ZOOM); update(); } } void FITSView::drawOverlay(QPainter * painter) { painter->setRenderHint(QPainter::Antialiasing, Options::useAntialias()); if (trackingBoxEnabled && getCursorMode() != FITSView::scopeCursor) drawTrackingBox(painter); if (!markerCrosshair.isNull()) drawMarker(painter); if (showCrosshair) drawCrosshair(painter); if (showObjects) drawObjectNames(painter); if (showEQGrid) drawEQGrid(painter); if (showPixelGrid) drawPixelGrid(painter); if (markStars) drawStarCentroid(painter); } void FITSView::updateMode(FITSMode fmode) { mode = fmode; } void FITSView::drawMarker(QPainter * painter) { painter->setPen(QPen(QColor(KStarsData::Instance()->colorScheme()->colorNamed("TargetColor")), 2)); painter->setBrush(Qt::NoBrush); float pxperdegree = (currentZoom / ZOOM_DEFAULT) * (57.3 / 1.8); float s1 = 0.5 * pxperdegree; float s2 = pxperdegree; float s3 = 2.0 * pxperdegree; float x0 = markerCrosshair.x() * (currentZoom / ZOOM_DEFAULT); float y0 = markerCrosshair.y() * (currentZoom / ZOOM_DEFAULT); float x1 = x0 - 0.5 * s1; float y1 = y0 - 0.5 * s1; float x2 = x0 - 0.5 * s2; float y2 = y0 - 0.5 * s2; float x3 = x0 - 0.5 * s3; float y3 = y0 - 0.5 * s3; //Draw radial lines painter->drawLine(QPointF(x1, y0), QPointF(x3, y0)); painter->drawLine(QPointF(x0 + s2, y0), QPointF(x0 + 0.5 * s1, y0)); painter->drawLine(QPointF(x0, y1), QPointF(x0, y3)); painter->drawLine(QPointF(x0, y0 + 0.5 * s1), QPointF(x0, y0 + s2)); //Draw circles at 0.5 & 1 degrees painter->drawEllipse(QRectF(x1, y1, s1, s1)); painter->drawEllipse(QRectF(x2, y2, s2, s2)); } void FITSView::drawStarCentroid(QPainter * painter) { float const ratio = currentZoom / ZOOM_DEFAULT; if (showStarsHFR) { QFont painterFont; // If we need to print the HFR out, give an arbitrarily sized font to the painter painterFont.setPointSizeF(painterFont.pointSizeF() * 3 * ratio); painter->setFont(painterFont); } painter->setPen(QPen(Qt::red, 2)); QFontMetrics const fontMetrics = painter->fontMetrics(); QRect const boundingRect(0, 0, painter->device()->width(), painter->device()->height()); foreach (auto const &starCenter, imageData->getStarCenters()) { int const x1 = std::round((starCenter->x - starCenter->width / 2.0f) * ratio); int const y1 = std::round((starCenter->y - starCenter->width / 2.0f) * ratio); int const w = std::round(starCenter->width * ratio); // Draw a circle around the detected star painter->drawEllipse(x1, y1, w, w); if (showStarsHFR) { // Ask the painter how large will the HFR text be QString const hfr = QString("%1").arg(starCenter->HFR, 0, 'f', 2); QSize const hfrSize = fontMetrics.size(Qt::TextSingleLine, hfr); // Store the HFR text in a rect QPoint const hfrBottomLeft(x1+w+5, y1+w/2); QRect const hfrRect(hfrBottomLeft.x(), hfrBottomLeft.y() - hfrSize.height(), hfrSize.width(), hfrSize.height()); // Render the HFR text only if it can be displayed entirely if (boundingRect.contains(hfrRect)) { painter->setPen(QPen(Qt::red, 3)); painter->drawText(hfrBottomLeft, hfr); painter->setPen(QPen(Qt::red, 2)); } } } } void FITSView::drawTrackingBox(QPainter * painter) { painter->setPen(QPen(Qt::green, 2)); if (trackingBox.isNull()) return; int x1 = trackingBox.x() * (currentZoom / ZOOM_DEFAULT); int y1 = trackingBox.y() * (currentZoom / ZOOM_DEFAULT); int w = trackingBox.width() * (currentZoom / ZOOM_DEFAULT); int h = trackingBox.height() * (currentZoom / ZOOM_DEFAULT); painter->drawRect(x1, y1, w, h); } /** This Method draws a large Crosshair in the center of the image, it is like a set of axes. */ void FITSView::drawCrosshair(QPainter * painter) { float scale = (currentZoom / ZOOM_DEFAULT); QPointF c = QPointF((qreal)image_width / 2 * scale, (qreal)image_height / 2 * scale); float midX = (float)image_width / 2 * scale; float midY = (float)image_height / 2 * scale; float maxX = (float)image_width * scale; float maxY = (float)image_height * scale; float r = 50 * scale; painter->setPen(QPen(QColor(KStarsData::Instance()->colorScheme()->colorNamed("TargetColor")))); //Horizontal Line to Circle painter->drawLine(0, midY, midX - r, midY); //Horizontal Line past Circle painter->drawLine(midX + r, midY, maxX, midY); //Vertical Line to Circle painter->drawLine(midX, 0, midX, midY - r); //Vertical Line past Circle painter->drawLine(midX, midY + r, midX, maxY); //Circles painter->drawEllipse(c, r, r); painter->drawEllipse(c, r / 2, r / 2); } /** This method is intended to draw a pixel grid onto the image. It first determines useful information from the image. Then it draws the axes on the image if the crosshairs are not displayed. Finally it draws the gridlines so that there will be 4 Gridlines on either side of the axes. Note: This has to start drawing at the center not at the edges because the center axes must be in the center of the image. */ void FITSView::drawPixelGrid(QPainter * painter) { float scale = (currentZoom / ZOOM_DEFAULT); double width = image_width * scale; double height = image_height * scale; double cX = width / 2; double cY = height / 2; double deltaX = width / 10; double deltaY = height / 10; //draw the Axes painter->setPen(QPen(Qt::red)); painter->drawText(cX - 30, height - 5, QString::number((int)((cX) / scale))); painter->drawText(width - 30, cY - 5, QString::number((int)((cY) / scale))); if (!showCrosshair) { painter->drawLine(cX, 0, cX, height); painter->drawLine(0, cY, width, cY); } painter->setPen(QPen(Qt::gray)); //Start one iteration past the Center and draw 4 lines on either side of 0 for (int x = deltaX; x < cX - deltaX; x += deltaX) { painter->drawText(cX + x - 30, height - 5, QString::number((int)((cX + x) / scale))); painter->drawText(cX - x - 30, height - 5, QString::number((int)((cX - x) / scale))); painter->drawLine(cX - x, 0, cX - x, height); painter->drawLine(cX + x, 0, cX + x, height); } //Start one iteration past the Center and draw 4 lines on either side of 0 for (int y = deltaY; y < cY - deltaY; y += deltaY) { painter->drawText(width - 30, cY + y - 5, QString::number((int)((cY + y) / scale))); painter->drawText(width - 30, cY - y - 5, QString::number((int)((cY - y) / scale))); painter->drawLine(0, cY + y, width, cY + y); painter->drawLine(0, cY - y, width, cY - y); } } bool FITSView::imageHasWCS() { if (imageData != nullptr) return imageData->hasWCS(); return false; } void FITSView::drawObjectNames(QPainter * painter) { painter->setPen(QPen(QColor(KStarsData::Instance()->colorScheme()->colorNamed("FITSObjectLabelColor")))); float scale = (currentZoom / ZOOM_DEFAULT); foreach (FITSSkyObject * listObject, imageData->getSkyObjects()) { painter->drawRect(listObject->x() * scale - 5, listObject->y() * scale - 5, 10, 10); painter->drawText(listObject->x() * scale + 10, listObject->y() * scale + 10, listObject->skyObject()->name()); } } /** This method will paint EQ Gridlines in an overlay if there is WCS data present. It determines the minimum and maximum RA and DEC, then it uses that information to judge which gridLines to draw. Then it calls the drawEQGridlines methods below to draw gridlines at those specific RA and Dec values. */ void FITSView::drawEQGrid(QPainter * painter) { float scale = (currentZoom / ZOOM_DEFAULT); if (imageData->hasWCS()) { wcs_point * wcs_coord = imageData->getWCSCoord(); if (wcs_coord != nullptr) { int size = image_width * image_height; double maxRA = -1000; double minRA = 1000; double maxDec = -1000; double minDec = 1000; for (int i = 0; i < (size); i++) { double ra = wcs_coord[i].ra; double dec = wcs_coord[i].dec; if (ra > maxRA) maxRA = ra; if (ra < minRA) minRA = ra; if (dec > maxDec) maxDec = dec; if (dec < minDec) minDec = dec; } auto minDecMinutes = (int)(minDec * 12); //This will force the Dec Scale to 5 arc minutes in the loop auto maxDecMinutes = (int)(maxDec * 12); auto minRAMinutes = (int)(minRA / 15.0 * 120.0); //This will force the scale to 1/2 minutes of RA in the loop from 0 to 50 degrees auto maxRAMinutes = (int)(maxRA / 15.0 * 120.0); double raConvert = 15 / 120.0; //This will undo the calculation above to retrieve the actual RA. double decConvert = 1.0 / 12.0; //This will undo the calculation above to retrieve the actual DEC. if (maxDec > 50 || minDec < -50) { minRAMinutes = (int)(minRA / 15.0 * 60.0); //This will force the scale to 1 min of RA from 50 to 80 degrees maxRAMinutes = (int)(maxRA / 15.0 * 60.0); raConvert = 15 / 60.0; } if (maxDec > 80 || minDec < -80) { minRAMinutes = (int)(minRA / 15.0 * 30); //This will force the scale to 2 min of RA from 80 to 85 degrees maxRAMinutes = (int)(maxRA / 15.0 * 30); raConvert = 15 / 30.0; } if (maxDec > 85 || minDec < -85) { minRAMinutes = (int)(minRA / 15.0 * 6); //This will force the scale to 10 min of RA from 85 to 89 degrees maxRAMinutes = (int)(maxRA / 15.0 * 6); raConvert = 15 / 6.0; } if (maxDec >= 89.25 || minDec <= -89.25) { minRAMinutes = (int)(minRA / 15); //This will force the scale to whole hours of RA in the loop really close to the poles maxRAMinutes = (int)(maxRA / 15); raConvert = 15; } painter->setPen(QPen(Qt::yellow)); QPointF pixelPoint, imagePoint, pPoint; //This section draws the RA Gridlines for (int targetRA = minRAMinutes; targetRA <= maxRAMinutes; targetRA++) { painter->setPen(QPen(Qt::yellow)); double target = targetRA * raConvert; if (eqGridPoints.count() != 0) eqGridPoints.clear(); double increment = std::abs((maxDec - minDec) / 100.0); //This will determine how many points to use to create the RA Line for (double targetDec = minDec; targetDec <= maxDec; targetDec += increment) { SkyPoint pointToGet(target / 15.0, targetDec); bool inImage = imageData->wcsToPixel(pointToGet, pixelPoint, imagePoint); if (inImage) { QPointF pt(pixelPoint.x() * scale, pixelPoint.y() * scale); eqGridPoints.append(pt); } } if (eqGridPoints.count() > 1) { for (int i = 1; i < eqGridPoints.count(); i++) painter->drawLine(eqGridPoints.value(i - 1), eqGridPoints.value(i)); QPointF pt = getPointForGridLabel(); if (pt.x() != -100) { if (maxDec > 50 || maxDec < -50) painter->drawText(pt.x(), pt.y(), QString::number(dms(target).hour()) + "h " + QString::number(dms(target).minute()) + '\''); else painter->drawText(pt.x() - 20, pt.y(), QString::number(dms(target).hour()) + "h " + QString::number(dms(target).minute()) + "' " + QString::number(dms(target).second()) + "''"); } } } //This section draws the DEC Gridlines for (int targetDec = minDecMinutes; targetDec <= maxDecMinutes; targetDec++) { if (eqGridPoints.count() != 0) eqGridPoints.clear(); double increment = std::abs((maxRA - minRA) / 100.0); //This will determine how many points to use to create the Dec Line double target = targetDec * decConvert; for (double targetRA = minRA; targetRA <= maxRA; targetRA += increment) { SkyPoint pointToGet(targetRA / 15, targetDec * decConvert); bool inImage = imageData->wcsToPixel(pointToGet, pixelPoint, imagePoint); if (inImage) { QPointF pt(pixelPoint.x() * scale, pixelPoint.y() * scale); eqGridPoints.append(pt); } } if (eqGridPoints.count() > 1) { for (int i = 1; i < eqGridPoints.count(); i++) painter->drawLine(eqGridPoints.value(i - 1), eqGridPoints.value(i)); QPointF pt = getPointForGridLabel(); if (pt.x() != -100) painter->drawText(pt.x(), pt.y(), QString::number(dms(target).degree()) + "° " + QString::number(dms(target).arcmin()) + '\''); } } //This Section Draws the North Celestial Pole if present SkyPoint NCP(0, 90); bool NCPtest = imageData->wcsToPixel(NCP, pPoint, imagePoint); if (NCPtest) { bool NCPinImage = (pPoint.x() > 0 && pPoint.x() < image_width) && (pPoint.y() > 0 && pPoint.y() < image_height); if (NCPinImage) { painter->fillRect(pPoint.x() * scale - 2, pPoint.y() * scale - 2, 4, 4, KStarsData::Instance()->colorScheme()->colorNamed("TargetColor")); painter->drawText(pPoint.x() * scale + 15, pPoint.y() * scale + 15, i18nc("North Celestial Pole", "NCP")); } } //This Section Draws the South Celestial Pole if present SkyPoint SCP(0, -90); bool SCPtest = imageData->wcsToPixel(SCP, pPoint, imagePoint); if (SCPtest) { bool SCPinImage = (pPoint.x() > 0 && pPoint.x() < image_width) && (pPoint.y() > 0 && pPoint.y() < image_height); if (SCPinImage) { painter->fillRect(pPoint.x() * scale - 2, pPoint.y() * scale - 2, 4, 4, KStarsData::Instance()->colorScheme()->colorNamed("TargetColor")); painter->drawText(pPoint.x() * scale + 15, pPoint.y() * scale + 15, i18nc("South Celestial Pole", "SCP")); } } } } } bool FITSView::pointIsInImage(QPointF pt, bool scaled) { float scale = (currentZoom / ZOOM_DEFAULT); if (scaled) return pt.x() < image_width * scale && pt.y() < image_height * scale && pt.x() > 0 && pt.y() > 0; else return pt.x() < image_width && pt.y() < image_height && pt.x() > 0 && pt.y() > 0; } QPointF FITSView::getPointForGridLabel() { float scale = (currentZoom / ZOOM_DEFAULT); //These get the maximum X and Y points in the list that are in the image QPointF maxXPt(image_width * scale / 2, image_height * scale / 2); for (auto &p : eqGridPoints) { if (p.x() > maxXPt.x() && pointIsInImage(p, true)) maxXPt = p; } QPointF maxYPt(image_width * scale / 2, image_height * scale / 2); for (auto &p : eqGridPoints) { if (p.y() > maxYPt.y() && pointIsInImage(p, true)) maxYPt = p; } QPointF minXPt(image_width * scale / 2, image_height * scale / 2); for (auto &p : eqGridPoints) { if (p.x() < minXPt.x() && pointIsInImage(p, true)) minXPt = p; } QPointF minYPt(image_width * scale / 2, image_height * scale / 2); for (auto &p : eqGridPoints) { if (p.y() < minYPt.y() && pointIsInImage(p, true)) minYPt = p; } //This gives preference to points that are on the right hand side and bottom. //But if the line doesn't intersect the right or bottom, it then tries for the top and left. //If no points are found in the image, it returns a point off the screen //If all else fails, like in the case of a circle on the image, it returns the far right point. if (image_width * scale - maxXPt.x() < 10) { return QPointF( image_width * scale - 50, maxXPt.y() - 10); //This will draw the text on the right hand side, up and to the left of the point where the line intersects } if (image_height * scale - maxYPt.y() < 10) return QPointF( maxYPt.x() - 40, image_height * scale - 10); //This will draw the text on the bottom side, up and to the left of the point where the line intersects if (minYPt.y() * scale < 30) return QPointF( minYPt.x() + 10, 20); //This will draw the text on the top side, down and to the right of the point where the line intersects if (minXPt.x() * scale < 30) return QPointF( 10, minXPt.y() + 20); //This will draw the text on the left hand side, down and to the right of the point where the line intersects if (maxXPt.x() == image_width * scale / 2 && maxXPt.y() == image_height * scale / 2) return QPointF(-100, -100); //All of the points were off the screen return QPoint(maxXPt.x() - 40, maxXPt.y() - 10); } void FITSView::setFirstLoad(bool value) { firstLoad = value; } QPixmap &FITSView::getTrackingBoxPixmap(uint8_t margin) { if (trackingBox.isNull()) return trackingBoxPixmap; int x1 = (trackingBox.x() - margin) * (currentZoom / ZOOM_DEFAULT); int y1 = (trackingBox.y() - margin) * (currentZoom / ZOOM_DEFAULT); int w = (trackingBox.width() + margin * 2) * (currentZoom / ZOOM_DEFAULT); int h = (trackingBox.height() + margin * 2) * (currentZoom / ZOOM_DEFAULT); trackingBoxPixmap = image_frame->grab(QRect(x1, y1, w, h)); return trackingBoxPixmap; } void FITSView::setTrackingBox(const QRect &rect) { if (rect != trackingBox) { trackingBox = rect; updateFrame(); if(showStarProfile) viewStarProfile(); } } void FITSView::resizeTrackingBox(int newSize) { int x = trackingBox.x() + trackingBox.width() / 2; int y = trackingBox.y() + trackingBox.height() / 2; int delta = newSize / 2; setTrackingBox(QRect( x - delta, y - delta, newSize, newSize)); } bool FITSView::isCrosshairShown() { return showCrosshair; } bool FITSView::isEQGridShown() { return showEQGrid; } bool FITSView::areObjectsShown() { return showObjects; } bool FITSView::isPixelGridShown() { return showPixelGrid; } void FITSView::toggleCrosshair() { showCrosshair = !showCrosshair; updateFrame(); } void FITSView::toggleEQGrid() { showEQGrid = !showEQGrid; if (!imageData->isWCSLoaded() && !wcsWatcher.isRunning()) { QFuture future = QtConcurrent::run(imageData, &FITSData::loadWCS); wcsWatcher.setFuture(future); return; } if (image_frame != nullptr) updateFrame(); } void FITSView::toggleObjects() { showObjects = !showObjects; if (!imageData->isWCSLoaded() && !wcsWatcher.isRunning()) { QFuture future = QtConcurrent::run(imageData, &FITSData::loadWCS); wcsWatcher.setFuture(future); return; } if (image_frame != nullptr) updateFrame(); } void FITSView::toggleStars() { toggleStars(!markStars); if (image_frame != nullptr) updateFrame(); } void FITSView::toggleStarProfile() { #ifdef HAVE_DATAVISUALIZATION showStarProfile = !showStarProfile; if(showStarProfile && trackingBoxEnabled) viewStarProfile(); if(toggleProfileAction) toggleProfileAction->setChecked(showStarProfile); if(showStarProfile) { //The tracking box is already on for Guide and Focus Views, but off for Normal and Align views. //So for Normal and Align views, we need to set up the tracking box. if(mode == FITS_NORMAL || mode == FITS_ALIGN) { setCursorMode(selectCursor); connect(this, SIGNAL(trackingStarSelected(int, int)), this, SLOT(move3DTrackingBox(int, int))); trackingBox = QRect(0, 0, 128, 128); setTrackingBoxEnabled(true); if(starProfileWidget) connect(starProfileWidget, SIGNAL(sampleSizeUpdated(int)), this, SLOT(resizeTrackingBox(int))); } if(starProfileWidget) connect(starProfileWidget, SIGNAL(rejected()), this, SLOT(toggleStarProfile())); } else { //This shuts down the tracking box for Normal and Align Views //It doesn't touch Guide and Focus Views because they still need a tracking box if(mode == FITS_NORMAL || mode == FITS_ALIGN) { if(getCursorMode() == selectCursor) setCursorMode(dragCursor); disconnect(this, SIGNAL(trackingStarSelected(int, int)), this, SLOT(move3DTrackingBox(int, int))); setTrackingBoxEnabled(false); if(starProfileWidget) disconnect(starProfileWidget, SIGNAL(sampleSizeUpdated(int)), this, SLOT(resizeTrackingBox(int))); } if(starProfileWidget) { disconnect(starProfileWidget, SIGNAL(rejected()), this, SLOT(toggleStarProfile())); starProfileWidget->close(); starProfileWidget = nullptr; } emit starProfileWindowClosed(); } updateFrame(); #endif } void FITSView::move3DTrackingBox(int x, int y) { int boxSize = trackingBox.width(); QRect starRect = QRect(x - boxSize / 2, y - boxSize / 2, boxSize, boxSize); setTrackingBox(starRect); } void FITSView::viewStarProfile() { #ifdef HAVE_DATAVISUALIZATION if(!trackingBoxEnabled) { setTrackingBoxEnabled(true); setTrackingBox(QRect(0, 0, 128, 128)); } if(!starProfileWidget) { starProfileWidget = new StarProfileViewer(this); //This is a band-aid to fix a QT bug with createWindowContainer //It will set the cursor of the Window containing the view that called the Star Profile method to the Arrow Cursor //Note that Ekos Manager is a QDialog and FitsViewer is a KXmlGuiWindow QWidget * superParent = this->parentWidget(); while(superParent->parentWidget() != 0 && !superParent->inherits("QDialog") && !superParent->inherits("KXmlGuiWindow")) superParent = superParent->parentWidget(); superParent->setCursor(Qt::ArrowCursor); //This is the end of the band-aid connect(starProfileWidget, SIGNAL(rejected()), this, SLOT(toggleStarProfile())); if(mode == FITS_ALIGN || mode == FITS_NORMAL) { starProfileWidget->enableTrackingBox(true); imageData->setStarAlgorithm(ALGORITHM_CENTROID); connect(starProfileWidget, SIGNAL(sampleSizeUpdated(int)), this, SLOT(resizeTrackingBox(int))); } } QList starCenters = imageData->getStarCentersInSubFrame(trackingBox); if(starCenters.size() == 0) { // FIXME, the following does not work anymore. //imageData->findStars(&trackingBox, true); // FIXME replacing it with this imageData->findStars(ALGORITHM_CENTROID, trackingBox); starCenters = imageData->getStarCentersInSubFrame(trackingBox); } starProfileWidget->loadData(imageData, trackingBox, starCenters); starProfileWidget->show(); starProfileWidget->raise(); if(markStars) updateFrame(); //this is to update for the marked stars #endif } void FITSView::togglePixelGrid() { showPixelGrid = !showPixelGrid; updateFrame(); } int FITSView::findStars(StarAlgorithm algorithm, const QRect &searchBox) { int count = 0; if(trackingBoxEnabled) count = imageData->findStars(algorithm, trackingBox); else count = imageData->findStars(algorithm, searchBox); return count; } void FITSView::toggleStars(bool enable) { markStars = enable; if (markStars && !imageData->areStarsSearched()) { QApplication::setOverrideCursor(Qt::WaitCursor); emit newStatus(i18n("Finding stars..."), FITS_MESSAGE); qApp->processEvents(); int count = findStars(); if (count >= 0 && isVisible()) emit newStatus(i18np("1 star detected.", "%1 stars detected.", count), FITS_MESSAGE); QApplication::restoreOverrideCursor(); } } void FITSView::processPointSelection(int x, int y) { - //if (mode != FITS_GUIDE) - //return; - - //image_data->getCenterSelection(&x, &y); - - //setGuideSquare(x,y); emit trackingStarSelected(x, y); } void FITSView::processMarkerSelection(int x, int y) { markerCrosshair.setX(x); markerCrosshair.setY(y); updateFrame(); } void FITSView::setTrackingBoxEnabled(bool enable) { if (enable != trackingBoxEnabled) { trackingBoxEnabled = enable; //updateFrame(); } } void FITSView::wheelEvent(QWheelEvent * event) { //This attempts to send the wheel event back to the Scroll Area if it was taken from a trackpad //It should still do the zoom if it is a mouse wheel if (event->source() == Qt::MouseEventSynthesizedBySystem) { QScrollArea::wheelEvent(event); } else { QPoint mouseCenter = getImagePoint(event->pos()); if (event->angleDelta().y() > 0) ZoomIn(); else ZoomOut(); event->accept(); cleanUpZoom(mouseCenter); } } /** This method is intended to keep key locations in an image centered on the screen while zooming. If there is a marker or tracking box, it centers on those. If not, it uses the point called viewCenter that was passed as a parameter. */ void FITSView::cleanUpZoom(QPoint viewCenter) { int x0 = 0; int y0 = 0; double scale = (currentZoom / ZOOM_DEFAULT); if (!markerCrosshair.isNull()) { x0 = markerCrosshair.x() * scale; y0 = markerCrosshair.y() * scale; } else if (trackingBoxEnabled) { x0 = trackingBox.center().x() * scale; y0 = trackingBox.center().y() * scale; } else { x0 = viewCenter.x() * scale; y0 = viewCenter.y() * scale; } ensureVisible(x0, y0, width() / 2, height() / 2); updateMouseCursor(); } /** This method converts a point from the ViewPort Coordinate System to the Image Coordinate System. */ QPoint FITSView::getImagePoint(QPoint viewPortPoint) { QWidget * w = widget(); if (w == nullptr) return QPoint(0, 0); double scale = (currentZoom / ZOOM_DEFAULT); QPoint widgetPoint = w->mapFromParent(viewPortPoint); QPoint imagePoint = QPoint(widgetPoint.x() / scale, widgetPoint.y() / scale); return imagePoint; } void FITSView::initDisplayImage() { if (imageData->channels() == 1) { rawImage = QImage(image_width, image_height, QImage::Format_Indexed8); rawImage.setColorCount(256); for (int i = 0; i < 256; i++) rawImage.setColor(i, qRgb(i, i, i)); } else { rawImage = QImage(image_width, image_height, QImage::Format_RGB32); } } /** The Following two methods allow gestures to work with trackpads. Specifically, we are targeting the pinch events, so that if one is generated, Then the pinchTriggered method will be called. If the event is not a pinch gesture, then the event is passed back to the other event handlers. */ bool FITSView::event(QEvent * event) { if (event->type() == QEvent::Gesture) return gestureEvent(dynamic_cast(event)); return QScrollArea::event(event); } bool FITSView::gestureEvent(QGestureEvent * event) { if (QGesture * pinch = event->gesture(Qt::PinchGesture)) pinchTriggered(dynamic_cast(pinch)); return true; } /** This Method works with Trackpads to use the pinch gesture to scroll in and out It stores a point to keep track of the location where the gesture started so that while you are zooming, it tries to keep that initial point centered in the view. **/ void FITSView::pinchTriggered(QPinchGesture * gesture) { if (!zooming) { zoomLocation = getImagePoint(mapFromGlobal(QCursor::pos())); zooming = true; } if (gesture->state() == Qt::GestureFinished) { zooming = false; } zoomTime++; //zoomTime is meant to slow down the zooming with a pinch gesture. if (zoomTime > 10000) //This ensures zoomtime never gets too big. zoomTime = 0; if (zooming && (zoomTime % 10 == 0)) //zoomTime is set to slow it by a factor of 10. { if (gesture->totalScaleFactor() > 1) ZoomIn(); else ZoomOut(); } cleanUpZoom(zoomLocation); } /*void FITSView::handleWCSCompletion() { //bool hasWCS = wcsWatcher.result(); if(imageData->hasWCS()) this->updateFrame(); emit wcsToggled(imageData->hasWCS()); }*/ void FITSView::syncWCSState() { bool hasWCS = imageData->hasWCS(); bool wcsLoaded = imageData->isWCSLoaded(); if (hasWCS && wcsLoaded) this->updateFrame(); emit wcsToggled(hasWCS); if (toggleEQGridAction != nullptr) toggleEQGridAction->setEnabled(hasWCS); if (toggleObjectsAction != nullptr) toggleObjectsAction->setEnabled(hasWCS); if (centerTelescopeAction != nullptr) centerTelescopeAction->setEnabled(hasWCS); } void FITSView::createFloatingToolBar() { if (floatingToolBar != nullptr) return; floatingToolBar = new QToolBar(this); auto * eff = new QGraphicsOpacityEffect(this); floatingToolBar->setGraphicsEffect(eff); eff->setOpacity(0.2); floatingToolBar->setVisible(false); floatingToolBar->setStyleSheet( "QToolBar{background: rgba(150, 150, 150, 210); border:none; color: yellow}" "QToolButton{background: transparent; border:none; color: yellow}" "QToolButton:hover{background: rgba(200, 200, 200, 255);border:solid; color: yellow}" "QToolButton:checked{background: rgba(110, 110, 110, 255);border:solid; color: yellow}"); floatingToolBar->setFloatable(true); floatingToolBar->setIconSize(QSize(25, 25)); //floatingToolBar->setMovable(true); QAction * action = nullptr; floatingToolBar->addAction(QIcon::fromTheme("zoom-in"), i18n("Zoom In"), this, SLOT(ZoomIn())); floatingToolBar->addAction(QIcon::fromTheme("zoom-out"), i18n("Zoom Out"), this, SLOT(ZoomOut())); floatingToolBar->addAction(QIcon::fromTheme("zoom-fit-best"), i18n("Default Zoom"), this, SLOT(ZoomDefault())); floatingToolBar->addAction(QIcon::fromTheme("zoom-fit-width"), i18n("Zoom to Fit"), this, SLOT(ZoomToFit())); floatingToolBar->addSeparator(); action = floatingToolBar->addAction(QIcon::fromTheme("crosshairs"), i18n("Show Cross Hairs"), this, SLOT(toggleCrosshair())); action->setCheckable(true); action = floatingToolBar->addAction(QIcon::fromTheme("map-flat"), i18n("Show Pixel Gridlines"), this, SLOT(togglePixelGrid())); action->setCheckable(true); toggleStarsAction = floatingToolBar->addAction(QIcon::fromTheme("kstars_stars"), i18n("Detect Stars in Image"), this, SLOT(toggleStars())); toggleStarsAction->setCheckable(true); #ifdef HAVE_DATAVISUALIZATION toggleProfileAction = floatingToolBar->addAction(QIcon::fromTheme("star-profile", QIcon(":/icons/star_profile.svg")), i18n("View Star Profile"), this, SLOT(toggleStarProfile())); toggleProfileAction->setCheckable(true); #endif if (mode == FITS_NORMAL || mode == FITS_ALIGN) { floatingToolBar->addSeparator(); toggleEQGridAction = floatingToolBar->addAction(QIcon::fromTheme("kstars_grid"), i18n("Show Equatorial Gridlines"), this, SLOT(toggleEQGrid())); toggleEQGridAction->setCheckable(true); toggleEQGridAction->setEnabled(false); toggleObjectsAction = floatingToolBar->addAction(QIcon::fromTheme("help-hint"), i18n("Show Objects in Image"), this, SLOT(toggleObjects())); toggleObjectsAction->setCheckable(true); toggleEQGridAction->setEnabled(false); centerTelescopeAction = floatingToolBar->addAction(QIcon::fromTheme("center_telescope", QIcon(":/icons/center_telescope.svg")), i18n("Center Telescope"), this, SLOT(centerTelescope())); centerTelescopeAction->setCheckable(true); centerTelescopeAction->setEnabled(false); } } /** This methood either enables or disables the scope mouse mode so you can slew your scope to coordinates just by clicking the mouse on a spot in the image. */ void FITSView::centerTelescope() { if (imageHasWCS()) { if (getCursorMode() == FITSView::scopeCursor) { setCursorMode(lastMouseMode); } else { lastMouseMode = getCursorMode(); setCursorMode(FITSView::scopeCursor); } updateFrame(); } updateScopeButton(); } void FITSView::updateScopeButton() { if (centerTelescopeAction != nullptr) { if (getCursorMode() == FITSView::scopeCursor) { centerTelescopeAction->setChecked(true); } else { centerTelescopeAction->setChecked(false); } } } /** This method just verifies if INDI is online, a telescope present, and is connected */ bool FITSView::isTelescopeActive() { #ifdef HAVE_INDI if (INDIListener::Instance()->size() == 0) { return false; } foreach (ISD::GDInterface * gd, INDIListener::Instance()->getDevices()) { INDI::BaseDevice * bd = gd->getBaseDevice(); if (gd->getType() != KSTARS_TELESCOPE) continue; if (bd == nullptr) continue; return bd->isConnected(); } return false; #else return false; #endif } void FITSView::setStarsEnabled(bool enable) { markStars = enable; if (floatingToolBar != nullptr) { foreach (QAction * action, floatingToolBar->actions()) { if (action->text() == i18n("Detect Stars in Image")) { action->setChecked(markStars); break; } } } } void FITSView::setStarsHFREnabled(bool enable) { showStarsHFR = enable; }