switch (ffmt)
      {
        case dsm_image_format_JPG:
          { int kind;
            pimg = uint16_image_io_jpeg_fread(rd, &kind);
            /* At present, {uint16_image_io_jpeg_read} can suport only two kinds. */
            /* Revise the next line if {uint16_image_io_jpeg_read} is expanded to support other kinds. */
            assert((kind == JCS_GRAYSCALE) || (kind == JCS_RGB));
            /* Use the client-given gamma for decoding: */
            gamma_file = NAN;
            bias_file = NAN;
          }
          break;

        case dsm_image_format_PNG:
          { bool_t verbose_png = TRUE;
            pimg = uint16_image_io_png_fread (rd, &gamma_file, verbose_png);
            gamma_file = 1/gamma_file; /* PNG "gAMA" gives the encoding gamma, not decoding. */
            bias_file = NAN;
          } 
          break;

        case dsm_image_format_PNM:
          { pimg = uint16_image_read_pnm_file(rd);
            /* Parameters that approximate the standard PNM ITU-R BT.709 decoding: */
            gamma_file = 1.0/0.450; 
            bias_file = 0.0327;
          }
          break;

        default:
          assert(FALSE);
      }



vector<struct sub_slope> visual_rythm (
     dsm_image_t *img, 
     struct dsm_vs_data *vs, 
     struct dsm_common_data *cd, 
     vector<KLT_FeatureList> &features, 
     vector<plate> plates_gt
   ) {

   vector<struct sub_slope> slopes;

   int imw = vs->img_width;
   int img->ny = vs->img_height;
   int vrw = cd->vrw;
   int vrh = cd->vrh;

   /*Extract visual rythm row and store into {cd->visual_rythm}.*/
   extract_rythm_row (img, imw, img->ny, vs->start, vs->xreduc, vs->xnlines, vs->ynlines, vs->vect_Y);
   
   for (int x = 0; x < vrw; x++) {
      cd->visual_rythm[vrh * vrw + x] = quantize(vs->vect_Y[x]);
   } 

   estimate_background_gain_and_weight_of_row (vs, cd);
 
   /*Image binarization: */
   for (int x = 0; x < vrw; x++) {
      double F = (double)(cd->visual_background[vrh * vrw + x]);
      double Y = (double)(cd->visual_rythm[vrh * vrw + x]);
      double coeff = cd->visual_gain[vrh];
      if ( fabs(Y - coeff * F) < 3 * vs->sigma_bg) {
         /*Background*/
         cd->visual_segmentation[vrh * vrw + x] = 0;
      }
      else { 
         /*Foreground*/
         cd->visual_segmentation[vrh * vrw + x] = 255;
      }
   }

   if ((vrh % 300) == 0) {
      write_uchar_to_pgm (cd->visual_rythm, vrw, vrh, "%s/snapshot/visual_rythm.pgm", cd->output_path);
      //write_uchar_to_pgm (cd->visual_update, cd->vrw, cd->vrh, "%s/snapshot/visual_update.pgm", cd->output_path);
      write_uchar_to_pgm (cd->visual_weights, vrw, vrh, "%s/snapshot/visual_weights.pgm", cd->output_path);
      write_uchar_to_pgm (cd->visual_background, vrw, vrh, "%s/snapshot/visual_background.pgm", cd->output_path);
      write_uchar_to_pgm (cd->visual_segmentation, vrw, vrh, "%s/snapshot/visual_segmentation.pgm", cd->output_path);

      char filename[512];
      snprintf(filename, sizeof(filename), "%s/snapshot/coeff_%05ld.txt", cd->output_path, cd->iframe);
      FILE *CF = fopen (filename, "w");
      for (int y = 0; y < vrh; y++) {
         fprintf(CF, "%f\n", cd->visual_gain[y]);
      }
      fclose(CF);
   }

   cd->vrh++;

   return slopes;
}

