/* Copyright (C) 2001-2026 Artifex Software, Inc.

   All Rights Reserved.

   This software is provided AS-IS with no warranty, either express or

   implied.

   This software is distributed under license and may not be copied,

   modified or distributed except as expressly authorized under the terms

   of the license contained in the file LICENSE in this distribution.

   Refer to licensing information at http://www.artifex.com or contact

   Artifex Software, Inc.,  39 Mesa Street, Suite 108A, San Francisco,

   CA 94129, USA, for further information.

*/

/* Ordered Dither Screen Creation Tool. */

#include <stdlib.h>

#ifdef GS_LIB_BUILD

#define LIB_BUILD

#include "std.h"

#include "string_.h"

#include "gsmemory.h"

#include "math_.h"

#include "gp.h"

#   define ALLOC(mem, size) (gs_alloc_bytes((gs_memory_t *)mem, size, "gen_ordered"))

#   define FREE(mem, ptr)   gs_free_object((gs_memory_t *)mem, ptr, "gen_ordered")

#   define PRINTF(mem, str) outprintf((gs_memory_t *)mem, str)

#   define PRINTF2(mem, str, v1, v2) outprintf((gs_memory_t *)mem, str, v1, v2)

#   define PRINTF4(mem, str, v1, v2, v3, v4) outprintf((gs_memory_t *)mem, str, v1, v2, v3, v4)

#   define PRINTF7(mem, str, v1, v2, v3, v4, v5, v6, v7) outprintf((gs_memory_t *)mem, str, v1, v2, v3, v4, v5, v6, v7)

#   define EPRINTF(mem, str) errprintf((gs_memory_t *)mem, str)

#   define EPRINTF1(mem, str, v1) errprintf((gs_memory_t *)mem, str, v1)

#   define EPRINTF3(mem, str, v1, v2, v3) errprintf((gs_memory_t *)mem, str, v1, v2, v3)

#   define _FILE gp_file

#   define FOPEN(mem, v1, v2) gp_fopen(mem, v1, v2)

#   define FCLOSE(fid) gp_fclose(fid)

#   define FPRINTF(fid, str) gp_fprintf(fid, str)

#   define FPRINTF1(fid, str, v1) gp_fprintf(fid, str, v1)

#   define FPRINTF2(fid, str, v1, v2) gp_fprintf(fid, str, v1, v2)

#   define FPRINTF6(fid, str, v1, v2, v3, v4, v5, v6) gp_fprintf(fid, str, v1, v2, v3, v4, v5, v6)

#   define FWRITE(v1, v2, v3, fid) gp_fwrite(v1, v2, v3, fid)

#endif /* defined GS_LIB_BUILD */

#ifndef LIB_BUILD

#include <stdio.h>

#include <string.h>

#include <memory.h>

#include <sys/stat.h>

#include <math.h>

typedef unsigned char byte;

#define false 0

#define true 1

#ifndef __cplusplus

    typedef int bool;

#endif      /* __cpluplus */

/* Needed if standalone (main) */

#   define ALLOC(mem, size) (malloc(size))

#   define FREE(mem, ptr)   (free(ptr))

#   define PRINTF(mem, str) printf(str)

#   define PRINTF2(mem, str, v1, v2) printf(str, v1, v2)

#   define PRINTF4(mem, str, v1, v2, v3, v4) printf(str, v1, v2, v3, v4)

#   define PRINTF7(mem, str, v1, v2, v3, v4, v5, v6, v7) printf(str, v1, v2, v3, v4, v5, v6, v7)

#   define EPRINTF(mem, str) fprintf(stderr, str)

#   define EPRINTF1(mem, str, v1) fprintf(stderr, str, v1)

#   define EPRINTF3(mem, str, v1, v2, v3) fprintf(stderr, str, v1, v2, v3)

#   define _FILE FILE

#   define FOPEN(mem, v1, v2) fopen(v1, v2)

#   define FCLOSE(fid) fclose(fid)

#   define FPRINTF(fid, str) fprintf(fid, str)

#   define FPRINTF1(fid, str, v1) fprintf(fid, str, v1)

#   define FPRINTF2(fid, str, v1, v2) fprintf(fid, str, v1, v2)

#   define FPRINTF6(fid, str, v1, v2, v3, v4, v5, v6) fprintf(fid, str, v1, v2, v3, v4, v5, v6)

#   define FWRITE(v1, v2, v3, fid) fwrite(v1, v2, v3, fid)

#endif /* ndef LIB_BUILD */

#include "gen_ordered.h"

typedef struct htsc_point_s {

    double x;

    double y;

} htsc_point_t;

typedef struct htsc_threshpoint {

    int x;

    int y;

    int value;

    int index;

    double dist_to_center;

} htsc_threshpoint_t;

typedef struct htsc_vertices_s {

    htsc_point_t lower_left;

    htsc_point_t upper_left;

    htsc_point_t upper_right;

    htsc_point_t lower_right;

} htsc_vertices_t;

