Commit b1c82686 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Merge branch 'registration' into 'main' 

registration app

See merge request !1
parents 6caecf7e 9385775f

Dockerfile

0 → 100644
+10 −0
Original line number Original line Diff line number Diff line 
FROM node:16-alpine

RUN apk add --no-cache python3

RUN apk add --no-cache make

RUN apk add --no-cache g++

RUN apk add --no-cache gcc

WORKDIR /app

COPY . .

RUN npm install

EXPOSE 8080

CMD ["node","app.js"]
 
 No newline at end of file

addon/common.h

0 → 100644
+9 −0
Original line number Original line Diff line number Diff line 
#pragma once



#include <vector>



using Array1D = std::vector<double>;

using Array2D = std::vector<Array1D>;



using Array1D_int = std::vector<int>;

using Array2D_int = std::vector<Array1D_int>;

addon/distFunkce.h

0 → 100644
+101 −0
Original line number Original line Diff line number Diff line 
#pragma once



#include <cmath>



#include "common.h"



void distance( Array2D& u0, int n1, int n2, int N1, int N2, int rozsah, double h1, double h2, int pocetIteraci)

{

	float e=0.0001;

	float t= 5.0e-7;

	float rezid=100;



	// N1 radku na N2 sloupcu

	Array2D u(N1, Array1D(N2));

	Array2D v(N1, Array1D(N2));



	for (int i=n1; i<=n1+rozsah; i++){

		for (int j=n2; j<=n2+rozsah; j++){

			u[i][j] = u0[i][j];

		}

	}





	int n=0;

	while(n<pocetIteraci){

		n=n+1;



		for (int i=n1+1; i<n1 + rozsah; i++)

		for (int j=n2+1; j<n2 + rozsah; j++)

		{

			double clen1 = 0.0;

			double clen2 = 0.0;

			double clen3 = 0.0;

			double clen4 = 0.0;



			double S = u0[i][j]/(sqrt(e*e + u0[i][j]*u0[i][j]));



			if(S>0){

				clen1 = ((u[i][j] - u[i][j-1])/h1);

				if (clen1 <0){ clen1 = 0.0; };

				clen2 = ((u[i][j+1] - u[i][j])/h1);

				if (clen2 >0){ clen2 = 0.0; };

				if(clen1 < -clen2){clen1 = clen2;}



				clen3 = ((u[i][j] - u[i-1][j])/h2);

				if (clen3 <0){ clen3 = 0.0; };

				clen4 = ((u[i+1][j] - u[i][j])/h2);

				if (clen4 >0){ clen4 = 0.0; };

				if(clen3 < -clen4){clen3 = clen4;}

			}

			else{

				clen1 = ((u[i][j+1] - u[i][j])/h1);

				if (clen1 <0){ clen1 = 0.0; };

				clen2 = ((u[i][j] - u[i][j-1])/h1);

				if (clen2 >0){ clen2 = 0.0; };

				if(clen1 < -clen2){clen1 = clen2;}



				clen3 = ((u[i+1][j] - u[i][j])/h2);

				if (clen3 <0){ clen3 = 0.0; };

				clen4 = ((u[i][j] - u[i-1][j])/h2);

				if (clen4 >0){ clen4 = 0.0; };

				if(clen3 < -clen4){clen3 = clen4;}

			}



			v[i][j] = u[i][j] + t*S*( 1.0 - sqrt(clen1*clen1 + clen3*clen3));

		}





		for (int j=n2+1; j<=n2 + rozsah; j++){	v[n1][j]=u[n1+1][j];}

		for (int j=n2+1; j<=n2 + rozsah; j++){	v[n1 + rozsah][j]=u[n1 + rozsah-1][j];}



		for (int i=n1+1; i<=n1 + rozsah; i++){	v[i][n2]=u[i][n2+1];}

		for (int i=n1+1; i<=n1 + rozsah; i++){	v[i][n2 + rozsah]=u[i][n2 + rozsah-1];}



		v[n1][n2]=u[n1][n2+1];

		v[n1][n2 + rozsah]=u[n1][n2 + rozsah-1];



		v[n1 + rozsah][n2]=u[n1 + rozsah -1][n2+1];

		v[n1 + rozsah][n2 + rozsah]=u[n1 + rozsah-1][n2 + rozsah-1];





		rezid=0.0;

		for (int i=n1; i<=n1 + rozsah; i++){

			for (int j=n2; j<=n2 + rozsah; j++){

				rezid = rezid +(v[i][j] - u[i][j])*(v[i][j] - u[i][j]);

				u[i][j] = v[i][j];

			}

		}

		rezid = sqrt(rezid);



		// vypis z GUI nema smysl

		//if(n%10000==0){	cout << rezid << endl;}

	}



	for (int i=n1; i<=n1+rozsah; i++){

		for (int j=n2; j<=n2+rozsah; j++){

			u0[i][j] = v[i][j];

		}

	}



}

