/* @(#)scanima.c 17.1.1.1 (ESO-DMD) 01/25/02 17:11:17 */ /*=========================================================================== 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 Massachusetts Ave, Cambridge, MA 02139, USA. Correspondence 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 program SCANIMA version 1.00 940214 K. Banse ESO - Garching .KEYWORDS graphics, reconstruction .PURPOSE scan a 2-dim image which presents the graph of a 1-dim function and create a 1-dim image or a table as result .ALGORITHM get bottom/top pixel coords. in each column of input image (of pixel > threshold) .INPUT/OUTPUT the following keys are used: in_a/c/1/60 input frame out_a/c/1/60 result image or table p3/c/1/80 threshold or lo-,hiband p4/c/1/10 option: bottom/top p5/c/1/80 x1g,x1r,x2g,x2r p6/c/1/80 y1g,y1r,y2g,y2r p7/c/1/80 xa,ya:xb,yb p8/c/1/80 wsize,xpix1,ypix1 .VERSIONS 010710 last modif -------------------------------------------------- */ #include static int kstart, sublo[2], subhi[2]; static int wsize, xpix1, ypix1; static char optio[8]; void main() { int uni, iav, stat, jj; int tid, xcol, ycol; int imno, imnoc; int nulo, sizec, chunk; int naxis, npix[2], naxisc, npixc[2]; int ibuf[3], ll, ref_flag; register int nr; char *cpntr; char line[80], linea[80]; char frame[64], framec[64]; char cunit[64], ident[72]; char *osmmget(); double step[2], start[2]; double stepc[2], startc[2]; double ref[4], xg[2], xr[2], yg[2], yr[2]; float *pntr, *pntrc, *fpntr; float cuts[4], faux[5]; void mima(); /* set up MIDAS environment + enable automatic error abort */ SCSPRO("scanima"); for (nr=0; nr<2; nr++) { startc[nr] = 0.0; stepc[nr] = 1.0; npixc[nr] = 1; npix[nr] = 1; start[nr] = 0.0; step[nr] = 1.0; sublo[nr] = 0; } ref_flag = 0; /* get input frame, result frame, threshold(s), option and references */ stat = SCKGETC("IN_A",1,60,&iav,frame); stat = SCKGETC("OUT_A",1,60,&iav,framec); stat = SCFOPN(frame,D_R4_FORMAT,0,F_IMA_TYPE,&imno); stat = SCDRDI(imno,"NAXIS",1,1,&iav,&naxis,&uni,&nulo); stat = SCDRDI(imno,"NPIX",1,2,&iav,npix,&uni,&nulo); stat = SCDRDD(imno,"START",1,2,&iav,start,&uni,&nulo); stat = SCDRDD(imno,"STEP",1,2,&iav,step,&uni,&nulo); stat = SCDGETC(imno,"CUNIT",1,64,&iav,cunit); stat = SCDRDR(imno,"LHCUTS",1,2,&iav,cuts,&uni,&nulo); subhi[0] = npix[0] - 1; subhi[1] = npix[1] - 1; stat = SCKGETC("P3",1,80,&iav,line); if (line[0] == '+') /* no threshold given */ { faux[0] = (cuts[0] + cuts[1]) / 2; faux[4] = 1; } else { iav = CGN_CNVT(line,2,2,npixc,faux,step); if (iav < 1) SCETER(101,"invalid `threshold value' ..."); faux[4] = 2; } stat = SCKGETC("P4",1,4,&iav,line); CGN_UPCOPY(optio,line,4); jj = CGN_INDEXS(framec,",t"); /* check, if we should make a table */ if (jj < 0) jj = CGN_INDEXS(framec,",T"); if (jj > 0) { framec[jj] = '\0'; optio[4] = 'T'; } else optio[4] = 'I'; stat = SCKGETC("p5",1,80,&iav,line); /* get x,y reference values */ stat = SCKGETC("p6",1,80,&iav,linea); if ((line[0] != '+') && (linea[0] != '+')) { ref_flag = 1; iav = CGN_CNVT(line,4,4,npixc,cuts,ref); if (iav < 1) SCETER(101,"invalid `x-ref. values' ..."); xg[0] = ref[0]; xr[0] = ref[1]; xg[1] = ref[2]; xr[1] = ref[3]; iav = CGN_CNVT(linea,4,4,npixc,cuts,ref); if (iav < 1) SCETER(101,"invalid `y-ref. values' ..."); yg[0] = ref[0]; yr[0] = ref[1]; yg[1] = ref[2]; yr[1] = ref[3]; } stat = SCKGETC("p7",1,80,&iav,line); /* get subframe coords */ /* check if we have extraction coords. */ if (line[0] != '+') { if (line[0] != '[') { jj = strlen(line); for (nr=jj; nr>0; nr--) line[nr] = line[nr-1]; line[0] = '['; line[jj+1] = '\0'; } stat = Convcoo(1,imno,line,2,&jj,sublo,subhi); if (stat != 0) SCETER(102,"invalid coordinates for subframe boundaries..."); } if ((optio[1] == 'F') || (optio[1] == 'D')) { xpix1 = 0; ypix1 = 0; wsize = subhi[1] - sublo[1] + 1; stat = SCKGETC("P8",1,80,&iav,line); /* get wsize,xpix1,ypix1 */ if (line[0] != '+') { iav = CGN_CNVT(line,1,3,ibuf,cuts,ref); if (iav >= 1) { wsize = ibuf[0]; if (wsize < 1) SCETER(105,"invalid window size ..."); if (iav == 3) { xpix1 = ibuf[1] - 1; ypix1 = ibuf[2] - 1; if ((xpix1 < sublo[0]) || (xpix1 > subhi[0])) SCETER(106,"start x-pixel outside subframe ..."); if ((ypix1 < sublo[1]) || (ypix1 > subhi[1])) SCETER(106,"start y-pixel outside subframe ..."); } } } } /* get virtual memory buffers */ chunk = npix[0]*npix[1]; ll = chunk * sizeof(float); cpntr = osmmget((unsigned int)ll); if (cpntr == (char *) 0) SCETER(103,"could not allocate virtual memory..."); else pntr = (float *) cpntr; stat = SCFGET(imno,1,chunk,&iav,(char *)pntr); /* the result frame may be created larger than necessary... */ if (optio[4] == 'I') { stat = SCFCRE(framec,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,npix[0],&imnoc); stat = SCFMAP(imnoc,F_O_MODE,1,npix[0],&iav,&cpntr); } else { stat = SCFCRE("for_table",D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,npix[0],&imnoc); stat = SCFMAP(imnoc,F_X_MODE,1,npix[0],&iav,&cpntr); } pntrc =(float *) cpntr; fpntr = pntrc; /* initialize */ for (nr=0; nr 0) /* x-start has been shifted */ startc[0] += (kstart*stepc[0]); if (optio[4] == 'T') { stat = TCTINI(framec,F_RECORD,F_O_MODE,2,sizec,&tid); stat = TCCINI(tid,D_R8_FORMAT,1,"E22.12"," ","X",&xcol); stat = TCCINI(tid,D_R4_FORMAT,1,"F15.8"," ","Y",&ycol); *xr = startc[0]; fpntr = pntrc; for (nr=1; nr<(sizec+1); nr++) { if (*fpntr >= 0.0) { stat = TCEWRD(tid,nr,1,xr); stat = TCEWRR(tid,nr,2,fpntr); } *xr += stepc[0]; fpntr++; } } else { stat = SCDWRI(imnoc,"NAXIS",&naxisc,1,1,&uni); stat = SCDWRI(imnoc,"NPIX",npixc,1,naxisc,&uni); stat = SCDWRD(imnoc,"START",startc,1,naxisc,&uni); stat = SCDWRD(imnoc,"STEP",stepc,1,naxisc,&uni); strcpy(ident,"new frame "); stat = SCDWRC(imnoc,"IDENT",1,ident,1,72,&uni); jj = (naxisc+1) * 16; stat = SCDWRC(imnoc,"CUNIT",1,cunit,1,jj,&uni); mima(pntrc,sizec,cuts); cuts[2] = cuts[0]; cuts[3] = cuts[1]; cuts[1] = 0.0; cuts[2] = 0.0; stat = SCDWRR(imnoc,"LHCUTS",cuts,1,4,&uni); stat = CGN_DSCUPD(imno,imnoc," "); } SCSEPI(); } /* */ void mima(pntri,npix,cuts) float *pntri, *cuts; int npix; { register int nr; float fmi, fma; fmi = *pntri++; fma = fmi; for (nr=1; nr fma) fma = *pntri; pntri ++; } cuts[0] = fmi; cuts[1] = fma; } /* */ int get_graph(fx,pntri,npix,pntro) float *fx, *pntri, *pntro; int *npix; { register int nrx, nry; int end, mm, kk, kmax, count; int my, ky, ws2, fy, limi; char cbuf[64]; float val, valx, *lpntr, *inp, *inq, *outq, fmax; val = *fx; lpntr = pntro; /* ---------------------------------------- first we work on the fixed WINDOW option ---------------------------------------- */ if (optio[1] == 'F') { ws2 = wsize / 2; /* and split for BOTTOM or TOP direction */ if (optio[0] == 'B') { inp = pntri + (sublo[1]*npix[0]); /* final split for Maximum or Value option */ if (optio[2] == 'M') /* search for max value */ { if (xpix1 < 1) /* no start pixel given */ { fmax = val*0.99; fy = -1; for (nrx=sublo[0]; nrx<=subhi[0]; nrx++) { outq = lpntr + nrx; inq = inp + nrx; for (nry=sublo[1]; nry<= subhi[1]; nry++) { if (*inq > fmax) { fy = nry; fmax = *inq; } inq += npix[0]; /* move to next line */ } if (fy > -1) { *outq = fy; xpix1 = nrx; ypix1 = fy; goto win_1; } } SCETER(110,"no start pixel for window option found... "); } win_1: outq = lpntr + xpix1 + 1; /* point to next element */ for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { fmax = val*0.99; fy = -1; my = ypix1 - ws2; if (my < sublo[1]) my = sublo[1]; ky = my + wsize; if (ky > subhi[1]) ky = subhi[1] + 1; inq = pntri + nrx + (my*npix[0]); /* point to low window start */ for (nry=my; nry fmax) { fy = nry; fmax = *inq; } inq += npix[0]; } if (fy > -1) { *outq = fy; ypix1 = fy; } outq ++; } } else { if (xpix1 < 1) /* no start pixel given */ { for (nrx=sublo[0]; nrx<=subhi[0]; nrx++) { outq = lpntr + nrx; inq = inp + nrx; for (nry=sublo[1]; nry<= subhi[1]; nry++) { if (*inq >= val) { *outq = nry; xpix1 = nrx; ypix1 = nry; goto win_11; } inq += npix[0]; /* move to next line */ } } SCETER(110,"no start pixel for window option found... "); } win_11: outq = lpntr + xpix1 + 1; /* point to next element */ /* here we either have single value or band */ if (fx[4] > 1.0) { valx = fx[1]; for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { my = ypix1 - ws2; if (my < sublo[1]) my = sublo[1]; ky = my + wsize; if (ky > subhi[1]) ky = subhi[1] + 1; inq = pntri + nrx + (my*npix[0]); /* point to low w_start */ for (nry=my; nry= val) && (*inq <= valx)) { *outq = nry; ypix1 = nry; break; } inq += npix[0]; } outq ++; } } else { for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { my = ypix1 - ws2; if (my < sublo[1]) my = sublo[1]; ky = my + wsize; if (ky > subhi[1]) ky = subhi[1] + 1; inq = pntri + nrx + (my*npix[0]); /* point to low w_start */ for (nry=my; nry= val) { *outq = nry; ypix1 = nry; break; } inq += npix[0]; } outq ++; } } } } else { inp = pntri + (subhi[1]*npix[0]); /* final split for Maximum or Value option */ if (optio[2] == 'M') /* search for max value */ { if (xpix1 < 1) /* no start pixel given */ { fmax = val*0.99; fy = -1; for (nrx=sublo[0]; nrx<=subhi[0]; nrx++) { outq = lpntr + nrx; inq = inp + nrx; for (nry=subhi[1]; nry>= sublo[1]; nry--) { if (*inq > fmax) { fy = nry; fmax = *inq; } if (fy > -1) { *outq = fy; xpix1 = nrx; ypix1 = fy; goto win_2; } inq -= npix[0]; /* move to previous line */ } } SCETER(110,"no start pixel for window option found... "); } win_2: outq = lpntr + xpix1 + 1; /* point to next element */ for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { fmax = val*0.99; fy = -1; my = ypix1 + ws2; if (my > subhi[1]) my = subhi[1]; ky = my - wsize; if (ky < sublo[1]) ky = sublo[1] - 1; inq = pntri + nrx + (my*npix[0]); /* point to hi window start */ for (nry=my; nry>ky; nry--) { if (*inq > fmax) { fy = nry; fmax = *inq; } inq -= npix[0]; } if (fy > -1) { *outq = fy; ypix1 = fy; } outq ++; } } else { if (xpix1 < 1) /* no start pixel given */ { for (nrx=sublo[0]; nrx<=subhi[0]; nrx++) { outq = lpntr + nrx; inq = inp + nrx; for (nry=subhi[1]; nry>= sublo[1]; nry--) { if (*inq >= val) { *outq = nry; xpix1 = nrx; ypix1 = nry; goto win_22; } inq -= npix[0]; /* move to previous line */ } } SCETER(110,"no start pixel for window option found... "); } win_22: outq = lpntr + xpix1 + 1; /* point to next element */ /* here we either have single value or band */ if (fx[4] > 1.0) { valx = fx[1]; for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { my = ypix1 + ws2; if (my > subhi[1]) my = subhi[1]; ky = my - wsize; if (ky < sublo[1]) ky = sublo[1] - 1; inq = pntri + nrx + (my*npix[0]); /* point to hi w_start */ for (nry=my; nry>ky; nry--) { if ((*inq >= val) && (*inq <= valx)) { *outq = nry; ypix1 = nry; break; } inq -= npix[0]; } outq ++; } } else { for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { my = ypix1 + ws2; if (my > subhi[1]) my = subhi[1]; ky = my - wsize; if (ky < sublo[1]) ky = sublo[1] - 1; inq = pntri + nrx + (my*npix[0]); /* point to hi w_start */ for (nry=my; nry>ky; nry--) { if (*inq >= val) { *outq = nry; ypix1 = nry; break; } inq -= npix[0]; } outq ++; } } } } goto common_work; } /* ----------------------------------------- else we work on the dynamic WINDOW option ----------------------------------------- */ if (optio[1] == 'D') { ws2 = wsize / 2; /* and split for BOTTOM or TOP direction */ if (optio[0] == 'B') { inp = pntri + (sublo[1]*npix[0]); if (xpix1 < 1) /* no start pixel given */ { for (nrx=sublo[0]; nrx<=subhi[0]; nrx++) { outq = lpntr + nrx; inq = inp + nrx; for (nry=sublo[1]; nry<= subhi[1]; nry++) { if (*inq >= val) { *outq = nry; xpix1 = nrx; ypix1 = nry; goto win_111; } inq += npix[0]; /* move to next line */ } } SCETER(110,"no start pixel for window option found... "); } win_111: outq = lpntr + xpix1 + 1; /* point to next element */ /* here we either have single value or band */ if (fx[4] > 1.0) { valx = fx[1]; for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { for (ky=0; ky= val) && (*inq <= valx)) { win_1140: *outq = my; ypix1 = my; if (limi == 1) goto win_115; my --; if (my < sublo[1]) goto win_115; inq -= npix[0]; if ((*inq >= val) && (*inq <= valx)) goto win_1140; else goto win_115; } my = ypix1 + ky; if (my > subhi[1]) my = subhi[1]; inq = pntri + nrx + (my*npix[0]); /* move to upper end */ if ((*inq >= val) && (*inq <= valx)) { *outq = my; ypix1 = my; goto win_115; } } win_115: outq ++; } } else { for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { for (ky=0; ky= val) { win_1150: *outq = my; ypix1 = my; if (limi == 1) goto win_116; my --; if (my < sublo[1]) goto win_116; inq -= npix[0]; if (*inq >= val) goto win_1150; else goto win_116; } my = ypix1 + ky; if (my > subhi[1]) my = subhi[1]; inq = pntri + nrx + (my*npix[0]); /* move to upper end */ if (*inq >= val) { *outq = my; ypix1 = my; goto win_116; } } win_116: outq ++; } } } else /* first look up then down */ { inp = pntri + (subhi[1]*npix[0]); if (xpix1 < 1) /* no start pixel given */ { for (nrx=sublo[0]; nrx<=subhi[0]; nrx++) { outq = lpntr + nrx; inq = inp + nrx; for (nry=subhi[1]; nry>= sublo[1]; nry--) { if (*inq >= val) { *outq = nry; xpix1 = nrx; ypix1 = nry; goto win_222; } inq -= npix[0]; /* move to previous line */ } } SCETER(110,"no start pixel for window option found... "); } win_222: outq = lpntr + xpix1 + 1; /* point to next element */ /* here we either have single value or band */ if (fx[4] > 1.0) { valx = fx[1]; for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { for (ky=0; ky subhi[1]) { limi = 1; my = subhi[1]; } else limi = 0; inq = pntri + nrx + (my*npix[0]); /* move to upper end */ if ((*inq >= val) && (*inq <= valx)) { win_2240: *outq = my; ypix1 = my; if (limi == 1) goto win_225; my++ ; if (my > subhi[1]) goto win_225; inq += npix[0]; if ((*inq >= val) && (*inq <= valx)) goto win_2240; else goto win_225; } my = ypix1 - ky; if (my < sublo[1]) my = sublo[1]; inq = pntri + nrx + (my*npix[0]); /* move to lower end */ if ((*inq >= val) && (*inq <= valx)) { *outq = my; ypix1 = my; goto win_225; } } win_225: outq ++; } } else { for (nrx=xpix1+1;nrx<=subhi[0]; nrx++) { for (ky=0; ky subhi[1]) { limi = 1; my = subhi[1]; } else limi = 0; inq = pntri + nrx + (my*npix[0]); /* move to upper end */ if (*inq >= val) { win_2250: *outq = my; ypix1 = my; if (limi == 1) goto win_226; my++ ; if (my > subhi[1]) goto win_226; inq += npix[0]; if (*inq >= val) goto win_2250; else goto win_226; } my = ypix1 - ky; if (my < sublo[1]) my = sublo[1]; inq = pntri + nrx + (my*npix[0]); /* move to loer end */ if (*inq >= val) { *outq = my; ypix1 = my; goto win_226; } } win_226: outq ++; } } } goto common_work; } /* --------------------------------- else we work on the GLOBAL option --------------------------------- */ inq = pntri + (sublo[1]*npix[0]); inp = inq; /* and split for BOTTOM or TOP direction */ if (optio[0] == 'B') { /* split for Maximum or Value option */ if (optio[2] == 'M') /* search for max value */ { inp = pntri + (sublo[1]*npix[0]); /* point to bottom line */ for (nrx=sublo[0]; nrx fmax) { fy = nry; fmax = *inq; } inq += npix[0]; } if (fy > -1) *outq = fy; } } else { if (fx[4] > 1.0) { valx = fx[1]; for (nry=sublo[1]; nry<=subhi[1]; nry++) { inq = inp + sublo[0]; for (nrx=sublo[0]; nrx= val) && (*inq <= valx)) { outq = lpntr + nrx; if (*outq < 0.0) *outq = nry; /* store lowest y-pos */ } inq ++; } inp += npix[0]; } } else { for (nry=sublo[1]; nry<=subhi[1]; nry++) { inq = inp + sublo[0]; for (nrx=sublo[0]; nrx= val) { outq = lpntr + nrx; if (*outq < 0.0) *outq = nry; /* store lowest y-pos */ } inq ++; } inp += npix[0]; } } } } else { /* split for Maximum or Value option */ if (optio[2] == 'M') /* search for max value */ { inp = pntri + (subhi[1]*npix[0]); /* point to top line */ inq = inp + sublo[0]; outq = lpntr + sublo[0]; for (nrx=sublo[0]; nrx=sublo[1]; nry--) { if (*inq > fmax) { fy = nry; fmax = *inq; } inq -= npix[0]; } if (fy > -1) *outq = fy; inp ++; outq ++; } } else { if (fx[4] > 1.0) { valx = fx[1]; for (nry=sublo[1]; nry<=subhi[1]; nry++) { inq = inp + sublo[0]; for (nrx=sublo[0]; nrx= val) && (*inq <= valx)) { outq = lpntr + nrx; *outq = nry; /* always store highest y-position */ } inq ++; } inp += npix[0]; } } else { for (nry=sublo[1]; nry<=subhi[1]; nry++) { inq = inp + sublo[0]; for (nrx=sublo[0]; nrx= val) { outq = lpntr + nrx; *outq = nry; /* always store highest y-position */ } inq ++; } inp += npix[0]; } } } } common_work: kstart = -1; /* get first pixel above threshold */ outq = pntro; for (nrx=0; nrx (-1.0)) { kstart = nrx; break; } outq ++; } if (kstart < 0) return (-1); outq = pntro + npix[0] - 1; for (nrx=npix[0]; nrx>kstart; nrx--) { if (*outq > (-1.0)) { end = nrx-1; break; } outq --; } if (kstart == end) { *pntro = *(pntro+kstart); return (1); /* just 1 pixel... */ } count = end - kstart + 1; /* no. of pixels */ /* move to front, if necessary */ if (kstart > 0) { inq = pntro; outq = pntro + kstart; for (nrx=0; nrx kmax) kmax = kk; kk = 0; } inq ++; } if (kk > kmax) kmax = kk; /* don't forget last */ if (mm > 0) /* we must have at least 3 pixels */ { (void) sprintf(cbuf,"%d holes of max. %d size detected ...",mm,kmax); SCTPUT(cbuf); mm = 0; inq = pntro + 1; while(inq < (pntro+count)) { if (*inq < 0.0) { outq = inq + 1; /* inq-1 is != -1 ... */ if (*outq < 0.0) { loop: inq ++ ; outq = inq + 1; /* inq-1 is != -1 ... */ if (*outq < 0.0) goto loop; else goto next; } *inq = (*(inq-1)+*outq)/2; /* interpolate */ } next: inq ++; } } return (count); }