typedef struct htsc_dot_shape_search_s {

    double norm;

    int index_x;

    int index_y;

} htsc_dot_shape_search_t;

typedef struct htsc_matrix_s {

   htsc_vector_t row[2];

} htsc_matrix_t;

typedef struct htsc_dither_pos_s {

    htsc_point_t *point;

    int number_points;

    int *locations;

} htsc_dither_pos_t;

static void htsc_determine_cell_shape(int *x, int *y, int *v, int *u,

                               int *N, htsc_param_t params, void *mem);

static double htsc_spot_value(spottype_t spot_type, double x, double y);

static int htsc_getpoint(htsc_dig_grid_t *dig_grid, int x, int y);

static void htsc_setpoint(htsc_dig_grid_t *dig_grid, int x, int y, int value);

static int  htsc_create_dot_mask(htsc_dig_grid_t *dig_grid, int x, int y, int u, int v,

                       double screen_angle, htsc_vertices_t vertices);

static void htsc_find_bin_center(htsc_dig_grid_t *dot_grid, htsc_vector_t *bin_center);

static int htsc_sumsum(htsc_dig_grid_t dig_grid);

static void htsc_create_dot_profile(htsc_dig_grid_t *dig_grid, int N, int x, int y,

                            int u, int v, double horiz_dpi,

                            double vert_dpi, htsc_vertices_t vertices,

                            htsc_point_t *one_index, spottype_t spot_type,

                            htsc_matrix_t trans_matrix);

static int htsc_allocate_supercell(htsc_dig_grid_t *super_cell, int x, int y, int u,

                           int v, int target_size, bool use_holladay_grid,

                           htsc_dig_grid_t dot_grid, int N, int *S, int *H, int *L);

static void htsc_tile_supercell(htsc_dig_grid_t *super_cell, htsc_dig_grid_t *dot_grid,

                 int x, int y, int u, int v, int N);

int create_2d_gauss_filter(double *filter, int x_size, int y_size,

                            double stdvalx, double stdvaly, gs_memory_t *mem);

static int htsc_create_holladay_mask(htsc_dig_grid_t super_cell, int H, int L,

                               double gamma, htsc_dig_grid_t *final_mask);

static int htsc_create_dither_mask(htsc_dig_grid_t super_cell,

                             htsc_dig_grid_t *final_mask, int verbose,

                             int num_levels, int y, int x, double vert_dpi,

                             double horiz_dpi, int N, double gamma,

                             htsc_dig_grid_t dot_grid, htsc_point_t one_index,

                             gs_memory_t *mem);

static int htsc_create_nondithered_mask(htsc_dig_grid_t super_cell, int H, int L,

                          double gamma, htsc_dig_grid_t *final_mask);

static int htsc_gcd(int a, int b);

static int  htsc_lcm(int a, int b);

static int htsc_matrix_inverse(htsc_matrix_t matrix_in, htsc_matrix_t *matrix_out);

static void htsc_matrix_vector_mult(htsc_matrix_t matrix_in, htsc_vector_t vector_in,

                   htsc_vector_t *vector_out);

static int htsc_mask_to_tos(htsc_dig_grid_t *final_mask);

#if RAW_SCREEN_DUMP

static int htsc_dump_screen(htsc_dig_grid_t *dig_grid, char filename[], gs_memory_t *mem);

static int htsc_dump_float_image(double *image, int height, int width, double max_val,

                      char filename[], gs_memory_t *mem);

static int htsc_dump_byte_image(byte *image, int height, int width, double max_val,

                      char filename[], gs_memory_t *mem);

#endif

/* Initialize default values */

void htsc_set_default_params(htsc_param_t *params)

{

    params->scr_ang = 0;

    params->targ_scr_ang = 0;

    params->targ_lpi = 75;

    params->vert_dpi = 300;

    params->horiz_dpi = 300;

    params->targ_quant_spec = false;

    params->targ_quant = 256;

    params->targ_size = 1;

    params->targ_size_spec = false;

    params->spot_type = CIRCLE;

    params->holladay = false;

    params->gamma = 1.0;

    params->output_format = OUTPUT_TOS;

    params->verbose = 0;

}

int

htsc_gen_ordered(htsc_param_t params, int *S, htsc_dig_grid_t *final_mask, gs_memory_t *mem)