addon/filtr.h

0 → 100644
+358 −0
Original line number Original line Diff line number Diff line 
#pragma once



#include <cmath>

#include "common.h"



int max_obr(Array2D& u, int N1, int N2, int* mi, int* mj)

{

    int maxV = u[0][0];



    for (int i=0; i<N1; i++)

    {

        for (int j=0; j<N2; j++)

        {

             if(maxV<u[i][j])

             {

                 maxV = u[i][j];

                 *mi = i;

                 *mj = j;

             }

        }

    }

    return maxV;

}





void ekvalizace(Array2D& u, int N1, int N2)

{

    // EKVALIZACE : priprava



    

    

    int maxVal,mi,mj;



    

    for(int i =0;i<3;i++)

    {

        maxVal = max_obr(u,N1,N2,&mi,&mj);

        u[mi][mj] = 0;

    }

    

    if(maxVal<=300)

    {

    

        // MSVC does not support variable-length arrays

//        int sum[maxVal] = {};

//        double ratio[maxVal];

        Array1D_int sum(maxVal);

        Array1D ratio(maxVal);



        int cumsum = 0;



        //Ekvalizace:

        

        for (int i=0; i<N1; i++)

        {

            for (int j=0; j<N2; j++)

            {

                sum[int(u[i][j])] ++;

            }

        }

        

        for (int i=0; i<= maxVal; i++)

        {

            cumsum = cumsum + sum[i];

            ratio[i] = float(cumsum)/float(N1*N2);

        }

        

        for (int i=0; i<N1; i++)

        {

            for (int j=0; j<N2; j++)

            {

                u[i][j] = ceil(ratio[int(u[i][j])]*maxVal);

            }

        }

        

    }

    else

    {

        

        maxVal = 300;

        

        // MSVC does not support variable-length arrays

//        int sum[maxVal] = {};

//        double ratio[maxVal];

        Array1D_int sum(maxVal);

        Array1D ratio(maxVal);



        int cumsum = 0;

        int overMax = 0;

        

        //Ekvalizace:



        for (int i=0; i<N1; i++)

        {

            for (int j=0; j<N2; j++)

            {

                if(u[i][j]>maxVal)

                {

                    u[i][j] = maxVal;

                }

                sum[int(u[i][j])] ++;

            }

        }

        

        for (int i=0; i<= maxVal; i++)

        {

            cumsum = cumsum + sum[i];

            ratio[i] = float(cumsum)/float(N1*N2);

        }

        

        for (int i=0; i<N1; i++)

        {

            for (int j=0; j<N2; j++)

            {

                if(u[i][j]<=maxVal)

                    u[i][j] = ceil(ratio[int(u[i][j])]*maxVal);

                else

                {

                    u[i][j]=maxVal;

                }

            }

        }

        

        

    }

}









void ekvalizaceVyrez(Array2D_int& u, int N1, int N2, int n1, int n2, int rozsah)