void adjust_weights_old (
     int vrw, 
     double Y[],     /* Frame row values [0..255], with {vrw} elements.*/
     double F[],     /* Background row values [0..255], with {vrw} elements.*/
     int maxiter, 
     int rclean,     /* Radius for weight clean-up. */ 
     double *sigma,  /* (IN/OUT) deviation of background errors. */
     double *coeff,  /* (IN/OUT) camera gain of this frame. */
     double W[]      /* (OUT) Fitting weights [0..1]. */
   ) {

   int iter = 0;

   do {

      /*Recompute the weights based on a Gaussian function:*/
      for (int x = 0; x < vrw; x++) {
         double d = (Y[x] - (*coeff) * F[x])/(*sigma);
         W[x] = exp(-d*d/2);
      }
     
      /*Clean the small background bits.*/
      if (rclean > 0) {
         erode_weigths  (W, vrw, rclean);
         dilate_weigths (W, vrw, rclean);
         dilate_weigths (W, vrw, rclean);
         erode_weigths  (W, vrw, rclean);
      } 

      /*Normalizing weights*/
      double max_w = 0.0;
      for (int x = 0; x < vrw; x++) { 
        if (W[x] > max_w) {
           max_w = W[x];
        }  
      }
      for (int x = 0; x < vrw; x++) { W[x] = (W[x]/max_w); }
       
      /*Estimate the camera gain by linear regression between background and rythm. */
 
      /* gsl_fit_wmul (const double * x, const size_t xstride, const double * w,
      const size_t wstride, const double * y, const size_t ystride, size_t n, 
      double *c1, double *cov11, double *sumsq). This function computes the best-fit
      linear regression coefficient c1 of the model Y = c_1 X for the weighted
      datasets (x, y), two vectors of length n with strides xstride and ystride. The
      vector w, of length n and stride wstride, specifies the weight of each datapoint. 
      The weight is the reciprocal of the variance for each datapoint in y. The variance 
      of the parameter c1 is computed using the weights and returned via the parameter 
      cov11. The weighted sum of squares of the residuals from the best-fit line, 
      \chi^2, is returned in chisq. Y: dependent variable. X: independent variable.*/
      double cov, sumsq;

      int fstride = 1, ystride = 1, wstride = 1;

      gsl_fit_wmul (F, fstride, W, wstride, Y, ystride, vrw, coeff, &cov, &sumsq);

      iter++;

   } while (iter < maxiter);

}

void adjust_weights_old2 (
     int vrw, 
     double Y[],     /* Frame row values [0..255], with {vrw} elements.*/
     double F[],     /* Background row values [0..255], with {vrw} elements.*/
     int maxiter, 
     int rclean,     /* Radius for weight clean-up. */ 
     double *sigma,  /* (IN/OUT) deviation of background errors. */
     double *coeff,  /* (IN/OUT) camera gain of this frame. */
     double W[]      /* (OUT) Fitting weights [0..1]. */
   ) {

   int iter = 0;

   do {

      /*Normalizing weights*/
      double max_w = 0.0;
      for (int x = 0; x < vrw; x++) { 
        if (W[x] > max_w) {
           max_w = W[x];
        }  
      }
      for (int x = 0; x < vrw; x++) { W[x] = (W[x]/max_w); }
       
      /*Estimate the camera gain by linear regression between background and rythm. */
 
      /* gsl_fit_wmul (const double * x, const size_t xstride, const double * w,
      const size_t wstride, const double * y, const size_t ystride, size_t n, 
      double *c1, double *cov11, double *sumsq). This function computes the best-fit
      linear regression coefficient c1 of the model Y = c_1 X for the weighted
      datasets (x, y), two vectors of length n with strides xstride and ystride. The
      vector w, of length n and stride wstride, specifies the weight of each datapoint. 
      The weight is the reciprocal of the variance for each datapoint in y. The variance 
      of the parameter c1 is computed using the weights and returned via the parameter 
      cov11. The weighted sum of squares of the residuals from the best-fit line, 
      \chi^2, is returned in chisq. Y: dependent variable. X: independent variable.*/
      double cov, sumsq;

      int fstride = 1, ystride = 1, wstride = 1;

      gsl_fit_wmul (F, fstride, W, wstride, Y, ystride, vrw, coeff, &cov, &sumsq);

      iter++;

   } while (iter < maxiter);

}

/*Estimates row {cd->vrh} of the background {cd->visual_background},
  the camera gain {cd->visual_gain}, and the background indicator
  weight {cd->visual_weights} from row {cd->vrh} of the visual rythm,
  and the previous {vs->ynlines_bg} rows of these images, with EMA decay constant {vs->tau_bg}.
  If {cd->vrh == 0}, uses {vs->bg_row} as the initial guess for the background.*/
void estimate_background_gain_and_weight_of_row (
     struct dsm_vs_data *vs, 
     struct dsm_common_data *cd)
{
   int my = vs->ynlines_bg;
   double sigma = vs->sigma_bg;
   double tau = vs->tau_bg;
   double beta = vs->beta_bg;
   double rclean = vs->cleanup_bg;
   bool use_new = false;
   int maxiter = 30;

