PRIR FRAKTALE

/************************************************************
 * -- NIE USUWAJ TEJ INFORMACJI Z PROGRAMU ---------------- *
 ************************************************************
 * -- Program powstał na bazie kodu źródłowego ------------ *
 * -- udostępnionego studentom na potrzeby przedmiotu ----- *
 * -- Programowanie Równoległe i Rozproszone -------------- *
 * -- Copyright (c) 2013 Politechnika Śląska w Gliwicach -- *
 * -- Wydział Elektryczny --------------------------------- *
 * -- Radosław Sokół, Marcin Połomski --------------------- *
 ************************************************************/

#include <assert.h>
#include <windows.h>
#include <stdio.h>
#include <exception>
#include <omp.h>

#include "Complex.h"

TCHAR NazwaAplikacji[] = TEXT("FRAKTALE");
TCHAR NazwaKlasy[] = TEXT("OKNOGLOWNE");

const int ImageWidth = 1024;
const int ImageHeight = 768;
char *Obraz = NULL;

unsigned char NrAlgorytmu = 0;
double CzasPrzetwarzania = 0.0;
TCHAR tAlgorytm[64];
unsigned char tadl = 0;
const int L_WATKOW = 8;

struct RGB { unsigned char B, G, R; };

inline unsigned int RoundUp4(const unsigned int i)
{
    return (i+3)&(~3);
}

void WyczyscObraz(RGB* Raster, int Width, int Height)
{
    RGB *Punkt;
    for (int y = 0; y < Height; ++y)
    {
        Punkt = Raster + RoundUp4(y * Width * sizeof(RGB)) / sizeof(RGB);
        for (int x = 0; x < Width; ++x)
        {
            Punkt->R = 0;
            Punkt->G = 0;
            Punkt->B = 0;
            ++Punkt;
        }
    }
}

// Mandelbrot =====================================================================================

int IterateMandelbrot(Complex P, int maxIter)
{
    Complex z = P;
    int cnt = 0;
    while((z.x * z.x + z.y * z.y <= 4)&& (cnt < maxIter)){
        z = ComplexSquare(z);
        z.x += P.x;
        z.y += P.y;
        cnt++;
    }
    return cnt;
}

void Mandelbrot(RGB* Raster, int Width, int Height, int MaxIter, float RealMin, float RealMax, float ImagMin, float ImagMax)
{
    RGB *Punkt;
    float realInc = (RealMax - RealMin) / Width;
    float imagInc = (ImagMax - ImagMin) / Height;
    Complex z;
    int x, y;
    int cnt;

    for (y = 0, z.y = ImagMin; y < Height; ++y, z.y += imagInc)
    {
        Punkt = Raster + RoundUp4(y * Width * sizeof(RGB)) / sizeof(RGB);

        for (x = 0, z.x = RealMin; x < Width; ++x, z.x += realInc)
        {
            cnt = IterateMandelbrot(z, MaxIter);

            if(cnt == MaxIter)
            {
                Punkt->R = 255;
                Punkt->G = 255;
                Punkt->B = 255;
            }
            else
            {
                Punkt->R = 0;
                Punkt->G = cnt;
                Punkt->B = 0;
            }

            ++Punkt;
        }
    }
}

void MandelbrotParallel(RGB* Raster, int Width, int Height, int MaxIter, float RealMin, float RealMax, float ImagMin, float ImagMax)
{
    RGB *Punkt;
    float realInc = (RealMax - RealMin) / Width;
    float imagInc = (ImagMax - ImagMin) / Height;
    Complex z;
    int x, y;
    int cnt;

#pragma omp parallel for private(y, Punkt, x, cnt, z) num_threads(L_WATKOW)
        for (y = 0 ; y < Height; ++y )
        {
        z.y=ImagMin+imagInc*y;
            Punkt = Raster + RoundUp4(y * Width * sizeof(RGB)) / sizeof(RGB);
            for (x = 0, z.x = RealMin; x < Width; ++x, z.x += realInc)
            {
                cnt = IterateMandelbrot(z, MaxIter);

                if(cnt == MaxIter)
                {
                    Punkt->R = 255;
                    Punkt->G = 255;
                    Punkt->B = 255;
                }
                else
                {
                    Punkt->R = 0;
                    Punkt->G = cnt;
                    Punkt->B = 0;
                }
                ++Punkt;
            }
        }
}

