/* * makeit.c * * Main program. * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * * This file part of: SExtractor * * Copyright: (C) 1993-2013 Emmanuel Bertin -- IAP/CNRS/UPMC * * License: GNU General Public License * * SExtractor 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. * SExtractor 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 SExtractor. If not, see . * * Last modified: 13/02/2013 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ #ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include #include #include #include #include "define.h" #include "globals.h" #include "prefs.h" #include "fits/fitscat.h" #include "assoc.h" #include "back.h" #include "check.h" #include "fft.h" #include "field.h" #include "filter.h" #include "growth.h" #include "interpolate.h" #include "pattern.h" #include "psf.h" #include "profit.h" #include "som.h" #include "weight.h" #include "xml.h" static int selectext(char *filename); time_t thetimet, thetimet2; extern profitstruct *theprofit,*thedprofit, *thepprofit,*theqprofit; extern char profname[][32]; double dtime; /******************************** makeit *************************************/ /* Manage the whole stuff. */ void makeit() { checkstruct *check; picstruct *dfield, *field,*pffield[MAXFLAG], *wfield,*dwfield; catstruct *imacat; tabstruct *imatab; patternstruct *pattern; static time_t thetime1, thetime2; struct tm *tm; unsigned int modeltype; int nflag[MAXFLAG], nparam2[2], i, nok, ntab, next, ntabmax, forcextflag, nima0,nima1, nweight0,nweight1, npsf0,npsf1, npat,npat0; /* Install error logging */ error_installfunc(write_error); /* Processing start date and time */ dtime = counter_seconds(); thetimet = time(NULL); tm = localtime(&thetimet); sprintf(prefs.sdate_start,"%04d-%02d-%02d", tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday); sprintf(prefs.stime_start,"%02d:%02d:%02d", tm->tm_hour, tm->tm_min, tm->tm_sec); NFPRINTF(OUTPUT, ""); QPRINTF(OUTPUT, "----- %s %s started on %s at %s with %d thread%s\n\n", BANNER, MYVERSION, prefs.sdate_start, prefs.stime_start, prefs.nthreads, prefs.nthreads>1? "s":""); /* Initialize globals variables */ initglob(); NFPRINTF(OUTPUT, "Setting catalog parameters"); readcatparams(prefs.param_name); useprefs(); /* update things accor. to prefs parameters */ /* Check if a specific extension should be loaded */ if ((nima0=selectext(prefs.image_name[0])) != RETURN_ERROR) { forcextflag = 1; ntabmax = next = 1; } else forcextflag = 0; /* Do the same for other data (but do not force single extension mode) */ nima1 = selectext(prefs.image_name[1]); nweight0 = selectext(prefs.wimage_name[0]); nweight1 = selectext(prefs.wimage_name[1]); if (prefs.dpsf_flag) { npsf0 = selectext(prefs.psf_name[0]); npsf1 = selectext(prefs.psf_name[1]); } else npsf0 = selectext(prefs.psf_name[0]); for (i=0; inparam, 1); changecatparamarrays("VECTOR_MODELERR", &theprofit->nparam, 1); nparam2[0] = nparam2[1] = theprofit->nparam; changecatparamarrays("MATRIX_MODELERR", nparam2, 2); if (prefs.dprof_flag) thedprofit = profit_init(thedpsf, modeltype); if (prefs.pattern_flag) { npat0 = prefs.prof_disk_patternvectorsize; if (npat0size[2]; changecatparamarrays("DISK_PATTERN_VECTOR", &npat, 1); } if (FLAG(obj2.prof_disk_patternmodvector)) { npat = pattern->ncomp*pattern->nfreq; changecatparamarrays("DISK_PATTERNMOD_VECTOR", &npat, 1); } if (FLAG(obj2.prof_disk_patternargvector)) { npat = pattern->ncomp*pattern->nfreq; changecatparamarrays("DISK_PATTERNARG_VECTOR", &npat, 1); } pattern_end(pattern); } QPRINTF(OUTPUT, "Fitting model: "); for (i=0; inprof; i++) { if (i) QPRINTF(OUTPUT, "+"); QPRINTF(OUTPUT, "%s", theprofit->prof[i]->name); } QPRINTF(OUTPUT, "\n"); if (FLAG(obj2.prof_concentration)|FLAG(obj2.prof_concentration)) { thepprofit = profit_init(thepsf, MODEL_DIRAC); theqprofit = profit_init(thepsf, MODEL_EXPONENTIAL); } #else error(EXIT_FAILURE, "*Error*: model-fitting is not supported in this build.\n", " Please check your configure options"); #endif } if (prefs.filter_flag) { NFPRINTF(OUTPUT, "Reading detection filter"); getfilter(prefs.filter_name); /* get the detection filter */ } if (FLAG(obj2.sprob)) { NFPRINTF(OUTPUT, "Initializing Neural Network"); neurinit(); NFPRINTF(OUTPUT, "Reading Neural Network Weights"); getnnw(); } if (prefs.somfit_flag) { int margin; thesom = som_load(prefs.som_name); if ((margin=(thesom->inputsize[1]+1)/2) > prefs.cleanmargin) prefs.cleanmargin = margin; if (prefs.somfit_vectorsize>thesom->neurdim) { prefs.somfit_vectorsize = thesom->neurdim; sprintf(gstr,"%d", prefs.somfit_vectorsize); warning("Dimensionality of the SOM-fit vector limited to ", gstr); } } /* Prepare growth-curve buffer */ if (prefs.growth_flag) initgrowth(); /* Allocate memory for multidimensional catalog parameter arrays */ alloccatparams(); useprefs(); if (!(imacat = read_cat(prefs.image_name[0]))) error(EXIT_FAILURE, "*Error*: cannot open ", prefs.image_name[0]); close_cat(imacat); imatab = imacat->tab; if (!forcextflag) { ntabmax = imacat->ntab; /*-- Compute the number of valid input extensions */ next = 0; for (ntab = 0 ; ntabntab; ntab++, imatab = imatab->nexttab) { /*---- Check for the next valid image extension */ if ((imatab->naxis < 2) || !strncmp(imatab->xtension, "BINTABLE", 8) || !strncmp(imatab->xtension, "ASCTABLE", 8)) if (!imatab->isTileCompressed) continue; next++; } } thecat.next = next; /*-- Init the CHECK-images */ if (prefs.check_flag) { checkenum c; NFPRINTF(OUTPUT, "Initializing check-image(s)"); for (i=0; inexttab) { printf("matab->naxis %d\ imatab->xtension = %s compressed? %d\n", imatab->naxis, imatab->xtension, imatab->isTileCompressed); /*-- Check for the next valid image extension */ if (!forcextflag && ((imatab->naxis < 2) || !strncmp(imatab->xtension, "BINTABLE", 8) || !strncmp(imatab->xtension, "ASCTABLE", 8))) if (!imatab->isTileCompressed) continue; nok++; /*-- Initial time measurement*/ time(&thetime1); thecat.currext = nok; dfield = field = wfield = dwfield = NULL; if (prefs.dimage_flag) { /*---- Init the Detection and Measurement-images */ dfield = newfield(prefs.image_name[0], DETECT_FIELD, nima0<0? ntab:nima0); field = newfield(prefs.image_name[1], MEASURE_FIELD, nima1<0? ntab:nima1); if ((field->width!=dfield->width) || (field->height!=dfield->height)) error(EXIT_FAILURE, "*Error*: Frames have different sizes",""); /*---- Prepare interpolation */ if (prefs.dweight_flag && prefs.interp_type[0] == INTERP_ALL) init_interpolate(dfield, -1, -1); if (prefs.interp_type[1] == INTERP_ALL) init_interpolate(field, -1, -1); } else { field = newfield(prefs.image_name[0], DETECT_FIELD | MEASURE_FIELD, nima0<0? ntab:nima0); /*---- Prepare interpolation */ if ((prefs.dweight_flag || prefs.weight_flag) && prefs.interp_type[0] == INTERP_ALL) init_interpolate(field, -1, -1); /* 0.0 or anything else */ } /*-- Init the WEIGHT-images */ if (prefs.dweight_flag || prefs.weight_flag) { weightenum wtype; PIXTYPE interpthresh; if (prefs.nweight_type>1) { /*------ Double-weight-map mode */ if (prefs.weight_type[1] != WEIGHT_NONE) { /*-------- First: the "measurement" weights */ wfield = newweight(prefs.wimage_name[1],field,prefs.weight_type[1], (nima1<0 && prefs.image_name[1])? ntab : (nweight1<0?1:nweight1)); wtype = prefs.weight_type[1]; interpthresh = prefs.weight_thresh[1]; /*-------- Convert the interpolation threshold to variance units */ weight_to_var(wfield, &interpthresh, 1); wfield->weight_thresh = interpthresh; if (prefs.interp_type[1] != INTERP_NONE) init_interpolate(wfield, prefs.interp_xtimeout[1], prefs.interp_ytimeout[1]); } /*------ The "detection" weights */ if (prefs.weight_type[0] != WEIGHT_NONE) { interpthresh = prefs.weight_thresh[0]; if (prefs.weight_type[0] == WEIGHT_FROMINTERP) { dwfield=newweight(prefs.wimage_name[0],wfield,prefs.weight_type[0], nima0<0? ntab : (nweight0<0? 1 :nweight0)); weight_to_var(wfield, &interpthresh, 1); } else { dwfield = newweight(prefs.wimage_name[0], dfield?dfield:field, prefs.weight_type[0], nima0<0? ntab : (nweight0<0?1:nweight0)); weight_to_var(dwfield, &interpthresh, 1); } dwfield->weight_thresh = interpthresh; if (prefs.interp_type[0] != INTERP_NONE) init_interpolate(dwfield, prefs.interp_xtimeout[0], prefs.interp_ytimeout[0]); } } else { /*------ Single-weight-map mode */ wfield = newweight(prefs.wimage_name[0], dfield?dfield:field, prefs.weight_type[0], nima0<0? ntab : (nweight0<0?1:nweight0)); wtype = prefs.weight_type[0]; interpthresh = prefs.weight_thresh[0]; /*------ Convert the interpolation threshold to variance units */ weight_to_var(wfield, &interpthresh, 1); wfield->weight_thresh = interpthresh; if (prefs.interp_type[0] != INTERP_NONE) init_interpolate(wfield, prefs.interp_xtimeout[0], prefs.interp_ytimeout[0]); } } /*-- Init the FLAG-images */ for (i=0; iwidth!=field->width) || (pffield[i]->height!=field->height)) error(EXIT_FAILURE, "*Error*: Incompatible FLAG-map size in ", prefs.fimage_name[i]); } /*-- Compute background maps for `standard' fields */ QPRINTF(OUTPUT, dfield? "Measurement image:" : "Detection+Measurement image: "); makeback(field, wfield, prefs.wscale_flag[1]); QPRINTF(OUTPUT, (dfield || (dwfield&&dwfield->flags^INTERP_FIELD))? "(M) " "Background: %-10g RMS: %-10g / Threshold: %-10g \n" : "(M+D) " "Background: %-10g RMS: %-10g / Threshold: %-10g \n", field->backmean, field->backsig, (field->flags & DETECT_FIELD)? field->dthresh: field->thresh); if (dfield) { QPRINTF(OUTPUT, "Detection image: "); makeback(dfield, dwfield? dwfield : (prefs.weight_type[0] == WEIGHT_NONE?NULL:wfield), prefs.wscale_flag[0]); QPRINTF(OUTPUT, "(D) " "Background: %-10g RMS: %-10g / Threshold: %-10g \n", dfield->backmean, dfield->backsig, dfield->dthresh); } else if (dwfield && dwfield->flags^INTERP_FIELD) { makeback(field, dwfield, prefs.wscale_flag[0]); QPRINTF(OUTPUT, "(D) " "Background: %-10g RMS: %-10g / Threshold: %-10g \n", field->backmean, field->backsig, field->dthresh); } /*-- For interpolated weight-maps, copy the background structure */ if (dwfield && dwfield->flags&(INTERP_FIELD|BACKRMS_FIELD)) copyback(dwfield->reffield, dwfield); if (wfield && wfield->flags&(INTERP_FIELD|BACKRMS_FIELD)) copyback(wfield->reffield, wfield); /*-- Prepare learn and/or associations */ if (prefs.assoc_flag) init_assoc(field); /* initialize assoc tasks */ /*-- Update the CHECK-images */ if (prefs.check_flag) for (i=0; i1) { if (prefs.psf_flag) { /*------ Read other PSF extensions */ NFPRINTF(OUTPUT, "Reading PSF information"); psf_end(thepsf, thepsfit); if (prefs.dpsf_flag) { psf_end(thedpsf, thedpsfit); thedpsf = psf_load(prefs.psf_name[0], nok); thepsf = psf_load(prefs.psf_name[1], nok); } else thepsf = psf_load(prefs.psf_name[0], nok); } if (prefs.prof_flag) { /*------ Create profiles at full resolution */ profit_end(theprofit); theprofit = profit_init(thepsf, modeltype); if (prefs.dprof_flag) { profit_end(thedprofit); thedprofit = profit_init(thedpsf, modeltype); } if (prefs.pattern_flag) { pattern = pattern_init(theprofit, prefs.pattern_type, npat0); pattern_end(pattern); } if (FLAG(obj2.prof_concentration)|FLAG(obj2.prof_concentration)) { profit_end(thepprofit); profit_end(theqprofit); thepprofit = profit_init(thepsf, MODEL_DIRAC); theqprofit = profit_init(thepsf, MODEL_EXPONENTIAL); } } } /*-- Initialize PSF contexts and workspace */ if (prefs.psf_flag) { psf_readcontext(thepsf, field); psf_init(); if (prefs.dpsf_flag) { psf_readcontext(thepsf, dfield); psf_init(); } } /*-- Copy field structures to static ones (for catalog info) */ if (dfield) { thefield1 = *field; thefield2 = *dfield; } else thefield1 = thefield2 = *field; if (wfield) { thewfield1 = *wfield; thewfield2 = dwfield? *dwfield: *wfield; } else if (dwfield) thewfield2 = *dwfield; reinitcat(field); /*-- Start the extraction pipeline */ NFPRINTF(OUTPUT, "Scanning image"); scanimage(field, dfield, pffield, prefs.nimaflag, wfield, dwfield); /*-- Finish the current CHECK-image processing */ if (prefs.check_flag) for (i=0; itm_year+1900, tm->tm_mon+1, tm->tm_mday); sprintf(prefs.stime_end,"%02d:%02d:%02d", tm->tm_hour, tm->tm_min, tm->tm_sec); prefs.time_diff = counter_seconds() - dtime; /* Write XML */ if (prefs.xml_flag) write_xml(prefs.xml_name); endcat((char *)NULL); if (prefs.xml_flag || prefs.cat_type==ASCII_VO) end_xml(); return; } /******************************** initglob ***********************************/ /* Initialize a few global variables */ void initglob() { int i; for (i=0; i<37; i++) { ctg[i] = cos(i*PI/18); stg[i] = sin(i*PI/18); } return; } /****** selectext ************************************************************ PROTO int selectext(char *filename) PURPOSE Return the user-selected extension number [%d] from the file name. INPUT Filename character string. OUTPUT Extension number, or RETURN_ERROR if nos extension specified. NOTES The bracket and its extension number are removed from the filename if found. AUTHOR E. Bertin (IAP) VERSION 08/10/2007 ***/ static int selectext(char *filename) { char *bracl,*bracr; int next; if (filename && (bracl=strrchr(filename, '['))) { *bracl = '\0'; if ((bracr=strrchr(bracl+1, ']'))) *bracr = '\0'; next = strtol(bracl+1, NULL, 0); return next; } return RETURN_ERROR; } /****** write_error ******************************************************** PROTO int write_error(char *msg1, char *msg2) PURPOSE Manage files in case of a catched error INPUT a character string, another character string OUTPUT RETURN_OK if everything went fine, RETURN_ERROR otherwise. NOTES -. AUTHOR E. Bertin (IAP) VERSION 14/07/2006 ***/ void write_error(char *msg1, char *msg2) { char error[MAXCHAR]; sprintf(error, "%s%s", msg1,msg2); if (prefs.xml_flag) write_xmlerror(prefs.xml_name, error); /* Also close existing catalog */ endcat(error); end_xml(); return; }