{

    // EKVALIZACE : priprava

    int maxVal,mi,mj;

    

    maxVal = u[0][0];



    for (int i=0; i<N1; i++)

    {

        for (int j=0; j<N2; j++)

        {

             if(maxVal<u[i][j])

             {

                 maxVal = u[i][j];

             }

        }

    }



    // MSVC does not support variable-length arrays

//    int sum[maxVal] = {};

//    double ratio[maxVal];

    Array1D_int sum(maxVal + 1);

    Array1D ratio(maxVal + 1);



    int cumsum = 0;



    //Ekvalizace:

    if(maxVal<=200)

    {

        for (int i=n1; i<=n1+rozsah; i++)

        {

            for (int j=n2; j<n2 + rozsah; j++)

            {

                sum[int(u[i][j])] ++;

            }

        }



        for (int i=0; i<= maxVal; i++)

        {

            cumsum = cumsum + sum[i];

            ratio[i] = float(cumsum)/float((rozsah+1)*(rozsah+1));

        }

        for (int i=n1; i<=n1+rozsah; i++)

        {

            for (int j=n2; j<n2 + rozsah; j++)

            {

                u[i][j] = ceil(ratio[int(u[i][j])]*maxVal);

            }

        }

        

    }

    else

    {

        

        maxVal = 200;

        

        // MSVC does not support variable-length arrays

//        int sum[maxVal] = {};

//        double ratio[maxVal];

        Array1D_int sum(maxVal + 1);

        Array1D ratio(maxVal + 1);



        int cumsum = 0;

                

        //Ekvalizace:

        for (int i=n1; i<=n1+rozsah; i++)

        {

            for (int j=n2; j<n2 + rozsah; j++)

            {

                if(u[i][j]>maxVal)

                {

                    u[i][j] = 0;

                }

                sum[int(u[i][j])] ++;

            }

        }



        for (int i=0; i<= maxVal; i++)

        {

            cumsum = cumsum + sum[i];

            ratio[i] = float(cumsum)/float((rozsah+1)*(rozsah+1));

        }

        for (int i=n1; i<=n1+rozsah; i++)

        {

            for (int j=n2; j<n2 + rozsah; j++)

            {

                u[i][j] = ceil(ratio[int(u[i][j])]*maxVal);

            }

        }

    }

}









// u: obrazek

// N1: pocet radku

// N2: pocet sloupcu

void filtr(Array2D& u, int N1, int N2)

{

	double t=1.5e-8;

	float e=0.0001;



    // MSVC does not support variable-length arrays

//	float v[N1][N2];		// N1 radku na N2 sloupcu

	Array2D v(N1, Array1D(N2));

	float h1=1.0/fmax(N1, N2);

	float h2=h1;

    

        

    // Difuze:

    for (int n=1; n<=60; n++)

    {



        #pragma omp parallel for collapse(2) schedule(static)

        for (int i=1; i<N1-1;i++){

                for (int j=1; j<N2-1; j++){

                    double E=pow(((u[i+1][j]-u[i][j])/h1),2)+pow(((u[i][j+1]+u[i+1][j+1]-u[i][j-1]-u[i+1][j-1])/(4*h2)),2);

                    double N=pow(((u[i+1][j-1]+u[i+1][j]-u[i-1][j-1]-u[i-1][j])/(4*h1)),2)+pow(((u[i][j]-u[i][j-1])/h2),2);

                    double W=pow(((u[i][j]-u[i-1][j])/h1),2)+pow(((u[i-1][j+1]+u[i][j+1]-u[i-1][j-1]-u[i][j-1])/(4*h2)),2);

                    double S=pow(((u[i+1][j]+u[i+1][j+1]-u[i-1][j]-u[i-1][j+1])/(4*h1)),2)+pow(((u[i][j+1]-u[i][j])/h2),2);



                    double g1=sqrt(1/(e*e+E));

                    double g2=sqrt(1/(e*e+N));

                    double g3=sqrt(1/(e*e+W));

                    double g4=sqrt(1/(e*e+S));



                    double norma=(sqrt(E)+sqrt(N)+sqrt(W)+sqrt(S))/4;





                    v[i][j]=u[i][j]+ norma*t*((1/(h1*h1))*(g1*(u[i+1][j]-u[i][j])-g3*(u[i][j]-u[i-1][j]))

                            +(1/(h2*h2))*(-g2*(u[i][j]-u[i][j-1])+g4*(u[i][j+1]-u[i][j])));

                        }

                    }



        // okraje

            #pragma omp parallel for schedule(static)

            for (int j=1; j<N2-1; j++) {

                v[0][j]=u[1][j];

                v[N1-1][j]=u[N1-2][j];

            }

            #pragma omp parallel for schedule(static)

            for (int i=0; i<N1; i++) {

                v[i][0]=u[i][1];

                v[i][N2-1]=u[i][N2-2];

            }

            // (rohy neni treba resit samostatne, kdyz se predchozi cyklus projede od 0 do N1-1)



        // zamena matic

            #pragma omp parallel for collapse(2) schedule(static)

            for (int k=0; k<N1;k++){

                for (int l=0; l<N2; l++){

                    u[k][l]=v[k][l];

                        }

                    }



    }

	

 

}