// Julia ==========================================================================================

int IterateJulia(Complex zInit, int maxIter, Complex c)
{
    Complex z = zInit;
    int cnt = 0;
    while((z.x * z.x + z.y * z.y <= 4)&& (cnt < maxIter)){
        z = ComplexSquare(z);
        z.x += c.x;
        z.y += c.y;
        cnt++;
    }
    return cnt;
}

void Julia(RGB* Raster, int Width, int Height, int MaxIter, float RealMin, float RealMax, float ImagMin, float ImagMax, Complex c)
{
    RGB *Punkt;
    float realInc = (RealMax - RealMin) / Width;
    float imagInc = (ImagMax - ImagMin) / Height;
    Complex z;
    int x, y;
    int cnt;

    for (y = 0, z.y = ImagMin; y < Height; ++y, z.y += imagInc)
    {
        Punkt = Raster + RoundUp4(y * Width * sizeof(RGB)) / sizeof(RGB);

        for (x = 0, z.x = RealMin; x < Width; ++x, z.x += realInc)
        {
            cnt = IterateJulia(z, MaxIter, c);

            if(cnt == MaxIter)
            {
                Punkt->R = 255;
                Punkt->G = 255;
                Punkt->B = 255;
            }
            else
            {
                Punkt->R = 0;
                Punkt->G = cnt;
                Punkt->B = 0;
            }

            ++Punkt;
        }
    }
}

void JuliaParallel(RGB* Raster, int Width, int Height, int MaxIter, float RealMin, float RealMax, float ImagMin, float ImagMax, Complex c)
{
    RGB *Punkt;
    float realInc = (RealMax - RealMin) / Width;
    float imagInc = (ImagMax - ImagMin) / Height;
    Complex z;
    int x, y;
    int cnt;

    #pragma omp parallel for private(x, y, z, cnt, Punkt) num_threads(L_WATKOW)
    for (y = 0; y < Height; ++y)
    {
        z.y = ImagMin+imagInc*y;
        Punkt = Raster + RoundUp4(y * Width * sizeof(RGB)) / sizeof(RGB);

        for (x = 0, z.x = RealMin; x < Width; ++x, z.x += realInc)
        {
            cnt = IterateJulia(z, MaxIter, c);

            if(cnt == MaxIter)
            {
                Punkt->R = 255;
                Punkt->G = 255;
                Punkt->B = 255;
            }
            else
            {
                Punkt->R = 0;
                Punkt->G = cnt;
                Punkt->B = 0;
            }
            ++Punkt;
        }
        z.y += imagInc;
    }
}


static void GenerujFraktal()
{
    assert(sizeof(RGB) == 3);
    BITMAPFILEHEADER *bfh = reinterpret_cast<BITMAPFILEHEADER*>(Obraz);
    BITMAPINFO *bi = reinterpret_cast<BITMAPINFO*>(Obraz + sizeof(BITMAPFILEHEADER));
    RGB *Raster = reinterpret_cast<RGB*>(Obraz + bfh->bfOffBits);

    int Width = bi->bmiHeader.biWidth;
    int Height = bi->bmiHeader.biHeight;

    // Stała C1 uzywane podczas generowania zbioru Julii
    Complex c1;
    c1.x = 0.281f; c1.y = 0.011f;

    // Stała C2 uzywane podczas generowania zbioru Julii
    Complex c2;
    c2.x = -0.810f; c2.y = 0.179f;

    switch (NrAlgorytmu)
    {
        case 0:
            // Mandelbrot Sekwencyjny ----------------------------------------------------------------------
            tadl = wsprintf(tAlgorytm, TEXT(". . ."));
            WyczyscObraz(Raster, Width, Height);
            break;

        case 1:
            tadl = wsprintf(tAlgorytm, TEXT("Mandlebrot Sekwencyjny"));
            Mandelbrot(Raster, Width, Height, 500, -2.0f, 0.8f, -1.0f, 1.0f);
            break;

        case 2:
            tadl = wsprintf(tAlgorytm, TEXT("Julia Sekwencyjny dla C1"));
            Julia(Raster, Width, Height, 500, -1.0f, 1.0f, -1.0f, 1.0f, c1);
            break;

        case 3:
            tadl = wsprintf(tAlgorytm, TEXT("Julia Sekwencyjny dla C2"));
            Julia(Raster, Width, Height, 500, -1.0f, 1.0f, -0.8f, 0.8f, c2);
            break;

        case 4:
            {
                tadl = wsprintf(tAlgorytm, TEXT("Mandelbrot Zrównoleglony"));
                MandelbrotParallel(Raster, Width, Height, 500, -2.0f, 0.8f, -1.0f, 1.0f);
            }
            break;

        case 5:
            {
                tadl = wsprintf(tAlgorytm, TEXT("Julia Zrównoleglony dla C1"));
                JuliaParallel(Raster, Width, Height, 500, -1.0f, 1.0f, -1.0f, 1.0f, c1);
            }
            break;

        case 6:
            {
                tadl = wsprintf(tAlgorytm, TEXT("Julia Zrównoleglony dla C2"));
                JuliaParallel(Raster, Width, Height, 500, -1.0f, 1.0f, -0.8f, 0.8f, c2);
            }
            break;

        default:
            tadl = wsprintf(tAlgorytm, TEXT(". . ."));
            WyczyscObraz(Raster, Width, Height);
            break;
    }
}

