/* @(#)trace2d.c 17.1.1.1 (ESO-DMD) 01/25/02 17:39:38 */ /*=========================================================================== 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 ===========================================================================*/ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ .COPYRIGHT: Copyright (c) 1994 European Southern Observatory, all rights reserved .IDENTIFIER TRACE1D .LANGUAGE C .AUTHOR K. Banse IPG-ESO Garching .KEYWORDS Image Display, graphic trace .PURPOSE extract a rotated subimage from a displayed image .ALGORITHM use two cursors to define parallel lines .INPUT/OUTPUT call as stat = TRACE2D( imf, stepln, outfra ) input: int imf : id of the input image double stepln : step size along the trace char *outfra : name of the trace .RETURNS status: 0 = ok 1 = invalid corners .ENVIRONment MIDAS #include Prototypes for MIDAS interfaces #include Global variables for DISPLAY interfaces .VERSIONS 1.00 940630 from TRACE.FOR R.M.van Hees 010423 last modif ------------------------------------------------------------*/ /* Define _POSIX_SOURCE to indicate that this is a POSIX program */ #define _POSIX_SOURCE 1 #include #include #include #include #include #ifndef TRUE #define TRUE 1 #define FALSE 0 #endif #ifndef PI #define PI 3.14159265358979325e0 #endif #define HALFPI (PI / 2) #define D_LARGE 1.0e+15 #define CFORM 8 /* cursor shape (= x-hair) */ #define COLOR 2 /* cursor color */ #define POLL_POS 2 /* continuous cursor position */ #define ENTR_POS 1 /* cursor position after ENTER */ static int ffelem = 1; /* */ /*++++++++++++++++++++++++++++++ .IDENTIFIER SWAP_RECT .PURPOSE Get the point of a rectangle in the right order .INPUT/OUTPUT input: double rm : slope of the first line in/output: double *beta : angle of the second line struct Point rect[4]: points of the rectangle .RETURNS nothing .COMMENTS static function --------------------------------*/ #ifdef __STDC__ static void SWAP_RECT( double rm, double *beta, struct Point *rect ) #else static void SWAP_RECT( rm, beta, rect ) double rm, *beta; struct Point *rect; #endif { if ( rm < 0.0 ) { (*beta) += HALFPI; if ( rect[0].px < rect[3].px ) { if ( rect[0].px > rect[1].px ) { SWAP_PNT( 0, 1, rect ); SWAP_PNT( 2, 3, rect ); } SWAP_PNT( 1, 3, rect ); } else { if ( rect[0].px < rect[1].px ) { SWAP_PNT( 0, 1, rect ); SWAP_PNT( 2, 3, rect ); } SWAP_PNT( 0, 2, rect ); } } else { if ( rect[0].px < rect[3].px ) { if ( rect[0].px > rect[1].px ) { SWAP_PNT( 0, 1, rect ); SWAP_PNT( 2, 3, rect ); } } else { if ( rect[0].px < rect[1].px ) { SWAP_PNT( 0, 3, rect ); SWAP_PNT( 1, 2, rect ); } else { SWAP_PNT( 0, 2, rect ); SWAP_PNT( 1, 3, rect ); } } } } /* */ /*++++++++++++++++++++++++++++++ .IDENTIFIER CROSS_PNT .PURPOSE determine cross point (x,y) .ALGORITHM use (y - y1) / (x - x1) = rrm (y - y2) / (x - x2) = rm and set them equal .INPUT/OUTPUT input: double rm1 : slope of the first line struct Point *pnt1 : point on the first line double rm2 : slope of the second line struct Point *pnt2 : point on the second line output: struct Point *cross : point where both lines intersect .RETURNS nothing .COMMENTS static function --------------------------------*/ #ifdef __STDC__ static void CROSS_PNT( double rm1, const struct Point *pnt1, double rm2, const struct Point *pnt2, struct Point *cross ) #else static void CROSS_PNT( rm1, pnt1, rm2, pnt2, cross ) double rm1, rm2; struct Point *pnt1, *pnt2, *cross; #endif { double dummy; dummy = pnt1->px * rm2 - pnt2->px * rm1 + pnt2->py - pnt1->py; cross->px = CGN_NINT( dummy / (rm2 - rm1) ); cross->py = CGN_NINT( rm1 * (cross->px - pnt2->px) + pnt2->py); } /* */ int TRACE2D( imf, stepln, outfra ) char *outfra; int imf; double *stepln; { register int nn; int actvals, enter, imno, iter, knul, unit, idum[2], ipix[2], npix[2], kk, disp[4], xfig[5], yfig[5]; float ca, sa, *pntrX, *pntrY, *pntrS, *p_img, cuts[4], image[4], xa, ya, xb, yb, pxend[4], pixels[6]; double angle, beta, dx, dy, rm, rrm, pxstep, step[2], startS[2], stepS[2], dd1[4], dd2[4], dd3[4]; /* synchronize with fp2wc */ char *cpntr, actio, ident[33], cunit[49], output[81]; struct Point lpnts[1], rect[4]; /* initialized variables */ int first = TRUE; static int bright, draw, erase_flag, mxpnts, old_disp[4]; static float nulval; static struct Point cpos[2]; static char *info_mes1 = "Use both cursors and press enter to fix baseline of rectangle.", *info_mes2 = "Move cursor and press enter to determine the line parallel to the baseline", *info_mes3 = "and the one closing side of the rectangle", *info_mes4 = "Move cursor again and press enter to close the rectangle", *info_mes5 = "Finally we extract the data inside the rotated rectangle"; if ( first ) { (void) SCKRDR( "NULL", 1, 1, &actvals, &nulval, &unit, &knul ); mxpnts = 10000; /* flags for IDI interfaces */ if ( IDINUM == 11 ) { erase_flag = 99; draw = 99; } else { erase_flag = 0; draw = 200; } bright = 255; } /* setup the cursors */ CURS_SETUP(3,CFORM,COLOR); actio = '2'; /* -> EXTRACT/ROTATED ... */ /* read cursors continuously + update graph accordingly */ iter = 1; (void) SCTDIS( info_mes1, 80 ); do { if (GET_CPOS(0,&actio,POLL_POS,3,&iter,&enter,cpos) != 0) return 0; /* EXIT */ if ( cpos->px > cpos[1].px ) SWAP_PNT( 0, 1, cpos ); /* set the angles */ dx = cpos[1].px - cpos->px; dy = cpos[1].py - cpos->py; rm = dy / dx; /* calculate end points of line through CPOS */ END_PNT( rm, cpos, disp ); /* erase old vector + draw new vector */ nn = 0; while ( nn < 4 && disp[nn] == old_disp[nn] ) nn++; if ( nn != 4 ) { if ( ! first ) /* erase old vector + draw new one */ (void) IIGPLY_C( QDSPNO, QOVCH, old_disp, old_disp+2, 2, erase_flag, 1 ); else { first = FALSE; (void) IIMCMY_C( QDSPNO, &QOVCH, 1, 0 ); } (void) IIGPLY_C( QDSPNO, QOVCH, disp, disp+2, 2, draw, 1 ); /* store new vector into old one */ for ( nn = 0; nn < 4; nn++ ) old_disp[nn] = disp[nn]; } } while ( ! enter ); angle = atan2( dy, dx ); if ( IDINUM == 11 ) { (void) IIGPLY_C( QDSPNO, QOVCH, disp, disp+2, 2, erase_flag, 1 ); (void) IIGPLY_C( QDSPNO, QOVCH, disp, disp+2, 2, bright, 1 ); } /* now continue with cursor 0 */ if ( IDINUM == 11 ) (void) strcpy(output,"Now use only the mouse to continue:"); else (void) strcpy(output,"Now use only cursor 0 to continue:"); (void) SCTDIS(output, 0); (void) SCTDIS(info_mes2,0); (void) SCTDIS(info_mes3,0); iter = 1; if (GET_CPOS(0,&actio,ENTR_POS,0,&iter,idum,cpos) != 0) return 0; /* EXIT */ /* first point of the rectangle */ rect[1].px = cpos->px; rect[1].py = cpos->py; /* Determine line at 90 degrees angle */ beta = angle + HALFPI; rrm = tan( beta ); lpnts->px = disp[0]; /* a point on the previous line */ lpnts->py = disp[2]; CROSS_PNT( rm, cpos, rrm, lpnts, rect ); /* rect[0] is a intersection */ END_PNT( rrm, cpos, disp ); /* the whole line */ (void) IIGPLY_C( QDSPNO, QOVCH, disp, disp+2, 2, bright, 1 ); /* draw it */ /* Now get parallel line */ END_PNT( rm, cpos, disp ); (void) IIGPLY_C( QDSPNO, QOVCH, disp, disp+2, 2, bright, 1 ); /* draw it */ /* Get final cursor to close rectangle */ (void) SCTDIS( info_mes4, 0 ); iter = 1; if (GET_CPOS(0,&actio,ENTR_POS,0,&iter,idum,cpos) != 0) return 0; /* EXIT */ END_PNT( rrm, cpos, disp ); (void) IIGPLY_C( QDSPNO, QOVCH, disp, disp+2, 2, bright, 1 ); /* draw it */ /* And determine cross points with the two lines */ CROSS_PNT( rm, cpos, rrm, rect, rect+3 ); CROSS_PNT( rm, cpos, rrm, rect+1, rect+2 ); /* clear all graphs + show only rectangle */ (void) SCTDIS( info_mes5, 0 ); IIMCMY_C( QDSPNO, &QOVCH, 1, 0 ); /* but first put corner points into right order */ SWAP_RECT( rm, &beta, rect ); for ( nn = 0; nn < 4; nn++ ) { xfig[nn] = rect[nn].px; yfig[nn] = rect[nn].py; } xfig[4] = rect->px; yfig[4] = rect->py; (void) IIGPLY_C( QDSPNO, QOVCH, xfig, yfig, 5, bright, 1 ); /* Read the input image descriptor */ (void) SCDRDI(imf,"NPIX",1,2,&actvals,npix,&unit,&knul); (void) SCDRDD(imf,"STEP",1,2,&actvals,step,&unit,&knul); (void) SCDRDR(imf,"LHCUTS",1,2,&actvals,cuts,&unit,&knul); (void) SCDGETC(imf,"IDENT",1,32,&actvals,ident); (void) SCDGETC(imf,"CUNIT",1,48,&actvals,cunit); if (ZPLANE != 0) /* 3-dim frame was loaded */ { dd1[2] = ZPLANE; ffelem = ((ZPLANE-1)*npix[0]*npix[1]) + 1; } else ffelem = 1; if (Pixconv("INIT",imf,dd1,dd2,dd3) > 0) /* init wc conversion */ SCETER(69,"initialization of world coord. conversion failed ..."); /* now move from screen space to image channel space */ for (nn=0; nn<4; nn++) Sc2ch(1,xfig+nn,yfig+nn); /* now convert the corner points from screen pixels to real pixel no's + world coordinates */ dd1[0] = (double)xfig[1]; dd1[1] = (double)yfig[1]; if ( Pixconv("IRW",imf,dd1,dd2,dd3) != 0 ) { SCTPUT( "TRACE2D: error in start pixel conversion..." ); return 1; /* invalid corners, EXIT */ } else { image[0] = (float)dd2[0]; image[2] = (float)dd2[1]; startS[0] = dd3[0]; startS[1] = dd3[1]; } dd1[0] = (double)xfig[2]; dd1[1] = (double)yfig[2]; if ( Pixconv("IRW",imf,dd1,dd2,dd3) != 0 ) { SCTPUT("TRACE2D: error in end pixel conversion..."); return 1; /* invalid corners, EXIT */ } else { pxend[0] = (float)dd2[0]; pxend[2] = (float)dd2[1]; } dd1[0] = (double)*xfig; dd1[1] = (double)*yfig; if ( Pixconv("IRW",imf,dd1,dd2,dd3) != 0 ) { SCTPUT("TRACE2D: error in corner pixel conversion..."); return 1; /* invalid corners, EXIT */ } else { image[1] = (float)dd2[0]; image[3] = (float)dd2[1]; } pixels[0] = (double)xfig[3]; pixels[1] = (double)yfig[3]; if ( Pixconv("IRW",imf,dd1,dd2,dd3) != 0 ) { SCTPUT("TRACE2D: error in corner pixel conversion..."); return 1; /* invalid corners, EXIT */ } else { pxend[1] = (float)dd2[0]; pxend[3] = (float)dd2[1]; } /* pull out lines from area inside rectangle */ if ( fabs( stepln[1] ) <= 1.0e-9 ) { stepS[1] = step[1]; pxstep = 1.0; } else { stepS[1] = stepln[1]; pxstep = stepS[1] / step[1]; } pntrX = (float *) malloc((unsigned int)(mxpnts * sizeof(float))); pntrY = (float *) malloc((unsigned int)(mxpnts * sizeof(float))); xa = image[0]; ya = image[2]; xb = image[1]; yb = image[3]; ipix[1] = Cpixlin(xa,ya,xb,yb,pxstep,mxpnts,pntrX,pntrY); /* now extract the indices of the base line from the frame */ if ( fabs( *stepln ) <= 1.0e-9 ) { *stepS = *step; pxstep = 1.0; } else { *stepS = *stepln; pxstep = *stepS / *step; } /* pull out max. IPIXL pixels per line along base line */ image[1] = pxend[0]; image[3] = pxend[2]; xb = image[1]; yb = image[3]; ipix[0] = Cpixlin(xa,ya,xb,yb,pxstep,mxpnts,pntrX,pntrY); /* * create output file + get z-values for pixels recorded in XINDX + YINDX */ (void) SCIPUT(outfra,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,2,ipix, startS,stepS,ident,cunit,&cpntr,&imno); pntrS = (float *) cpntr; /* and pull out the z-values */ ca = (float) cos( beta ); sa = (float) sin( beta ); kk = npix[0]*npix[1]; p_img = (float *) malloc((unsigned int)(kk * sizeof(float))); (void) SCFGET(imf,ffelem,kk,&actvals,(char *) p_img); for (nn=0; nn