{

    int num_levels;

    int x=0, y=0, v=0, u=0, N=0;

    htsc_vertices_t vertices;

    htsc_vector_t bin_center;

    htsc_point_t one_index = { 0., 0. };

    htsc_dig_grid_t dot_grid;

    htsc_dig_grid_t super_cell;

    int code;

    int H, L;

    htsc_matrix_t trans_matrix, trans_matrix_inv;

    int min_size;

    dot_grid.data = NULL;

    dot_grid.memory = final_mask->memory;

    super_cell.data = NULL;

    super_cell.memory = final_mask->memory;

    final_mask->data = NULL;

    params.targ_scr_ang = params.targ_scr_ang % 90;

    params.scr_ang = params.targ_scr_ang;

    /* Get the vector values that define the small cell shape */

    htsc_determine_cell_shape(&x, &y, &v, &u, &N, params, final_mask->memory);

    /* Figure out how many levels to dither across. */

    if (params.targ_quant_spec) {

        num_levels = ROUND((double) params.targ_quant / (double)N);

    } else {

        num_levels = 1;

    }

    if (num_levels < 1) num_levels = 1;

    if (num_levels == 1) {

        if (params.verbose > 0)

            PRINTF(final_mask->memory, "No additional dithering , creating minimal sized periodic screen\n");

        params.targ_size = 1;

    }

    /* Lower left of the cell is at the origin.  Define the other vertices */

    vertices.lower_left.x = 0;

    vertices.lower_left.y = 0;

    vertices.upper_left.x = x;

    vertices.upper_left.y = y;

    vertices.upper_right.x = x + u;

    vertices.upper_right.y = y + v;

    vertices.lower_right.x = u;

    vertices.lower_right.y = v;

    /* Create the matrix that is used to get us correctly into the dot shape function */

    trans_matrix.row[0].xy[0] = u;

    trans_matrix.row[0].xy[1] = x;

    trans_matrix.row[1].xy[0] = v;

    trans_matrix.row[1].xy[1] = y;

    code = htsc_matrix_inverse(trans_matrix, &trans_matrix_inv);

    if (code < 0) {

        EPRINTF(final_mask->memory, "ERROR! Singular Matrix Inversion!\n");

        return -1;

    }

    /* Create a binary mask that indicates where we need to define the dot turn

       on sequence or dot profile */

    code = htsc_create_dot_mask(&dot_grid, x, y, u, v, params.scr_ang, vertices);

    if (code < 0) {

        return code;

    }

#if RAW_SCREEN_DUMP

    code = htsc_dump_screen(&dot_grid, "mask", mem);

    if (code < 0)

        return code;

#endif

    /* A sanity check */

    if (htsc_sumsum(dot_grid) != -N) {

        EPRINTF(final_mask->memory, "ERROR! grid size problem!\n");

        return -1;

    }

    /* From the binary mask, find the center point.  This is needed to remove

       ambiguity during the TOS calculation from the dot profile.  We want to

       turn on those dots that are closests to the center first when there

       are ties */

    htsc_find_bin_center(&dot_grid, &bin_center);

    /* Now actually determine the turn on sequence */

    htsc_create_dot_profile(&dot_grid, N, x, y, u, v, params.horiz_dpi,

                            params.vert_dpi, vertices, &one_index,

                            params.spot_type, trans_matrix_inv);

#if RAW_SCREEN_DUMP

    code = htsc_dump_screen(&dot_grid, "dot_profile", mem);

    if (code < 0)

        return code;

#endif

    /* Allocate super cell */

    code = htsc_allocate_supercell(&super_cell, x, y, u, v, params.targ_size,

                            params.holladay, dot_grid, N, S, &H, &L);

    if (code < 0) {

        EPRINTF(final_mask->memory, "ERROR! grid size problem!\n");

        return -1;

    }

    /* Make a warning about large requested quantization levels with no -s set */

    if (params.targ_size == 1 && num_levels > 1) {

        min_size = (int)ceil((double)params.targ_quant / (double)N);

        EPRINTF1(final_mask->memory, "To achieve %d quantization levels with the desired lpi,\n", params.targ_quant);

        EPRINTF1(final_mask->memory, "it is necessary to specify a SuperCellSize (-s) of at least %d.\n", min_size);

        EPRINTF(final_mask->memory, "Note that an even larger size may be needed to reduce pattern artifacts.\n");

        EPRINTF(final_mask->memory, "Because no SuperCellSize was specified, the minimum possible size\n");

        EPRINTF(final_mask->memory, "that can approximate the requested angle and lpi will be created\n");

        EPRINTF1(final_mask->memory, "with %d quantization levels.\n", N);

    }

    /* Go ahead and fill up the super cell grid with our growth dot values */

    htsc_tile_supercell(&super_cell, &dot_grid, x, y, u, v, N);

#if RAW_SCREEN_DUMP

    code = htsc_dump_screen(&super_cell, "super_cell_tiled", mem);

    if (code < 0)

        return code;

#endif

    /* If we are using the Holladay grid (non dithered) then we are done. */

    if (params.holladay) {

        code = htsc_create_holladay_mask(super_cell, H, L, params.gamma, final_mask);

    } else {

        if ((super_cell.height == dot_grid.height &&

            super_cell.width == dot_grid.width) || num_levels == 1) {

            code = htsc_create_nondithered_mask(super_cell, H, L, params.gamma, final_mask);

        } else {

            /* Dont allow nonsense settings */

            if (num_levels * N > super_cell.height * super_cell.width) {

                EPRINTF3(final_mask->memory,

                         "Notice, %d quantization levels not possible with super cell of %d by %d\n",

                         num_levels, super_cell.height, super_cell.width);

                num_levels = ROUND((super_cell.height * super_cell.width) / (double)N);

                EPRINTF1(final_mask->memory, "Reducing dithering quantization to %d\n",

                         num_levels);

                EPRINTF1(final_mask->memory, "For an effective quantization of %d\n",

                         super_cell.height * super_cell.width);

            }

            code = htsc_create_dither_mask(super_cell, final_mask, params.verbose, num_levels,

                                           y, x, params.vert_dpi, params.horiz_dpi, N,

                                           params.gamma, dot_grid, one_index, mem);

        }

    }

    final_mask->bin_center = bin_center;

    /* Now if the requested format is turn-on-sequence, convert the "data" */

    if (code == 0 && params.output_format == OUTPUT_TOS) {

        code = htsc_mask_to_tos(final_mask);

    }

    /* result in in final_mask, clean up working arrays allocated. */

    FREE(final_mask->memory, dot_grid.data);

    FREE(final_mask->memory, super_cell.data);

    return code;

}

