/* * check.c * * Manage "check-images". * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * * This file part of: SExtractor * * Copyright: (C) 1993-2011 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: 26/02/2013 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ #ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include #include #include #include "define.h" #include "globals.h" #include "fits/fitscat.h" #include "fitswcs.h" #include "check.h" /********************************* addcheck **********************************/ /* Add a PSF to a CHECK-image (with a multiplicative factor). Outside boundaries are taken into account. */ void addcheck(checkstruct *check, float *psf, int w,int h, int ix,int iy, float amplitude) { PIXTYPE *pix; int x,y, xmin,xmax,ymin,ymax,w2, dwpsf; /* Set the image boundaries */ w2 = w; ymin = iy-h/2; ymax = ymin + h; if (ymin<0) { psf -= ymin*w; ymin = 0; } if (ymax>check->height) ymax = check->height; xmin = ix-w/2; xmax = xmin + w; if (xmax>check->width) { w2 -= xmax-check->width; xmax = check->width; } if (xmin<0) { psf -= xmin; w2 += xmin; xmin = 0; } dwpsf = w-w2; /* Subtract the right pixels to the destination */ for (y=ymin; ypix+y*check->width+xmin; for (x=w2; x--;) *(pix++) += amplitude**(psf++); } return; } /****** addcheck_resample ***************************************************** PROTO void addcheck_resample(checkstruct *check, float *thumb, int w, int h, int ix,int iy, float zoom, float amplitude) PURPOSE Add a resampled thumbnail to a check image (with a multiplicative factor). INPUT Pointer to the check-image, pointer to the thumbnail, thumbnail width, thumbnail height, thumbnail center x coordinate, thumbnail center y coordinate, zoom factor, flux scaling factor. OUTPUT -. NOTES Outside boundaries are taken into account. AUTHOR E. Bertin (IAP) VERSION 26/02/2013 ***/ void addcheck_resample(checkstruct *check, float *thumb, int w, int h, int ix, int iy, float step2, float amplitude) { float interpm[CHECKINTERPW*CHECKINTERPW], *pix2, *mask,*maskt,*pix12, *pixin,*pixin0, *pixout,*pixout0, xc1,xc2,yc1,yc2, xs1,ys1, x1,y1, x,y, dxm,dym, val, norm, dx,dy; int i,j,k,n,t, *start,*startt, *nmask,*nmaskt, ixs2,iys2, ix2,iy2, dix2,diy2, nx2,ny2, iys1a, ny1, hmw,hmh, ixs,iys, ix1,iy1, w2,h2; pix2 = check->pix; w2 = check->width; h2 = check->height; dx = 0.5*(1.0 - step2); dy = 0.5*(1.0 - step2); xc1 = (float)(w/2); /* Im1 center x-coord*/ xc2 = (float)ix; /* Im2 center x-coord*/ xs1 = xc1 + dx - xc2*step2; /* Im1 start x-coord */ if ((int)xs1 >= w) return; ixs2 = 0; /* Int part of Im2 start x-coord */ if (xs1<0.0) { dix2 = (int)(1-xs1/step2); /*-- Simply leave here if the images do not overlap in x */ if (dix2 >= w2) return; ixs2 += dix2; xs1 += dix2*step2; } nx2 = (int)((w-1-xs1)/step2+1);/* nb of interpolated Im2 pixels along x */ if (nx2>(ix2=w2-ixs2)) nx2 = ix2; if (nx2<=0) return; yc1 = (float)(h/2); /* Im1 center y-coord */ yc2 = (float)iy; /* Im2 center y-coord */ ys1 = yc1 + dy - yc2*step2; /* Im1 start y-coord */ if ((int)ys1 >= h) return; iys2 = 0; /* Int part of Im2 start y-coord */ if (ys1<0.0) { diy2 = (int)(1-ys1/step2); /*-- Simply leave here if the images do not overlap in y */ if (diy2 >= h2) return; iys2 += diy2; ys1 += diy2*step2; } ny2 = (int)((h-1-ys1)/step2+1);/* nb of interpolated Im2 pixels along y */ if (ny2>(iy2=h2-iys2)) ny2 = iy2; if (ny2<=0) return; /* Set the yrange for the x-resampling with some margin for interpolation */ iys1a = (int)ys1; /* Int part of Im1 start y-coord with margin */ hmh = CHECKINTERPW/2 - 1; /* Interpolant start */ if (iys1a<0 || ((iys1a -= hmh)< 0)) iys1a = 0; ny1 = (int)(ys1+ny2*step2)+CHECKINTERPW-hmh; /* Interpolated Im1 y size */ if (ny1>h) /* with margin */ ny1 = h; /* Express everything relative to the effective Im1 start (with margin) */ ny1 -= iys1a; ys1 -= (float)iys1a; /* Allocate interpolant stuff for the x direction */ QMALLOC(mask, float, nx2*CHECKINTERPW); /* Interpolation masks */ QMALLOC(nmask, int, nx2); /* Interpolation mask sizes */ QMALLOC(start, int, nx2); /* Int part of Im1 conv starts */ /* Compute the local interpolant and data starting points in x */ hmw = CHECKINTERPW/2 - 1; x1 = xs1; maskt = mask; nmaskt = nmask; startt = start; for (j=nx2; j--; x1+=step2) { ixs = (ix1=(int)x1) - hmw; dxm = ix1 - x1 - hmw; /* starting point in the interpolation func */ if (ixs < 0) { n = CHECKINTERPW+ixs; dxm -= (float)ixs; ixs = 0; } else n = CHECKINTERPW; if (n>(t=w-ixs)) n=t; *(startt++) = ixs; *(nmaskt++) = n; norm = 0.0; for (x=dxm, i=n; i--; x+=1.0) norm += (*(maskt++) = CHECKINTERPF(x)); norm = norm>0.0? 1.0/norm : 1.0; maskt -= n; for (i=n; i--;) *(maskt++) *= norm; } QCALLOC(pix12, float, nx2*ny1); /* Intermediary frame-buffer */ /* Make the interpolation in x (this includes transposition) */ pixin0 = thumb+iys1a*w; pixout0 = pix12; for (k=ny1; k--; pixin0+=w, pixout0++) { maskt = mask; nmaskt = nmask; startt = start; pixout = pixout0; for (j=nx2; j--; pixout+=ny1) { pixin = pixin0+*(startt++); val = 0.0; for (i=*(nmaskt++); i--;) val += *(maskt++)**(pixin++); *pixout = val; } } /* Reallocate interpolant stuff for the y direction */ QREALLOC(mask, float, ny2*CHECKINTERPW); /* Interpolation masks */ QREALLOC(nmask, int, ny2); /* Interpolation mask sizes */ QREALLOC(start, int, ny2); /* Int part of Im1 conv starts */ /* Compute the local interpolant and data starting points in y */ hmh = CHECKINTERPW/2 - 1; y1 = ys1; maskt = mask; nmaskt = nmask; startt = start; for (j=ny2; j--; y1+=step2) { iys = (iy1=(int)y1) - hmh; dym = iy1 - y1 - hmh; /* starting point in the interpolation func */ if (iys < 0) { n = CHECKINTERPW+iys; dym -= (float)iys; iys = 0; } else n = CHECKINTERPW; if (n>(t=ny1-iys)) n=t; *(startt++) = iys; *(nmaskt++) = n; norm = 0.0; for (y=dym, i=n; i--; y+=1.0) norm += (*(maskt++) = CHECKINTERPF(y)); norm = norm>0.0? 1.0/norm : 1.0; maskt -= n; for (i=n; i--;) *(maskt++) *= norm; } /* Make the interpolation in y and transpose once again */ pixin0 = pix12; pixout0 = pix2+ixs2+iys2*w2; for (k=nx2; k--; pixin0+=ny1, pixout0++) { maskt = mask; nmaskt = nmask; startt = start; pixout = pixout0; for (j=ny2; j--; pixout+=w2) { pixin = pixin0+*(startt++); val = 0.0; for (i=*(nmaskt++); i--;) val += *(maskt++)**(pixin++); *pixout += amplitude*val; } } /* Free memory */ free(pix12); free(mask); free(nmask); free(start); return; } /********************************* blankcheck *******************************/ /* Blank a part of the CHECK-image according to a mask. */ void blankcheck(checkstruct *check, PIXTYPE *mask, int w,int h, int xmin,int ymin, PIXTYPE val) { PIXTYPE *pixt; int x,y, xmax,ymax,w2,wc; /* Don't go further if out of frame!! */ if (xmin+w<0 || xmin>=check->width || ymin+h<0 || ymin>=check->height) return; /* Set the image boundaries */ w2 = w; ymax = ymin + h; if (ymin<0) { mask -= ymin*w; ymin = 0; } if (ymax>check->height) ymax = check->height; xmax = xmin + w; if (xmax>check->width) { w2 -= xmax - check->width; xmax = check->width; } if (xmin<0) { mask += -xmin; w2 -= -xmin; xmin = 0; } w -= w2; wc = check->width; ymin = ymin*wc+xmin; ymax = ymax*wc+xmin; /* Blank the right pixels in the image */ for (y=ymin; ypix + y; for (x=w2; x--; pixt++) if (*(mask++) > -BIG) *pixt = val; } return; } /******************************** initcheck **********************************/ /* initialize check-image. */ checkstruct *initcheck(char *filename, checkenum check_type, int next) { catstruct *cat; checkstruct *check; QCALLOC(check, checkstruct, 1); check->type = check_type; check->next = next; cat = check->cat = new_cat(1); strcpy(cat->filename, filename); if (next>1) /*-- Create a "pure" primary HDU */ { init_cat(cat); addkeywordto_head(cat->tab, "NEXTEND ", "Number of extensions"); fitswrite(cat->tab->headbuf, "NEXTEND ", &next, H_INT, T_LONG); if (open_cat(cat, WRITE_ONLY) != RETURN_OK) error(EXIT_FAILURE,"*Error*: cannot open for writing ", filename); save_head(cat, cat->tab); remove_tabs(cat); } else open_cat(cat, WRITE_ONLY); return check; } /******************************** reinitcheck ********************************/ /* initialize check-image (for subsequent writing). */ void reinitcheck(picstruct *field, checkstruct *check) { catstruct *cat; tabstruct *tab; wcsstruct *wcs; char *fitshead; PIXTYPE *ptrf; double dval; int i; cat = check->cat; /* Inherit the field FITS header */ remove_tabs(cat); copy_tab_fromptr(field->tab, cat, 0); tab = cat->tab; tab->cat = cat; if (check->next<=1) prim_head(tab); check->y = 0; fitshead = tab->headbuf; /* Neutralize possible scaling factors */ tab->bscale = 1.0; tab->bzero = 0.0; fitswrite(fitshead, "BSCALE ", &tab->bscale, H_FLOAT, T_DOUBLE); fitswrite(fitshead, "BZERO ", &tab->bzero, H_FLOAT, T_DOUBLE); fitswrite(fitshead, "BITSGN ", &tab->bitsgn, H_INT, T_LONG); if (tab->compress_type != COMPRESS_NONE) { tab->compress_type = COMPRESS_NONE; fitswrite(fitshead, "IMAGECOD", "NONE", H_STRING, T_STRING); } fitswrite(fitshead, "ORIGIN ", BANNER, H_STRING, T_STRING); switch(check->type) { case CHECK_IDENTICAL: case CHECK_BACKGROUND: case CHECK_FILTERED: case CHECK_SUBTRACTED: tab->bitpix = BP_FLOAT; tab->bytepix = 4; tab->bitsgn = 0; tab->naxisn[0] = check->width = field->width; tab->naxisn[1] = check->height = field->height; check->npix = field->npix; QMALLOC(ptrf, PIXTYPE, check->width); check->pix = (void *)ptrf; save_head(cat, cat->tab); break; case CHECK_BACKRMS: case CHECK_SUBOBJECTS: tab->bitpix = BP_FLOAT; tab->bytepix = 4; tab->bitsgn = 0; tab->naxisn[0] = check->width = field->width; tab->naxisn[1] = check->height = field->height; check->npix = field->npix; QMALLOC(check->pix, PIXTYPE, check->width); save_head(cat, cat->tab); /*---- Allocate memory for replacing the blanked pixels by 0 */ if (!check->line) QMALLOC(check->line, PIXTYPE, field->width); break; case CHECK_OBJECTS: case CHECK_APERTURES: case CHECK_PSFPROTOS: case CHECK_SUBPSFPROTOS: case CHECK_PCPROTOS: case CHECK_SUBPCPROTOS: case CHECK_PCOPROTOS: case CHECK_PROFILES: case CHECK_SUBPROFILES: case CHECK_SPHEROIDS: case CHECK_SUBSPHEROIDS: case CHECK_DISKS: case CHECK_SUBDISKS: case CHECK_PATTERNS: tab->bitpix = -32; tab->bytepix = 4; tab->bitsgn = 0; tab->naxisn[0] = check->width = field->width; tab->naxisn[1] = check->height = field->height; check->npix = field->npix; check->overlay = 30*field->backsig; QCALLOC(check->pix, PIXTYPE, check->npix); save_head(cat, cat->tab); break; case CHECK_SEGMENTATION: tab->bitpix = BP_LONG; tab->bytepix = 4; tab->bitsgn = 1; tab->naxisn[0] = check->width = field->width; tab->naxisn[1] = check->height = field->height; check->npix = field->npix; QCALLOC(check->pix, ULONG, check->npix); save_head(cat, cat->tab); break; case CHECK_MASK: case CHECK_SUBMASK: tab->bitpix = BP_BYTE; tab->bytepix = 1; tab->bitsgn = 1; tab->naxisn[0] = check->width = field->width; tab->naxisn[1] = check->height = field->height; check->npix = field->npix; save_head(cat, cat->tab); /*---- Allocate memory */ if (!check->line) QMALLOC(check->line, FLAGTYPE, check->width); break; case CHECK_ASSOC: tab->bitpix = BP_FLOAT; tab->bytepix = 4; tab->bitsgn = 0; tab->naxisn[0] = check->width = field->width; tab->naxisn[1] = check->height = field->height; check->npix = field->npix; QMALLOC(check->pix, PIXTYPE, check->npix); /*---- Initialize the pixmap to IEEE NaN */ memset(check->pix, 0xFF, check->npix*sizeof(LONG)); save_head(cat, cat->tab); break; case CHECK_MINIBACKGROUND: case CHECK_MINIBACKRMS: tab->bitpix = BP_FLOAT; tab->bytepix = 4; tab->bitsgn = 0; tab->naxisn[0] = check->width = field->nbackx; tab->naxisn[1] = check->height = field->nbacky; /*---- Scale the WCS information if present */ if ((wcs=field->wcs)) { dval = wcs->cdelt[0]*field->backw; fitswrite(fitshead, "CDELT1 ", &dval, H_EXPO, T_DOUBLE); dval = wcs->cdelt[1]*field->backh; fitswrite(fitshead, "CDELT2 ", &dval, H_EXPO, T_DOUBLE); dval = (wcs->crpix[0]-0.5)/field->backw + 0.5; fitswrite(fitshead, "CRPIX1 ", &dval, H_EXPO, T_DOUBLE); dval = (wcs->crpix[1]-0.5)/field->backh + 0.5; fitswrite(fitshead, "CRPIX2 ", &dval, H_EXPO, T_DOUBLE); dval = wcs->cd[0]*field->backw; fitswrite(fitshead, "CD1_1 ", &dval, H_EXPO, T_DOUBLE); dval = wcs->cd[1]*field->backh; fitswrite(fitshead, "CD1_2 ", &dval, H_EXPO, T_DOUBLE); dval = wcs->cd[wcs->naxis]*field->backw; fitswrite(fitshead, "CD2_1 ", &dval, H_EXPO, T_DOUBLE); dval = wcs->cd[wcs->naxis+1]*field->backh; fitswrite(fitshead, "CD2_2 ", &dval, H_EXPO, T_DOUBLE); } check->npix = check->width*check->height; QMALLOC(check->pix, PIXTYPE, check->npix); if (check->type==CHECK_MINIBACKRMS) memcpy(check->pix, field->sigma, check->npix*sizeof(float)); else memcpy(check->pix, field->back, check->npix*sizeof(float)); save_head(cat, cat->tab); write_body(cat->tab, check->pix, check->npix); free(check->pix); break; case CHECK_MAPSOM: tab->bitpix = BP_FLOAT; tab->bytepix = 4; tab->bitsgn = 0; tab->naxisn[0] = check->width = field->width; tab->naxisn[1] = check->height = field->height; check->npix = field->npix; QMALLOC(ptrf, PIXTYPE, check->npix); check->pix = (void *)ptrf; for (i=check->npix; i--;) *(ptrf++) = -10.0; save_head(cat, cat->tab); break; case CHECK_OTHER: tab->bitpix = BP_FLOAT; tab->bytepix = 4; tab->bitsgn = 0; tab->naxisn[0] = check->width; tab->naxisn[1] = check->height; check->npix = check->width*check->height; QCALLOC(check->pix, PIXTYPE, check->npix); save_head(cat, cat->tab); break; default: error(EXIT_FAILURE, "*Internal Error* in ", "reinitcheck()!"); } return; } /******************************** writecheck *********************************/ /* Write ONE line of npix pixels of a check-image. */ void writecheck(checkstruct *check, PIXTYPE *data, int w) { if (check->type == CHECK_APERTURES || check->type == CHECK_SUBPSFPROTOS || check->type == CHECK_SUBPCPROTOS || check->type == CHECK_PCOPROTOS || check->type == CHECK_SUBPROFILES || check->type == CHECK_SUBSPHEROIDS || check->type == CHECK_SUBDISKS || check->type == CHECK_OTHER) { memcpy((PIXTYPE *)check->pix + w*(check->y++), data, w*sizeof(PIXTYPE)); return; } else if (check->type == CHECK_SUBOBJECTS) { int i; PIXTYPE *pixt; pixt = (PIXTYPE *)check->line; for (i=w; i--; data++) *(pixt++) = (*data>-BIG)? *data:0.0; write_body(check->cat->tab, (PIXTYPE *)check->line, w); } else if (check->type == CHECK_MASK) { int i; FLAGTYPE *pixt; pixt = (FLAGTYPE *)check->line; for (i=w; i--;) *(pixt++) = (*(data++)>-BIG)?0:1; write_ibody(check->cat->tab, (FLAGTYPE *)check->line, w); } else if (check->type == CHECK_SUBMASK) { int i; FLAGTYPE *pixt; pixt = (FLAGTYPE *)check->line; for (i=w; i--;) *(pixt++) = (*(data++)>-BIG)?1:0; write_ibody(check->cat->tab, (FLAGTYPE *)check->line, w); } else write_body(check->cat->tab, data, w); return; } /********************************* reendcheck ********************************/ /* Finish current check-image. */ void reendcheck(picstruct *field, checkstruct *check) { catstruct *cat; cat = check->cat; switch(check->type) { case CHECK_MINIBACKGROUND: case CHECK_MINIBACKRMS: return; case CHECK_IDENTICAL: case CHECK_BACKGROUND: case CHECK_BACKRMS: case CHECK_FILTERED: case CHECK_SUBTRACTED: free(check->pix); free(check->line); check->line = NULL; pad_tab(cat, check->npix*sizeof(PIXTYPE)); break; case CHECK_OBJECTS: case CHECK_APERTURES: case CHECK_PSFPROTOS: case CHECK_SUBPSFPROTOS: case CHECK_PCPROTOS: case CHECK_SUBPCPROTOS: case CHECK_PCOPROTOS: case CHECK_PROFILES: case CHECK_SUBPROFILES: case CHECK_SPHEROIDS: case CHECK_SUBSPHEROIDS: case CHECK_DISKS: case CHECK_SUBDISKS: case CHECK_ASSOC: case CHECK_PATTERNS: case CHECK_MAPSOM: case CHECK_OTHER: write_body(cat->tab, check->pix, check->npix); free(check->pix); pad_tab(cat, check->npix*sizeof(PIXTYPE)); break; case CHECK_SEGMENTATION: write_ibody(cat->tab, check->pix, check->npix); free(check->pix); pad_tab(cat, check->npix*sizeof(FLAGTYPE)); break; case CHECK_MASK: case CHECK_SUBMASK: { int y; for (y=field->ymin; yymax; y++) writecheck(check, &PIX(field, 0, y), field->width); free(check->line); check->line = NULL; pad_tab(cat, check->npix*sizeof(unsigned char)); break; } case CHECK_SUBOBJECTS: { int y; for (y=field->ymin; yymax; y++) writecheck(check, &PIX(field, 0, y), field->width); free(check->pix); free(check->line); check->line = NULL; pad_tab(cat, check->npix*sizeof(PIXTYPE)); break; } default: error(EXIT_FAILURE, "*Internal Error* in ", "endcheck()!"); } return; } /********************************* endcheck **********************************/ /* close check-image. */ void endcheck(checkstruct *check) { free_cat(&check->cat,1); free(check); return; }