simcirc.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpccsv.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Draw circle in x,y matrix with 1 inside circle and 0 outside of it",
  "for simple simulations, SW testing, and creating ROI masks.",
  "Matrix size and circle diameter must be given in pixels.",
  " ",
  "Usage: @P [Options] dimension diameter outputfile",
  " ",
  "Options:",
  " -format=<<TSV>|<CSV>|<float>|<short>>",
  "     Matrix data can be saved as tab separated values (TSV, default),",
  "     comma separated values (CSV), or as binary flat files with either",
  "     floats or short ints.",
  " -diam2=<2nd diameter>",
  "     If the second diameter is specified, then a ring is created.",
  " -iring=<value inside hole>",
  "     With option -diam2, the integer value for hole inside ring;",
  "     by default 0.",
  " -segm=<Number of segments>",
  "     Divide circle or ring into given number of segments with increasing",
  "     value in each.",
  " -inv",
  "     Write 0 inside the circle or ring, and 1 outside of it.",
  "     Not applicable with options -iring and -segm.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example #1:",
  "  @P 256 224 circle.txt",
  "Example #2:",
  "  @P -format=float 256 64 circle.bin",
  "  flat2nii circle.bin 256 256 1 1 circle",
  " ",
  "See also: flat2nii, asc2flat, simiart, simboxes, simelli, imgadd, imglkup",
  " ",
  "Keywords: simulation, image, mask, software testing",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int       ai, help=0, version=0, verbose=1;
  char     *cptr, outfile[FILENAME_MAX];
  int       dim=-1;
  int       csv_separator=2; // 0=comma, 1=semi-colon, 2=tab, 3=space
  int       output_format=0; // 0=ASCII, 1=float binary, 2=short int binary
  double    diameter=-1, diameter2=-1;
  int       invert=0; // 0 or 1
  int       segments=1; // divide circle or ring into nr of segments
  int       iring=0; // value inside ring hole
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  outfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "DIAM2=", 6)==0) {
      diameter2=atofVerified(cptr+6); 
      if(!isnan(diameter2) && diameter2>0.0) continue;
    } else if(strncasecmp(cptr, "SEGM=", 5)==0) {
      if(atoiCheck(cptr+5, &segments)==0 && segments>0) continue;
    } else if(strncasecmp(cptr, "IRING=", 6)==0) {
      if(atoiCheck(cptr+6, &iring)==0 && iring>=0) continue;
    } else if(strncasecmp(cptr, "F=", 2)==0) {
      cptr+=2;
      if(!strcasecmp(cptr, "TSV")) {csv_separator=2; output_format=0; continue;}
      if(!strcasecmp(cptr, "CSV")) {csv_separator=0; output_format=0; continue;}
      if(!strcasecmp(cptr, "FLOAT")) {output_format=1; continue;}
      if(!strcasecmp(cptr, "SHORT")) {output_format=2; continue;}
    } else if(strncasecmp(cptr, "FORMAT=", 7)==0) {
      cptr+=7;
      if(!strcasecmp(cptr, "TSV")) {csv_separator=2; output_format=0; continue;}
      if(!strcasecmp(cptr, "CSV")) {csv_separator=0; output_format=0; continue;}
      if(!strcasecmp(cptr, "FLOAT")) {output_format=1; continue;}
      if(!strcasecmp(cptr, "SHORT")) {output_format=2; continue;}
    } else if(strncasecmp(cptr, "INVERT", 3)==0) {
      invert=1; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}  
  
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) {
    if(atoiCheck(argv[ai], &dim) || dim<3) {
      fprintf(stderr, "Error: invalid dimension.\n");
      return(1);  
    }
    ai++;
  }
  if(ai<argc) {
    if(atofCheck(argv[ai], &diameter) || diameter<1.0 || diameter==diameter2) {
      fprintf(stderr, "Error: invalid diameter.\n");
      return(1);  
    }
    ai++;
  }
  if(ai<argc) {strlcpy(outfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing or wrong? */
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(invert && iring>0) {
    fprintf(stderr, "Error: options invert and iring cannot be combined.\n");
    return(1);
  }
  if(invert && segments>1) {
    fprintf(stderr, "Error: options invert and segm cannot be combined.\n");
    return(1);
  }
 
  /* Make sure that optional diameter2 is larger than diameter */
  if(diameter2>0.0 && diameter>diameter2) {
    double v; v=diameter2; diameter2=diameter; diameter=v;
  }
  
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("outfile := %s\n", outfile);
    printf("dim := %d\n", dim);
    printf("diameter := %g\n", diameter);
    if(diameter2>0.0) printf("diameter2 := %g\n", diameter2);
    if(segments>1) printf("segments := %d\n", segments);
    if(iring!=0) printf("iring := %d\n", iring);
    printf("csv_separator := %d\n", csv_separator);
    printf("output_format := %d\n", output_format);
    printf("invert := %d\n", invert);
  }
 
 
  /*
   *  Allocate memory for int list
   */
  int *ilist, *lptr;
  ilist=(int*)malloc(dim*dim*sizeof(int));
  if(ilist==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    return(2);
  }
 
 
  /*
   *  Make the matrix data as a list of integers
   */
 
  int i=0;
  if(diameter2<=0.0) {
    int x, y;
    double cx, cy, dx, dy, d, radius;
    /*
     *  Write 1 inside circle, 0 outside
     */
    i=roundf(diameter);
    cx=cy=0.5*(double)dim;
    if((dim%2==0 && i%2==0) || (dim%2!=0 && i%2!=0)) {cx+=0.5; cy+=0.5;} 
    //printf("center=%g,%g\n", cx,cy);
    radius=0.5*diameter;
    lptr=ilist; i=0;
    for(y=1; y<=dim; y++) for(x=1; x<=dim; x++) {
      dx=cx-(double)x; dy=cy-(double)y;
      d=hypot(dx, dy); 
      if(d<=radius) {
        *lptr=1;
        if(segments>1 && (dx!=0.0 || dy!=0.0)) {
          double ra=M_PI+atan2(dx,dy);
          int s=1+ra*(double)segments/(2.0*M_PI);
          *lptr=s;
        }
        i++;
      } else
        *lptr=0;
      lptr++;
      if(verbose>10) printf("d(%d,%d) := %g\n", x, y, d);
    }
  } else { 
    int x, y;
    double cx, cy, dx, dy, d, radius, radius2;
    /*
     *  Write 1 inside a ring, 0 outside
     */
    i=roundf(diameter);
    cx=cy=0.5*(double)dim;
    if((dim%2==0 && i%2==0) || (dim%2!=0 && i%2!=0)) {cx+=0.5; cy+=0.5;} 
    //printf("center=%g,%g\n", cx,cy);
    radius=0.5*diameter; radius2=0.5*diameter2;
    lptr=ilist; i=0;
    for(y=1; y<=dim; y++) for(x=1; x<=dim; x++) {
      dx=cx-(double)x; dy=cy-(double)y;
      d=hypot(dx, dy);
      if(d>radius2) {
        *lptr=0;
      } else if(d<=radius) {
        *lptr=iring;
      } else {
        *lptr=1+iring;
        if(segments>1 && (dx!=0.0 || dy!=0.0)) {
          double ra=M_PI+atan2(dx,dy);
          int s=1+iring+ra*(double)segments/(2.0*M_PI);
          *lptr=s;
        }
        i++;
      }
      lptr++;
      if(verbose>10) printf("d(%d,%d) := %g\n", x, y, d);
    }
  }
  if(i==0) {
    fprintf(stderr, "Error: with these settings all matrix values are 0.\n");
    free(ilist); return(4);
  }
  /* Invert ones and zeroes if requested */
  if(invert!=0) {
    lptr=ilist;
    for(int i=0; i<dim*dim; i++, lptr++) if(*lptr==0) *lptr=1; else *lptr=0;
  }
 
 
  /*
   *  Write ints into CSV or TSV file, if requested
   */
  if(output_format==0) {
    int x, y, ret;
    CSV csv; csvInit(&csv);
    if(csvAllocate(&csv, dim*dim)!=TPCERROR_OK) {
      fprintf(stderr, "Error: out of memory.\n");
      free(ilist); return(6);
    }
    lptr=ilist; ret=TPCERROR_OK;
    for(y=0; y<dim && ret==TPCERROR_OK; y++) {
      for(x=0; x<dim && ret==TPCERROR_OK; x++, lptr++) {
        if(x==0 && y>0) ret=csvPutInt(&csv, *lptr, 1);
        else ret=csvPutInt(&csv, *lptr, 0);
      }
    }
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot fill CSV data.\n");
      free(ilist); csvFree(&csv); return(7);
    }  
    free(ilist);
    /* Set the field separator */
    if(csv_separator==0) csv.separator=',';
    else if(csv_separator==1) csv.separator=';';
    else if(csv_separator==2) csv.separator='\t';
    else csv.separator=' ';
    /* Open file */
    FILE *fp;
    if(strcasecmp(outfile, "stdout")==0) {
      fp=stdout;
    } else {
      fp=fopen(outfile, "w");
      if(fp==NULL) {
        fprintf(stderr, "Error: cannot open file for write.\n");
        csvFree(&csv); return(10);
      }
    }
    /* Write matrix data */
    ret=csvWrite(&csv, 0, fp, &status);
    csvFree(&csv); fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      if(verbose>1) printf("ret := %d\n", ret);
      return(11);
    }
    // We are ready
    return(0);
  }
 
  /*
   *  Write matrix data as flat binary file
   */
  /* Open file */
  FILE *fp;
  fp=fopen(outfile, "wb");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for write.\n");
    free(ilist); return(10);
  }
  /* Write data */
  lptr=ilist; 
  if(output_format==1) { // binary floats
    float f;
    for(int i=0; i<dim*dim; i++, lptr++) {
      f=(float)*lptr;
      if(fwrite(&f, sizeof(float), 1, fp) != 1) {
        fprintf(stderr, "Error: cannot write in %s\n", outfile);
        fclose(fp); free(ilist); return(12);
      }
    }
  } else if(output_format==2) { // binary short int
    short int s;
    for(i=0; i<dim*dim; i++, lptr++) {
      s=(short int)*lptr;
      if(fwrite(&s, sizeof(short int), 1, fp) != 1) {
        fprintf(stderr, "Error: cannot write in %s\n", outfile);
        fclose(fp); free(ilist); return(12);
      }
    }
  }
  fclose(fp); free(ilist);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/