fitsdata.cpp 56.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/***************************************************************************
                          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.                  *
 ***************************************************************************/

Jasem Mutlaq's avatar
Jasem Mutlaq committed
20
#include <config-kstars.h>
21
#include "fitsdata.h"
22

23 24
#include <cmath>
#include <cstdlib>
25
#include <climits>
26

27
#include <QApplication>
Jasem Mutlaq's avatar
Jasem Mutlaq committed
28
#include <QLocale>
29
#include <QFile>
30
#include <QProgressDialog>
Jasem Mutlaq's avatar
Jasem Mutlaq committed
31

32
#include <KMessageBox>
33
#include <KLocalizedString>
34

Jasem Mutlaq's avatar
Jasem Mutlaq committed
35 36 37 38 39 40
#ifdef HAVE_WCSLIB
#include <wcshdr.h>
#include <wcsfix.h>
#include <wcs.h>
#endif

41
#include "ksutils.h"
42
#include "Options.h"
43

44
#define ZOOM_DEFAULT	100.0
45 46 47 48
#define ZOOM_MIN	10
#define ZOOM_MAX	400
#define ZOOM_LOW_INCR	10
#define ZOOM_HIGH_INCR	50
49

50
const int MINIMUM_ROWS_PER_CENTER=3;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
51

52 53
#define DIFFUSE_THRESHOLD       0.15

54
#define MAX_EDGE_LIMIT  10000
Jasem Mutlaq's avatar
Jasem Mutlaq committed
55 56
#define LOW_EDGE_CUTOFF_1   50
#define LOW_EDGE_CUTOFF_2   10
57 58
#define MINIMUM_EDGE_LIMIT  2
#define SMALL_SCALE_SQUARE  256
Jasem Mutlaq's avatar
Jasem Mutlaq committed
59

60 61
bool greaterThan(Edge *s1, Edge *s2)
{
62 63
    //return s1->width > s2->width;
    return s1->sum > s2->sum;
64 65
}

66
FITSData::FITSData(FITSMode fitsMode)
Jasem Mutlaq's avatar
Jasem Mutlaq committed
67
{
68
    channels = 0;
69 70
    image_buffer = NULL;
    bayer_buffer = NULL;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
71
    wcs_coord    = NULL;
72
    fptr = NULL;
73
    maxHFRStar = NULL;
74
    darkFrame = NULL;
75
    tempFile  = false;
76
    starsSearched = false;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
77
    HasWCS = false;
78
    HasDebayer=false;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
79
    mode = fitsMode;
80 81 82 83

    debayerParams.method  = DC1394_BAYER_METHOD_NEAREST;
    debayerParams.filter  = DC1394_COLOR_FILTER_RGGB;
    debayerParams.offsetX  = debayerParams.offsetY = 0;
84 85
}

86
FITSData::~FITSData()
87
{
88 89
    int status=0;

90 91
    clearImageBuffers();

Jasem Mutlaq's avatar
Jasem Mutlaq committed
92 93
    if (starCenters.count() > 0)
        qDeleteAll(starCenters);
94

95
    delete[] wcs_coord;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
96

97
    if (fptr)
98
    {
99
        fits_close_file(fptr, &status);
100 101 102 103 104

        if (tempFile)
             QFile::remove(filename);

    }
105
}
106