static void Tworz(const HWND Okno)
{
    BITMAPFILEHEADER   bmfHeader;
    BITMAPINFOHEADER   bi;

    bi.biSize = sizeof(BITMAPINFOHEADER);
    bi.biWidth = ImageWidth;
    bi.biHeight = ImageHeight;
    bi.biPlanes = 1;
    bi.biBitCount = 24;
    bi.biCompression = BI_RGB;
    bi.biSizeImage = 0;
    bi.biXPelsPerMeter = 0;
    bi.biYPelsPerMeter = 0;
    bi.biClrUsed = 0;
    bi.biClrImportant = 0;

    DWORD dwBmpSize = ImageWidth * ImageHeight * bi.biBitCount/8;
    HANDLE hDIB = GlobalAlloc(GHND, dwBmpSize);
    char *lpbitmap = (char *)GlobalLock(hDIB);

    DWORD dwSizeOfDIB = dwBmpSize + sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER);
    bmfHeader.bfOffBits = (DWORD)sizeof(BITMAPFILEHEADER) + (DWORD)sizeof(BITMAPINFOHEADER);
    bmfHeader.bfSize = dwSizeOfDIB;
    bmfHeader.bfType = 0x4D42;

    Obraz = new char[dwSizeOfDIB];
    CopyMemory(Obraz, &bmfHeader, sizeof(BITMAPFILEHEADER));
    CopyMemory(Obraz + sizeof(BITMAPFILEHEADER), &bi, sizeof(BITMAPINFOHEADER));
    CopyMemory(Obraz + sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER), lpbitmap, dwBmpSize);

    GlobalUnlock(hDIB);
    GlobalFree(hDIB);
}

static void Rysuj(const HWND Okno)
{
    RECT Rozmiar;
    TCHAR Tekst[128];
    PAINTSTRUCT PS;
    unsigned int dl;
    GetClientRect(Okno, &Rozmiar);
    const HDC DC = BeginPaint(Okno, &PS);
    if (Obraz != 0) {
        BITMAPFILEHEADER *bfh = reinterpret_cast<BITMAPFILEHEADER*>(Obraz);
        BITMAPINFO *bi = reinterpret_cast<BITMAPINFO*>(Obraz + sizeof(BITMAPFILEHEADER));
        SetStretchBltMode(DC, HALFTONE);
        StretchDIBits(DC, 0, 0, Rozmiar.right, Rozmiar.bottom, 0, 0, bi->bmiHeader.biWidth, bi->bmiHeader.biHeight,
            static_cast<void*>(Obraz + bfh->bfOffBits), bi, DIB_RGB_COLORS, SRCCOPY);
        dl = swprintf_s(Tekst, TEXT("Czas obliczeń: %.2lf s"), CzasPrzetwarzania);
        TextOut(DC, 5, 5, Tekst, dl);
        TextOut(DC, 5, 25, tAlgorytm, tadl);
    }
    EndPaint(Okno, &PS);
}

