20 #include <config-kstars.h>
27 #include <QApplication>
29 #include <QProgressDialog>
30 #include <KMessageBox>
36 #include <getwcstab.h>
41 #define ZOOM_DEFAULT 100.0
44 #define ZOOM_LOW_INCR 10
45 #define ZOOM_HIGH_INCR 50
49 #define JM_LOWER_LIMIT 5
50 #define JM_UPPER_LIMIT 400
52 #define LOW_EDGE_CUTOFF_1 50
53 #define LOW_EDGE_CUTOFF_2 10
70 starsSearched =
false;
81 if (starCenters.count() > 0)
82 qDeleteAll(starCenters);
88 fits_close_file(fptr, &status);
91 QFile::remove(filename);
98 int status=0, nulval=0, anynull=0;
99 long fpixel[2], nelements, naxes[2];
100 char error_status[512];
102 qDeleteAll(starCenters);
107 progress->setLabelText(i18n(
"Please hold while loading FITS file..."));
108 progress->setWindowTitle(i18n(
"Loading FITS"));
112 progress->setValue(30);
117 fits_close_file(fptr, &status);
120 QFile::remove(filename);
123 filename = inFilename;
125 if (filename.contains(
"/tmp/"))
130 filename.remove(
"file://");
133 if (fits_open_image(&fptr, filename.toAscii(), READONLY, &status))
135 fits_report_error(stderr, status);
136 fits_get_errstatus(status, error_status);
138 KMessageBox::error(0, i18n(
"Could not open file %1 (fits_get_img_param). Error %2", filename, QString::fromUtf8(error_status)), i18n(
"FITS Open"));
143 if (progress->wasCanceled())
147 progress->setValue(40);
150 if (fits_get_img_param(fptr, 2, &(
stats.bitpix), &(
stats.ndim), naxes, &status))
152 fits_report_error(stderr, status);
153 fits_get_errstatus(status, error_status);
156 KMessageBox::error(0, i18n(
"FITS file open error (fits_get_img_param): %1", QString::fromUtf8(error_status)), i18n(
"FITS Open"));
163 KMessageBox::error(0, i18n(
"1D FITS images are not supported in KStars."), i18n(
"FITS Open"));
168 if (fits_get_img_type(fptr, &data_type, &status))
170 fits_report_error(stderr, status);
171 fits_get_errstatus(status, error_status);
174 KMessageBox::error(0, i18n(
"FITS file open error (fits_get_img_type): %1", QString::fromUtf8(error_status)), i18n(
"FITS Open"));
179 if (progress->wasCanceled())
183 progress->setValue(60);
185 stats.dim[0] = naxes[0];
186 stats.dim[1] = naxes[1];
188 delete (image_buffer);
191 image_buffer =
new float[
stats.dim[0] *
stats.dim[1]];
193 if (image_buffer == NULL)
195 qDebug() <<
"Not enough memory for image_buffer";
200 if (progress->wasCanceled())
202 delete (image_buffer);
208 progress->setValue(70);
214 qApp->processEvents();
216 if (fits_read_2d_flt(fptr, 0, nulval, naxes[0], naxes[0], naxes[1], image_buffer, &anynull, &status))
218 fprintf(stderr,
"fits_read_pix error\n");
219 fits_report_error(stderr, status);
225 if (progress->wasCanceled())
227 delete (image_buffer);
235 progress->setValue(80);
240 qApp->processEvents();
247 progress->setValue(90);
252 if (progress->wasCanceled())
254 delete (image_buffer);
260 progress->setValue(100);
262 starsSearched =
false;
271 long fpixel[2], nelements;
280 if (fits_create_file(&new_fptr, newFilename.toAscii(), &status))
282 fits_report_error(stderr, status);
287 if (fits_copy_file(fptr, new_fptr, 1, 1, 1, &status))
289 fits_report_error(stderr, status);
294 if (fits_close_file(fptr, &status))
296 fits_report_error(stderr, status);
302 QFile::remove(filename);
306 filename = newFilename;
311 if (fits_write_pix(fptr, TFLOAT, fpixel, nelements, image_buffer, &status))
313 fits_report_error(stderr, status);
320 if (fits_update_key(fptr, TDOUBLE,
"DATAMIN", &(
stats.min),
"Minimum value", &status))
322 fits_report_error(stderr, status);
327 if (fits_update_key(fptr, TDOUBLE,
"DATAMAX", &(
stats.max),
"Maximum value", &status))
329 fits_report_error(stderr, status);
334 if (fits_write_date(fptr, &status))
336 fits_report_error(stderr, status);
340 QString history = QString(
"Modified by KStars on %1").arg(QDateTime::currentDateTime().toString(
"yyyy-MM-ddThh:mm:ss"));
342 if (fits_write_history(fptr, history.toAscii(), &status))
344 fits_report_error(stderr, status);
352 int FITSImage::calculateMinMax(
bool refresh)
354 int status, nfound=0;
359 if (refresh ==
false)
361 if (fits_read_key_dbl(fptr,
"DATAMIN", &(
stats.min), NULL, &status) ==0)
364 if (fits_read_key_dbl(fptr,
"DATAMAX", &(
stats.max), NULL, &status) == 0)
376 for (
int i=0; i < npixels; i++)
378 if (image_buffer[i] <
stats.min)
stats.min = image_buffer[i];
379 else if (image_buffer[i] >
stats.max)
stats.max = image_buffer[i];
389 calculateMinMax(refresh);
395 if (refresh && markStars)
397 starsSearched =
false;
401 double FITSImage::average()
405 int width =
stats.dim[0];
406 int height =
stats.dim[1];
408 if (!image_buffer)
return -1;
410 for (
int i= 0 ; i < height; i++)
413 for (
int j= 0; j < width; j++)
414 sum += image_buffer[row+j];
417 return (sum / (width * height ));
426 int width =
stats.dim[0];
427 int height =
stats.dim[1];
429 if (!image_buffer)
return -1;
431 for (
int i= 0 ; i < height; i++)
434 for (
int j= 0; j < width; j++)
436 lsum += (image_buffer[row + j] -
stats.average) * (image_buffer[row+j] -
stats.average);
440 return (sqrt(lsum/(width * height - 1)));
456 if (fits_hdr2str(fptr, 0, NULL, 0, &header, &nkeys, &status))
458 fits_report_error(stderr, status);
463 recordList = QString(header);
470 bool FITSImage::checkCollision(
Edge* s1,
Edge*s2)
474 int diff_x=s1->
x - s2->
x;
475 int diff_y=s1->
y - s2->
y;
477 dis = abs( sqrt( diff_x*diff_x + diff_y*diff_y));
501 while (initStdDev >= 1)
503 threshold =
stats.stddev* initStdDev;
506 qDebug() <<
"The threshold level is " << threshold << endl;
510 for (
int i=0; i <
stats.dim[1]; i++)
514 for(
int j=0; j <
stats.dim[0]; j++)
516 pixVal = image_buffer[j+(i*
stats.dim[0])] -
stats.min;
520 if ( pixVal >= threshold || (sum > 0 && badPix <= 2))
522 if (pixVal < threshold)
537 if (pixelRadius >= (minEdgeWidth - (
MINIMUM_STDVAR - initStdDev)))
541 float center = avg/sum;
546 int i_center = round(center);
553 newEdge->
val = image_buffer[i_center+(i*
stats.dim[0])] -
stats.min;
554 newEdge->
width = pixelRadius;
558 edges.append(newEdge);
577 qDebug() <<
"Total number of edges found is: " << edges.count() << endl;
599 for (
int i=0; i < edges.count(); i++)
602 qDebug() <<
"# " << i <<
" Edge at (" << edges[i]->x <<
"," << edges[i]->y <<
") With a value of " << edges[i]->val <<
" and width of "
603 << edges[i]->width <<
" pixels. with sum " << edges[i]->sum << endl;
607 if (edges[i]->scanned == 1)
610 qDebug() <<
"Skipping check for center " << i <<
" because it was already counted" << endl;
616 qDebug() <<
"Invetigating edge # " << i <<
" now ..." << endl;
622 cen_v = edges[i]->sum;
623 cen_w = edges[i]->width;
632 for (
int j=0; j < edges.count();j++)
634 if (edges[j]->scanned)
637 if (checkCollision(edges[j], edges[i]))
639 if (edges[j]->sum >= cen_v)
641 cen_v = edges[j]->sum;
642 cen_w = edges[j]->width;
645 edges[j]->scanned = 1;
648 avg_x += edges[j]->x * edges[j]->val;
649 avg_y += edges[j]->y * edges[j]->val;
650 sum += edges[j]->val;
668 qDebug() <<
"center_count: " << cen_count <<
" and initstdDev= " << initStdDev <<
" and limit is "
669 << cen_limit << endl;
672 if (cen_limit < 1 || (cen_w > (0.2 *
stats.dim[0])))
677 if (cen_count >= cen_limit)
683 rCenter->
x = avg_x/sum;
684 rCenter->
y = avg_y/sum;
686 width_sum += rCenter->
width;
691 rCenter->
width = cen_w;
694 qDebug() <<
"Found a real center with number " << rc_index <<
"with (" << rCenter->
x <<
"," << rCenter->
y <<
")" << endl;
696 qDebug() <<
"Profile for this center is:" << endl;
697 for (
int i=edges[rc_index]->width/2; i >= -(edges[rc_index]->width/2) ; i--)
698 qDebug() <<
"#" << i <<
" , " << image_buffer[(int) round(rCenter->
x-i+(rCenter->
y*
stats.dim[0]))] -
stats.min << endl;
707 cen_x = (
int) round(rCenter->
x);
708 cen_y = (int) round(rCenter->
y);
710 if (cen_x < 0 || cen_x >
stats.dim[0] || cen_y < 0 || cen_y >
stats.dim[1])
716 for (
int k=rCenter->
width/2; k >= -(rCenter->
width/2) ; k--)
717 FSum += image_buffer[cen_x-k+(cen_y*
stats.dim[0])] -
stats.min;
723 TF = image_buffer[cen_y *
stats.dim[0] + cen_x] -
stats.min;
725 int pixelCounter = 1;
728 for (
int k=1; k < rCenter->
width/2; k++)
730 TF += image_buffer[cen_y *
stats.dim[0] + cen_x + k] -
stats.min;
731 TF += image_buffer[cen_y *
stats.dim[0] + cen_x - k] -
stats.min;
736 qDebug() <<
"Stopping at TF " << TF <<
" after #" << k <<
" pixels." << endl;
745 rCenter->
HFR = pixelCounter * (HF / TF);
750 qDebug() <<
"HFR for this center is " << rCenter->
HFR <<
" pixels and the total flux is " << FSum << endl;
752 starCenters.append(rCenter);
757 if (starCenters.count() > 1 && mode !=
FITS_FOCUS)
759 float width_avg = width_sum / starCenters.count();
761 float lsum =0, sdev=0;
763 foreach(
Edge *center, starCenters)
764 lsum += (center->
width - width_avg) * (center->
width - width_avg);
766 sdev = (sqrt(lsum/(starCenters.count() - 1))) * 4;
769 foreach(
Edge *center, starCenters)
770 if (center->
width > sdev)
771 starCenters.removeOne(center);
789 if (starCenters.size() == 0)
796 for (
int i=0; i < starCenters.count() ; i++)
798 if (starCenters[i]->val > maxVal)
801 maxVal = starCenters[i]->val;
806 return starCenters[maxIndex]->HFR;
813 for (
int i=0; i < starCenters.count() ; i++)
815 avgHFR += starCenters[i]->val * starCenters[i]->HFR;
816 FSum += starCenters[i]->val;
822 return (avgHFR / FSum);
831 if (type ==
FITS_NONE || histogram == NULL)
835 float val=0,bufferVal =0;
838 image = image_buffer;
840 int width =
stats.dim[0];
841 int height =
stats.dim[1];
852 for (
int i=0; i < height; i++)
853 for (
int j=0; j < width; j++)
855 bufferVal = image[i * width + j];
856 if (bufferVal < min) bufferVal =
min;
857 else if (bufferVal > max) bufferVal =
max;
858 image[i * width + j] = bufferVal;
864 coeff = max / log(1 + max);
866 for (
int i=0; i < height; i++)
867 for (
int j=0; j < width; j++)
869 bufferVal = image[i * width + j];
870 if (bufferVal < min) bufferVal =
min;
871 else if (bufferVal > max) bufferVal =
max;
872 val = (coeff * log(1 + bufferVal));
873 if (val < min) val =
min;
874 else if (val > max) val =
max;
875 image[i * width + j] = val;
880 coeff = max / sqrt(max);
882 for (
int i=0; i < height; i++)
883 for (
int j=0; j < width; j++)
885 bufferVal = (int) image[i * width + j];
886 if (bufferVal < min) bufferVal =
min;
887 else if (bufferVal > max) bufferVal =
max;
888 val = (int) (coeff * sqrt(bufferVal));
889 image[i * width + j] = val;
901 for (
int i=0; i < height; i++)
902 for (
int j=0; j < width; j++)
904 bufferVal = image[i * width + j];
905 if (bufferVal < min) bufferVal =
min;
906 else if (bufferVal > max) bufferVal =
max;
907 image[i * width + j] = bufferVal;
921 for (
int i=0; i < height; i++)
922 for (
int j=0; j < width; j++)
924 bufferVal = image[i * width + j];
925 if (bufferVal < min) bufferVal =
min;
926 else if (bufferVal > max) bufferVal =
max;
927 image[i * width + j] = bufferVal;
934 QVarLengthArray<int, INITIAL_MAXIMUM_WIDTH> cumulativeFreq = histogram->
getCumulativeFreq();
935 coeff = 255.0 / (height * width);
937 for (
int i=0; i < height; i++)
938 for (
int j=0; j < width; j++)
940 bufferVal = (int) (image[i * width + j] - min) * histogram->
getBinWidth();
942 if (bufferVal >= cumulativeFreq.size())
943 bufferVal = cumulativeFreq.size()-1;
945 val = (int) (coeff * cumulativeFreq[bufferVal]);
947 image[i * width + j] = val;
954 for (
int i=0; i < height; i++)
955 for (
int j=0; j < width; j++)
957 bufferVal = image[i * width + j];
958 if (bufferVal < min) bufferVal =
min;
959 else if (bufferVal > max) bufferVal =
max;
960 image[i * width + j] = bufferVal;
975 for (
int i=0; i <
stats.dim[0]*
stats.dim[1]; i++)
977 image_buffer[i] -= dark_buffer[i];
978 if (image_buffer[i] < 0)
987 if (histogram == NULL)
990 if (starsSearched ==
false)
992 qDeleteAll(starCenters);
1002 starsSearched =
true;
1004 return starCenters.count();
1010 if (starCenters.count() == 0)
1018 foreach(
Edge *center, starCenters)
1019 if (checkCollision(pEdge, center))
1029 void FITSImage::checkWCS()
1035 int nkeyrec, nreject, nwcs, stat[2];
1036 double imgcrd[2], phi, pixcrd[2], theta, world[2];
1037 struct wcsprm *wcs=0;
1041 if (fits_hdr2str(fptr, 1, NULL, 0, &header, &nkeyrec, &status))
1043 fits_report_error(stderr, status);
1047 if ((status = wcspih(header, nkeyrec, WCSHDR_all, -3, &nreject, &nwcs, &wcs)))
1049 fprintf(stderr,
"wcspih ERROR %d: %s.\n", status, wcshdr_errmsg[status]);
1062 if (wcs->crpix[0] == 0)
1067 if ((status = wcsset(wcs)))
1069 fprintf(stderr,
"wcsset ERROR %d: %s.\n", status, wcs_errmsg[status]);
1075 wcs_coord =
new wcs_point[width*height];
1079 for (
int i=0; i < height; i++)
1081 for (
int j=0; j < width; j++)
1086 if ((status = wcsp2s(wcs, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0])))
1088 fprintf(stderr,
"wcsp2s ERROR %d: %s.\n", status,
1089 wcs_errmsg[status]);
#define LOW_EDGE_CUTOFF_1
bool greaterThan(Edge *s1, Edge *s2)
int getFITSRecord(QString &recordList, int &nkeys)
void findCentroid(int initStdDev=MINIMUM_STDVAR, int minEdgeWidth=MINIMUM_PIXEL_RANGE)
int saveFITS(const QString &filename)
struct FITSImage::@1 stats
void applyFilter(FITSScale type, float *image=NULL, int min=-1, int max=-1)
void getCenterSelection(int *x, int *y)
FITSImage(FITSMode mode=FITS_NORMAL)
double getHFR(HFRType type=HFR_AVERAGE)
#define LOW_EDGE_CUTOFF_2
void subtract(float *darkFrame)
bool loadFITS(const QString &filename, QProgressDialog *progress=NULL)
void setFITSMinMax(double newMin, double newMax)
void calculateStats(bool refresh=false)
QVarLengthArray< int, INITIAL_MAXIMUM_WIDTH > getCumulativeFreq()
const int MINIMUM_ROWS_PER_CENTER