107
bool FITSData::loadFITS (const QString &inFilename, bool silent)
108
{
109
    int status=0, anynull=0;
110
    long naxes[3];
111
    char error_status[512];
112
    QString errMessage;
113

Jasem Mutlaq's avatar
Jasem Mutlaq committed
114 115 116
    qDeleteAll(starCenters);
    starCenters.clear();

117
    if (fptr)
118 119
    {

120 121
        fits_close_file(fptr, &status);

122 123 124 125 126 127 128 129 130 131 132 133 134 135
        if (tempFile)
            QFile::remove(filename);
    }

    filename = inFilename;

    if (filename.contains("/tmp/"))
        tempFile = true;
    else
        tempFile = false;

    filename.remove("file://");


Jasem Mutlaq's avatar
Jasem Mutlaq committed
136
    if (fits_open_image(&fptr, filename.toLatin1(), READONLY, &status))
137 138 139
    {
        fits_report_error(stderr, status);
        fits_get_errstatus(status, error_status);
140
        errMessage = i18n("Could not open file %1. Error %2", filename, QString::fromUtf8(error_status));
141
        if (silent == false)
142
            KMessageBox::error(0, errMessage, i18n("FITS Open"));
143
        if (Options::fITSLogging())
144
            qDebug() << errMessage;
145
        return false;
146 147
    }

148
    if (fits_get_img_param(fptr, 3, &(stats.bitpix), &(stats.ndim), naxes, &status))
149 150
    {
        fits_report_error(stderr, status);
Jasem Mutlaq's avatar
Jasem Mutlaq committed
151
        fits_get_errstatus(status, error_status);
152
        errMessage = i18n("FITS file open error (fits_get_img_param): %1", QString::fromUtf8(error_status));
153
        if (silent == false)
154
            KMessageBox::error(0, errMessage, i18n("FITS Open"));
155
        if (Options::fITSLogging())
156
            qDebug() << errMessage;
157
        return false;
158 159
    }

160 161
    if (stats.ndim < 2)
    {
162
        errMessage = i18n("1D FITS images are not supported in KStars.");
163
        if (silent == false)
164
            KMessageBox::error(0, errMessage, i18n("FITS Open"));
165
        if (Options::fITSLogging())
166
            qDebug() << errMessage;
167 168 169
        return false;
    }

170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
    switch (stats.bitpix)
    {
        case 8:
            data_type = TBYTE;
            break;
        case 16:
            data_type = TUSHORT;
            break;
        case 32:
            data_type = TINT;
            break;
        case -32:
             data_type = TFLOAT;
             break;
        case 64:
             data_type = TLONGLONG;
             break;
        case -64:
             data_type = TDOUBLE;
        default:
190
            errMessage = i18n("Bit depth %1 is not supported.", stats.bitpix);
191
            if (silent == false)
192
                KMessageBox::error(NULL, errMessage, i18n("FITS Open"));
193
            if (Options::fITSLogging())
194
                qDebug() << errMessage;
195 196 197 198
            return false;
            break;
    }

199 200 201
    if (stats.ndim < 3)
        naxes[2] = 1;

202 203
    if (naxes[0] == 0 || naxes[1] == 0)
    {
204
        errMessage = i18n("Image has invalid dimensions %1x%2", naxes[0], naxes[1]);
205
        if (silent == false)
206
            KMessageBox::error(0, errMessage, i18n("FITS Open"));
207
        if (Options::fITSLogging())
208
            qDebug() << errMessage;
209 210 211
        return false;
    }

212 213
    stats.width  = naxes[0];
    stats.height = naxes[1];
214
    stats.samples_per_channel   = stats.width*stats.height;
215

216 217 218
    clearImageBuffers();

    channels = naxes[2];
219

220
    image_buffer = new float[stats.samples_per_channel * channels];
221 222
    if (image_buffer == NULL)
    {
223
       qDebug() << "FITSData: Not enough memory for image_buffer channel. Requested: " << stats.samples_per_channel * channels * sizeof(float) << " bytes.";
224 225
       clearImageBuffers();
       return false;
226
    }
227

228 229 230
    rotCounter=0;
    flipHCounter=0;
    flipVCounter=0;
231
    long nelements = stats.samples_per_channel * channels;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
232

233 234 235 236
    if (fits_read_img(fptr, TFLOAT, 1, nelements, 0, image_buffer, &anynull, &status))
    {
        char errmsg[512];
        fits_get_errstatus(status, errmsg);
237
        errMessage = i18n("Error reading image: %1", QString(errmsg));
238
        if (silent == false)
239
            KMessageBox::error(NULL, errMessage, i18n("FITS Open"));
240
        fits_report_error(stderr, status);
241
        if (Options::fITSLogging())
242
            qDebug() << errMessage;
243 244 245
        return false;
    }

246 247 248 249 250
    if (darkFrame != NULL)
        subtract(darkFrame);

    if (darkFrame == NULL)
        calculateStats();
Jasem Mutlaq's avatar
Jasem Mutlaq committed
251

Jasem Mutlaq's avatar
Jasem Mutlaq committed
252
    if (mode == FITS_NORMAL)
Jasem Mutlaq's avatar
Jasem Mutlaq committed
253
        checkWCS();
254

255 256 257
    if (checkDebayer())
        debayer();

258 259
    starsSearched = false;

260
    return true;
261

262 263
}

