/*----------------------------------------------------------------------PROGRAM: modelwarp_mc.cDATE:    6/2/94AUTHOR:  Baback Moghaddam, baback@media.mit.edu------------------------------------------------------------------------   This routine reads in an ASCII file of the format   filename  x1 x2 x3 x4 y1 y2 y3 y4   .   .   .    and warps the images according to the affine transform   between the feature locations indicated and those in    the model.bf file   It then masks the the warped files and normalized its contrast   it the writes the results to output directory with the same names   ---------------------------------------------------------------------- */#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>#include <float.h>#include "util.h"#include "io.h"#include "matrix.h"#include "affine.h"#define MAX_CHARS  512/* ----------- Command-Line Parsing Stuff ------- */extern int optind;extern char *optarg;char *progname;                /* used to store the name of the program  */char comline[MAX_CHARS];       /* used to store the entire command line  */#define OPTIONS "i:l:m:o:x:y:ra:b:w:p:s:zfv:hc:d:e:"char *usage = "\t-i indir -l list -m model.bf -o outdir\n\\t\t\t[-r] [-a pt_1] [-b pt_2] [-x out_xdim] [-y out_ydim]\n\\t\t\t[-w mask.bf] [-p pedestal] [-s scale] [-f] [-z]\n\\t\t\t[-v bgvalue] [-h] [-c pt_1] [-d pt_2] [-e pt_3]\n";char *help="\Correspondance-Based Affine Warping (+ Masking & Contrast)\n\n\-i indir    \t input directory\n\-l list     \t ASCII list of filenames and their features locations\n\            \t with the following line format:\n\n\            \t   filename  i1 i2 i3 ... iN j1 j2 j3 ... jN\n\n\-m model.bf \t BF-format file of reference point coordinates\n\            \t with the following 2-by-N matrix format:\n\n\            \t   i1 i2 i3 ... iN\n\            \t   j1 j2 j3 ... jN\n\n\-o outdir   \t output directory\n\-r          \t apply rigid transform (based on default correspondances)\n\-a rigid_1  \t rigid transform feature #1 (default = 1)\n\-b rigid_2  \t rigid transform feature #2 (default = 2)\n\-x out_xdim \t output x dimension (1st n columns) of warped image\n\-y out_ydim \t output y dimension (1st n rows)    of warped image\n\-w mask.bf  \t mask file\n\-p pedestal \t pedestal for contrast normalization\n\-s scale    \t scale    for contrast normalization\n\-f          \t write float output image (default is same as input)\n\-z          \t contrast statistics computed over non-zero part of image\n\-v          \t pixel value in void regions (default = 0)\n\-h          \t apply selective 3pt affine warp\n\-c pt_1     \t selective affine warp feature #1 (default = 1)\n\-d pt_2     \t selective affine warp feature #2 (default = 2)\n\-e pt_3     \t selective affine warp feature #3 (default = 3)\n";/* --- command-line defaults --- */int do_mask      = 0;int do_contrast  = 0;int contrast_nonzero_flag   = 0;int out_float_flag = -1;  /* -1 = unspecified *//* ---- function prototypes ----- */int contrast(float **image_in, float **image_out, 	     int nrow, int ncol,	     float pedestal, float scale,	     float *mean, float *sigma);/*----------------------------------------------------------------------*/main(int argc, char **argv){   register int i,j,k,l,ii,jj;  int f, c, nframe, nfeatures, sets, bytes_pixel;  int nrow, ncol, out_nrow, out_ncol, M, N;  char command[MAX_CHARS],indir[MAX_CHARS],listfile[MAX_CHARS], \    infile[MAX_CHARS], modelfile[MAX_CHARS],outdir[MAX_CHARS], \      filename[MAX_CHARS], maskfile[MAX_CHARS], line[MAX_CHARS];  float **P, **Q, **Ps, **Qs, **W, **Winv;   float **image_in, **image_out, **mask;  unsigned char **uchar_image;  float contrast_mean, contrast_sigma;  int   contrast_npixels;  float fval1, fval2, fval;  FILE *fp, *fp2;         /* for output values dump */  /* required input flags */    int errflag   = 0;  int inflag    = 0;  int listflag  = 0;  int modelflag = 0;  int outflag   = 0;  int y_flag    = 0;  int x_flag    = 0;    /*--- command line defaults ---- */  int rigid_select = 0;  int rigid_pt_1 = 1;  int rigid_pt_2 = 2;  int affine_select = 0;  int affine_select_pt[8] = {0, 1, 2, 3, 4, 5, 6, 7};    float contrast_pedestal = -1;  /* unspecified values */  float contrast_scale    = -1;  float bgvalue = 0.0;  progname = argv[0];  for (i=0; i<argc; i++)    strcat(comline, argv[i]),strcat(comline, " ");      /* ----------------------  Command Line Parse ------------------------ */    while ((c = getopt(argc, argv, OPTIONS)) != EOF)    switch (c) {          case 'i':      strcpy(indir, optarg);      inflag = 1;      break;    case 'l':      strcpy(listfile, optarg);      listflag = 1;      break;          case 'm':      strcpy(modelfile, optarg);      modelflag = 1;      break;    case 'o':      strcpy(outdir, optarg);      outflag = 1;      break;    case 'y':      out_nrow = atoi(optarg);      y_flag = 1;      break;    case 'x':      out_ncol = atoi(optarg);      x_flag = 1;      break;    case 'r':      rigid_select = 1;      break;    case 'a':      rigid_pt_1 = atoi(optarg);      rigid_select = 1;      break;    case 'b':      rigid_pt_2 = atoi(optarg);      rigid_select = 1;      break;    case 'w':      strcpy(maskfile, optarg);      do_mask = 1;      break;    case 'p':      contrast_pedestal = atof(optarg);      break;    case 's':      contrast_scale = atof(optarg);      break;    case 'z':      contrast_nonzero_flag = 1;      break;    case 'f':      out_float_flag = 1;      break;    case 'v':      bgvalue = atof(optarg);      break;    case 'h':      affine_select = 1;      break;    case 'c':      affine_select_pt[1] = atoi(optarg);      break;    case 'd':      affine_select_pt[2] = atoi(optarg);      break;    case 'e':      affine_select_pt[3] = atoi(optarg);      break;    case '?':      errflag = 1;      break;    }    /* command line error check */    if (errflag || !inflag || !listflag || !modelflag) {    fprintf(stderr,"\nUSAGE: %s %s\n%s\n", progname, usage, help);    exit(1);  }  if (contrast_pedestal>0 & contrast_scale>0)    do_contrast = 1;  /* ----  read indir descriptor file -------- */    read_descriptor(indir, &nframe, &sets, &bytes_pixel, &ncol, &nrow);  if (sets>1)     myerror("Input files must be single-set DAT files!");  if (y_flag==0)     out_nrow = nrow;  if (x_flag==0)    out_ncol = ncol;  if (out_float_flag == -1)    out_float_flag = (bytes_pixel==4 ? 1:0);  /* ---- read model file ---- */  P = read_BIN(modelfile, &M, &nfeatures);  if (affine_select==1 && nfeatures<3) {    fprintf(stderr,	    "ERROR: selective affine warp needs atleast a 3pt-model\n");    exit(3);  }  if (affine_select==1) {    for (i=1; i<=3; i++)      if ((affine_select_pt[i] > nfeatures) || (affine_select_pt[i]<1)) {	fprintf(stderr,		"ERROR: affine_select_pt[%d]=%d is invalid for nfeatures=%d\n",		i, affine_select_pt[i], nfeatures);	exit(4);      }  }  Q = matrix(1, M, 1, nfeatures);  if (affine_select==1) {    Ps = matrix(1, 2, 1, 3);    Qs = matrix(1, 2, 1, 3);  }  W    = matrix(1, 3, 1, 3);  Winv = matrix(1, 3, 1, 3);  image_in   = matrix(1, nrow, 1, ncol);  image_out  = matrix(1, nrow, 1, ncol);  if (bytes_pixel == 1)    uchar_image = cmatrix(1, out_nrow, 1, out_ncol);  printf("\nNo. of model features = %d\n\n", nfeatures);  for (i=1; i<=M; i++) {    for (j=1; j<=nfeatures; j++)      printf("%3d ", (int) P[i][j]);    printf("\n");  }  printf("\n");     /* ---- read mask file (if called for) ---- */  if (do_mask) {    int m,n;    mask = read_BIN(maskfile, &m, &n);    if (m!=out_ncol || n!=out_nrow) {      fprintf(stderr,	      "ERROR: Mask file dimensions %d-by-%d are invalid\n",m,n);      exit(10);    }  }  /* ---- loop over input list file and warp ------- */  if ((fp = fopen(listfile, "r")) == NULL) {    fprintf(stderr,"ERROR: Could not open input file %s \n\n", listfile);    exit(1);  }  nframe = 0;  while (strlen(fgets(line, MAX_CHARS, fp))) {    if (strncmp(line, "#", 1) != 0 && strlen(line)>1) {      nframe++;          /* --- read filename and reference coordinates --- */            strcpy(infile, strtok(line, " \t"));      for (i=1; i<=nfeatures; i++)	Q[1][i] = atof(strtok(NULL, " \t"));      for (i=1; i<=nfeatures; i++)	Q[2][i] = atof(strtok(NULL, " \t"));            fprintf(stdout, "%s  ", infile);      for (i=1; i<=nfeatures; i++)	fprintf(stdout,"%3d ",  (int) Q[1][i]);      for (i=1; i<=nfeatures; i++)	fprintf(stdout,"%3d ",  (int) Q[2][i]);      fprintf(stdout,"\n");            fprintf(stdout,"Affine matrix: \n");      if (rigid_select == 1)	rigid(Q, P, nfeatures, W, rigid_pt_1, rigid_pt_2);      else if (affine_select == 0)	affine(Q, P, nfeatures, W);      else {	for (i=1; i<=3; i++) {	  Ps[1][i] = P[1][affine_select_pt[i]];	  Ps[2][i] = P[2][affine_select_pt[i]];	  Qs[1][i] = Q[1][affine_select_pt[i]];	  Qs[2][i] = Q[2][affine_select_pt[i]];	}	affine(Qs, Ps, 3, W);      }      for (i=1; i<=3; i++) {	for (j=1; j<=3; j++)	  fprintf(stdout,"%+3.1f ", W[i][j]);	fprintf(stdout,"\n");      }      matrix_inverse(W, 3, Winv);                  /* --- read file  --- */            sprintf(filename,"%s/%s", indir, infile);      if (bytes_pixel == 4)	read_RAW_float(filename, image_in, nrow, ncol);      else	read_RAW(filename, image_in, nrow, ncol);      fprintf(stdout, "Read %d-by-%d %s file %s\n",	      nrow, ncol, (bytes_pixel==4 ? "float":"uchar"), infile);      fprintf(stdout,"Affine matrix: \n");        for (i=1; i<=3; i++) {	for (j=1; j<=3; j++)	  fprintf(stdout,"%+3.3f ", W[i][j]);	fprintf(stdout,"\n");      }                /* --- apply warp to image --- */            printf("bgvalue = %f\n", bgvalue);       i = affine_warp(image_in, image_out, nrow, ncol, W, bgvalue);      fprintf(stdout,"Applied affine warp ... %d lowpass iterations\n", i);      if (bytes_pixel == 1)	for (i=1; i<=out_nrow; i++)	  for (j=1; j<=out_ncol; j++)	    image_out[i][j] = (int) image_out[i][j];      /* --- apply mask --- */      if (do_mask) {	if (bytes_pixel == 1) {	  for (i=1; i<=out_nrow; i++)	    for (j=1; j<=out_ncol; j++) 	      image_out[i][j] = (int) (image_out[i][j] * mask[i][j]);	  	}	if (bytes_pixel == 4) {	  for (i=1; i<=out_nrow; i++)	    for (j=1; j<=out_ncol; j++)	      image_out[i][j] *= mask[i][j];	}	fprintf(stdout,"Applied mask ...\n");      }            /* --- normalize contrast ---- */            if (do_contrast) {	contrast_npixels = contrast(image_out, image_in, 				    out_nrow, out_ncol,				    contrast_pedestal, contrast_scale,				    &contrast_mean, &contrast_sigma);	fprintf(stdout,		"Contrast pixels = %d\nContrast mean   = %3.2f\nContrast sigma  = %3.2f\n",		contrast_npixels, contrast_mean, contrast_sigma);	for (i=1; i<=out_nrow; i++)	  for (j=1; j<=out_ncol; j++)	   image_out[i][j] = image_in[i][j];      }      /* --- write output image --- */          if (out_float_flag == 0) {	for (i=1; i<=out_nrow; i++) 	  for (j=1; j<=out_ncol; j++)	    if (image_out[i][j]>255)	      uchar_image[i][j] = 255;	    else if (image_out[i][j]<0)	      uchar_image[i][j] = 0;	    else 	      uchar_image[i][j] = (int) (image_out[i][j] + 0.5);      }      sprintf(filename,"%s/%s", outdir, infile);      if (out_float_flag == 1)	write_RAW_float(filename, image_out, out_nrow, out_ncol);      else	write_RAW(filename, uchar_image, out_nrow, out_ncol);      fprintf(stdout,"Wrote %d-by-%d %s output file %s\n\n", 	      out_nrow, out_ncol, 	      (out_float_flag==1 ? "float":"uchar"), filename);                }  }  /* --- write output dir descriptor --- */    write_descriptor(outdir, nframe, out_ncol, out_nrow, bytes_pixel, comline);  free_matrix(P, 1, M, 1, nfeatures);  free_matrix(Q, 1, M, 1, nfeatures);  free_matrix(W, 1, 3, 1, 3);  free_matrix(Winv, 1, 3, 1, 3);  free_matrix(image_in, 1, nrow, 1, ncol);  free_matrix(image_out, 1, nrow, 1, ncol);  if (bytes_pixel == 1)    free_cmatrix(uchar_image, 1, out_nrow, 1, out_ncol);  if (do_mask == 1)    free_matrix(mask, 1, out_nrow, 1, out_ncol);  if (affine_select) {    free_matrix(Ps, 1, 2, 1, 3);    free_matrix(Qs, 1, 2, 1, 3);  }  fclose(fp);  return 0;}/* ====================================================================== */int contrast(float **image_in, float **image_out, 	      int nrow, int ncol,	      float pedestal, float scale,	      float *mean, float *sigma){  register int i,j;  float sum;  float my_mean, my_sigma;  int count;  /* --- compute statistics ---- */  sum = count = 0;  for (i=1; i<=nrow; i++)    for (j=1; j<=ncol; j++) {      if (contrast_nonzero_flag == 1) {	if (image_in[i][j] != 0) {	  sum += image_in[i][j];	  count ++;	}      }      else {	sum += image_in[i][j];	count ++;      }    }  my_mean = sum/count;  sum = count = 0;  for (i=1; i<=nrow; i++)    for (j=1; j<=ncol; j++) {      if (contrast_nonzero_flag == 1) {	if (image_in[i][j] != 0) {	  sum += SQR(image_in[i][j] - my_mean);	  count ++;	}      }      else {	sum += SQR(image_in[i][j] - my_mean);	count ++;      }    }  my_sigma = sqrt(sum/(count-1));  /* --- remap input image ---- */  for (i=1; i<=nrow; i++)    for (j=1; j<=ncol; j++)       if (contrast_nonzero_flag == 1) {	if (image_in[i][j] != 0)	  image_out[i][j] = pedestal + 	    scale * (image_in[i][j]-my_mean)/my_sigma;	else	  image_out[i][j] = 0;      }  *mean = my_mean;  *sigma = my_sigma;  return count;}