static void Znak(const HWND Okno, const TCHAR Znak)
{
    if (Znak >= TEXT('0') && Znak <= TEXT('9')) {
        SYSTEMTIME S1, S2;
        NrAlgorytmu = Znak - TEXT('0');
        SetWindowText(Okno, TEXT("Przetwarzanie..."));
        GetSystemTime(&S1);
        if (NrAlgorytmu >= 0 && NrAlgorytmu <= 9)
        {
            GenerujFraktal();
        }
        GetSystemTime(&S2);
        CzasPrzetwarzania = 3600.0 * (S2.wHour - S1.wHour) + 60.0 * (S2.wMinute - S1.wMinute)
            + (S2.wSecond - S1.wSecond) + 0.001 * (S2.wMilliseconds - S1.wMilliseconds);
        InvalidateRect(Okno, NULL, FALSE);
        SetWindowText(Okno, NazwaAplikacji);
        UpdateWindow(Okno);
    }

    if(Znak == 27)
    {
        DestroyWindow(Okno);
    }
}

static LRESULT CALLBACK FunkcjaOkienkowa(HWND Okno, UINT Komunikat, WPARAM wParam, LPARAM lParam)
{
    switch (Komunikat) {
        case WM_CREATE:
            Tworz(Okno);
            tadl = wsprintf(tAlgorytm, TEXT("..."));
            break;
        case WM_CHAR:
            Znak(Okno, static_cast<TCHAR>(wParam));
            break;
        case WM_PAINT:
            Rysuj(Okno);
            break;
        case WM_DESTROY:
            delete [] Obraz;
            PostQuitMessage(0);
            break;
        default:
            return DefWindowProc(Okno, Komunikat, wParam, lParam);
    }
    return 0;
}

static bool RejestrujKlasy()
{
    WNDCLASSEX wc;
    wc.cbSize = sizeof(WNDCLASSEX);
    wc.cbClsExtra = wc.cbWndExtra = 0;
    wc.hbrBackground = (HBRUSH) (1 + COLOR_WINDOW);
    wc.hCursor = LoadCursor(NULL, IDC_ARROW);
    wc.hIcon = LoadIcon(NULL, IDI_APPLICATION);
    wc.hIconSm = LoadIcon(NULL, IDI_APPLICATION);
    wc.hInstance = GetModuleHandle(NULL);
    wc.lpfnWndProc = &FunkcjaOkienkowa;
    wc.lpszClassName = NazwaKlasy;
    wc.lpszMenuName = NULL;
    wc.style = CS_HREDRAW | CS_VREDRAW;
    return (RegisterClassEx(&wc) != 0);
}

static void WyrejestrujKlasy()
{
    UnregisterClass(NazwaKlasy, GetModuleHandle(NULL));
}

int WINAPI WinMain(HINSTANCE Instancja, HINSTANCE Poprzednia, LPSTR Parametry, int Widocznosc)
{
    // Zarejestruj klasę. Protestuj, jeżeli wystąpił błąd.
    if (!RejestrujKlasy()) {
        MessageBox(NULL, TEXT("Nie udało się zarejestrować klasy okna!"),
            NazwaAplikacji, MB_ICONSTOP | MB_OK);
        return 1;
    }
    // Stwórz główne okno. Również protestuj, jeżeli wystąpił błąd.
    HWND GlowneOkno = CreateWindowEx(WS_EX_APPWINDOW | WS_EX_CLIENTEDGE,
        NazwaKlasy, NazwaAplikacji, WS_OVERLAPPEDWINDOW,
        CW_USEDEFAULT, CW_USEDEFAULT, CW_USEDEFAULT, CW_USEDEFAULT,
        NULL, NULL, Instancja, NULL);
    if (GlowneOkno == NULL) {
        MessageBox(NULL, TEXT("Nie udało się stworzyć głównego okna!"),
            NazwaAplikacji, MB_ICONSTOP | MB_OK);
        return 2;
    }
    // Wyświetl i uaktualnij nowo stworzone okno.
    ShowWindow(GlowneOkno, Widocznosc);
    UpdateWindow(GlowneOkno);
    // Główna pętla komunikatów wątku.
    MSG Komunikat;
    while (GetMessage(&Komunikat, NULL, 0, 0) > 0) {
        TranslateMessage(&Komunikat);
        DispatchMessage(&Komunikat);
    }
    // Zwolnij pamięć klas i zakończ proces.
    WyrejestrujKlasy();
    return static_cast<int>(Komunikat.wParam);
}