/* comparison for use in qsort */

static int

compare(const void * a, const void * b)

{

    const htsc_threshpoint_t *val_a = a;

    const htsc_threshpoint_t *val_b = b;

    double cost = val_a->value - val_b->value;

    /* If same value, use distance to center for decision */

    if (cost == 0)

        cost = val_a->dist_to_center - val_b->dist_to_center;

    if (cost == 0)

        return 0;

    if (cost < 0)

        return -1;

    return 1;

}

static int

htsc_mask_to_tos(htsc_dig_grid_t *final_mask)

{

    int width = final_mask->width;

    int height = final_mask->height;

    htsc_vector_t center = final_mask->bin_center;

    int *buff_ptr = final_mask->data;

    int x, y, k = 0;

    int count;

    htsc_threshpoint_t *values;

    int *tos;

    if (check_int_multiply(height, width, &count) < 0)

        return -1;

    if (check_int_multiply(count, sizeof(htsc_threshpoint_t), &x) < 0)

        return -1;

    values = (htsc_threshpoint_t *) ALLOC(final_mask->memory,

                                          sizeof(htsc_threshpoint_t) * width * height);

    if (values == NULL) {

        EPRINTF(final_mask->memory, "ERROR! malloc failure in htsc_mask_to_tos!\n");

        return -1;

    }

    if (check_int_multiply(count, 2, &x) < 0)

        return -1;

    if (check_int_multiply(x, sizeof(int), &x) < 0)

        return -1;

    tos = (int *) ALLOC(final_mask->memory, sizeof(int) * 2 * height * width);

    if (tos == NULL) {

        FREE(final_mask->memory, values);

        EPRINTF(final_mask->memory, "ERROR! malloc failure in htsc_mask_to_tos!\n");

        return -1;

    }

    /* Do a sort on the values and then output the coordinates */

    /* First get a list made with the unsorted values and coordinates */

    for (y = 0; y < height; y++) {

        for ( x = 0; x < width; x++ ) {

            values[k].value = *buff_ptr;

            values[k].x = x;

            values[k].y = y;

            values[k].index = k;

            values[k].dist_to_center = (x - center.xy[0]) * (x - center.xy[0]) +

                                       (y - center.xy[1]) * (y - center.xy[1]);

            buff_ptr++;

            k = k + 1;

        }

    }

#if RAW_SCREEN_DUMP

    EPRINTF(final_mask->memory, "Unsorted\n");

    for (k = 0; k < count; k++) {

        EPRINTF(final_mask->memory, "Index %d : x = %d y = %d dist = %4.2lf value = %d \n",

               values[k].index, values[k].x, values[k].y, values[k].dist_to_center, values[k].value);

    }

#endif

        /* Sort */

    qsort(values, (size_t)height * width, sizeof(htsc_threshpoint_t), compare);

#if RAW_SCREEN_DUMP

    EPRINTF(final_mask->memory, "Sorted\n");

    for (k = 0; k < count; k++) {

        EPRINTF(final_mask->memory, "Index %d : x = %d y = %d dist = %4.2lf value = %d \n",

                values[k].index, values[k].x, values[k].y, values[k].dist_to_center, values[k].value);

    }

#endif

    FREE(final_mask->memory, final_mask->data);

    final_mask->data = tos;

    buff_ptr = tos;

    for (k=0; k < count; k++) {

        *buff_ptr++ = values[count - 1 - k].x;

        *buff_ptr++ = values[count - 1 - k].y;

    }

    FREE(final_mask->memory, values);

    return 0;

}

static void

