215 lines
7.9 KiB
C
215 lines
7.9 KiB
C
#include "tricorn_V3.h"
|
|
#include "color_table.h"
|
|
#include "converge_helper.h"
|
|
|
|
static void* compute(void * args);
|
|
static size_t get_thread_count();
|
|
|
|
// Argument structs passed to the threads for the caluclations
|
|
typedef struct{
|
|
unsigned n;
|
|
_Bool color;
|
|
_Bool flipped;
|
|
float res;
|
|
size_t width;
|
|
unsigned char* img;
|
|
|
|
} const_args;
|
|
typedef struct {
|
|
float start_r;
|
|
float start_i;
|
|
size_t start_row;
|
|
size_t stop_row;
|
|
const_args* consts;
|
|
} compute_args;
|
|
|
|
|
|
|
|
|
|
/**
|
|
* @brief Generates a Tricorn fractal image using an optimized SIMD + Threading implementation.
|
|
* @param start The starting complex coordinate (bottom-left or top-left depending on height sign).
|
|
* @param width The width of the target image in pixels.
|
|
* @param height The height of the target image. A negative height reverses the Y-axis calculation.
|
|
* @param res The step size (resolution) per pixel in the complex plane.
|
|
* @param n The maximum number of mathematical iterations per pixel.
|
|
* @param color If true, generates a 24-bit RGB image. If false, writes a grayscale image.
|
|
* @param img Pointer to the pre-allocated memory buffer where the bitmap data will be stored.
|
|
*/
|
|
void tricorn_V3(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img){
|
|
|
|
const size_t min_rows_per_thread = 10;
|
|
const _Bool flipped = height < 0;
|
|
const size_t abs_height = flipped ? -height : height;
|
|
|
|
|
|
// Calculating the threadcount so that it is limited to a reasonable distribution where thread starting doesnt outway calculation cost
|
|
size_t thread_count = get_thread_count();
|
|
if(abs_height >= min_rows_per_thread)
|
|
while((abs_height/thread_count) < min_rows_per_thread) thread_count--;
|
|
else thread_count = 1;
|
|
|
|
const size_t y_step = (abs_height+thread_count-1)/thread_count;
|
|
const float start_r = crealf(start);
|
|
const float directional_res = (flipped ? -1.0f: 1.0f)*res;
|
|
const float start_i = height > 0 ? cimag(start) : cimag(start) + ((abs_height - 1) * res);
|
|
|
|
const_args* consts = malloc(sizeof(*consts));
|
|
consts->width = width;
|
|
consts->flipped = flipped;
|
|
consts->color = color;
|
|
consts->res = res;
|
|
consts->img = img;
|
|
consts->n = n;
|
|
|
|
pthread_t threads[thread_count];
|
|
for(size_t i = 0; i < thread_count; i++){
|
|
compute_args* args = malloc(sizeof(*args));
|
|
args->start_r =start_r;
|
|
args->start_i =start_i+i*directional_res*y_step;
|
|
args->start_row = i*y_step;
|
|
args->stop_row = abs_height < (args->start_row+y_step) ? abs_height: (args->start_row+y_step);
|
|
args->consts = consts;
|
|
if(pthread_create(&threads[i],NULL,compute,args)){
|
|
compute(args);
|
|
free(args);
|
|
args = NULL;
|
|
}
|
|
}
|
|
for(size_t i = 0; i < thread_count; i++){
|
|
pthread_join(threads[i],NULL);
|
|
}
|
|
|
|
free(consts);
|
|
consts = NULL;
|
|
return;
|
|
}
|
|
|
|
static void* compute(void * vargs){
|
|
|
|
compute_args *args = vargs;
|
|
const_args *c = args->consts;
|
|
vargs = NULL;
|
|
|
|
|
|
const int bytes_per_color = c->color ? 3 : 1;
|
|
// Absolute height for iteration purposes
|
|
const size_t simd_width = (c->width / 4) * 4;
|
|
|
|
float c_start_real = args->start_r;
|
|
__m128 four = _mm_set1_ps(4.0f);
|
|
__m128 minus_two = _mm_set1_ps(-2.0f);
|
|
|
|
float c_imag = args->start_i;
|
|
__m128 c_imag_vector = _mm_set1_ps(c_imag); // Initialize imagiinary vector which will be acumulated
|
|
|
|
float four_res = 4 * c->res;
|
|
float vertical_step = (c->flipped? -1.0f: 1.0f) * c->res;
|
|
__m128 vertical_step_vector = _mm_set1_ps(vertical_step);
|
|
|
|
for(size_t y = args->start_row; y < args->stop_row;y++) {
|
|
size_t x = 0;
|
|
// Initialize real vector which will be accumulated
|
|
__m128 c_real_vector = _mm_setr_ps(c_start_real,c_start_real + 1.0f*c->res,c_start_real + 2.0f*c->res, c_start_real + 3.0f*c->res);
|
|
for(; x < simd_width;x += 4) { //SIMD Loop
|
|
// Initialize start values for the calculation
|
|
__m128 real = _mm_setzero_ps();
|
|
__m128 imag = _mm_setzero_ps();
|
|
__m128i iterations = _mm_setzero_si128();
|
|
__m128i active_mask = _mm_set1_epi32(-1); // Active mask will track which pixels have digerved - 0 for diverged, 1 - for contained, at start 0xFFFF
|
|
|
|
for(unsigned i = 0; i < c->n;i++) {
|
|
// Calculate squared values of real and imaginary
|
|
__m128 real2 = _mm_mul_ps(real,real);
|
|
__m128 imag2 = _mm_mul_ps(imag,imag);
|
|
|
|
// Detect which have diverged
|
|
__m128 unbounded = _mm_cmple_ps(_mm_add_ps(real2,imag2),four);
|
|
__m128i unbounded_i = _mm_castps_si128(unbounded);
|
|
|
|
// Update the active mask
|
|
active_mask = _mm_and_si128(active_mask,unbounded_i);
|
|
|
|
// If all have diverged break
|
|
if(_mm_test_all_zeros(active_mask,active_mask)) {
|
|
break;
|
|
}
|
|
|
|
// Add the iterations for each pixel only if it hasn't diverged
|
|
iterations = _mm_sub_epi32(iterations,active_mask);
|
|
|
|
imag = _mm_add_ps(c_imag_vector,_mm_mul_ps(minus_two, _mm_mul_ps(real,imag)));
|
|
real = _mm_add_ps(_mm_sub_ps(real2,imag2),c_real_vector);
|
|
}
|
|
|
|
// Calculate the value for the pixel
|
|
// value = 255 - 4 * (i % 64)
|
|
__m128i mod64 = _mm_and_si128(iterations, _mm_set1_epi32(63)); // mod is realised with a bitmask
|
|
__m128i times4 = _mm_slli_epi32(mod64, 2); // then multiplication by 4 is shift by 2
|
|
__m128i values = _mm_sub_epi32(_mm_set1_epi32(255), times4); // subtraction from 255
|
|
|
|
// pixels which have not escaped need to be colored black so value = 0
|
|
__m128i max_iter = _mm_set1_epi32(c->n);
|
|
__m128i is_max = _mm_cmpeq_epi32(iterations, max_iter);
|
|
values = _mm_andnot_si128(is_max, values);
|
|
|
|
// Convert vector to an array of ints
|
|
unsigned int value_array[4];
|
|
_mm_storeu_si128((__m128i*)value_array,values);
|
|
|
|
// Write the 4 values in the image buffer
|
|
for(unsigned i = 0 ; i < 4;i++) {
|
|
unsigned int value = value_array[i];
|
|
size_t position = (y*c->width + (x+i)) * bytes_per_color;
|
|
|
|
if(!c->color) {
|
|
// If grayscale just take value
|
|
c->img[position] = value;
|
|
}
|
|
else {
|
|
// If colored use the lookup table to get the RGB values
|
|
c->img[position] = color_lookup_table[value][0];
|
|
c->img[position + 1] = color_lookup_table[value][1];
|
|
c->img[position + 2] = color_lookup_table[value][2];
|
|
}
|
|
}
|
|
// Accumulate the real part
|
|
c_real_vector = _mm_add_ps(c_real_vector,_mm_set1_ps(four_res));
|
|
}
|
|
|
|
// Follow up loop for the remaining pixels
|
|
// Calculates them scalarly
|
|
for(;x < c->width; x++) {
|
|
float c_real = c_start_real + (x*c->res);
|
|
unsigned int value = converge(c_real,c_imag,c->n);
|
|
|
|
size_t position = (y*c->width + x) * bytes_per_color;
|
|
|
|
if(!c->color) {
|
|
c->img[position] = value;
|
|
}
|
|
else {
|
|
c->img[position] = color_lookup_table[value][0];
|
|
c->img[position + 1] = color_lookup_table[value][1];
|
|
c->img[position + 2] = color_lookup_table[value][2];
|
|
}
|
|
}
|
|
// Accumulate the imaginary part
|
|
c_imag += vertical_step;
|
|
c_imag_vector = _mm_add_ps(c_imag_vector,vertical_step_vector);
|
|
}
|
|
free(args);
|
|
args = NULL;
|
|
return NULL;
|
|
}
|
|
|
|
|
|
static inline size_t get_thread_count(){
|
|
#ifdef WIN32
|
|
SYSTEM_INFO sysinfo;
|
|
GetSystemInfo(&sysinfo);
|
|
return sysinfo.dwNumberOfProcessors;
|
|
#else
|
|
return sysconf(_SC_NPROCESSORS_ONLN);
|
|
#endif
|
|
} |