264
int FITSData::saveFITS( const QString &newFilename )
265
{
266
    int status=0, exttype=0;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
267
    long nelements;
268 269
    fitsfile *new_fptr;

270
    nelements = stats.samples_per_channel * channels;
271 272

    /* Create a new File, overwriting existing*/
Jasem Mutlaq's avatar
Jasem Mutlaq committed
273
    if (fits_create_file(&new_fptr, newFilename.toLatin1(), &status))
274 275
    {
        fits_report_error(stderr, status);
276
        return status;
277 278
    }

279 280 281 282 283 284
    if (fits_movabs_hdu(fptr, 1, &exttype, &status))
    {
        fits_report_error(stderr, status);
        return status;
    }

285 286 287
    if (fits_copy_file(fptr, new_fptr, 1, 1, 1, &status))
    {
        fits_report_error(stderr, status);
288
        return status;
289 290
    }

291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311
    if (HasDebayer)
    {
        if (fits_close_file(fptr, &status))
        {
            fits_report_error(stderr, status);
            return status;
        }

        if (tempFile)
        {
            QFile::remove(filename);
            tempFile = false;
        }

        filename = newFilename;

        fptr = new_fptr;

        return 0;
    }

312 313 314 315
    /* close current file */
    if (fits_close_file(fptr, &status))
    {
        fits_report_error(stderr, status);
316
        return status;
317
    }    
318

319 320
    status=0;

321 322 323 324 325 326 327 328
    if (tempFile)
    {
        QFile::remove(filename);
        tempFile = false;
    }

    filename = newFilename;

329 330
    fptr = new_fptr;

331 332 333 334 335 336
    if (fits_movabs_hdu(fptr, 1, &exttype, &status))
    {
        fits_report_error(stderr, status);
        return status;
    }

337
    /* Write Data */
Jasem Mutlaq's avatar
Jasem Mutlaq committed
338
    if (fits_write_img(fptr, TFLOAT, 1, nelements, image_buffer, &status))
339 340
    {
        fits_report_error(stderr, status);
341
        return status;
342 343 344 345 346 347 348 349
    }

    /* Write keywords */

    // Minimum
    if (fits_update_key(fptr, TDOUBLE, "DATAMIN", &(stats.min), "Minimum value", &status))
    {
        fits_report_error(stderr, status);
350
        return status;
351 352 353 354 355 356
    }

    // Maximum
    if (fits_update_key(fptr, TDOUBLE, "DATAMAX", &(stats.max), "Maximum value", &status))
    {
        fits_report_error(stderr, status);
357
        return status;
358 359
    }

360
    // NAXIS1
361
    if (fits_update_key(fptr, TUSHORT, "NAXIS1", &(stats.width), "length of data axis 1", &status))
362 363 364 365 366 367
    {
        fits_report_error(stderr, status);
        return status;
    }

    // NAXIS2
368
    if (fits_update_key(fptr, TUSHORT, "NAXIS2", &(stats.height), "length of data axis 2", &status))
369 370 371 372 373
    {
        fits_report_error(stderr, status);
        return status;
    }

374 375 376 377
    // ISO Date
    if (fits_write_date(fptr, &status))
    {
        fits_report_error(stderr, status);
378
        return status;
379 380
    }

381
    QString history = QString("Modified by KStars on %1").arg(QDateTime::currentDateTime().toString("yyyy-MM-ddThh:mm:ss"));
382
    // History
Jasem Mutlaq's avatar
Jasem Mutlaq committed
383
    if (fits_write_history(fptr, history.toLatin1(), &status))
384 385
    {
        fits_report_error(stderr, status);
386
        return status;
387 388
    }

389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
    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 || mirror)
        rotWCSFITS(rot, mirror);

    rotCounter=flipHCounter=flipVCounter=0;