   int vr_nx = cd->vr_nx; /*Width of visual rythm.*/
   int vr_ny = cd->vr_ny; /*Height of visual rythm.*/
   
   /*Get rythm, and guess initial background and gain of row {vr_ny}: */
   double coeff = (vr_ny == 0 ? 1.0 : cd->visual_gain[vr_ny-1]);
   for (int x = 0; x < vr_nx; x++) {
      vs->vect_Y[x] = (double)cd->visual_rythm[vr_ny * vr_nx + x]; 
      vs->vect_F[x] = (vr_ny == 0 ? vs->bg_row[x] : (double)cd->visual_background[(vr_ny - 1) * vr_nx + x]);
   }

   if (use_new) {
     adjust_weights_new (vr_nx, vs->vect_Y, vs->vect_F, maxiter, rclean, &sigma, &coeff, vs->vect_W);
   }
   else {
     adjust_weights_old (vr_nx, vs->vect_Y, vs->vect_F, maxiter, rclean, &sigma, &coeff, vs->vect_W);
   }
   cd->visual_gain[vr_ny] = coeff; 

   if (true) { //if (sigma < 1) {
      char filename[512];
      snprintf(filename, sizeof(filename), "%s/snapshot/fit_%05ld.txt", cd->output_path, cd->iframe);
      FILE *file = fopen(filename,"w");
      //fprintf(file, "#Sigma = %f\n", sigma);
      //fprintf(file, "#Coeff = %f\n", coeff);
      for (int x = 0; x < vr_nx; x++) {
         fprintf(file, "%6.2f  %6.2f %8.6f %f %f\n", vs->vect_F[x], vs->vect_Y[x], vs->vect_W[x], coeff, sigma);
      }
      fclose(file);
   }

   /*Convert weights to pixels:*/
   for (int x = 0; x < vr_nx; x++) {
      cd->visual_weights[vr_ny * vr_nx + x] = quantize(vs->vect_W[x] * 255);
   }

   /*Recompute the background based on the rythm and weights: */

   /*Update background images if the line is totally background: */
   double w = image_minb (cd->visual_weights, vr_nx, vr_ny - vs->band_ny, vr_ny) / 255.0;
   printf("Novo peso: %f - %f\n", w, image_minb (cd->visual_weights, vr_nx, vr_ny - vs->band_ny, vr_ny));