// u: obrazek

// N1: pocet radku

// N2: pocet sloupcu

void filtrZaver(Array2D_int& u, int N1, int N2)

{

	double t=1.0e-10;

	float e=0.0001;



    // MSVC does not support variable-length arrays

//	float v[N1][N2];		// N1 radku na N2 sloupcu

	Array2D v(N1, Array1D(N2));

	float h1=1.0/fmax(N1, N2);

	float h2=h1;

    

        

    // Difuze:

    for (int n=1; n<=30; n++)

    {

        #pragma omp parallel for collapse(2) schedule(static)

        for (int i=1; i<N1-1;i++){

                for (int j=1; j<N2-1; j++){

                    double E=pow(((u[i+1][j]-u[i][j])/h1),2)+pow(((u[i][j+1]+u[i+1][j+1]-u[i][j-1]-u[i+1][j-1])/(4*h2)),2);

                    double N=pow(((u[i+1][j-1]+u[i+1][j]-u[i-1][j-1]-u[i-1][j])/(4*h1)),2)+pow(((u[i][j]-u[i][j-1])/h2),2);

                    double W=pow(((u[i][j]-u[i-1][j])/h1),2)+pow(((u[i-1][j+1]+u[i][j+1]-u[i-1][j-1]-u[i][j-1])/(4*h2)),2);

                    double S=pow(((u[i+1][j]+u[i+1][j+1]-u[i-1][j]-u[i-1][j+1])/(4*h1)),2)+pow(((u[i][j+1]-u[i][j])/h2),2);



                    double g1=sqrt(1/(e*e+E));

                    double g2=sqrt(1/(e*e+N));

                    double g3=sqrt(1/(e*e+W));

                    double g4=sqrt(1/(e*e+S));



                    double norma=(sqrt(E)+sqrt(N)+sqrt(W)+sqrt(S))/4;





                    v[i][j]=u[i][j]+ norma*t*((1/(h1*h1))*(g1*(u[i+1][j]-u[i][j])-g3*(u[i][j]-u[i-1][j]))

                            +(1/(h2*h2))*(-g2*(u[i][j]-u[i][j-1])+g4*(u[i][j+1]-u[i][j])));

                        }

                    }



        // okraje

            #pragma omp parallel for schedule(static)

            for (int j=1; j<N2-1; j++) {

                v[0][j]=u[1][j];

                v[N1-1][j]=u[N1-2][j];

            }

            #pragma omp parallel for schedule(static)

            for (int i=0; i<N1; i++) {

                v[i][0]=u[i][1];

                v[i][N2-1]=u[i][N2-2];

            }

            // (rohy neni treba resit samostatne, kdyz se predchozi cyklus projede od 0 do N1-1)



        // zamena matic

            #pragma omp parallel for collapse(2) schedule(static)

            for (int k=0; k<N1;k++){

                for (int l=0; l<N2; l++){

                    u[k][l]=v[k][l];

                        }

                    }



    }

	

 

}