404
    return status;
405 406
}

407 408
void FITSData::clearImageBuffers()
{
409
        delete[] image_buffer;
410
        image_buffer=NULL;
411
        delete [] bayer_buffer;
412
        bayer_buffer=NULL;
413
}
Jasem Mutlaq's avatar
Jasem Mutlaq committed
414

415
int FITSData::calculateMinMax(bool refresh)
416
{
417 418 419 420
    int status, nfound=0;

    status = 0;

Jasem Mutlaq's avatar
Jasem Mutlaq committed
421
    if (refresh == false)
422
    {
423
        if (fits_read_key_dbl(fptr, "DATAMIN", &(stats.min[0]), NULL, &status) ==0)
424 425
            nfound++;

426
        if (fits_read_key_dbl(fptr, "DATAMAX", &(stats.max[0]), NULL, &status) == 0)
427 428
            nfound++;

429
        // If we found both keywords, no need to calculate them, unless they are both zeros
430
        if (nfound == 2 && !(stats.min[0] == 0 && stats.max[0] ==0))
431 432
            return 0;
    }
433

434 435
    stats.min[0]= 1.0E30;
    stats.max[0]= -1.0E30;
436

437 438 439 440 441 442 443
    stats.min[1]= 1.0E30;
    stats.max[1]= -1.0E30;

    stats.min[2]= 1.0E30;
    stats.max[2]= -1.0E30;

    if (channels == 1)
444
    {
445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472
        for (unsigned int i=0; i < stats.samples_per_channel; i++)
        {
            if (image_buffer[i] < stats.min[0]) stats.min[0] = image_buffer[i];
            else if (image_buffer[i] > stats.max[0]) stats.max[0] = image_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 (image_buffer[i] < stats.min[0])
                stats.min[0] = image_buffer[i];
            else if (image_buffer[i] > stats.max[0])
                stats.max[0] = image_buffer[i];

            if (image_buffer[i+g_offset] < stats.min[1])
                stats.min[1] = image_buffer[i+g_offset];
            else if (image_buffer[i+g_offset] > stats.max[1])
                stats.max[1] = image_buffer[i+g_offset];

            if (image_buffer[i+b_offset] < stats.min[2])
                stats.min[2] = image_buffer[i+b_offset];
            else if (image_buffer[i+b_offset] > stats.max[2])
                stats.max[2] = image_buffer[i+b_offset];
        }
473
    }
474

Jasem Mutlaq's avatar
Jasem Mutlaq committed
475
    //qDebug() << "DATAMIN: " << stats.min << " - DATAMAX: " << stats.max;
476
    return 0;
477 478
}

479

480
void FITSData::calculateStats(bool refresh)
481
{
482
    // Calculate min max
Jasem Mutlaq's avatar
Jasem Mutlaq committed
483
    calculateMinMax(refresh);
484

485
    // Get standard deviation and mean in one run
486
    runningAverageStdDev();
487

488
    stats.SNR = stats.mean[0] / stats.stddev[0];
489

490 491 492 493
    if (refresh && markStars)
        // Let's try to find star positions again after transformation
        starsSearched = false;

494 495
}