htsc_determine_cell_shape(int *x_out, int *y_out, int *v_out,

                          int *u_out, int *N_out, htsc_param_t params, void *mem)

{

    double x = 0., y = 0., v = 0., u = 0., N = 0.;

    double frac, scaled_x;

    double ratio;

    const double  pi = 3.14159265358979323846f;

    double true_angle, lpi;

    double prev_lpi, max_lpi = 0.;

    bool use = false;

    double x_use = 0.,y_use = 0.;

    /* Go through and find the rational angle options that gets us to the

       best LPI.  Pick the one that is just over what is requested.

       That is really our limiting factor here.  Pick it and

       then figure out how much dithering we need to do to get to the proper

       number of levels.  */

    frac = tan( params.scr_ang * pi / 180.0 );

    ratio = frac * params.horiz_dpi / params.vert_dpi;

    scaled_x = params.horiz_dpi / params.vert_dpi;

    /* The minimal step is in x */

    prev_lpi = 0;

    if (ratio < 1 && ratio != 0) {

        if (params.verbose > 0) {

            PRINTF(mem, "x\ty\tu\tv\tAngle\tLPI\tLevels\n");

            PRINTF(mem, "-----------------------------------------------------------\n");

        }

        for (x = 1; x < 11; x++) {

            x_use = x;

            y=ROUND((double) x_use / ratio);

            true_angle = 180.0 * atan(((double) x_use / params.horiz_dpi) / ( (double) y / params.vert_dpi) ) / pi;

            lpi = 1.0/( sqrt( ((double) y / params.vert_dpi) * ( (double) y / params.vert_dpi) +

                                    ( (double) x_use / params.horiz_dpi) * ((double) x_use / params.horiz_dpi) ));

            v = -x_use / scaled_x;

            u = y * scaled_x;

            N = y *u - x_use * v;

            if (prev_lpi == 0) {

                prev_lpi = lpi;

                if (params.targ_lpi > lpi) {

                    EPRINTF(mem, "Warning this lpi is not achievable!\n");

                    EPRINTF(mem, "Resulting screen will be poorly quantized\n");

                    EPRINTF(mem, "or completely stochastic!\n");

                    use = true;

                }

                max_lpi = lpi;

            }

            if (prev_lpi >= params.targ_lpi && lpi < params.targ_lpi) {

                if (prev_lpi == max_lpi) {

                    EPRINTF(mem, "Notice lpi is at the maximimum level possible.\n");

                    EPRINTF(mem, "This may result in poor quantization. \n");

                }

                /* Reset these to previous x */

                x_use = x - 1;

                y=ROUND((double) x_use / ratio);

                true_angle =

                    180.0 * atan(((double) x_use / params.horiz_dpi) / ( (double) y / params.vert_dpi) ) / pi;

                lpi =

                    1.0/( sqrt( ((double) y / params.vert_dpi) * ( (double) y / params.vert_dpi) +

                                        ( (double) x_use / params.horiz_dpi) * ((double) x_use / params.horiz_dpi) ));

                v = -x_use / scaled_x;

                u = y * scaled_x;

                N = y *u - x_use * v;

                use = true;

            }

            if (use == true) {

                if (prev_lpi == max_lpi) {

                    /* Print out the final one that we will use */

                    if (params.verbose > 0)

                        PRINTF7(mem, "%3.0lf\t%3.0lf\t%3.0lf\t%3.0lf\t%3.1lf\t%3.1lf\t%3.0lf\n",

                                x_use,y,u,v,true_angle,lpi,N);

                }

                break;

            }

            if (params.verbose > 0)

                PRINTF7(mem, "%3.0lf\t%3.0lf\t%3.0lf\t%3.0lf\t%3.1lf\t%3.1lf\t%3.0lf\n",

                        x_use,y,u,v,true_angle,lpi,N);

            prev_lpi = lpi;

        }

        x = x_use;

    }

    if (ratio >= 1) {

        /* The minimal step is in y */

        if (params.verbose > 0) {

            PRINTF(mem, "x\ty\tu\tv\tAngle\tLPI\tLevels\n");

            PRINTF(mem, "-----------------------------------------------------------\n");

        }

        for (y = 1, lpi = 99999999; lpi > params.targ_lpi; y++) {

            y_use = y;

            x = ROUND(y_use * ratio);

            /* compute the true angle */

            true_angle = 180.0 * atan((x / params.horiz_dpi) / (y_use / params.vert_dpi)) / pi;

            lpi = 1.0 / sqrt( (y_use / params.vert_dpi) * (y_use / params.vert_dpi) +

                                (x / params.horiz_dpi) * (x / params.horiz_dpi));

            v = ROUND(-x / scaled_x);

            u = ROUND(y_use * scaled_x);

            N = y_use * u - x * v;

            if (prev_lpi == 0) {

                prev_lpi = lpi;

                if (params.targ_lpi > lpi) {

                    EPRINTF(mem, "Warning this lpi is not achievable!\n");

                    EPRINTF(mem, "Resulting screen will be poorly quantized\n");

                    EPRINTF(mem, "or completely stochastic!\n");

                    use = true;

                }

                max_lpi = lpi;

            }

            if (prev_lpi >= params.targ_lpi && lpi < params.targ_lpi) {

                if (prev_lpi == max_lpi) {

                    EPRINTF(mem, "Warning lpi will be slightly lower than target.\n");

                    EPRINTF(mem, "An increase will result in poor \n");

                    EPRINTF(mem, "quantization or a completely stochastic screen!\n");

                } else {

                    /* Reset these to previous x */

                    y_use = y - 1;

                    x = ROUND(y_use * ratio);

                    /* compute the true angle */

                    true_angle = 180.0 * atan((x / params.horiz_dpi) / (y_use / params.vert_dpi)) / pi;

                    lpi = 1.0 / sqrt( (y_use / params.vert_dpi) * (y_use / params.vert_dpi) +

                                        (x / params.horiz_dpi) * (x / params.horiz_dpi));

                    v = ROUND(-x / scaled_x);

                    u = ROUND(y_use * scaled_x);

                    N = y_use * u - x * v;

                }

                use = true;

            }

            if (use == true) {

                if (prev_lpi == max_lpi) {

                    /* Print out the final one that we will use */

                    PRINTF7(mem, "%3.0lf\t%3.0lf\t%3.0lf\t%3.0lf\t%3.1lf\t%3.1lf\t%3.0lf\n",

                                    x,y_use,u,v,true_angle,lpi,N);

                }

                break;

            }

            if (params.verbose > 0)

                PRINTF7(mem, "%3.0lf\t%3.0lf\t%3.0lf\t%3.0lf\t%3.1lf\t%3.1lf\t%3.0lf\n",

                        x,y_use,u,v,true_angle,lpi,N);

            prev_lpi = lpi;

        }

        y = y_use;

    }

    if (ratio == 0) {

        /* 0 degrees */

        if (scaled_x >= 1) {

            if (params.verbose > 0) {

                PRINTF(mem, "x\ty\tu\tv\tAngle\tLPI\tLevels\n");

                PRINTF(mem, "-----------------------------------------------------------\n");

            }

            for (y = 1, lpi=9999999; lpi > params.targ_lpi; y++) {

                y_use = y;

                x = ROUND( y_use * ratio );

                v = ROUND(-x / scaled_x);

                u = ROUND(y_use * scaled_x);

                N = y_use * u - x * v;

                true_angle = 0;

                lpi = 1.0/(double) sqrt( (double) ((y_use / params.vert_dpi) *

                    (y_use / params.vert_dpi) + (x / params.horiz_dpi) * (x / params.horiz_dpi)) );

                if (prev_lpi == 0) {

                    prev_lpi = lpi;

                    if (params.targ_lpi > lpi) {

                        EPRINTF(mem, "Warning this lpi is not achievable!\n");

                        EPRINTF(mem, "Resulting screen will be poorly quantized\n");

                        EPRINTF(mem, "or completely stochastic!\n");

                        use = true;

                    }

                    max_lpi = lpi;

                }

                if (prev_lpi >= params.targ_lpi && lpi < params.targ_lpi) {

                    if (prev_lpi == max_lpi) {

                        EPRINTF(mem, "Warning lpi will be slightly lower than target.\n");

                        EPRINTF(mem, "An increase will result in poor \n");

                        EPRINTF(mem, "quantization or a completely stochastic screen!\n");

                    } else {

                        /* Reset these to previous x */

                        y_use = y - 1;

                        x = ROUND( y_use * ratio );

                        v = ROUND(-x / scaled_x);

                        u = ROUND(y_use * scaled_x);

                        N = y_use * u - x * v;

                        true_angle = 0;

                        lpi = 1.0/(double) sqrt( (double) ((y_use / params.vert_dpi) *

                            (y_use / params.vert_dpi) + (x / params.horiz_dpi) * (x / params.horiz_dpi)) );

                    }

                    use = true;

                }

                if (use == true) {

                    if (prev_lpi == max_lpi) {

                        /* Print out the final one that we will use */

                        if (params.verbose > 0)

                            PRINTF7(mem, "%3.0lf\t%3.0lf\t%3.0lf\t%3.0lf\t%3.1lf\t%3.1lf\t%3.0lf\n",

                                    x,y_use,u,v,true_angle,lpi,N);

                    }

                    break;

                }

                if (params.verbose > 0)

                    PRINTF7(mem, "%3.0lf\t%3.0lf\t%3.0lf\t%3.0lf\t%3.1lf\t%3.1lf\t%3.0lf\n",

                            x,y_use,u,v,true_angle,lpi,N);

                prev_lpi = lpi;

            }

            y = y_use;

        } else {

        if (params.verbose > 0) {

            PRINTF(mem, "x\ty\tu\tv\tAngle\tLPI\tLevels\n");

            PRINTF(mem, "-----------------------------------------------------------\n");

        }

            for (x = 1; x < 11; x++) {

                x_use = x;

                y = ROUND(x_use * ratio);

                true_angle = 0;

                lpi = 1.0/( sqrt( (y / params.vert_dpi) * (y / params.vert_dpi) +

                            (x_use / params.horiz_dpi) * (x_use / params.horiz_dpi) ));

                v = ROUND( -x_use / scaled_x);

                u = ROUND( y * scaled_x);

                N = y  *u - x_use * v;

                if (prev_lpi == 0) {

                    prev_lpi = lpi;

                    if (params.targ_lpi > lpi) {

                        EPRINTF(mem, "Warning this lpi is not achievable!\n");

                        EPRINTF(mem, "Resulting screen will be poorly quantized\n");

                        EPRINTF(mem, "or completely stochastic!\n");

                        use = true;

                    }

                    max_lpi = lpi;

                }

                if (prev_lpi > params.targ_lpi && lpi < params.targ_lpi) {

                    if (prev_lpi == max_lpi) {

                        EPRINTF(mem, "Warning lpi will be slightly lower than target.\n");

                        EPRINTF(mem, "An increase will result in poor \n");

                        EPRINTF(mem, "quantization or a completely stochastic screen!\n");

                    } else {

                        /* Reset these to previous x */

                        x_use = x - 1;

                        y = ROUND(x_use * ratio);

                        true_angle = 0;

                        lpi = 1.0/( sqrt( (y / params.vert_dpi) * (y / params.vert_dpi) +

                                    (x_use / params.horiz_dpi) * (x_use / params.horiz_dpi) ));

                        v = ROUND( -x_use / scaled_x);

                        u = ROUND( y * scaled_x);

                        N = y  *u - x_use * v;

                    }

                    use = true;

                }

                if (use == true) {

                    if (prev_lpi == max_lpi) {

                        /* Print out the final one that we will use */

                        if (params.verbose > 0)

                            PRINTF7(mem, "%3.0lf\t%3.0lf\t%3.0lf\t%3.0lf\t%3.1lf\t%3.1lf\t%3.0lf\n",

                                    x_use,y,u,v,true_angle,lpi,N);

                    }

                    break;

                }

                if (params.verbose > 0)

                    PRINTF7(mem, "%3.0lf\t%3.0lf\t%3.0lf\t%3.0lf\t%3.1lf\t%3.1lf\t%3.0lf\n",

                            x_use,y,u,v,true_angle,lpi,N);

                prev_lpi = lpi;

            }

            x = x_use;

        }

    }

    *x_out = (int)x;

    *y_out = (int)y;

    *v_out = (int)v;

    *u_out = (int)u;

    *N_out = (int)N;

}