addon/mri.cpp

0 → 100644
+193 −0
Original line number Original line Diff line number Diff line 
#include <napi.h>

#include <stdio.h>

#include <iostream>

#include "filtr.h"

#include "segmentFunkce.h"

#include "distFunkce.h"

#include "optfFunkce.h"



Napi::ArrayBuffer level_set(const Napi::CallbackInfo & info) {

    Napi::Env env = info.Env();

    float polomer_vnejsi = info[0].As<Napi::Number>();

    float polomer_vnitrni = info[1].As<Napi::Number>();

    int stred1 = info[2].As<Napi::Number>();

    int stred2 = info[3].As<Napi::Number>();

    float K1 = info[4].As<Napi::Number>();

    float K2 = info[5].As<Napi::Number>();

    const int N1 = info[6].As<Napi::Number>();

    const int N2 = info[7].As<Napi::Number>();



    Napi::ArrayBuffer buf = info[8].As<Napi::ArrayBuffer>();

    double * data = reinterpret_cast<double*>(buf.Data());



    Array2D u(N1, Array1D(N2));

    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            u[i][j] = data[i * N2 + j];



    filtr(u, N1, N2);



    const int n1 = stred1 - polomer_vnejsi - 15; // svisle

    const int n2 = stred2 - polomer_vnejsi - 15; // vodorovne

    const int rozsah = 2*(polomer_vnejsi + 15);

    const double h1 = 1.0 / std::max(N1, N2);

    const double h2 = h1;



    Array2D_int A(N1, Array1D_int(N2));

    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            A[i][j] = std::round(u[i][j]);



    Array2D grad(N1, Array1D(N2));

    Array2D u_vnejsi(N1, Array1D(N2));

    Array2D u_vnitrni(N1, Array1D(N2));

    Array2D v(N1, Array1D(N2));



//filtrovani druha faze ________________________



    ekvalizaceVyrez(A,N1,N2,n1,n2,rozsah);

    filtrZaver(A,N1,N2);

//_____________________________________________-



    // vnejsi

//    float K = 1.5e-6;

    getGrad(A, grad, N1, N2, h1, h2, K1);



    inicializace(u_vnejsi, stred1, stred2, polomer_vnejsi, N1, N2, h1);



    segment(u_vnejsi, v, grad, n1, n2, rozsah, h1, h2, 6000);

    distance(u_vnejsi, n1, n2, N1, N2, rozsah, h1, h2, 3.0e4);



    // modifikace - vnejsi -> vnitrni

    for (int i=n1; i<=n1 + rozsah; i++)

    for (int j=n2; j<=n2 + rozsah; j++)

        u_vnitrni[i][j] = u_vnejsi[i][j] + (polomer_vnejsi - polomer_vnitrni) * h1;





    // vnitrni

//    K = 4.0e-6;

    getGrad(A, grad, N1, N2, h1, h2, K2);

    segment(u_vnitrni, v, grad, n1, n2, rozsah, h1, h2, 1200);

    // now u_vnitrni == v



    // mezikruzi

    for (int i=n1; i<=n1 + rozsah; i++)

    for (int j=n2; j<=n2 + rozsah; j++)

        if(u_vnitrni[i][j]<0)

            v[i][j] = -u_vnitrni[i][j];     // vnitrek

        else

            v[i][j] = u_vnejsi[i][j];       // vnejsek



    distance(v, n1, n2, N1, N2, rozsah, h1, h2, 3.0e4);



    // normalizace pro vystup

    #pragma omp parallel for collapse(2) schedule(static)

    for (int i=0; i<N1; i++)

    for (int j=0; j<N2; j++) {

        if(v[i][j]> 10*h1)

            v[i][j] = 10*h1;

        if(v[i][j]< -10*h1)

            v[i][j] = -10*h1;

    }



    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            data[i * N2 + j] = v[i][j];



    return Napi::ArrayBuffer::New(env, buf.Data(), buf.ByteLength());

}



