src/snark/contur.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/contur.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  A PROCEDURE TO CONTOUR A RECONSTRUCTION INTO A BLACK-WHITE PICTURE
  11  */
  12
  13 #include <cstdio>
  14
  15 #include "blkdta.h"
  16 #include "geom.h"
  17 #include "math.h" //added. hstau
  18 #include "consts.h"
  19 #include "uiod.h"
  20
  21 #include "contur.h"
  22
  23 void contur(REAL* a, INTEGER m, INTEGER n, REAL thresh, REAL w1, REAL w2, REAL w3)
  24 {
  25         INTEGER i, last = m * n;
  26         REAL t = thresh;
  27
  28         if (w3 >= Consts.zero)
  29         {
  30                 REAL min, max, amin, amax, ai, at;
  31
  32                 if ((GeoPar.aveden <= w1) || (GeoPar.aveden >= w2))
  33                 {
  34                         fprintf(output,
  35                                         "\n **** invalid parameters to contour - (aveden <= w1 or w2 <= aveden)");
  36                         fprintf(output, "\n **** contour failed to continue");
  37                         return;
  38                 };
  39
  40                 // find min and max
  41                 min = a[0];
  42                 max = a[0];
  43                 for (i = 1; i < last; i++)
  44                 {
  45                         ai = a[i];
  46                         if (ai < min)
  47                                 min = ai;
  48                         if (ai > max)
  49                                 max = ai;
  50                 };
  51
  52                 if ((max - min) < Consts.zero)
  53                 {
  54                         fprintf(output, "\n **** picture is a constant");
  55                         fprintf(output, "\n **** contour failed to continue");
  56                         return;
  57                 };
  58
  59                 // perform a bisection search
  60                 min = min - 0.00001;
  61                 amin = w2;
  62                 amax = w1;
  63                 while ((amin != GeoPar.aveden) && (amax != GeoPar.aveden) && (max - min > 0.00001))
  64                 {
  65                         t = (min + max) / 2.0;
  66                         at = 0;
  67                         for (i = 0; i < last; i++)
  68                         {
  69                                 if (a[i] <= t)
  70                                         at += w1;
  71                                 else
  72                                         at += w2;
  73                         };
  74                         at /= last;
  75                         if (at < GeoPar.aveden)
  76                         {
  77                                 max = t;
  78                                 amax = at;
  79                         }
  80                         else
  81                         {
  82                                 min = t;
  83                                 amin = at;
  84                         }
  85                 };
  86
  87                 if (amin - GeoPar.aveden < GeoPar.aveden - amax)
  88                         t = min;
  89                 else
  90                         t = max;
  91
  92                 fprintf(output, "\n         contour threshold = %10.4f", t);
  93         };
  94
  95         for (i = 0; i < last; i++)
  96         {
  97                 if (a[i] <= t)
  98                         a[i] = w1;
  99                 else
 100                         a[i] = w2;
 101         };
 102
 103         return;
 104 }