   for (int x = 0; x < vr_nx; x++) {
      double sum_wef = 0.0; 
      double sum_we = 0.0; 
      double sum_e = 0.0; 
      for (int y = 0; y <= my; y++) {
         int ky = vr_ny - y;
         if (ky >= 0) {
            double w = (double)cd->visual_weights[ky * vr_nx + x]/255.0;
            double e = exp(-y/tau); /*exponential moving average decay weight.*/
            double f = cd->visual_rythm[ky * vr_nx + x] / cd->visual_gain[ky];
            sum_wef += w * e * f; 
            sum_we += w * e; 
            sum_e += e; 
         } 
      }  
      double fold = vs->vect_F[x];
      double fnew = (sum_we == 0 ? fold : sum_wef/sum_we);
      double wnew = sum_we/sum_e;
      assert ((wnew >= 0) && (wnew <= 1));
      double fest = fold;
      if (w > 0.05) {
      if (wnew >= beta) {
         fest = fnew; //wnew * fnew + (1 - wnew) * fold;
      } }
      int fint = (int)(fest + 0.5);
      if (fint < 0) { fint = 0; }
      else if (fint > 255) { fint = 255; }
      //cd->visual_background[vr_ny * vr_nx + x] = (unsigned char)fint; 
      cd->visual_background[vr_ny * vr_nx + x] = (vr_ny == 0 ? vs->bg_row[x] : (double)cd->visual_background[(vr_ny - 1) * vr_nx + x]); 
   }
}

   int sum_y = 0;
   for (int x = 0; x < vr_nx; x++) {
      if (cd->visual_update.p[vr_ny * vr_nx + x] > 20) {
         sum_y++;
      }
   }
   if (sum_y == 0) {
      if ((vs->end - vs->beg) > 4) {
          // int *prof = (int *)malloc(vr_nx * sizeof(int));
          for (int x = 0; x < vr_nx; x++) {
             int sum = 0;
             for (int y = vs->beg; y <= vs->end; y++) {
                 if (cd->visual_update.p[y * vr_nx + x] > 20) {
                    sum++;
                 }
             }
             prof[x] = sum;
             /*if (sum > 3) {
                cd->visual_update[(vr_ny - sum) * vr_nx + x] = 255;
             }*/
             /*else {
                cd->visual_update[vr_ny * vr_nx + x] = 255;
             }*/  
          }      
          //paint_chrom_rect (img, nx, ny, 0, y0-300, nx, y0, 128, 0); 
         // window_smooth (prof, vr_nx, 1, 9);
         // for (int x = 0; x < vr_nx; x++) {
           //printf("prof: %d\n", prof[x]);
           //cd->visual_update[(vr_ny - prof[x]) * vr_nx + x] = 255;
         // }
         // free(prof);
      }
      vs->beg = vs->end = vr_ny;
   }
   else {
      //cd->visual_update[vr_ny * vr_nx + 5] = 255;
      //cd->visual_update[vr_ny * vr_nx + 10] = sum_y;
      vs->end++;
   }

  /*Given a frame {img}, extracts row 
    {y0} by averaging {ny} rows around it and {nx} pixels
    around each pixel, and subsamples each row by taking
    one every {xreduc} pixels. Stores the result in the 
    vector {R[0..imw/xreduc-1]},
    as doubles in the range [0..255].*/
  void extract_rythm_row (
			  dsm_image_t *img,
			  int y0,
			  int xreduc,
			  int nx, 
			  int ny,
			  double *R) {

    double *wx = dsm_compute_hahn_weights (nx); 
    double *wy = dsm_compute_hahn_weights (ny); 

    assert((nx % 2) == 1);
    int rx = nx/2;

    assert((ny % 2) == 1);
    int ry = ny/2;

    int wnew = img->nx / xreduc;

    for (int xnew = 0; xnew < wnew; xnew++) {
      int xold = xnew * xreduc + xreduc/2;
      assert (xold < img->nx);
      double sum_wp = 0.0;
      double sum_w = 0.0;
      for (int ix = -rx; ix <= +rx; ix++) {
	for (int iy = -ry; iy <= +ry; iy++) {
	  int kx = xold + ix; 
	  int ky = y0 + iy; 
	  if ((kx >= 0) && (kx < img->nx) && (ky >= 0) && (ky < img->ny)) {
	    double p = (double)img->p[ky * img->nx + kx];
	    double w = wx[ix + rx] * wy[iy + ry];
	    sum_wp += w * p;
	    sum_w += w;
	  }
	}
      }
      R[xnew] = (sum_w == 0 ? 0 : sum_wp/sum_w);
    }
    free(wx);
    free(wy);
  }

  void adjust_weights_new 
    ( int vrw, 
      double Y[],     /* Frame row values [0..255], with {vrw} elements.*/
      double F[],     /* Background row values [0..255], with {vrw} elements.*/
      int maxiter, 
      int rclean,     /* Radius for weight clean-up. */ 
      double *sigma,  /* (IN/OUT) deviation of background errors. */
      double *coeff,  /* (IN/OUT) camera gain of this frame. */
      double W[]      /* (OUT) Fitting weights [0..1]. */
    ) 
  {
  
    bool verbose = false;
 
    double avg_bad = 127, dev_bad = 1000.0;

    double *Wa = (double *)malloc(vrw * sizeof(double));

    for (int i = 0; i < vrw; i++) {
      Wa[i] = 1.0;
    } 

    lsq_robust_fit (vrw, F, Y, Wa, maxiter, rclean, coeff, sigma, &avg_bad, &dev_bad, W, verbose);

    free(Wa);
  }


void window_smooth (int *prof, int size, int iter, int wsize) {

    int *temp = (int *)malloc(size * sizeof(int));
    for (int j = 0; j < size; j++) {
      temp[j] = 0;
    }
    while (iter > 0) {
      for (int j = wsize / 2; j < size - wsize / 2; j++) {
	int sum = 0;
	for (int i = - wsize / 2; i <= + wsize / 2; i++) {
	  sum += prof[j + i];
	}
	temp[j] = sum/wsize;
      }
      for (int j = 0; j < size; j++) {
	prof[j] = temp[j];
      }
      iter--;
    }
    free(temp);
  }

      snprintf(filename, sizeof(filename), "%s/snapshot/frame_background_band_%05ld.pgm", cd->output_path, cd->iframe);
      write_uchar_to_pgm (vs->background_band.p, nx_band, ny_band, filename);
    
      snprintf(filename, sizeof(filename), "%s/snapshot/frame_background_band_%05ld.pgm", cd->output_path, cd->iframe);
      write_uchar_to_pgm (vs->background_band.p, nx_band, ny_band, filename);