Napi::ArrayBuffer OPTF(const Napi::CallbackInfo & info) {

    Napi::Env env = info.Env();

    float polomer = info[0].As<Napi::Number>();

    int stred1 = info[1].As<Napi::Number>();

    int stred2 = info[2].As<Napi::Number>();

    const int N1 = info[3].As<Napi::Number>();

    const int N2 = info[4].As<Napi::Number>();



    Napi::ArrayBuffer buf = info[5].As<Napi::ArrayBuffer>();

    double * data = reinterpret_cast<double*>(buf.Data());



    Array2D J(N1, Array1D(N2));

    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            J[i][j] = data[i * N2 + j];



    Array2D I(N1, Array1D(N2));

    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            I[i][j] = data[i * N2 + j + N1 * N2];



    Array2D_int T(N1, Array1D_int(N2));

    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            T[i][j] = std::round(data[i * N2 + j + N1 * N2 * 2]);



    const int n1 = stred1 - polomer - 15; // svisle

    const int n2 = stred2 - polomer - 15; // vodorovne

    const int rozsah = 2*(polomer + 15);

    const double h1 = 1.0 / std::max(N1, N2);

    const double h2 = h1;



    Array2D u(N1, Array1D(N2));

    Array2D v(N1, Array1D(N2));



    // pomocna pole

    Array2D_int p(N1, Array1D_int(N2));



    for (int i = 0; i < N1; i++)

    for (int j = 0; j < N2; j++) {

        u[i][j] = 0.0;

        v[i][j] = 0.0;

        J[i][j] = std::max(N1, N2) * J[i][j];

        I[i][j] = std::max(N1, N2) * I[i][j];

    }



    for (int i = n1 + 1; i < n1 + rozsah; i++)

    for (int j = n2 + 1; j < n2 + rozsah; j++){

        u[i][j] =  0.1 * h1;

        v[i][j] = -0.1 * h1;

    }



    const float gama = 4.5;   // -> minimalizace

    const float alfa = 0.75;  // -> brightness const.

    const float beta = 4.5;   // -> vyhlazení



    const float e = 1.0e-3;

    const float t = 7.0e-8;



    const int pocet_iteraci = 6.0e3;



    optf(u, v, I, J, n1, n2, N1, N2, rozsah, h1, h2, pocet_iteraci, e, t, alfa, beta, gama);



    #pragma omp parallel for collapse(2) schedule(static)

    for (int i = 0; i < N1; i++)

    for (int j = 0; j < N2; j++) {

        const int ra = i - std::round(v[i][j] / h2);

        const int sl = j - std::round(u[i][j] / h1);

        if (ra >= 0 && ra < N1 && sl >= 0 && sl < N2)

            p[i][j] = T[ra][sl];

        else

            p[i][j] = 0;

    }



    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            data[i * N2 + j] = p[i][j];

    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            data[i * N2 + j + N1 * N2] = u[i][j];

    for (int i = 0; i < N1; i++)

        for (int j = 0; j < N2; j++)

            data[i * N2 + j + N1 * N2 * 2] = v[i][j];



    return Napi::ArrayBuffer::New(env, buf.Data(), buf.ByteLength());

}



Napi::Object InitAll(Napi::Env env, Napi::Object exports) {

  exports.Set("levelSet", Napi::Function::New(env, level_set));

  exports.Set("optf", Napi::Function::New(env, OPTF));

  return exports;

}



NODE_API_MODULE(mriAddon, InitAll)
 
 No newline at end of file