496 497 498 499 500 501
void FITSData::runningAverageStdDev()
{
    int m_n = 2;
    double m_oldM=0, m_newM=0, m_oldS=0, m_newS=0;
    m_oldM = m_newM = image_buffer[0];

502
    for (unsigned int i=1; i < stats.samples_per_channel; i++)
503 504 505 506 507 508 509 510 511 512 513
    {
        m_newM = m_oldM + (image_buffer[i] - m_oldM)/m_n;
        m_newS = m_oldS + (image_buffer[i] - m_oldM) * (image_buffer[i] - m_newM);

        m_oldM = m_newM;
        m_oldS = m_newS;
        m_n++;
    }

    double variance = m_newS / (m_n -2);

514 515
    stats.mean[0] = m_newM;
    stats.stddev[0]  = sqrt(variance);
516 517
}

518
void FITSData::setMinMax(double newMin,  double newMax, uint8_t channel)
519
{
520 521
    stats.min[channel] = newMin;
    stats.max[channel] = newMax;
522 523
}

524
int FITSData::getFITSRecord(QString &recordList, int &nkeys)
525
{
526 527 528 529 530 531 532 533 534 535 536 537 538 539 540
    char *header;
    int status=0;

    if (fits_hdr2str(fptr, 0, NULL, 0, &header, &nkeys, &status))
    {
        fits_report_error(stderr, status);
        free(header);
        return -1;
    }

    recordList = QString(header);

    free(header);

    return 0;
541
}
542

543
bool FITSData::checkCollision(Edge* s1, Edge*s2)
544 545 546 547 548 549
{
    int dis; //distance

    int diff_x=s1->x - s2->x;
    int diff_y=s1->y - s2->y;

550
    dis = std::abs( sqrt( diff_x*diff_x + diff_y*diff_y));
551 552 553 554 555 556 557 558 559 560
    dis -= s1->width/2;
    dis -= s2->width/2;

    if (dis<=0) //collision
    return true;

    //no collision
    return false;
}

Jasem Mutlaq's avatar
Jasem Mutlaq committed
561 562