static void

htsc_create_dot_profile(htsc_dig_grid_t *dig_grid, int N, int x, int y, int u, int v,

                        double horiz_dpi, double vert_dpi, htsc_vertices_t vertices,

                        htsc_point_t *one_index, spottype_t spot_type, htsc_matrix_t trans_matrix)

{

    int done, dot_index, hole_index, count, index_x, index_y;

    htsc_dot_shape_search_t dot_search;

    int k, val_min;

    double dist;

    int j;

    htsc_vector_t vector_in, vector_out;

    double x_offset = (x % 2 == 0 ? 0.5 : 0);

    double y_offset = (y % 2 == 0 ? 0.5 : 0);

    done = 0;

    dot_index = 1;

    hole_index = N;

    count = 0;

    val_min=MIN(0,v);

    while (!done) {

        /* First perform a search for largest dot value for those remaining

           dots */

        index_x = 0;

        dot_search.index_x = 0;

        dot_search.index_y = 0;

        dot_search.norm = -100000000; /* Hopefully the dot func is not this small */

        for (k = 0; k < x + u; k++) {

            index_y = 0;

            for (j = val_min; j < y; j++) {

                if ( htsc_getpoint(dig_grid, index_x, index_y) == -1 ) {

                    /* For the spot function we want to make sure that

                       we are properly adjusted to be in the range from

                       -1 to +1.  j and k are moving in the transformed

                       (skewed/rotated) space. First untransform the value */

                    vector_in.xy[0] = k + x_offset;

                    vector_in.xy[1] = j + y_offset;

                    htsc_matrix_vector_mult(trans_matrix, vector_in,

                                            &vector_out);

                    /* And so now we are in the range 0, 1 get us to -.5, .5 */

                    vector_out.xy[0] = 2.0 * vector_out.xy[0] - 1.0;

                    vector_out.xy[1] = 2.0 * vector_out.xy[1] - 1.0;

                    dist = htsc_spot_value(spot_type, vector_out.xy[0],

                                           vector_out.xy[1]);

                    if (dist > dot_search.norm) {

                        dot_search.norm = dist;

                        dot_search.index_x = index_x;

                        dot_search.index_y = index_y;

                    }

                }

                index_y++;

            }

            index_x++;

        }

        /* Assign the index for this position */

        htsc_setpoint(dig_grid, dot_search.index_x, dot_search.index_y,

                      dot_index);

        dot_index++;

        count++;

        if (count == N) {

            done = 1;

            break;

        }

        /* The ones position for the dig_grid is located at the first dot_search

           entry.  We need this later so grab it now */

        if (count == 1) {

            one_index->x = dot_search.index_x;

            one_index->y = dot_search.index_y;

        }

        /* Now search for the closest one to a vertex (of those remaining).

           and assign the current largest index */

        index_x = 0;

        dot_search.index_x = 0;

        dot_search.index_y = 0;

        dot_search.norm = 10000000000;  /* Or this large */

        for (k = 0; k < x + u; k++) {

            index_y = 0;

            for (j = val_min; j < y; j++) {

                if ( htsc_getpoint(dig_grid, index_x, index_y) == -1 ) {

                    /* For the spot function we want to make sure that

                       we are properly adjusted to be in the range from

                       -1 to +1.  j and k are moving in the transformed

                       (skewed/rotated) space. First untransform the value */

                    vector_in.xy[0] = k + x_offset;

                    vector_in.xy[1] = j + y_offset;

                    htsc_matrix_vector_mult(trans_matrix, vector_in,

                                            &vector_out);

                    /* And so now we are in the range 0, 1 get us to -.5, .5 */

                    vector_out.xy[0] = 2.0 * vector_out.xy[0] - 1.0;

                    vector_out.xy[1] = 2.0 * vector_out.xy[1] - 1.0;

                    dist = htsc_spot_value(spot_type, vector_out.xy[0],

                                           vector_out.xy[1]);

                    if (dist < dot_search.norm) {

                        dot_search.norm = dist;

                        dot_search.index_x = index_x;

                        dot_search.index_y = index_y;

                    }

                }

                index_y++;

            }

            index_x++;

        }

        /* Assign the index for this position */

        htsc_setpoint(dig_grid, dot_search.index_x, dot_search.index_y, hole_index);

        hole_index--;

        count++;

        if (count == N) {

            done = 1;

            break;

        }

    }

}

