Menu
I want to simulate the flow around objects in two dimensions. Therefore I wrote a program in C which uses the Navier-Stokes equations to describe the motion of fluids. Now I came to the point where I actually want more than just placing a rectangle in the simulation domain. To draw such a rectangle I just do something like:
for(int i=start_x; i<end_x; i++)
for(int j=start_y; j<end_y; j++)
M[i][j] = 1; // barrier cell = 1
Doing this I get a nice rectangle. No surprise. But what would be an approach if I want to simulate the flow around a circle, a cross, a triangle, a wing profile or any other arbitrary polygon? Is there an easy way to draw such 2D objects in a matrix M of size m x n?
I just found an easy way to draw almost any shape I want. The answer of @Nominal Animal inspired me to find this solution. I just use a .png file and convert it to a .pgm file using the command convert picture.png picture.pgm (using Linux). In my code I only need a few more lines:
FILE *pgmFile;
pgmFile = fopen("picture.pgm", "r");
for(int i=0; i<1024; i++){
for(int j=0; j<1024; j++){
int d = fgetc(pgmFile);
if(d < 255){
M[i][j] = 1; // barrier cell = 1
}
}
}
fclose(pgmFile);
Here I use a picture of 1024 x 1024 pixels. If the value of the pixel is smaller than 255 (not white) than I set the pixel of M[i][j] to 1. Here is a result I made with the Stack Overflow logo (flux is coming from the left):
Velocity plot, Re = 20000 (Reynolds number)
935
There might be more efficient ways of doing this, but here's one way.
Define a function in C using the equation of the polygon you wish to draw. The function is defined such that it accepts a point coordinates, and returns whether the point lies inside the polygon or not. For example, for a circle, the function could accept the point (x,y), the centre (x0,y0), and the radius r, and return (x-x0)^2 + (y-y0)^2 - r^2 < 0. Let this function be f.
Determine the bounding box rectangle of the polygon, if possible, or else, the smallest rectangle you can determine which completely encloses the polygon. This will give you a rectangular matrix.
Now, iterate over the points in the rectangular matrix. For each point, call the function you previously defined. Assign the coordinate a 1 if it returns True, and 0 if it returns False. This will construct the polygon.
Suppose you want to draw a circle with centre (x0,y0), radius r, then you can use:
int f(int i, int j, int x0, int y0, int r)
{
return pow((i-x0),2) + pow((j-y0),2) - pow(r,2) < 0;
}
for(int i = x0-r; i <= x0 + r; i++)
{
for(int j = y0-r; j <= y0 + r; j++)
{
if(f(i,j,x0,y0,r))
{
M[i][j] = 1;
}
else
{
M[i][j] = 0;
}
}
}
The problem at hand boils down to rasterisation (Wikipedia); and to scan line conversion (siggraph.org) in particular.
The siggraph.org article contains detailed explanations on how to draw straight lines, circles and ellipses, and convex and concave polygons.
However, this is a problem that has already been solved a large number of times. While OP could certainly implement the necessary primitives (lines, ellipses, triangles, polygons), there is a much simpler approach.
I suggest that OP implements a simple NetPBM format reader for P5 (binary grayscale pixmap) format, and the netpbm tools (from netpbm package in Linux distributions and BSD variants; see the Netpbm home page for other systems) to convert any image to an easy-to-read PGM (P5) file, where each pixel corresponds to one element in OP's matrix.
That way, one can use e.g. Inkscape to draw the system using vector graphics, rasterize it at any size (by e.g. exporting as PNG image), convert to PGM (P5) format using the netpbm tools (pngtopnm or anytopnm, followed by ppmtopgm), and read the file. Indeed, in POSIX.1 systems (just about everywhere except in windows), one can use popen("anytopnm path-to-file | pnmtopng", "r") (or slightly more complicated two-fork() piped solution) to read any pixmap image in PGM (P5) format.
Alternatively, one could consider using e.g. the ImageMagick library to read just about any format pixmap images (JPEG, GIF, PNG etc.).
Personally, both as a developer and user (although note that I'm explicitly a non-Windows user; haven't used Microsoft products in over a decade), I would prefer the netpbm approach. The program, say mysim, would use e.g. /usr/lib/mysim/read-image shell script (or program in Windows, perhaps macs; or, if defined, the script or program defined by the MYSIM_READ_IMAGE environment variable), to read an image specified on the command line, emitting it in PGM (P5) format. The main program would simply read the output of the helper.
This way, if a user needs special handling for input files, they can trivially copy the existing script, modify it to suit their own needs, and install somewhere under their own home directory (or globally, or even replace the existing one, if it is used by all users anyway).
The program can use either popen() or fork()+execv() to execute the script, with the input filename as a command-line parameter, and read the output in the parent process to construct the initial matrix.
I prefer this approach over the image library approach for a number of reasons. First, it is more modular, allowing the user to override the image reading mechanism and manipulate it if necessary. (In my experience, such overrides are not very often needed, but when they are, they are extremely useful, and definitely overall worth it.) Second, the image processing (which in many cases is quite complex) is done in a separate process, which means all memory (for code and data) needed to read and decipher the image are released when the image has been fully read. Third, this approach follows the Unix philosophy and the KISS principle, which have a proven track record of guiding the development of robust and useful tools.
Here is an example program that reads a binary PBM, PGM, or PPM file (NetPBM P4, P5, and P6 formats, respectively) from the standard input, into a matrix structure, filling the matrix with 0 or 1 (based on the colors or grayscale values read from the image). For ease of testing, the program outputs the matrix to standard output in PGM (P5) format.
The program follows the format specification in the NetPBM manual pages (for PBM (P4), PGM (P5), and PPM (P6) formats, respectively). The Wikipedia article on NetPBM formats currently show examples with invalid comments (between the header and the data). The NetPBM manual pages state that the final header value is followed by a single whitespace character, and not a comment. (If a comment may follow the final header value, it is impossible to know whether a # (binary 0x23 = 35) in the binary data starts a comment, or is an actual data value.)
This is explicitly in the public domain, or, equivalently, licensed under the Creative Commons CC0 license. This means you are completely free to use the below code any way and anywhere you like, even in commercial projects, but that there are no guarantees: if it breaks, or breaks something, or sets your hair on fire, you get to keep all the pieces and only blame yourself.
That said, it is only lightly tested, so if you find a bug in it, let me know in the comments so I can verify and fix.
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <errno.h>
/* Matrix to read data into */
typedef struct {
int rows;
int cols;
long rowstride;
long colstride;
unsigned char *data; /* data[row*rowstride + col*colstride] */
} matrix;
#define MATRIX_INIT { 0, 0, 0, 0, NULL }
/* NetPBM (binary) formats supported */
#define PNM_PBM 4
#define PNM_PGM 5
#define PNM_PPM 6
/* Error codes from pnm_*() functions */
#define PNM_EOF -1
#define PNM_INVALID -2
#define PNM_OVERFLOW -3
/* This helper function returns the NetPBM file identifier;
PNM_PBM, PNM_PGM, PNM_PPM, or PNM_INVALID if unsupported.
*/
static int pnm_type(FILE *in)
{
/* First character must be 'P'. */
if (getc(in) != 'P')
return PNM_INVALID;
/* Second character determines the type. */
switch (getc(in)) {
case '4': return PNM_PBM;
case '5': return PNM_PGM;
case '6': return PNM_PPM;
default: return PNM_INVALID;
}
}
/* This helper function reads a number from a NetPBM header,
correctly handling comments. Since all numbers in NetPBM
headers are nonnegative, this function returns negative
when an error occurs:
-1: Premature end of input
-2: Value is too large (int overflow)
-3: Invalid input (not a NetPBM format file)
*/
static int pnm_value(FILE *in)
{
int c;
/* Skip leading whitespace and comments. */
c = getc(in);
while (c == '\t' || c == '\n' || c == '\v' ||
c == '\f' || c == '\r' || c == ' ' || c == '#')
if (c == '#') {
while (c != EOF && c != '\n')
c = getc(in);
} else
c = getc(in);
if (c == EOF)
return PNM_EOF;
if (c >= '0' && c <= '9') {
int value = 0;
while (c >= '0' && c <= '9') {
const int oldvalue = value;
value = 10*value + (c - '0');
if ((int)(value / 10) != oldvalue)
return PNM_OVERFLOW;
c = getc(in);
}
/* Do not consume the separator. */
if (c != EOF)
ungetc(c, in);
/* Success. */
return value;
}
return PNM_INVALID;
}
/* This helper function consumes the single newline
following the final value in the header.
Returns 0 if success, PNM_INVALID otherwise.
*/
static int pnm_newline(FILE *in)
{
int c;
c = getc(in);
if (c == '\r')
c = getc(in);
if (c == '\n')
return 0;
return PNM_INVALID;
}
static void pnm_matrix_free(matrix *to)
{
if (to) {
free(to->data);
to->rows = 0;
to->cols = 0;
to->rowstride = 0;
to->colstride = 0;
to->data = NULL;
}
}
static int pnm_matrix_init(matrix *to, int rows, int cols)
{
size_t cells, bytes;
if (rows < 1 || cols < 1)
return PNM_INVALID;
cells = (size_t)rows * (size_t)cols;
if ((size_t)(cells / (size_t)rows) != (size_t)cols ||
(size_t)(cells / (size_t)cols) != (size_t)rows)
return PNM_OVERFLOW;
bytes = cells * sizeof to->data[0];
if ((size_t)(bytes / sizeof to->data[0]) != cells)
return PNM_OVERFLOW;
to->data = malloc(bytes);
if (!to->data)
return PNM_OVERFLOW;
to->rows = rows;
to->cols = cols;
/* Default to a row-major data order. */
to->colstride = 1L;
to->rowstride = cols;
return 0;
}
static int pnm_p4_matrix(FILE *in, matrix *to)
{
int rows, cols, result, r, c, byte = 0;
cols = pnm_value(in);
if (cols < 1)
return PNM_INVALID;
rows = pnm_value(in);
if (rows < 1)
return PNM_INVALID;
if (pnm_newline(in))
return PNM_INVALID;
result = pnm_matrix_init(to, rows, cols);
if (result)
return result;
for (r = 0; r < rows; r++) {
const long ri = r * to->rowstride;
for (c = 0; c < cols; c++) {
const long i = ri + c * to->colstride;
switch (c & 7) {
case 0:
byte = getc(in);
if (byte == EOF) {
pnm_matrix_free(to);
return PNM_INVALID;
}
to->data[i] = !!(byte & 128);
break;
case 1:
to->data[i] = !!(byte & 64);
break;
case 2:
to->data[i] = !!(byte & 32);
break;
case 3:
to->data[i] = !!(byte & 16);
break;
case 4:
to->data[i] = !!(byte & 8);
break;
case 5:
to->data[i] = !!(byte & 4);
break;
case 6:
to->data[i] = !!(byte & 2);
break;
case 7:
to->data[i] = !!(byte & 1);
break;
}
}
}
return 0;
}
static int pnm_p5_matrix(FILE *in, matrix *to)
{
int rows, cols, max, r, c, result;
cols = pnm_value(in);
if (cols < 1)
return PNM_INVALID;
rows = pnm_value(in);
if (rows < 1)
return PNM_INVALID;
max = pnm_value(in);
if (max < 1 || max > 65535)
return PNM_INVALID;
if (pnm_newline(in))
return PNM_INVALID;
result = pnm_matrix_init(to, rows, cols);
if (result)
return result;
if (max < 256) {
const int limit = (max + 1) / 2;
int val;
for (r = 0; r < rows; r++) {
const long ri = r * to->rowstride;
for (c = 0; c < cols; c++) {
const long i = ri + c * to->colstride;
val = getc(in);
if (val == EOF) {
pnm_matrix_free(to);
return PNM_INVALID;
}
to->data[i] = (val < limit);
}
}
} else {
const int limit = (max + 1) / 2;
int val, low;
for (r = 0; r < rows; r++) {
const long ri = r * to->rowstride;
for (c = 0; c < cols; c++) {
const long i = ri + c * to->colstride;
val = getc(in);
low = getc(in);
if (val == EOF || low == EOF) {
pnm_matrix_free(to);
return PNM_INVALID;
}
val = 256*val + low;
to->data[i] = (val < limit);
}
}
}
return 0;
}
static int pnm_p6_matrix(FILE *in, matrix *to)
{
int rows, cols, max, r, c, result;
cols = pnm_value(in);
if (cols < 1)
return PNM_INVALID;
rows = pnm_value(in);
if (rows < 1)
return PNM_INVALID;
max = pnm_value(in);
if (max < 1 || max > 65535)
return PNM_INVALID;
if (pnm_newline(in))
return PNM_INVALID;
result = pnm_matrix_init(to, rows, cols);
if (result)
return result;
if (max < 256) {
const int limit = 128 * max;
int val, rval, gval, bval;
for (r = 0; r < rows; r++) {
const long ri = r * to->rowstride;
for (c = 0; c < cols; c++) {
const long i = ri + c * to->colstride;
rval = getc(in);
gval = getc(in);
bval = getc(in);
if (rval == EOF || gval == EOF || bval == EOF) {
pnm_matrix_free(to);
return PNM_INVALID;
}
val = 54 * rval
+ 183 * gval
+ 19 * bval;
to->data[i] = (val < limit);
}
}
} else {
const int limit = 128 * max;
int val, rhi, rlo, ghi, glo, bhi, blo;
for (r = 0; r < rows; r++) {
const long ri = r * to->rowstride;
for (c = 0; c < cols; c++) {
const long i = ri + c * to->colstride;
rhi = getc(in);
rlo = getc(in);
ghi = getc(in);
glo = getc(in);
bhi = getc(in);
blo = getc(in);
if (rhi == EOF || rlo == EOF ||
ghi == EOF || glo == EOF ||
bhi == EOF || blo == EOF) {
pnm_matrix_free(to);
return PNM_INVALID;
}
val = 54 * (rhi*256 + rlo)
+ 183 * (ghi*256 + glo)
+ 19 * (bhi*256 + blo);
to->data[i] = (val < limit);
}
}
}
return 0;
}
int pnm_matrix(FILE *in, matrix *to)
{
/* If the matrix is specified, initialize it. */
if (to) {
to->rows = 0L;
to->cols = 0L;
to->rowstride = 0L;
to->colstride = 0L;
to->data = NULL;
}
/* Sanity checks on parameters. */
if (!to || !in || ferror(in))
return PNM_INVALID;
switch (pnm_type(in)) {
case PNM_PBM: return pnm_p4_matrix(in, to);
case PNM_PGM: return pnm_p5_matrix(in, to);
case PNM_PPM: return pnm_p6_matrix(in, to);
default: return PNM_INVALID;
}
}
int main(void)
{
int r, c;
matrix m = MATRIX_INIT;
if (pnm_matrix(stdin, &m)) {
fprintf(stderr, "Cannot parse standard input.\n");
return EXIT_FAILURE;
}
fprintf(stderr, "Read %d rows, %d columns, from standard input.\n", m.rows, m.cols);
/* For ease of debugging, we output the matrix as a PGM file. */
printf("P5\n%d %d\n255\n", m.cols, m.rows);
for (r = 0; r < m.rows; r++)
for (c = 0; c < m.cols; c++)
if (m.data[r * m.rowstride + c * m.colstride] == 0)
putchar(255); /* White */
else
putchar(0); /* Black */
return EXIT_SUCCESS;
}
Note that I did not verify whether the bit/grayscale/color conversion is the right way with respect to how OP intends to use the matrix. (That is, whether "white" or light colours should yield a 0 or a 1 in the matrix.) If you need to invert it for PBM images, use !(byte & NUMBER) instead. If you need to invert it for PGM or PPM images, use (val >= limit) instead.
The program should be valid C (even down to C89), and compile on any architecture. On silly architectures like Windows, you might have to open/reopen the standard input in "binary mode" (include b in the fopen() flags), as they otherwise may mangle the input.
On Linux, I compiled and tested the program (example.c) with
gcc -Wall -O2 example.c -o example
./example < inputfile.pbm > result-pbm.pgm
./example < inputfile.pgm > result-pgm.pgm
./example < inputfile.ppm > result-ppm.pgm
I prefer to use 'sin' and 'cos' to make a circle. If you want some special shape like Oval. You can use the 'fac_specialX' and 'fac_specialY' make it different. If 'fac_specialX' and 'fac_specialY' not a fixed value(maybe be changed each time in loop), can make the shape more special.(or just try to modify a part of the circle array)
int r=10;// radius
int x0=25,y0=25; // center
int circle_points = 300; // accuracy --> higher cause better quality but slow
int circleX[circle_points]; // circle array
int circleY[circle_points]; // circle array
// #define PI 3.1415926
double fac_angle = ( 2*PI ) / circle_points;
// normal circle : fac_specialX & fac_specialY set 1
// Oval : fac_specialX --> higher cause longer in X
// fac_specialY --> higher cause longer in Y
double fac_specialX = 0.5;
double fac_specialY = 1.5;
// Calculate the coordinates
for(int i=0 ; i<circle_points ; i++) {
// #include <math.h> ->> sin cos
circleX[i] = x0 + r*sin( (i+1)*fac_angle )*fac_specialX;
circleY[i] = y0 + r*cos( (i+1)*fac_angle )*fac_specialY;
// set the ponts in M array
M[ circleY[i] ][ circleX[i] ] = 1;
}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
