/* @(#)mo_fitoff.c 17.1.1.1 (ESO-DAG) 01/25/02 17:49:06 /*=========================================================================== Copyright (C) 1995 European Southern Observatory (ESO) 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. This program 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 this program; if not, write to the Free Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, MA 02139, USA. Corresponding concerning ESO-MIDAS should be addressed as follows: Internet e-mail: midas@eso.org Postal address: European Southern Observatory Data Management Division Karl-Schwarzschild-Strasse 2 D 85748 Garching bei Muenchen GERMANY ===========================================================================*/ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ .IDENTIFICATION: routine mo_fitoff .KEYWORDS: bulk data frames, mosaicing, offset .PURPOSE: Finds the relative offsets in background intensity between every pair of overlapping frames in a mosaic. Construct a mosaic of a number of image frames, the result will a frame with nx subframes in x and ny subframes in y. The subframes must have the same size in x and y. Blank subframes can be indicated on the command line. .ALGORITHM: Find the best fitting background offsets for each frame in a mosaic. This routine uses a chi-square minimization to try to minimize the residual errors between relative offsets. The routine findoff must be run first to create a file that contains the offsets between each pair of frames. This routine has three inputs: input_array containing the output from mo_offset; exclude_array with subraster not be be included; output_table with the resulting offset for each subraster; Besides the output table, this routine will output the offset for each frame and the average of the square of residuals for each frame. At the end it will print the mean of these and the standard deviation. Bad frames will have large residuals and can be excluded from the fit by adding them to the exclude file and rerunning the routine. .VERSION: 950920 RHW Created; original Mike Regan (U of MD,mregan@astro.umd.edu -----------------------------------------------------------------------------*/ /* * Define _POSIX_SOURCE to indicate * that this is a POSIX program */ #define _POSIX_SOURCE 1 /* * definition of the used functions */ #include #include #include #include #include #include #include #include void MO_GETOFF(); #define NPT 1000 #define MA 120 #define SPREAD 0.001 /* ++++++++++++++++++++++++++++++++++++++++++++++++++ * here starts the code of the function --------------------------------------------------- */ MO_FITOFF(npt, nover, y, x, iref, jref, excount, exclud, deltai) int npt; int nover; float *y; float *x; int *iref; int *jref; int excount; int *exclud; float *deltai; { int idum=(-911); int count; int i, io, ip; int iran; int itst; int j,k; int mfit; int numsub; int exclude; float alamda; float chisq,ochisq; float avgoff; float sum; float stdev,avgerr; int newindex[300]; int xcount; int ex[300]; int *lista; float *sig; float *cerror; float *error; float **covar; float **alpha; float **diff; char line[80]; static float a[MA+1]; /* Allocated memory */ lista = (int *) osmmget(npt * sizeof(int)); error = (float *) osmmget(SQR(NPT)/2 * sizeof(float)); cerror = (float *) osmmget(SQR(NPT)/2 * sizeof(float)); sig = (float *) osmmget(SQR(NPT)/2 * sizeof(float)); covar = (float **) osmmget(npt * sizeof(float)); alpha = (float **) osmmget(npt * sizeof(float)); for (i = 1; i <= npt; i++) { covar[i] = (float *) osmmget(npt * sizeof(float)); alpha[i] = (float *) osmmget(npt * sizeof(float)); } diff = (float **) osmmget((npt+1) * sizeof(float)); for (i = 0; i <= npt; i++) diff[i] = (float *) osmmget((npt+1) * sizeof(float)); /* Initialize the array */ for (i=1; i<=npt; i++) for (j=i+1; j<=npt; j++) diff[i][j] = 0.0; /* Stores the input in the arrays */ count = 0; for (ip = 1; ip <= nover; ip++) { i = iref[ip]; j = jref[ip]; sig[ip] = 1.0; diff[i][j] = y[ip]; cerror[ip] = 0.0; count++; } SCTPUT(" "); sprintf(line, "Pairs of subrasters before editing = %d",count); SCTPUT(line); if (excount > 0) { for (iran = 1; iran <= excount; iran++) ex[iran] = exclud[iran-1]; sprintf(line, "Number of subrasters deleted = %d", excount); SCTPUT(line); numsub = 0; for (i = 1; i <= npt; i++) { exclude = 0; for (j = 1; j <= excount; j++) { if (i == ex[j]) { exclude = 1; numsub++; } } if (exclude == 1) newindex[i] = 0; else newindex[i] = i-numsub; } for (i=1; i<=npt; i++) for (j=i+1; j<=npt; j++) diff[newindex[i]][newindex[j]] = diff[i][j]; npt = npt - excount; count = 0; xcount = 1; for (i=1; i<=npt; i++){ for (j=i+1; j<=npt; j++){ if (diff[i][j] != 0.0){ y[count+1] = diff[i][j]; x[count+1] = xcount; count++; } xcount++; } } } sprintf(line, "Pairs of subrasters after deletion = %d",count); SCTPUT(line); SCTPUT(" "); mfit = npt; for (i = 1; i <= mfit; i++) lista[i]=i; alamda = -1; for (i = 1; i <= mfit; i++) a[i]=0; MO_MRQMIN(x,y,sig,count,a,npt,lista,mfit,covar,alpha,&chisq,MO_GETOFF,&alamda); k = 1; itst = 0; while (itst < 3) { sprintf(line, "%s %2d %17s %10.4f %10s %9.2e","Iteration #",k, "chi-squared:",chisq/count,"alamda:",alamda); SCTPUT(line); k++; ochisq=chisq; MO_MRQMIN(x,y,sig,count,a,npt,lista,mfit,covar,alpha,&chisq, MO_GETOFF,&alamda); if (chisq > ochisq) itst = 0; else if (fabs(ochisq-chisq) < 0.01) itst++; } alamda=0.0; MO_MRQMIN(x,y,sig,count,a,npt,lista,mfit,covar,alpha,&chisq,MO_GETOFF,&alamda); for (i=1; i<=npt; i++){ for (j=i+1; j<=npt; j++){ if (diff[i][j] != 0.0){ error[i] = error[i] + SQR((diff[i][j]-(a[j]-a[i]))); cerror[i]++; cerror[j]++; error[j] = error[j] + SQR((diff[i][j]-(a[j]-a[i]))); } } } sum = 0; for (i=1; i<=npt; i++) sum+=a[i]; avgoff = sum/(float) npt; SCTPUT(" "); sprintf(line, "Average offset = %f",avgoff); SCTPUT(line); sum = 0; numsub = 0; for (i=1; i<=npt+excount; i++) { exclude=0; for (j=1; j<=excount; j++){ if (i == ex[j]) { exclude=1; numsub++; } } if (exclude==0) { a[i-numsub] = a[i-numsub]-avgoff; sprintf(line, "frame %d offset %6.1f error = %6.2f",i, a[i-numsub],error[i-numsub]/cerror[i-numsub]); SCTPUT(line); sum+=error[i-numsub]/cerror[i-numsub]; } } avgerr = sum/(float) npt; sum = 0; for (i=1; i<=npt; i++) sum += SQR(avgerr-(error[i]/cerror[i])); stdev = sqrt(sum/(npt-1)); sprintf(line, "Average error= %f , Std Dev = %f",avgerr,stdev); SCTPUT(line); sum = 0; count = 1; numsub = 0; a[0] = 0; for (io=0; io