/* This creates a mask for creating the dot shape */

static int

htsc_create_dot_mask(htsc_dig_grid_t *dot_grid, int x, int y, int u, int v,

                 double screen_angle, htsc_vertices_t vertices)

{

    int k,j,val_min,index_x,index_y;

    double slope1, slope2;

    int t1,t2,t3,t4,point_in;

    double b3, b4;

    htsc_point_t test_point;

    if (screen_angle != 0) {

        slope1 = (double) y / (double) x;

        slope2 = (double) v / (double) u;

        val_min=MIN(0,v);

        dot_grid->height = abs(val_min) + y;

        dot_grid->width = x + u;

        if (check_int_multiply(dot_grid->height, dot_grid->width, &t1) < 0)

            return -1;

        if (check_int_multiply(t1, sizeof(int), &t1) < 0)

            return -1;

        dot_grid->data =

    (int *) ALLOC(dot_grid->memory, (size_t)dot_grid->height * dot_grid->width * sizeof(int));

        if (dot_grid->data == NULL)

            return -1;

        memset(dot_grid->data,0,(size_t)dot_grid->height * dot_grid->width * sizeof(int));

        index_x = 0;

        for (k = 0; k < (x+u); k++) {

            index_y=0;

            for (j = val_min; j < y; j++) {

                test_point.x = k + 0.5;

                test_point.y = j + 0.5;

                /* test 1 and 2 */

                t1 = (slope1 * test_point.x >= test_point.y);

                t2 = (slope2 * test_point.x <= test_point.y);

                /* test 3 */

                b3 = vertices.upper_left.y - slope2 * vertices.upper_left.x;

                t3 = ((slope2 * test_point.x + b3) > test_point.y);

                /* test 4 */

                b4 = vertices.lower_right.y - slope1 * vertices.lower_right.x;

                t4=((slope1 * test_point.x + b4) < test_point.y);

                point_in = (t1 && t2 && t3 && t4);

                if (point_in) {

                    htsc_setpoint(dot_grid, index_x, index_y, -1);

                }

                index_y++;

            }

            index_x++;

        }

    } else {

        /* All points are valid */

        dot_grid->height = y;

        dot_grid->width = u;

        if (check_int_multiply(y, u, &t1) < 0)

            return -1;

        if (check_int_multiply(sizeof(int), t1, &t1) < 0)

            return -1;

        dot_grid->data = (int *) ALLOC(dot_grid->memory, (size_t)y * u * sizeof(int));

        if (dot_grid->data == NULL)