Drop density estimator

First of all, it does not look too pretty. But more importantly now we
can reduce the bucket size to four (instead of five) and vectorize that
stuff.

9 files changed, 22 insertions, 545 deletions

diff --git a/filters.c b/filters.c
index 8bb9aa4..392c7eb 100644
--- a/filters.c
+++ b/filters.c
@@ -273,143 +273,6 @@ int flam3_create_spatial_filter(flam3_frame *spec, int field, double **filter) {
    return (fwidth);
 }
 
-flam3_de_helper flam3_create_de_filters(double max_rad, double min_rad, double curve, int ss) {
-
-   flam3_de_helper de;
-   double comp_max_radius, comp_min_radius;
-   double num_de_filters_d;
-   int num_de_filters,de_max_ind;
-   int de_row_size, de_half_size;
-   int filtloop;
-   int keep_thresh=100;
-
-   de.kernel_size=-1;
-
-   if (curve <= 0.0) {
-      fprintf(stderr,"estimator curve must be > 0\n");
-      return(de);
-   }
-
-   if (max_rad < min_rad) {
-      fprintf(stderr,"estimator must be larger than estimator_minimum.\n");
-      fprintf(stderr,"(%f > %f) ? \n",max_rad,min_rad);
-      return(de);
-   }
-
-   /* We should scale the filter width by the oversample          */
-   /* The '+1' comes from the assumed distance to the first pixel */
-   comp_max_radius = max_rad*ss + 1;
-   comp_min_radius = min_rad*ss + 1;
-
-   /* Calculate how many filter kernels we need based on the decay function */
-   /*                                                                       */
-   /*    num filters = (de_max_width / de_min_width)^(1/estimator_curve)    */
-   /*                                                                       */
-   num_de_filters_d = pow( comp_max_radius/comp_min_radius, 1.0/curve );
-   if (num_de_filters_d>1e7) {
-      fprintf(stderr,"too many filters required in this configuration (%g)\n",num_de_filters_d);
-      return(de);
-   }
-   num_de_filters = (int)ceil(num_de_filters_d);
-         
-   /* Condense the smaller kernels to save space */
-   if (num_de_filters>keep_thresh) { 
-      de_max_ind = (int)ceil(DE_THRESH + pow(num_de_filters-DE_THRESH,curve))+1;
-      de.max_filtered_counts = (int)pow( (double)(de_max_ind-DE_THRESH), 1.0/curve) + DE_THRESH;
-   } else {
-      de_max_ind = num_de_filters;
-      de.max_filtered_counts = de_max_ind;
-   }
-
-   /* Allocate the memory for these filters */
-   /* and the hit/width lookup vector       */
-   de_row_size = (int)(2*ceil(comp_max_radius)-1);
-   de_half_size = (de_row_size-1)/2;
-   de.kernel_size = (de_half_size+1)*(2+de_half_size)/2;
-
-   de.filter_coefs = (double *) calloc (de_max_ind * de.kernel_size,sizeof(double));
-   de.filter_widths = (double *) calloc (de_max_ind,sizeof(double));
-
-   /* Generate the filter coefficients */
-   de.max_filter_index = 0;
-   for (filtloop=0;filtloop<de_max_ind;filtloop++) {
-
-      double de_filt_sum=0.0, de_filt_d;
-      double de_filt_h;
-      int dej,dek;
-      double adjloop;
-      int filter_coef_idx;
-
-      /* Calculate the filter width for this number of hits in a bin */
-      if (filtloop<keep_thresh)
-         de_filt_h = (comp_max_radius / pow(filtloop+1,curve));
-      else {
-         adjloop = pow(filtloop-keep_thresh,(1.0/curve)) + keep_thresh;
-         de_filt_h = (comp_max_radius / pow(adjloop+1,curve));
-      }
-
-      /* Once we've reached the min radius, don't populate any more */
-      if (de_filt_h <= comp_min_radius) {
-         de_filt_h = comp_min_radius;
-         de.max_filter_index = filtloop;
-      }
-
-      de.filter_widths[filtloop] = de_filt_h;
-
-      /* Calculate norm of kernel separately (easier) */
-      for (dej=-de_half_size; dej<=de_half_size; dej++) {
-         for (dek=-de_half_size; dek<=de_half_size; dek++) {
-            
-            de_filt_d = sqrt( (double)(dej*dej+dek*dek) ) / de_filt_h;
-
-            /* Only populate the coefs within this radius */
-            if (de_filt_d <= 1.0) {
-
-               /* Gaussian */
-               de_filt_sum += flam3_spatial_filter(flam3_gaussian_kernel,
-                        flam3_spatial_support[flam3_gaussian_kernel]*de_filt_d);
-
-               /* Epanichnikov */
-//             de_filt_sum += (1.0 - (de_filt_d * de_filt_d));
-            }
-         }
-      }
-
-      filter_coef_idx = filtloop*de.kernel_size;
-
-      /* Calculate the unique entries of the kernel */
-      for (dej=0; dej<=de_half_size; dej++) {
-         for (dek=0; dek<=dej; dek++) {
-            de_filt_d = sqrt( (double)(dej*dej+dek*dek) ) / de_filt_h;
-
-            /* Only populate the coefs within this radius */
-            if (de_filt_d>1.0)
-               de.filter_coefs[filter_coef_idx] = 0.0;
-            else {
-
-               /* Gaussian */
-               de.filter_coefs[filter_coef_idx] = flam3_spatial_filter(flam3_gaussian_kernel,
-                        flam3_spatial_support[flam3_gaussian_kernel]*de_filt_d)/de_filt_sum; 
-                            
-               /* Epanichnikov */
-//             de_filter_coefs[filter_coef_idx] = (1.0 - (de_filt_d * de_filt_d))/de_filt_sum;
-            }
-                  
-            filter_coef_idx ++;
-         }
-      }
-
-      if (de.max_filter_index>0)
-         break;
-   }
-
-   if (de.max_filter_index==0)
-      de.max_filter_index = de_max_ind-1;
-
-   
-   return(de);
-}                       
-
 double flam3_create_temporal_filter(int numsteps, int filter_type, double filter_exp, double filter_width,
                                     double **temporal_filter, double **temporal_deltas) {
 
diff --git a/filters.h b/filters.h
index c7c18a2..1f7ac16 100644
--- a/filters.h
+++ b/filters.h
@@ -21,21 +21,10 @@
 
 #include "private.h"
 
-#define DE_THRESH   100
-
-typedef struct {
-   int max_filtered_counts;
-   int max_filter_index;
-   int kernel_size;
-   double *filter_widths;
-   double *filter_coefs;
-} flam3_de_helper;
-
 extern double flam3_spatial_support[flam3_num_spatialfilters];
 
 double flam3_spatial_filter(int knum, double x);
 int flam3_create_spatial_filter(flam3_frame *spec, int field, double **filter);
-flam3_de_helper flam3_create_de_filters(double max_rad, double min_rad, double curve, int ss);
 double flam3_create_temporal_filter(int numsteps, int filter_type, double filter_exp, double filter_width,
                                     double **temporal_filter, double **temporal_deltas);
  
diff --git a/flam3-genome.c b/flam3-genome.c
index 2ba04d1..a6a8611 100644
--- a/flam3-genome.c
+++ b/flam3-genome.c
@@ -54,9 +54,6 @@ void test_cp(flam3_genome *cp) {
    cp->sample_density = 1;
    cp->nbatches = 1;
    cp->ntemporal_samples = 1;
-   cp->estimator = 0.0;
-   cp->estimator_minimum = 0.0;
-   cp->estimator_curve = 0.6;
 }
 
 flam3_genome *string_to_cp(char *s, int *n, randctx * const rc) {
diff --git a/flam3.c b/flam3.c
index 11e04e0..2380c6f 100644
--- a/flam3.c
+++ b/flam3.c
@@ -1141,10 +1141,6 @@ void clear_cp(flam3_genome *cp, int default_flag) {
        cp->spatial_filter_radius = 0.5;
        cp->zoom = 0.0;
        cp->sample_density = 1;
-       /* Density estimation stuff defaulting to ON */
-       cp->estimator = 9.0;
-       cp->estimator_minimum = 0.0;
-       cp->estimator_curve = 0.4;
        cp->gam_lin_thresh = 0.01;
 //       cp->motion_exp = 0.0;
        cp->nbatches = 1;
@@ -1170,9 +1166,6 @@ void clear_cp(flam3_genome *cp, int default_flag) {
        cp->width = -1;
        cp->height = -1;
        cp->sample_density = -1;
-       cp->estimator = -1;
-       cp->estimator_minimum = -1;
-       cp->estimator_curve = -1;
        cp->gam_lin_thresh = -1;
 //       cp->motion_exp = -999;
        cp->nbatches = 0;
@@ -1404,12 +1397,6 @@ void flam3_apply_template(flam3_genome *cp, flam3_genome *templ) {
    }
    if (templ->height > 0)
       cp->height = templ->height;
-   if (templ->estimator >= 0)
-      cp->estimator = templ->estimator;
-   if (templ->estimator_minimum >= 0)
-      cp->estimator_minimum = templ->estimator_minimum;
-   if (templ->estimator_curve >= 0)
-      cp->estimator_curve = templ->estimator_curve;
    if (templ->gam_lin_thresh >= 0)
       cp->gam_lin_thresh = templ->gam_lin_thresh;
    if (templ->nbatches>0)
@@ -1559,8 +1546,6 @@ void flam3_print(FILE *f, flam3_genome *cp, char *extra_attributes, int print_ed
       fprintf(f, " highlight_power=\"%g\"", cp->highlight_power);
       
    fprintf(f, " vibrancy=\"%g\"", cp->vibrancy);
-   fprintf(f, " estimator_radius=\"%g\" estimator_minimum=\"%g\" estimator_curve=\"%g\"",
-      cp->estimator, cp->estimator_minimum, cp->estimator_curve);
    fprintf(f, " gamma_threshold=\"%g\"", cp->gam_lin_thresh);
    
    if (flam3_palette_mode_step == cp->palette_mode)
@@ -3353,8 +3338,6 @@ typedef float abucket_float[4];
    if (UINT_MAX - dest > delta) dest += delta; else dest = UINT_MAX; \
 } while (0)
 #define iter_thread iter_thread_float
-#define de_thread_helper de_thread_helper_33
-#define de_thread de_thread_33
 #include "rect.c"
 #undef iter_thread
 #undef render_rectangle
@@ -3363,8 +3346,6 @@ typedef float abucket_float[4];
 #undef abucket
 #undef bump_no_overflow
 #undef abump_no_overflow
-#undef de_thread_helper
-#undef de_thread
 
 int flam3_render(flam3_frame *spec, void *out,
         int field, int nchan, int trans, stat_struct *stats) {
diff --git a/flam3.h b/flam3.h
index b841956..0107672 100644
--- a/flam3.h
+++ b/flam3.h
@@ -488,12 +488,6 @@ typedef struct {
    int nbatches;
    int ntemporal_samples;
 
-   /* Density estimation parameters for blurring low density hits */
-   double estimator;             /* Filter width for bin with one hit */
-   double estimator_curve;              /* Exponent on decay function ( MAX / a^(k-1) ) */
-   double estimator_minimum;         /* Minimum filter width used -
-                                    forces filter to be used of at least this width on all pts */
-
    /* XML Edit structure */
    xmlDocPtr edits;
 
diff --git a/interpolation.c b/interpolation.c
index 09c3a7a..528e312 100644
--- a/interpolation.c
+++ b/interpolation.c
@@ -438,9 +438,6 @@ void flam3_interpolate_n(flam3_genome *result, int ncp,
    INTERP(rotate);
    INTERI(nbatches);
    INTERI(ntemporal_samples);
-   INTERP(estimator);
-   INTERP(estimator_minimum);
-   INTERP(estimator_curve);
    INTERP(gam_lin_thresh);
    
    /* Interpolate the chaos array */
diff --git a/palettes.c b/palettes.c
index adcc242..816f3a2 100644
--- a/palettes.c
+++ b/palettes.c
@@ -376,7 +376,6 @@ static double try_colors(flam3_genome *g, int color_resolution) {
 
     g->sample_density = 1;
     g->spatial_oversample = 1;
-    g->estimator = 0.0;
     
     /* Scale the image so that the total number of pixels is ~10000 */
     pixtotal = g->width * g->height;    
@@ -446,7 +445,6 @@ static double try_colors(flam3_genome *g, int color_resolution) {
     g->pixels_per_unit = saved.pixels_per_unit;
     g->nbatches = saved.nbatches;
     g->ntemporal_samples = saved.ntemporal_samples;
-    g->estimator = saved.estimator;
 
     /* Free xform storage */
     clear_cp(&saved,flam3_defaults_on);
diff --git a/parser.c b/parser.c
index 1952dec..d27bbe1 100644
--- a/parser.c
+++ b/parser.c
@@ -471,12 +471,6 @@ int parse_flame_element(xmlNode *flame_node, flam3_genome *loc_current_cp,
          cp->vibrancy = flam3_atof(att_str);
       } else if (!xmlStrcmp(cur_att->name, (const xmlChar *)"hue")) {
          cp->hue_rotation = fmod(flam3_atof(att_str), 1.0);
-      } else if (!xmlStrcmp(cur_att->name, (const xmlChar *)"estimator_radius")) {
-         cp->estimator = flam3_atof(att_str);
-      } else if (!xmlStrcmp(cur_att->name, (const xmlChar *)"estimator_minimum")) {
-         cp->estimator_minimum = flam3_atof(att_str);
-      } else if (!xmlStrcmp(cur_att->name, (const xmlChar *)"estimator_curve")) {
-         cp->estimator_curve = flam3_atof(att_str);
       } else if (!xmlStrcmp(cur_att->name, (const xmlChar *)"gamma_threshold")) {
          cp->gam_lin_thresh = flam3_atof(att_str);
       } else if (!xmlStrcmp(cur_att->name, (const xmlChar *)"soloxform")) {
diff --git a/rect.c b/rect.c
index 5a4482b..7a64cbd 100644
--- a/rect.c
+++ b/rect.c
@@ -36,201 +36,6 @@
 #define FUSE_28 100
 #define WHITE_LEVEL 255
 
-
-typedef struct {
-
-   bucket *b;
-   abucket *accumulate;
-   int width, height, oversample;
-   flam3_de_helper *de;
-   double k1,k2;
-   double curve;
-   int last_thread;
-   int start_row, end_row;
-   flam3_frame *spec;
-   int *aborted;
-   int progress_size;
-   
-} de_thread_helper;
-
-static void de_thread(void *dth) {
-
-   de_thread_helper *dthp = (de_thread_helper *)dth;
-   int oversample = dthp->oversample;   
-   int ss = (int)floor(oversample / 2.0);
-   int scf = !(oversample & 1);
-   double scfact = pow(oversample/(oversample+1.0), 2.0);
-   int wid=dthp->width;
-   int hig=dthp->height;
-   int ps =dthp->progress_size;
-   int str = (oversample-1)+dthp->start_row;
-   int enr = (oversample-1)+dthp->end_row;
-   int i,j;
-   time_t progress_timer=0;
-   struct timespec pauset;
-   int progress_count = 0;
-   int pixel_num;
-   
-   pauset.tv_sec = 0;
-   pauset.tv_nsec = 100000000;
-         
-   /* Density estimation code */         
-   for (j = str; j < enr; j++) {
-      for (i = oversample-1; i < wid-(oversample-1); i++) {
-
-         int ii,jj;
-         double f_select=0.0;
-         int f_select_int,f_coef_idx;
-         int arr_filt_width;
-         bucket *b;
-         double c[4],ls;
-               
-         b = dthp->b + i + j*wid;
-         
-         /* Don't do anything if there's no hits here */
-         if (b[0][4] == 0 || b[0][3] == 0)
-            continue;
-
-         /* Count density in ssxss area   */
-         /* Scale if OS>1 for equal iters */
-         for (ii=-ss; ii<=ss; ii++) {
-            for (jj=-ss; jj<=ss; jj++) {
-               b = dthp->b + (i + ii) + (j + jj)*wid;
-               f_select += b[0][4]/255.0;
-            }
-         }
-               
-         if (scf)
-            f_select *= scfact;
-                  
-         if (f_select > dthp->de->max_filtered_counts)
-            f_select_int = dthp->de->max_filter_index;                  
-         else if (f_select<=DE_THRESH)
-            f_select_int = (int)ceil(f_select)-1;
-         else
-            f_select_int = (int)DE_THRESH +
-               (int)floor(pow(f_select-DE_THRESH,dthp->curve));
-
-         /* If the filter selected below the min specified clamp it to the min */
-         if (f_select_int > dthp->de->max_filter_index)
-            f_select_int = dthp->de->max_filter_index;
-
-         /* We only have to calculate the values for ~1/8 of the square */
-         f_coef_idx = f_select_int*dthp->de->kernel_size;
-
-         arr_filt_width = (int)ceil(dthp->de->filter_widths[f_select_int])-1;
-
-         b = dthp->b + i + j*wid;
-
-         for (jj=0; jj<=arr_filt_width; jj++) {
-            for (ii=0; ii<=jj; ii++) {
-            
-               /* Skip if coef is 0 */
-               if (dthp->de->filter_coefs[f_coef_idx]==0.0) {
-                  f_coef_idx++;
-                  continue;
-               }
-                     
-               c[0] = (double) b[0][0];
-               c[1] = (double) b[0][1];
-               c[2] = (double) b[0][2];
-               c[3] = (double) b[0][3];
-
-               ls = dthp->de->filter_coefs[f_coef_idx]*(dthp->k1 * log(1.0 + c[3] * dthp->k2))/c[3];
-
-               c[0] *= ls;
-               c[1] *= ls;
-               c[2] *= ls;
-               c[3] *= ls;
-
-               if (jj==0 && ii==0) {
-                  add_c_to_accum(dthp->accumulate,i,ii,j,jj,wid,hig,c);
-               }
-               else if (ii==0) {
-                  add_c_to_accum(dthp->accumulate,i,jj,j,0,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,-jj,j,0,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,0,j,jj,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,0,j,-jj,wid,hig,c);
-               } else if (jj==ii) {
-                  add_c_to_accum(dthp->accumulate,i,ii,j,jj,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,-ii,j,jj,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,ii,j,-jj,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,-ii,j,-jj,wid,hig,c);
-               } else {
-                  add_c_to_accum(dthp->accumulate,i,ii,j,jj,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,-ii,j,jj,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,ii,j,-jj,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,-ii,j,-jj,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,jj,j,ii,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,-jj,j,ii,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,jj,j,-ii,wid,hig,c);
-                  add_c_to_accum(dthp->accumulate,i,-jj,j,-ii,wid,hig,c);
-               }
-
-               f_coef_idx++;
-
-            }
-         }
-      }
-            
-      pixel_num = (j-str+1)*wid;
-
-      if (dthp->last_thread) {
-         /* Standard progress function */
-         if (dthp->spec->verbose && time(NULL) != progress_timer) {
-            progress_timer = time(NULL);
-            fprintf(stderr, "\rdensity estimation: %d/%d          ", j-str, enr-str);
-            fflush(stderr);
-         }
-         
-      }
-      /* Custom progress function */
-      if (dthp->spec->progress && pixel_num > progress_count + ps) {
-      
-         progress_count = ps * floor(pixel_num/(double)ps);
-      
-         if (dthp->last_thread) {
-         
-            int rv = (*dthp->spec->progress)(dthp->spec->progress_parameter,
-				      100*(j-str)/(double)(enr-str), 1, 0);
-				      
-            if (rv==2) { /* PAUSE */
-            
-               *(dthp->aborted) = -1;
-                              
-               do {
-				   nanosleep(&pauset,NULL);
-                  rv = (*dthp->spec->progress)(dthp->spec->progress_parameter,
-                     100*(j-str)/(double)(enr-str), 1, 0);
-               } while (rv==2);
-
-               *(dthp->aborted) = 0;
-
-            }
-
-	        if (rv==1) {
-			   *(dthp->aborted) = 1;
-               pthread_exit((void *)0);
-            }
-         } else {
-
-            if (*(dthp->aborted)<0) {
-            
-                do {
-				   nanosleep(&pauset,NULL);
-                } while (*(dthp->aborted)<0);                            
-            }
-            
-            if (*(dthp->aborted)>0) pthread_exit((void *)0);
-         }
-      }
-
-   }
-
-     pthread_exit((void *)0);
-
-}
-
 static void iter_thread(void *fth) {
    double sub_batch;
    int j;
@@ -487,7 +292,6 @@ static int render_rectangle(flam3_frame *spec, void *out,
    int vib_gam_n = 0;
    time_t progress_began=0;
    int verbose = spec->verbose;
-   int gnm_idx,max_gnm_de_fw,de_offset;
    flam3_genome cp;
    unsigned short *xform_distrib;
    flam3_iter_constants fic;
@@ -588,32 +392,6 @@ static int render_rectangle(flam3_frame *spec, void *out,
    */
    gutter_width = (filter_width - oversample) / 2;
 
-   /* 
-      Check the size of the density estimation filter.
-      If the 'radius' of the density estimation filter is greater than the          
-      gutter width, we have to pad with more.  Otherwise, we can use the same value.
-   */
-   max_gnm_de_fw=0;
-   for (gnm_idx = 0; gnm_idx < spec->ngenomes; gnm_idx++) {
-      int this_width = (int)ceil(spec->genomes[gnm_idx].estimator) * oversample;
-      if (this_width > max_gnm_de_fw)
-         max_gnm_de_fw = this_width;
-   }
-
-   /* Add OS-1 for the averaging at the edges, if it's > 0 already */
-   if (max_gnm_de_fw>0)
-      max_gnm_de_fw += (oversample-1);
-
-   /* max_gnm_de_fw is now the number of pixels of additional gutter      */
-   /* necessary to appropriately perform the density estimation filtering */
-   /* Check to see if it's greater than the gutter_width                  */
-   if (max_gnm_de_fw > gutter_width) {
-      de_offset = max_gnm_de_fw - gutter_width;
-      gutter_width = max_gnm_de_fw;
-   } else
-      de_offset = 0;
-
-
    /* Allocate the space required to render the image */
    fic.height = oversample * image_height + 2 * gutter_width;
    fic.width  = oversample * image_width  + 2 * gutter_width;
@@ -644,45 +422,11 @@ static int render_rectangle(flam3_frame *spec, void *out,
 
    /* Batch loop - outermost */
    for (batch_num = 0; batch_num < nbatches; batch_num++) {
-      double de_time;
       double sample_density=0.0;
       double k1, area, k2;
-      flam3_de_helper de;
-
-      de_time = spec->time + temporal_deltas[batch_num*ntemporal_samples];
 
       memset((char *) buckets, 0, sizeof(bucket) * nbuckets);
 
-      /* interpolate and get a control point                      */
-      /* ONLY FOR DENSITY FILTER WIDTH PURPOSES                   */
-      /* additional interpolation will be done in the temporal_sample loop */
-      flam3_interpolate(spec->genomes, spec->ngenomes, de_time, 0, &cp);
-
-      /* if instructed to by the genome, create the density estimation */
-      /* filter kernels.  Check boundary conditions as well.           */
-      if (cp.estimator < 0.0 || cp.estimator_minimum < 0.0) {
-         fprintf(stderr,"density estimator filter widths must be >= 0\n");
-         return(1);
-      }
-
-      if (spec->bits <= 32) {
-         if (cp.estimator > 0.0) {
-            fprintf(stderr, "warning: density estimation disabled with %d bit buffers.\n", spec->bits);
-            cp.estimator = 0.0;
-         }
-      }
-
-      /* Create DE filters */
-      if (cp.estimator > 0.0) {
-         de = flam3_create_de_filters(cp.estimator,cp.estimator_minimum,
-                                      cp.estimator_curve,oversample);
-         if (de.kernel_size<0) {
-            fprintf(stderr,"de.kernel_size returned 0 - aborting.\n");
-            return(1);
-         }
-      } else
-         de.max_filter_index = 0;
-      
       /* Temporal sample loop */
       for (temporal_sample_num = 0; temporal_sample_num < ntemporal_samples; temporal_sample_num++) {
 
@@ -864,111 +608,31 @@ static int render_rectangle(flam3_frame *spec, void *out,
       printf("k1=%f,k2=%15.12f\n",k1,k2);
 #endif
 
-      if (de.max_filter_index == 0) {
-
-         for (j = 0; j < fic.height; j++) {
-            for (i = 0; i < fic.width; i++) {
-
-               abucket *a = accumulate + i + j * fic.width;
-               bucket *b = buckets + i + j * fic.width;
-               double c[4], ls;
-
-               c[0] = (double) b[0][0];
-               c[1] = (double) b[0][1];
-               c[2] = (double) b[0][2];
-               c[3] = (double) b[0][3];
-               if (0.0 == c[3])
-                  continue;
-
-               ls = (k1 * log(1.0 + c[3] * k2))/c[3];
-               c[0] *= ls;
-               c[1] *= ls;
-               c[2] *= ls;
-               c[3] *= ls;
-
-               abump_no_overflow(a[0][0], c[0]);
-               abump_no_overflow(a[0][1], c[1]);
-               abump_no_overflow(a[0][2], c[2]);
-               abump_no_overflow(a[0][3], c[3]);
-            }
-         }
-      } else {
-      
-         de_thread_helper *deth;
-         int de_aborted=0;
-         int myspan = (fic.height-2*(oversample-1)+1);
-         int swath = myspan/(spec->nthreads);
-                  
-         /* Create the de helper structures */
-         deth = (de_thread_helper *)calloc(spec->nthreads,sizeof(de_thread_helper));
-         
-         for (thi=0;thi<(spec->nthreads);thi++) {
-         
-            /* Set up the contents of the helper structure */
-            deth[thi].b = buckets;
-            deth[thi].accumulate = accumulate;
-            deth[thi].width = fic.width;
-            deth[thi].height = fic.height;
-            deth[thi].oversample = oversample;
-            deth[thi].progress_size = spec->sub_batch_size/10;
-            deth[thi].de = &de;
-            deth[thi].k1 = k1;
-            deth[thi].k2 = k2;
-            deth[thi].curve = cp.estimator_curve;
-            deth[thi].spec = spec;
-            deth[thi].aborted = &de_aborted;
-            if ( (spec->nthreads)>myspan) { /* More threads than rows */
-               deth[thi].start_row=0;
-               if (thi==spec->nthreads-1) {
-                  deth[thi].end_row=myspan;
-                  deth[thi].last_thread=1;
-               } else {
-                  deth[thi].end_row=-1;
-                  deth[thi].last_thread=0;
-               }
-            } else { /* Normal case */
-               deth[thi].start_row=thi*swath;
-               deth[thi].end_row=(thi+1)*swath;
-               if (thi==spec->nthreads-1) {
-                  deth[thi].end_row=myspan;
-                  deth[thi].last_thread=1;
-               } else {
-                  deth[thi].last_thread=0;
-               }
-            }
-         }
-
-         /* Let's make some threads */
-         myThreads = (pthread_t *)malloc(spec->nthreads * sizeof(pthread_t));
-
-         pthread_attr_init(&pt_attr);
-         pthread_attr_setdetachstate(&pt_attr,PTHREAD_CREATE_JOINABLE);
+      for (j = 0; j < fic.height; j++) {
+         for (i = 0; i < fic.width; i++) {
 
-         for (thi=0; thi <spec->nthreads; thi ++)
-            pthread_create(&myThreads[thi], &pt_attr, (void *)de_thread, (void *)(&(deth[thi])));
+            abucket *a = accumulate + i + j * fic.width;
+            bucket *b = buckets + i + j * fic.width;
+            double c[4], ls;
 
-         pthread_attr_destroy(&pt_attr);
+            c[0] = (double) b[0][0];
+            c[1] = (double) b[0][1];
+            c[2] = (double) b[0][2];
+            c[3] = (double) b[0][3];
+            if (0.0 == c[3])
+               continue;
 
-         /* Wait for them to return */
-         for (thi=0; thi < spec->nthreads; thi++)
-            pthread_join(myThreads[thi], NULL);
-         
-         free(myThreads);            
+            ls = (k1 * log(1.0 + c[3] * k2))/c[3];
+            c[0] *= ls;
+            c[1] *= ls;
+            c[2] *= ls;
+            c[3] *= ls;
 
-         free(deth);
-                  
-         if (de_aborted) {
-            if (verbose) fprintf(stderr, "\naborted!\n");
-            goto done;
+            abump_no_overflow(a[0][0], c[0]);
+            abump_no_overflow(a[0][1], c[1]);
+            abump_no_overflow(a[0][2], c[2]);
+            abump_no_overflow(a[0][3], c[3]);
          }
-
-      } /* End density estimation loop */
-
-
-      /* If allocated, free the de filter memory for the next batch */
-      if (de.max_filter_index > 0) {
-         free(de.filter_coefs);
-         free(de.filter_widths);
       }
 
    }
@@ -1050,9 +714,9 @@ static int render_rectangle(flam3_frame *spec, void *out,
       }
 
       /* Apply the spatial filter */
-      y = de_offset;
+      y = 0;
       for (j = 0; j < image_height; j++) {
-         x = de_offset;
+         x = 0;
          for (i = 0; i < image_width; i++) {
             int ii, jj,rgbi;
             void *p;