1 /*
2 * usimage.cpp
3 * PHD Guiding
4 *
5 * Created by Craig Stark.
6 * Copyright (c) 2006-2010 Craig Stark.
7 * All rights reserved.
8 *
9 * This source code is distributed under the following "BSD" license
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions are met:
12 * Redistributions of source code must retain the above copyright notice,
13 * this list of conditions and the following disclaimer.
14 * Redistributions in binary form must reproduce the above copyright notice,
15 * this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 * Neither the name of Craig Stark, Stark Labs nor the names of its
18 * contributors may be used to endorse or promote products derived from
19 * this software without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 */
34
35 #include "phd.h"
36 #include "image_math.h"
37
38 #include <algorithm>
39
Init(const wxSize & size)40 bool usImage::Init(const wxSize& size)
41 {
42 // Allocates space for image and sets params up
43 // returns true on error
44
45 unsigned int prev = NPixels;
46 NPixels = size.GetWidth() * size.GetHeight();
47 Size = size;
48 Subframe = wxRect(0, 0, 0, 0);
49 MinADU = MaxADU = MedianADU = 0;
50
51 if (NPixels != prev)
52 {
53 delete[] ImageData;
54
55 if (NPixels)
56 {
57 ImageData = new unsigned short[NPixels];
58 if (!ImageData)
59 {
60 NPixels = 0;
61 return true;
62 }
63 }
64 else
65 ImageData = nullptr;
66 }
67
68 return false;
69 }
70
SwapImageData(usImage & other)71 void usImage::SwapImageData(usImage& other)
72 {
73 unsigned short *t = ImageData;
74 ImageData = other.ImageData;
75 other.ImageData = t;
76 }
77
CalcStats()78 void usImage::CalcStats()
79 {
80 if (!ImageData || !NPixels)
81 return;
82
83 MinADU = 65535; MaxADU = 0;
84 FiltMin = 65535; FiltMax = 0;
85
86 if (Subframe.IsEmpty())
87 {
88 // full frame, no subframe
89
90 for (unsigned int i = 0; i < NPixels; i++)
91 {
92 unsigned short d = ImageData[i];
93 if (d < MinADU) MinADU = d;
94 if (d > MaxADU) MaxADU = d;
95 }
96
97 unsigned short *tmpdata = new unsigned short[NPixels];
98
99 Median3(tmpdata, ImageData, Size, wxRect(Size));
100
101 const unsigned short *src = tmpdata;
102 for (unsigned int i = 0; i < NPixels; i++)
103 {
104 unsigned short d = *src++;
105 if (d < FiltMin) FiltMin = d;
106 if (d > FiltMax) FiltMax = d;
107 }
108
109 memcpy(tmpdata, ImageData, NPixels * sizeof(unsigned short));
110 std::nth_element(tmpdata, tmpdata + NPixels / 2, tmpdata + NPixels);
111 MedianADU = tmpdata[NPixels / 2];
112
113 delete[] tmpdata;
114 }
115 else
116 {
117 // Subframe
118
119 unsigned int pixcnt = Subframe.width * Subframe.height;
120 unsigned short *tmpdata = new unsigned short[pixcnt];
121
122 unsigned short *dst;
123
124 dst = tmpdata;
125 for (int y = 0; y < Subframe.height; y++)
126 {
127 const unsigned short *src = ImageData + Subframe.x + (Subframe.y + y) * Size.GetWidth();
128 for (int x = 0; x < Subframe.width; x++)
129 {
130 unsigned short d = *src;
131 if (d < MinADU) MinADU = d;
132 if (d > MaxADU) MaxADU = d;
133 *dst++ = *src++;
134 }
135 }
136
137 dst = new unsigned short[pixcnt];
138
139 Median3(dst, tmpdata, Subframe.GetSize(), wxRect(Subframe.GetSize()));
140
141 const unsigned short *src = dst;
142 for (unsigned int i = 0; i < pixcnt; i++)
143 {
144 unsigned short d = *src++;
145 if (d < FiltMin) FiltMin = d;
146 if (d > FiltMax) FiltMax = d;
147 }
148
149 memcpy(dst, tmpdata, pixcnt * sizeof(unsigned short));
150 std::nth_element(dst, dst + pixcnt / 2, dst + pixcnt);
151 MedianADU = dst[pixcnt / 2];
152
153 delete[] dst;
154 delete[] tmpdata;
155 }
156 }
157
CopyToImage(wxImage ** rawimg,int blevel,int wlevel,double power)158 bool usImage::CopyToImage(wxImage **rawimg, int blevel, int wlevel, double power)
159 {
160 wxImage *img = *rawimg;
161
162 if (!img || !img->Ok() || (img->GetWidth() != Size.GetWidth()) || (img->GetHeight() != Size.GetHeight()) ) // can't reuse bitmap
163 {
164 delete img;
165 img = new wxImage(Size.GetWidth(), Size.GetHeight(), false);
166 }
167
168 unsigned char *ImgPtr = img->GetData();
169 unsigned short *RawPtr = ImageData;
170
171 if (power == 1.0 || blevel >= wlevel)
172 {
173 float range = (float) wxMax(1, wlevel); // Go 0-max
174 for (unsigned int i = 0; i < NPixels; i++, RawPtr++)
175 {
176 float d;
177 if (*RawPtr >= range)
178 d = 255.0;
179 else
180 d = ((float) (*RawPtr) / range) * 255.0;
181
182 *ImgPtr++ = (unsigned char) d;
183 *ImgPtr++ = (unsigned char) d;
184 *ImgPtr++ = (unsigned char) d;
185 }
186 }
187 else
188 {
189 float range = (float) (wlevel - blevel);
190 for (unsigned int i = 0; i < NPixels; i++, RawPtr++ )
191 {
192 float d;
193 if (*RawPtr <= blevel)
194 d = 0.0;
195 else if (*RawPtr >= wlevel)
196 d = 255.0;
197 else
198 {
199 d = ((float) (*RawPtr) - (float) blevel) / range;
200 d = pow(d, (float) power) * 255.0;
201 }
202 *ImgPtr++ = (unsigned char) d;
203 *ImgPtr++ = (unsigned char) d;
204 *ImgPtr++ = (unsigned char) d;
205 }
206 }
207
208 *rawimg = img;
209 return false;
210 }
211
BinnedCopyToImage(wxImage ** rawimg,int blevel,int wlevel,double power)212 bool usImage::BinnedCopyToImage(wxImage **rawimg, int blevel, int wlevel, double power)
213 {
214 wxImage *img;
215 unsigned char *ImgPtr;
216 unsigned short *RawPtr;
217 int x,y;
218 float d;
219 //, s_factor;
220 int full_xsize, full_ysize;
221 int use_xsize, use_ysize;
222
223 full_xsize = Size.GetWidth();
224 full_ysize = Size.GetHeight();
225 use_xsize = full_xsize;
226 use_ysize = full_ysize;
227 if (use_xsize % 2) use_xsize--;
228 if (use_ysize % 2) use_ysize--;
229 // Binsize = 2;
230
231 img = *rawimg;
232 if ((!img->Ok()) || (img->GetWidth() != (full_xsize/2)) || (img->GetHeight() != (full_ysize/2)) ) { // can't reuse bitmap
233 delete img; // Clear out current image if it exists
234 img = new wxImage(full_xsize/2, full_ysize/2, false);
235 }
236 ImgPtr = img->GetData();
237 RawPtr = ImageData;
238 // s_factor = (((float) Max - (float) Min) / 255.0);
239 float range = (float) (wlevel - blevel);
240
241 if ((power == 1.0) || (range == 0.0)) {
242 range = wlevel; // Go 0-max
243 if (range == 0.0) range = 0.001;
244 for (y=0; y<use_ysize; y+=2) {
245 for (x=0; x<use_xsize; x+=2) {
246 RawPtr = ImageData + x + y*full_xsize;
247 d = (float) (*RawPtr + *(RawPtr+1) + *(RawPtr+full_xsize) + *(RawPtr+1+full_xsize)) / 4.0;
248 d = (d / range) * 255.0;
249 if (d < 0.0) d = 0.0;
250 else if (d > 255.0) d = 255.0;
251 *ImgPtr = (unsigned char) d;
252 ImgPtr++;
253 *ImgPtr = (unsigned char) d;
254 ImgPtr++;
255 *ImgPtr = (unsigned char) d;
256 ImgPtr++;
257 }
258 }
259 }
260 else {
261 for (y=0; y<use_ysize; y+=2) {
262 for (x=0; x<use_xsize; x+=2) {
263 RawPtr = ImageData + x + y*full_xsize;
264 d = (float) (*RawPtr + *(RawPtr+1) + *(RawPtr+full_xsize) + *(RawPtr+1+full_xsize)) / 4.0;
265 d = (d - (float) blevel) / range ;
266 if (d < 0.0) d= 0.0;
267 else if (d > 1.0) d = 1.0;
268 d = pow(d, (float) power) * 255.0;
269 *ImgPtr = (unsigned char) d;
270 ImgPtr++;
271 *ImgPtr = (unsigned char) d;
272 ImgPtr++;
273 *ImgPtr = (unsigned char) d;
274 ImgPtr++;
275 }
276 }
277 }
278 *rawimg = img;
279 return false;
280 }
281
InitImgStartTime()282 void usImage::InitImgStartTime()
283 {
284 ImgStartTime = wxDateTime::UNow();
285 }
286
Save(const wxString & fname,const wxString & hdrNote) const287 bool usImage::Save(const wxString& fname, const wxString& hdrNote) const
288 {
289 bool bError = false;
290
291 try
292 {
293 fitsfile *fptr; // FITS file pointer
294 int status = 0; // CFITSIO status value MUST be initialized to zero!
295
296 PHD_fits_create_file(&fptr, fname, true, &status);
297
298 long fsize[] = {
299 (long) Size.GetWidth(),
300 (long) Size.GetHeight(),
301 };
302 fits_create_img(fptr, USHORT_IMG, 2, fsize, &status);
303
304 FITSHdrWriter hdr(fptr, &status);
305
306 float exposure = (float) ImgExpDur / 1000.0;
307 hdr.write("EXPOSURE", exposure, "Exposure time in seconds");
308
309 if (ImgStackCnt > 1)
310 hdr.write("STACKCNT", (unsigned int) ImgStackCnt, "Stacked frame count");
311
312 if (!hdrNote.IsEmpty())
313 hdr.write("USERNOTE", hdrNote.utf8_str(), 0);
314
315 hdr.write("DATE", wxDateTime::UNow(), wxDateTime::UTC, "file creation time, UTC");
316 hdr.write("DATE-OBS", ImgStartTime, wxDateTime::UTC, "Image capture start time, UTC");
317 hdr.write("CREATOR", wxString(APPNAME _T(" ") FULLVER).c_str(), "Capture software");
318 hdr.write("PHDPROFI", pConfig->GetCurrentProfile().c_str(), "PHD2 Equipment Profile");
319
320 if (pCamera)
321 {
322 hdr.write("INSTRUME", pCamera->Name.c_str(), "Instrument name");
323 unsigned int b = pCamera->Binning;
324 hdr.write("XBINNING", b, "Camera X Bin");
325 hdr.write("YBINNING", b, "Camera Y Bin");
326 hdr.write("CCDXBIN", b, "Camera X Bin");
327 hdr.write("CCDYBIN", b, "Camera Y Bin");
328 float sz = b * pCamera->GetCameraPixelSize();
329 hdr.write("XPIXSZ", sz, "pixel size in microns (with binning)");
330 hdr.write("YPIXSZ", sz, "pixel size in microns (with binning)");
331 unsigned int g = (unsigned int) pCamera->GuideCameraGain;
332 hdr.write("GAIN", g, "PHD Gain Value (0-100)");
333 unsigned int bpp = pCamera->BitsPerPixel();
334 hdr.write("CAMBPP", bpp, "Camera resolution, bits per pixel");
335 }
336
337 if (pPointingSource)
338 {
339 double ra, dec, st;
340 bool err = pPointingSource->GetCoordinates(&ra, &dec, &st);
341 if (!err)
342 {
343 hdr.write("RA", (float) (ra * 360.0 / 24.0), "Object Right Ascension in degrees");
344 hdr.write("DEC", (float) dec, "Object Declination in degrees");
345
346 {
347 int h = (int) ra;
348 ra -= h;
349 ra *= 60.0;
350 int m = (int) ra;
351 ra -= m;
352 ra *= 60.0;
353 hdr.write("OBJCTRA", wxString::Format("%02d %02d %06.3f", h, m, ra).c_str(), "Object Right Ascension in hms");
354 }
355
356 {
357 int sign = dec < 0.0 ? -1 : +1;
358 dec *= sign;
359 int d = (int) dec;
360 dec -= d;
361 dec *= 60.0;
362 int m = (int) dec;
363 dec -= m;
364 dec *= 60.0;
365 hdr.write("OBJCTDEC", wxString::Format("%c%d %02d %06.3f", sign < 0 ? '-' : '+', d, m, dec).c_str(), "Object Declination in dms");
366 }
367 }
368
369 PierSide p = pPointingSource->SideOfPier();
370 if (p != PierSide::PIER_SIDE_UNKNOWN)
371 hdr.write("PIERSIDE", (unsigned int) p, "Side of Pier 0=East 1=West");
372 }
373
374 float sc = (float) pFrame->GetCameraPixelScale();
375 hdr.write("SCALE", sc, "Image scale (arcsec / pixel)");
376 hdr.write("PIXSCALE", sc, "Image scale (arcsec / pixel)");
377 hdr.write("PEDESTAL", (unsigned int) Pedestal, "dark subtraction bias value");
378 hdr.write("SATURATE", (1U << BitsPerPixel) - 1, "Data value at which saturation occurs");
379
380 const PHD_Point& lockPos = pFrame->pGuider->LockPosition();
381 if (lockPos.IsValid())
382 {
383 hdr.write("PHDLOCKX", (float) lockPos.X, "PHD2 lock position x");
384 hdr.write("PHDLOCKY", (float) lockPos.Y, "PHD2 lock position y");
385 }
386
387 if (!Subframe.IsEmpty())
388 {
389 hdr.write("PHDSUBFX", (unsigned int) Subframe.x, "PHD2 subframe x");
390 hdr.write("PHDSUBFY", (unsigned int) Subframe.y, "PHD2 subframe y");
391 hdr.write("PHDSUBFW", (unsigned int) Subframe.width, "PHD2 subframe width");
392 hdr.write("PHDSUBFH", (unsigned int) Subframe.height, "PHD2 subframe height");
393 }
394
395 long fpixel[3] = { 1, 1, 1 };
396 fits_write_pix(fptr, TUSHORT, fpixel, NPixels, ImageData, &status);
397
398 PHD_fits_close_file(fptr);
399
400 bError = status ? true : false;
401 }
402 catch (const wxString& Msg)
403 {
404 POSSIBLY_UNUSED(Msg);
405 bError = true;
406 }
407
408 return bError;
409 }
410
fhdr_int(fitsfile * fptr,const char * key,int * val)411 static bool fhdr_int(fitsfile *fptr, const char *key, int *val)
412 {
413 char *k = const_cast<char *>(key);
414 int status = 0;
415 fits_read_key(fptr, TINT, k, val, nullptr, &status);
416 return status == 0;
417 }
418
Load(const wxString & fname)419 bool usImage::Load(const wxString& fname)
420 {
421 bool bError = false;
422
423 try
424 {
425 if (!wxFileExists(fname))
426 {
427 pFrame->Alert(_("File does not exist - cannot load ") + fname);
428 throw ERROR_INFO("File does not exist");
429 }
430
431 int status = 0; // CFITSIO status value MUST be initialized to zero!
432 fitsfile *fptr; // FITS file pointer
433 if (!PHD_fits_open_diskfile(&fptr, fname, READONLY, &status))
434 {
435 int hdutype;
436 if (fits_get_hdu_type(fptr, &hdutype, &status) || hdutype != IMAGE_HDU)
437 {
438 pFrame->Alert(_("FITS file is not of an image: ") + fname);
439 throw ERROR_INFO("Fits file is not an image");
440 }
441
442 // Get HDUs and size
443 int naxis = 0;
444 fits_get_img_dim(fptr, &naxis, &status);
445 long fsize[3];
446 fits_get_img_size(fptr, 2, fsize, &status);
447 int nhdus = 0;
448 fits_get_num_hdus(fptr, &nhdus, &status);
449 if ((nhdus != 1) || (naxis != 2)) {
450 pFrame->Alert(wxString::Format(_("Unsupported type or read error loading FITS file %s"), fname));
451 throw ERROR_INFO("unsupported type");
452 }
453 if (Init((int) fsize[0], (int) fsize[1]))
454 {
455 pFrame->Alert(wxString::Format(_("Memory allocation error loading FITS file %s"), fname));
456 throw ERROR_INFO("Memory Allocation failure");
457 }
458 long fpixel[3] = { 1, 1, 1 };
459 if (fits_read_pix(fptr, TUSHORT, fpixel, (int)(fsize[0] * fsize[1]), nullptr, ImageData, nullptr, &status)) { // Read image
460 pFrame->Alert(wxString::Format(_("Error reading data from FITS file %s"), fname));
461 throw ERROR_INFO("Error reading");
462 }
463
464 char *key = const_cast<char *>("EXPOSURE");
465 float exposure;
466 status = 0;
467 fits_read_key(fptr, TFLOAT, key, &exposure, nullptr, &status);
468 if (status == 0)
469 ImgExpDur = (int) (exposure * 1000.0);
470
471 int stackcnt;
472 if (fhdr_int(fptr, "STACKCNT", &stackcnt))
473 ImgStackCnt = stackcnt;
474
475 int pedestal;
476 if (fhdr_int(fptr, "PEDESTAL", &pedestal))
477 Pedestal = (unsigned short) pedestal;
478
479 int saturate;
480 if (fhdr_int(fptr, "SATURATE", &saturate))
481 BitsPerPixel = saturate > 255 ? 16 : 8;
482
483 wxRect subf;
484 bool ok = fhdr_int(fptr, "PHDSUBFX", &subf.x);
485 if (ok) ok = fhdr_int(fptr, "PHDSUBFY", &subf.y);
486 if (ok) ok = fhdr_int(fptr, "PHDSUBFW", &subf.width);
487 if (ok) ok = fhdr_int(fptr, "PHDSUBFH", &subf.height);
488 if (ok) Subframe = subf;
489
490 PHD_fits_close_file(fptr);
491
492 CalcStats();
493 }
494 else
495 {
496 pFrame->Alert(wxString::Format(_("Error opening FITS file %s"), fname));
497 throw ERROR_INFO("error opening file");
498 }
499 }
500 catch (const wxString& Msg)
501 {
502 POSSIBLY_UNUSED(Msg);
503 bError = true;
504 }
505
506 return bError;
507 }
508
CopyFrom(const usImage & src)509 bool usImage::CopyFrom(const usImage& src)
510 {
511 if (Init(src.Size))
512 return true;
513 memcpy(ImageData, src.ImageData, NPixels * sizeof(unsigned short));
514 return false;
515 }
516
Rotate(double theta,bool mirror)517 bool usImage::Rotate(double theta, bool mirror)
518 {
519 wxImage *pImg = 0;
520
521 CalcStats();
522
523 CopyToImage(&pImg, MinADU, MaxADU, 1.0);
524
525 wxImage mirrored = *pImg;
526
527 if (mirror)
528 {
529 mirrored = pImg->Mirror(false);
530 }
531
532 wxImage rotated = mirrored.Rotate(theta, wxPoint(0,0));
533
534 CopyFromImage(rotated);
535
536 delete pImg;
537
538 return false;
539 }
540
CopyFromImage(const wxImage & img)541 bool usImage::CopyFromImage(const wxImage& img)
542 {
543 Init(img.GetSize());
544
545 const unsigned char *pSrc = img.GetData();
546 unsigned short *pDest = ImageData;
547
548 for (unsigned int i = 0; i < NPixels; i++)
549 {
550 *pDest++ = ((unsigned short) *pSrc) << 8;
551 pSrc += 3;
552 }
553
554 return false;
555 }
556