/*** Find center of stars and calculate Half Flux Radius */
563
void FITSData::findCentroid(int initStdDev, int minEdgeWidth)
Jasem Mutlaq's avatar
Jasem Mutlaq committed
564
{
565
    double threshold=0;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
566 567
    double avg = 0;
    double sum=0;
568
    double min=0;
569 570
    double noiseAvg=0,noiseSum=0;
    int pixelRadius =0, noisePixelRadius=0;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
571
    int pixVal=0;
572 573
    int noisePix=0;
    int minimumEdgeCount = MINIMUM_EDGE_LIMIT;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
574

575
    int noisePixLimit=0;
576
    double JMIndex = histogram->getJMIndex();
577

Jasem Mutlaq's avatar
Jasem Mutlaq committed
578 579
    QList<Edge*> edges;

580
    if (JMIndex < DIFFUSE_THRESHOLD)
581 582
    {
            minEdgeWidth = JMIndex*35+1;
583
            minimumEdgeCount=minEdgeWidth-1;
584 585 586 587
    }
    else
    {
            minEdgeWidth =6;
588
            minimumEdgeCount=4;
589 590
    }

591 592
    while (initStdDev >= 1)
    {
593
        minEdgeWidth--;
594
        minimumEdgeCount--;
595 596 597

        if (minEdgeWidth < 3)
            minEdgeWidth = 3;
598 599
        if (minimumEdgeCount < 3)
            minimumEdgeCount=3;
600

601
       if (JMIndex < DIFFUSE_THRESHOLD)
602
       {
603 604 605
           //threshold = stats.max[0] - stats.stddev[0]* (MINIMUM_STDVAR - initStdDev +1);
           // 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);
606

607
           min =stats.min[0];
608

609
           noisePixLimit=minEdgeWidth*0.5;
610 611 612
       }
       else
       {
613 614 615
           //threshold = (stats.max[0] - stats.min[0])/2.0 + stats.min[0]  + stats.stddev[0]* (MINIMUM_STDVAR - initStdDev);
           //if ( (stats.max[0] - stats.min[0])/2.0 > (stats.mean[0]+stats.stddev[0]*5))
           threshold = stats.mean[0]+stats.stddev[0]*initStdDev*(0.95 - (MINIMUM_STDVAR - initStdDev) * 0.05);
616
           min = stats.min[0];
617
           noisePixLimit =2;
618

619
       }
620

621
        if (Options::fITSLogging())
622 623 624 625
        {
            qDebug() << "SNR: " << stats.SNR;
            qDebug() << "The threshold level is " << threshold << " minimum edge width" << minEdgeWidth << " minimum edge limit " << minimumEdgeCount;
        }
626

627
        threshold -= stats.min[0];
628

629 630
    int subX, subY, subW, subH;

Jasem Mutlaq's avatar
Jasem Mutlaq committed
631 632 633 634 635 636 637 638 639
    if (mode == FITS_GUIDE)
    {
        subX = stats.width/10;
        subY = stats.height/10;
        subW = stats.width - subX;
        subH = stats.height - subY;
    }
    else
    {
Jasem Mutlaq's avatar
Jasem Mutlaq committed
640 641
        subX = 0;
        subY = 0;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
642 643 644
        subW = stats.width;
        subH = stats.height;
    }
645

Jasem Mutlaq's avatar
Jasem Mutlaq committed
646
    // Detect "edges" that are above threshold
647
    for (int i=subY; i < subH; i++)
Jasem Mutlaq's avatar
Jasem Mutlaq committed
648 649 650
    {
        pixelRadius = 0;

651
        for(int j=subX; j < subW; j++)
Jasem Mutlaq's avatar
Jasem Mutlaq committed
652
        {
653
            pixVal = image_buffer[j+(i*stats.width)] - min;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
654 655

            // If pixel value > threshold, let's get its weighted average
656
            if ( pixVal >= threshold || (sum > 0 && noisePix <= noisePixLimit))
Jasem Mutlaq's avatar
Jasem Mutlaq committed
657
            {
658
                if (pixVal < threshold)
659 660 661 662 663 664 665 666 667 668 669 670 671 672
                {
                    noisePix++;
                    noiseAvg += j * pixVal;
                    noiseSum += pixVal;
                    noisePixelRadius++;
                    continue;
                }
                else if (noisePix)
                {
                   avg += noiseAvg;
                   sum += noiseSum;
                   pixelRadius += noisePixelRadius;
                   noisePix=noiseAvg=noiseSum=noisePixelRadius=0;
                }
673 674 675 676

                avg += j * pixVal;
                sum += pixVal;
                pixelRadius++;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
677
            }
Yuri Chornoivan's avatar
Yuri Chornoivan committed
678
            // Value < threshold but avg exists
Jasem Mutlaq's avatar
Jasem Mutlaq committed
679 680
            else if (sum > 0)
            {
681

Jasem Mutlaq's avatar
Jasem Mutlaq committed
682
                // We found a potential centroid edge
683
                if (pixelRadius >= (minEdgeWidth - (MINIMUM_STDVAR - initStdDev)))
Jasem Mutlaq's avatar
Jasem Mutlaq committed
684
                {
685
                    float center = avg/sum + 0.5;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
686 687
                    if (center > 0)
                    {
688
                        int i_center = floor(center);
Jasem Mutlaq's avatar
Jasem Mutlaq committed
689

Jasem Mutlaq's avatar
Jasem Mutlaq committed
690
                        Edge *newEdge = new Edge();
Jasem Mutlaq's avatar
Jasem Mutlaq committed
691

Jasem Mutlaq's avatar
Jasem Mutlaq committed
692
                        newEdge->x          = center;
693
                        newEdge->y          = i + 0.5;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
694
                        newEdge->scanned    = 0;
695
                        newEdge->val        = image_buffer[i_center+(i*stats.width)] - min;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
696 697 698 699 700 701
                        newEdge->width      = pixelRadius;
                        newEdge->HFR        = 0;
                        newEdge->sum        = sum;

                        edges.append(newEdge);
                    }
Jasem Mutlaq's avatar
Jasem Mutlaq committed
702 703 704 705

                }

                // Reset
706 707
                avg= sum = pixelRadius=0;
                noisePix = noiseAvg = noiseSum = noisePixelRadius = 0;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
708 709 710 711
            }
         }
     }

712
    if (Options::fITSLogging())
713
        qDebug() << "Total number of edges found is: " << edges.count();
Jasem Mutlaq's avatar
Jasem Mutlaq committed
714

715 716 717 718 719 720 721
    // In case of hot pixels
    if (edges.count() == 1 && initStdDev > 1)
    {
        initStdDev--;
        continue;
    }

722 723
    if (edges.count() >= MAX_EDGE_LIMIT)
    {
724
        if (Options::fITSLogging())
725 726
            qDebug() << "Too many edges, aborting... " << edges.count();

727 728 729 730
        qDeleteAll(edges);
        return;
    }

731
    if (edges.count() >= minimumEdgeCount)
732 733 734 735 736 737 738
        break;

      qDeleteAll(edges);
      edges.clear();
      initStdDev--;
    }

Jasem Mutlaq's avatar
Jasem Mutlaq committed
739 740 741 742
    int cen_count=0;
    int cen_x=0;
    int cen_y=0;
    int cen_v=0;
743 744
    int cen_w=0;
    int width_sum=0;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
745

746
    // Let's sort edges, starting with widest
747
    qSort(edges.begin(), edges.end(), greaterThan);
748

Jasem Mutlaq's avatar
Jasem Mutlaq committed
749 750 751
    // Now, let's scan the edges and find the maximum centroid vertically
    for (int i=0; i < edges.count(); i++)
    {
752
        if (Options::fITSLogging())
753 754
            qDebug() << "# " << 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;
755

Jasem Mutlaq's avatar
Jasem Mutlaq committed
756 757 758
        // If edge scanned already, skip
        if (edges[i]->scanned == 1)
        {
759
            if (Options::fITSLogging())
760
                qDebug() << "Skipping check for center " << i << " because it was already counted";
Jasem Mutlaq's avatar
Jasem Mutlaq committed
761 762 763
            continue;
        }

764
        if (Options::fITSLogging())
765
            qDebug() << "Invetigating edge # " << i << " now ...";
766

Jasem Mutlaq's avatar
Jasem Mutlaq committed
767 768 769
        // Get X, Y, and Val of edge
        cen_x = edges[i]->x;
        cen_y = edges[i]->y;
770
        cen_v = edges[i]->sum;
771
        cen_w = edges[i]->width;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
772

773 774 775
        float avg_x = 0;
        float avg_y = 0;

776
        sum = 0;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
777 778 779
        cen_count=0;

        // Now let's compare to other edges until we hit a maxima
780
        for (int j=0; j < edges.count();j++)
Jasem Mutlaq's avatar
Jasem Mutlaq committed
781 782 783 784
        {
            if (edges[j]->scanned)
                continue;

785 786 787 788 789
            if (checkCollision(edges[j], edges[i]))
            {
                if (edges[j]->sum >= cen_v)
                {
                    cen_v = edges[j]->sum;
790
                    cen_w = edges[j]->width;
791 792 793 794 795
                }

                edges[j]->scanned = 1;
                cen_count++;

796 797 798 799
                avg_x += edges[j]->x * edges[j]->val;
                avg_y += edges[j]->y * edges[j]->val;
                sum += edges[j]->val;

800 801 802
                continue;
            }

Jasem Mutlaq's avatar
Jasem Mutlaq committed
803 804
        }

805 806
        int cen_limit = (MINIMUM_ROWS_PER_CENTER - (MINIMUM_STDVAR - initStdDev));

Jasem Mutlaq's avatar
Jasem Mutlaq committed
807 808 809 810 811 812 813
        if (edges.count() < LOW_EDGE_CUTOFF_1)
        {
            if (edges.count() < LOW_EDGE_CUTOFF_2)
                cen_limit = 1;
            else
                cen_limit = 2;
        }
814

815
    if (Options::fITSLogging())
816
        qDebug() << "center_count: " << cen_count << " and initstdDev= " << initStdDev << " and limit is " << cen_limit;
817

818
        if (cen_limit < 1)
819 820
            continue;

Jasem Mutlaq's avatar
Jasem Mutlaq committed
821
        // If centroid count is within acceptable range
822
        //if (cen_limit >= 2 && cen_count >= cen_limit)
823
        if (cen_count >= cen_limit)
Jasem Mutlaq's avatar
Jasem Mutlaq committed
824 825 826 827
        {
            // We detected a centroid, let's init it
            Edge *rCenter = new Edge();

828 829
            rCenter->x = avg_x/sum;
            rCenter->y = avg_y/sum;
830 831
            width_sum += rCenter->width;
            rCenter->width = cen_w;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
832

833
           if (Options::fITSLogging())
834
            qDebug() << "Found a real center with number with (" << rCenter->x << "," << rCenter->y << ")";
Jasem Mutlaq's avatar
Jasem Mutlaq committed
835 836 837 838 839 840

            // Calculate Total Flux From Center, Half Flux, Full Summation
            double TF=0;
            double HF=0;
            double FSum=0;

841 842
            cen_x = (int) floor(rCenter->x);
            cen_y = (int) floor(rCenter->y);
843

844
            if (cen_x < 0 || cen_x > stats.width || cen_y < 0 || cen_y > stats.height)
845 846 847
                continue;


848 849 850
            // Complete sum along the radius
            //for (int k=0; k < rCenter->width; k++)
            for (int k=rCenter->width/2; k >= -(rCenter->width/2) ; k--)
851
                FSum += image_buffer[cen_x-k+(cen_y*stats.width)] - min;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
852 853 854 855 856

            // Half flux
            HF = FSum / 2.0;

            // Total flux starting from center
857
            TF = image_buffer[cen_y * stats.width + cen_x] - min;
Jasem Mutlaq's avatar
Jasem Mutlaq committed
858 859 860 861 862 863 864 865

            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)
                {
866
                    if (Options::fITSLogging())
867
                        qDebug() << "Stopping at TF " << TF << " after #" << k << " pixels.";
Jasem Mutlaq's avatar
Jasem Mutlaq committed
868 869 870
                    break;
                }

871 872 873
                TF += image_buffer[cen_y * stats.width + cen_x + k] - min;
                TF += image_buffer[cen_y * stats.width + cen_x - k] - min;

Jasem Mutlaq's avatar
Jasem Mutlaq committed
874 875 876 877 878 879 880 881
                pixelCounter++;
            }

            // Calculate weighted Half Flux Radius
            rCenter->HFR = pixelCounter * (HF / TF);
            // Store full flux
            rCenter->val = FSum;

882
            if (Options::fITSLogging())
883 884 885
                qDebug() << "HFR for this center is " << rCenter->HFR << " pixels and the total flux is " << FSum;

            starCenters.append(rCenter);
886

Jasem Mutlaq's avatar
Jasem Mutlaq committed
887 888 889
        }
    }

890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910
    if (starCenters.count() > 1 && mode != FITS_FOCUS)
    {
        float width_avg = width_sum / starCenters.count();

        float lsum =0, sdev=0;

        foreach(Edge *center, starCenters)
            lsum += (center->width - width_avg) * (center->width - width_avg);

        sdev = (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;

    }

Jasem Mutlaq's avatar
Jasem Mutlaq committed
911 912
    // Release memory
    qDeleteAll(edges);
913

Jasem Mutlaq's avatar
Jasem Mutlaq committed
914 915
}

916
double FITSData::getHFR(HFRType type)
Jasem Mutlaq's avatar
Jasem Mutlaq committed
917
{