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